HW 4 Solution¶

In [1]:
#Importing numpy and matplotlib
import numpy as np
import matplotlib.pyplot as plt

#For Curve Fitting
from scipy.optimize import curve_fit  

#For reading urls
import urllib

plt.rc('figure', figsize = (6, 3.5)) # Reduces overall size of figures

Problem 1¶

In [2]:
link = 'https://www.eg.bucknell.edu/~phys310/hw/assignments/fitting_1/fit_1.dat'
fh = urllib.request.urlopen(link)
y1, u1 = np.loadtxt(fh, unpack=True)
x1 = np.arange(len(y1))

#We perform a weighted fit of the data
def cons(x,c):
    return 0*x+c

p1,pcov1 = curve_fit(cons,x1,y1, sigma=u1, absolute_sigma=True)

#Plot Data with Results
plt.figure(0)
plt.errorbar(x1,y1,yerr = u1, fmt = 'ko',label = 'data')
plt.plot(x1,cons(x1,*p1), label = 'fit')
plt.xlabel('trial');
plt.ylabel('y');
plt.legend()

#Print Value of variable
y1values,y1err = p1[0],np.sqrt(pcov1[0][0])
print('y = '+ str(np.round(y1values,3)) + '+/-' + str(np.round(y1err,3)))

#Determine Chi-Squared
chi2_1 = np.sum(((y1values- y1)/u1)**2)
print('The Chi-Squared value is '+ str(np.round(chi2_1,4)))
print('The Reduced Chi-Squared value is '+ str(np.round(chi2_1/(len(y1)--len(p1)),4)))

y = 14.469+/-0.059
The Chi-Squared value is 14.669
The Reduced Chi-Squared value is 0.6985
In [3]:
## Computing Weighted Sum

weights = 1/(y1err**2) #form a weighted average with the given data
yWeightedAvg = np.sum(weights*y1values)/np.sum(weights)
alpha_y = np.sqrt(1/np.sum(weights))
yWeightedAvg, alpha_y
Out[3]:
(14.468891058040507, 0.05916554751147506)

Note That this yields the same value as the weighted fit.

Problem 2¶

In [4]:
link = 'https://www.eg.bucknell.edu/~phys310/hw/assignments/fitting_2/fit_2.dat'
fh = urllib.request.urlopen(link)
x2,y2, u2 = np.loadtxt(fh, unpack=True)


#We perform a weighted fit of the data
def flin(x,m,b):
    return m*x+b


p2,pcov2 = curve_fit(flin,x2,y2, sigma=u2, absolute_sigma=True)

#Plot Data with Results

plt.figure(1)
plt.errorbar(x2,y2,yerr = u2, fmt = 'ro',label = 'data')
plt.plot(x2,flin(x2,*p2), label = 'fit')
plt.xlabel('x');
plt.ylabel('y');
plt.legend()

#Print Value of variable
m,m_err = p2[0],np.sqrt(pcov2[0][0])
b,b_err = p2[1],np.sqrt(pcov2[1][1])

print('slope = '+ str(np.round(m,3)) + '+/-' + str(np.round(m_err,3)))
print('intercept = '+ str(np.round(b,3)) + '+/-' + str(np.round(b_err,3)))

slope = 2.626+/-0.155
intercept = 1.799+/-1.858
In [5]:
plt.figure(2)
res2 = (flin(x2,*p2)-y2)/u2
plt.errorbar(x2,res2,np.ones(len(x2)),fmt = 'ro')
#plt.errorbar(x2,u2*res2,u2,fmt = 'ko')
plt.axhline(y = 1,linestyle = 'dashed',label = '+/- 1 sigma')
plt.axhline(y = -1,linestyle = 'dashed')
plt.axhline(y = 0)
plt.ylabel('Norm Residuals')
plt.xlabel('x')
plt.legend()
Out[5]:
<matplotlib.legend.Legend at 0x28cedb30a10>
In [6]:
#Counting the number of points within +/- 1 standard error

n1sig = np.sum(np.abs(res2)<1)

print(str(n1sig)+' points are withing one standard deviation')

15 points are withing one standard deviation
In [7]:
##Compute Chi-Squared
#Determine Chi-Squared
chi2_2 = np.sum(res2**2)
print('The Chi-Squared value is '+ str(np.round(chi2_2,4)))
print('The Reduced Chi-Squared value is '+ str(np.round(chi2_2/(len(y2)-len(p2)),4)))

The Chi-Squared value is 16.6917
The Reduced Chi-Squared value is 0.9273
In [ ]:

Problem 3¶

In [8]:
link = 'https://www.eg.bucknell.edu/~phys310/hw/assignments/fitting_3/fit_3.dat'
fh = urllib.request.urlopen(link)
x3,y3, u3 = np.loadtxt(fh, unpack=True)


#We perform a weighted fit of the data
def fSin2(x,a1,a2):
    return a1*np.sin(2*np.pi*x)+a2*np.sin(4*np.pi*x)


p3,pcov3 = curve_fit(fSin2,x3,y3, sigma=u3, absolute_sigma=True)

#Plot Data with Results

plt.figure(1)
plt.errorbar(x3,y3,yerr = u3, fmt = 'mo',label = 'data')
plt.plot(x3,fSin2(x3,*p3), label = 'fit')
plt.xlabel('x');
plt.ylabel('y');
plt.legend()
plt.title('$y = a_1 sin(2\pi x)+ a_2 sin(4\pi x)$')


#Print Value of variable
a1,a1_err = p3[0],np.sqrt(pcov3[0][0])
a2,a2_err = p3[1],np.sqrt(pcov3[1][1])

print('a1 = '+ str(np.round(a1,3)) + '+/-' + str(np.round(a1_err,3)))
print('a2 = '+ str(np.round(a2,3)) + '+/-' + str(np.round(a2_err,3)))

a1 = 1.965+/-0.08
a2 = 0.879+/-0.08

Note that this is a linear fit.

In [9]:
plt.figure(4)
res3 = (fSin2(x3,*p3)-y3)/u3
plt.errorbar(x3,res3,np.ones(len(x3)),fmt = 'mo')
#plt.errorbar(x2,u2*res2,u2,fmt = 'ko')
plt.axhline(y = 1,linestyle = 'dashed',label = '+/- 1 sigma')
plt.axhline(y = -1,linestyle = 'dashed')
plt.axhline(y = 0)
plt.ylabel('Norm Residuals')
plt.xlabel('x')
plt.legend()
Out[9]:
<matplotlib.legend.Legend at 0x28cedbd0a10>
In [10]:
##Compute Chi-Squared
#Determine Chi-Squared
chi2_3 = np.sum(res3**2)
print('The Chi-Squared value is '+ str(np.round(chi2_3,4)))
print('The Reduced Chi-Squared value is '+ str(np.round(chi2_3/(len(y3)-len(p3)),4)))

The Chi-Squared value is 36.6146
The Reduced Chi-Squared value is 0.7472
In [ ]:
In [ ]: