[R] Derivative of a splinefun function.

[Martin Maechler]

>>>>> "BeT" == Berwin A Turlach <berwin at maths.uwa.edu.au>
>>>>>     on Sat, 18 Mar 2006 17:08:47 +0800 writes:

    BeT> G'day Rolf,
>>>>> "RT" == Rolf Turner <rolf at math.unb.ca> writes:

    RT> Is there a way of calculating the derivative of a function
    RT> returned by splinefun()?
    BeT> Yes, of course. :)
    BeT> As somebody else once said: "this is R, anything can be done, the
    BeT> question is just how easy it is to do so".

    BeT> It may depend on what kind of spline you are fitting with
    BeT> splinefun()... 

    RT> Such a function is a cubic spline, whence it has a calculable
    RT> derivative, but is there a (simple) way of getting at it?
    BeT> Depends on your definition of simple. ;-))

    RT> One workaround that I have thought of is to take a fine grid
    RT> of points, evaluate the function returned by splinefun() at
    RT> these points, put an interpolating spline through these points
    RT> using smooth.spline() with spar=0, and then calculate the
    RT> derivative with predict(...,deriv=1).
    BeT> Why easy if one can also make it complicated?  There is the function
    BeT> interpSpline() in the package splines for which, according to its help
    BeT> page, the predict method also works.  So using smooth.spline with
    BeT> spar=0 does indeed look a bit kludgy. :)

    RT> Is there a better way?
    BeT> splinefun() essentially returns a function that contains one object in
    BeT> its environment, namely z, which contains all the necessary
    BeT> information for calculating the spline and whose body is just a call
    BeT> to an internal C function.

    BeT> The components y, b, c and d of the object z contain the piecewise
    BeT> polynomial representation between knots of the spline if you choose
    BeT> the methods "fmm" or "natural".  If you use the method "periodic" then
    BeT> the component e of the object z is used too, but I don't know whether
    BeT> it is just used to hold intermediate results while calculating the
    BeT> spline or whether it is actually necessary for calculating the
    BeT> spline.  (An analysis of the C code that is called to calculate the
    BeT> spline and to evaluate the spline should answer that question but
    BeT> until now I never bothered of doing this since I didn't have use for
    BeT> periodic splines yet.)

    BeT> Thus, since the spline is stored in a piecewise polynomial fashion, it
    BeT> shouldn't be hard to define a function that calculates the derivative
    BeT> of the function returned by splinefun().  If my calculus is correct,
    BeT> the following code should construct such a function:

    >> par(mfrow=c(1,2))
    >> n <- 9
    >> x <- 1:n
    >> y <- rnorm(n)
    >> f <- splinefun(x, y)
    >> ls(envir = environment(f))
    > [1] "z"
    >> splinecoef <- eval(expression(z), envir = environment(f))
    >> curve(f(x), 1, 10, col = "green", lwd = 1.5)
    >> points(splinecoef, col = "purple", cex = 2)

    > ##
    > ## construct a function that calculates the derivative of f()
    > ##
    >> g <- f
    >> environment(g) <- new.env(parent=baseenv())
    >> z <- get("z", envir=environment(f)) 
    >> z$y <- z$b
    >> z$b <- 2*z$c
    >> z$c <- 3*z$d
    >> z$d <- rep(0, z$n)
    >> assign("z", z, envir=environment(g))
    >> curve(g(x), 1, 10, col = "green", lwd = 1.5)

    BeT> As indicated above, I would only vouch for this method if splinefun()
    BeT> is called with method="fmm" or method="natural".  And the second
    BeT> caveat is, of course, that you are assuming that I get my calculus
    BeT> correct on a Saturday morning which might be a bit much to ask. ;-))

It does work with "periodic" as well.  The "e" component could
have been zapped (in splinefun before returning).

Berwin, do you feel like making this into patch for splinefun(),
which in the future would return a   function(x, deriv = 0)
and would also work for deriv \in {1,2,3}
(for 'deriv = 3' it could return an object of class "stepfun" !) ?

    BeT> Hope that helps.

I'm pretty sure it did.  Well done for early Saturday morning!
(I've checked numerically with my  sfsmisc::D1ss() function that
 your result is correct)

Regards,
Martin