[R] Derivative of a splinefun function.
Martin Maechler
maechler at stat.math.ethz.ch
Sat Mar 18 14:45:28 CET 2006
>>>>> "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
More information about the R-help
mailing list