regressão polinomial
Pessoal tenho o seginte conjunto de dados: f=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25.2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24.4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920) w=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45.0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44.4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000) se vocês plotarem eles, obteram uma curva, gostaria de obter a função desta curva ou uma regressão polinomial de grau elevado dela já tentei usar a função lm e não deu um bom resultado, quando chego em graus elevados e calculo utilizando os coeficientes obtidos os valores dão 10 vezes maiores (não sei bem porque) mas quando ploto ela aparece bonitinha no lugar que deveria aprecer em outro caso tentei na raça mesmo por multiplicação matricial mas também deu errado pois os valores da multiplicação da inversa são muito pequenos alguém pode me dar uma sugestão ou uma luz do que eu estou fazendo de errado? Tito Conte
Certamente ficará mas facil o objetivo para quem for responder se voce mostrar o que está tentando fazer. Envie um CMR (leia o guia de postagem da lista) On Tue, 24 Jul 2012, Tito Conte wrote:
Pessoal tenho o seginte conjunto de dados:
f=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25. 2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24 .4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920) w=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45. 0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44 .4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000) se vocês plotarem eles, obteram uma curva, gostaria de obter a função desta curva ou uma regressão polinomial de grau elevado dela
já tentei usar a função lm e não deu um bom resultado, quando chego em graus elevados e calculo utilizando os coeficientes obtidos os valores dão 10 vezes maiores (não sei bem porque) mas quando ploto ela aparece bonitinha no lugar que deveria aprecer
em outro caso tentei na raça mesmo por multiplicação matricial mas também deu errado pois os valores da multiplicação da inversa são muito pequenos
alguém pode me dar uma sugestão ou uma luz do que eu estou fazendo de errado?
Tito Conte
Segue o código
# pontos
y=c(-26.0000,-25.9862,-25.**9343,-25.8822,-25.8433,-25.** 8054,-25.7948,-25.7872,-25.**7668,-25.7284,-25.7015,-25.** 6282,-25.4612,-25.3016,-25. 2564,-25.2412,-25.2412,-25.**2232,-25.0869,-25.0000,-25.** 0000,-24.9856,-24.9397,-24.**8976,-24.8533,-24.6587,-24.** 6373,-24.5740,-24.5406,-24 .4379,-24.3934,-24.3628,-24.**3277,-24.3042,-24.2972,-24.** 2973,-24.3020,-24.3068,-24.**3068,-24.2956,-24.2920)
x=c(-45.8529,-45.8302,-45.**7575,-45.6680,-45.6107,-45.**
5050,-45.4574,-45.4058,-45.**3538,-45.2913,-45.2634,-45.** 2264,-45.1667,-45.0719,-45. 0213,-45.0000,-45.0000,-44.**9748,-44.8726,-44.8405,-44.** 8405,-44.8351,-44.8259,-44.**7972,-44.7764,-44.6867,-44.** 6674,-44.5912,-44.5587,-44 .4927,-44.4515,-44.4115,-44.**3314,-44.2664,-44.2171,-44.** 1729,-44.1288,-44.0820,-44.**0509,-44.0118,-44.0000)
# regressão polinomial m=lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)) # criar polinomio de ordem 5 # extraindo os coeficientes b0=coefficient(m)[1] b1=coefficient(m)[2] b2=coefficient(m)[3] b3=coefficient(m)[4] b4=coefficient(m)[5] b5=coefficient(m)[6] # recalculando m_teste=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5 A idéia seria que os pontos de m_teste dessem semelhantes a y mas isso não acontece
O seu modelo m retorna NA nos dois ultimos coeficientes, portanto o m_teste será NA. Além disso, coefficient() é uma função de algum pacote específico? O padrão é coef(). Vc tentou predict(m) # ignora os coeficientes com NA plot(x,y) lines(x, predict(m)) ? --- Fernando Mayer Universidade Federal de Santa Catarina - UFSC Departamento de Ecologia e Zoologia - ECZ/CCB URL: http://sites.google.com/site/fernandomayer e-mail: fernandomayer [@] gmail.com 2012/7/24 Tito Conte <tito.conte@gmail.com>:
Segue o código
# pontos
y=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25.
2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24
.4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920)
x=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45.
0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44
.4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000)
# regressão polinomial m=lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)) # criar polinomio de ordem 5
# extraindo os coeficientes b0=coefficient(m)[1] b1=coefficient(m)[2] b2=coefficient(m)[3] b3=coefficient(m)[4] b4=coefficient(m)[5] b5=coefficient(m)[6]
# recalculando m_teste=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5
A idéia seria que os pontos de m_teste dessem semelhantes a y mas isso não acontece
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
Sem testar código segue um comentário genérico: Estes polinomios de alto grau podem gerar variáveis de valores muito elevados e com isto instabiliades numéricas nos calculos matriciais tente usar polinomios ortonormalizados fazendo: y ~ poly(x, 5) Os coeficientes serão diferentes mas preicoes serão as mesmas On Tue, 24 Jul 2012, Tito Conte wrote:
Segue o código # pontos y=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25. 2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24 .4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920)
x=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45. 0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44 .4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000)
# regressão polinomial m=lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)) # criar polinomio de ordem 5
# extraindo os coeficientes b0=coefficient(m)[1] b1=coefficient(m)[2] b2=coefficient(m)[3] b3=coefficient(m)[4] b4=coefficient(m)[5] b5=coefficient(m)[6]
# recalculando m_teste=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5
A idéia seria que os pontos de m_teste dessem semelhantes a y mas isso não acontece
as dicas foram uteis e muito obrigado, muito embora quando eu dou o comando
r=lm(y~poly(x,5)) e depois dou a equação pelos coeficientes b0=r$coefficient[1] b1=r$coefficient[2] b2=r$coefficient[3] b3=r$coefficient[4] b4=r$coefficient[5] b5=r$coefficient[6] m_teste=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5 o m_teste extrapola de novo será que estou extraindo os coeficientes de forma correta? e aplicando a equação de forma correta? vocês podem me dar outra luz? obrigado
O problema que voce relata é pq o polinomio nao foi ajustado com x, x^2 ... x^5 e sim pelo geado pelo poly() POnrtanto seus coneficientes nao sao compatívies com as potencias de x na forma "bruta" Obtenbha os falores ajustados por fitted(r) e prediçoes em outros pontos usando predict() para fazer a predicão On Wed, 25 Jul 2012, Tito Conte wrote:
as dicas foram uteis e muito obrigado, muito embora quando eu dou o comando
r=lm(y~poly(x,5))
e depois dou a equação pelos coeficientes
b0=r$coefficient[1] b1=r$coefficient[2] b2=r$coefficient[3] b3=r$coefficient[4] b4=r$coefficient[5] b5=r$coefficient[6]
m_teste=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5
o m_teste extrapola de novo
será que estou extraindo os coeficientes de forma correta? e aplicando a equação de forma correta?
vocês podem me dar outra luz?
obrigado
Tito seu exemplo ilustra um o problema numérico que pode ocorrer em ajustes. O perigo é que o resultado é prouzido e se não avaliao o erro pode passar desapercebido e a curva de ajuste erraa pode acabar sendo usada. Veja o efeito no gráfico. Teoricamente os ajustes deveriam ser iguais, entretanto numericamente não são. m=lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)) # criar polinomio de ordem 5 coef(m) m1 <- lm(y ~ poly(x, 5)) ## polinomio e ordem 5 com variáveis ortonormalizadas coef(m1) xpred <- data.frame(x = seq(min(x), max(x), len=200)) plot(x, y) lines(xpred$x, predict(m, newdata=xpred)) ## ajuste/predicao com erro numérico lines(xpred$x, predict(m1, newdata=xpred), col=2) ## adequado On Tue, 24 Jul 2012, Tito Conte wrote:
Segue o código # pontos y=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25. 2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24 .4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920)
x=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45. 0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44 .4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000)
# regressão polinomial m=lm(y~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)) # criar polinomio de ordem 5
# extraindo os coeficientes b0=coefficient(m)[1] b1=coefficient(m)[2] b2=coefficient(m)[3] b3=coefficient(m)[4] b4=coefficient(m)[5] b5=coefficient(m)[6]
# recalculando m_teste=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5
A idéia seria que os pontos de m_teste dessem semelhantes a y mas isso não acontece
Paulo agora outra dúvida, tenho dez mil pontos, como comparar os que estão
a direita e a esquerda dessa curva (como se fosse um dentro ou fora do sistema), já que não consigo extrair o coeficiente de da função poly? Necessito dos dados numericamente pois estes pontos serão submetidos a novos tratamentos posteriormente
Seja mais claro na sua colocação e forneça um CMR ilustrando com comentários os pontos em dúvida. À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Walmes vou criar um exemplo para facilitar o entendimento o que eu preciso é comparar pontos amostrados com os resultados de um modelo de interpolação polinomial por exemplo #minhas amostras y=seq(1:9) x=seq(3:11) #pontos do modelo xm=seq(2,11) ym=seq(0:8) #tirando os valores do modelo da reta r=lm(ym~poly(xm,1)) #extraindo coeficientes a=r$coefficients[1] b=r$coefficients[0] validacao=vector() #comparando cada ponto obtido pelo modelo ajustado com a amostra for(i in c(1:length(y)){ # calculando valor do modelo m=a*x[i]+b # comparando modelo com amostra ifelse(m<y[i],a='fora',a='dentro') validacao[i]=a } é algo deste estilo porém com um polinomio de grau elevado, Quando ploto os resultados eles estão bonitos, mas quando calculo ele como o método abaixo os valores extrapolam muito da realidade e não entendo o porque, veja o exemplo aqui y=c(-26.0000,-25.9862,-25.**9343,-25.8822,-25.8433,-25.**
8054,-25.7948,-25.7872,-25.**7668,-25.7284,-25.7015,-25.** 6282,-25.4612,-25.3016,-25. 2564,-25.2412,-25.2412,-25.**2232,-25.0869,-25.0000,-25.** 0000,-24.9856,-24.9397,-24.**8976,-24.8533,-24.6587,-24.** 6373,-24.5740,-24.5406,-24 .4379,-24.3934,-24.3628,-24.**3277,-24.3042,-24.2972,-24.** 2973,-24.3020,-24.3068,-24.**3068,-24.2956,-24.2920)
x=c(-45.8529,-45.8302,-45.**7575,-45.6680,-45.6107,-45.**
5050,-45.4574,-45.4058,-45.**3538,-45.2913,-45.2634,-45.** 2264,-45.1667,-45.0719,-45. 0213,-45.0000,-45.0000,-44.**9748,-44.8726,-44.8405,-44.** 8405,-44.8351,-44.8259,-44.**7972,-44.7764,-44.6867,-44.** 6674,-44.5912,-44.5587,-44 .4927,-44.4515,-44.4115,-44.**3314,-44.2664,-44.2171,-44.** 1729,-44.1288,-44.0820,-44.**0509,-44.0118,-44.0000)
r=lm(y~poly(x,5)) b0=r$coefficient[1] b1=r$coefficient[2] b2=r$coefficient[3] b3=r$coefficient[4] b4=r$coefficient[5] b5=r$coefficient[6] modelo=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5 modelo veja que os resultados estão muito fora. mas deveriam ser os mesmos, como proceder?
Voce parece estar complicando o fácil. A sua logica de comparar o dado com o ajustao aplicando a equacao caso a caso pode ser alterada por algo mais simples O objheto com modelo retorna residuos, valores preditos me portanto voce pode obter isto diretamente bastando tomar os residuos que sao positivos ou negativos. alem disto o loop for é desnecessário resid(r) ## mostra residuos ifelse(resid(r) > 0, "fora", "dentro") On Wed, 25 Jul 2012, Tito Conte wrote:
Walmes vou criar um exemplo para facilitar o entendimento o que eu preciso é comparar pontos amostrados com os resultados de um modelo de interpolação polinomial
por exemplo
#minhas amostras y=seq(1:9) x=seq(3:11)
#pontos do modelo xm=seq(2,11) ym=seq(0:8)
#tirando os valores do modelo da reta r=lm(ym~poly(xm,1))
#extraindo coeficientes a=r$coefficients[1] b=r$coefficients[0]
validacao=vector() #comparando cada ponto obtido pelo modelo ajustado com a amostra for(i in c(1:length(y)){ # calculando valor do modelo m=a*x[i]+b # comparando modelo com amostra ifelse(m<y[i],a='fora',a='dentro') validacao[i]=a }
é algo deste estilo porém com um polinomio de grau elevado, Quando ploto os resultados eles estão bonitos, mas quando calculo ele como o método abaixo os valores extrapolam muito da realidade e não entendo o porque, veja o exemplo aqui
y=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25. 2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24 .4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920)
x=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45. 0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44 .4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000)
r=lm(y~poly(x,5))
b0=r$coefficient[1] b1=r$coefficient[2] b2=r$coefficient[3] b3=r$coefficient[4] b4=r$coefficient[5] b5=r$coefficient[6]
modelo=b0+b1*x+b2*x^2+b3*x^3+b4*x^4+b5*x^5
modelo
veja que os resultados estão muito fora. mas deveriam ser os mesmos, como proceder?
Paulo mas o polinomio corresponde a uma borda por isso dentro fora o trabalho que estou realizando corresponde a saber se as partículas jogadas ao acaso em uma determinada área ultrapassam ou não as bordas. Os dados que passei para vocês corresponde a uma dessas bordas, Porém os dados não são responsáveis pela construção de nenhuma dessas bordas
Paulo e Walmes, encontrei a seguinte forma para resolver o problema, gostaria que dessem uma olhada #dados do limite da borda direitra y=c(-26.0000,-25.9862,-25.9343,-25.8822,-25.8433,-25.8054,-25.7948,-25.7872,-25.7668,-25.7284,-25.7015,-25.6282,-25.4612,-25.3016,-25.2564,-25.2412,-25.2412,-25.2232,-25.0869,-25.0000,-25.0000,-24.9856,-24.9397,-24.8976,-24.8533,-24.6587,-24.6373,-24.5740,-24.5406,-24.4379,-24.3934,-24.3628,-24.3277,-24.3042,-24.2972,-24.2973,-24.3020,-24.3068,-24.3068,-24.2956,-24.2920) x=c(-45.8529,-45.8302,-45.7575,-45.6680,-45.6107,-45.5050,-45.4574,-45.4058,-45.3538,-45.2913,-45.2634,-45.2264,-45.1667,-45.0719,-45.0213,-45.0000,-45.0000,-44.9748,-44.8726,-44.8405,-44.8405,-44.8351,-44.8259,-44.7972,-44.7764,-44.6867,-44.6674,-44.5912,-44.5587,-44.4927,-44.4515,-44.4115,-44.3314,-44.2664,-44.2171,-44.1729,-44.1288,-44.0820,-44.0509,-44.0118,-44.0000) # partícula 1 amostra=c(x=-45.3,y=-25.111) # inserir dados das amostras no limite x[length(x)+1]=amostra['x'] y[length(y)+1]=amostra['y'] #criar polinomio m=lm(y~poly(x,3)) xpred <- data.frame(x = seq(min(x), max(x), len=250)) #plotar dados plot(y,x) lines(xpred$x, predict(m, newdata=xpred)) if(abs(resid(m[length(m)]))<0.1 or resid(m[length(m)]))<0 ) i='dentro' else i='fora'
se "direita"e "esquerda" se referem a residuaos posotivos e negativos use os residuos do modelo ajustado resid(OBJ do Modelo) On Wed, 25 Jul 2012, Tito Conte wrote:
Paulo agora outra dúvida, tenho dez mil pontos, como comparar os que estão a direita e a esquerda dessa curva (como se fosse um dentro ou fora do sistema), já que não consigo extrair o coeficiente de da função poly? Necessito dos dados numericamente pois estes pontos serão submetidos a novos tratamentos posteriormente
participantes (4)
-
Fernando Mayer -
Paulo Justiniano -
Tito Conte -
Walmes Zeviani