gsl_nls: GSL Nonlinear Least Squares fitting

In gslnls: GSL Nonlinear Least-Squares Fitting

GSL Nonlinear Least Squares fitting

Description

Determine the nonlinear least-squares estimates of the parameters of a nonlinear model using the gsl_multifit_nlinear module present in the GNU Scientific Library (GSL).

Usage

gsl_nls(fn, ...)

## S3 method for class 'formula'
gsl_nls(
  fn,
  data = parent.frame(),
  start,
  algorithm = c("lm", "lmaccel", "dogleg", "ddogleg", "subspace2D"),
  control = gsl_nls_control(),
  jac = NULL,
  fvv = NULL,
  trace = FALSE,
  subset,
  weights,
  na.action,
  model = FALSE,
  ...
)

## S3 method for class 'function'
gsl_nls(
  fn,
  y,
  start,
  algorithm = c("lm", "lmaccel", "dogleg", "ddogleg", "subspace2D"),
  control = gsl_nls_control(),
  jac = NULL,
  fvv = NULL,
  trace = FALSE,
  weights,
  ...
)

Arguments

fn	a nonlinear model defined either as a two-sided formula including variables and parameters, or as a function returning a numeric vector, with first argument the vector of parameters to be estimated. See the individual method descriptions below.
data	an optional data frame in which to evaluate the variables in fn if defined as a formula. Can also be a list or an environment, but not a matrix.
y	numeric response vector if fn is defined as a function, equal in length to the vector returned by evaluation of the function fn.
start	a named list or named numeric vector of starting estimates. start is only allowed to be missing if fn is a selfStart model. If fn is a formula, a naive guess for start is tried, but this should not be relied on.
algorithm	character string specifying the algorithm to use. The following choices are supported: "lm" Levenberg-Marquardt algorithm (default) "lmaccel" Levenberg-Marquardt algorithm with geodesic acceleration. Can be faster than "lm" but less stable. Stability is controlled by the avmax parameter in control, setting avmax to zero is analogous to not using geodesic acceleration. "dogleg" Powell's dogleg algorithm "ddogleg" Double dogleg algorithm, an improvement over "dogleg" by including information about the Gauss-Newton step while the iteration is still far from the minimum. "subspace2D" 2D generalization of the dogleg algorithm. This method searches a larger subspace for a solution, it can converge more quickly than "dogleg" on some problems.
control	an optional list of control parameters to tune the least squares iterations. See gsl_nls_control for the available control parameters and their default values.
jac	either NULL (default) or a function returning the n by p dimensional Jacobian matrix of the nonlinear model fn, where n is the number of observations and p the number of parameters. If a function, the first argument must be the vector of parameters of length p. If NULL, the Jacobian is computed internally using a finite difference approximations. Can also be TRUE, in which case jac is derived symbolically with deriv, this only works if fn is defined as a (non-selfstarting) formula. If fn is a selfStart model, the Jacobian specified in the "gradient" attribute of the self-start model is used instead.
fvv	either NULL (default) or a function returning an n dimensional vector containing the second directional derivatives of the nonlinear model fn, with n the number of observations. This argument is only used if geodesic acceleration is enabled (algorithm = "lmaccel"). If a function, the first argument must be the vector of parameters of length p and the second argument must be the velocity vector also of length p. If NULL, the second directional derivative vector is computed internal using a finite difference approximation. Can also be TRUE, in which case fvv is derived symbolically with deriv, this only works if fn is defined as a (non-selfstarting) formula. If the model function in fn also returns a "hessian" attribute (similar to the "gradient" attribute in a selfStart model), this Hessian matrix is used to evaluate the second directional derivatives instead.
trace	logical value indicating if a trace of the iteration progress should be printed. Default is FALSE. If TRUE, the residual (weighted) sum-of-squares and the current parameter estimates are printed after each iteration.
subset	an optional vector specifying a subset of observations to be used in the fitting process. This argument is only used if fn is defined as a formula.
weights	an optional numeric vector of (fixed) weights. When present, the objective function is weighted least squares.
na.action	a function which indicates what should happen when the data contain NAs. The default is set by the na.action setting of options, and is na.fail if that is unset. The 'factory-fresh' default is na.omit. Value na.exclude can be useful. This argument is only used if fn is defined as a formula.
model	a logical value. If TRUE, the model frame is returned as part of the object. Defaults to FALSE. This argument is only used if fn is defined as a formula.
...	additional arguments passed to the calls of fn, jac and fvv if defined as functions.

Value

If fn is a formula returns a list object of class nls. If fn is a function returns a list object of class gsl_nls. See the individual method descriptions for the structures of the returned lists and the generic functions applicable to objects of both classes.

Methods (by class)

• formula: If fn is a formula, the returned list object is of classes gsl_nls and nls. Therefore, all generic functions applicable to objects of class nls, such as anova, coef, confint, deviance, df.residual, fitted, formula, logLik, nobs, predict, print, profile, residuals, summary, vcov and weights are also applicable to the returned list object. In addition, a method confintd is available for inference of derived parameters.

• function: If fn is a function, the first argument must be the vector of parameters and the function should return a numeric vector containing the nonlinear model evaluations at the provided parameter and predictor or covariate vectors. In addition, the argument y needs to contain the numeric vector of observed responses, equal in length to the numeric vector returned by fn. The returned list object is (only) of class gsl_nls. Although the returned object is not of class nls, the following generic functions remain applicable for an object of class gsl_nls: anova, coef, confint, deviance, df.residual, fitted, formula, logLik, nobs, predict, print, residuals, summary, vcov and weights. In addition, a method confintd is available for inference of derived parameters.

