Markov連鎖モンテカルロ法

同義語:マルコフ連鎖モンテカルロ法Markov chain Monte CarloMCMC

概要

Markov連鎖モンテカルロ法(Markov chain Monte Carlo、MCMC)とは、正規化の定数が分からない分布 $\pi$ について、$\pi$ を定常分布にもつ Markov 連鎖を走らせ、その時間平均で期待値 $E_\pi[h]$ を見積もる計算法である。Metropolis–Hastings 法は提案した候補を確率 $\min\bigl(1,\frac{\pi(y)Q(y,x)}{\pi(x)Q(x,y)}\bigr)$ で受理し、比 $\pi(y)/\pi(x)$ だけを使って詳細釣り合いを満たす連鎖を作る。状態が有限個で $\pi>0$ なら、提案が既約で支えが対称なとき時間平均は確率 $1$ で $E_\pi[h]$ に収束し、自己ループがあれば $X_n$ の分布も指数的な速さで $\pi$ に近づく。Gibbs サンプリングはその特別な場合である。

$$\newcommand{C}[0]{\mathbb{C}} \newcommand{div}[0]{\mathbin{÷}} \newcommand{N}[0]{\mathbb{N}} \newcommand{Q}[0]{\mathbb{Q}} \newcommand{R}[0]{\mathbb{R}} \newcommand{Z}[0]{\mathbb{Z}} $$

前提知識: Markov連鎖, Monte Carlo法, 確率漸化式と定常分布, エルゴード定理

Monte Carlo法 では、分布 $\pi$ に従う独立な標本 $Y_1,\dots,Y_n$ を作り、期待値 $E_\pi[h]=\sum_x\pi(x)h(x)$ を標本平均で見積もる。ところが、統計物理の Ising 模型($N$ 個のスピンの配置 $2^N$ 通りの上の分布)や Bayes 統計の事後分布のように、$\pi(x)=\tilde\pi(x)/Z$ の 比 $\tilde\pi(y)/\tilde\pi(x)$ は簡単に計算できても、正規化の定数 $Z=\sum_x\tilde\pi(x)$ が膨大な和で計算できず、$\pi$ から直接に標本を作れないことが多い。
Markov連鎖モンテカルロ法(MCMC)は、$\pi$ を定常分布にもつ Markov連鎖 $X_0,X_1,X_2,\dots$ を作って走らせ、独立な標本の代わりに連鎖の値の時間平均
$$ \frac1n\sum_{m=0}^{n-1}h(X_m) $$
で $E_\pi[h]$ を見積もる方法である。連鎖の作り方には $Z$ が要らない。この記事では、状態の集合 $S$ が有限で、すべての $x\in S$ で $\pi(x)>0$ の場合に、連鎖を作る 2 つの標準的な方法(Metropolis–Hastings 法 と Gibbs サンプリング)を定義し、それらが $\pi$ を定常分布にもつこと(thm-mcmc-mh・prop-mcmc-gibbs)と、時間平均が $E_\pi[h]$ に収束し、さらに $X_n$ の分布が $\pi$ に指数的な速さで近づくこと(thm-mcmc-main)を証明する。

定義と使える条件

$S$ を有限集合、$\pi$ を $S$ 上の分布とする。$S$ 上の推移行列 $P=(P(x,y))_{x,y\in S}$(各行の成分が $0$ 以上で和が $1$ の行列)について、分布を行ベクトルとみて $\pi P=\pi$ となるとき、$\pi$ を $P$ の 定常分布 という(Markov連鎖 の定義「定常分布」)。

詳細釣り合い

推移行列 $P$ と分布 $\pi$ が
$$ \pi(x)P(x,y)=\pi(y)P(y,x)\qquad(x,y\in S) $$
を満たすとき、$P$ は $\pi$ について 詳細釣り合い(detailed balance)を満たす、または $\pi$ について 可逆 であるという。

詳細釣り合いなら定常

$P$ が $\pi$ について詳細釣り合いを満たすならば、$\pi$ は $P$ の定常分布である。

各 $y$ について $(\pi P)(y)=\sum_x\pi(x)P(x,y)=\sum_x\pi(y)P(y,x)=\pi(y)\sum_xP(y,x)=\pi(y)$ である。$\square$

