ctest/matplotlib_fit_hist.py

   1 #! /usr/bin/env python
   2
   3 import matplotlib
   4 #matplotlib.use('GTK')
   5 matplotlib.use('Agg')
   6 import matplotlib.pyplot as plt
   7 import numpy as np
   8 from scipy.optimize import curve_fit
   9 from math import sqrt
  10 import sys
  11
  12 def prob_T(x,a):
  13     g=np.float128(105.)
  14     Tr=np.float128(300.)
  15     return a * (  x**((g-2)/2) * np.exp( -g*x/(2*Tr) ) )
  16
  17 #x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
  18 #x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
  19 #x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
  20 if (sys.argv[1] == '1L2Y_NH_GB000.stat' or sys.argv[1] == '1L2Y_NH_GB.stat'):
  21  x,y= np.loadtxt(sys.argv[1],usecols=(0,11),skiprows=10000,unpack=True)
  22 else:
  23  x,y= np.loadtxt(sys.argv[1],usecols=(0,10),skiprows=10000,unpack=True)
  24
  25 h,bin=np.histogram(y,bins=50,density=True)
  26
  27 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
  28 plt.xlim(min(bin), max(bin))
  29 plt.ylim(0,max(h))
  30 plt.ylabel('probality')
  31 plt.xlabel('temperature')
  32
  33 center = (bin[:-1] + bin[1:]) / 2
  34 #print bin
  35 #print center
  36 popt, pcov = curve_fit(prob_T, center, h, p0=10E-107)
  37 xfine = np.linspace(min(bin), max(bin), 100)
  38 #print popt
  39
  40 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
  41 print '%15.10f' % chi_squared
  42
  43 #print xfine
  44 #print prob_T(xfine,popt[0])
  45 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
  46 plt.savefig(sys.argv[1]+'.png')
  47 #plt.show()