DifferentialWeight: dist.py.txt

File dist.py.txt, 4.8 KB (added by trac, 7 years ago)
Line 
1import math
2import numpy
3
4def read_file(filepos,ratio=1):
5    out=[]
6    for line in file(filepos):
7        data=line.split()
8        out.append(ratio*float(data[1]))
9    return out
10
11def read_correlation(filepos):
12
13    text=file(filepos).read()
14    text=text[:-1].replace('\n',';')
15    out=numpy.matrix(text)
16    return out
17
18
19class Chi_carre:
20   
21    def __init__(self,expected,variance,correlation,ratio=1):
22        self.expected=read_file(expected,ratio)
23        self.variance=read_file(variance,math.sqrt(ratio))
24        self.correlation=read_correlation(correlation)
25        self.len=len(self.expected)
26       
27    def compute_from_value(self,dico,mode=1):
28        total=0
29        if 1==mode: # standard chi-carre
30            for key in range(0,self.len):
31                total+=((dico[key]-self.expected[key])/self.variance[key])**2       
32            return total/len(self.expected)
33        if 2==mode: #test ofcoherence
34            for key in range(0,self.len):
35                prov=0
36                fact1=(dico[key]-self.expected[key])/self.variance[key]
37                for key2 in range(0,self.len):
38                    prov+=fact1*(dico[key2]-self.expected[key2])/(self.variance[key2]*self.correlation[key,key2])
39                total+=prov
40            total=math.fabs(total)/(len(self.expected)**2)
41            return total
42        if 3==mode:
43            inv_mat=self.correlation.I
44            total=0
45            for key in range(0,self.len):
46                prov_num=0
47                prov_den=0
48                for key2 in range(0,self.len):
49                    prov_num+=inv_mat[key,key2]*(dico[key]-self.expected[key])
50                    prov_den=inv_mat[key,key2]*self.expected[key2]
51                total+=(prov_num/prov_den)**2
52            return total/float(self.len**2)
53        if 4==mode:
54            inv_mat=self.correlation.I
55            total=0
56            for key in range(0,self.len):
57                prov_num=0
58                prov_den=self.variance[key]
59                for key2 in range(0,self.len):
60                    prov_num+=inv_mat[key,key2]*(dico[key]-self.expected[key])
61                total+=(prov_num/prov_den)**2
62            return total/float(self.len**2)
63        if 5==mode:
64            cov=numpy.matrix(self.correlation)
65            for key in range(0,self.len):
66                for key2 in range(0,self.len):
67                    cov[key,key2]=cov[key,key2]*self.variance[key]*self.variance[key2]
68            inv_mat=cov.I
69            total=0
70            for key in range(0,self.len):
71                prov=0
72                for key2 in range(0,self.len):
73                    prov+=(dico[key]-self.expected[key])*inv_mat[key,key2]*(dico[key2]-self.expected[key2])
74                    total+=(dico[key]-self.expected[key])*inv_mat[key,key2]*(dico[key2]-self.expected[key2])
75                print prov
76            return total/float(self.len)
77
78    def compute_from_file(self,filepos,mode=1):
79        dico=read_file(filepos)
80        return self.compute_from_value(dico,mode=mode)
81
82class p_value:
83
84    def __init__(self):
85        self.p_value=[value/10.0 for value in range(0,26)]
86        self.ntot=0
87        self.value={}
88        for value in self.p_value:
89            self.value[value]=0
90
91    def read_file(self,filename):
92
93        for line in file(filename):
94            self.add_value(float(line))
95
96    def add_value(self,value):
97            self.ntot+=1
98            for pvalue in self.p_value:
99                if value>pvalue:
100                    self.value[pvalue]+=1
101                else:
102                    break
103
104    def compute_result(self):
105        self.result()
106        for value in self.p_value:
107            self.value[value]/=float(self.ntot)
108            self.value[value]*=100
109
110    def result(self):
111        text=''
112        for pvalue in self.p_value:
113            text+= str(pvalue)+' '+str(self.value[pvalue])+'\n'
114        return text
115
116
117def compute_pvalue():
118    stat_obj=[]
119    range_test=1
120    for i in range(0,range_test): stat_obj.append(p_value())
121#    obj1=Chi_carre('value','error','correlation')
122    obj2=Chi_carre('dist_value_200','dist_error_200','dist_correlation_200')
123
124    ff=open('chi_carre','w')
125    for i in range(0,10000):
126        for j in range(0,range_test):
127#            chi_carre=obj1.compute_from_file('10bin_'+str(i),mode=5)
128#            stat_obj[0].add_value(chi_carre)
129            chi_carre=obj2.compute_from_file('dist_'+str(i),mode=5)
130            stat_obj[0].add_value(chi_carre)
131            ff.writelines(str(chi_carre)+'\n')
132 
133    for j in range(0,range_test):
134        stat_obj[j].compute_result()
135        text=stat_obj[j].result()
136        print text
137    ff.close()
138    ff=open('pvalue','w')
139    ff.writelines(text)
140    ff.close()
141
142if '__main__'==__name__:
143    obj2=Chi_carre('MW_mean','MW_error','MW_correlation',1000)
144    chi_carre=obj2.compute_from_file('MW_SM_event',mode=5)
145    print chi_carre