詳細釣り合いは、$\pi$ に従って選んだ状態から 1 歩進むとき、「$x$ から $y$ へ」と「$y$ から $x$ へ」の確率の流れが等しいことを言っている。定常であるには、各状態について入る流れの総和と出る流れの総和がつり合えばよいので、詳細釣り合いは定常であるための十分条件であって、必要条件ではない(ex-mcmc-gibbs の巡回走査)。

Metropolis–Hastings 法

すべての $x\in S$ で $\pi(x)>0$ とし、$S$ 上の推移行列 $Q$(提案行列)を 1 つとる。$x\ne y$ で $Q(x,y)>0$ のとき
$$ \alpha(x,y):=\min\Bigl(1,\ \frac{\pi(y)\,Q(y,x)}{\pi(x)\,Q(x,y)}\Bigr) $$
(受理確率)とおき、推移行列 $P$ を
$$ P(x,y):=\begin{cases}Q(x,y)\,\alpha(x,y)&(x\ne y,\ Q(x,y)>0),\\ 0&(x\ne y,\ Q(x,y)=0),\\ 1-\sum_{z\ne x}P(x,z)&(y=x)\end{cases} $$
と定める。$P$ を推移行列とする Markov 連鎖を作ることを Metropolis–Hastings 法 という。$Q$ が対称($Q(x,y)=Q(y,x)$)のときは $\alpha(x,y)=\min(1,\pi(y)/\pi(x))$ で、これを Metropolis 法 という。

連鎖を 1 歩進めるには、現在の状態 $x$ から $Q(x,\cdot)$ に従って候補 $y$ を選び、確率 $\alpha(x,y)$ で $y$ へ移り(受理)、そうでなければ $x$ にとどまる(棄却)。これは上の $P$ に従って 1 歩進むことと同じである。受理確率には比 $\pi(y)/\pi(x)=\tilde\pi(y)/\tilde\pi(x)$ しか現れないので、$Z$ は要らない(Mac03 §29.4、pp. 365–366 の (29.31))。
使える条件は次のとおりである。

  1. $S$ は有限集合で、すべての $x\in S$ で $\pi(x)>0$ である($\alpha$ の分母が $0$ にならない)。
  2. 提案の支えが対称である:$Q(x,y)>0$ と $Q(y,x)>0$ は同値(ex-mcmc-counterexamples の 2 行目)。
  3. $Q$ は 既約 である:どの $x,y$ についても、$Q$ の正の成分をたどって $x$ から $y$ へ行ける(ex-mcmc-counterexamples の 1 行目)。
  4. $X_n$ の分布の収束(thm-mcmc-main の 3)では、さらに $P(z,z)>0$ となる状態 $z$ が 1 つあることを仮定する。これは十分条件で、必要条件ではない($S=\{1,2,3\}$、$\pi$ が一様で、ほかの 2 つの状態のどちらかを確率 $\frac12$ ずつ選ぶ提案なら $P(z,z)=0$ だが、$P^2$ の成分はすべて正である)。この条件を外すと収束しない例がある(ex-mcmc-counterexamples の 4 行目)。

公式の表

方法提案受理確率 $\alpha(x,y)$使える条件根拠
Metropolis–Hastings 法一般の $Q$$\min\bigl(1,\frac{\pi(y)Q(y,x)}{\pi(x)Q(x,y)}\bigr)$条件 1–3thm-mcmc-mh
Metropolis 法対称な $Q$$\min\bigl(1,\frac{\pi(y)}{\pi(x)}\bigr)$条件 1–3thm-mcmc-mh
独立型$Q(x,y)=q(y)$$\min\bigl(1,\frac{\pi(y)q(x)}{\pi(x)q(y)}\bigr)$$q>0$ が $S$ の全体でthm-mcmc-mh
Barker 法一般の $Q$$\frac{\pi(y)Q(y,x)}{\pi(x)Q(x,y)+\pi(y)Q(y,x)}$条件 1–3下の注意
Gibbs サンプリング(座標 $i$ の更新)条件付き分布 $\pi(\cdot\mid x_{-i})$$1$直積の集合、$\pi>0$prop-mcmc-gibbs
詳細釣り合い ⇒ 定常———prop-mcmc-db
定常な行列の積・平均$P_1,P_2$ が $\pi$ を定常にもつ——$\pi P_1P_2=\pi$、$\pi\frac{P_1+P_2}2=\pi$

