MCMCについて

MCMC について 水産資源学におけるベイズ統計の応用ワークショップ 2007 年 8 月 2-3 日, 中央水研 遠洋水産研究所外洋資源部 鯨類管理研究室 岡村寛

事後分布からサンプルする Pr(θ 1, θ 2, θ 3, ) θ 1 は重要. あとは必要だけど直接的じゃないというようなとき, P(θ 1 x)= P(θ 1,, θ n x)dθ 2 dθ n を計算したい. モンテカルロ近似 E π (f(x)) = f(x)π(x)dx を計算したい (1/N)Σf(x (i) ) E π (f(x)), x (1),, x (N) ~ π(x) π(x) からデータを取り出せれば良い. しかし, ある確率分布からランダムにデータを取り出すのは一般に難しい

MCMC でない方法

逆関数法 まず分布関数 F(x) = Pr(X x) の逆関数を求める. 一様分布 U[0,1] から乱数 u を発生.F -1 (u) は分布関数 F(x) の確率分布からとったランダムサンプルとなる. なぜなら, Pr(F -1 (U) x) = Pr(U F(x)) = F(x) だから,F -1 (U) は X に対応している. 例 : 指数分布密度関数 f(x) = 1/2 exp(-x/2) 分布関数 F(x) = Pr(X x) = 1-exp(-x/2) U = F(v) v = - 2log(1-U)

逆関数法の例 R プログラム : u1 <- runif(100000,0,1) par(mfrow=c(1,1),mar=c(5,5,5,5)) plot(density(- 2*log(u1)),type="l",lty=1,lwd=3,col=2, xlim=c(0,25),ylim=c(0,0.5),xlab="",yla b="density",main="exponential distribution",cex.lab=2,cex.main=2) x1 <- seq(0,30,by=0.01) lines(x1,dgamma(x1,1,1/2),lty=2,lwd =2,col=1)

グリッド法 ( 復習 ) 1 つ 2 つのパラメータ, 細かい分析をする必要がなければ有効 ( 楽チン ) 1. パラメータのグリッドを用意する 2. 各値に事前確率を割り当てる ( 一様分布なら必要なし ) 3. 事後確率 ( 尤度 事前確率 ) を計算. 4. 事後確率に従ってグリッドの値を抽出する 事後サンプル

2 変数の場合 上の事後確率行列を M とする. 1. まず x の周辺事後分布から x を B 個サンプルする 2. 次に上でサンプルした x の条件付確率 p(y x) にしたがって y を B 個サンプルする 3. 得られた (x, y) は上の事後確率を持つ分布からのサンプルである

R プログラム x <- c(1,2,3); y <- c(1,2,3,4); B <- 1000 M <- matrix(c(0.017,0.067,0.017,0.083,0.2,0.03 3,0.083,0.117,0.150,0.1,0.133,0),nrow=3,n col=4) nx <- sample(nrow(m),b,replace=t,prob=rowsu ms(m)) x1 <- x[nx] f1 <- function(x) sample(ncol(m),1,replace=t,prob=m[x,]) ny <- apply(as.matrix(nx), 1, f1) y1 <- y[ny] x2 <- x1+rnorm(length(x1),0,0.05) y2 <- y1+rnorm(length(y1),0,0.05) par(mfrow=c(1,1),mar=c(5,5,5,5)) plot(y2,-x2,xlab="y",ylab="x",cex.lab=2)

シミュレーションによる資源管理方式 の頑健性評価 さまざまな不確実性 t-1 年末の個体数プロダクションモデル t 年はじめの個体数管理方式によって漁獲 t 年末の個体数 加入 自然死亡の変動自然災害による大量死個体群動態モデルの不確実性資源評価の誤り混獲 投棄 密漁 虚偽報告 100 年後の個体群の減少率, 平均漁獲数などから, 管理方式の違いによる影響を評価

