「データ解析のための統計モデリング入門」読書ノート 8章 MCMCとメトロポリス法

8章の図を再現する

今回は，本書の8章「マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル」の解析を実際に自分でやってみて，図を再現することを目的とする．基本的に本書の作図を真似するということを目標にして調整しているが，細部で少し違いがあるかもしれないので注意していただきたい．

問題設定とサンプルデータ

問題設定は6章と同じく，上限があるカウントデータを用いる．植物個体から8個種子を取ってきた時の種子の生存数が二項分布に従うとした上で，そのときの尤度を最大化し，生存確率であるパラメータqを推定することを目標とする．

まずはサンプルデータと二項分布の対数尤度関数の設定．

1
2

data <- c(4,3,4,5,5,2,3,1,4,0,1,5,5,6,5,4,4,5,3,4)
loglik <- function(q){data*log(q)+(8-data)*log(1-q)+log(choose(8,data))}

図 8.2・8.3

生存確率qと対数尤度の関係を図示する．8章の解析の大きな特徴としては，この後に出てくるメトロポリス法の説明のために，本来は連続変数である生存確率qを0から1までの0.01刻みの離散値として扱っている点にある．以下の図では，それらの生存確率に対応する対数尤度がマル印で示されてる．

また，解析的に求まる最大対数尤度とその時の生存確率q=0.45625を縦の赤線とバツ印で表示してある．今回は生存確率を離散値として表現しているため，最大対数尤度に対応する対数尤度の値は図示されていない．この図の中では，対数尤度の最大値を取るときの生存確率はq=0.46となっている．

1
2
3
4
5

a <- seq(0.01, 0.99, by=0.01)
plot(a, apply(matrix(a),1,function(z){sum(loglik(z))}),
     xlim=c(0.2,0.7), ylim=c(-60,-35), xlab="q", ylab="log L(q)")
abline(v=0.45625, lwd=2, col="red", lty=2)
points(0.45625, sum(loglik(0.45625)), pch=4, ,cex=2, col="red")

図 8.8

MCMCのアルゴリズムの一つであるメトロポリス法を実際に動かしてサンプリングをしてみる．

メトロポリス法の実装は本書にコード例が無かったので自作している．このコードで使っているアルゴリズムに関しては，本書P.176で述べられている手順を使っているわけではなく，スタンダードなメトロポリス法の手法を用いて先に対数尤度比を求めて条件分岐を作っている．細かい部分で違いはあれど基本的にやっていることは同じで作図もうまくいっているので，たぶん問題はないだろう．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

# Metropolis method
niter <- 100000
x <- 0.30
for(i in 2:niter){
  z <- sample(c(-0.01,0.01),1)
  r <- min(1, exp(loglik(x[i-1]+z)-loglik(x[i-1])))

  if(runif(1) < r){
    x[i] <- x[i-1] + z
  }else{
    x[i] <- x[i-1]
  }
}

# Fig. 8.8
par(mfrow=c(3,2))
plot(1:100,x[1:100], type="l", xlab="Iterations", main="100 MC steps")
plot(density(x[1:100]), xlim=c(0,1), main="")
plot(1:1000,x[1:1000], type="l", xlab="Iterations", main="1000 MC steps")
plot(density(x[1:1000]), xlim=c(0,1), main="")
plot(1:niter,x, type="l", xlab="Iterations", main="100000 MC steps")
plot(density(x), xlim=c(0,1),  main="")

図 8.9

メトロポリス法では，ステップ数を大きくしていくとサンプリングの値は定常分布に収束する．しかし，ステップ数が小さい段階では初期値に影響されるため，Burn-inとしてある程度のステップまでの結果を捨てるということが行われる．このように初期値が影響する様子を見るために，以下の図では初期値を0.1から0.9までの0.1刻みとした上でメトロポリス法を実行し，その結果を重ねあわせている．

これを見ると，ステップスを重ねる毎に初期値の影響が少なくなっていくことがわかるが，では実際にどのあたりまでを捨てるべきかという疑問については図を見ても判断できない．Burn-inの範囲設定は色々と手法があるものの，ある程度は実験者の裁量に任されているようだ．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

palette <- c("#000000","#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# http://wiki.stdout.org/rcookbook/Graphs/Colors%20(ggplot2)/

plot(1:500, rep(0,500), ylim=c(0,1), type="n", xlab="Iterations", ylab="q", main="500 MC steps")
for(init in seq(1,9)){
  niter <- 500
  x <- init*0.1
  for(i in 2:niter){
    z <- sample(c(-0.01,0.01),1)
    r <- min(1, exp(loglik(x[i-1]+z)-loglik(x[i-1])))

    if(runif(1) < r){
      x[i] <- x[i-1] + z
    }else{
      x[i] <- x[i-1]
    }
  }
  lines(x, col=palette[init], ylim=c(0,1), lwd=2)
}

参考