独立型は $Q(y,x)/Q(x,y)=q(x)/q(y)$ を代入したものである。Barker 法の受理確率 $\beta(x,y)$ では、$\pi(x)Q(x,y)\beta(x,y)=\frac{\pi(x)Q(x,y)\,\pi(y)Q(y,x)}{\pi(x)Q(x,y)+\pi(y)Q(y,x)}$ が $x,y$ について対称なので、詳細釣り合いが成り立つ。

主定理と証明

Metropolis–Hastings 法の詳細釣り合い

def-mcmc-mh の $P$ は推移行列であり、$\pi$ について詳細釣り合いを満たす。とくに $\pi$ は $P$ の定常分布である。

$x\ne y$ について $0\le P(x,y)\le Q(x,y)$ なので、$\sum_{z\ne x}P(x,z)\le\sum_{z\ne x}Q(x,z)\le1$ であり、$P(x,x)\ge0$ である。各行の和は定義から $1$ である。よって $P$ は推移行列である。
詳細釣り合い:$x=y$ では自明なので $x\ne y$ とする。$Q(x,y)>0$ のとき
$$ \pi(x)P(x,y)=\pi(x)Q(x,y)\min\Bigl(1,\frac{\pi(y)Q(y,x)}{\pi(x)Q(x,y)}\Bigr)=\min\bigl(\pi(x)Q(x,y),\ \pi(y)Q(y,x)\bigr) $$
である。$Q(x,y)=0$ のときは $\pi(x)P(x,y)=0$ で、これも $\min\bigl(\pi(x)Q(x,y),\pi(y)Q(y,x)\bigr)=0$ に等しい。どちらの場合も右辺は $x$ と $y$ について対称なので、$\pi(x)P(x,y)=\pi(y)P(y,x)$ である。最後の主張は prop-mcmc-db による。$\square$

この証明は、使える条件の 2 と 3 を使っていない。どんな提案行列でも $\pi$ は定常になるが、それだけでは $\pi$ 以外の定常分布もありうる(ex-mcmc-counterexamples)。次の 2 つの補題で、連鎖が $\pi$ に収束するための条件を確かめる。

既約性の受け継ぎ

使える条件 1–3 のもとで、def-mcmc-mh の $P$ は既約である。

$x\ne y$ で $Q(x,y)>0$ なら、条件 2 により $Q(y,x)>0$ で、$\pi(y)Q(y,x)>0$ だから $\alpha(x,y)>0$、よって $P(x,y)>0$ である。すなわち $Q$ の正の成分(対角成分を除く)は $P$ でも正である。条件 3 により $x$ から $y$ への $Q$ の正の成分の道があり、対角成分を通る段を除いても道のままなので、それは $P$ の正の成分の道でもある。$\square$

1 か所の自己ループ

$S$ を $d$ 個の元からなる有限集合、$P$ を既約な推移行列とし、ある $z\in S$ で $P(z,z)>0$ とする。このとき $k:=2(d-1)$($d=1$ なら $k:=1$)とおくと、$P^k$ のすべての成分は正である。

$d=1$ なら $P=(1)$ である。$d\ge2$ とし、$x,y\in S$ を任意にとる。既約性により $x$ から $z$ への正の成分の道があり、同じ状態を 2 度通らないものを選べば、その長さ $\ell_1$ は $d-1$ 以下である。同様に $z$ から $y$ への長さ $\ell_2\le d-1$ の道がある($x=z$ なら $\ell_1=0$ とし、$z=y$ なら $\ell_2=0$ とする)。$j:=k-\ell_1-\ell_2\ge0$ 回だけ $z$ にとどまる段をはさむと、長さちょうど $k$ の $x$ から $y$ への道ができ、その確率は $P(z,z)^j>0$ を含む正の数の積である。$P^k(x,y)$ はこの道の確率以上なので正である。$\square$

