PYTHON 30
###### Datafit.py Guest on 2nd May 2021 02:28:27 AM
1. #-------------------------
2. # datafit.py
3. # read data file
4. #------------------------
5. import matplotlib.pyplot as plt
6. from numpy import *
7. #------------------------------------------
8. def linfit(x,y,wt):
9. # Function to fit a line
10. # Inputs -
11. #   x - Independent variable
12. #   y - Dependent variable
13. #   wt - weight
14. # Outputs -
15. #   a_fit - Fit parameters; a(1) is intercept, a(2) is slope
16. #   sig_a - Estimated error in the parameters a()
17. #-----------------------------------------
18.     N =   size(x)
19.     s =   sum(wt)
20.     sx=dot(x,wt)
21.     sy=dot(y,wt)
22.     sxy=0.0;  sxx=0.0
23.     for k in num:
24.         sxy=sxy+x[k]*y[k]*wt[k]
25.         sxx=sxx+x[k]*x[k]*wt[k]
26. # % matrix method
27.     alpha = [[s, sx], [ sx, sxx]];
28.     betas  = [sy, sxy];
29. #
30.     errmat = linalg.inv(alpha);
31.     param  = dot(betas,errmat);
32. #  the fitted parameters
33.     a_fit = param
34. #    print ' a_fit  ' %a_fit
35. #    return a_fit
36. #  calculate errors from error matrix
37.     sig_a0 = sqrt(errmat[[0],[0]])
38.     sig_a1 = sqrt(errmat[[1],[1]])
39.     sig_a=[0,0]
40.     sig_a[0]=sig_a0
41.     sig_a[1]=sig_a1
42. #    print 'sig ', sig_a[0],sig_a[1]
43.     return a_fit,sig_a
44. #---------------------------------------
45. #   calculate chi sq ----------------
46. def chi2(y,yy,wtt):
47.     ylen=len(y)
48.     chs=0.0
49.     for iy in range(0,ylen):
50.         chs=chs+((y[iy]-yy[iy])**2)*wtt[iy]
51. #    print 'chi sq %8.2f' %chs
52.     return chs
53. #-----------------------------
54. # read data file
55. num=[] ; x=[]; errdat=[]; errors=[]; wtt=[]
56. f=open("lindata.txt")
57. k=0
58. for line in f:
59.     win=line.split()
60.     xin=float(win[0])
61.     datin=float(win[1])
62.     erin=float(win[2])
63.     x.append(xin)
64.     errdat.append(datin)
65.     errors.append(erin)
66. # weight factor
67.     wtt.append((1.0/erin)**2)
68.     num.append(k)
69.     k=k+1
70. #--------------------------------
71. #for k in num:
72. #    print '%d  %4.2f  %6.2f  %6.2f  %6.4f' %(k,x[k],errdat[k],errors[k],wtt[k])
73. #----------------------------------
74. # fit the line
75. a_fit,sig_a=linfit(x,errdat,wtt)
76. print "a_fit %6.2f  %6.2f" %(a_fit[0],a_fit[1])
77. print "sig_a %6.2f  %6.2f" %(sig_a[0],sig_a[1])
78. #-----------------------------------
79. # calculate chi sq for best fit
80. xlen=len(x)
81. yy=zeros(xlen)
82. inter=a_fit[0]
83. slpe=a_fit[1]
84. for iiy in range(0,xlen):
85.     yy[iiy]=inter + slpe*x[iiy]
86. chs=chi2(errdat,yy,wtt)
87. print ' %5d points      chi sq %6.2f' %(xlen,chs)
88. #-----------------------------------
89. # chi sq as a function of slope
90. slope=zeros(80); chst=zeros(80)
91. yyy=zeros(xlen)
92. for k in range(0,79):
93.     slope[k]=90.5+0.25*k
94.     for iiy in range(0,xlen):
95.         yyy[iiy]=inter+ slope[k]*x[iiy]
96.     chst[k]=chi2(errdat,yyy,wtt)
97. #-----------------------------------------
98. plt.figure()
99. # ---------  plot chi sq as a function of slope
100. #plt.subplot(211)
101. plt.plot(slope,chst,lw=0)
102. plt.errorbar(slope,chst,0,fmt='o')
103. plt.xlabel("slope")
104. plt.ylabel("chisq")
105. plt.ylim(chs-1,chs+4)
106. plt.xlim(slpe-4,slpe+4)
107. plt.grid(True)
108. #plt.title("%4d points   chisq %6.2f    slope %7.2f     inter %7.2f " %(xlen,chs,a_fit[1],a_fit[0]))
109. #-------- plot data and fitted line -----------
110. plt.figure()
111. #plt.subplot(212)
112. plt.plot(x,errdat,lw=0)
113. plt.plot(x,yy,lw=1)
114. plt.errorbar(x,errdat,errors,fmt='o')
115. plt.xlabel("x")
116. plt.ylabel("data")
117. plt.ylim(0,1100)
118. plt.grid(True)
119. plt.title("%4d points    chisq %6.2f    slope %7.2f   intercept %7.2f " %(xlen,chs,a_fit[1],a_fit[0]))
120. plt.show()

Paste-bin is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Recent Pastes

Raw Paste

or to edit or fork this paste. It's free.