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

    
data=np.load('decay_osc.npz')
tdata=data['tdata']
ydata=data['ydata']
dy=data['dy']



In [32]:

    
tdata,ydata,dy









    Out[32]:





(array([  0.23181023,   0.75756308,   0.79449188,   1.06440841,
          1.15442   ,   1.2074788 ,   1.6082512 ,   1.84020235,
          2.26847064,   2.4663936 ,   2.73670264,   3.51277747,
          3.69390206,   4.40311262,   4.92725839,   5.01612571,
          5.45756459,   5.77519885,   5.92089063,   6.34168734,
          7.25390092,   7.53691783,   7.5695311 ,   7.90738928,
          9.14672657,   9.91487324,  10.00612156,  10.90520466,
         11.11871571,  11.2010993 ,  11.54757171,  11.73261207,
         12.33314215,  12.53762032,  12.95405221,  13.2536048 ,
         13.26900451,  14.07723698,  14.38998527,  14.56464101,
         14.97249975,  16.09455364,  16.86871453,  17.01532558,
         18.03660694,  18.76147469,  19.04467933,  19.55174508,
         19.59947156,  19.69287049]),
 array([-5.11971604, -2.84019899, -4.55222367, -1.06305079, -1.09953019,
        -1.92218542, -0.59020311,  1.85089586,  3.16111471,  3.06495644,
         4.48901465,  4.51430527,  3.3223624 ,  3.060024  , -1.29224873,
        -0.3399654 , -3.40969819, -2.56486305, -2.37543706, -2.48478812,
        -2.81948031, -0.70505326,  1.34483188,  0.40966188,  3.42459268,
         0.69676068,  2.01185751,  0.16295044, -0.76996363, -1.27652668,
        -0.53415007, -1.07417858, -2.1761134 , -2.2696592 , -1.75609407,
        -1.85163761, -1.34630634,  0.70453535,  0.52977223,  1.59203359,
         2.77592308,  1.01095984,  0.46912102,  1.23047778, -1.16710647,
        -0.445895  , -2.57912718,  0.42990156, -0.561772  , -0.04998557]),
 array([-0.36496298,  0.5273403 , -1.31654655,  1.11718996,  0.7022342 ,
        -0.3474073 , -0.74961891,  0.74378254,  0.60086554,  0.01513279,
         0.99358238,  1.23501009,  0.38040594,  2.08020409, -0.64090611,
         0.56560325, -1.44522866, -0.11280122,  0.21087276,  0.16255611,
        -1.45250123,  0.02860089,  2.00304072,  0.28861505,  1.49829406,
        -0.93999144,  0.47556696,  0.01131728, -0.56794271, -0.9436255 ,
         0.29211413, -0.03459027, -0.75892383, -0.84310407, -0.48875094,
        -0.82467465, -0.33400402,  0.77781567,  0.23315726,  1.1089334 ,
         1.94636771,  0.08477602,  0.10017441,  0.99299024, -0.60108347,
         0.31703138, -1.84874594,  0.97018015, -0.04638158,  0.4140329 ]))



In [34]:

    
plt.plot(tdata,ydata)









    Out[34]:





[<matplotlib.lines.Line2D at 0x7ffba368f6a0>]



In [35]:

    
plt.errorbar?



In [36]:

    
plt.errorbar(tdata,ydata,dy,fmt='k.')









    Out[36]:





<Container object of 3 artists>



In [ ]:

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

    
def model(t,A,o,l,d):
    return A*np.exp(-l*t)*np.cos(o*t)+d



In [47]:

    
theta_best,theta_cov=opt.curve_fit(model,tdata,ydata,np.array((6,1,1,0)),dy,absolute_sigma=True)
print('A = {0:.3f} +/- {1:.3f}'.format(theta_best[0], np.sqrt(theta_cov[0,0])))
print('omega = {0:.3f} +/- {1:.3f}'.format(theta_best[1], np.sqrt(theta_cov[1,1])))
print('lambda = {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
omega = 1.001 +/- 0.001
lambda = 0.094 +/- 0.003
delta = 0.027 +/- 0.014



In [53]:

    
tfit=np.linspace(0,20,100)
A,o,l,d=theta_best
yfit=A*np.exp(-l*tfit)*np.cos(o*tfit)+d
plt.plot(tfit,yfit)
plt.plot(tdata,ydata,'k.')
plt.xlabel('time')
plt.ylabel('y')
plt.title('Decaying Oscillator')
plt.axhline(0,color='lightgray')









    Out[53]:





<matplotlib.lines.Line2D at 0x7ffba341cda0>



In [ ]:

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