Metropolis–Hastings 法では、$Q(z,z)>0$ となる $z$ があるか、棄却が正の確率で起こる組($Q(z,y)>0$ かつ $\pi(y)Q(y,z)<\pi(z)Q(z,y)$)があれば、$P(z,z)\ge Q(z,z)+Q(z,y)\bigl(1-\alpha(z,y)\bigr)>0$ となって、この補題の仮定が満たされる。

Markov 連鎖モンテカルロ法の収束

$S$ を $d$ 個の元からなる有限集合、$\pi$ を $S$ 上の分布で $\pi>0$、$P$ を $\pi$ を定常分布にもつ既約な推移行列とし、$(X_n)$ を推移行列 $P$ の Markov 連鎖(初期分布は任意)とする。$h\colon S\to\mathbb R$ とし、$\|h\|_\infty:=\max_x\lvert h(x)\rvert$ とおく。

  1. $\pi$ は $P$ のただ 1 つの定常分布である。
  2. ほとんど確実に $\displaystyle\lim_{n\to\infty}\frac1n\sum_{m=0}^{n-1}h(X_m)=E_\pi[h]$ である。
  3. さらにある $z$ で $P(z,z)>0$ ならば、$k$ を lem-mcmc-aperiodic のとおりとし、$\delta:=\min_{x,y}P^k(x,y)>0$ とおくと、すべての $n$ で
    $$ \bigl\lvert E[h(X_n)]-E_\pi[h]\bigr\rvert\le2\,\|h\|_\infty\,(1-d\,\delta)^{\lfloor n/k\rfloor} $$
    である。

1 は Markov連鎖 の定理「既約な有限連鎖の定常分布」による。2 は、既約な有限連鎖について時間平均が定常分布での期待値に収束することを述べた エルゴード定理 の定理「有限 Markov 連鎖のエルゴード定理」そのものである(初期分布は任意でよく、周期的でもよい)。
3:lem-mcmc-aperiodic により $P^k$ のすべての成分は $\delta>0$ 以上で、$P$ は 確率漸化式と定常分布 の定義「正則な確率行列」の意味で正則である。$\mu$ を $X_0$ の分布とすると $X_n$ の分布は $\mu P^n$ であり、同じ記事の定理「正則な連鎖の収束」により、1 で示した唯一の定常分布 $\pi$ について
$$ \sum_{y\in S}\bigl\lvert(\mu P^n)(y)-\pi(y)\bigr\rvert\le2(1-d\,\delta)^{\lfloor n/k\rfloor} $$
である。したがって
$$ \bigl\lvert E[h(X_n)]-E_\pi[h]\bigr\rvert=\Bigl\lvert\sum_y\bigl((\mu P^n)(y)-\pi(y)\bigr)h(y)\Bigr\rvert\le\|h\|_\infty\sum_y\bigl\lvert(\mu P^n)(y)-\pi(y)\bigr\rvert $$
から主張が従う。$\square$

Metropolis–Hastings 法では、thm-mcmc-mh と lem-mcmc-irreducible により、使える条件 1–3 のもとで thm-mcmc-main の 1・2 が成り立ち、条件 4 も満たせば 3 も成り立つ。2 は、時間平均の 一致性(標本の数を増やすと正しい値に近づくこと)を保証する。3 は、最初の $b$ 歩を捨ててから平均をとる バーンイン(burn-in)の考え方の根拠である。ただし評価の係数 $1-d\delta$ は $d$ が大きいとき $1$ にとても近くなりうる粗い上界で、実際の収束は評価よりずっと速いことがある(ex-mcmc-three)。

Gibbs サンプリング

状態が多変数 $x=(x_1,\dots,x_r)$ で、全体の分布からは標本を作れないが、1 つの変数の条件付き分布からは作れる、という場合に使う方法である。$S=S_1\times\cdots\times S_r$(各 $S_i$ は有限集合)とし、すべての $x\in S$ で $\pi(x)>0$ とする。$x$ から第 $i$ 成分を除いたものを $x_{-i}$ と書き、$\pi_{-i}(x_{-i}):=\sum_{z\in S,\ z_{-i}=x_{-i}}\pi(z)$ とおく。

