LMFIT  Contents  Download  Develop  
NonLinear LeastSquares Minimization and CurveFitting for Python  Introduction  Parameters  Models 
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 or comparing results from related fits.
The minimize() function is a wrapper around Minimizer for running an optimization problem. It takes an objective function (the function that calculates the array to be minimized), a Parameters object, and several optional arguments. See Writing a Fitting Function for details on writing the objective.
find values for the params so that the sumofsquares of the array returned from function is minimized.
Parameters: 


Returns:  MinimizerResult instance, which will contain the optimized parameter, and several goodnessoffit statistics. 
On output, the params will be unchanged. The bestfit values, and where appropriate, estimated uncertainties and correlations, will all be contained in the returned MinimizerResult. See MinimizerResult – the optimization result for further details.
For clarity, it should be emphasized that this function is simply a wrapper around Minimizer that runs a single fit, implemented as:
fitter = Minimizer(fcn, params, fcn_args=args, fcn_kws=kws,
iter_cb=iter_cb, scale_covar=scale_covar, **fit_kws)
return fitter.minimize(method=method)
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:
calculate objective residual to be minimized from parameters.
Parameters: 


Returns:  residual array (generally datamodel) to be minimized in the leastsquares sense. 
Return type:  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 LevenbergMarquardt 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 leastsquares 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.e10:
period = sign(period)*1.e10
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.e10)
but putting this directly in the function with:
if abs(period) < 1.e10:
period = sign(period)*1.e10
is also a reasonable approach. Similarly, one could place bounds on the decay parameter to take values only between pi/2 and pi/2.
By default, the LevenbergMarquardt 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 wellbehaved for most curvefitting needs, and making it easy to estimate uncertainties for and correlations between pairs of fit variables, as discussed in MinimizerResult – the optimization result.
Alternative algorithms can also be used by providing the method keyword to the minimize() function or Minimizer.minimize() class as listed in the Table of Supported Fitting Methods.
Table of Supported Fitting Method, eithers:
Fitting Method method arg to minimize() or Minimizer.minimize() LevenbergMarquardt leastsq NelderMead nelder LBFGSB lbfgsb Powell powell Conjugate Gradient cg NewtonCG newton COBYLA cobyla Truncated Newton tnc Dogleg dogleg Sequential Linear Squares Programming slsqp Differential Evolution differential_evolution
Note
The objective function for the LevenbergMarquardt 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 LevenbergMarquardt method is the method used. Many of the fit statistics and estimates for uncertainties in parameters discussed in MinimizerResult – the optimization result are done only for this method.
An optimization with minimize() or Minimizer.minimize() will return a MinimizerResult object. This is an otherwise plain container object (that is, with no methods of its own) that simply holds the results of the minimization. These results will include several pieces of informational data such as status and error messages, fit statistics, and the updated parameters themselves.
Importantly, the parameters passed in to Minimizer.minimize() will be not be changed. To to find the bestfit values, uncertainties and so on for each parameter, one must use the MinimizerResult.params attribute.
the Parameters actually used in the fit, with updated values, stderr and correl.
ordered list of variable parameter names used in optimization, and useful for understanding the the values in init_vals and covar.
number of function evaluations
boolean (True/False) for whether fit succeeded.
boolean (True/False) for whether uncertainties were estimated.
message about fit success.
integer error value from scipy.optimize.leastsq() (leastsq only).
message from scipy.optimize.leastsq() (leastsq only).
number of variables in fit \(N_{\rm varys}\)
number of data points: \(N\)
degrees of freedom in fit: \(N  N_{\rm varys}\)
residual array, return value of func(): \({\rm Resid}\)
chisquare: \(\chi^2 = \sum_i^N [{\rm Resid}_i]^2\)
reduced chisquare: \(\chi^2_{\nu}= {\chi^2} / {(N  N_{\rm varys})}\)
Akaike Information Criterion statistic (see below)
Bayesian Information Criterion statistic (see below).
Table of Fit Results: These values, including the standard GoodnessofFit statistics, are all attributes of the MinimizerResult object returned by minimize() or Minimizer.minimize().
Attribute Name  Description / Formula 

