ギブスサンプリングを今更理解する

趣旨

MCMCの一つであるギブスサンプリングを理解する。

MCMCって?

前の記事を参考にしてください。

usapyoi.hatenablog.com


では解説する。

ギブスサンプリングとは

MCMCの一つで、Metropolis-Hastingsアルゴリズムの特殊な例と言える。

問題設定

今、区間 (0,1)で一様に分布する擬似乱数を発生させるアルゴリズムが得られている。これを元に、 p(\boldsymbol{z})に従う L個のサンプル \boldsymbol{z}^{(1)}, \boldsymbol{z}^{(2)}, \cdots, \boldsymbol{z}^{(L)}を抽出したい。ここで、 \boldsymbol{z} D次元変数 \boldsymbol{z} = (z_1, z_2, \cdots, z_D)であるとする。

ギブスサンプリング

1. 手順

  1.  \boldsymbol{z} = (z_1, z_2, \cdots, z_D)に初期値 \boldsymbol{z}^{(0)} = (z_1^{(0)}, z_2^{(0)}, \cdots, z_D^{(0)})を与える。
  2.  \tau = 1, 2, \cdots, Tについて、以下のSTEP 3, 4, 5を行う。
  3.  d=1, 2, \cdots, Dについて、以下のSTEP 4, 5を行う。
  4.  z_d^{(\tau)} p(z_d | \boldsymbol{z}_{\setminus d}^{(\tau-1)})からサンプリングする.。ただし \boldsymbol{z}_{\setminus d} z_d以外の \boldsymbol{z}の要素。
  5. 得られた z_d^{(\tau)} \boldsymbol{z}_{\setminus i}^{(\tau-1)}を組み合わせたもの \boldsymbol{z}^*をサンプルとして保存する。

この手順によって得られる \boldsymbol{z} DTなので、 T = L/Dとすれば最短のイテレーションで済む。ただし、初期値依存性を排除するために、STEP 5を行わない、いわゆる「burn-in」を設けることが必要。

2. なぜ上手くいくのか

ギブスサンプリングはMetropolis-Hastingsアルゴリズムの特殊な例である、と先に述べた。そうであることを示せば、あとは以前の記事で紹介した通りうまくいくことが示される。

2-1. 復習 Metropolis-Hastingsアルゴリズムの手順

現在の状態が \boldsymbol{z}^{(\tau)}であるとき、分布 q_k(\boldsymbol{z} | \boldsymbol{z}^{(\tau)})からサンプル \boldsymbol{z}^*を抽出し、確率 A_k (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})に従って採択する。ただし、


\begin{align}
A_k (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) = \min  \left( 1, \frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \right)
\end{align}

また、状態 \boldsymbol{z}^{(\tau)}の更新則は、


\boldsymbol{z}^{(\tau + 1)} =
\left\{
    \begin{array}{cc}
        \boldsymbol{z}^* & \text{ if } \boldsymbol{z}^* \text{ is accepted} \\
        \boldsymbol{z}^{(\tau)} & \text{ otherwise}
    \end{array}
\right.

で与えられる。提案分布はこのマルコフ連鎖が既約で非周期的であれば良い。十分条件として、任意の \boldsymbol{z}_1, \boldsymbol{z}_2に対して q(\boldsymbol{z}_2 | \boldsymbol{z}_1) > 0であることが挙げられる。

2-2. ギブスサンプリングとの対応関係

ギブスサンプリングでは、目的の(正規化していない)分布 \tilde{p}(\boldsymbol{z})や提案分布 q(\boldsymbol{z})は、どちらも目的分布 p(\boldsymbol{z})になっている。また、ギブスサンプリングにおけるサンプル \boldsymbol{z}^*の生成では、一つの要素 z_kのみが新たに生成され、残る要素 \boldsymbol{z}_{\setminus k}^*は直前の状態、つまり \boldsymbol{z}^{(\tau)}を流用する。まとめると、


\begin{align}
\tilde{p}(\boldsymbol{z}^*) &= p(\boldsymbol{z}^*) = p(z_k^* | \boldsymbol{z}_{\setminus k}^{(\tau)}) p(\boldsymbol{z}_{\setminus k}^{(\tau)}) \\\
\tilde{p}(\boldsymbol{z}^{(\tau)}) &= p(\boldsymbol{z}^{(\tau)}) = p(z_k^{(\tau)} | \boldsymbol{z}_{\setminus k}^{(\tau)}) p(\boldsymbol{z}_{\setminus k}^{(\tau)}) \\\
q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)}) &= p (z_k^*| \boldsymbol{z}_{\setminus k}^{(\tau)}) \\\
q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*) &= p (z_k^{(\tau)}| \boldsymbol{z}_{\setminus k}^*) \\\
\end{align}

この条件下でサンプリングを行った時、サンプル \boldsymbol{z}^*の採択確率は


\frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} = 1

となる。よって、ギブスサンプリングは、Metropolis-Hastingsアルゴリズムの、サンプルが必ず採択されるような特殊な例であると言える。

まとめ

ギブスサンプリングについて勉強した。言うほど難しくはないっすね。PRMLではギブスサンプリングの欠点や変種などを紹介していたが、ここでは触れない。

おまけ

割と最近までギブソンサンプリングだと思っていたことを告白します。