Gibbs サンプリング

$i=1,\dots,r$ について、推移行列 $G_i$ を
$$ G_i(x,y):=\begin{cases}\dfrac{\pi(y)}{\pi_{-i}(x_{-i})}&(y_{-i}=x_{-i}),\\[2mm] 0&(\text{それ以外})\end{cases} $$
と定める。$G_i$ は、第 $i$ 成分だけを、ほかの成分を固定した条件付き分布 $\pi(\cdot\mid x_{-i})$ に従って引き直す操作である。成分を順に更新する $G:=G_1G_2\cdots G_r$ を 巡回走査、更新する成分を一様に選ぶ $R:=\frac1r\sum_{i=1}^rG_i$ を 無作為走査 の Gibbs サンプリング という。

Gibbs サンプリングの性質

上の設定で次が成り立つ。

  1. 各 $G_i$ は推移行列で、$\pi$ について詳細釣り合いを満たす。$G_i$ は、提案行列を $G_i$ とした Metropolis–Hastings 法で受理確率がつねに $1$ のものに等しい。
  2. 無作為走査 $R$ は $\pi$ について詳細釣り合いを満たし、既約で、すべての $x$ で $R(x,x)>0$ である。
  3. 巡回走査 $G$ は $\pi$ を定常分布にもち、$G$ のすべての成分は正である。
    したがって $R$ と $G$ のどちらにも thm-mcmc-main の 1–3 が当てはまる。

1:$\sum_yG_i(x,y)=\sum_{y_{-i}=x_{-i}}\pi(y)/\pi_{-i}(x_{-i})=1$ である。$y_{-i}=x_{-i}$ なら $\pi_{-i}(x_{-i})=\pi_{-i}(y_{-i})$ なので $\pi(x)G_i(x,y)=\pi(x)\pi(y)/\pi_{-i}(x_{-i})$ は $x,y$ について対称であり、$y_{-i}\ne x_{-i}$ なら両辺とも $0$ である。$G_i(x,y)>0$ のとき $\frac{\pi(y)G_i(y,x)}{\pi(x)G_i(x,y)}=\frac{\pi(y)\pi(x)}{\pi(x)\pi(y)}=1$ なので、受理確率は $1$ であり、def-mcmc-mh の $P$ は $G_i$ そのものになる。
2:詳細釣り合いの式は $P$ について 1 次なので、$G_i$ の平均 $R$ も満たす。$u$ と $v$ が第 $i$ 成分だけ違うとき $R(u,v)\ge\frac1rG_i(u,v)=\frac{\pi(v)}{r\,\pi_{-i}(u_{-i})}>0$ であり、$x$ から $y$ へは成分を 1 つずつ $y$ の値に変えていけば行けるので、$R$ は既約である。$R(x,x)\ge\frac1rG_1(x,x)>0$ である。
3:prop-mcmc-db により $\pi G_i=\pi$ なので $\pi G=\pi G_1G_2\cdots G_r=\pi$ である。$x,y\in S$ について、$x^{(0)}:=x$、$x^{(i)}:=(y_1,\dots,y_i,x_{i+1},\dots,x_r)$ とおくと、$x^{(i)}$ と $x^{(i-1)}$ は第 $i$ 成分以外が等しいので
$$ G(x,y)\ge\prod_{i=1}^rG_i\bigl(x^{(i-1)},x^{(i)}\bigr)>0 $$
である($x^{(r)}=y$)。成分がすべて正なので既約で、$G(x,x)>0$ である。$\square$

巡回走査 $G$ は一般に詳細釣り合いを満たさないが、$\pi$ は定常である(ex-mcmc-gibbs)。Gibbs サンプリングは、変数どうしの相関が強いと 1 成分ずつの更新で状態があまり動かず、収束が遅くなる(ex-mcmc-gibbs の後半)。Gibbs サンプリングを Metropolis 法の一種とみる見方と、詳細釣り合いの用語は Mac03 §29.5–29.6(pp. 370–374。詳細釣り合いは p. 374 の (29.46))にある。

