From: <berwin_at_maths.uwa.edu.au>

Date: Tue, 07 Oct 2008 13:35:04 +0200 (CEST)

[1] 0.5682923 0.5956102 0.6229280 0.6502459 0.6775638 0.7048816 [7] 0.7321995 0.7595174 0.7868352 0.8141531 0.8414710

[1] 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787 [7] 0.2731787 0.2731787 0.2731787 0.2731787 0.2731787

[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

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

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

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

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