Introduction to Scipy: Fitting data

We have talked about the Numpy and Matplotlib libraries, but there is a third library that is invaluable for Scientific Analysis: Scipy. Scipy is basically a very large library of functions that you can use for scientific analysis. A good place to start to find out about the top-level scientific functionality in Scipy is the Documentation.

Examples of the functionality include:

  • Integration (scipy.integrate)
  • Optimization/Fitting (scipy.optimize)
  • Interpolation (scipy.interpolate)
  • Fourier Transforms (scipy.fftpack)
  • Signal Processing (scipy.signal)
  • Linear Algebra (scipy.linalg)
  • Spatial data structures and algorithms (scipy.spatial)
  • Statistics (scipy.stats)
  • Multi-dimensional image processing (scipy.ndimage)

and so on.

This week, we will take a look at how to fit models to data. When analyzing scientific data, fitting models to data allows us to determine the parameters of a physical system (assuming the model is correct).

There are a number of routines in Scipy to help with fitting, but we will use the simplest one, curve_fit, which is imported as follows:

In [1]:
import numpy as np
from scipy.optimize import curve_fit

The full documentation for the curve_fit is available here, and we will look at a simple example here, which involves fitting a straight line to a dataset.

We first create a fake dataset with some random noise:

In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
In [3]:
x = np.random.uniform(0., 100., 100)
y = 3. * x + 2. + np.random.normal(0., 10., 100)
plt.plot(x, y, '.')
Out[3]:
[<matplotlib.lines.Line2D at 0x1186f73c8>]

Let's now imagine that this is real data, and we want to determine the slope and intercept of the best-fit line to the data. We start off by definining a function representing the model:

In [4]:
def line(x, a, b):
    return a * x + b

The arguments to the function should be x, followed by the parameters. We can now call curve_fit to find the best-fit parameters using a least-squares fit:

In [5]:
popt, pcov = curve_fit(line, x, y)

The curve_fit function returns two items, which we can popt and pcov. The popt argument are the best-fit paramters for a and b:

In [6]:
popt
Out[6]:
array([ 2.99069214,  2.45335655])

which is close to the initial values of 3 and 2 used in the definition of y.

The reason the values are not exact is because there are only a limited number of random samples, so the best-fit slope is not going to be exactly those used in the definition of y. The pcov variable contains the covariance matrix, which indicates the uncertainties and correlations between parameters. This is mostly useful when the data has uncertainties.

Let's now try and fit the data assuming each point has a vertical error (standard deviation) of +/-10:

In [7]:
e = np.repeat(10., 100)
plt.errorbar(x, y, yerr=e, fmt="none")
Out[7]:
<Container object of 3 artists>
In [8]:
popt, pcov = curve_fit(line, x, y, sigma=e)
In [9]:
popt
Out[9]:
array([ 2.99069216,  2.45335518])

Now pcov will contain the true variance and covariance of the parameters, so that the best-fit parameters are:

In [10]:
print("a =", popt[0], "+/-", pcov[0,0]**0.5)
print("b =", popt[1], "+/-", pcov[1,1]**0.5)

a = 2.99069215825 +/- 0.0327811500383
b = 2.45335518448 +/- 1.87958698376

We can now plot the best-fit line:

In [11]:
plt.errorbar(x, y, yerr=e, fmt="none")
xfine = np.linspace(0., 100., 100)  # define values to plot the function for
plt.plot(xfine, line(xfine, popt[0], popt[1]), 'r-')
Out[11]:
[<matplotlib.lines.Line2D at 0x11868e4a8>]

You should now be able to fit simple models to datasets! Note that for more complex models, more sophisticated techniques may be required for fitting, but curve_fit will be good enough for most simple cases.

Note that there is a way to simplify the call to the function with the best-fit parameters, which is:

line(x, *popt)

The * notation will expand a list of values into the arguments of the function. This is useful if your function has more than one or two parameters. Hence, you can do:

In [12]:
plt.errorbar(x, y, yerr=e, fmt="none")
plt.plot(xfine, line(xfine, *popt), 'r-')
Out[12]:
[<matplotlib.lines.Line2D at 0x11b9ca4e0>]

Important Note: the way curve_fit determines the uncertainty is to actually renormalize the errors so that the reduced $\chi^2$ value is one, so the magnitude of the errors doesn't matter, only the relative errors. In some fields of science (such as astronomy) we do not renormalize the errors, so for those cases you can specify absolute_sigma=True in order to preserve the original errors.

Exercise 1

In the following code, we generate some random data points:

In [13]:
x = np.random.uniform(0., 10., 100)
y = np.polyval([1, 2, -3], x) + np.random.normal(0., 10., 100)
e = np.random.uniform(5, 10, 100)

Fit a line and a parabola to it and overplot the two models on top of the data:

In [14]:
# your solution here

Exercise 2

As before, we use the data/munich_temperatures_average_with_bad_data.txt file, which gives the temperature in Munich every day for several years:

In [15]:
# The following code reads in the file and removes bad values
import numpy as np
date, temperature = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt', unpack=True)
keep = np.abs(temperature) < 90
date = date[keep]
temperature = temperature[keep]

Fit the following function to the data:

$$f(t) = a~\cos{(2\pi t + b)} + c$$

where $t$ is the time in years. Make a plot of the data and the best-fit model in the range 2008 to 2012. What are the best-fit values of the parameters? What is the overall average temperature in Munich, and what are the typical daily average values predicted by the model for the coldest and hottest time of year? What is the meaning of the b parameter, and does its value make sense?

In [16]:
# your solution here