Skip to content

polynomial issue

4 messages · Ivan Cabrera, Ben Bolker, David Winsemius +1 more

Original
Ivan Cabrera
# 
Hi all, I'm new to any programming and I've been doing some stuff in R  
the last couple of month, I hope you can help me here?

I'm using lm() to do a 12 degree polynomial fit to a series of data  
x,y. Then I built a polynomial using the coefficients given by coef()  
function (replacing any NA by a 0), and evaluate it in x.

 > x<-data[,1]; y <- data[,2]
 > pol <- paste(paste("I(x",1:12,sep="^"),")",sep="")
 > form <- as.formula(paste("y~", paste(pol,collapse="+")))
 > p12 <- lm(form)
 > coef(p12)
   (Intercept)        I(x^1)        I(x^2)        I(x^3)         
I(x^4)        I(x^5)        I(x^6)        I(x^7)        I(x^8)         
I(x^9)       I(x^10)
  1.549683e+05 -2.801410e+02  2.238321e-01 -1.037342e-04  3.071376e-08  
-6.013217e-12  7.735856e-16 -6.168457e-20  2.493099e-24            NA  
-3.184072e-33
       I(x^11)       I(x^12)
            NA  3.973316e-42
 > pl12 <- coef(p12)[1] + coef(p12)[2]*(x^1) + coef(p12)[3]*(x^2) +  
coef(p12)[4]*(x^3) + coef(p12)[5]*(x^4) + coef(p12)[6]*(x^5) +  
coef(p12)[7]*(x^6) + coef(p12)[8]*(x^7) + coef(p12)[9]*(x^8) +  
coef(p12)[11]*(x^10) + coef(p12)[13]*(x^12)

My problem comes when I replace the coef(p12)[ ] with their actual  
values. I get completely different values of the polynomial?

 > pl12.1 <- (1.549683e+05) + (-2.801410e+02)*(x^1) +  
(2.238321e-01)*(x^2) + (-1.037342e-04)*(x^3) + (3.071376e-08)*(x^4) +  
(-6.013217e-12)*(x^5) + (7.735856e-16)*(x^6) + (-6.168457e-20)*(x^7) +  
(2.493099e-24)*(x^8) + (-3.184072e-33)*(x^10) + (3.973316e-42)*(x^12)
 > rle(pl12==pl12.1)
Run Length Encoding
   lengths: int 3412
   values : logi FALSE

and not even the shape of the functions are alike when i plot them.

can anyone give me a clue of what's happening and how to fix it?

Thanks,
Ivan
Ben Bolker
#
On 10-11-05 07:05 PM, Ivan Cabrera wrote:
This is not appropriate for this list (MacOS special interest group),
please post it to the general R help list (r-help at r-project.org).

  cheers
    Ben Bolker
David Winsemius
#
On Nov 5, 2010, at 7:05 PM, Ivan Cabrera wrote:

            

      
      
      
    

        
This does really look more like an Rhelp question.

            

      
      
      
    

        
I am wondering what would happen if you tried:

fit12 <- lm(y ~ poly(x), degree=12)

p12 <- predict(fit12)

plot(p12)

(Untested in absence of reproducible data.)

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ken Beath
    

    

      
      #
    

  


  

      
On 06/11/2010, at 12:30 PM, David Winsemius wrote:

        

            

      
      
      
    

        
Yes,  a numerical problem. The two NA in the coefs make this a certainty. Using (and understanding poly) will solve the problem. However in most cases using splines will produce something closer to the result you want. While a polynomial will fit the points well, it is likely to do ridiculous things between the points.

Further questions should be directed to Rhelp or a general stats list.

Ken