階層モデルの構造｜ベイズ統計

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

📎 前提：収束診断　|　関連：正規正規モデル

要点（BLUF）

階層モデルは、グループ構造を持つデータを「各グループのパラメータ $\theta_j$ が**共通の親分布（ハイパープライアー）**から生まれた」とみなして結びます。
3つの極：完全プーリング（全グループ同一）／非プーリング（各グループ独立）／部分プーリング（階層・両者の中間に収縮）。階層モデルは部分プーリングを自動で実現します。
部分プーリングは精度の加重平均で、不確実な（データが少ない・誤差が大きい）グループほど全体平均へ強く収縮します。グループ間で情報を「借り合う」ので、小グループの推定が安定します。

1. グループ構造のデータと3つの選択肢

学校ごとのテスト効果、店舗ごとの転換率、選手ごとの打率——グループに分かれたデータはいたるところにあります。各グループ $j$ の観測平均 $y_j$ と、その不確実性（標準誤差） $\sigma_j$ があるとき、グループ平均 $\theta_j$ をどう推定するか。素朴な選択肢は両極端です。

完全プーリング：グループ差を無視し、全データをまとめて1つの $\theta$ を推定。→ グループの個性を消す。
非プーリング：各グループを完全に独立に $\theta_j=y_j$ 。→ データの少ないグループの推定が暴れる。

flowchart TD
  H["親分布 θ_j ~ N(μ, τ²)<br/>（全グループ共通のハイパープライアー）"] --> T1["θ_1"]
  H --> T2["θ_2"]
  H --> Tj["… θ_J"]
  T1 --> Y1["データ y_1 ~ N(θ_1, σ_1²)"]
  T2 --> Y2["データ y_2 ~ N(θ_2, σ_2²)"]
  Tj --> Yj["データ y_J ~ N(θ_J, σ_J²)"]

2. 部分プーリング＝階層モデル：親分布で結ぶ

階層モデルは中間を取ります。各 $\theta_j$ が共通の親分布から生まれたとモデル化し、親分布のパラメータ（全体平均 $\mu$ 、グループ間ばらつき $\tau$ ）もデータから推定します。

y_j\sim\mathcal N(\theta_j,\sigma_j^2),\qquad \theta_j\sim\mathcal N(\mu,\tau^2)

$\tau\to0$ なら全 $\theta_j$ が $\mu$ に一致＝完全プーリング、 $\tau\to\infty$ なら親が無情報で $\theta_j=y_j$ ＝非プーリング。有限の $\tau$ がその中間＝部分プーリングです。 $\mu,\tau$ を固定すれば、各 $\theta_j$ の事後平均は分散既知の正規–正規（正規正規モデル）そのもので、精度の加重平均になります。

\hat\theta_j=(1-B_j)\,y_j+B_j\,\mu,\qquad B_j=\frac{\sigma_j^2}{\sigma_j^2+\tau^2}\ (\text{収縮率})

3つを同じデータで比べます（ $\mu,\tau$ は周辺尤度最大化＝経験ベイズで推定。詳細は経験ベイズ）。

import numpy as np
from scipy import optimize

# J=8 グループ：観測平均 y_j と既知の標準誤差 σ_j（精密な群～粗い群）
rng = np.random.default_rng(2)
J = 8
sigma = np.array([2, 2, 3, 3, 4, 5, 8, 12.0])
theta_true = rng.normal(50.0, 10.0, J)            # 真のグループ平均（μ=50, τ=10）
y = rng.normal(theta_true, sigma)

# 完全プーリング：全グループ同一 → 精度加重平均
mu_complete = np.sum(y/sigma**2) / np.sum(1/sigma**2)
# 非プーリング：各グループ独立 → そのまま y_j
theta_nopool = y.copy()
# 部分プーリング：μ,τ を周辺尤度最大化（経験ベイズ）
def neg_logmarg(p):
    mu, log_tau = p; v = sigma**2 + np.exp(2*log_tau)
    return 0.5*np.sum(np.log(2*np.pi*v) + (y-mu)**2/v)
res = optimize.minimize(neg_logmarg, x0=[y.mean(), np.log(8)], method="Nelder-Mead")
mu_hat, tau_hat = res.x[0], np.exp(res.x[1])
theta_partial = (y/sigma**2 + mu_hat/tau_hat**2) / (1/sigma**2 + 1/tau_hat**2)

print(f"推定 μ={mu_hat:.2f}, τ={tau_hat:.2f}  / 完全プーリング={mu_complete:.2f}")
print(f"{'群':<4}{'σ_j':>5}{'非プーリングy_j':>14}{'部分プーリング':>14}{'収縮率B':>9}")
for j in range(J):
    B = sigma[j]**2/(sigma[j]**2 + tau_hat**2)
    print(f"{j+1:<4}{sigma[j]:>5.0f}{theta_nopool[j]:>14.2f}{theta_partial[j]:>14.2f}{B:>9.2f}")

