ギブスサンプリング（Gibbs Sampling）は、MCMC（マルコフ連鎖モンテカルロ）法の代表的なアルゴリズムの1つです。ベイズ推定で事後分布からのサンプリングが解析的にできない場合に、近似的にサンプルを得る手法として広く使われています。
今回はギブスサンプリングの仕組みについて、できるだけ分かりやすく解説していきます。
本記事の内容
- MCMC法の概要とギブスサンプリングの位置づけ
- ギブスサンプリングのアルゴリズム
- 2次元ガウス分布を例にした導出
- Python での実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
MCMC法とは
MCMC（Markov Chain Monte Carlo）法は、複雑な確率分布からサンプルを生成するための手法群です。
ベイズ推定では、事後分布 $p(\theta|\mathcal{D})$ から期待値などの統計量を計算したい場面が多くあります。しかし、事後分布が複雑な場合、解析的に積分を計算することが困難です。
$$ E[\theta|\mathcal{D}] = \int \theta \cdot p(\theta|\mathcal{D}) \, d\theta $$
MCMC法は、目的の分布をマルコフ連鎖の定常分布として持つようなサンプル列を生成し、この期待値をモンテカルロ近似で計算します。
ギブスサンプリングの基本的な考え方
ギブスサンプリングの核心的なアイデアは、多変量の同時分布からの直接サンプリングが難しい場合でも、各変数の条件付き分布からのサンプリングは容易であることを利用する点にあります。
$D$ 次元の確率変数 $\bm{\theta} = (\theta_1, \theta_2, \dots, \theta_D)$ の同時分布 $p(\bm{\theta})$ からサンプルを得たい場合、ギブスサンプリングは以下のようにします。
ギブスサンプリングのアルゴリズム
- 初期値 $\bm{\theta}^{(0)} = (\theta_1^{(0)}, \theta_2^{(0)}, \dots, \theta_D^{(0)})$ を設定する
- $t = 1, 2, \dots$ について以下を繰り返す： – $\theta_1^{(t)} \sim p(\theta_1 | \theta_2^{(t-1)}, \theta_3^{(t-1)}, \dots, \theta_D^{(t-1)})$ – $\theta_2^{(t)} \sim p(\theta_2 | \theta_1^{(t)}, \theta_3^{(t-1)}, \dots, \theta_D^{(t-1)})$ – $\vdots$ – $\theta_D^{(t)} \sim p(\theta_D | \theta_1^{(t)}, \theta_2^{(t)}, \dots, \theta_{D-1}^{(t)})$
各ステップでは、1つの変数だけを更新し、それ以外の変数は直前の値を使います。更新する際は、他の変数を条件として与えた完全条件付き分布からサンプリングします。
2次元ガウス分布の例
最もわかりやすい例として、2次元ガウス分布を考えましょう。
$$ \begin{equation} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} \sim \mathcal{N}\left(\begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix}, \begin{pmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}\right) \end{equation} $$
ここで $\rho$ は相関係数です。
完全条件付き分布の導出
多変量ガウス分布の条件付き分布の公式より、
$$ \begin{align} p(x_1|x_2) &= \mathcal{N}\left(\mu_1 + \rho\frac{\sigma_1}{\sigma_2}(x_2 – \mu_2), \, (1-\rho^2)\sigma_1^2\right) \\ p(x_2|x_1) &= \mathcal{N}\left(\mu_2 + \rho\frac{\sigma_2}{\sigma_1}(x_1 – \mu_1), \, (1-\rho^2)\sigma_2^2\right) \end{align} $$
ギブスサンプリングでは、(2) と (3) を交互に使ってサンプルを生成します。
Python での実装
2次元ガウス分布に対するギブスサンプリングを実装しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
def gibbs_sampling_2d_gaussian(mu1, mu2, sigma1, sigma2, rho, n_samples, burn_in=500):
"""2次元ガウス分布に対するギブスサンプリング"""
samples = np.zeros((n_samples + burn_in, 2))
# 初期値
x1, x2 = 0.0, 0.0
for t in range(n_samples + burn_in):
# x1をサンプリング（x2で条件付け）
cond_mean1 = mu1 + rho * (sigma1 / sigma2) * (x2 - mu2)
cond_var1 = (1 - rho**2) * sigma1**2
x1 = np.random.normal(cond_mean1, np.sqrt(cond_var1))
# x2をサンプリング（x1で条件付け）
cond_mean2 = mu2 + rho * (sigma2 / sigma1) * (x1 - mu1)
cond_var2 = (1 - rho**2) * sigma2**2
x2 = np.random.normal(cond_mean2, np.sqrt(cond_var2))
samples[t] = [x1, x2]
# バーンイン期間を除去
return samples[burn_in:]
# パラメータ設定
mu1, mu2 = 0, 0
sigma1, sigma2 = 1.0, 1.0
rho = 0.8
n_samples = 5000
np.random.seed(42)
samples = gibbs_sampling_2d_gaussian(mu1, mu2, sigma1, sigma2, rho, n_samples)
# 理論分布
mean = [mu1, mu2]
cov = [[sigma1**2, rho*sigma1*sigma2],
[rho*sigma1*sigma2, sigma2**2]]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# 左: サンプリングの軌跡（最初の100ステップ）
ax1 = axes[0]
trace = samples[:100]
ax1.plot(trace[:, 0], trace[:, 1], 'b-', alpha=0.3, linewidth=0.5)
ax1.scatter(trace[:, 0], trace[:, 1], c=range(100), cmap='viridis', s=10)
ax1.set_xlabel('x1', fontsize=12)
ax1.set_ylabel('x2', fontsize=12)
ax1.set_title('Gibbs Sampling Trace (first 100)', fontsize=12)
ax1.grid(True, alpha=0.3)
# 中央: サンプルの散布図
ax2 = axes[1]
ax2.scatter(samples[:, 0], samples[:, 1], alpha=0.1, s=5, color='blue')
# 理論的な等高線
x_grid = np.linspace(-3.5, 3.5, 100)
y_grid = np.linspace(-3.5, 3.5, 100)
X, Y = np.meshgrid(x_grid, y_grid)
pos = np.dstack((X, Y))
rv = multivariate_normal(mean, cov)
ax2.contour(X, Y, rv.pdf(pos), levels=5, colors='red', linewidths=1.5)
ax2.set_xlabel('x1', fontsize=12)
ax2.set_ylabel('x2', fontsize=12)
ax2.set_title('Gibbs Samples vs True Contours', fontsize=12)
ax2.grid(True, alpha=0.3)
# 右: 周辺分布の比較
ax3 = axes[2]
ax3.hist(samples[:, 0], bins=50, density=True, alpha=0.5, color='blue', label='x1 (Gibbs)')
ax3.hist(samples[:, 1], bins=50, density=True, alpha=0.5, color='green', label='x2 (Gibbs)')
x_range = np.linspace(-4, 4, 200)
from scipy.stats import norm
ax3.plot(x_range, norm.pdf(x_range, mu1, sigma1), 'b-', linewidth=2, label='x1 (Theory)')
ax3.plot(x_range, norm.pdf(x_range, mu2, sigma2), 'g-', linewidth=2, label='x2 (Theory)')
ax3.set_xlabel('Value', fontsize=12)
ax3.set_ylabel('Density', fontsize=12)
ax3.set_title('Marginal Distributions', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
plt.suptitle(f'Gibbs Sampling from 2D Gaussian (rho={rho})', fontsize=14)
plt.tight_layout()
plt.show()
左のグラフでは、ギブスサンプリングが水平・垂直に交互にサンプルを生成する特徴的な軌跡が見えます。中央のグラフでは、サンプルが理論的な等高線と一致していることが確認できます。
ギブスサンプリングの特徴
メリット
- 完全条件付き分布が既知の場合、実装が容易
- メトロポリス・ヘイスティングス法のような棄却ステップが不要（常にサンプルが受理される）
- 共役事前分布を使ったベイズモデルと相性が良い
デメリット
- 相関が強い場合、サンプル間の自己相関が高くなり、収束が遅くなる
- 完全条件付き分布が解析的に求められない場合は直接適用できない
まとめ
本記事では、ギブスサンプリングについて解説しました。
- ギブスサンプリングはMCMC法の一種で、完全条件付き分布から交互にサンプリングする手法
- 多変量分布の直接サンプリングが困難でも、条件付き分布が簡単なら適用可能
- 2次元ガウス分布の例では、条件付き分布が1次元ガウス分布になるため簡単に実装できる
- バーンイン期間を設けることで、初期値の影響を除去する
次のステップとして、以下の記事も参考にしてください。