モデル比較：Formal Analysis

「一般式と Bayes Factor」からの続き。

ここでは、説明のために極端な例を想定する。

設定例：コインを製造する工場が二つある。一方は裏に偏ったコインを製造して偏りの分布を ω₁ = 0.25、もう一方は表に偏った ω₂ = 0.75 という分布。このどちらかの工場で製造されたコイン投げの結果（表 9 回、裏 3 回）から、それがどちらの工場で作られたかの確率を見積もる。

左は図 10.2 で、二つのモデルの階層型モデル。二つのモデルとは、設定例の ω₁ と ω₂ のこと。左の点線枠が ω₁（裏に偏ったモデル）、右が ω₂（表に偏ったモデル）。

y_i のデータから、m = 1 , m = 2 の事後確率 p(m = 1), p(m = 2) を求める。事前確率は 50/50（どちらか不明）。

前回取り上げた左式 10.5 から、Bayes Factor を求めるため、各モデルの p(D | m) の値を算出する。

表の回数 z, 施行回数 N の Bernoulli 尤度関数から、Bayes' rule の分母 the marginal likelihood は次式 10.6 

これを R で関数定義すると

　pD <- function(z,N,a,b) { beta(z+a, N-z+b) / beta(a,b) }

この式ではアンダーフロー (underflow) の懸念があるので、自然対数で回避する

　pD <- function(z,N,a,b) { exp(lbeta(z+a, N-z+b) - lbeta(a,b)) }

a, b は図 10.2 の A₁, B₁, A₂, B₂ に対応して各モデル ω₁ = 0.25,  ω₂ = 0.75 毎の θ ~ beta(θ | ω(k-2)+1, (1-ω)(k-2)+1) より（ここでは k = 12 とした）

　θ ~ beta(θ | a1, b1) = beta(θ | 0.25(12-2)+1, (1-0.25)(12-2)+1) = beta(θ | 3.5, 8.5)
　θ ~ beta(θ | a2, b2) = beta(θ | 0.75(12-2)+1, (1-0.75)(12-2)+1) = beta(θ | 8.5, 3.5)

定義した pD 関数から Bayes Factor の分子と分母（各モデルの marginal likelihood）を求める

> (m1<-pD(z=6,N=9,a=3.5,b=8.5))
[1] 0.0004993439
> (m2<-pD(z=6,N=9,a=8.5,b=3.5))
[1] 0.002338556
> m1/m2
[1] 0.2135266

表と裏の偏りの確率が 50/50 なら、この結果から裏の偏りの odds は劣って、4.68 (1/0.213 or m2/m1) 表の方が勝っている

> m2/m1
[1] 4.683258

また式 10.5 で表すと

　p(m=1 | D) / p(m=2 | D) = (0.0004993439 / 0.002338556) * (0.5/0.5) = 0.213

m₁ か m₂ のいずれかなので、p(m=1 | D) + p(m=2 | D) = 1 の関係から

　p(m=1 | D) / p(m=2 | D) = p(m=1 | D) / (1 - p(m=1 | D)) = 0.213

これを変形して

　p(m=1 | D) = 0.213 / (1 + 0.213) = 17.6%

よって

　p(m=2 | D) = 1 - 17.6% = 82.4%


結論として、9 回の試行で 6 回の表のデータから、表の偏りの事前確率 p(m=2 | D) = 50% は p(m=2 | D) = 82.4% の事後確率へと「reallocated beliefs 確信を再編成した」。

注意すべきは、ここで行ったのはモデル比較であって「どちらのモデルの可能性が高いか」で、「コインの偏り分布 p(θ | D, m) ではない」ということ。つまり、どちらの工場がそのコインを製造したかの見積りであって、コインの偏りの見積りでなない。


k の影響

次は「実習 10.1」で、同様にどちらかに偏ったコインから、どちらの工場であるかの確率を求める。違いは観測データが「10 回の試行で 7 回の表」、そして k の値

本解答のための R スクリプトは末尾に記した。なお、コード中の

　p1gD <- BF12xPriorOdds / (1.0+BF12xPriorOdds)

は、先の p(m=1 | D) = 0.213 / (1 + 0.213) に対応し、0.213 が BF12xPriorOdds 。

(A) ω₁  = 0.25,  ω₂  = 0.75, k = 6 の事後確率は（どちらの工場で製造されたか）？
> show(BF12); show(p1gD); show(p2gD)
[1] 0.3333333
[1] 0.25
[1] 0.75

裏の偏りの工場：表の偏りの工場 = 25% : 75%

(B) ω₁  = 0.25,  ω₂  = 0.75, k = 202 の事後確率は（どちらの工場で製造されたか）？
> show(BF12); show(p1gD); show(p2gD)
[1] 0.01621596
[1] 0.0159572
[1] 0.9840428

裏の偏りの工場：表の偏りの工場 = 1.6% : 98.4%

(C) A と B で事後確率の違いの理由は？
データが「10 回中 7 回表」なので、「表に偏った工場」なのだが、どちらの工場であるかの確証は (A), (B) のように算出することで明確になる。この違いの理由は k で、k = 202 の (B) はデータが (A) よりも決定的な要因となった。

k = 202 の (B) は、コインが 0.25 か 0.75 の「どちらかの値に極めて近い」という強い確証。片や k = 6 の (A) では 0.25 か 0.75 の「どちらかもしれないが良く分からない」という感じで確証は弱い。そして、データから (B) は更に確証を得て、50% : 50% の事前確率から「1.6% : 98.4% に確率が更新」ということ。

注意点を繰り返すが、これは「どちらのモデルの可能性が高いか」のモデル比較であって、「コインの偏りの分布 p(θ | D, m) ではない」。(B) の場合「98.4% 表に偏っている」ではなく、「98.4% の確率で、表に偏った工場で製造されたコイン」ということ。

「Grid Approximation」に続く。


R コードサンプル

# The function for computing the marginal likelihood
pD <- function(z,N,a,b) {exp(lbeta(z+a,N-z+b)-lbeta(a,b))}
omega1 <- 0.25
omega2 <- 0.75
kappa <- 6 # 6 for (A), 202 for (B)
a1 <- omega1*(kappa-2) + 1
b1 <- (1-omega1)*(kappa-2) + 1
a2 <- omega2*(kappa-2) + 1
b2 <- (1-omega2)*(kappa-2) + 1
# data
z <- 7; N <- 10
# The marginal likelihood
pDg1 <- pD(z,N,a1,b1)
pDg2 <- pD(z,N,a2,b2)
# Bayes' factor
BF12 <- pDg1/pDg2
# Prior probs
p1 <- 0.5; p2 <- 1-p1
# Post probs
BF12xPriorOdds <- (pDg1/pDg2)*(p1/p2)
p1gD <- BF12xPriorOdds / (1.0+BF12xPriorOdds)
p2gD <- 1.0 - p1gD
# Results
show(BF12); show(p1gD); show(p2gD)