| 1 | import math
|
|---|
| 2 | import numpy
|
|---|
| 3 |
|
|---|
| 4 | def 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 |
|
|---|
| 11 | def 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 |
|
|---|
| 19 | class 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 |
|
|---|
| 82 | class 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 |
|
|---|
| 117 | def 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 |
|
|---|
| 142 | if '__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
|
|---|