シミュレーションを用いる利点 様々な不確実性を取り込んだ上で影響評価が可能 これまでの水産資源評価 管理は, 不確実性の大きさのため, 意思決定が不可能であった. 分からないから何もしない から 分からない中でどうするか へ

シミュレーションモデル N t+1 = [N t + r N t {1-(N t /K) z }]exp(u t ), u t ~N(0,s 2 ) 今後 5 年ごとに調査を行って新しい個体数推定値を得る 翌年から最新の個体数推定値に従って捕獲枠を決定.5 年ごとに見直しをする ( 前年に新しい個体数推定値が得られることを想定 ) 捕獲の方法として, 個体数推定値の 2% をとる場合,PBR でとる場合を考える

PBR Potential Biological Removal PBR = 1/2R max N min F R N min : 最小個体数推定値 ( 下限 20% 点 ) N min = Nexp(-0.84CV) R max : 純増加率 ( 鯨類 0.04, 鰭脚類 0.12) F R : 回復係数 (0.1 for endangered, 0.5 for depleted)

有効なデータ 捕獲統計 個体数推定値

コンディショニング 実際のデータに合うような環境収容力のサンプルを取り出す ( ここでグリッド法を使う ) それから将来に対するシミュレーションを行う 環境収容力 K のサンプルの取り出し方 K にある範囲のグリッドデータを割り当てるそれぞれの K に対して, 既存の個体数推定値に基づく尤度を計算尤度の大きさに従って K をサンプルプログラムを配布するので詳細はそれを見てください

コンディショニング結果

管理方式の評価 2% 捕獲 PBR(FR=1)

頑健性テスト PBR(FR=1) PBR(FR=0.5) 追加分散 = プロセスエラーがある場合の結果

Rejection サンプリング ある確率分布 f(x) から乱数を発生させるのが難しいとき, まず発生させやすい確率分布 g(x) から乱数 x* を発生させて, さらに一様分布 U[0,1] から値 u をとり, u < f(x*)/{cg(x*)} ならば,x* をとっておき, u f(x*)/{cg(x*)} ならば x* を捨てる. 結果, 残った x* は f(x) からのランダムサンプル.

正規分布からのサンプリング cg(x) = 10/60 f(x * ) N(x 1,3) からサンプルすることを考える. x=1 のとき, 高さ 0.133 なので, 一様分布 U[-30,30] からデータ x * を発生させてやると,c = 10 とすれば,10/60 = 0.167 となり, 一様分布 u ~ U[0,1] と N(x * 1,3)/(10/60) を比較することになる. x * ~ U[-30,30]

R プログラム x1 <- runif(100000,-30,30) u <- runif(100000) ac1 <- which(u < dnorm(x1,1,3)/(10/60)) x1 <- x1[ac1] plot(density(x1),col=2,lty=1,lwd=4, xlim=c(-15,15),ylim=c(0,0.15), main="n(1,3)",xlab="x",ylab=" Density") x2 <- seq(-15,15,by=0.1) lines(x2,dnorm(x2,1,3),col=1,lty= 2,lwd=3)

インポータンスサンプリング Q. 対数正規分布 LN(1.1, 0.6) の平均を求めよ. f(x) からのサンプルが難しいとき, サンプルしやすい g(x) からサンプルして, E(x) = (1/n)Σ{f(x*)/g(x*)}x* と計算する. アルゴリズム : この重みをインポータンス比と呼ぶ x1 <- runif(100000,0,60) mean(dlnorm(x1,1.1,0.6)*x1*60) = 3.587 真の平均 exp(1.1+0.6 2 /2) = 3.597

SIR 法 Sampling-Importance Resampling 法 f(x) からのサンプルは困難,g(x) からは容易 1. g(x) からM 個のサンプルx 1*,, x M* を発生 2. w * = f(x*)/g(x*) を計算 3. x 1*,, x M* からw 1*,, w M* に比例する確率でサン プルをm 個とってくる 4. 得られたm 個のサンプルはf(x) からの事後サンプル

