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

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

--MP_/kvy20nVajVG/n.8m=_ZjLAX

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

```

> 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

Perhaps I should have thought a bit more about this. For a natural spline c is zero, and d 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) )
```

 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

--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) ) + res[ind] <- ifelse(deriv == 1, z\$y, 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--

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