ベイズの基礎：サンプリング利用

「課題21-2 Testing Prior（2つの確率は同じと仮定）」からの続き。


ここでは課題 21-2 の観測データ、薬Aでは 8 人の患者中 4 人が完治、薬Bでは 8 人の患者中 2 人が完治、を使用する。

完治する確率 p_A, p_B は共に不明なので、一様分布を事前確率とする

PRIOR(p_A) = 1, PRIOR(p_B) = 1

二つの確率 p_A, p_B は独立の関係なので

PRIOR(p_A, p_B) = 1

尤度は

LIKELIHOOD = p_A⁴ (1 - p_A)⁴  × p_B²(1 - p_B)⁶

よって事後確率は

POST = PRIOR × LIKELIHOOD = p_A⁴ (1 - p_A)⁴  × p_B²(1 - p_B)⁶

s を 完治、f を未完治とすると、s_A = 4, f_A = 4, s_B = 2, f_B = 6 、事後確率 p_A は beta(s_A + 1, f_A + 1), p_B は beta(s_B + 1, f_B + 1) 、よって

　p_A⁴ (1 - p_A)⁴  × から beta(5, 5) , p_B²(1 - p_B)⁶ から beta(3, 7)


サンプリング（乱数生成）

薬A, B のそれぞれの事後確率が beta(5, 5), beta(4, 7) と算出できた。この二つの分布の違いが、薬A, B の効用の違いとなる。そこで、二つの事後確率からそれぞれ 1,000 件の乱数を生成して比較する。

> pA <- rbeta(1000,5,5)
> pB <- rbeta(1000,3,7)
> df<-data.frame(pA,pB)
> plot(df)
> abline(0,1,col="red")
> text(0.3,0.6,"pA<pB")
> text(0.9,0.6,"pA>pB")

赤い線は y = x の直線で、この線の下にある p_A > p_B の数が目視からも多いと確認できる。この線を境界に二分した区間にあるプロット数は

> nrow(subset(df,pA < pB))
[1] 160
> nrow(subset(df,pA > pB))
[1] 840

つまり 16 : 84 で薬Aが上回っている。一般的な認識では、16 : 84 から薬Aが効果がありそうに見えるが、そうだろうか？ HDI でどう判断できるか？

p_A と  p_B  の差 p_A - p_B 、すなわち「薬AがBより上回った度数」をヒストグラムにしたのが左図。

> d <- df$pA - df$pB
> h<-hist(d,xlab="pA-pB",ylab="",main="",
breaks=seq(-0.6,1.0,by=0.2))

各度数は以下の通り

> h$counts
[1]   1  20 139 333 340 139  28   0
> sum(h$counts)
[1] 1000

ヒストグラムで、度数が高い 4 区間 [ -0.2, 0.6 ] の度数とその合計は

> h$counts[3:6]
[1] 139 333 340 139
> sum(h$counts[3:6])
[1] 951

95.1% の確率区間にマイナス区間 [ -0.2, 0 ] が含まれていることから、薬Aの効能が高いと確信をもって結論はできない。

ちなみに 840 である p_A > p_B の個数は、以下でも求まる

> sum(h$counts[4:7])
[1] 840

「84% が上回っているから...」という理由は、95% 信頼区間においては通用しない。


テキスト終了

この課題、図らずも同時期に取組んでいたのランダムウォークと重なる。ランダムウォークの事後確率からサンプリングの採取は、今回のベータ関数の乱数生成に対応する。Beyesian の学習の王道の流れを感じてしまう（笑）。

さて、この投稿で本テキストに関する投稿は終わる。現在は未定だが、「ベイズの基礎」として良書の候補もあるので、面白いトピックが見つかり次第、投稿は続けていくつもり。

「目次『Workshop Statistics: Discovery with Data, A Bayesian Approach』」に続く。