「Rで学ぶベイズ統計学入門」をNumpy, Scipyで書き換えるでござるの巻
Rのpdisc()関数をpythonで書いた。
参考"yuji_nkのメモ"
"ハリ・セルダンになりたくて"
from numpy import * def pdisc(p,prior, data): p1=[] s=data[0] f=data[1] p1 = p + 0.5*(p==0) - 0.5*(p==1) p1=array(p1) like= s * log(p1) + f * log(1.0-p1) for i in range(len(p1)): if (p[i] ==0 and s>0) or (p[i] == 1 and f > 0): like[i] = -999.0 like = exp(like - max(like)) product = like * prior post = product/sum(product) return post
以下は実行結果
>>> p=linspace(5,95,10) >>> P=p*0.01 >>> prior=array([ 0.03460208, 0.1799308 , 0.27681661, 0.24913495, 0.15916955, ... 0.07266436, 0.02422145, 0.00346021, 0. , 0. ]) >>> data=array([11,16]) >>> pdisc(P,prior,data) array([ 1.45171410e-08, 2.25602200e-03, 1.29134895e-01, 4.76790956e-01, 3.33834971e-01, 5.58782010e-02, 2.09829685e-03, 6.64274607e-06, 0.00000000e+00, 0.00000000e+00])