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]:

    
data = np.load('decay_osc.npz')



In [3]:

    
tdata = data['tdata']
ydata = data['ydata']
dy = data['dy']



In [4]:

    
plt.figure(figsize=(8,6))
plt.plot(tdata, ydata)
plt.xlabel('Time');
plt.ylabel('Y Data');
plt.title('Oscillation Data');
plt.tick_params(axis='x',top='off',direction='out');
plt.tick_params(axis='y',right='off',direction='out');
plt.errorbar(tdata, ydata, dy,fmt='.k', ecolor='lightgray');



In [5]:

    
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 [6]:

    
def decay_model(t,A,Y,w,d):
    return A*np.exp(-Y*t)*np.cos(w*t)+d



In [7]:

    
theta_best, theta_cov = opt.curve_fit(decay_model, tdata, ydata,absolute_sigma=True)



In [8]:

    
A = theta_best[0]
Y = theta_best[1]
w = theta_best[2]
d = theta_best[3]



In [9]:

    
print('A = {0:.3f} +/- {1:.3f}'.format(A, np.sqrt(theta_cov[0,0])))
print('λ = {0:.3f} +/- {1:.3f}'.format(Y, np.sqrt(theta_cov[1,1])))
print('ω = {0:.3f} +/- {1:.3f}'.format(w, np.sqrt(theta_cov[2,2])))
print('δ = {0:.3f} +/- {1:.3f}'.format(d, np.sqrt(theta_cov[3,3])))









    



A = -5.285 +/- 0.520
λ = 0.071 +/- 0.014
ω = -1.012 +/- 0.009
δ = 0.229 +/- 0.147



In [12]:

    
tfit = np.linspace(0,20)



In [17]:

    
plt.figure(figsize=(8,6))
plt.xlim(0,20)
plt.ylim(-6,6)
plt.scatter(tdata, ydata,color = 'black')
plt.plot(tfit,decay_model(tfit,A,Y,w,d))
plt.xlabel('Time');
plt.ylabel('Y Data');
plt.title('Oscillation Curve Fit');
plt.tick_params(axis='x',top='off',direction='out');
plt.tick_params(axis='y',right='off',direction='out');



In [11]:

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