Here’s a common thing scientists need to do, and it’s easy to accomplish in python. Suppose that you have a data set consisting of temperature vs time data for the cooling of a cup of coffee. We’ll start by importing the needed libraries and defining a fitting function:

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit def fitFunc(t, a, b, c): return a*np.exp(-b*t) + c

Now we create some fake data as numpy arrays and add some noise; the fake data will be called *noisy:*

t = np.linspace(0,4,50) temp = fitFunc(t, 2.5, 1.3, 0.5) noisy = temp + 0.25*np.random.normal(size=len(temp))

The scipy.optimize module contains a least squares curve fit routine that requires as input a user-defined fitting function (in our case *fitFunc* ), the x-axis data (in our case, *t*) and the y-axis data (in our case, *noisy*). The curve_fit routine returns an array of fit parameters, and a matrix of covariance data (the square root of the diagonal values are the 1-sigma uncertainties on the fit parameters—provided you have a reasonable fit in the first place.):

fitParams, fitCovariances = curve_fit(fitFunc, t, noisy) print fitParams print fitCovariance

output:

[ 2.42573207 1.38417604 0.57396876] [[ 0.01738482 0.00815523 -0.00077816] [ 0.00815523 0.02509109 0.00606179] [-0.00077816 0.00606179 0.00292134]]

Now we plot the data points with error bars, plot the best fit curve, and label the axes:

plt.ylabel('Temperature (C)', fontsize = 16) plt.xlabel('time (s)', fontsize = 16) plt.xlim(0,4.1) # plot the data as red circles with vertical errorbars plt.errorbar(t, noisy, fmt = 'ro', yerr = 0.2) # now plot the best fit curve and also +- 1 sigma curves # (the square root of the diagonal covariance matrix # element is the uncertianty on the fit parameter.) sigma = [fitCovariance[0,0], \ fitCovariance[1,1], \ fitCovariance[2,2] \ ] plt.plot(t, fitFunc(t, fitParams[0], fitParams[1], fitParams[2]),\ t, fitFunc(t, fitParams[0] + sigma[0], fitParams[1] - sigma[1], fitParams[2] + sigma[2]),\ t, fitFunc(t, fitParams[0] - sigma[0], fitParams[1] + sigma[1], fitParams[2] - sigma[2])\ ) # save plot to a file savefig('dataFitted.png', bbox_inches=0)

The output file looks like this:

You can see a the curve fitting routine as a python script, and you can see an ipython notebook here.

I am using your script to fit this function: exp(-b*t)*(c-a)**-1. When running the script it fails because of the (c-a)**-1. I assume that (c-a)**-1 is causing 1/0 when varying a and c. Is there a way to fit this function?

Sure, but because the function blows up when , I would suggest taking the natural logarithm of both sides to get

Then, once you obtain the slope () and the intercept () you will have to do a quick calculation to find . I’m not sure that you will be able to

obtain and individually, however.

Without loss of generality, why not simplify the function to (1/a)*exp(-b*t).

Then everything should work.

Thank you, a very helpful code segment!