計算の手順

Metropolis–Hastings 法で $E_\pi[h]$ を見積もる手順は次のとおりである。

  1. 目標の分布を、定数倍を除いた形 $\tilde\pi$ で書く。すべての状態で $\tilde\pi>0$ であることを確かめる(使える条件 1)。
  2. 提案行列 $Q$ を選ぶ。近くの状態へ動く対称な提案(Metropolis 法)が基本である。支えが対称で(条件 2)、既約である(条件 3)ことを確かめる。
  3. 初期値 $X_0$ を選び、$m=0,1,2,\dots$ について次を繰り返す:$Q(X_m,\cdot)$ に従って候補 $y$ を引き、$[0,1]$ の一様乱数 $U$ を引いて、$U\le\alpha(X_m,y)$ なら $X_{m+1}:=y$、そうでなければ $X_{m+1}:=X_m$ とする。棄却したときも $X_m$ をもう一度数える。
  4. 最初の $b$ 歩を捨て、$\hat h:=\frac1n\sum_{m=b}^{b+n-1}h(X_m)$ を見積もりとする(thm-mcmc-main の 2 により、$b$ を固定して $n\to\infty$ とすれば $\hat h\to E_\pi[h]$)。
  5. 結果を確かめる。受理された割合が $0$ や $1$ に近すぎないか、異なる初期値から始めた複数の連鎖が同じ値に近づくかを見る。連鎖の値は互いに独立でないので、独立な標本の公式 $\sigma/\sqrt n$ で誤差を見積もってはいけない(Monte Carlo法 の定理「推定量の基本性質」は独立な標本についての主張である)。
    Gibbs サンプリングでは、3 の代わりに、成分を 1 つ選んで(または順に)その条件付き分布から引き直す。

例題

3 状態の Metropolis 法

$S=\{1,2,3\}$、$\tilde\pi=(1,2,3)$($\pi=(\frac16,\frac13,\frac12)$)とし、提案は「ほかの 2 つの状態のどちらかを確率 $\frac12$ ずつ」とする。
使う公式:def-mcmc-mh、thm-mcmc-main。
答え:推移行列は
$$ P=\begin{pmatrix}0&\frac12&\frac12\\ \frac14&\frac14&\frac12\\ \frac16&\frac13&\frac12\end{pmatrix} $$
で、$\pi P=\pi$ である。$X_0=1$ から始めると $X_n$ の分布は $(0,\frac12,\frac12)$、$(\frac5{24},\frac7{24},\frac12)$、$(\frac5{32},\frac{11}{32},\frac12)$ と進み、$\pi$ との差の絶対値の和は $\frac13,\frac1{12},\frac1{48},\dots$ と毎回ちょうど $\frac14$ 倍になる($n=6$ で約 $3.26\times10^{-4}$)。

計算と評価の比較を開く

$1\to2$ の受理確率は $\min(1,\frac21)=1$ で $P(1,2)=\frac12$、$2\to1$ は $\min(1,\frac12)$ で $P(2,1)=\frac14$、$3\to1$ は $\frac13$ で $P(3,1)=\frac16$、$3\to2$ は $\frac23$ で $P(3,2)=\frac13$ である。対角成分は行の和が $1$ になるように決まる。$P$ の固有値は $1,-\frac14,0$ で、差が毎回 $\frac14$ 倍になるのは固有値 $-\frac14$ による。thm-mcmc-main の 3 の証明で使った 確率漸化式と定常分布 の定理「正則な連鎖の収束」を直接当てると、$P^2$ の成分がすべて正なので(最小は $\frac7{48}$)、$k=2$、$1-3\delta=\frac9{16}$ として差は $2(\frac9{16})^{\lfloor n/2\rfloor}$ 以下と評価される。実際の差はこの上界よりずっと小さい。

2 項分布を目標にした酔歩

