LMFIT Contents Download Develop Non-Linear Least-Squares Minimization and Curve-Fitting for Python Introduction Parameters Models

#### Previous topic

Parameter and Parameters

#### Next topic

Modeling Data and Curve Fitting

# Performing Fits, Analyzing Outputs¶

As shown in the previous chapter, a simple fit can be performed with the minimize() function. For more sophisticated modeling, the Minimizer class can be used to gain a bit more control, especially when using complicated constraints.

## The minimize() function¶

The minimize function takes a objective function (the function that calculates the array to be minimized), a Parameters ordered dictionary, and several optional arguments. See Writing a Fitting Function for details on writing the function to minimize.

minimize(function, params[, args=None[, kws=None[, method='leastsq'[, scale_covar=True[, iter_cb=None[, **leastsq_kws]]]]]])

find values for the params so that the sum-of-squares of the array returned from function is minimized.

Parameters: function (callable.) – function to return fit residual. See Writing a Fitting Function for details. params (dict or Parameters.) – a Parameters dictionary. Keywords must be strings that match [a-z_][a-z0-9_]* and is not a python reserved word. Each value must be Parameter. args (tuple) – arguments tuple to pass to the residual function as positional arguments. kws (dict) – dictionary to pass to the residual function as keyword arguments. method (string (default leastsq)) – name of fitting method to use. See Choosing Different Fitting Methods for details scale_covar (bool (default True)) – whether to automatically scale covariance matrix (leastsq only) iter_cb (callable or None) – function to be called at each fit iteration leastsq_kws (dict) – dictionary to pass to scipy.optimize.leastsq(). Minimizer object, which can be used to inspect goodness-of-fit statistics, or to re-run fit.

On output, the params will be updated with best-fit values and, where appropriate, estimated uncertainties and correlations. See Goodness-of-Fit and estimated uncertainty and correlations for further details.

If provided, the iter_cb function should take arguments of params, iter, resid, *args, **kws, where params will have the current parameter values, iter the iteration, resid the current residual array, and *args and **kws as passed to the objective function.

## Writing a Fitting Function¶

An important component of a fit is writing a function to be minimized – the objective function. Since this function will be called by other routines, there are fairly stringent requirements for its call signature and return value. In principle, your function can be any python callable, but it must look like this:

func(params, *args, **kws):

calculate objective residual to be minimized from parameters.

Parameters: params (dict) – parameters. args – positional arguments. Must match args argument to minimize() kws – keyword arguments. Must match kws argument to minimize() residual array (generally data-model) to be minimized in the least-squares sense. numpy array. The length of this array cannot change between calls.

A common use for the positional and keyword arguments would be to pass in other data needed to calculate the residual, including such things as the data array, dependent variable, uncertainties in the data, and other data structures for the model calculation.

The objective function should return the value to be minimized. For the Levenberg-Marquardt algorithm from leastsq(), this returned value must be an array, with a length greater than or equal to the number of fitting variables in the model. For the other methods, the return value can either be a scalar or an array. If an array is returned, the sum of squares of the array will be sent to the underlying fitting method, effectively doing a least-squares optimization of the return values.

Since the function will be passed in a dictionary of Parameters, it is advisable to unpack these to get numerical values at the top of the function. A simple way to do this is with Parameters.valuesdict(), as with:

def residual(pars, x, data=None, eps=None):
# unpack parameters:
#  extract .value attribute for each parameter
parvals = pars.valuesdict()
period = parvals['period']
shift = parvals['shift']
decay = parvals['decay']

if abs(shift) > pi/2:
shift = shift - sign(shift)*pi

if abs(period) < 1.e-10:
period = sign(period)*1.e-10

model = parvals['amp'] * sin(shift + x/period) * exp(-x*x*decay*decay)

if data is None:
return model
if eps is None:
return (model - data)
return (model - data)/eps


In this example, x is a positional (required) argument, while the data array is actually optional (so that the function returns the model calculation if the data is neglected). Also note that the model calculation will divide x by the value of the ‘period’ Parameter. It might be wise to ensure this parameter cannot be 0. It would be possible to use the bounds on the Parameter to do this:

params['period'] = Parameter(value=2, min=1.e-10)


but putting this directly in the function with:

if abs(period) < 1.e-10:
period = sign(period)*1.e-10


is also a reasonable approach. Similarly, one could place bounds on the decay parameter to take values only between -pi/2 and pi/2.

## Choosing Different Fitting Methods¶

