Model Comparison and Hierachical Modeling

• 階層モデリングの延長線
• これまでは一つのモデル内の連続型パラメーターを推定するために事前分布・超事前分布
• 複数のモデルを用意し、超事前分布をモデルのインデックスとして考える。

• あるデータが生成される過程は一つではない。= 複数のモデルが存在しうる
• モデル1とモデル2、どっちが正しいか
• データが得られた時、あるモデルが正しい確率は？

\[ p(m|D) = \frac{p(D|m)p(m)}{p(D)} \]

• モデル・インデックスを除くパラメーターを周辺化する。

\[ \frac{p(m = 1|D)}{p(m = 2|D)} = \underbrace{\frac{p(D|m = 1)}{p(D|m = 2)}}_{\textrm{Bayes Factor(BF)}} \frac{p(m = 1)}{p(m = 2)} \]

• \( p(m = 1) \):モデル1の確率(事前分布)
• \( p(m = 2) \):モデル2の確率(事前分布)
• \( p(D|m = 1) \):モデル1で今回のデータが得られる確率
• \( p(D|m = 2 \):モデル2で今回のデータが得られる確率

• ２つの工場で作られたコイン

  • 工場1(\( m = 1 \)): 裏が出やすい(tail-biased)
    • バイアスは\( \omega = 0.25, \kappa = 12 \)のベータ分布に従う
  • 工場2(\( m = 2 \)): 表が出やすい(head-biased)
    • バイアスは\( \omega = 0.75, \kappa = 12 \)のベータ分布に従う

• どっちの工場のコインかに関する事前分布は\( p(m = 1) = p(m = 2) = 0.5 \)とする。

• データ: コインを9回(\( N \))投げ、6回表(\( z \))

コイン投げの結果が得られるプロセスは以下のとおりと仮定する。

\[ \begin{eqnarray} Y_i & \sim & \textrm{Bernoulli}(\theta) \\ \theta & \sim & \textrm{Beta}(\omega(\kappa - 2) + 1, (1 - \omega)(\kappa - 2) + 1) \\ \omega & \sim & \textrm{Categorical}( K = (0.25, 0.75), p = (p(m=1), p(m=2)))\\ p(m = 1) & = & p(m = 2) = 0.5 \\ \kappa & = & 12 \end{eqnarray} \]

• コイン投げのデータ(\( D \))が得られてから、工場\( m \)のコインである確率は

\[ p(m|D) = \frac{p(D|\theta, m)p(\theta|m)p(m)}{p(D)} \]

• これを\( \theta \)に対して周辺化すると

\[ \begin{eqnarray} p(m|D) & = & \frac{p(D|m)p(m)}{p(D)} \notag \\ p(D|m) & = & p(D|\theta, m)p(\theta|m) \notag \end{eqnarray} \]

• \( p(m) \)は0.5であり、\( p(D) \)は後でどうせなくなるので関心があるのは\( p(D|m) \)のみ。

• 事前のベータ分布の\( a \)と\( b \)が与えられている場合の\( p(D) = p(N, z) \)は(教科書 p.132)

\[ p(z, N) = \frac{B(z + a, N - z + b)}{B(a, b)} \]

• \( B(a, b) \)はベータ分布でなくベータ関数

手順

• \( p(\theta|m = 1) \)だと、\( \omega = 0.25, \kappa = 12 \)なので

\[ \begin{eqnarray} a & = & \omega (\kappa - 2) + 1 = 3.5\\ b & = & (1 - \omega) (\kappa - 2) + 1 = 8.5\\ p(D|m) & = & \frac{B(z + a, N - z + b)}{B(a, b)} \\ p(D|m) & = & \frac{B(6 + 3.5, 3 + 8.5)}{B(3.5, 8.5)} \end{eqnarray} \]

• \( m = 2 \)でも同様の手順で計算したら

\[ p(D|m = 2) = \frac{B(6 + 8.5, 3 + 3.5)}{B(8.5, 3.5)} \]

\( p(D|m = 1) \)と\( p(D|m = 2) \)は

pDm1 <- beta(9.5, 11.5) / beta(3.5, 8.5)
pDm2 <- beta(14.5, 6.5) / beta(8.5, 3.5)

したがってBayes Factor(BF)は

BF <- pDm1 / pDm2
print(BF)

[1] 0.2135266

ここに\( \omega \)の事前分布の比で重みを付ける。

事前分布は\( p(m = 1) = p(m = 2) = 0.5 \)なのでこれを掛けるとPosterior oddsが得られる。

(pm1D <- (pDm1 * 0.5)/(pDm2 * 0.5))

[1] 0.2135266

(pm2D <- (pDm2 * 0.5)/(pDm1 * 0.5))

[1] 4.683258

pm2Dが3以上なので、pm1Dよりは良いモデルと評価できよう。

一応、各モデルが合ってる確率も計算できる。

たとえば\( \frac{p(m = 1|D)}{p(m = 2|D)} \)が0.2135266ということは

\[ \frac{p(m = 1 | D)}{1 - p(m = 1 | D)} = 0.2135266 \]

を意味する。(\( p(m = 2|D) = 1 - p(m = 1|D) \)だから)

したがって、\( p(m = 1|D) = \frac{0.213}{1.213} \)

print(0.213/1.213)

[1] 0.1755977

約17.6%である。

教科書の範囲は超えるが、もし工場1の人員が工場2より二倍で、最新設備を備えていると仮定する。

適当に工場1のコイン生産量が工場2より2倍なら

(BF * 2); (1/BF * 0.5); (0.4270531/1.4270531)

[1] 0.4270531

[1] 2.341629

[1] 0.2992552

工場1である確率が約17.6%から約30%まで上がった。

まず、データ、\( \omega \)、\( \kappa \)を予め用意する。grid.gapはグリッドの刻み単位を意味する。

grid.gap  <- 0.01
omega     <- c(0.25, 0.75)
kappa     <- 12
coin.data <- c(1, 1, 1, 1, 1, 1, 0, 0, 0)

\( \omega \)と\( \kappa \)が与えられた時にベータ分布のパラメータを計算する関数(clac.shape)を作成する。

calc.shape <- function(omega, kappa){
    alpha <- omega * (kappa - 2) + 1
    beta  <- (1 - omega) * (kappa - 2) + 1
    return(list(alpha, beta))
}

#⊂_ヽ 
#　 ＼＼  Λ_Λ 
#　　 ＼ ('ㅅ')  オメガ！ 
#　　　 >　⌒ヽ 
#　　　/ 　 へ＼ 
#　　 /　　/　＼＼ 
#　　 ﾚ　ノ　　 ヽ_つ 
#　　/　/ カッパ！
#　 /　/| 
#　(　(ヽ 
#　|　|、＼ 
#　| 丿 ＼ ⌒) 
#　| |　　) / 
#(`ノ )　　Lﾉ

グリッドを作成する。

Grid.df <- expand.grid(omega, kappa, seq(0, 1, grid.gap))
names(Grid.df) <- c("omega", "kappa", "theta")

グリッドの内の\( \theta \), \( \omega \)と\( \kappa \)を使って事前確率を計算する。

Grid.df <- Grid.df %>% mutate(Prior = dbeta(theta, 
                                            calc.shape(omega, kappa)[[1]],
                                            calc.shape(omega, kappa)[[2]]))

グリッドの中身を見る。

head(Grid.df)

  omega kappa theta        Prior
1  0.25    12  0.00 0.000000e+00
2  0.75    12  0.00 0.000000e+00
3  0.25    12  0.01 7.936872e-03
4  0.75    12  0.01 8.345905e-13
5  0.25    12  0.02 4.160600e-02
6  0.75    12  0.02 1.472908e-10

\( \theta \)の周辺分布を計算する。

Mar.T <- Grid.df %>% group_by(theta) %>% 
    summarise(MT = sum(Prior))

\( \theta \)の周辺分布を表示する。

plot(x = seq(0, 1, grid.gap), y = Mar.T$MT, type = "l",
     lwd = 2, xlab = expression(theta), 
     ylab = expression(paste("p(", theta, ")", sep = "")),
     main = "Prior", cex.lab = 1.25)

\( \omega = 0.25 \), \( \omega = 0.75 \)の場合の\( \theta \)の事前分布も一緒に見てみよう。

lines(x = seq(0, 1, grid.gap), y = filter(Grid.df, omega == 0.25)$Prior, 
      col = "red", lty = 2, lwd = 2)
lines(x = seq(0, 1, grid.gap), y = filter(Grid.df, omega == 0.75)$Prior, 
      col = "blue", lty = 2, lwd = 2)
legend("topright", col = c("red", "blue", "black"), lty = c(2, 2, 1), lwd = rep(2, 3),
       legend = c(expression(paste(omega, "= 0.25")),
                  expression(paste(omega, "= 0.75")),
                  "Marginal"))

尤度を計算する。

Grid.df <- Grid.df %>% mutate(Likelihood = dbinom(sum(coin.data), length(coin.data), theta))

事前分布と掛け算したら事後分布(正規化前)

Grid.df <- Grid.df %>% mutate(Post = Likelihood * Prior)

計算結果を見る

head(Grid.df)

  omega kappa theta        Prior   Likelihood         Post
1  0.25    12  0.00 0.000000e+00 0.000000e+00 0.000000e+00
2  0.75    12  0.00 0.000000e+00 0.000000e+00 0.000000e+00
3  0.25    12  0.01 7.936872e-03 8.150512e-11 6.468957e-13
4  0.75    12  0.01 8.345905e-13 8.150512e-11 6.802339e-23
5  0.25    12  0.02 4.160600e-02 5.059848e-09 2.105201e-10
6  0.75    12  0.02 1.472908e-10 5.059848e-09 7.452690e-19

\( \omega = 0.25 \)の場合の尤度は

plot(x = seq(0, 1, grid.gap), y = filter(Grid.df, omega == 0.25)$Likelihood, type = "l", lwd = 2, xlab = expression(theta), ylab = "Likelihood", cex.lab = 2, ylim = c(0, 1))

\( \omega = 0.75 \)の場合の尤度は

plot(x = seq(0, 1, grid.gap), y = filter(Grid.df, omega == 0.75)$Likelihood, type = "l", lwd = 2, xlab = expression(theta), ylab = "Likelihood", cex.lab = 2, ylim = c(0, 1))

• どっちもかわらない。尤度そのものは\( m \)とは無関係のようだ。

• 何か変わるんだとしたら、事前分布(\( p(m) \))と掛け算した事後分布だろう。早速見てみる。

モデル1(\( \omega = 0.25 \))の場合の事後分布は

plot(x = seq(0, 1, grid.gap), y = filter(Grid.df, omega == 0.25)$Post, 
     xlim = c(0, 1), type = "l", lwd = 2, cex.lab = 2,
     xlab = expression(theta),
     ylab = expression(paste("P(", theta, "|", omega, " = 0.25, D)", sep = "")))

モデル2(\( \omega = 0.75 \))の場合の事後分布は

par(mar = c(5.1, 5.1, 2.1, 2.1))
plot(x = seq(0, 1, grid.gap), y = filter(Grid.df, omega == 0.75)$Post, 
     xlim = c(0, 1), type = "l", lwd = 2, cex.lab = 2,
     xlab = expression(theta),
     ylab = expression(paste("P(", theta, "|", omega, " = 0.25, D)", sep = "")))

par(mar = c(5.1, 4.1, 4.1, 2.1))

比較してみる

周辺分布として見ると

• しかし、ここでは\( \theta \)よりはどっちの\( \omega \)が良いモデルかが見たいだけだ。

MO <- Grid.df %>% group_by(omega) %>% summarize(MO <- sum(Post))
par(mar = c(5.1, 5.1, 2.1, 2.1))
plot(x = c(0.25, 0.75), y = MO$MO, xlim = c(0, 1), cex.lab = 2,
     type = "h", lwd = 5, col = "skyblue", ylim = c(0, 25),
     xlab = expression(omega), xaxt = "n",
     ylab = expression(paste("P(", omega, "|D)", sep = "")))
axis(1, at = c(0.25, 0.75), labels = c(0.25, 0.75))

par(mar = c(5.1, 4.1, 4.1, 2.1))

• データ\( D \)が得られた場合、\( \omega \)である確率は

# P(m = 1|D)
print(MO$MO[1])

NULL

# P(m = 2|D)
print(MO$MO[2])

NULL

# Bayes Factor
print(MO$MO[1]/MO$MO[2])

numeric(0)

• Formal solutionの時の同じ値が得られた。

• 教科書ではJAGSを使うが、対象外だからやめよう
• これまでの例で\( P(\omega) \)はカテゴリカル分布であったが、実はStanでカテゴリカルなパラメーターを直接指定するのは出来ない。
  • Stanの問題というよりはHMCの限界かなと
• 潜在離散変数を使えば迂回してできるらしい
• 例としては下のリンクを参照

• 一般的に考えられるのは比較したモデルの中で最も妥当なモデル(\( m = b \))を用いること

\[ p(\hat{y}|D, m = b) = \int d \theta_b p_b(\hat{y}|\theta_b, m = b)p_b(\theta_b|D, m = b) \]

• しかし、データが階層的なプロセスで生成されたものなら、他のモデルも考慮に入れる必要がある。
• \( p(\theta | D, m = 1) \)と\( p(\theta | D, m = 2) \)に\( p(m) \)で重みをかけたものの和 \( \rightarrow \) model averaging

\[ \begin{eqnarray} P(\hat{y}|D) & = & \sum_{m} p(\hat{y}|D, m)p(m|D) \notag \\ & = & \sum_{m} \int d \theta_m p_m(\hat{y}|D_m, m)p_m(\theta_m|D, m)p(m|D) \notag \end{eqnarray} \]

• 複雑(complex)なモデルは単純(simple)なモデルより手元のデータをよりよく表現できる。
• しかし、ノイズのあるデータの場合、常に複雑なモデルが良いとは言えない
• Bayesian model comparisonではパラメータの事前分布によって複雑なモデルが好まれる状況を補正できる。

コインの例

• Simple model: \( \omega_s = 0.5, \kappa_s = 1000 \) \( \leftarrow \) 尖ってる
• Complex model: \( \omega_c = 0.5, \kappa_c = 2 \) \( \leftarrow \) 一様分布

例1) コインを20回(\( N \))投げて15回(\( z \))表が出た場合

