Model Comparison and Hierachical Modeling

Takaaki Asano, Jaehyun Song
2016-07-08

Basic

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

General Formular and the Bayes Factor

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

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

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

General Formular and the Bayes Factor

\[ \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で今回のデータが得られる確率

Example: Two Factories of Coins

  • 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 \))

Solution by foraml Analysis: Step 1

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

\[ \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} \]

Solution by foraml Analysis: Step 2

  • コイン投げのデータ(\( 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) \)のみ。

Solution by foraml Analysis: Step3

  • 事前のベータ分布の\( 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) \)はベータ分布でなくベータ関数

Solution by foraml Analysis: Step4

手順

  • \( 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)} \]

Solution by foraml Analysis: Step 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 \)の事前分布の比で重みを付ける。

Solution by foraml Analysis: Step 6

事前分布は\( 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よりは良いモデルと評価できよう。

Solution by foraml Analysis: Step 7

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

たとえば\( \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%である。

Solution by foraml Analysis: Step +

教科書の範囲は超えるが、もし工場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%まで上がった。

Solution by grid approximation: Step 1

まず、データ、\( \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)

Solution by grid approximation: Step 2

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

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

Solution by grid approximation: Step 3

#⊂_ヽ 
#  \\  Λ_Λ 
#   \ ('ㅅ')  オメガ! 
#    > ⌒ヽ 
#   /   へ\ 
#   /  / \\ 
#   レ ノ   ヽ_つ 
#  / / カッパ!
#  / /| 
# ( (ヽ 
# | |、\ 
# | 丿 \ ⌒) 
# | |  ) / 
#(`ノ )  Lノ

Solution by grid approximation: Step 4

グリッドを作成する。

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

Solution by grid approximation: Step 5

グリッドの中身を見る。

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

Solution by grid approximation: Step 6A

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

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

Solution by grid approximation: Step 6B

\( \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)

Solution by grid approximation: Step 6C

plot of chunk show_marginal_theta2

Solution by grid approximation: Step 6D

\( \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"))

Solution by grid approximation: Step 6E

plot of chunk show_marginal_theta4

Solution by grid approximation: Step 7

尤度を計算する。

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

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

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

Solution by grid approximation: Step 8

計算結果を見る

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

Solution by grid approximation: Step 9-1

\( \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))

plot of chunk grid_show_likelihood1

Solution by grid approximation: Step 9-2

\( \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))

plot of chunk grid_show_likelihood2

Solution by grid approximation: Step 9

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

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

Solution by grid approximation: Step 10-1

モデル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 = "")))

plot of chunk unnamed-chunk-3

Solution by grid approximation: Step 10-2

モデル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 = "")))

plot of chunk show_grid_post_plot2

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

Solution by grid approximation: Step 10-3

比較してみる

plot of chunk show_grid_post_plot3

Solution by grid approximation: Step 10-3

周辺分布として見ると

plot of chunk show_grid_post_plot4

Solution by grid approximation: Step 11

  • しかし、ここでは\( \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))

plot of chunk show_grid_marginal_omega

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

Solution by grid approximation: Step 12

  • データ\( 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の時の同じ値が得られた。

Solution by MCMC

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

Ito Hiroki - Dealing with latent discrete parameters in Stan [Slideshare]

潜在離散パラメータ [Github]

Prediction: Model Averaging

  • 一般的に考えられるのは比較したモデルの中で最も妥当なモデル(\( 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} \]

Model Complexity Naturally Accounted For (1)

  • 複雑(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 \) 一様分布

Model Complexity Naturally Accounted For (2)

例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 \))が良いモデルらしい。

Model Complexity Naturally Accounted For (3)

例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 \))がより評価される。

Caveats regarding nested model comparison

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

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

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

Extrem Sensitivity to Prior Distribution (1)

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

コインの例

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

Extrem Sensitivity to Prior Distribution (2)

例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

Extrem Sensitivity to Prior Distribution (3)

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

Priors of different models should be equally informed(1)

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

またコイン投げの例

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

Priors of different models should be equally informed(2)

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は安定する