[R-br] regressão polinomial
Paulo Justiniano
paulojus em leg.ufpr.br
Quarta Julho 25 12:04:01 BRT 2012
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
>
>
Mais detalhes sobre a lista de discussão R-br