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

**Douglas Bates
**
bates@stat.wisc.edu

*04 Sep 2000 16:53:14 -0500*

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