Confidence interval for exponential curve fit
confidence-interval, numpy, python, scipy
Solution
Gabriel's answer is incorrect. Here in red the 95% confidence band for his data as calculated by GraphPad Prism:
Background: the "confidence interval of a fitted curve" is typically called confidence band. For a 95% confidence band, one can be 95% confident that it contains the true curve. (This is different from prediction bands, shown above in gray. Prediction bands are about future data points. For more details, see, e.g., this page of the GraphPad Curve Fitting Guide.)
In Python, kmpfit can calculate the confidence band for non-linear least squares. Here for Gabriel's example:
from pylab import *
from kapteyn import kmpfit
x, y = np.loadtxt('_exp_fit.txt', unpack=True)
def model(p, x):
a, b, c = p
return a*np.exp(b*x)+c
f = kmpfit.simplefit(model, [.1, .1, .1], x, y)
print f.params
# confidence band
a, b, c = f.params
dfdp = [np.exp(b*x), a*x*np.exp(b*x), 1]
yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model)
scatter(x, y, marker='.', s=10, color='#0000ba')
ix = np.argsort(x)
for i, l in enumerate((upper, lower, yhat)):
plot(x[ix], l[ix], c='g' if i == 2 else 'r', lw=2)
show()
The `dfdp` are the partial derivatives ∂f/∂p of the model f = a*e^(b*x) + c with respect to each parameter p (i.e., a, b, and c). For background, see the kmpfit Tutorial or this page of the GraphPad Curve Fitting Guide. (Unlike my sample code, the kmpfit Tutorial does not use `confidence_band()` from the library but its own, slightly different, implementation.)
Finally, the Python plot matches the Prism one:
Problem
I'm trying to obtain a confidence interval on an exponential fit to some `x,y` data (available here). Here's the MWE I have to find the best exponential fit to the data: ``` from pylab import * from scipy.optimize import curve_fit # Read data. x, y = np.loadtxt('exponential_data.dat', unpack=True) def func(x, a, b, c): '''Exponential 3-param function.''' return a * np.exp(b * x) + c # Find best fit. popt, pcov = curve_fit(func, x, y) print popt # Plot data and best fit curve. scatter(x, y) x = linspace(11, 23, 100) plot(x, func(x, *popt), c='r') show() ``` which produces: How can I obtain the 95% (or some other value) confidence interval on this fit preferably using either pure `python`, `numpy` or `scipy` (which are the packages I already have installed)?