LMFIT  Contents  Download  Develop  
NonLinear LeastSquares Minimization and CurveFitting for Python  Introduction  Parameters  Models 
As shown in the introduction, 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 takes a function to minimize, a dictionary of Parameter , and several optional arguments. See Writing a Fitting Function for details on writing the function to minimize.
find values for the params so that the sumofsquares of the returned array from function is minimized.
Parameters: 


Returns:  Minimizer object, which can be used to inspect goodnessoffit statistics, or to rerun fit. 
On output, the params will be updated with bestfit values and, where appropriate, estimated uncertainties and correlations. See GoodnessofFit and estimated uncertainty and correlations for further details.
An important component of a fit is writing a function to be minimized in the leastsquares sense. 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 example would look like:
def residual(pars, x, data=None, eps=None):
# unpack parameters:
# extract .value attribute for each parameter
amp = pars['amp'].value
period = pars['period'].value
shift = pars['shift'].value
decay = pars['decay'].value
if abs(shift) > pi/2:
shift = shift  sign(shift)*pi
if abs(period) < 1.e10:
period = sign(period)*1.e10
model = 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 varied value of the ‘period’ Parameter. It might be wise to make sure 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 might be wiser to put this directly in the function with:
if abs(period) < 1.e10:
period = sign(period)*1.e10
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 GoodnessofFit and estimated uncertainty and correlations.
Alternative algorithms can also be used. These include simulated annealing which promises a better ability to avoid local minima, and BFGS, which is a modification of the quasiNewton method.
To select which of these algorithms to use, use 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() LevenbergMarquardt leastsq leastsq() Not available NelderMead nelder fmin() NelderMead LBFGSB lbfgsb lbfgsb() LBFGSB Simulated Annealing anneal anneal() Anneal Powell powell Powell Conjugate Gradient cg CG NewtonCG newton NewtonCG COBYLA cobyla COBYLA Sequential Linear Squares Programming slsqp SLSQP
Note
Use of scipy.optimize.minimize() requires scipy 0.11 or higher.
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
The LevenbergMarquardt method is by far the most tested fit method, and much of this documentation assumes that this is the method used. For example, many of the fit statistics and estimates for uncertainties in parameters discussed in GoodnessofFit and estimated uncertainty and correlations are done only for the leastsq method.
In particular, the simulated annealing method appears to not work correctly.... understanding this is on the ToDo list.
On a successful fit using the leastsq method, several goodnessoffit 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 GoodnessofFit 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. 
message  message about fit success. 
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  chisquare: \(\chi^2 = \sum_i^N [{\rm Resid}_i]^2\) 
redchi  reduced chisquare: \(\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 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.
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.
creates a Minimizer, for finegrain access to fitting methods and attributes.
Parameters: 


Returns:  Minimizer object, which can be used to inspect goodnessoffit statistics, or to rerun fit. 
The Minimizer object has a few public methods:
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 LBFGSB algorithm. Keywords will be passed directly to scipy.optimize.fmin_l_bfgs_b().
lbfgsb() arg Default Value Description factr 1000.0 approx_grad True calculate approximations of gradient maxfun 2000*(nvar+1) maximum number of function calls (nvar= # of variables)
perform fit with NelderMead downhill simplex algorithm. Keywords will be passed directly to scipy.optimize.fmin().
fmin() arg Default Value Description ftol 1.e4 function tolerance xtol 1.e4 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 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, 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.
generate and return text of report of bestfit values, uncertainties, and correlations from fit.
Parameters: 


If the first argument is a Minimizer object, as returned from minimize(), the report will include some goodnessoffit statistics.
print text of report from fit_report().
An example fit with an error report:
p_true = Parameters()
p_true.add('amp', value=14.0)
p_true.add('period', value=5.33)
p_true.add('shift', value=0.123)
p_true.add('decay', value=0.010)
def residual(pars, x, data=None):
amp = pars['amp'].value
per = pars['period'].value
shift = pars['shift'].value
decay = pars['decay'].value
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 = 2500
xmin = 0.
xmax = 250.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})
fit = residual(fit_params, x)
report_errors(fit_params)
would generate this report:
[[Variables]]
amp: 13.969724 +/ 0.050145 (0.36%) initial = 13.000000
decay: 0.009990 +/ 0.000042 (0.42%) initial = 0.020000
period: 5.331423 +/ 0.002788 (0.05%) initial = 2.000000
shift: 0.125333 +/ 0.004938 (3.94%) initial = 0.000000
[[Correlations]] (unreported correlations are < 0.100)
C(period, shift) = 0.800
C(amp, decay) = 0.576