12

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?

David LeBauer
  • 31,011
  • 31
  • 115
  • 189

3 Answers3

19

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!

Gavin Simpson
  • 170,508
  • 25
  • 396
  • 453
David LeBauer
  • 31,011
  • 31
  • 115
  • 189
2

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

Greg Snow
  • 48,497
  • 6
  • 83
  • 110
2

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.
IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • +1 thanks for the insight; however I am using a spline to approximate a computationally intensive model. – David LeBauer Dec 02 '10 at 22:16
  • 1
    I was afraid of that. I would advise running through the examples on the splinefun page to get a better idea how they behave with higher order (than 1) arguments to derive=. It made sense after I spent some time plotting but at first the results seemed fairly jarring. – IRTFM Dec 02 '10 at 22:32