$S=\{0,1,\dots,10\}$、$\tilde\pi(x)=\binom{10}x$ とし、提案は $x\pm1$ を確率 $\frac12$ ずつ($S$ の外に出る候補は棄却としてその場にとどまる、すなわち $Q(0,0)=Q(10,10)=\frac12$)とする。$E_\pi[X]=5$、$E_\pi[X^2]=27.5$ である。
使う公式:def-mcmc-mh(Metropolis 法)、thm-mcmc-main。
答え:$x\to x+1$ の受理確率は $\min\bigl(1,\frac{10-x}{x+1}\bigr)$、$x\to x-1$ は $\min\bigl(1,\frac x{11-x}\bigr)$ である。たとえば $P(3,4)=\frac12$、$P(3,2)=\frac3{16}$、$P(3,3)=\frac5{16}$ である。$X_0=0$ から始めると、$E[X_n]$ は $n=10$ で約 $3.6702$、$n=50$ で約 $4.99769$ であり、$n=100$ で分布と $\pi$ の差の絶対値の和は約 $3.7\times10^{-7}$ になる。乱数で連鎖を走らせ、最初の $1000$ 歩を捨てて $n=10^5$ 歩の平均をとった 1 回の実行では、$X$ の平均は約 $5.004$、$X^2$ の平均は約 $27.53$、受理の割合は約 $0.75$ であった。

Gibbs サンプリングの 2 つの走査

$S=\{0,1\}^2$ 上の $\pi(0,0)=\pi(1,1)=0.4$、$\pi(0,1)=\pi(1,0)=0.1$ とする。
使う公式:def-mcmc-gibbs、prop-mcmc-gibbs。
答え:第 1 成分の更新 $G_1$ は、$x_2=0$ のとき $x_1=0$ を確率 $0.8$、$x_2=1$ のとき $x_1=1$ を確率 $0.8$ で選ぶ。状態を $(0,0),(0,1),(1,0),(1,1)$ の順に並べると、巡回走査 $G=G_1G_2$ の第 1 行は $(\frac{16}{25},\frac4{25},\frac1{25},\frac4{25})$ である。$\pi G=\pi$ だが、
$$ \pi(0,0)\,G\bigl((0,0),(0,1)\bigr)=0.4\cdot\tfrac4{25}=\tfrac8{125}\ne\tfrac2{125}=0.1\cdot\tfrac4{25}=\pi(0,1)\,G\bigl((0,1),(0,0)\bigr) $$
なので詳細釣り合いは成り立たない。無作為走査 $R=\frac12(G_1+G_2)$ は詳細釣り合いを満たす。

相関と収束の速さを開く

$G$ の固有値は $1,0.36,0,0$ で、分布の $\pi$ との差はおよそ毎回 $0.36$ 倍になる。$\pi(0,0)=\pi(1,1)=c$、$\pi(0,1)=\pi(1,0)=\frac12-c$ として相関を強めると、$G$ の $1$ 以外の固有値は $c=0.45$ で $0.64$、$c=0.49$ で $0.9216$ と $1$ に近づき、収束が遅くなる。$0.36=0.6^2$、$0.64=0.8^2$、$0.9216=0.96^2$ で、$0.6,0.8,0.96$ はそれぞれ $x_1$ と $x_2$ の相関係数 $4c-1$ である。

使える条件と反例

外す条件反例成り立たなくなること
$Q$ が既約(条件 3)$S=\{1,2,3,4\}$、$\pi$ は一様、$Q$ は同じ組($\{1,2\}$ または $\{3,4\}$)の 2 状態のどちらかを確率 $\frac12$ ずつ選ぶ定常分布がただ 1 つであること、時間平均が $E_\pi[h]$ に収束すること
提案の支えが対称(条件 2)$S=\{1,2,3\}$、$\pi$ は一様、$Q$ は $1\to2\to3\to1$ と巡回させる$P$ が既約であること($P=I$ になる)
受理確率の $Q(y,x)/Q(x,y)$ の因子$S=\{1,2\}$、$\pi=(\frac12,\frac12)$、$Q(1,2)=1$、$Q(2,1)=Q(2,2)=\frac12$ で $\alpha=\min(1,\pi(y)/\pi(x))$ を使う$\pi$ が定常分布であること
$P(z,z)>0$ となる $z$(条件 4)$S=\{1,2\}$、$\pi=(\frac12,\frac12)$、$Q$ は必ず入れ替える$X_n$ の分布が $\pi$ に収束すること
反例の確かめ
4 つの確かめを開く

