MCMCの基礎(メトロポリス・ヘイスティングス法)

MCMC(マルコフ連鎖モンテカルロ法)とは,マルコフ連鎖を活用して目標となる確率分布に従うサンプルを得る手法のことです。

この説明だけではわかりにくいので,MCMCの意味を理解するためにMCMCの代表例であるメトロポリス・ヘイスティングス法(MH法)を解説します。

MCMCの包含関係

問題設定

問題設定
  • 確率関数が π(x)\pi(x) である確率分布に従うサンプルを生成したい
  • π(x)\pi(x) は計算できないが,π(x)\pi(x) に比例する関数 π~(x)\tilde{\pi}(x) の計算はできる
問題のイメージと重要性

xx のとりうる値が x=1,2,3x=1,2,3」で「π(x)\pi(x)x2x^2 に比例する」ような確率分布 π(x)\pi(x) に従うサンプルを生成したい状況を考えます。つまり π~(x)=x2(x=1,2,3)\tilde{\pi}(x)=x^2\:(x=1,2,3) です。

この場合は正規化定数が計算できて π(x)=112+22+32x2=114x2\pi(x)=\dfrac{1}{1^2+2^2+3^2}x^2=\dfrac{1}{14}x^2 なので直接 π(x)\pi(x) からサンプリングできますが,正規化定数が計算できない(π~(x)\tilde{\pi}(x) はわかるが π(x)\pi(x) はわからない)場合もあります。

例えば,ベイズ推定の文脈でこの問題設定は頻出です。→ベイズ推定の簡単な例と利点

P(XD)=P(DX)P(X)P(D)P(X\mid D)=\dfrac{P(D\mid X)P(X)}{P(D)} において,分子はわかるが分母が計算できない状況です。つまり π~=P(DX)P(X)\tilde{\pi}=P(D\mid X)P(X) が計算できることを利用して π=P(DX)P(X)P(D)\pi=\dfrac{P(D\mid X)P(X)}{P(D)} からのサンプルを得たい,という状況です。

メトロポリス・ヘイスティングス法

適当な点 x1x_1 から,ある規則に従って順に x2,x3,x4...x_2,x_3,x_4... を生成していきます。規則をうまく作って nn が十分大きいとき xnx_n が従う分布が π(x)\pi(x) になるようにします。

つまり,以下の手順です。

メトロポリス・ヘイスティングス法の概要
  1. 初期点 x1x_1 をランダムに(適当に)決める
  2. t=1,2,...t=1,2,... に対して,後述の「規則」に従って xtx_t から xt+1x_{t+1} を生成する
  3. 十分大きい nn に対して xnx_n を出力する

