Fitting data with SciPy

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:

dataFitted

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

 

 

 

 

About PaulNakroshis

Professor of Physics at the University of Southern Maine.
This entry was posted in curve fitting, data analysis, error bars, Matplotlib, Plotting and tagged , , . Bookmark the permalink.

2 Responses to Fitting data with SciPy

  1. Levi says:

    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?

    • PaulNakroshis says:

      Sure, but because the function blows up when c-a \rightarrow 0, I would suggest taking the natural logarithm of both sides to get

          \[\ln y = -bt -\ln(c-a).\]

      Then, once you obtain the slope (= -b) and the intercept (\ln(c-a)) you will have to do a quick calculation to find c-a. I’m not sure that you will be able to
      obtain c and a individually, however.

Leave a Reply

Your email address will not be published. Required fields are marked *

*