N <- 20; z <- 15
s.a <- 0.5 * (1000 - 2) + 1
s.b <- (1 - 0.5) * (1000 - 2) + 1
c.a <- 0.5 * (2 - 2) + 1
c.b <- (1 - 0.5) * (2 - 2) + 1
s.model <- beta(z + s.a, N - z + s.b) / beta(s.a, s.b)
c.model <- beta(z + c.a, N - z + c.b) / beta(c.a, c.b)
print(s.model/c.model)

[1] 0.3229023

複雑なモデル(\( \omega_c = 0.5, \kappa_c = 2 \))が良いモデルらしい。

例2) コインを20回(\( N \))投げて11回(\( z \))表が出た場合

N <- 20; z <- 11
s.a <- 0.5 * (1000 - 2) + 1
s.b <- (1 - 0.5) * (1000 - 2) + 1
c.a <- 0.5 * (2 - 2) + 1
c.b <- (1 - 0.5) * (2 - 2) + 1
s.model <- beta(z + s.a, N - z + s.b) / beta(s.a, s.b)
c.model <- beta(z + c.a, N - z + c.b) / beta(c.a, c.b)
print(s.model/c.model)

[1] 3.337148

ここではむしろ簡単なモデル(\( \omega_s = 0.5, \kappa_s = 1000 \))がより評価される。

