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

This is a low priority bug that has been around for a while, but I came across
it again while alpha testing 2.8.

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.

Here is an example:
> x <- 1:10
> y <- sin(x)
>
> splfun <- splinefun(x,y, method=3D'natural')
>
> # these should be linear (and are)
> splfun( seq(0,1, 0.1) )
 [1] 0.5682923 0.5956102 0.6229280 0.6502459 0.6775638 0.7048816 0.7321995
 [8] 0.7595174 0.7868352 0.8141531 0.8414710
>
> # these should all be the same
> splfun( seq(0,1, 0.1), deriv=3D1 )
 [1] -0.34287807 -0.22582729 -0.12109764 -0.02868912  0.05139826  0.11916451
 [7]  0.17460962  0.21773359  0.24853643  0.26701813 
0.27317870
>
> # these should all be 0
> splfun( seq(0,1, 0.1), deriv=3D2 )
 [1] 1.2321135 1.1089022 0.9856908 0.8624795 0.7392681 0.6160568 0.4928454
 [8] 0.3696341 0.2464227 0.1232114 0.0000000
> splfun( seq(0,1, 0.1), deriv=3D3 )
 [1] -1.232114 -1.232114 -1.232114 -1.232114 -1.232114 -1.232114 -1.232114
 [8] -1.232114 -1.232114 -1.232114 -1.232114
>
>
> # everything works like it should on the right end
> splfun( seq(10,11, 0.1) )
 [1] -0.5440211 -0.6444546 -0.7448881 -0.8453215 -0.9457550 -1.0461885
 [7] -1.1466220 -1.2470554 -1.3474889 -1.4479224
-1.5483559
> splfun( seq(10,11, 0.1), deriv=3D1 )
 [1] -1.004335 -1.004335 -1.004335 -1.004335 -1.004335 -1.004335 -1.004335
 [8] -1.004335 -1.004335 -1.004335 -1.004335
> splfun( seq(10,11, 0.1), deriv=3D2 )
 [1] 0 0 0 0 0 0 0 0 0 0 0
> splfun( seq(10,11, 0.1), deriv=3D3 )
 [1] 0 0 0 0 0 0 0 0 0 0 0
>

--please do not edit the information below--

Version:
 platform =3D i386-pc-mingw32
 arch =3D i386
 os =3D mingw32
 system =3D i386, mingw32
 status =3D alpha
 major =3D 2
 minor =3D 8.0
 year =3D 2008
 month =3D 10
 day =3D 01
 svn rev =3D 46589
 language =3D R
 version.string =3D R version 2.8.0 alpha (2008-10-01 r46589)

Windows XP (build 2600) Service Pack 3

Locale:
LC_COLLATE=3DEnglish_United States.1252;LC_CTYPE=3DEnglish_United
States.1252;LC_MONETARY=3DEnglish_United
States.1252;LC_NUMERIC=3DC;LC_TIME=3DEnglish_United States.1252

Search Path:
 .GlobalEnv, package:stats, package:graphics, package:grDevices, package:utils,
package:datasets, package:methods, Autoloads, package:base



--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111

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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: R-patch
URL:
<https://stat.ethz.ch/pipermail/r-devel/attachments/20081007/8f9eea45/attachment.pl>