出力：

推定 μ=49.39, τ=11.41  / 完全プーリング=47.02
群     σ_j     非プーリングy_j       部分プーリング     収縮率B
1       2         52.45         52.36     0.03
2       2         43.66         43.84     0.03
3       3         48.80         48.84     0.06
4       3         24.65         26.25     0.06
5       4         66.68         64.79     0.11
6       5         57.48         56.18     0.16
7       8         50.39         50.06     0.33
8      12         56.55         52.79     0.53

出力の意味：部分プーリングの推定は、非プーリング $y_j$ と全体平均 $\mu=49.4$ の間に必ず入ります。鍵は収縮率 $B_j$ がグループの不確実性で変わること。精密な群（ $\sigma=2$ 、群1・2）は $B\approx0.03$ でほとんど収縮せず $y_j$ のまま信頼されます。一方、粗い群（ $\sigma=12$ 、群8）は $B=0.53$ で全体平均へ半分以上引き戻され、 $56.55\to52.79$ に。データの乏しい群は「自分の値」より「全体の傾向」を借りるわけです。これにより、非プーリングなら暴れる小グループの推定が安定します。

3. 収縮の直観：情報を借り合う

部分プーリングが効く理由は、グループが互いに情報を貸し借りすることです。

全体平均 $\mu$ と群間ばらつき $\tau$ は全グループのデータから推定される。
各群はその「全体像」を事前として使い、自分のデータと混ぜる（精度の加重平均）。
データの少ない群ほど全体像に頼る（収縮が強い）、多い群ほど自前のデータを信じる。

結果として、非プーリングより外れにくく、完全プーリングより個性を保つ。この「いいとこ取り」が階層モデルの最大の価値です。収縮がなぜ推定を改善するのかの数理は収縮の数理で、 $\mu,\tau$ をどう推定するかは経験ベイズで深掘りします。実データでは $\theta_j,\mu,\tau$ を同時に MCMC（第4章）で推定するのが普通です。

非プーリング・部分プーリング・完全プーリングを並べて図示します。

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

rng = np.random.default_rng(2)
J = 8; sigma = np.array([2,2,3,3,4,5,8,12.0])
y = rng.normal(rng.normal(50.0, 10.0, J), sigma)
mu_complete = np.sum(y/sigma**2)/np.sum(1/sigma**2)
def neg_logmarg(p):
    v = sigma**2 + np.exp(2*p[1]); return 0.5*np.sum(np.log(2*np.pi*v)+(y-p[0])**2/v)
r = optimize.minimize(neg_logmarg, [y.mean(), np.log(8)], method="Nelder-Mead")
mu_hat, tau_hat = r.x[0], np.exp(r.x[1])
theta_partial = (y/sigma**2 + mu_hat/tau_hat**2)/(1/sigma**2 + 1/tau_hat**2)

g = np.arange(1, J+1)
plt.figure(figsize=(8, 4.5))
plt.errorbar(g, y, yerr=sigma, fmt="o", color="C0", label="非プーリング y_j ± σ_j", capsize=3)
plt.plot(g, theta_partial, "s", color="C1", label="部分プーリング（収縮後）")
plt.axhline(mu_complete, color="gray", ls="--", label=f"完全プーリング {mu_complete:.1f}")
for j in range(J):
    plt.plot([g[j], g[j]], [y[j], theta_partial[j]], color="C1", lw=0.8, alpha=0.6)
plt.xlabel("グループ"); plt.ylabel("推定値"); plt.legend(fontsize=9)
plt.title("収縮：不確実な群ほど全体平均へ引き戻される")
plt.tight_layout(); plt.show()

グラフの意味：青の点（非プーリング $y_j$ ）から橙の四角（部分プーリング）への縦線が収縮の量。誤差棒 $\sigma_j$ の大きい群（右ほど大きい）ほど縦線が長く、灰色の破線（完全プーリング）へ強く引き戻されます。誤差棒の小さい群はほとんど動きません。

⚠️ よくある誤解

「完全プーリングか非プーリングの二択」ではない：階層モデルが連続的な中間（部分プーリング）を $\tau$ で自動調整します。
「収縮はデータを歪める」ではない：収縮はむしろ平均二乗誤差を下げます（収縮の数理のジェームズ–スタイン）。特に小グループで有効。
「全グループ同じだけ収縮する」ではない：収縮率 $B_j=\sigma_j^2/(\sigma_j^2+\tau^2)$ はグループの不確実性で変わります。精密な群は収縮せず、粗い群は強く収縮。
「 $\tau$ は適当でよい」ではない： $\tau$ が収縮の強さを決める要。データから推定します（経験ベイズ＝経験ベイズ、完全ベイズなら MCMC）。

要点（BLUF）

1. グループ構造のデータと3つの選択肢

2. 部分プーリング＝階層モデル：親分布で結ぶ

3. 収縮の直観：情報を借り合う

⚠️ よくある誤解

関連ノート