• Nested modelとは Full modelに内包されるReduced model

• Full modelの方が複雑なので、一般的なモデル比較ではfull modelの方が勝る

• しかし、Bayesian model comparisonにおいては先と同様、Reduced modelによって手元のデータがちゃんと説明できればFull modelに勝ることもあり得る。

• Bayes Factorは「尤度 \( \times \) 事前分布」なので事前分布の影響に敏感

コインの例

• Simple model: 事前分布として\( \textrm{Beta}(500, 500) \)
• Complex model1: 事前分布として\( \textrm{Beta}(1, 1) \)
• Complex model2: 事前分布として\( \textrm{Beta}(0.01, 0.01) \)

例2) コインを100回(\( N \))投げて65回(\( z \))表が出た場合

N <- 100; z <- 65
s.model  <- exp(lbeta(z + 500, N - z + 500) - lbeta(500, 500))
c.model1 <- exp(lbeta(z + 1, N - z + 1) - lbeta(1, 1))
c.model2 <- exp(lbeta(z + 0.01, N - z + 0.01) - lbeta(0.01, 0.01))
print(s.model/c.model1)

[1] 0.125287

複雑なモデル1の方が好まれるが、同じデータで複雑モデル2と比べればどうだろうか。

print(s.model/c.model2)

[1] 5.728066

