Constrained Fitting of a Model to Data
Description
Fitting a model to data, with lower and/or upper bounds
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27  modFit(f, p, ..., lower = Inf, upper = Inf,
method = c("Marq", "Port", "Newton",
"NelderMead", "BFGS", "CG", "LBFGSB", "SANN",
"Pseudo", "bobyqa"), jac = NULL,
control = list(), hessian = TRUE)
## S3 method for class 'modFit'
summary(object, cov=TRUE,...)
## S3 method for class 'modFit'
deviance(object, ...)
## S3 method for class 'modFit'
coef(object, ...)
## S3 method for class 'modFit'
residuals(object, ...)
## S3 method for class 'modFit'
df.residual(object, ...)
## S3 method for class 'modFit'
plot(x, ask = NULL, ...)
## S3 method for class 'summary.modFit'
print(x, digits = max(3, getOption("digits")  3),
...)

Arguments
f 
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

p 
initial values for the parameters to be optimized over. 
... 
additional arguments passed to function 
lower 
lower bounds on the parameters; if unbounded set equal to

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

method 
The method to be used, one of "Marq", "Port", "Newton", "NelderMead", "BFGS", "CG", "LBFGSB", "SANN", "Pseudo", "bobyqa"  see details. 
jac 
A function that calculates the Jacobian; it should be
called as 
hessian 

control 
additional control arguments passed to the optimisation
routine  see details of 
object 
an object of class 
x 
an object of class 
digits 
number of significant digits in printout. 
cov 
when 
ask 
logical; if 
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
:
"NelderMead", the default from
optim
,"BFGS", a quasiNewton method,
"CG", a conjugategradient method,
"LBFGSB", constrained quasiNewton method,
"SANN", method of simulated annealing.
Or one of the following:
"Marq", the LevenbergMarquardt algorithm (
nls.lm
from packageminpack
)  the default. Note that this method is the only least squares method."Newton", a Newtontype algorithm (see
nlm
),"Port", the Port algorithm (see
nlminb
),"Pseudo", a pseudorandomsearch algorithm (see
pseudoOptim
),"bobyqa", derivativefree optimization by quadratic approximation from package
minqa
.
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.
For all methods except LBFGSB
, Port
, Pseudo
, and
bobyqa
that handle box constraints internally, 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 = 1/2 * (upper + lower) + (upper  lower) * 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 * (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
, Marq
or bobyqa
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 LevenbergMarquard 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
.
Value
a list of class modFit containing the results as returned from the called optimisation routines.
This includes the following:
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 
ms 
the mean squared residuals, i.e. 
var_ms 
the weighted and scaled variable mean squared residuals,
one value per observed variable; only when 
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. 
Note: this means that some return arguments of the original optimisation functions are renamed.
More specifically, "objective" and "counts" from routine
nlminb
(method = "Port") are renamed; "value" and
"counts"; "niter" and "minimum" from routine
nls.lm
(method=Marq) are renamed; "counts"
and "value"; "minimum" and "estimate" from routine nlm
(method = "Newton"
) are renamed.
The list returned by modFit
has a method for the
summary
, deviance
, coef
,
residuals
, df.residual
and
print.summary
– see note.
Note
The summary
method is based on an estimate of the
parameter covariance matrix.
In computing the covariance matrix of the fitted parameters, the problem is
treated as if it were a linear least squares problem, linearizing around
the parameter values that minimize Chi^2.
The covariance matrix is estimated as 1/(0.5*Hessian).
This computation relies on several things, i.e.:
the parameter values are located at the minimum (i.e. the fitting algorithm has converged).
the observations y_j are subject to independent errors whose variances are well estimated by 1 / (n  p) times the residual sum of squares (where n = number of data points, p = number of parameters).
the model is not too nonlinear.
This means that the estimated covariance (correlation) matrix and the confidence intervals derived from it may be worthless if the assumptions behind the covariance computation are invalid.
If in doubt about the validity of the summary computations, use Monte Carlo
fitting instead, or run a modMCMC
.
Other methods included are:

deviance
, which returns the model deviance, 
coef
, which extracts the values of the fitted parameters, 
residuals
,which extracts the model residuals, 
df.residual
which returns the residual degrees of freedom 
print.summary
, producing a nice printout of the summary.
Specifying a function to estimate the Jacobian matrix via argument
jac
may increase speed. The Jacobian is used in the methods
"Marq", "BFGS", "CG", "LBFGS", "Port", and is also used at the end,
to estimate the Hessian at the optimum.
Specification of the gradient
in routines "BFGS", "CG", "LBFGS"
from optim
and "port" from nlminb
is not allowed here.
Within modFit
, the gradient is rather estimated from the Jacobian
jac
and the function f
.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>,
Thomas Petzoldt <thomas.petzoldt@tudresden.de>
References
Bates, D., Mullen, K. D. Nash, J. C. and Varadhan, R. 2014. minqa: Derivativefree optimization algorithms by quadratic approximation. R package. https://cran.rproject.org/package=minqa
Gay, D. M., 1990. Usage Summary for Selected Optimization Routines. Computing Science Technical Report No. 153. AT&T Bell Laboratories, Murray Hill, NJ 07974.
Powell, M. J. D. (2009). The BOBYQA algorithm for bound constrained optimization without derivatives. Report No. DAMTP 2009/NA06, Centre for Mathematical Sciences, University of Cambridge, UK. http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2009_06.pdf
Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P., 2007. Numerical Recipes in C. Cambridge University Press.
Price, W.L., 1977. A Controlled Random Search Procedure for Global Optimisation. The Computer Journal, 20: 367370.
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
Please see also additional publications on the help pages of the individual algorithms.
See Also
constrOptim
for constrained optimization.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116  ## =======================================================================
## 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/N01) * 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/N01)*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$Nout) # 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))
## LevenbergMarquardt
print(system.time(
mf < modFit(f = FUN, p = runif(30), jac = JAC)
))
max(abs(mf$par  X))