R プログラム x1 <- runif(100000,0,60) n2 <- sample(length(x1), 1000, replace=t, prob=dlnorm(x1,1.1,0.6)) x2 <- x1[n2] par(mfrow=c(1,1),mar=c(5,5,5,5)) plot(density(x2),lty=1,lwd=3,col=2,xla b="x",ylab="density",main="lognormal distribution LN(1.1,0.6)",xlim=c(0,20),ylim=c(0,0.3 ),cex.lab=2,cex.main=1.5) x3 <- seq(0,20,by=0.01) lines(x3,dlnorm(x3,1.1,0.6),lty=2,lwd= 2,col=1)

MCMC

マルコフ連鎖とは Pr(x t+1 x t,x t-1,x t-2, ) = Pr(x t+1 x t ) 良い性質を持つマルコフ連鎖は定常分布に収束する ( 生態学の例では,Leslie 行列モデルが安定齢分布を持つ, というのを想像せよ ). この性質を使って, 定常分布が ( サンプルしたいけどしにくい )f(x) になるようにデータを生成するのがマルコフ連鎖モンテカルロ法 (MCMC) である. マルコフ性により得られるサンプルは独立でないことに注意 ( これまでと違うとこ ).

なぜ MCMC? Rejection サンプリング,SIR 法では, 効率良いサンプリングのために, サンプラー g(x) を上手に選んでやる必要がある 非常に複雑な分布 ( 非線形モデル ), たくさんのパラメータがあるとき, 適当な g(x) を見つけるのは難しい ( ときに不可能 ) マルコフ連鎖を使うのは, 暗闇でいきなり動くとタンスに頭をぶつけるが, 手探りで ( 危ない時にはそろそろ, 行ける時には大胆に ) 動けばより安全というイメージ ( 手探りで動くという部分がマルコフ連鎖 )

メトロポリス法 f(x) からサンプルしたい. 1. q(y x)=q(x y) となる分布 ( 対称サンプラーという ) を用意する. たとえば,N(y x,σ 2 ) は対称サンプラー. 2. 初期値 x 0 を決める. 3. q(y x 0 ) からサンプル y を取り出す. 4. 一様分布 U[0,1] から乱数 u を発生. 5. u < f(y)/f(x 0 ) なら,x 1 =y, そうでなければ x 1 =x 0. 6. この操作を繰り返して得た x 1,, x n は f(x) からのサンプルである. x n+1 x n x n+1

N(0,1) のサンプリング B <- 3000 x0 <- -2 e1 <- runif(b,-5,5) u1 <- runif(b,0,1) x1 <- rep(na,b+1) x1[1] <- x0 for (i in 1:B){ xs <- x1[i]+e1[i] x1[i+1] <- ifelse(u1[i] < dnorm(xs)/dnorm(x1[i]), xs, x1[i]) } plot 部分は省略

初期値の影響

メトロポリス ヘイスティングス (MH) 法 対称サンプラーでない場合を考えてみる. x q(x y) q(y x) y もし q(x y) > q(y x) ならば,x は y よりもずっとサンプルされやすくなる. これでは正しい f(x) からのサンプルにならない. インポータンスサンプリングと同様に重み付き分布の考え方を使って, 事後分布をインポータンス比で置き換えてやる. つまり,f(x)/q(x y) である.

MH 法のアルゴリズム f(x) からサンプルしたい. 1. ある条件付分布 q(y x) を用意する. 2. 初期値 x 0 を決める. 3. q(y x 0 ) からサンプル y を取り出す. 4. 一様分布 U[0,1] から乱数 u を発生. 5. u < {f(y)/q(y x 0 )}/{f(x 0 )/q(x 0 y)} = {f(y)q(x 0 y)}/{f(x 0 )q(y x 0 )} なら,x 1 = y, そうでなければ x 1 =x 0. 6. この操作を繰り返して得た x 1,, x n は f(x) からのサンプルである.

