MCMC（Markov Chain Monte Carlo、マルコフ連鎖モンテカルロ法）について

今日はMCMC法についての解説です。

メモ程度のものですが、ご参考になれば幸いです。

日本語の良本はこれ。

有名な解説論文：

Sampling-Based Approaches to Calculating Marginal Densities.

Gelfand AE and Afrian F. M. Smith.

Journal of the American Statistical Association, 85;410:398-409, 1990.

【概念】

Monte Carlo（モンテカルロ法）

• モンテカルロ：金持ちの町、F1もやってる
• 興味のある値を「頻度」を使って推定する
  • 円の面積を求める
    • 正方形の中に円を描いてランダムに点を打つ
    • 何回も繰り返して円の中に点が入った割合を求める
    • 円の面積の近似値とする
  • 昔は乱数を発生する事が一苦労だったのでとても時間がかかった
  • 「金持ちしか出来ないような方法」という意味でネーミングされた
  • 今はPCの性能が上がっているので、大規模データでなければそんなに時間はかからない
• モンテカルロ積分と表現される事もある

Markov Chain（マルコフ連鎖）

• 一般マルコフ：次に起こる事象の確率分布が、それ以前の事象に影響される
• 1次マルコフ：次に起こる事象の確率分布が、直前（今）の事象に影響される
  • だいたいの手法はこれ
  • 2次以上のマルコフを使う場面を見た事はない

MCMC

• 以上の2つをミックス

例1. 条件付き分布は分かっているけどパラメータを推定するために積分計算をして周辺分布を得ることが困難

• 階層モデル
• 混合効果モデル

例2. サンプルの数が少なくて最尤法ではパラメータの値が不安定になる

• ロジスティック回帰などで説明変数にカテゴリ変数や交互作用項が多い場合
• 方法
  • 初期値を与えてパラメータが収束するまでマルコフ連鎖を行う（Markov Chain）
  • 何回も繰り返しパラメータの分布を得る（Monte Carlo）
  • 中央値や平均値をパラメータの推定値とする

Bayes

• 事後確率 = 尤度×事前確率
  • 事前確率：前もって得ている知識を確率で表現したもの
  • 尤度：得られたデータを確率で表現したもの
  • 事後確率：事前確率を尤度によって更新した確率
• モデル選択の枠組みでは確率→オッズ、尤度→ベイズファクターとも呼ばれる
  • ベイズファクターはBIC（ベイズ情報量基準）に近似可能
• MCMCの「条件付き分布を使って分布を更新していく」というところがベイズ的

【計算方法】

• Gibbs sampling（ギブスサンプリング）が最も有名（提案者：Geman and Geman）
  • 画像修復、ニューラルネットワーク、システム系で使われていた

• 同時分布を計算できなくても、各パラメータの条件付き分布が分かっていれば良い
  • 階層モデルに便利

【収束判定】

• Raftery and Lewisの診断
  • 分位点を使った診断
  • 診断のためには一定数以上のMCMC連鎖が必要（Lower Bound）
    • 正の自己相関があると連鎖の必要数が増える
    • Dependenceが自己相関（autocorrelation）の指標
      • 5より大きいと自己相関が強すぎる（”stickiness”と表現）
      • 初期値が不適格であることが疑われる
  • q：推定したい分位点
  • N：収束に（？）必要な連鎖回数

• Gelman-Rubin統計量（2つ以上の連鎖が必要）
  • 連鎖同士の平均値と個々の連鎖の値を比較する
  • 統計量が1に近ければ収束しているとみなす

• Gewekeの診断
  • 連鎖内で適当な分位点同士の差の検定を行う

【小ネタ】

• 私が居た教室で2006年あたりの修士論文が、「混合効果モデルに対してSASでギブスサンプリングを実装した」というテーマが3つくらいありました。
• @doryokujin君の修士論文は「MCMCが急速混合するために必要な試行回数の研究と改善」という内容をPythonを使って実装していた。というものだったようです。
  • 彼は学部のときにはParticle Filter (非線型ベイジアンフィルタ)について研究していたようです（難しい事やってたんだな～笑）

つまり、、

MCMCをプログラミングで実装できれば修士号を取れる！？

かもしれません笑

library(MCMCpack)

y <- c(1,3,3,3,5)
x1 <- c(1,2,3,5,5)
x2 <- c(2,3,1,5,2)
x <- matrix(c(x1,x2),ncol=2,nrow=5)

#---burn-in→mcmc
result3.sim <- MCMCregress(y~x, data = parent.frame(), burnin = 1000, mcmc = 10000) #返り値：パラメータの行列
result4.sim <- MCMCregress(y~x, data = parent.frame(), burnin = 1000, mcmc = 1000, B0 = 0.1)
plot(result3.sim)
summary(result3.sim)	#Naive SEはパラメータの標準誤差の推定値ではない

#---収束判定
#Raftery and Lewisの診断
raftery.diag(result3.sim)

#Gelman-Rubin統計量（2つ以上の連鎖が必要）
result34 <- mcmc.list(result3.sim, result4.sim)
gelman.diag(result34)

#Gewekeの診断
geweke.diag(result3.sim)


#---example
line   <- list(X = c(-2,-1,0,1,2), Y = c(1,3,3,3,5))
posterior  <- MCMCregress(Y~X, data=line, verbose=1000)	#verbose: この数毎に値がプリントされる
plot(posterior)
raftery.diag(posterior)
summary(posterior)


#普通の線形回帰
summary(glm(y ~ x))
plot(data.frame(x1, x2, y))




#MCMCpackのロジスティック回帰
data(birthwt)
posterior <- MCMClogit(low~age+as.factor(race)+smoke, data=birthwt)
plot(posterior)
summary(posterior)

#普通のロジスティック回帰
summary(glm(low~age+as.factor(race)+smoke, family=binomial, data=birthwt))