numericDeriv: Evaluate Derivatives (Jacobian) Numerically

In ArkaB-DS/nlsj: JN radical alternative to port of stats::nls()

Description

numericDeriv numerically evaluates the Jacobian of an expression for residuals in a nonlinear least squares problem expressed in nls() style.

Usage

numericDeriv(expr, theta, rho = parent.frame(), dir = 1,
             eps = .Machine$double.eps ^ (1/if(central) 3 else 2), central = FALSE)

Arguments

expr	expression or call to be differentiated. Should evaluate to a numeric vector.
theta	character vector of names of numeric variables used in expr.
rho	environment containing all the variables needed to evaluate expr.
dir	numeric vector of directions, typically with values in -1, 1 to use for the finite differences; will be recycled to the length of theta.
eps	a positive number, to be used as unit step size h for the approximate numerical derivative (f(x+h)-f(x))/h or the central version, see central.
central	logical indicating if central divided differences should be computed, i.e., (f(x+h) - f(x-h)) / 2h . These are typically more accurate but need more evaluations of f().

Details

This is an all-R routine that replaces a front end to the C function numeric_deriv, which is described in Writing R Extensions. The description in Writing R Extensions does not quite match the actual code in base R, nor the present nlsj package, since the dir and central arguments are ignored.

The numeric variables must be of type double and not integer.

Value

The value of eval(expr, envir = rho) plus a matrix attribute "gradient". The columns of this matrix are the derivatives of the value with respect to the variables listed in theta. This is structured for use by nls() style Gauus-Newton programs. The term "gradient" is meaningful for derivatives of the individual functions in the residual vector, but together these become the Jacobian matrix of the nonlinear least squares problem.

Author(s)

Saikat DebRoy [e-mail]; tweaks and eps, central options by R Core Team.

Examples

myenv <- new.env()
myenv$mean <- 0.
myenv$sd   <- 1.
myenv$x    <- seq(-3., 3., length.out = 31)
nD <- numericDeriv(quote(pnorm(x, mean, sd)), c("mean", "sd"), myenv)
str(nD)

## Visualize :
require(graphics)
matplot(myenv$x, cbind(c(nD), attr(nD, "gradient")), type="l")
abline(h=0, lty=3)
## "gradient" is close to the true derivatives, you don't see any diff.:
curve( - dnorm(x), col=2, lty=3, lwd=2, add=TRUE)
curve(-x*dnorm(x), col=3, lty=3, lwd=2, add=TRUE)
##
## IGNORE_RDIFF_BEGIN
# shows 1.609e-8 on most platforms
all.equal(attr(nD,"gradient"),
          with(myenv, cbind(-dnorm(x), -x*dnorm(x))))
## IGNORE_RDIFF_END

ArkaB-DS/nlsj documentation built on Dec. 17, 2021, 9:43 a.m.