References

M. Galassi et al., GNU Scientific Library Reference Manual (3rd Ed.), ISBN 0954612078.

See Also

nls

https://www.gnu.org/software/gsl/doc/html/nls.html

Examples

# Example 1: exponential model
# (https://www.gnu.org/software/gsl/doc/html/nls.html#exponential-fitting-example)

## data
set.seed(1)
n <- 50
x <- (seq_len(n) - 1) * 3 / (n - 1)
f <- function(A, lam, b, x) A * exp(-lam * x) + b
y <- f(A = 5, lam = 1.5, b = 1, x) + rnorm(n, sd = 0.25)

## model fit
ex1_fit <- gsl_nls(
  fn = y ~ A * exp(-lam * x) + b,                        ## model formula
  data = data.frame(x = x, y = y),                       ## model fit data
  start = c(A = 0, lam = 0, b = 0)                       ## starting values
)
summary(ex1_fit)                                         ## model summary
predict(ex1_fit, interval = "prediction")                ## prediction intervals

## analytic Jacobian 1
gsl_nls(
  fn = y ~ A * exp(-lam * x) + b,                        ## model formula
  data = data.frame(x = x, y = y),                       ## model fit data
  start = c(A = 0, lam = 0, b = 0),                      ## starting values
  jac = function(par) with(as.list(par),                 ## jacobian
    cbind(A = exp(-lam * x), lam = -A * x * exp(-lam * x), b = 1)
  )
)

## analytic Jacobian 2
gsl_nls(
  fn = y ~ A * exp(-lam * x) + b,                        ## model formula
  data = data.frame(x = x, y = y),                       ## model fit data
  start = c(A = 0, lam = 0, b = 0),                      ## starting values
  jac = TRUE                                             ## automatic derivation
)

## self-starting model
gsl_nls(
  fn =  y ~ SSasymp(x, Asym, R0, lrc),                   ## model formula
  data = data.frame(x = x, y = y)                        ## model fit data
)

# Example 2: Gaussian function
# (https://www.gnu.org/software/gsl/doc/html/nls.html#geodesic-acceleration-example-2)

## data
set.seed(1)
n <- 300
x <- seq_len(n) / n
f <- function(a, b, c, x) a * exp(-(x - b)^2 / (2 * c^2))
y <- f(a = 5, b = 0.4, c = 0.15, x) * rnorm(n, mean = 1, sd = 0.1)

## Levenberg-Marquardt (default)
gsl_nls(
  fn = y ~ a * exp(-(x - b)^2 / (2 * c^2)),             ## model formula
  data = data.frame(x = x, y = y),                      ## model fit data
  start = c(a = 1, b = 0, c = 1),                       ## starting values
  trace = TRUE                                          ## verbose output
)

## Levenberg-Marquardt w/ geodesic acceleration 1
gsl_nls(
  fn = y ~ a * exp(-(x - b)^2 / (2 * c^2)),             ## model formula
  data = data.frame(x = x, y = y),                      ## model fit data
  start = c(a = 1, b = 0, c = 1),                       ## starting values
  algorithm = "lmaccel",                                ## algorithm
  trace = TRUE                                          ## verbose output
)

## Levenberg-Marquardt w/ geodesic acceleration 2
## second directional derivative
fvv <- function(par, v, x) {
  with(as.list(par), {
    zi <- (x - b) / c
    ei <- exp(-zi^2 / 2)
    2 * v[["a"]] * v[["b"]] * zi / c * ei + 2 * v[["a"]] * v[["c"]] * zi^2 / c * ei -
      v[["b"]]^2 * a / c^2 * (1 - zi^2) * ei -
      2 * v[["b"]] * v[["c"]] * a / c^2 * zi * (2 - zi^2) * ei -
      v[["c"]]^2 * a / c^2 * zi^2 * (3 - zi^2) * ei
  })
}

## analytic fvv 1
gsl_nls(
  fn = y ~ a * exp(-(x - b)^2 / (2 * c^2)),             ## model formula
  data = data.frame(x = x, y = y),                      ## model fit data
  start = c(a = 1, b = 0, c = 1),                       ## starting values
  algorithm = "lmaccel",                                ## algorithm
  trace = TRUE,                                         ## verbose output
  fvv = fvv,                                            ## analytic fvv
  x = x                                                 ## argument passed to fvv
)

## analytic fvv 2
gsl_nls(
  fn = y ~ a * exp(-(x - b)^2 / (2 * c^2)),             ## model formula
  data = data.frame(x = x, y = y),                      ## model fit data
  start = c(a = 1, b = 0, c = 1),                       ## starting values
  algorithm = "lmaccel",                                ## algorithm
  trace = TRUE,                                         ## verbose output
  fvv = TRUE                                            ## automatic derivation
)

# Example 3: Branin function
# (https://www.gnu.org/software/gsl/doc/html/nls.html#comparing-trs-methods-example)

## Branin model function
branin <- function(x) {
  a <- c(-5.1 / (4 * pi^2), 5 / pi, -6, 10, 1 / (8 * pi))
  f1 <- x[2] + a[1] * x[1]^2 + a[2] * x[1] + a[3]
  f2 <- sqrt(a[4] * (1 + (1 - a[5]) * cos(x[1])))
  c(f1, f2)
}

## Dogleg minimization w/ model as function
gsl_nls(
  fn = branin,                   ## model function
  y = c(0, 0),                   ## response vector
  start = c(x1 = 6, x2 = 14.5),  ## starting values
  algorithm = "dogleg"           ## algorithm
)

gslnls documentation built on Jan. 17, 2023, 5:15 p.m.