t(15) のサンプリング 独立サンプラー N(0,1.5 2 ) を使って,t 分布 t(15) からサンプルする. x0 <- -1 s1 <- 1.5 B <- 3000 x1 <- rep(na,(b+1)) x1[1] <- x0 y <- rnorm(b,0,s1) u1 <- runif(b) for (i in 1:B){ a1 <- min(1,dt(y[i],15)*dnorm(x1[i],0,s1)/(dt(x1[i], 15)*dnorm(y[i],0,s1))) x1[i+1] <- ifelse(u1[i] < a1, y[i], x1[i]) }

Single-component MH 法 多変量 f(x 1,x 2, ) の場合 ひとつの変数に対して順番にMH 法を繰り返す 1. x (0), y (0) を決定 2. f(x y (0) ) からのサンプルをMH 法で,x (1) =x * or x (0) 3. f(y x (1) ) からのサンプルをMH 法で,y (1) =y * or y (0) 4. この操作を繰り返す ここで,f(x y) = f(x,y)/ f(x,y)dx ( フル条件付分布 )

なぜうまくいくか ( 直感的に ) y x 4 y x 7 x y 1 x y 2 x y 6 禁ななめ移動

例題 データ x = (1,2,3) を得たとき, 正規分布 N(m,s 2 ) の事後分布は, 事前分布をP(m,s) 1/sとして, P(m,s x) (1/s 4 )exp{-{s 2 +3(m-mx) 2 }/(2s 2 )} となる ( ここで,mx = (1+2+3)/3,S 2 =Σ(x i -mx) 2 ). この P(m,s (1,2,3)) からサンプリングせよ

R プログラム x0 <- -1; y0 <- 1; B <- 3000 e1 <- runif(b,-7,7); u1 <- runif(b,0,1) ys <- rgamma(b,1,1); u2 <- runif(b,0,1) x1 <- y1 <- rep(na,b+1); x1[1] <- x0; y1[1] <- y0 xobs <- c(1,2,3) mx <- mean(xobs); S2 <- sum((xobs-mx)^2) n <- length(xobs); nu <- n-1 MH 法の場合, 分母の積分値がいらないのに注意 post1 <- function(m,s1) (1/s1^(n+1))*exp(-(S2+n*(m-mx)^2)/(2*s1^2)) for (i in 1:B){ xs <- x1[i]+e1[i] x1[i+1] <- ifelse(u1[i] < post1(xs,y1[i])/post1(x1[i],y1[i]), xs, x1[i]) y1[i+1] <- ifelse(u2[i] < post1(x1[i+1],ys[i])*dgamma(y1[i],2,0.5)/(post1(x1[i+1],y1[i])*dgamm a(ys[i],2,0.5)), ys[i], y1[i]) }

事後分布 実は, muの周辺事後分布は, t(mx, S 2 /n, n-1) (tは一般化 t 分布 ) s 2 の周辺事後分布は, IG(n/2, S 2 /2) (IGは逆ガンマ分布) に従うことが分かっている.

ギブスサンプリング Single-component MH 法でサンプラーをフル条件付分布にする. MH 法の採択確率は,f(x)q(y x)/f(y)q(x y) であるから, サンプルしたい目的の関数 f(x) とサンプルしやすい関数 q(x y) をどちらも同じフル条件付分布にすることになる. よって, 採択確率 =1 効率が良い

ギブスサンプリング 多変量関数 f(x 1,x 2,x 3, ) からサンプルひとつの変数に対して順番にフル条件付分布からサンプルする 1. x (0) 1, x (0) 2,x (0) 3, を決定 2. f(x 1 x (0) 2,x (0), 3 ) からサンプルx (1) 1 を得る. 3. f(x 2 x (1) 1,x (0), 3 ) からサンプルx (1) 2 を得る. 4. この操作を繰り返す

