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

From: <berwin_at_maths.uwa.edu.au>
Date: Tue, 07 Oct 2008 13:35:04 +0200 (CEST)


--MP_/Rxf/JAvsQx5JLkhZFc9Jmn4

Content-Type: text/plain; charset=US-ASCII
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

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[1] to zero (d[1] is usually not zero). This, of course, breaks down when derivatives are evaluated (in this case spline_eval should set c[1], b[1] or a[1] 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) )

 [1] 0.5682923 0.5956102 0.6229280 0.6502459 0.6775638 0.7048816  [7] 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 )

 [1] 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787  [7] 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 )

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

HTH. Cheers,         

        Berwin

--MP_/Rxf/JAvsQx5JLkhZFc9Jmn4

Content-Type: text/plain; name=R-patch
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment; filename=R-patch

Index: src/library/stats/R/splinefun.R


Index: src/library/stats/man/splinefun.Rd


 ## Manual spline evaluation --- demo the coefficients : -.x <- get("ux", envir = environment(f)) +.x <- splinecoef$x
 u <- seq(3,6, by = 0.25)
 (ii <- findInterval(u, .x))
 dx <- u - .x[ii]

--MP_/Rxf/JAvsQx5JLkhZFc9Jmn4--



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Tue 07 Oct 2008 - 11:46:14 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.

list of date sections of archive