PHYS 310: HW 4 Solutions

There are three problems in this set. You must download the data for these problems from the website http://www.eg.bucknell.edu/physics/ph310/fit1.html or Moodle.

In [1]:
import scipy as sp
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
%matplotlib inline
In [2]:
# Set figure defaults
plt.rc('figure', figsize = (6, 4.5))            
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('xtick', labelsize=14)
plt.rc('ytick', labelsize=14)

Start by downloading the data files into the working directory The loadtxt function loads the content of the data file into a Matrix. The unpack = 'True' option transposes the matrix so that each column is in a separate array

In [3]:
data1 = sp.loadtxt('hw4-1.dat', unpack = 'True')
data2 = sp.loadtxt('hw4-2.dat', unpack = 'True')
data3 = sp.loadtxt('hw4-3.dat', unpack = 'True')
In [4]:
#Uncomment to look at the data.
#data1
In [5]:
### Linear Model Fit FUnctions
## No need to edit this.
## The particular model is modified by editing the basis functions.

def func(x, a):
    '''Given basis functions and coefficients, returns value of linear function'''
    return sp.dot(basis(x).T, a)


def LinearModelFit(x, y, u):
    '''Performs linear least squares fit to a set of 2-d data with uncertainties
    
    x = array of x values [x0, x1, x2, ...]
    y = array of values of dependent variable [y0, y1, y2, ...]
    u = array of uncertainties for dependent variable [u0, u1, u2, ...]
    '''
    X = basis(x).T    # Basis functions evaluated at all x (the X_j(x_i)) of N.R.)
    W = sp.diag(1/u)  # Matrix with uncertainties on diagonal
    Xw = sp.dot(W,X)  # A_ij of Eq. (14.3.4)
    Yw = sp.dot(y,W)  # b_i of Eq. (14.3.5)
    fit = sp.linalg.lstsq(Xw,Yw)  # lstq returns: best values, chi2, ....
    covariance = sp.linalg.inv(sp.dot(Xw.T,Xw))
    uncertainty = sp.sqrt(sp.diag(covariance))
    return(fit[0], uncertainty,fit[1], covariance)

Problem 1

Use the “hw4-1.dat” file for this problem. These data come from twenty (20) experiments nominally measuring the same physical quantity. The data for each measurement (point) is on a single line in the file. The first number on each line is the value of the measured quantity, and the second is the uncertainty in the measurement.

(a) What value do you quote for this quantity based on the data? (Include an uncertainty.)

In [6]:
y,u = data1;
w = 1/u**2;
mean = sp.sum(y*w)/sp.sum(w);
wunc = 1./sp.sqrt(sp.sum(w))

#We can round the numbers to three significant digits using the sp.around() function.
# The str command casts the number to a string.
print('Weighted mean = ' + str(sp.around(mean,3))+ ' +/- '+ str(sp.around(wunc,3)))
Weighted mean = 14.469 +/- 0.059

(b) Each data point has its own uncertainty $\sigma_i$. How many of the data points lie within 1 $\sigma$ of the mean value you determined?

(c) What is the goodness-of-fit parameter $\chi^2$ for this data?

In [7]:
x = sp.arange(len(y))
ymean = sp.zeros(len(y))+mean;

plt.errorbar(x,y,u, fmt = 'ko')
plt.plot(x,ymean,'r-')
plt.xlabel('experiment number');

plt.ylabel('Result')
Out[7]:
<matplotlib.text.Text at 0xa3fe860>

In order to compute the number within one standard deviation of the weighted, we can calculate the normalized residuals (z) and count the number with magnitude < 1

The gooddess of fit chi squared is simply the sum of the square of residuals

In [8]:
z = (y - ymean)/u
n = sum(abs(z) < 1)

print(str(n) + ' data points are within 1 sigma of the weighted mean.' )

chisq = sum(z**2)

print('chisq =' + ' '+str(sp.around(chisq,1)))
15 data points are within 1 sigma of the weighted mean.
chisq = 14.7

Problem 2

Data for an experiment in which there is a suspected linear relationship between measured values of x and y. The data for each point is on a single line in the file. The first number on each line is the value of x, the second is the value of y, and the third is the uncertainty in y. Uncertainties in x are assumed to be negligible.

(a.) Using this data, what value do you quote for the slope and intercept? (Include uncertainties.)

(b.) Plot your residuals. How many of the data points lie within 1σ of the line you determined?

(c.) Use Eq. (5.9) from Hughes & Hase and calculate the goodness-of-fit parameter chi squared for this data?

