モデル比較：Pseudo-Priors 階層型MCMC

「非階層型 MCMC で p(D) 」からの続き。

左は図 10.2 、この図に従って実装した階層型モデルを使う。次は Jags-Ydich-Xnom1subj-MbernBetaModelComp.R のモデル定義。

model {
  for ( i in 1:N ) {
    y[i] ~ dbern( theta )
  }
  theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 ) 
  omega[1] <- .25
  omega[2] <- .75
  kappa <- 12
  m ~ dcat( mPriorProb[] )
  mPriorProb[1] <- .5
  mPriorProb[2] <- .5
}

theta ~ dbeta(omega[m]... から、θ のサンプリングは二つのモデルを「まとめて一緒」にしている。この実装で予想されるのは、「優位なモデルの方が多くサンプリングされる」ということ。

注意：この実装例は、psudo-prior の必要性の説明のためのもので、実用的ではない。

Jags-Ydich-Xnom1subj-MbernBetaModelComp.R の実行結果が左図で、本書の図 10.4 の下段に相当する（ 図 10.4 の上段は、コードのデータ部分をコメントアウトして prior を示したもの）。

先に指摘したように、この実装の問題点は「θ のサンプリング数の不平等さ」にある。Model Index のグラフから、モデル2が圧倒的に選択されているため、モデル間のサンプリング数の比率は

> length(thetaM1)/(length(thetaM1)+length(thetaM2))
[1] 0.17334
> length(thetaM2)/(length(thetaM1)+length(thetaM2))
[1] 0.82666

「モデル1：モデル2 = 17% : 83%」の大差となっている。つまり、このデータの場合はモデル2が「好ましいモデル」のため、「好ましくないモデル」の θ₁ は少ないサンプリング数の分布となっている。


Pseudo-Priors の導入

先の問題を解決するために導入するのが pseudo-priors（pseudo は「仮の」「擬似」という感じ）。以下は、pseudo-priors の実装スクリプト Jags-Ydich-Xnom1subj-MbernBetaModelCompPseudoPrior.R のモデル定義部分。

model {
  for ( i in 1:N ) {
    y[i] ~ dbern( theta )
  }
  theta <- equals(m,1)*theta1 + equals(m,2)*theta2
  theta1 ~ dbeta( omega1[m]*(kappa1[m]-2)+1 , (1-omega1[m])*(kappa1[m]-2)+1 ) 
  omega1[1] <- .10 # true prior value
  omega1[2] <- .40 # pseudo prior value
  kappa1[1] <- 20 # true prior value
  kappa1[2] <- 50 # pseudo prior value
theta2 ~ dbeta( omega2[m]*(kappa2[m]-2)+1 , (1-omega2[m])*(kappa2[m]-2)+1 ) 
  omega2[1] <- .70 # pseudo prior value
  omega2[2] <- .90 # true prior value
  kappa2[1] <- 50 # pseudo prior value
  kappa2[2] <- 20 # true prior value
  m ~ dcat( mPriorProb[] )
  mPriorProb[1] <- .5
  mPriorProb[2] <- .5
}

太字表記が pseudo-priors 。

pseudo-priors の求め方は二段階を要し、以下に各段階の実行例を示す。


pseudo-priors の求め方

段階#1：pseudo-priors に「true prior」を設定して実行

先のスクリプトのモデル定義で、以下のように true-priors と同じ値を pseudo-priors に設定する。

model {
  for ( i in 1:N ) {
    y[i] ~ dbern( theta )
  }
  theta <- equals(m,1)*theta1 + equals(m,2)*theta2
  theta1 ~ dbeta( omega1[m]*(kappa1[m]-2)+1 , (1-omega1[m])*(kappa1[m]-2)+1 ) 
  omega1[1] <- .10 # true prior value
  omega1[2] <- .10 # pseudo prior value
  kappa1[1] <- 20 # true prior value
  kappa1[2] <- 20 # pseudo prior value
theta2 ~ dbeta( omega2[m]*(kappa2[m]-2)+1 , (1-omega2[m])*(kappa2[m]-2)+1 ) 
  omega2[1] <- .90 # pseudo prior value
  omega2[2] <- .90 # true prior value
  kappa2[1] <- 20 # pseudo prior value
  kappa2[2] <- 20 # true prior value
  m ~ dcat( mPriorProb[] )
  mPriorProb[1] <- .5
  mPriorProb[2] <- .5
}

ここでの問題点は、実行結果のグラフの左側の Autocorrelation で ESS = 498.4 と低すぎる（MCMC のステップ数は 1000 回に設定）し、autocorrelation も発生している。また、Param. Values の図から、サンプリングの「動きの停滞」も見られる。

各 θ のサンプル数は

> length(theta1M1);length(theta1M2);length(theta2M1);length(theta2M2)

[1] 733

[1] 9267

[1] 733

[1] 9267

つまり、true-priors と pseudo-priors の合計のサンプリング数は二つのモデルで同じ。よって次の段階は、先の autocorrelation 等の問題解決と、true-priors と pseudo-priors の posteriors を近似させること。

段階#2：#1 の結果から true-priors の事後確率を「真似る」ような pseudo-priors を設定して実行

#1 の実行結果から、モデル1, 2 の pseudo-priors の事後確率は、それぞれ mode = 0.114, 0.885 と、pseudo-priors の 0.10, 0.90 に近似して、明らかに「データを表現していない」。

true-priors の事後確率を真似る ため、スクリプトの pseudo-priors を omega1[2] <- 0.40, omega2[1] <- 0.70 とする。これは、#1 の実行結果の true-priors の事後確率の mode = 0.414, 0.673 から。

そして pseudo-priors の k は

The k values were determined by remembering that the concentration parameter increases by N.

なので、スクリプトで kappa1[2] <- 50, kappa2[1] <- 50 に変更。

ESS = 9501.8 、autocorrelation も見られず、サンプリングの停滞もない。モデル比較は、実行結果の右上 p(m=1|D) = 0.0808, p(m=2|D) = 0.919 の通り。


あくまでもモデル比較

ところでこの結果から、モデル比較ではなく、各モデルのパラメータ θ₁, θ₂ を調べる場合、モデル1のサンプリング数は少なすぎる。

> length(theta1M1);length(theta1M2);length(theta2M1);length(theta2M2)

[1] 808

[1] 9192

[1] 808

[1] 9192

> length(theta1M1)/(length(theta1M1)+length(theta2M2))

[1] 0.0808

たかだか 800 のサンプリング数では少ない。対処方法は、単にステップ数を増やす（サンプリング時間は増える）などがあるが、詳細については本書に譲る。

「モデル比較の注意点」に続く。