Question

Sometimes, we simply want to use non-linear least squares to fit a function to data, perhaps to estimate paramters for a mechanistic or phenomenological model. The curve_fit function uses the quasi-Newton Levenberg-Marquadt aloorithm to perform such fits. Behind the scnees, curve_fit is just a wrapper around the leastsq function that we have already seen in a more conveneint format.

from scipy.optimize import curve_fit

def logistic4(x, a, b, c, d):
    """The four paramter logistic function is often used to fit dose-response relationships."""
    return ((a-d)/(1.0+((x/c)**b))) + d

nobs = 24
xdata = np.linspace(0.5, 3.5, nobs)
ptrue = [10, 3, 1.5, 12]
ydata = logistic4(xdata, *ptrue) + 0.5*np.random.random(nobs)

popt, pcov = curve_fit(logistic4, xdata, ydata)

perr = yerr=np.sqrt(np.diag(pcov))
print 'Param\tTrue\tEstim (+/- 1 SD)'
for p, pt, po, pe  in zip('abcd', ptrue, popt, perr):
    print '%s\t%5.2f\t%5.2f (+/-%5.2f)' % (p, pt, po, pe)

Param       True    Estim (+/- 1 SD)
a   10.00   10.26 (+/- 0.15)
b    3.00    3.06 (+/- 0.76)
c    1.50    1.62 (+/- 0.11)
d   12.00   12.41 (+/- 0.20)

x = np.linspace(0, 4, 100)
y = logistic4(x, *popt)
plt.plot(xdata, ydata, 'o')
plt.plot(x, y);