ぶっ細工ながらCMH統計量計算関数作成。 Yとyの値が近くなるとえっらく誤差が出ます。 作りかえんとね。 とりあえず、テンポラリあっぷ。 from inspect import * import decimal from math import * class MyPyStat: def MHmethod(self,tbls): self.MH_RES = [] Y = decimal.Decimal(0) w_s = decimal.Decimal(0) wy_s = decimal.Decimal(0) for pt_tbls in tbls : a = decimal.Decimal(pt_tbls[0][0]) b = decimal.Decimal(pt_tbls[0][1]) c = decimal.Decimal(pt_tbls[1][0]) d = decimal.Decimal(pt_tbls[1][1]) w = decimal.Decimal(1/((1/a)+(1/b)+(1/c)+(1/d))) print w r1 = a + b; r2 = c + d; c1 = a + c; c2 = b + d tot = r1 + r2 print (a*d)/(b*c) y = decimal.Decimal(str(log((a*d)/(b*c)))) print y w_s = w_s + w wy_s = wy_s + w*y tbls_prob = dict({'mfr':[r1,r2],'mfc':[c1,c2],'tot':tot, 'est':[[r1*c1/tot,r2*c1/tot],[r2*c1/tot,r2*c2/tot]], 'ratio':[[a/r1,b/r1],[c/r2,d/r2]], 'w':w,'y':y}) print tbls_prob self.MH_RES = self.MH_RES + [tbls_prob] Y = (wy_s)/w_s print Y self.MH_CHISQ = decimal.Decimal(0) for solo_mh in self.MH_RES: w=solo_mh['w'];y=solo_mh['y'] print 'w : ' + str(w);print 'y : ' + str(y) self.MH_CHISQ = self.MH_CHISQ + w*(y-Y)*(y-Y) def print_mbr(self,name): print self.__dict__[name] mps = MyPyStat(); mytbs = [[[1011,81],[390,77]],[[383,66],[365,123]]] mps.MHmethod(mytbs) mps.print_mbr("MH_RES") mps.print_mbr("MH_CHISQ")