In [9]:
x2,y2,u2 = data2

def basis(x):
    X1 = x
    X0 = 0.*X1 + 1. # Need array of len(x), thus the 0.*X1
    return sp.array([X0,X1])

f0,unc,f1,covm = LinearModelFit(x2,y2,u2)
b,a = f0;
siga,sigb = unc; 

print ('Linear Model: y = ax + b')
print(' a = ' + str(sp.around(a,2)) + ' +/-' + str(sp.around(siga,2)))
print(' b = ' + str(sp.around(b,2)) + ' +/-' + str(sp.around(sigb,2)))

yfit = f0[1]*x2 + f0[0];
plt.errorbar(x2,y2,u2, fmt = 'ko')
plt.plot(x2,yfit, 'r-');
plt.xlabel('x'); 
plt.ylabel('y')
plt.title('data2: data and fit')
Linear Model: y = ax + b
 a = 2.63 +/-1.86
 b = 1.8 +/-0.16
Out[9]:
<matplotlib.text.Text at 0xba8de80>
In [10]:
# Plot(normalized residuals)
z2 = (y2 - yfit)/u2
plt.errorbar(x2,z2,1, fmt = 'ko')
plt.xlabel('x')
plt.ylabel('Normalized Residuals')
plt.title('data2 residuals')

plt.plot((0,21), (0,0), 'r--')
plt.xlim(0, 21)


n2 = sum(abs(z2) < 1)
print(str(n2) + ' data points are within 1 sigma of the weighted mean.' )


chisq = sum(z2**2)
print('chisq =' + ' '+str(sp.around(chisq,1)))
print('This looks like a good fit')
15 data points are within 1 sigma of the weighted mean.
chisq = 16.7
This looks like a good fit

Problem 3

For the data in this set, there is reason to believe that that x and y are related by a function of the form y=a1sin(2πx)+a2sin(4πx). The data for each point is on a single line in the file. The first number on each line is the value of x, the second is the value of y, and the third is the uncertainty in y. Uncertainties in x are assumed to be negligible.

(a.) Fit the data to the assumed functional form to find values for a1 and a2.

(b.) Plot your residuals.

(c.) Use Eq. (5.9) from Hughes & Hase and calculate the goodness-of-fit parameter χ2 for this data?

(d.) Is your fit good?

(e.) Give your values for a1 and a2 (including uncertainties).

In [11]:
x3,y3,u3 = data3

## A linear function used to model the data
def basis(x):
    X1 = sp.sin(4*sp.pi*x)
    X0 = sp.sin(2*sp.pi*x)
    return sp.array([X0,X1])

f0,unc,f1,cov3 = LinearModelFit(x3,y3,u3)

a1,a2 = f0;
sig1,sig2 = unc; 

yfit3 =  a1*sp.sin(2*sp.pi*x3) + a2*sp.sin(4*sp.pi*x3)

print ('y = a1 sin(2 pi x) + a2 sin(2 pix)')
print(' a1 = ' + str(sp.around(a1,2)) + ' +/-' + str(sp.around(sig1,2)))
print(' a2 = ' + str(sp.around(a2,2)) + ' +/-' + str(sp.around(sig2,2)))


#yfit3 = SinModel(x3, *vals3)
plt.errorbar(x3,y3,u3, fmt = 'ko')
plt.plot(x3,yfit3, 'r-');
plt.xlabel('x'); 
plt.ylabel('y')
plt.title('data3: data and fit')
y = a1 sin(2 pi x) + a2 sin(2 pix)
 a1 = 1.97 +/-0.08
 a2 = 0.88 +/-0.08
Out[11]:
<matplotlib.text.Text at 0xbd5f5f8>
In [12]:
# Plot(residuals)
res3 = y3 - yfit3;
z3 = (y3 - yfit3)/u3
plt.errorbar(x3,res3,u3, fmt = 'ko')
plt.xlabel('x')
plt.ylabel('residuals = y - model')
plt.title('data3 residuals')

plt.plot((-0.1,1.1), (0,0), 'r-')
plt.xlim(-0.1, 1.1)

n3 = sum(abs(z3) < 1)
print(str(n3) + ' data points are within 1 sigma of the weighted mean.' )


chisq = sum(z3**2)
print('chisq =' + ' '+str(sp.around(chisq,1)))
print('This looks like a good fit')
38 data points are within 1 sigma of the weighted mean.
chisq = 36.6
This looks like a good fit