ギブスサンプリング

ギブスサンプリングをMH法の特殊例という観点で導入する。

状況設定メトロポリス・ヘイスティングス法からの導入受理確率を 1 にするアルゴリズムメリット

状況設定

次のような状況を考える。

確率分布 $P(\bm x)$ に従う乱数列 $\{\bm x^i\}$ を生成したい
目標となる確率分布 $P(\bm x)$ から直接乱数 $\bm x$ をサンプルすることは難しい
条件付き確率 $P(x_j | \bm x_{\setminus j})$ からサンプルすることは容易い

ここで $\bm x_{\setminus j}$ は、 $\bm x$ から $j$ 番目の要素を除外したような配列である。

\begin{aligned} \bm x \coloneqq \left[\begin{darray}{} x_1 \\ \vdots \\ x_{j-1} \\ x_j \\ x_{j+1} \\ \vdots \\ x_N \end{darray}\right] \implies \bm x_{\setminus j} \coloneqq \left[\begin{darray}{} x_1 \\ \vdots \\ x_{j-1} \\ x_{j+1} \\ \vdots \\ x_N \end{darray}\right] \end{aligned}

メトロポリス・ヘイスティングス法からの導入

M-H法では、以下のようなアルゴリズムを利用すればよいのであった。

Metropolis-Hastings Algorithm

\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 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_j | \bm x_{\setminus j}) && \text{if } \bm y_{\setminus j} = \bm x_{\setminus j} \\ & 0 && \text{otherwise} \end{aligned}\right. \end{aligned}

すなわち、 $j$ 番目以外の要素はすべて $\bm x$ のものと同じにし ( $\bm y_{\setminus j} = \bm x_{\setminus j})$ 、 $j$ 番目の要素は $P(y_j | \bm x_{\setminus j})$ に従うようにサンプリングする。このようにすると、提案分布 $Q(\bm y | \bm x)$ からサンプルされた $\bm y$ は必ず $\bm y_{\setminus j} = \bm x_{\setminus j}$ を満たす。これに注意すると、受理確率が必ず $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_j | \bm y_{\setminus j})}{P(\bm x) P(y_j | \bm x_{\setminus j})} \\ &\left|\small\quad\begin{aligned} P(\bm x) ={}& P(x_j, \bm x_{\setminus j}) \\ ={}& P(x_j | \bm x_{\setminus j}) P(\bm x_{\setminus j}) \end{aligned}\right. \\ ={}& \frac{P(y_j | \bm y_{\setminus j}) P(\bm y_{\setminus j}) P(x_j | \bm y_{\setminus j})}{P(x_j | \bm x_{\setminus j}) P(\bm x_{\setminus j}) P(y_j | \bm x_{\setminus j})} \\ &\left|\small\quad\begin{aligned} \bm y_{\setminus j} = \bm x_{\setminus j} \end{aligned}\right. \\ ={}& \frac{P(y_j | \bm x_{\setminus j}) P(\bm x_{\setminus j}) P(x_j | \bm x_{\setminus j})}{P(x_j | \bm x_{\setminus j}) P(\bm x_{\setminus j}) P(y_j | \bm x_{\setminus j})} \\ ={}& 1 \end{aligned}

アルゴリズム

以上の考察から、次のようなアルゴリズムを得る。

Algorithm: Gibbs sampler.
$\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 3. & \qquad \text{for } j \in \{1, 2, \dots, N\} : \\ 4. & \qquad \qquad x_j \sim Q(x_j | \bm x_{\setminus j})\\ 5. & \qquad \qquad X \gets X \cup \{\bm x\} \\ 6. & \operatorname{return} X \end{darray}$

このように、元の確率分布から直接サンプルする代わりに、条件付き確率分布から次々と各要素をサンプルすることによって乱数列を近似的に生成する方法をギブスサンプリング (Gibbs sampling; Gibbs sampler) という。

メリット

受理確率を1にすることのメリットは大きく2つある。

提案された乱数 $\bm y$ が棄却されることはないため、乱数列 $\{\bm x^t\}_{t=1}^n$ を長くする必要がない
メトロポリステストを実施する必要がないため、アルゴリズムは非常にシンプルなものとなる

もちろん、ギブスサンプリングが使えるのは条件付き確率分布 $P(x_j | \bm x_{\setminus j})$ からのサンプリングが容易な場合のみで、より一般のM-H法に比べれば限定的である。

とはいえ、この方法が使える場合には非常に効率のよい乱数サンプリングアルゴリズムとして重宝するのはいうまでもない。実際、さまざまな最適化アルゴリズムやシミュレーションアルゴリズムに応用されている。