By default, the Levenberg-Marquardt algorithm is used for fitting. While often criticized, including the fact it finds a local minima, this approach has some distinct advantages. These include being fast, and well-behaved for most curve-fitting needs, and making it easy to estimate uncertainties for and correlations between pairs of fit variables, as discussed in Goodness-of-Fit and estimated uncertainty and correlations.

Alternative algorithms can also be used by providing the method keyword to the minimize() function or use the corresponding method name from the Minimizer class as listed in the Table of Supported Fitting Methods.

Table of Supported Fitting Methods:

Fitting Method method arg to minimize() Minimizer method method arg to scalar_minimize()
Levenberg-Marquardt leastsq leastsq() Not available
L-BFGS-B lbfgsb lbfgsb() L-BFGS-B
Powell powell   Powell
Newton-CG newton   Newton-CG
COBYLA cobyla   COBYLA
COBYLA cobyla   COBYLA
Truncated Newton tnc   TNC
Trust Newton-CGn trust-ncg   trust-ncg
Dogleg dogleg   dogleg
Sequential Linear Squares Programming slsqp   SLSQP

Note

The objective function for the Levenberg-Marquardt method must return an array, with more elements than variables. All other methods can return either a scalar value or an array.

Warning

Much of this documentation assumes that the Levenberg-Marquardt method is the method used. Many of the fit statistics and estimates for uncertainties in parameters discussed in Goodness-of-Fit and estimated uncertainty and correlations are done only for this method.

## Goodness-of-Fit and estimated uncertainty and correlations¶

On a successful fit using the leastsq method, several goodness-of-fit statistics and values related to the uncertainty in the fitted variables will be calculated. These are all encapsulated in the Minimizer object for the fit, as returned by minimize(). The values related to the entire fit are stored in attributes of the Minimizer object, as shown in Table of Fit Results while those related to each fitted variables are stored as attributes of the corresponding Parameter.

