Fitting Models Exercise 2

Imports



In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt

Fitting a decaying oscillation

For this problem you are given a raw dataset in the file decay_osc.npz. This file contains three arrays:

tdata: an array of time values
ydata: an array of y values
dy: the absolute uncertainties (standard deviations) in y

Your job is to fit the following model to this data:

$$ y(t) = A e^{-\lambda t} \cos{\omega t + \delta} $$

First, import the data using NumPy and make an appropriately styled error bar plot of the raw data.



In [2]:

    
f = open('decay_osc.npz','r')
r = np.load('decay_osc.npz')
tdata = r['tdata']
ydata = r['ydata']
dy = r['dy']
f.close()



In [3]:

    
plt.figure(figsize=(7,5))
plt.errorbar(tdata, ydata, dy, fmt='.b', ecolor='gray')
plt.tick_params(right=False,top=False,direction='out')
plt.xlabel('t')
plt.ylabel('y');



In [4]:

    
assert True # leave this to grade the data import and raw data plot

Now, using curve_fit to fit this model and determine the estimates and uncertainties for the parameters:

Print the parameters estimates and uncertainties.
Plot the raw and best fit model.
You will likely have to pass an initial guess to curve_fit to get a good fit.
Treat the uncertainties in $y$ as absolute errors by passing absolute_sigma=True.



In [5]:

    
def model(t,A,lamb,omega,delta):
    return A*np.exp(-lamb*t)*np.cos(omega*t)+delta



In [6]:

    
theta_best, theta_cov = opt.curve_fit(model, tdata, ydata, sigma=dy, absolute_sigma=True)



In [7]:

    
print('A = {0:.3f} +/- {1:.3f}'.format(theta_best[0], np.sqrt(theta_cov[0,0])))
print('lamb = {0:.3f} +/- {1:.3f}'.format(theta_best[1], np.sqrt(theta_cov[1,1])))
print('omega = {0:.3f} +/- {1:.3f}'.format(theta_best[2], np.sqrt(theta_cov[2,2])))
print('delta = {0:.3f} +/- {1:.3f}'.format(theta_best[3], np.sqrt(theta_cov[3,3])))









    



A = -4.896 +/- 0.063
lamb = 0.094 +/- 0.003
omega = -1.001 +/- 0.001
delta = 0.027 +/- 0.014



In [8]:

    
A_fit = theta_best[0]
lamb_fit = theta_best[1]
omega_fit = theta_best[2]
delta_fit = theta_best[3]



In [9]:

    
t_fit = np.linspace(0,20,50)
y_fit = A_fit*np.exp(-lamb_fit*t_fit)*np.cos(omega_fit*t_fit)+delta_fit

plt.figure(figsize=(7,5))
plt.errorbar(tdata, ydata, dy, fmt='.b', ecolor='gray')
plt.plot(t_fit,y_fit,color='r')
plt.tick_params(right=False,top=False,direction='out')
plt.xlabel('t')
plt.title('Best Fit Curve')
plt.ylabel('y');



In [ ]:

    
assert True # leave this cell for grading the fit; should include a plot and printout of the parameters+errors