SPI Data Reading, Linear, and Non-linear Fit
Here is the program, extending last week's program, that reads in single photon interference data and performs a couple fits. The program reads in data and plots them media:spi_data.png. Then a linear fit is performed and a plot of the data plus this model is generated media:spi_data_linear_model.png. Finally a nonlinear fit is performed and a plot with the two models and the data is generated media:spi_data_linear_and_nonlinear_model.png. See comments in the
import numpy as np
import pylab as plt
#Read data
_volts, _A, _B, _C, _AB, _AC, _BC, _ABC = np.loadtxt('spi_indistinguishable.txt',\
unpack=True, skiprows=9, usecols=(2,3,4,5,6,7,8,9), \
delimiter=',')
print '_volts =', _volts
print '_ABC =', _ABC
print '_AC =', _AC
volts = np.arange(30,71)
print 'volts =', volts
N = len(volts)
AC = np.zeros(N)
AC_err = np.zeros(N)
nsamp = np.zeros(N)
print 'range(N) =', range(N)
for i in range(N):
#print "i =", i
#print "voltage =", volts[i]
indices = np.where(_volts == volts[i])
#print "indices =", indices
#print "_volts at indices =", _volts[indices]
#print "_AC at indices =", _AC[indices]
AC[i] = np.mean(_AC[indices])
AC_err[i] = np.std(_AC[indices])
nsamp[i] = len(_AC[indices])
plt.clf()
plt.errorbar(volts, AC, AC_err/np.sqrt(nsamp), fmt='r.', label='AC')
plt.xlabel('Piezoelectric Voltage (V)')
plt.ylabel('Count (per 0.01 s)')
plt.title('Indistiguishable Paths')
plt.legend(loc='lower left')
plt.savefig('spi_data.png')
#Linear fit y(V) = a*cos(2*pi*V/period) + b*sin(2*pi*V/period) + c
# = a*f1(V)+b*f2(V)+c*f3(V)
period = 9.
#Building projection matrix M
cos = np.cos(2*np.pi*volts/period)
sin = np.sin(2*np.pi*volts/period)
one = np.ones(len(volts))
MT = np.array([cos, sin, one])
M = np.transpose(MT)
print "Dim of M =", M.shape
#Building noise inverse matrix
Ninv = np.diag(1/(AC_err**2/nsamp))
print "Dim of Ninv =", Ninv.shape
#Compute the best fit parameters a = (alpha)^-1 beta
alpha = np.dot(MT,np.dot(Ninv,M)) #np.dot is matrix multiplication
print "Dim of alpha =", alpha.shape
beta = np.dot(MT,np.dot(Ninv,AC))
cov = np.linalg.inv(alpha)
params = np.dot(cov,beta)
print "Linear fit results "
print "a =", params[0], '+/-', np.sqrt(cov[0][0])
print "b =", params[1], '+/-', np.sqrt(cov[1][1])
print "c =", params[2], '+/-', np.sqrt(cov[2][2])
plt.plot(volts, params[0]*cos+params[1]*sin+params[2]*one, label='Linear Model')
plt.legend(loc='lower left')
plt.savefig('spi_data_and_linear_model.png')
#Try non-linear fit -- determine the frequency/period
# model: a*sin(2*pi*b*V+c)+d (4 parameters)
def y(V,a,b,c,d):
return a*np.sin(2*np.pi*b*V+c)+d
def dyda(V,a,b,c,d):
return np.sin(2*np.pi*b*V+c)
def dydb(V,a,b,c,d):
return a*2*np.pi*V*np.cos(2*np.pi*b*V+c)
def dydc(V,a,b,c,d):
return a*np.cos(2*np.pi*b*V+c)
def dydd(V,a,b,c,d):
return 1.
def alpha(lam,V,a,b,c,d,var):
alpha = np.zeros((4,4))
alpha[0][0] = np.sum(dyda(V,a,b,c,d)**2/var*(1+lam))
alpha[1][1] = np.sum(dydb(V,a,b,c,d)**2/var*(1+lam))
alpha[2][2] = np.sum(dydc(V,a,b,c,d)**2/var*(1+lam))
alpha[3][3] = np.sum(dydd(V,a,b,c,d)**2/var*(1+lam))
alpha[0,1] = np.sum(dyda(V,a,b,c,d)*dydb(V,a,b,c,d)/var)
alpha[0,2] = np.sum(dyda(V,a,b,c,d)*dydc(V,a,b,c,d)/var)
alpha[0,3] = np.sum(dyda(V,a,b,c,d)*dydd(V,a,b,c,d)/var)
alpha[1,2] = np.sum(dydb(V,a,b,c,d)*dydc(V,a,b,c,d)/var)
alpha[1,3] = np.sum(dydb(V,a,b,c,d)*dydd(V,a,b,c,d)/var)
alpha[2,3] = np.sum(dydc(V,a,b,c,d)*dydd(V,a,b,c,d)/var)
alpha[1,0] = alpha[0,1]
alpha[2,0] = alpha[0,2]
alpha[3,0] = alpha[0,3]
alpha[2,1] = alpha[1,2]
alpha[3,1] = alpha[1,3]
alpha[3,2] = alpha[2,3]
return alpha
def beta(V,a,b,c,d,data,var):
beta = np.zeros(4)
beta[0]=np.sum((data-y(V,a,b,c,d))*dyda(V,a,b,c,d)/var)
beta[1]=np.sum((data-y(V,a,b,c,d))*dydb(V,a,b,c,d)/var)
beta[2]=np.sum((data-y(V,a,b,c,d))*dydc(V,a,b,c,d)/var)
beta[3]=np.sum((data-y(V,a,b,c,d))*dydd(V,a,b,c,d)/var)
return beta
def chisq(V,data,var,a,b,c,d):
return np.sum((data-y(V,a,b,c,d))**2/var)
def marquardt(V, data, var, a0, b0, c0, d0, verbose=False):
a=a0
b=b0
c=c0
d=d0
lam = 0.001
i=0
while(True):
chisq0 = chisq(V,data,var,a,b,c,d)
da,db,dc,dd = np.dot(np.linalg.inv(alpha(lam,V,a,b,c,d,var)), \
beta(V,a,b,c,d,data,var))
chisq1 = chisq(V,data,var,a+da,b+db,c+dc,d+dd)
dchisq = chisq1-chisq0
if dchisq > 0: #Uphill, oh no!
lam *= 10
continue #go back to top of while loop
if dchisq < 0: #Downhill, life is good!
lam /= 10
a = a+da
b = b+db
c = c+dc
d = d+dd
i+=1
if verbose:
print i, chisq1
if abs(dchisq/chisq0) < 0.001:
break #leave the loop
return a, b, c, d, np.linalg.inv(alpha(lam,V,a,b,c,d,var))
print "Non-linear Marquardt Search (iteration, chi-sq)"
a,b,c,d,cov = marquardt(volts, AC, AC_err**2/nsamp, .5 , .1, 0.,4., verbose=True)
print "Non-linear fit results"
print "a =", a, '+/-', np.sqrt(cov[0][0])
print "b =", b, '+/-', np.sqrt(cov[1][1])
print "c =", c, '+/-', np.sqrt(cov[2][2])
print "d =", d, '+/-', np.sqrt(cov[3][3])
plt.plot(volts, a*np.sin(2*np.pi*b*volts+c)+d, label='Nonlinear Model')
plt.legend(loc='lower left')
plt.savefig("spi_data_linear_and_nonlinear_model.png")