In [34]:
# In this script we will fit data to a model that is linear in the fit parameters.
# The fit will be weighted by the measurement uncertainties.
import numpy as np
In [35]:
# Enter the data.
t = np.arange(2, 62, 4)
T = np.array([40.5, 37.8, 36.9, 34.7, 33.0, 32.4, 30.5, 29.3, 28.6, 27.0, 26.6, 25.7, 25.0, 24.1, 23.6])
errT = np.array([0.8, 0.8, 0.8, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.3, 0.3])
In [36]:
# Plot the data using errorbar(x,y,e).
import matplotlib.pyplot as plt
plt.errorbar(t, T, errT, fmt = 'ko', markersize = 5,\
                 linewidth = 1.8,\
                 markeredgecolor = 'b',\
                 markerfacecolor = (.49, 1, .63),\
                 capsize = 5)
plt.xlabel('time (hours)')
plt.ylabel('temperature (Celsius)');
In [37]:
# To do the actual fit, we will use the 'curve_fit()' function from, the
# 'SciPy' module.  This way of fitting is very nice because we
# use it for all types of fit models (linear, polynomial, line-in-parameter fits,
# and nonlinear fit).  It is capable of doing both unweighted
# and weighted fits and it will return uncertainties in the fit parameters.

# The first step is to define a function for the model that we will fit our
# data to.  In this case, the modeil is a quadratic.
def LinParamFunc(x, Trm, T0):
    y = Trm*(1 - np.exp(-x/50)) + T0*np.exp(-x/50)
    return y
In [38]:
# Here is the actual command to execute the fit.  At a minimum, 'curve_fit()' 
# requires as inputs the function that defines the model, the x-data, and 
# the y-data.  The statement below tells 'curve_fit()' to return a list of the 
# the best-fit parameters and the covariance matrix which will be used to 
# determine the error in the fit parameters.
from scipy.optimize import curve_fit
a_fit, cov = curve_fit(LinParamFunc, t, T)
In [39]:
# Extract the best-fit parameters.
print('The best-fit parameters are:\n Trm =', a_fit[0], '\n',\
      'T0 =', a_fit[1])

The best-fit parameters are:
 Trm = 15.358535838161755 
 T0 = 41.16050335848032
In [40]:
# The uncertainties of the best-fit parameters are determined from the 
# square roots of the diagonal elements of the covariance matrix.  We
# can select the diagonal elements using (Note the use of the unicode
# character \u0394 to print $\Delta$):
print('The errors in the parameters are:\n \u0394Trm =', np.sqrt(np.diag(cov))[0],\
      '\n', '\u0394T0 =', np.sqrt(np.diag(cov))[1])

The errors in the parameters are:
 ΔTrm = 0.24119673138027387 
 ΔT0 = 0.181411792462186
In [41]:
# Next, plot the data and best-fit curve on the same graph.

# First, reproduce the data plot that was generated above.
plt.errorbar(t, T, errT, fmt = 'ko', markersize = 5,\
                 linewidth = 1.8,\
                 markeredgecolor = 'b',\
                 markerfacecolor = (.49, 1, .63),\
                 capsize = 5)
plt.xlabel('time (hours)')
plt.ylabel('temperature (Celsius)')

# Next, use the function that we defined to generate the best-fit curve.
xx = np.arange(0, 70, 0.1)
plt.plot(xx, LinParamFunc(xx, a_fit[0], a_fit[1]), 'k-')
plt.axis((0, 65, 23, 42));
In [42]:
# All of this has produced an "unweighted" fit to the data.  To include 
# weights, all we need to do is include another option in 'curve_fit()'.
# Everything else is exactly the same!  The new option is 'sigma' and it is
# simply a list of the errors in the y-values.  Note that many fitting
# routines require you to provide the actual weights as 1/sigma^2.  That is
# not the case here.  You just have to provide the absolute y-uncertainties.
# Here's a block of code that does the fit, extracts the best fit function, the
# best-fit parameters, the uncertainties, and then plots the result.
a_fit, cov = curve_fit(LinParamFunc, t, T, sigma = errT)

# Extract the fit parameters and their uncertainties.
print('The best-fit parameters are:\
      \n Trm \u00B1 \u0394Trm =', a_fit[0], '\u00B1', np.sqrt(np.diag(cov))[0],\
     '\n T0 \u00B1 \u0394T0 =', a_fit[1], '\u00B1', np.sqrt(np.diag(cov))[1])

# Plot the data.
plt.errorbar(t, T, errT, fmt = 'ko', markersize = 5,\
                 linewidth = 1.8,\
                 markeredgecolor = 'b',\
                 markerfacecolor = 'r',\
                 capsize = 5)
plt.xlabel('time (hours)')
plt.ylabel('temperature (Celsius)')

# Plot the best-fit line.
xx = np.arange(0, 70, 0.1)
plt.plot(xx, LinParamFunc(xx, a_fit[0], a_fit[1]), 'k-')
plt.axis((0, 65, 23, 42));

The best-fit parameters are:      
 Trm ± ΔTrm = 15.515080640989071 ± 0.18888502589858366 
 T0 ± ΔT0 = 41.004498558138955 ± 0.21243658642593372