ギブスサンプリング｜シミュレーション

← シミュレーション一覧

🎓 レベル：発展　|　重要度：A（必須）

📎 前提：メトロポリス・ヘイスティングス法　|　ベイズ：ギブスサンプリング

要点（BLUF）

ギブスサンプリング：多変量分布を、各変数を順に**「他を固定したときの条件付き分布」**から引いて更新する MCMC。
提案が条件付き分布そのものなので受理率は常に100%——メトロポリス・ヘイスティングス（メトロポリス・ヘイスティングス法）の特殊形です。
条件付き分布が閉形式で書けるときに強力。相関0.8の二変量正規を交互更新で再現し、相関0.798を復元します。

1. アルゴリズム

$d$ 個の変数 $\mathbf{x} = (x_1, \dots, x_d)$ の同時分布 $\pi(\mathbf{x})$ からサンプルしたいとき、各変数を他のすべてを固定した条件付き分布から順に引きます。

$x_1^{(t+1)} \sim \pi(x_1 \mid x_2^{(t)}, x_3^{(t)}, \dots, x_d^{(t)})$
$x_2^{(t+1)} \sim \pi(x_2 \mid x_1^{(t+1)}, x_3^{(t)}, \dots, x_d^{(t)})$
… 最後の $x_d$ まで順に更新して1ステップ完了

最新の値を使いながら1変数ずつ更新するのがポイント（更新したらすぐ次の条件に反映）。これを繰り返すと、状態列は同時分布 $\pi(\mathbf{x})$ を定常分布とするマルコフ連鎖になります。

2. なぜ受理率100%か

ギブスは MH で提案分布を「条件付き分布 $q(x_i'|x_{-i}) = \pi(x_i'|x_{-i})$ 」に選んだ特殊ケースです。このとき受理確率を計算すると、 $\pi(\mathbf{x}) = \pi(x_i|x_{-i})\pi(x_{-i})$ を使って

\alpha = \min\!\left(1,\ \frac{\pi(x')\,q(x|x')}{\pi(x)\,q(x'|x)}\right) = \min\!\left(1,\ \frac{\pi(x_i'|x_{-i})\pi(x_{-i})\cdot \pi(x_i|x_{-i})}{\pi(x_i|x_{-i})\pi(x_{-i})\cdot \pi(x_i'|x_{-i})}\right) = \min(1, 1) = 1

きれいに約分されて常に1。提案が「その変数の正しい条件付き分布」なので、棄却する理由がない——だから100%受理で効率的です。条件付き分布から直接サンプルできること（閉形式や標準分布であること）が使える条件です。

3. 実装：相関0.8の二変量正規

二変量標準正規（相関 $\rho=0.8$ ）の条件付き分布は、 $x \mid y \sim \mathcal{N}(\rho y,\ 1-\rho^2)$ 、 $y \mid x \sim \mathcal{N}(\rho x,\ 1-\rho^2)$ と閉形式。これを交互に引きます。

import numpy as np

# 乱数シードを固定
rng = np.random.default_rng(42)

rho = 0.8
n = 200_000
x = np.empty(n); y = np.empty(n)
xc, yc = 0.0, 0.0
sd = np.sqrt(1 - rho**2)
for i in range(n):
    xc = rng.normal(rho * yc, sd)   # x | y からサンプル
    yc = rng.normal(rho * xc, sd)   # y | x からサンプル（更新後のxを使う）
    x[i] = xc; y[i] = yc

burn = 2000
print(f"平均  x={x[burn:].mean():.3f}  y={y[burn:].mean():.3f}  (目標 0, 0)")
print(f"標準偏差 x={x[burn:].std():.3f}  y={y[burn:].std():.3f}  (目標 1, 1)")
print(f"相関 corr(x,y) = {np.corrcoef(x[burn:], y[burn:])[0,1]:.3f}  (目標 0.8)")

出力：

平均  x=0.000  y=0.001  (目標 0, 0)
標準偏差 x=0.999  y=1.000  (目標 1, 1)
相関 corr(x,y) = 0.798  (目標 0.8)

出力の意味：条件付き正規から交互に引くだけで、二変量正規の平均（0,0）・標準偏差（1,1）・相関0.798（目標0.8）を見事に復元。同時分布を直接サンプルする式を持っていなくても、条件付き分布さえ書ければ同時分布が再現できる——これがギブスの威力です。ベイズの階層モデルでは各パラメータの条件付き分布が共役で書けることが多く、ギブスが定番になります（ベイズ）。

4. 強みと弱み

強み：受理率100%、チューニング不要（歩幅がない）、条件付きが共役なら高速。
弱み：条件付き分布が閉形式で書けないと使えない（書けなければ各ステップで MH を挟む＝Metropolis-within-Gibbs）。
弱み：変数間の相関が強いと収束が遅い。軸に沿ってしか動けないので、強相関の細長い分布をジグザグにしか進めません。相関0.8でも自己相関が残ります。相関を緩める再パラメータ化やブロック更新が対策です。

数式の直観的意味

ギブスは「一度に1軸ずつ、その軸の正しい断面分布に沿って動く」探索です。 $x$ を更新するときは $y$ を固定した「縦の断面」、次は $x$ を固定した「横の断面」——各断面では条件付き分布が正確に分かるので、提案がそのまま正解（受理率100%）。ただし軸に平行にしか動けないので、分布が斜めに細長い（強相関）と、対角線方向へ進むのに小刻みなジグザグを強いられ、収束が遅くなる。これがギブスの幾何的な限界で、MH のランダムウォーク（任意方向に提案できる）との使い分けの基準です。条件付きが書けるなら受理率100%のギブス、書けない・強相関ならブロック化や HMC を検討します。

⚠️ よくある誤解・落とし穴

「更新は古い値で一斉に」ではない：1変数更新したらすぐ次の条件に反映（最新値を使う）。一斉更新は別物（正しくない）。
「受理率100%だから自己相関もゼロ」ではない：受理率と自己相関は別。強相関分布では受理100%でも自己相関が大きい。
「条件付き分布が書けなくても使える」ではない：書けないと使えません。その場合は Metropolis-within-Gibbs。
「どの順で更新してもいい」ではない：系統スキャン（固定順）とランダムスキャンで挙動が変わります。固定順は厳密には可逆でない場合あり（定常分布は同じ）。
「ギブスは MH と無関係」ではない：MH の提案を条件付き分布にした特殊ケースです。

対応シミュレーション参照

本文の二変量正規ギブス（相関0.8・default_rng(42)、相関0.798復元）。応用はベイズへ。