SPI Data Reading, Linear, and Non-linear Fit
Jump to navigation
Jump to search
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. The text output of the program should look like this:
_volts = [ 30. 30. 30. ..., 70. 70. 70.] _ABC = [ 0. 0. 0. ..., 0. 0. 0.] _AC = [ 3. 4. 7. ..., 6. 3. 2.] volts = [30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70] range(N) = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40] Dim of M = (41, 3) Dim of Ninv = (41, 41) Dim of alpha = (3, 3) Linear fit results a = 0.374361111765 +/- 0.015423507935 b = 0.805224133015 +/- 0.0153865196485 c = 3.82291964389 +/- 0.0110404319537 Non-linear Marquardt Search (iteration, chi-sq) 1 4155.25421469 2 3428.96612097 3 3337.99353442 4 3231.38829482 5 2556.16118098 6 837.69911394 7 171.075641211 8 118.588679965 9 117.272843887 10 117.253784918 Non-linear fit results a = -0.889535513171 +/- 0.0154429366169 b = 0.112274935769 +/- 0.000238019986401 c = -3.0535050281 +/- 0.0751514516099 d = 3.824237946 +/- 0.0110433609376
Note the non-linear model will give a slightly different function from the linear model as the frequency of the interference pattern is being fit (and comes out to be 8.9V instead of 9V assumed in the linear fit). Here is the program.
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(1.,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")