※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。



ぶっ細工ながらCMH統計量計算関数作成。
Yとyの値が近くなるとえっらく誤差が出ます。
作りかえんとね。
とりあえず、テンポラリあっぷ。(最下段)

某本の例題なんだけど

chisq on Book : 0.896
chisq on my Python PG : 0.929

Rの出力
確率密度
> dchisq(0.929,df=1)
[1] 0.2601192
> dchisq(0.896,df=1)
[1] 0.2692726

確率累積分布(いわゆる確率)
> pchisq(0.929,df=1)
[1] 0.6648771
> pchisq(0.896,df=1)
[1] 0.6561435 

んーーー、これは誤差と考えていいか???
RでPythonと同じ計算をさせた結果、
Pythonと等しかった。
つまり、正確でないのは本の解答か?

(11/23)
decimal→floatに変えてみる。
各要素に使う型は何が良いのだろう?

 CMH統計量計算プログラム

 from inspect import *
 import decimal
 from math import *
 import scipy
 
 
 class MyPyStat:
     
     def MHmethod(self,tbls):
         
         self.MH_RES = []
         Y = float(0)
         w_s = float(0)
         wy_s = float(0)
         for pt_tbls in tbls :
             a = float(pt_tbls[0][0])
             b = float(pt_tbls[0][1])
             c = float(pt_tbls[1][0])
             d = float(pt_tbls[1][1])
             w = float(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 = float(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 = float(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)
             
     
   #いわゆるreflection
     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")