1 行目:$Q$ は対称で支えも対称で、$Q(z,z)=\frac12>0$ だが、既約でない。受理確率はすべて $1$ で $P=Q$ である。$\pi$ は定常だが、$\{1,2\}$ の上の一様分布も定常で、定常分布は 1 つでない。$X_0=1$ なら連鎖は $\{1,2\}$ から出ないので、$h=\mathbf 1_{\{3,4\}}$ の時間平均はつねに $0$ で、$E_\pi[h]=\frac12$ にならない。

2 行目:$Q$ は既約だが、$Q(1,2)>0$ で $Q(2,1)=0$ である。$P(x,y)=\min\bigl(Q(x,y),\pi(y)Q(y,x)/\pi(x)\bigr)$ が $x\ne y$ ですべて $0$ になり、$P=I$ で、連鎖は初期値から動かない。

3 行目:$\pi(1)=\pi(2)$ なので補正なしの受理確率はつねに $1$ で、$P=Q$ となる。$Q$ の定常分布は $(\frac13,\frac23)$ で、$\pi$ ではない。正しい受理確率では $\alpha(1,2)=\min(1,\frac{(1/2)(1/2)}{(1/2)\cdot1})=\frac12$、$\alpha(2,1)=1$ で、$P=\begin{pmatrix}\frac12&\frac12\\ \frac12&\frac12\end{pmatrix}$ となり、$\pi$ が定常である。

4 行目:受理確率はつねに $1$ で $P=\begin{pmatrix}0&1\\1&0\end{pmatrix}$ であり、既約で $\pi$ は定常だが、$P(z,z)=0$ である。$X_0=1$ なら $X_n$ の分布は $(1,0)$ と $(0,1)$ を交互にとり、収束しない。時間平均は $E_\pi[h]$ に収束する(thm-mcmc-main の 2 は周期的でも成り立つ)。

よくある誤り

誤り正しくはなぜ
棄却した候補を捨て、受理した状態だけを数える棄却したら現在の状態をもう一度数える受理した状態だけの列は $P$ の連鎖でなく、$\pi$ に従わない
非対称な提案で $\min(1,\pi(y)/\pi(x))$ を使う$Q(y,x)/Q(x,y)$ を掛けるex-mcmc-counterexamples の 3 行目
連鎖の値を独立な標本とみて誤差を $\sigma/\sqrt n$ とする自己相関を考えて誤差を見積もる連鎖の値は互いに依存する
詳細釣り合いは定常であるための必要条件だとする十分条件にすぎないex-mcmc-gibbs の巡回走査
$\pi$ が定常なら必ず $\pi$ に収束するとする既約性が要る。分布の収束にはさらに周期をなくす条件(たとえば自己ループ)が要るex-mcmc-counterexamples

補足

状態が実数のベクトルで $\pi$ が密度で与えられる場合も、提案の密度 $Q(x,y)$ を使って受理確率を同じ式で定め、Metropolis–Hastings 法と Gibbs サンプリングが使われる(Mac03 §29.4–29.5、pp. 365–371)。状態が可算無限個の場合には、既約で周期がなく定常分布をもつ連鎖について $X_n$ の分布の収束が知られている(Dur19 Theorem 5.6.6、p. 313。この記事では証明しない)。収束の速さは推移行列の固有値や、分布が $\pi$ に十分近づくまでの回数(混合時間)で調べられる。
MCMC は、Ising 模型などの統計物理の計算、Bayes 統計の事後分布の計算、組合せ的な対象の数え上げの近似に使われる。期待値の見積もりそのものについては Monte Carlo法 を参照のこと。

関連項目

参考文献

[1]
David J. C. MacKay, Information Theory, Inference, and Learning Algorithms, Cambridge University Press, 2003, §29.4 The Metropolis–Hastings method(pp. 365–366、受理確率 (29.31))、§29.5 Gibbs sampling(pp. 370–371)、§29.6 Terminology for Markov chain Monte Carlo methods(pp. 372–374、詳細釣り合い (29.46))

Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する