MCMCを今更理解する

趣旨

MCMC(マルコフ連鎖モンテカルロ)の理解を目的に、基本的なMetropolis-Hastingsアルゴリズムをまとめる。

なぜMCMC?

名前が可愛いから 情報学徒として当然だよね。

では解説する。

MCMCとは

サンプリングを行うアルゴリズムの群。MCMCという特定のアルゴリズムがあるわけではなく、あくまで大雑把な指針。

背景

期待値を評価したい。

\mathbb{E} \lbrack f \rbrack = \int f(\boldsymbol{z}) p(\boldsymbol{z}) d \boldsymbol{z}

しかし直接計算できないとする。

サンプリング法

サンプリング法とは、標本平均で期待値を近似しようというアプローチ。分布 p(\boldsymbol{z})に従う L個のサンプル \boldsymbol{z}^{(1)}, \boldsymbol{z}^{(2)}, \cdots, \boldsymbol{z}^{(L)}を抽出して、

\mathbb{E} \lbrack f \rbrack \approx \frac{1}{L} \sum_{l=1}^L f(\boldsymbol{z}^{(l)})

と近似する。このアプローチで最も重要なのは、どうやって分布 p(\boldsymbol{z})から独立にサンプリングするのか、という点。

問題設定

今、区間 (0,1)で一様に分布する擬似乱数を発生させるアルゴリズムが得られている。これを元に、 p(\boldsymbol{z})に従う L個のサンプル \boldsymbol{z}^{(1)}, \boldsymbol{z}^{(2)}, \cdots, \boldsymbol{z}^{(L)}を抽出したい。更に、任意の \boldsymbol{z}に対して p(\boldsymbol{z})は得られないが、ある未知の正規化定数 Z_pが存在して \tilde{p}(\boldsymbol{z}) = Z_p p(\boldsymbol{z})ならば容易に計算できるものとする。

  • 少し恣意的な問題設定に見えるかもしれないが、要するに \tilde{p}(\boldsymbol{z})は分かるけど Z_p = \int \tilde{p}(\boldsymbol{z}) d\boldsymbol{z}は分からないのだと考えれば良い。積分、難しいもんな。

棄却サンプリング

MCMCの前に、棄却サンプリングについて簡単に説明する。

1. 手順

  1. p(\boldsymbol{z})からサンプリングできないので、サンプリングが容易な分布 q(\boldsymbol{z})を用意する(提案分布という)
  2. 定数 kを任意の \boldsymbol{z}について kq(\boldsymbol{z}) \geq \tilde p(\boldsymbol{z})となるように定める
  3. \boldsymbol{z}^*を q(\boldsymbol{z})からサンプリングし、確率 \tilde{p}(\boldsymbol{z}^*) / \lbrace kq(\boldsymbol{z}^*) \rbraceで採択する( kの定義から、この値は必ず 1以下になる)

2. なぜ上手くいくのか

所望の分布 p(\boldsymbol{z})と区別するため、確率を示したい場合は \text{Pr}(\boldsymbol{z})と表記する。変数 aを、ステップ3において採択されれば a=1、棄却されれば a=0となるように定義する。この時、採択された \boldsymbol{z}の分布は \text{Pr}( \boldsymbol{z} | a=1)と書ける。ベイズの定理より、

\begin{align} \text{Pr}(\boldsymbol{z} | a=1) &= \frac{\text{Pr}(\boldsymbol{z}) \text{Pr}(a=1 | \boldsymbol{z})}{\text{Pr}(a=1)} \end{align}

と書ける。右辺を見ていこう。

\begin{align} \text{Pr}(\boldsymbol{z}) &= q(\boldsymbol{z}) \\\ \text{Pr}(a=1 | \boldsymbol{z}) &= \frac{\tilde{p}(\boldsymbol{z})}{kq(\boldsymbol{z})} \\\ &= \frac{Z_p p(\boldsymbol{z})}{kq(\boldsymbol{z})} \\\ \text{Pr}(a=1) &= \int \text{Pr}(a=1, \boldsymbol{z}) d \boldsymbol{z} \\\ &= \int \text{Pr}(\boldsymbol{z}) \text{Pr}(a=1 | \boldsymbol{z}) d \boldsymbol{z} \\\ &= \int q(\boldsymbol{z}) \frac{\tilde{p}(\boldsymbol{z})}{kq(\boldsymbol{z})} d \boldsymbol{z} \\\ &= \frac{Z_p}{k} \int p(\boldsymbol{z}) d \boldsymbol{z} \\\ &= \frac{Z_p}{k} \end{align}

