ギブスサンプリング
ギブスサンプリングにより2変量正規分布からのサンプリングを行った例。共分散を \sigma_{x_1 x_2} = 0.5 とした。
次のような状況を考える。
- 確率分布 P(\bm x) に従う乱数列 \{\bm x^t\}_{t=1}^T を生成したい
- 目標となる確率分布 P(\bm x) から直接乱数 \bm x をサンプルすることは難しい
- 条件付き確率 P(x_n | \bm x_{\setminus n}) からサンプルすることは容易い
ここで \bm x_{\setminus n} \in \R^{N-1} は、\bm x \in \R^N から n 番目の要素を除外したような配列である。
\begin{aligned}
\bm x \coloneqq \left[\begin{darray}{} x_1 \\ \vdots \\ x_{n-1} \\ x_n \\ x_{n+1} \\ \vdots \\ x_N \end{darray}\right]
\implies
\bm x_{\setminus n} \coloneqq \left[\begin{darray}{} x_1 \\ \vdots \\ x_{n-1} \\ x_{n+1} \\ \vdots \\ x_N \end{darray}\right]
\end{aligned}
メトロポリス・ヘイスティングス法からの導入
M-H法では、以下のようなアルゴリズムを利用すればよいのであった。
\begin{darray}{rl}
1. & X \gets \{\} \\
2. & \text{for } t \in \{1, 2, \dots, T\} : \\
3. & \qquad y \sim Q(y | x) \\
4. & \qquad r \sim {\cal U}(r | 0, 1) \\
5. & \qquad \text{if } A(y | x) \gt r : \\
6. & \qquad \qquad x \gets y \\
7. & \qquad X \gets X \cup \{x\} \\
8. & \operatorname{return} X
\end{darray}
ここで Q(\bm y | \bm x) は提案分布と呼ばれる確率分布、A(\bm y | \bm x) は受理確率と呼ばれる関数で、次のような式で定義される。
\begin{aligned}
A(\bm y | \bm x) = \frac{P(\bm y) Q(\bm x | \bm y)}{P(\bm x) Q(\bm y | \bm x)}
\end{aligned}
受理確率を 1 にする
問題は提案分布をどう定義するかであるが、ここでは次のようなものを使用する。
\begin{aligned}
Q(\bm y | \bm x) =
\left\lbrace\begin{aligned}
& P(y_n | \bm x_{\setminus n}) && \text{if } \bm y_{\setminus n} = \bm x_{\setminus n} \\
& 0 && \text{otherwise}
\end{aligned}\right.
\end{aligned}
すなわち、\bm y の n 番目以外の要素はすべて \bm x のものと同じにし (\bm y_{\setminus n} = \bm x_{\setminus n})、n 番目の要素は P(y_n | \bm x_{\setminus n}) に従うようにサンプリングする のである。
このようにすると、提案分布 Q(\bm y | \bm x) からサンプルされた \bm y は必ず \bm y_{\setminus n} = \bm x_{\setminus n} を満たす。したがって、次のように受理確率が必ず 1 になる。
\begin{aligned}
A(\bm y | \bm x)
={}& \frac{P(\bm y) Q(\bm x | \bm y)}{P(\bm x) Q(\bm y | \bm x)} \\
={}& \frac{P(\bm y) P(x_n | \bm y_{\setminus n})}{P(\bm x) P(y_n | \bm x_{\setminus n})} \\
&\left|\small\quad\begin{aligned}
P(\bm x)
={}& P(x_n, \bm x_{\setminus n}) \\
={}& P(x_n | \bm x_{\setminus n}) P(\bm x_{\setminus n}) \\
P(\bm y)
={}& P(y_n, \bm y_{\setminus n}) \\
={}& P(y_n | \bm y_{\setminus n}) P(\bm y_{\setminus n})
\end{aligned}\right. \\
={}& \frac{P(y_n | \bm y_{\setminus n}) P(\bm y_{\setminus n}) P(x_n | \bm y_{\setminus n})}{P(x_n | \bm x_{\setminus n}) P(\bm x_{\setminus n}) P(y_n | \bm x_{\setminus n})} \\
&\left|\small\quad\begin{aligned}
\bm y_{\setminus n} = \bm x_{\setminus n}
\end{aligned}\right. \\
={}& \frac{P(y_n | \bm x_{\setminus n}) P(\bm x_{\setminus n}) P(x_n | \bm x_{\setminus n})}{P(x_n | \bm x_{\setminus n}) P(\bm x_{\setminus n}) P(y_n | \bm x_{\setminus n})} \\
={}& 1
\end{aligned}
アルゴリズム
以上の考察から、次のようなアルゴリズムを得る。
\begin{darray}{rl}
1. & X \gets \{\} \\
2. & \text{for } t \in \{1, 2, \dots, T\} : \\
3. & \qquad \text{for } n \in \{1, 2, \dots, N\} : \\
4. & \qquad \qquad x_n \sim Q(x_n | \bm x_{\setminus n})\\
5. & \qquad \qquad X \gets X \cup \{\bm x\} \\
6. & \operatorname{return} X
\end{darray}
このように、元の確率分布から直接サンプルする代わりに、条件付き確率分布から次々と各要素をサンプルすることによって乱数列を近似的に生成する方法をギブスサンプリング (Gibbs sampler) という。
メリット
受理確率を1にすることのメリットは大きく2つある。
- 提案された乱数が棄却されることはないため、乱数列 \{\bm x^t\}_{t=1}^T を長くする必要がない
- メトロポリステストを実施する必要がないため、アルゴリズムは非常にシンプルなものとなる
もちろん、ギブスサンプリングが使えるのは条件付き確率分布 P(x_n | \bm x_{\setminus n}) からのサンプリングが容易な場合のみで、より一般のM-H法に比べれば限定的である。
とはいえ、この方法が使える場合には非常に効率のよい乱数サンプリングアルゴリズムとして重宝する。実際、さまざまな最適化アルゴリズムやシミュレーションアルゴリズムに応用されている。
Discussion