R/jacobian.R

In TIMP: Fitting Separable Nonlinear Models in Spectroscopy and Microscopy

jacobian <- function(func, nm, x, ...) {
  arglist <- list(...)
  arglist[[nm]] <- x
  f <- do.call(func, arglist)
  n <- length(x)
  eps <- sqrt(.Machine$double.eps)
  df <- matrix(NA, length(f), n)
  newarglist <- arglist
  for (i in 1:n) {
    xx <- abs(x[i])
    if (xx == 0) {
      delta <- eps
    } else {
      delta <- eps * xx
    }
    newx <- x
    newx[i] <- newx[i] + delta
    newarglist[[nm]] <- newx
    df[, i] <- (do.call(func, newarglist) - f) / delta
  }
  df
}