|
lmfit
|
lmfit -- Non-linear least-squares fit by Levenberg-Marquardt
method.
DESCRIPTION:
Implement Levenberg-Marquardt method to perform a non-linear least
squares fit to a function of an arbitrary number of parameters. The
function may be any non-linear function. If available, partial
derivatives can be calculated by the user function, else this routine
will estimate partial derivatives with a forward difference
approximation.
CATEGORY:
E2 - Curve and Surface Fitting.
SYNTAX:
result= lmfit(f, x, a, y, w, ...);
INPUTS:
F: The model function to fit. The function must be written as
described under RESTRICTIONS, below.
X: Anything useful for the model function, for instance: independent
variables, a complex structure of data or even nothing!. The
LMFIT routine does not manipulate or use values in X, it simply
passes X to the user-written function F.
A: A vector that contains the initial estimate for each parameter.
Y: Array of dependent variables (i.e., the data). Y can have any
geometry, but it must be the same as the result returned by F.
W: Optional weight, must be conformable with Y and all values of W
must be positive or null (default = 1.0). Data points with zero
weight are not fitted. Here are some examples:
- For no weighting (lest square fit): W = 1.0
- For instrumental weighting: W(i) = 1.0/Y(i)
- Gaussian noise: W(i) = 1.0/Var(Y(i))
OUTPUTS:
A: The vector of fitted parameters.
Returns a structure lmfit_result with fields:
NEVAL: (long) number of model function evaluations.
NITER: (long) number of iteration, i.e. successful CHI2 reductions.
NFIT: (long) number of fitted parameters.
NFREE: (long) number of degrees of freedom (i.e., number of valid
data points minus number of fitted parameters).
MONTE_CARLO: (long) number of Monte Carlo simulations.
CHI2_FIRST: (double) starting error value: CHI2=sum(W*(F(X,A)-Y)^2).
CHI2_LAST: (double) last best error value: CHI2=sum(W*(F(X,A)-Y)^2).
CONV: (double) relative variation of CHI2.
SIGMA: (double) estimated uniform standard deviation of data. If a
weight is provided, a value of SIGMA different from one
indicates that, if the model is correct, W should be
multiplied by 1/SIGMA^2. Computed so that
sum(W*(F(X,A)-Y)^2)/SIGMA^2=NFREE.
LAMBDA: (double) last value of LAMBDA.
STDEV: (pointer) standard deviation vector of the parameters.
STDEV_MONTE_CARLO: (pointer) standard deviation vector of the
parameters estimated by Monte Carlo simulations.
CORREL: (pointer) correlation matrice of the parameters.
KEYWORDS:
FIT: List of indices of parameters to fit, the others remaing
constant. The default is to tune all parameters.
CORREL: If set to a non zero and non-nil value, the correlation
matrice of the parameters is stored into LMFIT result.
STDEV: If set to a non zero and non-nil value, the standard deviation
vector of the parameters is stored into LMFIT result.
DERIV: When set to a non zero and non-nil value, indicates that the
model function F is able to compute its derivatives with respect
to the parameters (see RESTRICTIONS). By default, the
derivatives will be estimated by LMFIT using forward difference.
If analytical derivatives are available they should always be
used.
EPS: Small positive value used to estimate derivatives by forward
difference. Must be such that 1.0+EPS and 1.0 are numerically
different and should be about sqrt(machine_precision)/100
(default = 1e-6).
TOL: Stop criteria for the convergence (default = 1e-7). Should not
be smaller than sqrt(machine_precision). The routine returns
when the relative decrease of CHI2 is less than TOL in an
interation.
ITMAX: Maximum number of iterations. Default = 100.
GAIN: Gain factor for tuning LAMBDA (default = 10.0).
LAMBDA: Starting value for parameter LAMBDA (default = 1.0e-3).
MONTE_CARLO: Number of Monte Carlo simulations to perform to estimate
standard deviation of parameters (by default no Monte Carlo
simulations are undergone). May spend a lot of time if you use a
large number; but should not be too small!
GLOBAL VARIABLES:
None.
SIDE EFFECTS:
The values of the vector of parameters A are modified.
PROCEDURE:
The function to be fitted must be defined as follow:
func F(x, a) {....}
and returns a model with same shape as data Y. If you want to provide
analytic derivatives, F should be defined as:
func F(x, a, &grad, deriv=)
{
y= ...;
if (deriv) {
grad= ...;
}
return y;
}
Where X are the independent variables (anything the function needs to
compute synthetic data except the model parameters), A are the model
parameters, DERIV is a flag set to non-nil and non-zero if the
gradient is needed and the output gradient GRAD is a numberof(Y) by
numberof(A) array: GRAD(i,j) = derivative of ith data point model with
respect to jth parameter.
LMFIT tune parameters A so as to minimize: CHI2=sum(W*(F(X,A)-Y)^2).
The Levenberg-Marquardt method consists in varying between the
inverse-Hessian method and the steepest descent method where the
quadratic expansion of CHI2 does not yield a better model. The
initial guess of the parameter values should be as close to the
actual values as possible or the solution may not converge or may give
a wrong answer.
RESTRICTIONS:
Beware that the result does depend on your initial guess A. In the
case of numerous local minima, the only way to get the correct
solution is to start with A close enough to this solution.
The estimates of standard deviation of the parameters are rescaled
assuming that, for a correct model and weights, the expected value of
CHI2 should be of the order of NFREE=numberof(Y)-numberof(A) (LMFIT
actually compute NFREE from the number of valid data points and number
of fitted parameters). If you don't like this you'll have to rescale
the returned standard deviation to meet your needs (all necessary
information are in the structure returned by LMFIT).
EXAMPLE:
This example is from ODRPACK (version 2.01). The function to fit is
of the form:
f(x) = a1+a2*(exp(a3*x)-1.0)^2
Starting guess:
a= [1500.0, -50.0, -0.1];
Independent variables:
x= [ 0.0, 0.0, 5.0, 7.0, 7.5, 10.0,
16.0, 26.0, 30.0, 34.0, 34.5, 100.0];
Data:
y= [1265.0, 1263.6, 1258.0, 1254.0, 1253.0, 1249.8,
1237.0, 1218.0, 1220.6, 1213.8, 1215.5, 1212.0];
Function definition (without any optimization):
func foo(x, a, &grad, deriv=)
{
if (deriv)
grad= [array(1.0, dimsof(x)),
(exp(a(3)*x)-1.0)^2,
2.0*a(2)*x*exp(a(3)*x)*(exp(a(3)*x)-1.0)];
return a(1)+a(2)*(exp(a(3)*x)-1.0)^2;
}
Fitting this model by:
r= lmfit(foo, x, a, y, 1., deriv=1, stdev=1, monte_carlo=500, correl=1)
produces typically the following result:
a = [1264.84, -54.9987, -0.0829835]
r.neval = 12
r.niter = 6
r.nfit = 3
r.nfree = 9
r.monte_carlo = 500
r.chi2_first = 40.4383
r.chi2_last = 40.4383
r.conv = 3.84967e-09
r.sigma = 0.471764
r.lambda = 1e-09
*r.stdev = [1.23727, 1.78309, 0.00575123]
*r.stdev_monte_carlo = [1.20222, 1.76120, 0.00494790]
*r.correl = [[ 1.000, -0.418, -0.574],
[-0.418, 1.000, -0.340],
[-0.574, -0.340, 1.000]]
HISTORY:
- Basic ideas borrowed from "Numerical Recipes in C", CURVEFIT.PRO (an
IDL version by DMS, RSI, of the routine "CURFIT: least squares fit to
a non-linear function", Bevington, Data Reduction and Error Analysis
for the Physical Sciences) and ODRPACK ("Software for Weigthed
Orthogonal Distance Regression" freely available at: www.netlib.org).
- Added: fitting of a subset of the parameters, Monte-Carlo
simulations...
|