[R] Polynomial fitting

[David Winsemius]

On Nov 11, 2009, at 6:14 AM, Julia Cains wrote:

> Dear R helpers
>
>
> Suppose I have a following data
>
> y  <- c(9.21, 9.51, 9.73, 9.88, 10.12. 10.21)
>
> t  <- c(0, 0.25, 1, 3, 6, 12)
>
> I want to find out the polynomial which fits y in terms of t i.e. y  
> = f(t) some function of t.
>
> e.g.   y = bo + b1*t + (b2 * t^2) + (b3 * t^3) + ...... and so on.
>
> In Excel I have defined y as independent variable, then defined t,  
> t^2 and t^3 and using regression I could arrive at the equation
>
> y = 9.505799 + (0.191092 * t) - (0.0225 * t^2) + (0.001245 * t^3)
>
> However I feel this is wrong as I am trying to use linear regression  
> but here I am having polynomial in t.

Not sure what you see as wrong. the lm function will give you a least  
squares fit

 > y  <- c(9.21, 9.51, 9.73, 9.88, 10.12, 10.21)  # Note: fixed the  
erroneous period separator
 >
 > t  <- c(0, 0.25, 1, 3, 6, 12)
 > lm(y ~ t + I(t^2) +I(t^3) )

Call:
lm(formula = y ~ t + I(t^2) + I(t^3))

Coefficients:
(Intercept)            t       I(t^2)       I(t^3)
    9.339845     0.332386    -0.047351     0.002143

So those Excel results were only off by roughly 50% on the last two  
coefficients! You should learn to doubt Excel results for statistical  
or regression work. MS continues to ignore the numerous errors  
reported in its statistical routines. It is rather amazing that  
financial institutions continue to use it widely.

I also ran this in OpenOffice.org with its linest function and get  
these estimates which are obviously in reversed order:
0.00262
-0.06546
0.43248
9.30948

So now you have three estimates. I know which one I trust.

In case anyone has doubts then check the plots of predicted that I  
attach. Here's the code that produces it:

 > plot(t, predict(lm(y ~ t + I(t^2) +I(t^3) ) ) , type="l",  
ylim=c(min(y)-.05, max(y)+.05))
 > points(t, y)
 > OOpred <- 0.00262*t^3 -0.06546*t^2 +0.43248*t +9.30948
 > lines(t, OOpred, col="red")
 > Excelpred <- 9.505799 + (0.191092 * t) - (0.0225 * t^2) + (0.001245  
* t^3)
 > lines(t, Excelpred, col="blue")

-------------- next part --------------
A non-text attachment was scrubbed...
Name: polyfit.pdf
Type: application/pdf
Size: 94521 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091111/1567c202/attachment-0002.pdf>
-------------- next part --------------





>
> I am not that good in stats as well as in mathmatics.
>
> I request you to kindly help me as to how to express the 'y' in  
> polynomial in terms of t.
>
> Thanking you in advance
>
> Julia

David Winsemius, MD
Heritage Laboratories
West Hartford, CT