R devel - Oct 2008 - splinefun gives incorrect derivs when extrapolating to the (PR#13138)

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

=========================== 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

--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
==================================================================---
src/library/stats/R/splinefun.R	(revision 46635)
+++ src/library/stats/R/splinefun.R	(working copy)
@@ -38,7 +38,6 @@
 	    y <- as.vector(tapply(y,x,ties))# as.v: drop dim & dimn.
 	    x <- sort(ux)
 	    nx <- length(x)
-	    rm(ux)
 	} else {
 	    o <- order(x)
 	    x <- x[o]
@@ -93,7 +92,7 @@
 	    d=double(nx),
 	    e=double(if(iMeth == 1) nx else 0),
 	    PACKAGE="stats")
-    rm(x,y,nx,o,method,iMeth)
+    rm(x,y,nx,ux,o,method,iMeth)
     z$e <- NULL
     function(x, deriv = 0) {
 	deriv <- as.integer(deriv)
@@ -114,18 +113,31 @@
         ##           where dx := (u[j]-x[i]); i such that x[i] <= u[j] <=
x[i+1},
         ##                u[j]:= xout[j] (unless sometimes for periodic spl.)
         ##           and  d_i := d[i] unless for natural splines at left
-	.C("spline_eval",
-	   z$method,
-	   as.integer(length(x)),
-	   x=as.double(x),
-	   y=double(length(x)),
-	   z$n,
-	   z$x,
-	   z$y,
-	   z$b,
-	   z$c,
-	   z$d,
-	   PACKAGE="stats")$y
+	res <- .C("spline_eval",
+                  z$method,
+                  as.integer(length(x)),
+                  x=as.double(x),
+                  y=double(length(x)),
+                  z$n,
+                  z$x,
+                  z$y,
+                  z$b,
+                  z$c,
+                  z$d,
+                  PACKAGE="stats")$y
+
+        ## deal with points to the left of first knot if natural
+        ## splines are used  (Bug PR#13132)
+        if( deriv > 0 && z$method==2 && any(ind <-
x<=z$x[1]) ){
+          dx <- z$x[1] - x[ind]
+          switch(deriv,
+                 z$c[1] <- 0,  ## values were calculated above from
+                 z$b[1] <- 0,  ## z$d[1], but ought to be zero.
+                 z$y[1] <- 0)
+          res[ind] <- with(z, y[1] + dx*b[1])
+        }
+
+        res
     }
 }
 
Index: src/library/stats/man/splinefun.Rd
==================================================================---
src/library/stats/man/splinefun.Rd	(revision 46635)
+++ src/library/stats/man/splinefun.Rd	(working copy)
@@ -131,7 +131,7 @@
 par(op)
 
 ## 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--