この記事は,ベイズ統計学勉強会2020年春合宿で使用したスライドの紹介記事です。
『社会科学のためのベイズ統計モデリング』の第4章「MCMC」の内容をベースにMCMCの仕組みについて解説しているスライドです。StanやJAGSのような確率的プログラミング言語を使えば理屈を知らずともMCMCを簡単に実行するができますが,MCMCの仕組みを知っていれば,うまく収束しないときに解決策を見つけやすくなるかもしれません。以下,スライド内で紹介しているRコードを(コピペしやすいように)掲載しておきます。
まずはメトロポリス・アルゴリズムのコードから (スライド15枚目)。このコードは『社会科学のためのベイズ統計モデリング』のpp.52-53で紹介されているコードを改変し,メトロポリス・アルゴリズムの各ステップをより分かりやすくしています。
Rコード01: メトロポリス・アルゴリズム
## つのドリルデータ ##
data <- rep.int(c(1,0),times=c(18,46-18))
## 確率モデル (尤度) ##
lik <- function(x,theta) theta^sum(x)*(1-theta)^(length(x)-sum(x))
## 事前分布 ##
prior <- function(theta) dbeta(theta,1,1)
## メトロポリス・アルゴリズム ##
n_iter <- 1000 #反復回数
#Step0: 初期値を設定
theta0 <- 0.10
MCMCsample <- 1:(n_iter+1) #結果を格納する変数 (初期値を含める)
MCMCsample[1] <- theta0
for(i in 1:n_iter){
#Step1: 一様乱数からtheta1を生成
theta1 <- runif(1,0,1)
#Step2: 事後オッズrを計算
posterior0 <- lik(data,theta0)*prior(theta0)
posterior1 <- lik(data,theta1)*prior(theta1)
r <- posterior1 / posterior0
#Step3: 遷移確率alphaを計算
alpha <- min(1,r)
#Step4: 次のtheta0を決定
theta0 <-ifelse(alpha>=runif(1,0,1),theta1,theta0)
#MCMCサンプルとして保存
MCMCsample[i+1] <- theta0
}
このコードを実行することで,R上で(メトロポリス・アルゴリズムによる)MCMCを実行することができます。確率モデル(尤度関数)や事前分布,データを変えて遊んでみてください。
続いて,スライド28枚目にあるマルコフ連鎖のシミュレーション用のRコードを掲載します。
Rコード02: マルコフ連鎖のシミュレーション
## 時点tの最大数 ##
max_t <- 10
## 推移確率行列 ##
P <- matrix(c(1/5,3/5,1/5,
1/3,1/3,1/3,
1/2,1/4,1/4),nrow=3,byrow=T)
## 初期の確率ベクトル ##
pi_t <- matrix(c(1,0,0),ncol=3)
## 全時点の確率ベクトルを格納する箱 ##
pi_all <- data.frame(t=0:max_t,ramen=0,chahan=0,bread=0)
pi_all[1,-1] <- pi_t
## マルコフ連鎖開始 ##
for (i in 1:max_t){
pi_t <- pi_t %*% P
pi_all[i+1,-1] <- pi_t
}
このコードを走らせればマルコフ連鎖が定常分布に収束することを確認できます。初期値や推移確率行列などを変えてみても楽しいです。
このスライドで解説しているMCMCのアルゴリズム(メトロポリス・アルゴリズムとMHアルゴリズム)は非効率的なのであまり実用的ではありませんが,MCMCの基本的な仕組みとして知っていれば何かと役立つと思います。