Table of Fit Results: These values, including the standard Goodness-of-Fit statistics, are all attributes of the Minimizer object returned by minimize().
Minimizer Attribute Description / Formula
nfev number of function evaluations
success boolean (True/False) for whether fit succeeded.
errorbars boolean (True/False) for whether uncertainties were estimated.
ier integer error value from scipy.optimize.leastsq()
lmdif_message message from scipy.optimize.leastsq()
nvarys number of variables in fit $$N_{\rm varys}$$
ndata number of data points: $$N$$
nfree  degrees of freedom in fit: $$N - N_{\rm varys}$$
residual residual array (return of func(): $${\rm Resid}$$
chisqr chi-square: $$\chi^2 = \sum_i^N [{\rm Resid}_i]^2$$
redchi reduced chi-square: $$\chi^2_{\nu}= {\chi^2} / {(N - N_{\rm varys})}$$
var_map list of variable parameter names for rows/columns of covar
covar covariance matrix (with rows/columns using var_map

Note that the calculation of chi-square and reduced chi-square assume that the returned residual function is scaled properly to the uncertainties in the data. For these statistics to be meaningful, the person writing the function to be minimized must scale them properly.

After a fit using using the leastsq() method has completed successfully, standard errors for the fitted variables and correlations between pairs of fitted variables are automatically calculated from the covariance matrix. The standard error (estimated $$1\sigma$$ error-bar) go into the stderr attribute of the Parameter. The correlations with all other variables will be put into the correl attribute of the Parameter – a dictionary with keys for all other Parameters and values of the corresponding correlation.

In some cases, it may not be possible to estimate the errors and correlations. For example, if a variable actually has no practical effect on the fit, it will likely cause the covariance matrix to be singular, making standard errors impossible to estimate. Placing bounds on varied Parameters makes it more likely that errors cannot be estimated, as being near the maximum or minimum value makes the covariance matrix singular. In these cases, the errorbars attribute of the fit result (Minimizer object) will be False.

## Using the Minimizer class¶

For full control of the fitting process, you’ll want to create a Minimizer object, or at least use the one returned from the minimize() function.

class Minimizer(function, params, fcn_args=None, fcn_kws=None, iter_cb=None, scale_covar=True, **kws)

Parameters: function (callable.) – objective function to return fit residual. See Writing a Fitting Function for details. params (dict) – a dictionary of Parameters. Keywords must be strings that match [a-z_][a-z0-9_]* and is not a python reserved word. Each value must be Parameter. fcn_args (tuple) – arguments tuple to pass to the residual function as positional arguments. fcn_kws (dict) – dictionary to pass to the residual function as keyword arguments. iter_cb (callable or None) – function to be called at each fit iteration scale_covar – flag for automatically scaling covariance matrix and uncertainties to reduced chi-square (leastsq only) kws (dict) – dictionary to pass as keywords to the underlying scipy.optimize method. Minimizer object, which can be used to inspect goodness-of-fit statistics, or to re-run fit.

The Minimizer object has a few public methods:

leastsq(scale_covar=True, **kws)

perform fit with Levenberg-Marquardt algorithm. Keywords will be passed directly to scipy.optimize.leastsq(). By default, numerical derivatives are used, and the following arguments are set:

leastsq() arg Default Value Description
xtol 1.e-7 Relative error in the approximate solution
ftol 1.e-7 Relative error in the desired sum of squares
maxfev 2000*(nvar+1) maximum number of function calls (nvar= # of variables)
Dfun None function to call for Jacobian calculation
lbfgsb(**kws)

perform fit with L-BFGS-B algorithm. Keywords will be passed directly to scipy.optimize.fmin_l_bfgs_b().

lbfgsb() arg Default Value Description
factr 1000.0
maxfun 2000*(nvar+1) maximum number of function calls (nvar= # of variables)
fmin(**kws)

perform fit with Nelder-Mead downhill simplex algorithm. Keywords will be passed directly to scipy.optimize.fmin().

fmin() arg Default Value Description
ftol 1.e-4 function tolerance
xtol 1.e-4 parameter tolerance
maxfun 5000*(nvar+1) maximum number of function calls (nvar= # of variables)

perform fit with any of the scalar minimization algorithms supported by scipy.optimize.minimize().

scalar_minimize() arg Default Value Description
tol 1.e-7 fitting and parameter tolerance
hess None Hessian of objective function
prepare_fit(**kws)

prepares and initializes model and Parameters for subsequent fitting. This routine prepares the conversion of Parameters into fit variables, organizes parameter bounds, and parses, checks and “compiles” constrain expressions.

This is called directly by the fitting methods, and it is generally not necessary to call this function explicitly. An exception is when you would like to call your function to minimize prior to running one of the minimization routines, for example, to calculate the initial residual function. In that case, you might want to do something like:

myfit = Minimizer(my_residual, params,  fcn_args=(x,), fcn_kws={'data':data})

myfit.prepare_fit()
init = my_residual(p_fit, x)
pylab.plot(x, init, 'b--')

myfit.leastsq()


That is, this method should be called prior to your fitting function being called.

## Getting and Printing Fit Reports¶

fit_report(params, modelpars=None, show_correl=True, min_correl=0.1)

generate and return text of report of best-fit values, uncertainties, and correlations from fit.

Parameters: params – Parameters from fit, or Minimizer object as returned by minimize(). modelpars – Parameters with “Known Values” (optional, default None) show_correl – whether to show list of sorted correlations [True] min_correl – smallest correlation absolute value to show [0.1]

If the first argument is a Minimizer object, as returned from minimize(), the report will include some goodness-of-fit statistics.

report_fit(params, modelpars=None, show_correl=True, min_correl=0.1)

print text of report from fit_report().

An example fit with report would be

#!/usr/bin/env python
#<examples/doc_withreport.py>

from __future__ import print_function
from lmfit import Parameters, minimize, fit_report
from numpy import random, linspace, pi, exp, sin, sign

p_true = Parameters()

def residual(pars, x, data=None):
vals = pars.valuesdict()
amp =  vals['amp']
per =  vals['period']
shift = vals['shift']
decay = vals['decay']

if abs(shift) > pi/2:
shift = shift - sign(shift)*pi
model = amp * sin(shift + x/per) * exp(-x*x*decay*decay)
if data is None:
return model
return (model - data)

n = 1001
xmin = 0.
xmax = 250.0

random.seed(0)

noise = random.normal(scale=0.7215, size=n)
x     = linspace(xmin, xmax, n)
data  = residual(p_true, x) + noise

fit_params = Parameters()

out = minimize(residual, fit_params, args=(x,), kws={'data':data})

fit = residual(fit_params, x)
print(fit_report(fit_params))

#<end of examples/doc_withreport.py>


which would write out:

[[Variables]]
amp:      13.9121944 +/- 0.141202 (1.01%) (init= 13)
decay:    0.03264538 +/- 0.000380 (1.16%) (init= 0.02)
period:   5.48507044 +/- 0.026664 (0.49%) (init= 2)
shift:    0.16203677 +/- 0.014056 (8.67%) (init= 0)
[[Correlations]] (unreported correlations are <  0.100)
C(period, shift)             =  0.797
C(amp, decay)                =  0.582
`