modFit: Constrained Fitting of a Model to Data

Description

Fitting a model to data, with lower and/or upper bounds

Usage

modFit(f, p, ..., lower = -Inf, upper = Inf, method = c("Marq", "Port", "Newton", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Pseudo"), jac = NULL, control = list(), hessian = TRUE)
"summary"(object, cov=TRUE,...)
"deviance"(object,  ...)
"coef"(object, ...)
"residuals"(object, ...)
"df.residual"(object, ...)
"plot"(x, ask = NULL, ...)
"print"(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

a function to be minimized, with first argument the vector of parameters over which minimization is to take place. It should return either a vector of residuals (of model versus data) or an element of class modCost (as returned by a call to modCost.

initial values for the parameters to be optimized over.

...

additional arguments passed to function f (modFit) or passed to the methods.

lower

lower bounds on the parameters; if unbounded set equal to -Inf.

upper

upper bounds on the parameters; if unbounded set equal to Inf.

method

The method to be used, one of "Marq", "Port", "Newton", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Pseudo" - see details.

jac

A function that calculates the Jacobian; it should be called as jac(x, ...) and return the matrix with derivatives of the model residuals as a function of the parameters. Supplying the Jacobian can substantially improve performance; see last example.

hessian

TRUE if Hessian is to be estimated. Note that, if set to FALSE, then a summary cannot be estimated.

control

additional control arguments passed to the optimisation routine - see details of nls.lm ("Marq"), nlminb ("Port"), optim ("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN"), nlm ("Newton") or pseudoOptim("Pseudo").

object

an object of class modFit.

digits

number of significant digits in printout.

cov

when TRUE also calculates the parameter covariances.

ask

logical; if TRUE, the user is asked before each plot, if NULL the user is only asked if more than one page of plots is necessary and the current graphics device is set interactive, see par(ask=.) and dev.interactive.

Value

par: the best set of parameters found.
ssr: the sum of squared residuals, evaluated at the best set of parameters.
Hessian: A symmetric matrix giving an estimate of the Hessian at the solution found - see note.
residuals: the result of the last f evaluation; that is, the residuals.
ms: the mean squared residuals, i.e. ssr/length(residuals).
var_ms: the weighted and scaled variable mean squared residuals, one value per observed variable; only when f returns an element of class modCost; NA otherwise.
var_ms_unscaled: the weighted, but not scaled variable mean squared residuals
var_ms_unweighted: the raw variable mean squared residuals, unscaled and unweighted.
...: any other arguments returned by the called optimisation routine.

Details

Note that arguments after ... must be matched exactly.

The method to be used is specified by argument method which can be one of the methods from function optim:

"Nelder-Mead", the default from optim,
"BFGS", a quasi-Newton method,
"CG", a conjugate-gradient method,
"L-BFGS-B", constrained quasi-Newton method,
"SANN", method of simulated annealing.

Or one of the following:

"Marq", the Levenberg-Marquardt algorithm (nls.lm from package minpack) - the default. Note that this method is the only least squares method.
"Newton", a Newton-type algorithm (see nlm),
"Port", the Port algorithm (see nlminb),
"Pseudo", a pseudorandom-search algorithm (see pseudoOptim).

For difficult problems it may be efficient to perform some iterations with Pseudo, which will bring the algorithm near the vicinity of a (the) minimum, after which the default algorithm (Marq) is used to locate the minimum more precisely.

The implementation for the routines from optim differs from constrOptim which implements an adaptive barrier algorithm and which allows a more flexible implementation of linear constraints.

In modFit, bounds on parameters are imposed by a transformation of the parameters to be fitted.

In case both lower and upper bounds are specified, this is achieved by a tangens and arc tangens transformation.

This is, parameter values, $p'$, generated by the optimisation routine, and which are located in the range [-Inf, Inf] are transformed, before they are passed to f as:

$$p = (upper + lower)/2 + (upper - lower) \cdot \arctan(p')/\pi$$.

which maps them into the interval [lower, upper].

Before the optimisation routine is called, the original parameter values, as given by argument p are mapped from [lower,upper] to [-Inf, Inf] by:

$$p' = \tan(\pi/2 \cdot (2 p - upper - lower) / (upper - lower))$$

In case only lower or upper bounds are specified, this is achieved by a log transformation and a corresponding exponential back transformation.

In case parameters are transformed (all methods) or in case the method Port, Pseudo or Marq is selected, the Hessian is approximated as $2 * J^T * J$, where J is the Jacobian, estimated by finite differences.

This ignores the second derivative terms, but this is reasonable if the method has truly converged to the minimum.

Note that finite differences are not extremely precise.

In case the Levenberg-Marquard method (Marq) is used, and parameters are not transformed, 0.5 times the Hessian of the least squares problem is returned by nls.lm, the original Marquardt algorithm. To make it compatible, this value is multiplied with 2 and the TRUE Hessian is thus returned by modFit.

References

Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P., 2007. Numerical Recipes in C. Cambridge University Press.

Gay, D. M., 1990. Usage Summary for Selected Optimization Routines. Computing Science Technical Report No. 153. AT&T Bell Laboratories, Murray Hill, NJ 07974. http://netlib.bell-labs.com/cm/cs/cstr/153.pdf Soetaert, K. and Petzoldt, T., 2010. Inverse Modelling, Sensitivity and Monte Carlo Analysis in R Using Package FME. Journal of Statistical Software 33(3) 1--28. http://www.jstatsoft.org/v33/i03

Examples

Run this code


## =======================================================================
## logistic growth model
## =======================================================================
TT    <- seq(1, 60, 5)
N0    <- 0.1
r     <- 0.5
K     <- 100

## perturbed analytical solution
Data <- data.frame(
  time = TT,
     N = K / (1+(K/N0-1) * exp(-r*TT)) * (1 + rnorm(length(TT), sd = 0.01))
)

plot(TT, Data[,"N"], ylim = c(0, 120), pch = 16, col = "red",
     main = "logistic growth", xlab = "time", ylab = "N")


##===================================
## Fitted with analytical solution  #
##===================================

## initial "guess"
parms <- c(r = 2, K = 10, N0 = 5)

## analytical solution
model <- function(parms,time)
  with (as.list(parms), return(K/(1+(K/N0-1)*exp(-r*time))))

## run the model with initial guess and plot results
lines(TT, model(parms, TT), lwd = 2, col = "green")

## FITTING algorithm 1
ModelCost <- function(P) {
 out <- model(P, TT)
 return(Data$N-out)  # residuals
}

(Fita <- modFit(f = ModelCost, p = parms))

times <- 0:60
lines(times, model(Fita$par, times), lwd = 2, col = "blue")
summary(Fita)

##===================================
##  Fitted with numerical solution  #
##===================================

## numeric solution
logist <- function(t, x, parms) {
  with(as.list(parms), {
    dx <- r * x[1] * (1 - x[1]/K)
    list(dx)
  })
}

## model cost,
ModelCost2 <- function(P) {
 out <- ode(y = c(N = P[["N0"]]), func = logist, parms = P, times = c(0, TT))
 return(modCost(out, Data))  # object of class modCost
}

Fit <- modFit(f = ModelCost2, p = parms, lower = rep(0, 3),
              upper = c(5, 150, 10))

out <- ode(y = c(N = Fit$par[["N0"]]), func = logist, parms = Fit$par,
           times = times)

lines(out, col = "red", lty = 2)
legend("right", c("data", "original", "fitted analytical", "fitted numerical"),
       lty = c(NA, 1, 1, 2), lwd = c(NA, 2, 2, 1),
       col = c("red", "green", "blue", "red"), pch = c(16, NA, NA, NA))
summary(Fit)
plot(residuals(Fit))

## =======================================================================
## the use of the Jacobian
## =======================================================================

## We use modFit to solve a set of linear equations
A <- matrix(nr = 30, nc = 30, runif(900))
X <- runif(30)
B <- A %*% X

## try to find vector "X"; the Jacobian is matrix A

## Function that returns the vector of residuals
FUN <- function(x)
  as.vector(A %*% x - B)

## Function that returns the Jacobian
JAC <- function(x) A

## The port algorithm
print(system.time(
  mf <- modFit(f = FUN, p = runif(30), method = "Port")
))
print(system.time(
  mf <- modFit(f = FUN, p = runif(30), method = "Port", jac = JAC)
))
max(abs(mf$par - X))  # should be very small

## BFGS
print(system.time(
  mf <- modFit(f = FUN, p = runif(30), method = "BFGS")
))
print(system.time(
  mf <- modFit(f = FUN, p = runif(30), method = "BFGS", jac = JAC)
))
max(abs(mf$par - X))

## Levenberg-Marquardt
print(system.time(
  mf <- modFit(f = FUN, p = runif(30), jac = JAC)
))
max(abs(mf$par - X))

Run the code above in your browser using DataLab