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.

Raw Paste

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