lmcurve2 - Levenberg-Marquardt least-squares fit of a curve (t,y,dy) (obsolescent)
lmcurve2() is obsolescent, and may be withdrawn in lmfit 13. Its own outputs are consistent with each other, but the convention they follow is hard-wired and unstated, and it is not the convention that lmmin2(3) applies to its parerr. For the full account, see ERROR ESTIMATES below.
The recommended replacement is lmcurve3(3), which takes one additional argument, scale_cov, selecting the convention explicitly. Migration preserves all numbers exactly: insert scale_cov=0 after covar.
lmcurve2() itself is unchanged, in signature and in results, and remains available throughout lmfit 11 and 12.
#include <lmcurve2.h>
void lmcurve2( const int n_par, double *par, double *parerr, double *covar, const int m_dat, const double *t, const double *y, const double *dy, double (*f)( const double ti, const double *par ), const lm_control_struct *control, lm_status_struct *status);
extern const lm_control_struct lm_control_double;
extern const lm_control_struct lm_control_float;
extern const char *lm_infmsg[];
extern const char *lm_shortmsg[];
lmcurve2() wraps the more generic minimization function lmmin2(), for use in curve fitting.
lmcurve2() determines a vector par that minimizes the sum of squared elements of a residue vector r[i] := (y[i] - f(t[i];par)) / dy[i]. Typically, lmcurve2() is used to approximate a data set t,y,dy, where dy represents the standard deviation of empirical data y, by a parametric function f(ti;par). On success, par represents a local minimum, not necessarily a global one; it may depend on its starting value. Users must ensure that all dy[i] are positive.
Function arguments:
Number of free variables. Length of parameter vector par.
Parameter vector. On input, it must contain a reasonable guess. On output, it contains the solution found to minimize ||r||.
Parameter uncertainties vector. Array of length n_par or NULL. On output, unless it is NULL, it contains
parerr[j] = sqrt( covar[j*n_par+j] )i.e. the uncertainties for the weighted model, on the assumption that the dy are the standard deviations of the data. See ERROR ESTIMATES.
Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains the inverse curvature matrix
covar = (J^T J)^-1where J is the Jacobian at the solution, computed by forward differences. Since the residues are divided by dy here, this is the covariance matrix of the parameters, in their own units. See ERROR ESTIMATES.
Both parerr and covar are set to zero if no error estimate can be obtained, which is the case if the sum of squares underflows, if f cannot be evaluated during the finite differencing, or if the curvature matrix is singular.
Number of data points. Length of vectors t, y, dy. Must statisfy n_par <= m_dat.
Array of length m_dat. Contains the abcissae (time, or "x") for which function f will be evaluated.
Array of length m_dat. Contains the ordinate values that shall be fitted.
Array of length m_dat. Contains the standard deviations of the values y.
A user-supplied parametric function f(ti;par).
Parameter collection for tuning the fit procedure. In most cases, the default &lm_control_double is adequate. If f is only computed with single-precision accuracy, &lm_control_float should be used. Parameters are explained in lmmin2(3).
A record used to return information about the minimization process: For details, see lmmin2(3).
Parameter uncertainties are not defined absolutely. They are defined relative to an assumption about the errors of the data, and there are two such assumptions in common use; lmmin3(3) explains both. lmcurve2() implements the first one, and offers no way to request the second:
Trusted sigma, which lmcurve2() assumes. The dy are the standard deviations of the data. Since the residues are divided by dy, the covariance matrix is (J^T J)^-1, which is what covar contains, and parerr is the square root of its diagonal. Both are then in the units of the parameters.
Inferred sigma, which lmcurve2() cannot deliver. The dy are trusted only up to a common factor, which is deduced from the residues themselves under the assumption that the model is correct. Both outputs would then have to be multiplied by status.s2 resp. its square root -- status.s2 = fnorm^2/(m_dat-n_par) being the reduced chi-squared. A caller can do this by hand, but the function does not say that it has to be done, and the omission is silent: the numbers look plausible either way.
Two further points to be aware of:
The parerr of lmcurve2() does not mean the same as the parerr of lmmin2(3). The latter carries an additional factor sqrt(status.s2), because that function mixes the two conventions. Comparing the two functions on the same fit, or moving code from one to the other, changes the reported uncertainties by that factor.
If your dy are relative weights rather than standard deviations -- for instance all set to 1 because only the shape of the data is trusted -- then the trusted-sigma convention is the wrong one, and the uncertainties reported here are too small or too large by sqrt of the reduced chi-squared.
lmcurve3(3) makes the choice an explicit argument, and preserves the identity parerr[j] = sqrt(covar[j*n_par+j]) in both conventions.
Fit a data set y(x) with standard deviations dy(x) by a curve f(x;p):
#include "lmcurve2.h"
#include <stdio.h>
/* model function: a parabola */
double f( double t, const double *p )
{
return p[0] + p[1]*t + p[2]*t*t;
}
int main()
{
int n = 3; /* number of parameters in model function f */
double par[3] = { 100, 0, -10 }; /* really bad starting value */
double parerr[3];
double covar[3*3];
/* data points: a slightly distorted standard parabola */
int m = 9;
int i;
double t[9] = { -4., -3., -2., -1., 0., 1., 2., 3., 4. };
double y[9] = { 16.6, 9.9, 4.4, 1.1, 0., 1.1, 4.2, 9.3, 16.4 };
double dy[9] = { 4, 3, 2, 1, 2, 3, 4, 5, 6 };
lm_control_struct control = lm_control_double;
lm_status_struct status;
control.verbosity = 1;
printf( "Fitting ...\n" );
/* now the call to lmfit */
lmcurve2( n, par, parerr, covar, m, t, y, dy, f, &control, &status );
printf( "Results:\n" );
printf( "status after %d function evaluations:\n %s\n",
status.nfev, lm_infmsg[status.outcome] );
printf("obtained parameters:\n");
for ( i = 0; i < n; ++i)
printf(" par[%i] = %12g uncertainty = %12g\n", i, par[i], parerr[i]);
printf("obtained norm:\n %12g\n", status.fnorm );
printf("fitting data as follows:\n");
for ( i = 0; i < m; ++i)
printf(
" t[%1d]=%2g y=%5.1f+-%4.1f fit=%8.5f residue=%8.4f weighed=%8.4f\n",
i, t[i], y[i], dy[i], f(t[i],par), y[i] - f(t[i],par),
(y[i] - f(t[i],par))/dy[i] );
return 0;
}Copyright (C) 2009-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH
Software: FreeBSD License
Documentation: Creative Commons Attribution Share Alike
Homepage: https://jugit.fz-juelich.de/mlz/lmfit
Please send bug reports and suggestions to the author <j.wuttke@fz-juelich.de>.