モデル比較：非階層型MCMCでp(D)

「Grid Approximation」からの続き。

「Formal Analysis」のように、モデル比較では Bayes Factor を求める。Bayes Factor とは p(D | m=1) / p(D | m=2) で、p(D | m) を求める必要がある。それを前回までの方法で求めるには、ごく少数のパラメータ分析には対応可能だが、多数のパラメータを扱う複雑なモデルでは現実的ではない。

そこで、用いるのが MCMC 。ここでは、一つのモデルの p(D | m) を算出方法を検討する。モデルが一つなので、以降は p(D) と記す。


数式による考察

まず、確率分布を近似する関数の「基本原則」を決める。次は式 10.7

日本語訳するより、原文が分かりやすいので引用。

For any function f (θ), the integral of that function, weighted by the probability distribution p(θ), is approximated by the average of the function values at points sampled from the probability distribution.

つまり、確率分布のサンプリング値の平均で近似できる、ということ。

次の式は、本書に分かりやすい解説があるが、ここでは式だけを示す。

ここまでの目的は、モデルの the marginal likelihood である p(D) = ∫ dθp(D|θ)p(θ) , 事前確率の p(θ) 。以下は p(D) で

MCMC サンプリング利用のために、この式を変形する。次は式 10.8 で、Bayes' rule から始まる。

2 行目で両辺を ∫ dθh(θ) = 1 を掛けるのが特徴。

h(θ) については、本書もしくは「Metropolis Algorithm の実装例 」の「Evidence p(D) の算出」を参照。


R 実装例

ここでは、非階層型 MCMCで p(D) の算出を実装する。

本書提供の R スクリプト（Jags-Ydich-Xnom1subj-Mbernbeta.R）で、サンプリングする theta を以下のように求める。

model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dbern( theta )
  }

  #theta ~ dbeta(0.75*(12-2)+1, (1-0.75)*(12-2)+1) # ω=0.75, k=12
  theta ~ dbeta(0.25*(12-2)+1, (1-0.25)*(12-2)+1) # ω=0.25, k=12
}

p(D) の算出コードは末尾に記した。以下は、式 10.8 の最後の行を実装した R コード

oneOverPD <- mean(dbeta(theta,aPost,bPost)/
    (theta^6*(1- theta)^(9-6) * dbeta(theta,omega*(kappa-2)+1,(1-omega)*(kappa-2)+1)))
PD <- 1/oneOverPD


結果は
　ω = 0.75 の場合 0.002338531
　ω = 0.25 の場合 0.0004993379

これは「Formal Analysis」で求めた値とほぼ同値（サンプリングなのでイコールにはなりにくい）。

「Pseudo-Priors 階層型 MCMC」に続く。


R コードサンプル

kappa <- 12; omega <- 0.25 # 0.25 or 0.75
myData <- c(rep(0,9-6),rep(1,6)) # 9 filips with 6 heads
mcmcCoda <- genMCMC(data=myData, numSavedSteps = 10000)
# Convert from coda object to vector
theta <- as.matrix(mcmcCoda)[,"theta"]

# Compute mean and sd of MCMC values
meanTheat <- mean(theta)
sdTheta <- sd(theta)
# Convert to a,b shape params for use in h(theta) function
aPost <- meanTheat * (meanTheat * (1 - meanTheat) / sdTheta^2 -1)
bPost <- (1-meanTheat) * (meanTheat * (1 - meanTheat) / sdTheta^2 -1)

# Compute 1/p(D)
oneOverPD <- mean(dbeta(theta,aPost,bPost)/
    (theta^6*(1- theta)^(9-6) * dbeta(theta,omega*(kappa-2)+1,(1-omega)*(kappa-2)+1)))
PD <- 1/oneOverPD
show(PD)

/** tex **/