Scipy contains functions for fitting equations with Python, in its scipy.optimize module. The two main ones I’ve used in the past are leastsq and curve_fit, which in itself is a convenience wrapper around leastsq.

curve_fit

For this operation you require three (four) things:

• a function to fit of form f(x, *params)
• x data
• y data
• Optionally error data

You can also supply an initial guess with the p0 argument.

For example a simple linear equation

def f(x, a, b):
    '''
    Fit a linear function y = a * x + b
    '''
    return a * x  + b

Load the data to fit into x, y and e arrays

p0 = [1., 1.]
popt, pcov = curve_fit(f, x, y, sigma=e, p0=p0)

In this example popt contains the “optimal” results, and pcov the covariance array after fitting.

The problem

This method of working is very powerful but you cannot place limits on the extent of the input fitting parameters, and you are locked into using the leastsq underlying function with this nice interface. For example the function

$$
f(x) = x^a
$$

will explode and the fitting routine will complain when asked to fit around a = 0 when x = 0.

lmfit provides this functionality in a convenient object-oriented interface. The function to fit is phrased a little differently but the functionality is the same.

The easiest way is to write out the desired function as with curve_fit and include a “residuals” wrapper to calculate the normalised (if errors are given) distance away from the model

$$
R = \frac{y - m}{\sigma}
$$

where $R$ are the residuals, $y$ are the observed y values, $m$ the values as calculated by the function that’s being fitted, and $\sigma$ the uncertainties. These residuals tested so that the square of the value above is minimised. The functional form of the residuals function is a little different, and parameters must be accessed with a dictionary lookup and the attribute value.

import lmfit

def resids(params, x, y, e):
    '''
    Given the function f outlined above, return the normalised residuals
    '''
    a = params['a'].value
    b = params['b'].value

    model = f(x, a, b)
    return (y - model) / e

# generate the input parameters
params = lmfit.Parameters()
params.add('a', value=1., min=0, max=10)
params.add('b', value=1., min=0, max=10)

out = lmfit.minimize(resids, params, args=(x, y, e))

The above code alters the params object in place so the best fit parameters are given with

print "Best fit a: {}".format(params['a'].value)
print "Best fit b: {}".format(params['b'].value)

The uncertainties can be accessed with the stderr parameter. A nice report can be printed with lmfit.report_fit(params).

Crucially the initial values for the parameter guesses as well as bounds and (something I’ve not explored much) functional relationships between them. This is a simple way of applying some basic Bayesian priors to the parameter assumptions.