# Re: [Rd] splinefun gives incorrect derivs when extrapolating to the left (PR#13132)

From: Berwin A Turlach <berwin_at_maths.uwa.edu.au>
Date: Tue, 07 Oct 2008 19:31:03 +0800

G'day Greg,

On Mon, 6 Oct 2008 19:50:13 +0200 (CEST) Greg.Snow_at_imail.org wrote:

> This is a low priority bug that has been around for a while, but I
> came across it again while alpha testing 2.8.

Indeed, that bug must have been around since splinefun was changed to return a function with a deriv argument. Looks as if the person who produced the initial patch didn't really think through all possibilities. ;-)

> The resulting function for splinefun gives incorrect deriv values
> when x is less than the smallest x-value used to create the function
> (at least in one circumstance), but does the correct thing for x
> greater than the largest x-value.

In a nutshell, splinefun calls, after some initial processing, the C routine spline_coef to fit the desired spline. spline_coef returns coefficients for the following representation of the spline;

a[i] + dx*(b[i] + dx*(c[i] + dx* d[i])
where dx := u-x[i]; i such that x[i] <= u <= x[i+1]. This would evaluate the spline at u; the xs are the original knots.

I.e., the four coefficients returned for a knot correspond to the polynomial to the "right" of that knot. Hence no problem with evaluating derivatives to the right of the largest knot.

The routine returned by splinefun, calls the C function spline_eval to evaluate the spline at the desired points. If deriv > 0, then the above representation of the polynomials that make up the splines is differentiated within the R code to obtain the desired derivative of the spline. The problem is that spline_eval does not know which derivative is evaluated and assumes that the function itself is evaluated. Thus, for evaluation points left of the first knot it sets d to zero (d is usually not zero). This, of course, breaks down when derivatives are evaluated (in this case spline_eval should set c, b or a to zero depending on whether the 1, 2 or 3 derivative is evaluated).

The solution would be to either make spline_eval aware that a derivative is evaluated, but this would imply that that C function gets another parameter and the changes in the C code would probably not be very "elegant", or to have the R code of the function returned by splinefun check whether a natural spline was fitted and is now evaluated at points to the left of the first knot; with appropriate actions if this happens.

The attached patch (against the current SVN version of R) implements the latter strategy. With this patch applied, "make check FORCE=FORCE" passes on my machine. The version of R that is build seems to give the correct answer in your example:

```R> x <- 1:10
R> y <- sin(x)
R>
R> splfun <- splinefun(x,y, method='natural')
R>
R> # these should be linear (and are)
R> splfun( seq(0,1, 0.1) )
```

 0.5682923 0.5956102 0.6229280 0.6502459 0.6775638 0.7048816   0.7321995 0.7595174 0.7868352 0.8141531 0.8414710
```R>
R> # these should all be the same
R> splfun( seq(0,1, 0.1), deriv=1 )
```

 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787   0.2731787 0.2731787 0.2731787 0.2731787 0.2731787
```R>
R> # these should all be 0
R> splfun( seq(0,1, 0.1), deriv=2 )
```

 0 0 0 0 0 0 0 0 0 0 0
R> splfun( seq(0,1, 0.1), deriv=3 )
 0 0 0 0 0 0 0 0 0 0 0

HTH. Cheers,

Berwin

• Full address ============================= Berwin A Turlach Tel.: +65 6516 4416 (secr) Dept of Statistics and Applied Probability +65 6516 6650 (self) Faculty of Science FAX : +65 6872 3919 National University of Singapore 6 Science Drive 2, Blk S16, Level 7 e-mail: statba_at_nus.edu.sg Singapore 117546 http://www.stat.nus.edu.sg/~statba

R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Tue 07 Oct 2008 - 11:35:27 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 07 Oct 2008 - 12:30:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.