ギブスサンプリング フル条件付分布からのサンプリングが簡単な場合とても魅力的 フル条件付分布からのサンプリングが簡単ではない場合,Rejection サンプリングや MH 法と組み合わせる (MH 法と組み合わせたとき, MH within Gibbs と呼んだりする ) 明日使う WinBUGS は基本的に Gibbs Sampler を用いて事後分布からサンプルするソフトである

例題 2 変量正規分布 x x 1 2 1 ~ N, 2 1 0.7 0.7 1 からデータをサンプルする. この場合, x 1 x 2 ~ N(1+0.7(x 2-2), 1-0.7 2 ) x 2 x 1 ~ N(2+0.7(x 1-1), 1-0.7 2 ) となることが分かっている. フル条件付分布からのサンプルは通常の正規分布からのサンプル

R プログラム B <- 2000 x1 <- x2 <- rep(na,b+1) x1[1] <- -2; x2[1] <- 1 for (i in 1:B){ x1[i+1] <- rnorm(1,1+0.7*(x2[i]- 2),sqrt(1-0.7^2)) x2[i+1] <- rnorm(1,2+0.7*(x1[i+1]- 1),sqrt(1-0.7^2)) }

スライスサンプリング X~f(x) (X,V)~U[(x,v):0 v f(x)] これを使って, 1. u (0), x (0) を決めてやる 2. u (1) ~ U[0, f(x (0) )] 3. x (1) ~ U[x: f(x) u (1) ] 4. この操作を繰り返す (f(x) の代わりにf 1 (x)=cf(x) を用いてもOK)

N(0,1) からのサンプリング u (t+1) N(0,1) からのサンプリングを考える. これは exp(-x 2 /2) からのサンプリングを考えても同じことである. この場合, exp(-x 2 /2) u は,x 2-2log(u) となる. これが slice x (t) x (t+1)

R プログラム B <- 3000 x1 <- u1 <- rep(na,b+1) x1[1] <- 0.5 u1[1] <- 0.5 for (i in 1:B){ u1[i+1] <- runif(1,0,dnorm(x1[i],0,1)*sqrt(2*pi)) x1[i+1] <- runif(1,-sqrt(- 2*log(u1[i+1])),sqrt(-2*log(u1[i+1]))) }

Burn-in と Thinning Burn-in: 初期値が目的とする分布のはじっこの方にあるとき, 最初の方のサンプルはまだ定常分布に収束していないので捨ててやる必要がある Thinning: MCMC の結果得られたサンプルは互いに相関しており, 独立なサンプルとなっていない. 独立サンプルを得るためのひとつの (ad hoc) 方法として, ある間隔ごとにデータをとる. コンピュータメモリの節約にも.

初期値の影響 自己相関 初期値の影響が大きい 自己相関している

収束診断 MCMC 時系列が目的とする分布に収束している必要がある. 収束したかどうかの判定にはいくつかの方法が考えられている. よく使われるのは, 数本のMCMCを走らせて, それらがどれも同じようなところに分布するかどうかを見る方法である ( 詳細は明日のWinBUGSの話で見る ).MCMCではとにかく色々とプロットしてみることが重要である ( そんなときRが便利 ).

MCMC 注意 何でも OK な万能薬ではない パラメータが多いと計算時間がかかる パラメータ間の相関が大きいと計算時間がかかる log をとったりすると計算時間がかかる 初期値を変えて同じ結果になるか確かめろ, 計算時間がかかるけど うまくプログラムしないと計算時間がやたらとかかる

参考文献 平松一彦 2004. MCMC 入門. 水産資源管理談話会報 34: 72-76. 大森裕浩 2001. マルコフ連鎖モンテカルロ法の最近の展開. 日本統計学会誌 31: 305-344. 伊庭幸人 2005. マルコフ連鎖モンテカルロ法の基礎. 統計科学のフロンティア計算統計 II : 1-106.