[Rd] bug in spline()? (PR#653)

[Douglas Bates]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._