xtx_t から xt+1x_{t+1} を決める「規則」:
手順1. 適当に決めた分布(→補足) Q(xx)Q(x'\mid x) をもとに次の点 xt+1x_{t+1} の候補を決める。つまり確率分布 Q(xxt)Q(x'\mid x_t) からサンプルを生成し xt+1x_{t+1}' とする
手順2. 確率 min(1,π~(xt+1)Q(xtxt+1)π~(xt)Q(xt+1xt))\min\left(1,\dfrac{\tilde{\pi}(x_{t+1}')Q(x_t\mid x_{t+1}')}{\tilde{\pi}(x_t)Q(x_{t+1}'\mid x_t)}\right) で候補を採択する。つまり xt+1=xt+1x_{t+1}=x_{t+1}' として終了。そうでない場合は棄却,つまり xt+1x_{t+1}' を捨てて手順1に戻る

補足:提案分布

Q(xx)Q(x'\mid x) のことを提案分布と言います。提案分布は計算できる(サンプリングできる)分布なら何でも良いです。例えば(xx にすら依存しない) xx' についての一様分布とすればOKです。

※ただし,提案分布が目標分布 π(x)\pi(x') に近くなるように選べれば,上記の手順2で棄却される確率が低くなるので計算時間が短くなります。

うまくいく理由・証明

定理

QQ がゆるい条件を満たせば)メトロポリス・ヘイスティングス法で得られる点列 x1,x2,...,x_1,x_2,..., について,nn を大きくしていくと,xnx_n が従う分布は限りなく π(x)\pi(x) に近づく。

ただし,ゆるい条件とは「マルコフ連鎖が既約で非周期的になる」という条件です。例えば一様分布など任意の x,xx,x' に対して Q(xx)>0Q(x'\mid x)>0 ならこの条件を満たします。

この定理の証明はとてもおもしろいですが少し大変です。

定常分布・極限分布・詳細釣り合い条件・収束定理 で述べた「詳細釣り合い条件に関する定理」と「収束定理」を使います。

証明

メトロポリス・ヘイスティングス法において,xtx_t から xt+1x_{t+1} を生成する規則は tt によらない。また,xt+1x_{t+1}xtx_{t} のみから決まる。つまりマルコフ連鎖である。

「規則」において xx から yy へ遷移する確率を P(x,y)P(x,y) とおくと,任意の x,yx,y に対して π(x)P(x,y)=π(y)P(y,x)\pi(x)P(x,y)=\pi(y)P(y,x) を満たす(→後述の補足計算)。

つまり,詳細釣り合い条件を満たすので π(x)\pi(x) はこのマルコフ連鎖の定常分布である。

また,ゆるい条件よりマルコフ連鎖は既約で非周期的である。よって,収束定理により π(x)\pi(x) はこのマルコフ連鎖の極限分布でもある。つまり,x1x_1 をどの点にとっても nn を大きくしていくと,xnx_n が従う分布は限りなく π(x)\pi(x) に近づく。

※以下,補足計算:
手順2における採択確率を AA とおく。つまり, A(x,y)=min(1,π~(y)Q(xy)π~(x)Q(yx))A(x,y)=\min\left(1,\dfrac{\tilde{\pi}(y)Q(x\mid y)}{\tilde{\pi}(x)Q(y\mid x)}\right) とおく。遷移する確率の比は「QQ で候補に選ばれる」かつ「採択される」確率の比なので, P(x,y)P(y,x)=Q(yx)A(x,y)Q(xy)A(y,x)\dfrac{P(x,y)}{P(y,x)}=\dfrac{Q(y\mid x)A(x,y)}{Q(x\mid y)A(y,x)} となる。ところが,A(x,y)A(x,y)A(y,x)A(y,x)min\min の第二引数部分は互いに逆数なのでどちらかが 11 以上でもう片方が 11 以下であるので,結局どちらの場合も A(x,y)A(y,x)=π~(y)Q(xy)π~(x)Q(yx)\dfrac{A(x,y)}{A(y,x)}=\dfrac{\tilde{\pi}(y)Q(x\mid y)}{\tilde{\pi}(x)Q(y\mid x)} 以上より P(x,y)P(y,x)=π~(y)π~(x)=π(y)π(x)\dfrac{P(x,y)}{P(y,x)}=\dfrac{\tilde{\pi}(y)}{\tilde{\pi}(x)}=\dfrac{\pi(y)}{\pi(x)} 分母を払うと詳細釣り合い条件を得る。

補足

サンプルがたくさん欲しいとき

  • 十分大きい nn に対して xnx_n だけでなく xn+1,xn+2,...x_{n+1},x_{n+2},... なども選べばそれらはすべて π(x)\pi(x) に従うのでサンプルがたくさん得られます。
  • しかし,xnx_nxn+1x_{n+1} には相関があります。独立なサンプルをたくさん欲しい場合,十分大きい nnTT に対して xn,xn+T,xn+2T,...x_n,x_{n+T},x_{n+2T},... を選べばよいです。
  • nnTT を大きくすれば良いサンプルが得られますが,計算時間は長くなります。

メトロポリス法

Q(xx)=Q(xx)Q(x'\mid x)=Q(x\mid x') のとき,手順2の採択確率は QQ の部分が消えて min(1,π~(xt)π~(xt+1))\min\left(1,\dfrac{\tilde{\pi}(x_t)}{\tilde{\pi}(x_{t+1}')}\right) になります。この場合のメトロポリス・ヘイスティングス法のことをメトロポリス法と言います。

「MCMC」も「メトロポリス・ヘイスティングス」も名前がかっこいいです。