754540c8135c154576a2e693c50ba9ba50ff50ad
[unres.git] / 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     gg=np.float128(g)
14     aa=np.float128(10**(-gg-2)*a)
15     Tr=np.float128(300.)
16     return np.exp( np.log(aa) + (gg-2)/2*np.log(x) - gg*x/(2*Tr) )
17 #    return aa * (  x**((gg-2)/2) * np.exp( -gg*x/(2*Tr) ) )
18
19 #x,y= np.loadtxt('1L2Y_L_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
20 #x,y= np.loadtxt('1L2Y_NH_GB000.stat',usecols=(0,11),skiprows=10000,unpack=True)
21 #x,y= np.loadtxt('1L2Y_B_GB000.stat',usecols=(0,10),skiprows=30,unpack=True)
22 g=int(sys.argv[2])
23 if (sys.argv[1] == '1L2Y_NH_GB000.stat' or sys.argv[1] == '1L2Y_NH_GB.stat'):
24  x,y= np.loadtxt(sys.argv[1],usecols=(0,11),skiprows=10000,unpack=True)
25 else: 
26  x,y= np.loadtxt(sys.argv[1],usecols=(0,10),skiprows=10000,unpack=True)
27
28 h,bin=np.histogram(y,bins=50,density=True)
29
30 plt.bar(bin[:-1], h, width = bin[2]-bin[1])
31 plt.xlim(min(bin), max(bin))
32 plt.ylim(0,max(h))
33 plt.ylabel('probality')
34 plt.xlabel('temperature')
35
36 center = (bin[:-1] + bin[1:]) / 2
37 #print bin
38 #print center
39 popt, pcov = curve_fit(prob_T, center, h)
40 xfine = np.linspace(min(bin), max(bin), 100)
41 #print popt
42
43 chi_squared = np.sum((prob_T(center, *popt)-h)**2)
44 print '%15.10f' % chi_squared
45
46 #print xfine
47 #print prob_T(xfine,popt[0])
48 plt.plot(xfine,prob_T(xfine,popt[0]),'-',c='red')
49 plt.savefig(sys.argv[1]+'.png')
50 #plt.show()