secondDeriv function

Numeric second derivatives

Computes numeric second derivatives (hessian) of an arbitrary multidimensional function at a particular location.

secondDeriv(x, FUN, eps = 1e-08, ...)

Arguments

• x: The location (a vector) where the second derivatives of FUN are desired.

• FUN: An R function for which the second derivatives are sought. This must be a function of the form FUN <- function(x, ...){...} where x is a vector of variable parameters to FUN at which to evaluate the 2nd derivative, and ... are additional parameters needed to evaluate the function. FUN must return a single value (scalar), the height of the surface above x, i.e., FUN evaluated at x.

• eps: A vector of small relative distances to add to x

  when evaluating derivatives. This determines the '$dx$' of the numerical derivatives. That is, the function is evaluated at x, x+dx, and x+2*dx, where $dx$ = x*eps^0.25, in order to compute the second derivative. eps defaults to 1e-8 for all dimensions which equates to setting $dx$ to one percent of each x (i.e., by default the function is evaluate at x, 1.01*x and 1.02*x

  to compute the second derivative).

  One might want to change eps if the scale of dimensions in x varies wildly (e.g., kilometers and millimeters), or if changes between FUN(x) and FUN(x*1.01) are below machine precision. If length of eps is less than length of x, eps is replicated to the length of x.

• ...: Any arguments passed to FUN.

Details

This function uses the "5-point" numeric second derivative method advocated in numerous numerical recipe texts. During computation of the 2nd derivative, FUN must be capable of being evaluated at numerous locations within a hyper-ellipsoid with cardinal radii 2*x(eps)^0.25 = 0.02x at the default value of eps.

A handy way to use this function is to call an optimization routine like nlminb with FUN, then call this function with the optimized values (solution) and FUN. This will yield the hessian at the solution and this is can produce a better estimate of the variance-covariance matrix than using the hessian returned by some optimization routines. Some optimization routines return the hessian evaluated at the next-to-last step of optimization.

An estimate of the variance-covariance matrix, which is used in Rdistance, is solve(hessian) where hessian is secondDeriv(<parameter estimates>, <likelihood>).

Examples

func <- function(x){-x*x} # second derivative should be -2
secondDeriv(0,func)
secondDeriv(3,func)

func <- function(x){3 + 5*x^2 + 2*x^3} # second derivative should be 10+12x
secondDeriv(0,func)
secondDeriv(2,func)

func <- function(x){x[1]^2 + 5*x[2]^2} # should be rbind(c(2,0),c(0,10))
secondDeriv(c(1,1),func)