MCMCアルゴリズムと呼ぶ

8.1‐8.3 MCMC
FM13009 庭月野翔太
発表の流れ
1. 概要
• MCMC法
2. 例題：種子の生存確率（個体差なし）
3. 試行錯誤による最尤推定
4. MCMCアルゴリズムのひとつ：メトロポリス法
• メトロポリス法によるサンプリング
• マルコフ連鎖の定常分布
• 推定された定常分布は何を表しているのか
5. まとめ
概要
• GLMMはランダム効果の発生源が個体差のみだった
• 場所差などランダム効果の発生源が増えるにつれパ
ラメータの推定が困難になる
 最尤推定値の計算時間の増加，推定値の探索の困難化
複雑な統計モデルのあてはめで有効なのがMCMC法．
MCMC法
• 正しくはマルコフ連鎖モンテカルロ法
• 資料では，MCMC法による様々な手続きをMCMCア
ルゴリズム，データに対しアルゴリズムを適用させる
ことをMCMCサンプリング，サンプリングによって得ら
れたデータをMCMCアルゴリズムと呼ぶ
例題：種子の生存確率（個体差なし）
種子の生存確率の計算（1）
• 個体𝑖の生存数は
｛ｙ₁，ｙ₂…ｙ₂₀｝=｛4,3,4,5,5,2,3,1,4,0,1,5,5,6,5,4,4,5,3,4}
• 生存種子数𝑦𝑖 が二項分布にしたがい，種子生存確率
を𝑞 とすると，ある個体𝑖の種子数が𝑦𝑖 である確率は
種子の生存確率の計算（2）
• 種子生存確率𝐿 𝑞 の尤度は
両辺の対数をとると
となり，log𝐿 𝑞 を最大化するような𝑞 が最尤推定値
の𝑞^ となる．
種子の生存確率の計算（3）
• 例題条件下においての最尤推定値は下式の通り
• 以上より，生存確率は0.46ぐらいであると推定される
試行錯誤による最尤推定
目的：試行錯誤による数値的な最尤推定
• 以下のルールのもとで𝑞を変化させ，対数尤度が高く
なる𝑞^ を探索させる．
① パラメータ𝑞の初期値を選ぶ
② パラメータ𝑞 を増やすか減らすかをランダムに決める
③ あたらしいパラメータ𝑞新において尤度が大きくなる（当てはまりが良くなる）ならば𝑞 の値
を𝑞新に変更する．
𝑞の値が変化しなくなるまで②～③をくり返す．
ふらふら試行の錯誤による最尤推定（２）
• 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)) }
• 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")
•
参考URL： http://yagays.github.io/blog/2012/11/20/glm-mcmc-chp8/
ふらふら試行の錯誤による最尤推定（3）
• 「対数尤度の高い方向にだけ移動する」というルールで𝑞を変化さ
せると初期値に関係なく最尤推定値（ 𝑞 =0.46）に変化していく．
MCMCアルゴリズムのひとつ：メトロポリス法
• MCMCアルゴリズムのひとつであるメトロポリス法の手順．
① パラメータ𝑞の初期値を選ぶ
② パラメータ𝑞 を増やすか減らすかをランダムに決める
③ あたらしいパラメータ𝑞新において尤度が大きくなる（当てはまりが良くなる）ならば𝑞 の値
を𝑞新に変更する
④ パラメータ𝑞新で尤度が小さくなる（あてはまりが悪くなる）場合でも確率𝑟で𝑞 の値を𝑞新
に変更する．確率𝑟は尤度比𝑟 =
い．）．
𝐿
𝑞新
𝐿 𝑞
に等しい（尤度差が小さいほど移動しやす
前の状態に基づいて新しい状態を作り出すことをマルコフ連鎖，乱数を利用したアルゴリズ
ムはモンテカルロ法と呼ばれる．
メトロポリス法によるサンプリング（1）
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="")
メトロポリス法によるサンプリング（2）
• メトロポリス法をデータに適用し，各ステップごとの生存確率𝑞と対
数尤度log𝐿 𝑞 の変化の過程を示したのが下図である．
メトロポリス法によるサンプリング（3）
•
•
MCMCアルゴリズムの目的はステップ数とともに変化する値の生成である．
下図のヒストグラムの横に示されている曲線の確率分布は，例題の統計モデ
ルとメトロポリス法によって決まるマルコフ連鎖の「定常分布」である．
マルコフ連鎖の定常分布
•
定常分布とは変数𝑞 のマルコフ連鎖が一定の条件を満たしているときに，そ
のマルコフ連鎖から発生する𝑞の値が従う確率分布である．
•
定常分布を推定できるような標本集団を得るためには，十分な数のMCMCサ
ンプリングが必要となる．Cの状態においては，定常分布が推定できるような𝑞
の値のサンプルが得られたとする．
推定された定常分布は何を表しているのか（1）
•
この章で取り上げた例題の定常分布𝑝 𝑞 𝑌 は尤度𝐿 𝑞 に比例する確率分布
である．次のように定義される．
推定された定常分布は何を表しているのか（2）
•
•
推定された確率分布𝑝 𝑞 𝑌 は𝑞の尤もらしさである尤度に比例するためあ
るデータ𝑌に統計モデルをあてはめたときに𝑞がとる値の確率分布と解釈する
ことができる．
統計量の評価を行い，あてはめの結果として得られた𝑞の部分を要約すれば
統計モデルによるさまざまな予測ができる可能性がある．
まとめ
• MCMC法は多数のパラメータを推定する方法である
• 最尤推定法は尤度最大になるパラメータを探索する
最適化である
• MCMCアルゴリズムは定常分布からのランダムサン
プリングが目的である．この章の例題では定常分布
は尤度に比例する確率分布である

Download Report