nfev  number of function evaluations 
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 value of func(): \({\rm Resid}\) 
chisqr  chisquare: \(\chi^2 = \sum_i^N [{\rm Resid}_i]^2\) 
redchi  reduced chisquare: \(\chi^2_{\nu}= {\chi^2} / {(N  N_{\rm varys})}\) 
aic  Akaike Information Criterion statistic (see below) 
bic  Bayesian Information Criterion statistic (see below) 
var_names  ordered list of variable parameter names used for init_vals and covar 
covar  covariance matrix (with rows/columns using var_names 
init_vals  list of initial values for variable parameters 
Note that the calculation of chisquare and reduced chisquare 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\) errorbar) 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.
The MinimizerResult includes the tradtional chisquare and reduced chisquare statistics:
where \(r\) is the residual array returned by the objective function (likely to be (datamodel)/uncertainty for data modeling usages), \(N\) is the number of data points (ndata), and \(N_{\rm varys}\) is number of variable parameters.
Also included are the Akaike Information Criterion, and Bayesian Information Criterion statistics, held in the aic and bic attributes, respectively. These give slightly different measures of the relative quality for a fit, trying to balance quality of fit with the number of variable parameters used in the fit. These are calculated as
Generally, when comparing fits with different numbers of varying parameters, one typically selects the model with lowest reduced chisquare, Akaike information criterion, and/or Bayesian information criterion. Generally, the Bayesian information criterion is considered themost conservative of these statistics.
An iteration callback function is a function to be called at each iteration, just after the objective function is called. The iteration callback allows usersupplied code to be run at each iteration, and can be used to abort a fit.
usersupplied function to be run at each iteration
Parameters: 


Returns:  residual array (generally datamodel) to be minimized in the leastsquares sense. 
Return type:  None for normal behavior, any value like True to abort fit. 
Normally, the iteration callback would have no return value or return None. To abort a fit, have this function return a value that is True (including any nonzero integer). The fit will also abort if any exception is raised in the iteration callback. When a fit is aborted this way, the parameters will have the values from the last iteration. The fit statistics are not likely to be meaningful, and uncertainties will not be computed.
For full control of the fitting process, you’ll want to create a Minimizer object.
creates a Minimizer, for more detailed access to fitting methods and attributes.
Parameters: 


The Minimizer object has a few public methods:
perform fit using either leastsq() or scalar_minimize().
Parameters: 


Returns:  MinimizerResult object, containing updated parameters, fitting statistics, and information. 
Additonal keywords are passed on to the correspond leastsq() or scalar_minimize() method.
perform fit with LevenbergMarquardt 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.e7 Relative error in the approximate solution ftol 1.e7 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
perform fit with any of the scalar minimization algorithms supported by scipy.optimize.minimize().
scalar_minimize() arg Default Value Description method NelderMead fitting method tol 1.e7 fitting and parameter tolerance hess None Hessian of objective function
prepares and initializes model and Parameters for subsequent fitting. This routine prepares the conversion of Parameters into fit variables, organizes parameter bounds, and parses, “compiles” and checks constrain expressions. The method also creates and returns a new instance of a MinimizerResult object that contains the copy of the Parameters that will actually be varied in the fit.
This method is called directly by the fitting methods, and it is generally not necessary to call this function explicitly.
generate and return text of report of bestfit values, uncertainties, and correlations from fit.
Parameters: 


If the first argument is a Parameters object, goodnessoffit statistics will not be included.
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()
p_true.add('amp', value=14.0)
p_true.add('period', value=5.46)
p_true.add('shift', value=0.123)
p_true.add('decay', value=0.032)
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()
fit_params.add('amp', value=13.0)
fit_params.add('period', value=2)
fit_params.add('shift', value=0.0)
fit_params.add('decay', value=0.02)
out = minimize(residual, fit_params, args=(x,), kws={'data':data})
print(fit_report(out))
#<end of examples/doc_withreport.py>
which would write out:
[[Fit Statistics]]
# function evals = 85
# data points = 1001
# variables = 4
chisquare = 498.812
reduced chisquare = 0.500
[[Variables]]
amp: 13.9121944 +/ 0.141202 (1.01%) (init= 13)
period: 5.48507044 +/ 0.026664 (0.49%) (init= 2)
shift: 0.16203677 +/ 0.014056 (8.67%) (init= 0)
decay: 0.03264538 +/ 0.000380 (1.16%) (init= 0.02)
[[Correlations]] (unreported correlations are < 0.100)
C(period, shift) = 0.797
C(amp, decay) = 0.582
C(amp, shift) = 0.297
C(amp, period) = 0.243
C(shift, decay) = 0.182
C(period, decay) = 0.150