In [1]:
import numpy as np
from scipy import optimize
from scipy import stats
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
In [2]:
# Following is an Ipython magic command that puts figures in notebook.
%matplotlib notebook
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot params for new size
Define functions¶
In [3]:
def f(x,m,b):
'''Simple linear function with slope m and intercept b'''
return x*m + b
def chi2(x, y, u, m, b):
'''Chisquare as a function of data (x, y, and yerr=u), and model
parameters slope and intercept (m and b)'''
return np.sum((y - f(x, m, b))**2/u**2)
Linear fit to data for $m$ and $b$¶
Data to be fit:¶
In [4]:
# Or: data = np.loadtxt("file.dat")
# Format: [[x1,y1,u1], [x2,y2,u2], ... ] where u1 is uncertainty in y1
data = np.array([[1, 2.947032612427293, 0.5],
[2, 6.168779380682309, 0.5],
[3, 7.1618838821688, 0.5],
[4, 9.590549514954866, 0.5],
[5, 11.20657, 0.5]])
x, y, u = data.T
In [5]:
xc = np.linspace(0,6,201) # quasi-continuous set of x's for function plot
plt.figure()
plt.title("data",fontsize=14)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.axhline(0,color='magenta')
plt.xlim(0,6)
plt.errorbar(x,y,yerr=u,fmt='o');
Perform fit¶
In [6]:
popt, pcov = optimize.curve_fit(f, x, y, sigma=u, absolute_sigma=True)
slope = popt[0]
alpha_m = np.sqrt(pcov[0,0]) # Std error in slope
intercept = popt[1]
alpha_b = np.sqrt(pcov[1,1]) # Std error in intercept
print("slope =", slope,"+/-", alpha_m,"\n")
print("intercept =", intercept,"+/-", alpha_b,"\n")
print("covariance matrix =","\n",pcov,"\n")
pcov_data = pcov
print("chi2 =", chi2(x, y, u,*popt))
# print("reduced chi2 = chi2/(len(x)-len(popt)) =", chi2(x, y, u, *popt)/(len(x)-len(popt)))
a = chi2(x,y,u,*popt)
In [7]:
xc = np.linspace(0,6,201) # quasi-continuous set of x's function plot
plt.figure()
plt.title("data with best fit line",fontsize=14)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.axhline(0, color='magenta')
plt.xlim(0,6) # Pad x-range on plot
plt.errorbar(x, y, yerr=u, fmt='o');
plt.plot(xc ,f(xc, slope, intercept));
Residuals:¶
In [8]:
plt.figure()
plt.axhline(0,color='magenta')
plt.title('normalized residuals')
plt.xlabel('$x$')
plt.ylabel('normalized residuals')
plt.grid(True)
plt.errorbar(x,(f(x,slope,intercept)-y)/u,1,fmt='o')
plt.xlim(0,6);
What confidence can we have in the model?¶
If normalized residuals look ok:
- Assume the model, along with the best fit parameters are ok ("null hypothesis")
- Generate simulated data based on the model, best-fit parameters, and uncertainties. (Use same $x$ values).
- Fit each of the simulated data sets
- What's the probability of getting values of $\chi^2$ from the simulated data sets that are greater than the value given by the data that was actually measured. In other words, is our value of $\chi^2$ a reasonable result given the fluctuations we expect in our data?
In [9]:
n_sim = 10000 # Number of data sets to simulate
count = 0
for i in range(n_sim):
mSim = np.array([]) # Array for simulated slopes
bSim = np.array([]) # Array for simulated intercepts
chi_sim = np.array([]) # Array for chisq from simulated sets
for i in range(n_sim):
ySim = stats.norm.rvs(slope*x+intercept,u) # Generate simulated data set
popt, pcov = optimize.curve_fit(f, x, ySim, sigma=u, absolute_sigma=True) # Fit simulated data
bSim = np.append(bSim, popt[1]) # Record intercept
mSim = np.append(mSim, popt[0]) # Record slope
aa = chi2(x,ySim,u,*popt)
chi_sim = np.append(chi_sim, aa)
if aa > a:
count += 1
print(count/n_sim)
In [10]:
1 - stats.chi2.cdf(3.700621, df=3)
Out[10]:
Make "data" for contour plot¶
- Choose ranges of $m$ and $b$ for contour plot
- Make meshgrid of slope and intercept values
- Calculate values of $\chi^2$ at grid points
In [11]:
m = np.linspace(1.5, 2.5, 201)
b = np.linspace(0.5, 2.5, 201)
M, B = np.meshgrid(m, b, indexing='ij')
Z = np.zeros((len(m),len(b)))
for i in range(len(m)):
for j in range(len(b)):
Z[i,j] = chi2(x, y, u, m[i],b[j]) - chi2(x, y, u, *popt)
In [12]:
plt.figure()
CS = plt.contour(M, B, Z, levels=[1,2.7,4,9]) # contour plots of chi2 = 1, 2.7, 4, 9 above chi2_min
#
#SC = plt.scatter(mSim, bSim, s=0.01) # overlay scatterplot of simulted data
#
plt.xlabel('$m$')
plt.ylabel('$b$')
plt.grid()
plt.axhline(intercept)
#plt.axhline(intercept + alpha_b, linestyle='--')
#plt.axhline(intercept - alpha_b, linestyle='--')
plt.axvline(slope)
#plt.axvline(slope + alpha_m,linestyle='--')
#plt.axvline(slope - alpha_m,linestyle='--')
plt.clabel(CS, inline=1, fontsize=10);
Correlation between fit parameters¶
Note that the contours are tilted and squashed; this indicates that there is correlation between the two fit parameters, mainly caused because the y-intercept is not within the data range. These elongated contours are telling you that a small change in $m$ can, to some extent, be compensated for by a corresponding change in $b$.
Fits with correlated data aren't great, because they tend to lead to an overestimation of the uncertainties in the fit parameters. In this case, note tha thte uncertainty in $b$ is $\sim 0.5$ due to this correlation effect.
We can do better if we perform our linear fit on a dataset whose x-values are centered on $x=0$.
In [13]:
x_mean = np.mean(x) # determine the center of the distribution of x values
x_shift = x - x_mean # shift the array so that the x data are centered on x = 0
plt.figure()
plt.title("shifted data",fontsize=14)
plt.xlabel('$x_{shift}$')
plt.ylabel('$y$')
plt.axhline(0,color='magenta')
plt.xlim(0-x_mean,6-x_mean)
plt.errorbar(x_shift,y,yerr=u,fmt='o');
In [14]:
popt, pcov = optimize.curve_fit(f, x_shift, y, sigma=u, absolute_sigma=True)
slope = popt[0]
alpha_m = np.sqrt(pcov[0,0]) # Std error in slope
intercept = popt[1]
alpha_b = np.sqrt(pcov[1,1]) # Std error in intercept
print("slope =", slope,"+/-", alpha_m,"\n")
print("intercept =", intercept,"+/-", alpha_b,"\n")
print("covariance matrix =","\n",pcov,"\n")
pcov_data = pcov
print("chi2 =", chi2(x_shift, y, u,*popt))
# print("reduced chi2 = chi2/(len(x)-len(popt)) =", chi2(x, y, u, *popt)/(len(x)-len(popt)))
a = chi2(x_shift,y,u,*popt)
In [15]:
m = np.linspace(1.5, 2.5, 201)
b = np.linspace(6.5, 8.5, 201)
M, B = np.meshgrid(m, b, indexing='ij')
Z = np.zeros((len(m),len(b)))
for i in range(len(m)):
for j in range(len(b)):
Z[i,j] = chi2(x_shift, y, u, m[i],b[j]) - chi2(x_shift, y, u, *popt)
In [16]:
plt.figure()
CS = plt.contour(M, B, Z, levels=[1,2.7,4,9]) # contour plots of chi2 = 1, 2.7, 4, 9 above chi2_min
#
#SC = plt.scatter(mSim, bSim, s=0.01) # overlay scatterplot of simulted data
#
plt.xlabel('$m$')
plt.ylabel('$b$')
plt.grid()
plt.axhline(intercept)
#plt.axhline(intercept + alpha_b, linestyle='--')
#plt.axhline(intercept - alpha_b, linestyle='--')
plt.axvline(slope)
#plt.axvline(slope + alpha_m,linestyle='--')
#plt.axvline(slope - alpha_m,linestyle='--')
plt.clabel(CS, inline=1, fontsize=10);
Version information¶
%version_information is an IPython magic extension for showing version information for dependency modules in a notebook;
See
https://github.com/jrjohansson/version_information%version_informationis available on Bucknell computers on the linux network. You can easily install it on any computer.
In [17]:
%load_ext version_information
In [18]:
%version_information numpy, scipy, matplotlib
Out[18]:
In [ ]: