Example: Generating data and fitting a curve

In [5]:
#  Our first step is to import any required modules.
#  I will use a shortcut to import commonly-used modules in scientific computing:
%pylab inline   
#  A 'cleaner' way would be to import only those modules that we really need.
Populating the interactive namespace from numpy and matplotlib

\(1\). Generate a quadratic with random coefficients and plot.

In [4]:
coef = np.random.random(3)
print coef
x = np.linspace(-5, 5, 100)
y0 = coef[0]*x**2 + coef[1]*x + coef[2]
plt.plot(x, y0);
[ 0.94120775  0.31219185  0.87239774]

\(2\). Add in some Gaussian noise, representing random error in a measurement.

In [4]:
sigma = 1.0
y = y0 + sigma*np.random.normal(size=len(x))
plt.scatter(x, y);

\(3\). Save the data to a text file. Here we will use NumPy instead of the core library.

In [5]:
fname = "quadratic_with_noise.dat"
np.savetxt(fname, [x, y])

\(4\). Read the data from a text file (just to show we can), and plot once more.

In [6]:
x,y = np.loadtxt(fname)
plt.plot(x, y, 's');

\(5\). Now attempt to fit the data set using "curve_fit" from SciPy.

In [7]:
from scipy.optimize import curve_fit
def func(x, a, b, c):
    return a*x**2 + b*x + c
co, pcov = curve_fit(func, x, y)

\(6\). Check that the fitted coefficients are reasonable estimates.

In [9]:
for i in range(3):
    print "{}\t{}\n".format(coef[i], co[i])

ytrue = coef[0]*x**2 + coef[1]*x + coef[2]
yest = co[0]*x**2 + co[1]*x + co[2]
plt.plot(x, y, 's')
plt.plot(x, ytrue, '-')
plt.plot(x, yest, '+');
0.968082653695	0.974474964906

0.247884464147	0.313064145505

0.97281445601	1.05805371579