how can I evaluate the derivative of a spline function in R?

Question

R can generate a spline function using splinefun() in the splines library. However, I need to evaluate this function at its first and second derivatives. Is there a way to do this?

for example

library(splines)
x <- 1:10
y <- sin(pi/x) #just an example
f_of_x <- splinefun(x,y)

How can I evaluate f'(x) for a vector of x's?

Answer 1

It is very easy to do since the ability to evaluate the function at its derivatives is built in to the function!

 f_of_x(x, deriv = 1)

Thanks R-core!

Answer 2

You may also be interested in the TkSpline function in the TeachingDemos package which will plot the spline function along with its derivatives.

Answer 3

The use of the deriv = argument to splinefun is sensible, and it should be added that second and third derivatives are supposed to be available, but if you work through the examples you will realize that the linear approximations are jagged and or discontinuous at higher degrees.

In the situation in which you have an analytical expression there are some admittedly limited provisions for algorithmic differentiation. See the help(deriv) page for more details.

> deriv(~sin(pi/x), "x")
expression({
    .expr1 <- pi/x
    .value <- sin(.expr1)
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- -(cos(.expr1) * (pi/x^2))
    attr(.value, "gradient") <- .grad
    .value
})

And then constructing "by hand" a second function with that result. Or you could use the DD example provided on the help(deriv) page to automate the process a bit more:

 DD <- function(expr,name, order = 1) {
    if(order < 1) stop("'order' must be >= 1")
    if(order == 1) D(expr,name)
    else DD(D(expr, name), name, order - 1)
 }
 DD(expression(sin(pi/x)), "x", 2)
-(sin(pi/x) * (pi/x^2) * (pi/x^2) - cos(pi/x) * (pi * (2 * x)/(x^2)^2))
 DD(expression(sin(pi/x)), "x")
-(cos(pi/x) * (pi/x^2))
funD<- function(x){}
body(funD) <- DD(expression(sin(pi/x)), "x")
funD
   #function (x) 
     #-(cos(pi/x) * (pi/x^2))
funD(2)
#   [1] -4.809177e-17  as it should be at a maximum
funDD <- function(x){}
body(funDD) <- DD(expression(sin(pi/x)), "x", 2)
funDD(2)
#  [1] -0.6168503   as it should be at a maximum.