gerumaniumの日記

「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])

Permalink | コメント(2) | トラックバック(0) | 05:01