ベイズの基礎：課題21-2 Testing Prior（2つの確率は同じと仮定）

課題21-1 どちらの薬が効果的？ からの続き。

前回はもっとも単純な観測データを使った比較だったが、今回は複数の観測データを使う。

Activity 21-2: Is the New Drug Better? (cont)

課題：課題21-1 の二つの薬の効果が同じである可能性が高いとして、事前確率を左表にした。


この事前確率は testing prior と呼ぶ。つまり、「二つの割合は同じ」ということをテストするという意味。

　同じ  ：0.5
　Aが良：0.25
　Bが良：0.25

今回は以下の 3 つの実験データを使う。

　実験#1：薬A は 8 人中、4 人が完治、薬B は 8 人中、2 人が完治
　実験#2：薬A は 16 人中、8 人が完治、薬B は 16 人中、4 人が完治
　実験#3：薬A は 32 人中、16 人が完治、薬B は 32 人中、8 人が完治

考え方は、前回の一人の場合と同様に考える。例えば、1 を完治、0 を未完治として薬 (A, B) の実験結果を表すと

　実験#1： (1, 1) が 2 、(1, 0) が 2 、(0, 0) が 4
　実験#2： (1, 1) が 4 、(1, 0) が 4 、(0, 0) が 8
　実験#3： (1, 1) が 8 、(1, 0) が 8 、(0, 0) が 16

Testing Prior 表のセルの尤度は、例えば左上隅セル (0.3, 0.3) は

　実験#1：(0.3 × 0.3)^2 × [ 0.3 × (1 - 0.3) ]^2 × [ (1 - 0.3) × (1 - 0.3) ]^4
　実験#2：(0.3 × 0.3)^4 × [ 0.3 × (1 - 0.3) ]^4 × [ (1 - 0.3) × (1 - 0.3) ]^8
　実験#3：(0.3 × 0.3)^8 × [ 0.3 × (1 - 0.3) ]^8 × [ (1 - 0.3) × (1 - 0.3) ]^16

セル (0.3, 0.5) の場合は

　実験#1：(0.3 × 0.5)^2 × [ 0.3 × (1 - 0.5) ]^2 × [ (1 - 0.3) × (1 - 0.5) ]^4
　実験#2：(0.3 × 0.5)^4 × [ 0.3 × (1 - 0.5) ]^4 × [ (1 - 0.3) × (1 - 0.5) ]^8
　実験#3：(0.3 × 0.5)^8 × [ 0.3 × (1 - 0.5) ]^8 × [ (1 - 0.3) × (1 - 0.5) ]^16

これ以降の事後確率の算出は、末尾に自作した R の実装例を参考。もっとスマートな実装も可能だろうが、その点は今回無視して素直な実装を優先した。

次が実行結果

> posts(2,2,4)
          [,1]       [,2]        [,3]
[1,] 0.2703677 0.04987169 0.004560555
[2,] 0.2715236 0.20033966 0.009160106
[3,] 0.1351838 0.04987169 0.009121109
> posts(4,4,8)
          [,1]      [,2]         [,3]
[1,] 0.2377382 0.0161781 0.0001352866
[2,] 0.4795509 0.1305339 0.0005457836
[3,] 0.1188691 0.0161781 0.0002705732
> posts(8,8,16)
           [,1]         [,2]         [,3]
[1,] 0.10042475 0.0009300963 6.504025e-08
[2,] 0.81722606 0.0302753619 1.058556e-06
[3,] 0.05021238 0.0009300963 1.300805e-07

注意：テキストの実験#3 の(0.5, 0.7) の確率が 0.067 は誤りだろう。この確率は実験#1, #2より高くなることはない。

各実験で高い値を示す (0.5, 0.3) の確率は

　実験#1：(0.5, 0.3) = 0.2715
　実験#2：(0.5, 0.3) = 0.4795
　実験#3：(0.5, 0.3) = 0.8172

さて、実はこの課題では、複雑な算出方法は分かったが、肝心の分析結果はよく分からない。もちろんテキストにもないので、中途半場だがここで終わる。ただ次回、サンプリングを使った分析では結論は示せるので、それを参照して頂きたい。


保留にした課題21-3

課題21-3 は、9 × 9 の 81 セルの Testing Prior だが、どうしてもテキストのような値を導き出すことができなかった。ある程度の考察まではいったが、確信が持てなかったので掲載は取りやめた。

別の書籍等を使って今後もベイジアンには取組むので、この課題は今のところは保留とする。

サンプリング利用 に続く。


R 実装例

# Function for Likelihoods
likelihood <- function(pA, pB, nBothGood, nBetterA, nBothBad) {
　(pA*pB)^nBothGood * (pA*(1-pB))^nBetterA * ((1-pA)*(1-pB))^nBothBad
}

# Matrix of Priors
priors <- matrix(c(2/12,1/12,1/12,1/12,2/12,1/12,1/12,1/12,2/12), nrow=3, ncol=3)

# Function for Posteriors
# [Parameters]
#  nBothGood: The number of cases that both drugs are good.
#  nBetterA : The number of cases that drug A are better than drug B.
#  nBothBad : The number of cases that both drugs are bad.
posts <- function(nBothGood, nBetterA, nBothBad) {
     # likeVec: 事前確率行列との積に使うため、[1,1][2,1][3,1][1,2]...のように行から格納。
  # 1st Col
  likeVec <- likelihood(0.3, 0.3, nBothGood, nBetterA, nBothBad)
  likeVec <- c(likeVec, likelihood(0.5, 0.3, nBothGood, nBetterA, nBothBad))
  likeVec <- c(likeVec, likelihood(0.7, 0.3, nBothGood, nBetterA, nBothBad))
  # 2nd Col
  likeVec <- c(likeVec, likelihood(0.3, 0.5, nBothGood, nBetterA, nBothBad))
  likeVec <- c(likeVec, likelihood(0.5, 0.5, nBothGood, nBetterA, nBothBad))
  likeVec <- c(likeVec, likelihood(0.7, 0.5, nBothGood, nBetterA, nBothBad))
  # 3rd Col
  likeVec <- c(likeVec, likelihood(0.3, 0.7, nBothGood, nBetterA, nBothBad))
  likeVec <- c(likeVec, likelihood(0.5, 0.7, nBothGood, nBetterA, nBothBad))
  likeVec <- c(likeVec, likelihood(0.7, 0.7, nBothGood, nBetterA, nBothBad))

  total <- sum(priors*likeVec)
  priors*likeVec/total
}

posts(2,2,4)
posts(4,4,8)
posts(8,8,16)