[R] Fitting spline using Pspline
Ravi Varadhan
rvaradhan at jhmi.edu
Mon May 30 05:11:07 CEST 2011
Yes, you are right that the results of smooth.spline are slightly worse than that of sm.spline.
The Doppler function is "tricky". At small `x' values, it oscillates rapidly. Hence it is not surprising that the smoothers do not do as well.
Here is a noisy version of your Doppler function. I have also considered another smoother `glkerns'. As you can see, the smoothers do better for larger `x' than for small `x'. It is impossible to distinguish changes in function from noise.
require(pspline)
require(lokern)
x=array(0,1000)
y=array(0,1000)
for (i in 1:1000){
x[i] = i/1000
y[i] = (x[i]*(1-x[i]))^.5 * sin(2*pi*(1.05/(x[i]+.05)))
}
y <- y * (1 + rnorm(1000, 0, 0.2))
plot(x,y, cex=0.4, xlim=c(0,0.1))
fit = sm.spline(x, y, norder=2, cv=FALSE)
lines(fit$x,fit$y, col=2)
fit2 = smooth.spline(x, y, cv=FALSE)
lines(fit2$x,fit2$y, col=3)
fit3 = glkerns(x, y)
lines(fit3$x.out,fit3$est, col=4)
Ravi.
________________________________________
From: r-help-bounces at r-project.org [r-help-bounces at r-project.org] On Behalf Of guy33 [david.reshef at magd.ox.ac.uk]
Sent: Sunday, May 29, 2011 6:28 PM
To: r-help at r-project.org
Subject: Re: [R] Fitting spline using Pspline
Ravi,
Thanks so much! You're right, smooth.spline does work on larger n.
Although, for some reason it's results are different (slightly less good?,
but I'm not sure). For example, on the simple doppler function below,
sm.spline seems to be closer to the true function than smooth.spline:
x=array(0,1000)
y=array(0,1000)
for (i in 1:1000){
x[i] = i/1000
y[i] = (x[i]*(1-x[i]))^.5 * sin(2*pi*(1.05/(x[i]+.05)))
}
plot(x,y)
fit = sm.spline(x, y, norder=2, cv=FALSE)
lines(fit$x,fit$y)
fit2 = smooth.spline(x, y, cv=FALSE)
lines(fit2$x,fit2$y)
What do you make of that?
-guy33
--
View this message in context: http://r.789695.n4.nabble.com/Fitting-spline-using-Pspline-tp3559202p3559610.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.
More information about the R-help
mailing list