ギブスサンプリング

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：メトロポリスヘイスティングス　|　数理：MCMC（マルコフ連鎖モンテカルロ）（統計）・正規正規モデル

要点（BLUF）

ギブスサンプリングは、多次元の事後を「1成分ずつ、他を固定した完全条件付き分布から順に引く」MCMC。受容判定が要らず、提案は常に受容されます。
前提は「各成分の完全条件付き分布が既知でサンプリング可能」なこと。共役構造（第2章）があると条件付きが標準分布（正規・逆ガンマなど）になり、ギブスがそのまま回ります。
平均も分散も未知の正規モデルで実装： $\mu\mid\sigma^2$ は正規（正規正規モデルの精度の加重平均）、 $\sigma^2\mid\mu$ は逆ガンマ。2次元グリッドの真値と一致することを確かめます。

1. アイデア：1変数ずつの条件付きに分解する

$d$ 次元の同時事後 $p(\theta_1,\dots,\theta_d\mid D)$ を直接サンプリングするのは難しい。ギブスは「他の成分を今の値に固定したときの $\theta_j$ の条件付き分布」＝**完全条件付き分布（full conditional）**から、 $j=1,\dots,d$ と順に引きます。

\theta_j^{(t+1)}\sim p\big(\theta_j\mid \theta_{-j}^{(\text{最新})},\,D\big)\qquad(\theta_{-j}=\theta\ \text{から}\ \theta_j\ \text{を除いたもの})

これを1サイクルとして繰り返すと、点列が同時事後 $p(\theta\mid D)$ からの標本になります（なぜ正しいか＝完全条件付きから引く操作が受容確率1の MH に当たる、は MCMC（マルコフ連鎖モンテカルロ））。

2. 正規モデルの完全条件付き（共役が効く）

データ $x_1,\dots,x_n\sim\mathcal N(\mu,\sigma^2)$ で、 $\mu$ と $\sigma^2$ の両方を推定します。事前は $\mu\sim\mathcal N(\mu_0,1/\tau_0)$ 、 $\sigma^2\sim\mathrm{Inv\text{-}Gamma}(a_0,b_0)$ （独立）。同時事後は閉形式で出ませんが、各成分の完全条件付きは第2章の共役そのものです。

$\mu\mid\sigma^2,D$ ：分散 $\sigma^2$ を固定すれば「分散既知の正規–正規」。精度の加重平均で

\mu\mid\sigma^2,D\ \sim\ \mathcal N\!\left(\frac{\tau_0\mu_0+(n/\sigma^2)\bar x}{\tau_0+n/\sigma^2},\ \frac{1}{\tau_0+n/\sigma^2}\right)

$\sigma^2\mid\mu,D$ ：平均 $\mu$ を固定すれば「平均既知の正規–逆ガンマ」。平方和 $\sum(x_i-\mu)^2$ を使って

\sigma^2\mid\mu,D\ \sim\ \mathrm{Inv\text{-}Gamma}\!\left(a_0+\frac{n}{2},\ b_0+\frac12\sum_i (x_i-\mu)^2\right)

どちらも標準分布なので直接サンプリングでき、交互に引くだけ。これがギブスの気持ちよさです。

3. numpy で実装し、2次元グリッドの真値と一致させる

import numpy as np
from scipy import stats
from scipy.integrate import trapezoid          # numpy 2.0+ で np.trapz は廃止

rng = np.random.default_rng(3)
data = rng.normal(5.0, 2.0, size=30)             # 真は μ=5, σ²=4
n = len(data); xbar = data.mean()
mu0, s0 = 0.0, 10.0; tau0 = 1/s0**2              # 事前 μ~N(0,10²)
ag, bg = 2.0, 2.0                                # 事前 σ²~InvGamma(2,2)

def gibbs(N=60000, burn=5000, seed=0):
    rng = np.random.default_rng(seed)
    mu, s2 = xbar, data.var()
    M = np.empty(N); S = np.empty(N)
    for i in range(N):
        # μ | σ², data ~ Normal（精度の加重平均）
        tau_n = tau0 + n/s2
        mu_n  = (tau0*mu0 + (n/s2)*xbar)/tau_n
        mu = rng.normal(mu_n, np.sqrt(1/tau_n))
        # σ² | μ, data ~ InvGamma(a0+n/2, b0+½Σ(x-μ)²)
        a_n = ag + n/2
        b_n = bg + 0.5*np.sum((data - mu)**2)
        s2 = stats.invgamma(a=a_n, scale=b_n).rvs(random_state=rng)
        M[i], S[i] = mu, s2
    return M[burn:], S[burn:]

