BUG IN SPLINE()? Version R-1.0.1, system i486,linux If the spline(x,y,method="natural") function is given values outside the range of the data, it does not give a warning. Moreover, the extrapolated value reported is not the ordinate of the natural spline defined by (x,y). Example. Let x <- c(2,5,8,10) and y <- c(1.2266,-1.7606,-0.5051,1.0390). Then interpolate/extrapolate with a natural cubic spline at 1:10 using either spline(x,y,n=10,method="natural",xmin=1,xmax=10) or fn <- splinefun(x,y,method="natural") fn(1:10) This gives c(2.5366,1.2266,-0.0834,...,1.0390). I agree with all values except that at x=1. The interpolation formula in Green and Silverman's book on splines gives the derivative of a natural spline at the first knot as f'(x1)=(f2-f1)/(x2-x1) -(1/6)(x2-x1)d, where d=second derivative at the second knot. For the above example, d=0.7071, f'(x1)=-1.3493, and the natural spline interpolant has value 1.2266-(-1.3493)=2.5759 at x=1, not 2.5366. PS There are different e-mail addresses for bug reports given in the FAQ file and in the help information for bug.report. ************************************************ * I.White * * ICAPB, University of Edinburgh * * Ashworth Laboratories, West Mains Road * * Edinburgh EH9 3JT * * Fax: 0131 667 3210 Tel: 0131 650 5490 * * E-mail: i.m.s.white@ed.ac.uk * ************************************************ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
imsw@holyrood.ed.ac.uk writes:> This gives c(2.5366,1.2266,-0.0834,...,1.0390). I agree with all values > except that at x=1. The interpolation formula in Green and Silverman's > book on splines gives the derivative of a natural spline at the first > knot as f'(x1)=(f2-f1)/(x2-x1) -(1/6)(x2-x1)d, where d=second derivative > at the second knot. For the above example, d=0.7071, f'(x1)=-1.3493, > and the natural spline interpolant has value 1.2266-(-1.3493)=2.5759 > at x=1, not 2.5366.Also, plot(diff(fn(seq(1,3,.01)))) reveals that something is wrong At the other end, plot(diff(fn(seq(8,11,.01)))) looks even weirder, especially if one modifies to y <- c(1.2266,-1.7606,-0.5051,1.01)> PS There are different e-mail addresses for bug reports given in the > FAQ file and in the help information for bug.report.Hmm, and neither looks quite right! One would probably want the address to be in the r-project.org domain (in case we move it away from Copenhagen). "r-bugs@list.r-project.org", however, is none too logical and assumes that there's an r-bugs user or an alias in Zurich which forwards to biostat.dk (there is, it seems). Same thing with r-bugs@r-project.org, since the Wisconsin machine aliases r-bugs to r-bugs@lists.r-project.org. More logical would be to use r-bugs@r-project.org and have the forwarding there go to biostat.dk. Alternatively, it ought to work as r-bugs@bugs.r-project.org, but that doesn't work whenever a RedHat upgrade blows away the resending rules in my sendmail.cf... (the program in r-bugs' .forward needs to execute on a machine that matches the architecture of the binary!) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
imsw@holyrood.ed.ac.uk writes:> BUG IN SPLINE()? > > Version R-1.0.1, system i486,linux > > If the spline(x,y,method="natural") function is given values outside the > range of the data, it does not give a warning. Moreover, the extrapolated > value reported is not the ordinate of the natural spline defined by (x,y). > > Example. Let x <- c(2,5,8,10) and y <- c(1.2266,-1.7606,-0.5051,1.0390). > Then interpolate/extrapolate with a natural cubic spline at 1:10 using > either > > spline(x,y,n=10,method="natural",xmin=1,xmax=10) > > or > > fn <- splinefun(x,y,method="natural") > fn(1:10) > > This gives c(2.5366,1.2266,-0.0834,...,1.0390). I agree with all values > except that at x=1. The interpolation formula in Green and Silverman's > book on splines gives the derivative of a natural spline at the first > knot as f'(x1)=(f2-f1)/(x2-x1) -(1/6)(x2-x1)d, where d=second derivative > at the second knot. For the above example, d=0.7071, f'(x1)=-1.3493, > and the natural spline interpolant has value 1.2266-(-1.3493)=2.5759 > at x=1, not 2.5366.I think the `interpSpline' function in the splines package gives the desired natural interpolation spline. Perhaps we should change the underlying code for the spline function. > library(splines) > x <- c(2,5,8,10) > y <- c(1.2266,-1.7606,-0.5051,1.0390) > ns1 <- interpSpline(y ~ x) > plot(ns1) # gives a plot that looks reasonable > points(x,y) > predict(ns1, 1:10) $x [1] 1 2 3 4 5 6 7 8 9 10 $y [1] 2.5758923 1.2266000 -0.0834080 -1.1577100 -1.7606000 -1.7349409 [7] -1.2378717 -0.5051000 0.2669514 1.0390000 attr(,"class") [1] "xyVector" -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._