lmcurve3 - Levenberg-Marquardt least-squares fit of a curve (t,y,dy)
#include <lmcurve3.h>
void lmcurve3( const int n_par, double *par, double *parerr, double *covar, const int scale_cov, 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[];
lmcurve3() wraps the more generic minimization function lmmin3(), for use in curve fitting.
It replaces the obsolescent lmcurve2(3), from which it differs only by the additional argument scale_cov; with scale_cov=0 the results are the same.
lmcurve3() 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, lmcurve3() 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 the square roots of the diagonal elements of covar,
parerr[j] = sqrt( covar[j*n_par+j] )in the convention selected by scale_cov. This identity holds in both conventions.
Covariance matrix. Array of length n_par * n_par or NULL. On output, unless it is NULL, it contains
covar = (J^T J)^-1 if scale_cov == 0
covar = (J^T J)^-1 * status.s2 if scale_cov != 0where J is the Jacobian at the solution, computed by forward differences. 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.
Selects the convention in which parerr and covar are expressed. Use 0 if the dy are the standard deviations of the data, nonzero if they are trusted only up to a common factor, which is then inferred from the residues. See ERROR ESTIMATES.
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 lmmin3(3).
A record used to return information about the minimization process: For details, see lmmin3(3). Of particular interest here is status.s2 = status.fnorm^2/(m_dat-n_par), the reduced chi-squared, which is the factor between the two conventions below.
Parameter uncertainties are defined relative to an assumption about the errors of the data, and scale_cov selects that assumption. Write J for the Jacobian of the residue vector r at the solution, and s2 = status.s2 = fnorm^2/(m_dat-n_par).
covar = (J^T J)^-1
parerr[j] = sqrt( covar[j*n_par+j] )Since the residues are divided by dy, this is correct if the dy are the standard deviations of the data. Both outputs are then in the units of the parameters, and they do not depend on how well the model happens to fit. This is the convention of MINPACK and of GSL's gsl_multifit_covar(), and it is what lmcurve2(3) always did.
covar = (J^T J)^-1 * s2
parerr[j] = sqrt( covar[j*n_par+j] )Correct if the dy give only the relative weights of the data points, their common scale being unknown and to be deduced from the residues under the assumption that the model is correct. The typical case is all dy[i] equal, for instance all 1. This is what scipy.optimize.curve_fit does by default (absolute_sigma=False), and what gnuplot's fit and the summaries of R's nls report.
The rule of thumb: if you would be prepared to defend your dy as standard deviations, and to conclude from a reduced chi-squared far from 1 that the model or the errors are wrong, use 0. If your dy are weights, use nonzero.
Since status.s2 is reported either way, a program that needs both can convert instead of fitting twice: divide or multiply parerr by sqrt(status.s2).
For the relation to the obsolescent lmmin2(3) and lmcurve2(3), and for the history of these conventions in lmfit, see ERROR ESTIMATES in lmmin3(3).
Fit a data set y(x) with standard deviations dy(x) by a curve f(x;p):
#include "lmcurve3.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 */
/* dy are standard deviations, so the trusted-sigma
convention applies: scale_cov = 0 */
lmcurve3( n, par, parerr, covar, 0, 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-2026 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>.