mu_s, s2_s = gibbs()
print(f"ギブス:     E[μ]={mu_s.mean():.4f}  E[σ²]={s2_s.mean():.4f}  "
      f"μ95%=[{np.quantile(mu_s,0.025):.3f}, {np.quantile(mu_s,0.975):.3f}]")

# 真値：2次元グリッドで同時事後を数値積分して周辺平均を出す
mug = np.linspace(2, 8, 400); s2g = np.linspace(0.5, 12, 400)
MU, S2 = np.meshgrid(mug, s2g, indexing='ij')
logpost = (-0.5*n*np.log(S2)
           - 0.5*((data[:,None,None]-MU[None,:,:])**2).sum(0)/S2
           - 0.5*tau0*(MU-mu0)**2
           + (-(ag+1))*np.log(S2) - bg/S2)
post = np.exp(logpost - logpost.max())
post /= trapezoid(trapezoid(post, s2g, axis=1), mug)
Emu = trapezoid(trapezoid(post*MU, s2g, axis=1), mug)
Es2 = trapezoid(trapezoid(post*S2, s2g, axis=1), mug)
print(f"グリッド真値: E[μ]={Emu:.4f}  E[σ²]={Es2:.4f}")

出力：

ギブス:     E[μ]=5.1073  E[σ²]=5.1068  μ95%=[4.286, 5.923]
グリッド真値: E[μ]=5.1076  E[σ²]=5.1067

出力の意味：ギブスの事後平均 $E[\mu]=5.107$ 、 $E[\sigma^2]=5.107$ が、2次元グリッドで数値積分した真値とほぼ完全に一致します。同時事後は閉形式で書けないのに、完全条件付き（正規・逆ガンマ）を交互に引くだけで正しい事後を再現できました。共役（第2章）が「条件付きを標準分布にする」形でここに効いています。

4. ギブスの利点と前提

受容判定が不要：完全条件付きから直接引くので、提案は常に受容。MH の特殊ケース（受容確率がつねに1）と見なせます。棄却で足踏みしない分、効率がよいことが多い。
前提が厳しい：各成分の完全条件付きが既知でサンプリング可能でなければなりません。共役な階層モデルでは成立しますが、一般には崩れます。
メトロポリス内ギブス：条件付きから直接引けない成分があれば、その成分だけ MH を使います（メトロポリスヘイスティングス）。
相関に注意：成分どうしが強く相関すると、1軸ずつ動くギブスは進みが遅くなります（収束診断で診断。勾配を使う HMC が有利な場面がハミルトニアンモンテカルロとNUTS）。

サンプルが2次元事後をどう埋めるかを散布図で見ます。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用

rng = np.random.default_rng(3)
data = rng.normal(5.0, 2.0, size=30); n = len(data); xbar = data.mean()
mu0, tau0, ag, bg = 0.0, 1/100, 2.0, 2.0
mu, s2 = xbar, data.var()
M = np.empty(20000); S = np.empty(20000)
for i in range(20000):
    tau_n = tau0 + n/s2; mu = rng.normal((tau0*mu0+(n/s2)*xbar)/tau_n, np.sqrt(1/tau_n))
    s2 = stats.invgamma(a=ag+n/2, scale=bg+0.5*np.sum((data-mu)**2)).rvs(random_state=rng)
    M[i], S[i] = mu, s2

plt.figure(figsize=(6, 5))
plt.scatter(M[::5], S[::5], s=3, alpha=0.15)
plt.xlabel("μ"); plt.ylabel("σ²")
plt.title("ギブス標本が描く同時事後 p(μ, σ² | D)")
plt.tight_layout(); plt.show()

グラフの意味：点群は $\mu\approx5.1$ 、 $\sigma^2\approx5$ を中心に広がり、同時事後の形を埋めます。 $\mu$ の散らばりは $\sigma^2$ が大きい行ほど横に広い（分散が大きいと平均の不確実性も増す）——2変数の依存が散布図に表れます。

⚠️ よくある誤解

「ギブスは MH と別物」ではない：完全条件付きから引く＝受容確率1の MH。受容判定が要らないだけです。
「ギブスはいつでも使える」ではない：完全条件付きが既知でサンプリング可能という前提が要る。崩れたら成分ごとに MH を混ぜます（メトロポリス内ギブス）。
「条件付きが正規・逆ガンマになるのは偶然」ではない：共役（第2章）の帰結。 $\mu\mid\sigma^2$ は分散既知の正規–正規、 $\sigma^2\mid\mu$ は平均既知の正規–逆ガンマです。
「ギブスなら相関を気にしなくてよい」ではない：成分が強相関だと1軸ずつのギブスは混合が遅い。診断（収束診断）と再パラメータ化（階層モデルの実例と再パラメータ化）で対処します。