お、比較結果、複雑モデル1>単純モデル>複雑モデル2

• 一般的なベイズ推定の場合、ある程度のデータさえあれば、事前分布の不確実性はそこまで問題にならない(ロバストである)
• しかし、モデル比較においては重要

• どうする？
• データの一部を切り取って、その分を事前分布に渡し、残りでモデル比較

またコイン投げの例

• 先ほどの三つのモデルで10%を取った結果、10回投げ-6回表だとする。
• ならば、\( N = 100-10, z = 65-6 \)となる。
• とられた10%のデータは事前分布に渡す

N <- 100 - 10; z <- 65 - 6
s.model  <- exp(lbeta(z + 500 + 6, N - z + 500 + 10 - 6) - lbeta(500 + 6, 500 + 10 - 6))
c.model1 <- exp(lbeta(z + 1 + 6, N - z + 1 + 10 - 6) - lbeta(1 + 6, 1 + 10 - 6))
c.model2 <- exp(lbeta(z + 0.01 + 6, N - z + 0.01 + 10 - 6) - lbeta(0.01 + 6, 0.01 + 10 - 6))
print(s.model/c.model1)

[1] 0.05570509

print(s.model/c.model2)

[1] 0.05748123

どのケースにおいても複雑モデルがが単純モデルに勝った。

• 両モデルに同じ大きさのデータの一部(small amount of representative data)を渡すだけで、Bayes Factorは安定する