[R] Polynomial fitting

[Berry Boessenkool]

Hey Matti,

lm() _is_ a least square approximation.
Did you notivce in the poly-help you'll probably need to set raw=TRUE ?
Helped me a lot once I figured that one out...
check this:

a <- c(-5:10)  ;  b <- 5*a^3 + 2*a^2 - 7*a + 3  + rnorm(15, 0, 100) 
data.frame(a,b) 
modell <- lm(b ~  poly(a,3, raw=T))     ;    modell
plot(a,b)  ; lines(a, predict(modell))

now try this with raw=FALSE. The graph will look the same, but the coefficients are normalized and centered.
Notice that the value estimated for the intercept is quite different form the "real" intercept (3), as the random numbers added have great variation (sd=100).
The other coefficents are closer to the "real" ones. Note also, that they are printed in reverse order than specified in b <- ...

A way to avoid the whole polyfunction altogether is:

lm(b ~ I(a^3) + I(a^2) + a)

You decide what's more elegant...

Hope this helps,
Berry

-------------------------------------
Berry Boessenkool
Potsdam
-------------------------------------


> Date: Fri, 8 Jul 2011 21:12:46 +0300
> From: matti.jokipii at gmail.com
> To: Gerrit.Eichner at math.uni-giessen.de
> CC: r-help at r-project.org
> Subject: Re: [R] Polynomial fitting
> 
> Thank you Gerrit for the quick reply! And yes, i'm Matti.
> 
> I can get the coeffs now, though i'm not sure whether i'm doing 
> something wrong or whether poly is just not the right method for what 
> i'm trying to find. I will look into this more closely and give it 
> another try.
> 
> Is poly best for fitting on noisy data that's been generated by a 
> polynomial and not that good for approximating an arbitrary function? I 
> tried a least squares fitting with a web applet and got all exited 
> because the approximation looked quite promising. I understand that R is
> designed mainly for statistical computing and may not be the best tool 
> for my purposes. Before i look elsewhere i would like to ask if there is
> some other R method i should try, perhaps a least squares approximation?
> 
[[elided Hotmail spam]]
> 
> Matti Jokipii
> 
> 08.07.2011 08:25, Gerrit Eichner kirjoitti:
> > Hello, mfa (Matti?),
> >
> > if x and y contain the coordinates of your data points and k is the
> > wanted polynomial degree, then
> >
> > fit <- lm( y ~ poly( x, k))
> >
> > fits orthonormal polynomials up to degree k to your data. Using
> >
> > dummy.coef( fit)
> >
> > should give the coefficients you are interested in.
> >
> > Hth -- Gerrit
> >
> > On Thu, 7 Jul 2011, mfa wrote:
> >
> >> Hello,
> >>
> >> i'm fairly familiar with R and use it every now and then for math related
> >> tasks.
> >>
> >> I have a simple non polynomial function that i would like to approximate
> >> with a polynomial. I already looked into poly, but was unable to
> >> understand
> >> what to do with it. So my problem is this. I can generate virtually any
> >> number of datapoints and would like to find the coeffs a1, a2, ... up
> >> to a
> >> given degree for a polynomial a1x^1 + a2x^2 + ... that approximates my
> >> simple function. How can i do this with R?
> >>
[[elided Hotmail spam]]
> >>
> >> --
> >> View this message in context:
> >> http://r.789695.n4.nabble.com/Polynomial-fitting-tp3652816p3652816.html
> >> Sent from the R help mailing list archive at Nabble.com.
> >>
> >> ______________________________________________
> >> R-help at r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide
> >> http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >>
> > ---------------------------------------------------------------------
> > Dr. Gerrit Eichner Mathematical Institute, Room 212
> > gerrit.eichner at math.uni-giessen.de Justus-Liebig-University Giessen
> > Tel: +49-(0)641-99-32104 Arndtstr. 2, 35392 Giessen, Germany
> > Fax: +49-(0)641-99-32109 http://www.uni-giessen.de/cms/eichner
> > ---------------------------------------------------------------------
> >
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.