これらを代入すれば \text{Pr}(\boldsymbol{z} | a=1) = p(\boldsymbol{z})が得られる。

3. 欠点

せっかくサンプリングしたのに高い確率で棄却してしまうのは非効率的なので、 kはできるだけ小さくあって欲しい。そのためには、 q(\boldsymbol{z})は p(\boldsymbol{z})に限りなく近いものであって欲しい。そんなこと最初から出来たら苦労しない。

MCMC

MCMCでは、提案分布 q(\boldsymbol{z})を現在の状態に依存させることでサンプリングをより効率的に行うことを目指す。この時、提案分布は現在の状態 \boldsymbol{z}^{(\tau)}に基づいた分布 q(\boldsymbol{z} | \boldsymbol{z}^{(\tau)})として表現される。

MCMCの例: Metropolis-Hastingsアルゴリズム

1. 手順

現在の状態が \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. なぜ上手くいくのか

あんまりよく分かりません 既約で非周期的なマルコフ連鎖が定常分布を持つなら、必ずその定常分布に収束する(収束定理)。問題はその定常分布が所望の分布 p(\boldsymbol{z})かどうかということ。

分布 p(\boldsymbol{z})が定常分布であることを保証するための十分条件は、以下のような詳細釣り合い条件が挙げられる。

p(\boldsymbol{z}^{(\tau)}) \text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*) = p(\boldsymbol{z}^*) \text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})

ここで、 \text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)は状態が \boldsymbol{z}^{(\tau)}から \boldsymbol{z}^*に遷移する確率を指し、今回の設定では

\begin{align} \text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*) &= q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)}) A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) \end{align}

と表される。 \text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})も同様。この時、詳細釣り合い条件を少し変形させたもの

\begin{align} \frac{\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)}{\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})} &= \frac{p(\boldsymbol{z}^*)}{p(\boldsymbol{z}^{(\tau)})} \end{align}

が本当に成立するのか確かめるために、左辺について考える。上で定義した式を代入すれば、

\begin{align} \frac{\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)}{\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})} &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)}) A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*) A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} \\\ &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \frac{A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} \end{align}

と書ける。ここで、最右辺の右側の分数について更に考察する。分母分子について、

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

なのだが、よく見ると \minの第二項が逆数になっているだけである。つまり、どちらかが 1より大きく、どちらかが 1より小さくなることが保証される。どちらがどうなるにせよ、先ほどの分数について、

\begin{align} \frac{A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} = \frac{\tilde{p}(\boldsymbol{z}^*) q (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \end{align}

となることは簡単な計算で分かるだろう。これを代入すれば、

\begin{align} \frac{\text{Pr} (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)}{\text{Pr} (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})} &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \frac{A (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)})}{A (\boldsymbol{z}^{(\tau)}, \boldsymbol{z}^*)} \\\ &= \frac{q(\boldsymbol{z}^* | \boldsymbol{z}^{(\tau)})}{q(\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)} \frac{\tilde{p}(\boldsymbol{z}^*) q (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \\\ &= \frac{\tilde{p}(\boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)})} = \frac{p(\boldsymbol{z}^*)}{p(\boldsymbol{z}^{(\tau)})} \end{align}

はい示せた。

結局

MCMCではまず提案分布 q(\boldsymbol{z})を置き、何回かのサンプリングでそれを p(\boldsymbol{z})に近似する。十分に近似できたらサンプルとして収集すれば良い。蒸留酒の最初の部分捨てるみたいだね。

まとめ

MCMCのコンセプトを理解するために、まずは基本的なMetropolis-Hastingsアルゴリズムをまとめた。変種としてギブスサンプリングやスライスサンプリングなどがあるが、勘弁してください今後理解する予定です。

おまけ

えむしーえむしー \leftarrowかわいい