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

--MP_/kvy20nVajVG/n.8m=_ZjLAX
Content-Type: text/plain; charset=US-ASCII
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

On Tue, 7 Oct 2008 19:31:03 +0800
Berwin A Turlach <berwin at maths.uwa.edu.au> wrote:
> 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:
Perhaps I should have thought a bit more about this.  For a natural
spline c[1] is zero, and d[1] is typically not but for evaluations
left of the first knot it should be taken as zero.  So the attached
patch solves the problem in what some might consider a more elegant
manner. :)

With the patch "make check FORCE=FORCE" works on my machine and it
also solves 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_/kvy20nVajVG/n.8m=_ZjLAX
Content-Type: text/plain; name=R-patch1
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment; filename=R-patch1

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,25 @@
         ##           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]) )
+          res[ind] <- ifelse(deriv == 1, z$y[1], 0)
+
+        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_/kvy20nVajVG/n.8m=_ZjLAX--