ある学区の10の中学校で数学の試験を行ったとき、各学校の「本当の平均点」をどう推定すればよいでしょうか。素朴には各学校の標本平均を使えば良さそうですが、生徒数が少ない学校の推定はばらつきが大きくなります。一方、全学校のデータをプールして一つの平均を推定するのは、学校間の差を無視してしまいます。
この「個別推定とプール推定の間のちょうどよいバランス」を自動的に実現するのが階層ベイズモデル（hierarchical Bayesian model）です。各学校のパラメータが共通の分布（超母集団分布）から生成されるという構造を仮定することで、情報の適切な「借り入れ」（borrowing of strength）が実現されます。
階層ベイズモデルを理解すると、以下のような応用が可能になります。
- 教育統計: 学校・教師・学区レベルの効果の同時推定
- 臨床試験: 多施設試験における施設間変動の考慮
- スポーツ分析: 選手のパフォーマンスの推定（少ない出場機会のデータからの推定）
- 自然言語処理: トピックモデル（LDAなど）の構造
本記事の内容
- 階層ベイズモデルの構造と動機
- 縮小推定（shrinkage）の直感と理論
- 正規-正規モデルでの具体的な計算
- 経験ベイズ法との比較
- MCMCによるパラメータ推定
- Pythonによる実装と可視化
前提知識
階層ベイズモデルの構造
なぜ「階層」が必要なのか
統計的推定においてしばしば直面するジレンマがあります。たとえば、10の中学校の数学の平均点を推定する問題を考えましょう。
アプローチ1: 完全な個別推定 — 各学校のデータだけを使ってそれぞれの平均点を推定します。生徒が100人いる学校では精度の良い推定ができますが、生徒が5人しかいない学校では推定がばらつきます。極端な例では、たまたま成績の良い5人が集まった学校の平均点が過大評価されてしまいます。
アプローチ2: 完全なプール推定 — 全学校のデータをまとめて一つの平均を計算します。推定の精度は上がりますが、学校間の個性（学力レベルの違い）が完全に無視されます。
どちらのアプローチにも問題があります。直感的に望ましいのは、データが多い学校はそのデータを信頼し、データが少ない学校は全体の傾向を参考にする、というバランスです。階層ベイズモデルは、このバランスを数学的に最適な形で実現するのです。
この問題は野球の打率推定でも同じです。シーズン序盤で打席数が少ない選手の打率は大きくばらつきますが、「この選手はリーグの平均的な打者に似ているだろう」という事前知識を使えば、極端な推定を避けられます。打席が増えるにつれて、推定は実際のデータに近づいていきます。これがまさに階層ベイズの「情報の借り入れ」です。
2段階の確率モデル
階層ベイズモデルは、パラメータ自体が確率分布に従うという2段階（以上）の確率構造を持ちます。この構造を段階ごとに説明しましょう。
第1段階（データモデル）: 各グループ $j$ のデータ $y_{ij}$ はグループ固有のパラメータ $\theta_j$ に従います。
$$ y_{ij} | \theta_j \sim f(y | \theta_j), \quad i = 1, \dots, n_j $$
これは「各学校にはそれぞれ固有の真の平均点 $\theta_j$ があり、個々の生徒の成績はその周りにばらつく」ということを表しています。
第2段階（事前分布）: グループパラメータ $\theta_j$ は共通のハイパーパラメータ $\phi$ を持つ分布から生成されます。
$$ \theta_j | \phi \sim g(\theta | \phi), \quad j = 1, \dots, J $$
これは「各学校の真の平均点は、学区全体を特徴づける分布からサンプルされたものである」という仮定です。$\phi$ は分布の中心（学区全体の平均レベル）や広がり（学校間の差の大きさ）を制御します。
第3段階（ハイパー事前分布）: ハイパーパラメータ $\phi$ にも事前分布を置きます。
$$ \phi \sim h(\phi) $$
学区全体の平均レベルや学校間の差がどの程度であるかについて、弱い事前情報を導入します。
この多段階構造が「階層」の名前の由来です。パラメータのパラメータ（ハイパーパラメータ）にさらに分布を置くという入れ子構造が、モデルに柔軟性と頑健性を与えます。
階層モデルのグラフィカル表現
階層ベイズモデルの構造を理解するためには、グラフィカルモデル（有向非巡回グラフ、DAG）で表現すると分かりやすいです。ノードは確率変数を、矢印は条件付き依存関係を示します。
学校の例では、次のような依存構造になります。
$$ \phi \to \theta_1, \theta_2, \dots, \theta_J \to y_{11}, y_{12}, \dots, y_{Jn_J} $$
ハイパーパラメータ $\phi$ が各グループパラメータ $\theta_j$ に影響を与え、各 $\theta_j$ がそのグループのデータ $y_{ij}$ を生成します。重要なのは、$\theta_j$ たちは $\phi$ が与えられたもとで条件付き独立であるということです。しかし、$\phi$ を周辺化した（積分で除いた）場合、$\theta_j$ たちの間には$\phi$を介した間接的な依存関係が生じます。この依存関係が「情報の借り入れ」の源です。
この構造を数式で表すと、同時事後分布は次のようになります。
$$ p(\theta_1, \dots, \theta_J, \phi | \bm{y}) \propto \left[\prod_{j=1}^J \prod_{i=1}^{n_j} f(y_{ij}|\theta_j)\right] \left[\prod_{j=1}^J g(\theta_j|\phi)\right] h(\phi) $$
この同時分布の構造が、MCMCサンプリングや解析的な計算の出発点になります。
次に、階層ベイズモデルの最も重要な帰結である縮小推定を詳しく見ていきましょう。
縮小推定の直感と理論
縮小推定とは
階層ベイズモデルの最も重要な帰結は縮小推定（shrinkage estimation）です。各グループの推定値 $\hat{\theta}_j$ は、そのグループのデータから得られる推定値（個別推定）と全体平均の間のどこかに「縮小」されます。
この「縮小」の量はデータの量に依存します。
- データが少ないグループ: 全体平均に強く引き寄せられます。データが少ないので信頼性が低く、「他のグループと似ているだろう」という事前の信念が支配的になります。
- データが多いグループ: 個別推定に近い値を取ります。データが十分にあるので、そのグループ固有の情報が事前の信念に勝ります。
この適応的な重み付けが、階層ベイズモデルの最大の強みです。データの豊富さに応じて、全体の傾向と個別のデータの間のバランスを自動的に調整してくれるのです。
ジェームズ・スタインの推定量との関係
縮小推定の有用性は、統計学の歴史の中で衝撃的な発見として知られています。1961年にスタインとジェームズが証明したジェームズ・スタイン推定量（James-Stein estimator）は、3つ以上の正規分布の平均を同時に推定する場合、各グループの標本平均（最尤推定量）よりも全体を全体平均に向かって縮小した推定量の方が、平均二乗誤差の意味で常に優れていることを示しました。
この結果は「ジェームズ・スタインのパラドックス」と呼ばれ、当初は直感に反するものとして大きな議論を巻き起こしました。なぜ、野球選手の打率の推定が、気温の推定値から影響を受けるべきなのでしょうか。
階層ベイズモデルはこのパラドックスに自然な解釈を与えます。グループパラメータが共通の分布から生成されるという仮定のもとでは、縮小推定はベイズの定理の自然な帰結なのです。階層ベイズの観点では、ジェームズ・スタイン推定量は正規-正規階層モデルの事後平均の特殊ケースとして理解できます。
それでは、最も基本的な階層モデルである正規-正規モデルで、縮小推定の数学的な構造を見ていきましょう。
正規-正規階層モデル
モデルの定式化
最も基本的で分析的に扱いやすい階層モデルとして、正規-正規モデル（normal-normal hierarchical model）を考えます。
第1段階（データモデル）: 各グループ $j$ の標本平均 $\bar{y}_j$ は、真の平均 $\theta_j$ の周りに正規分布でばらつきます。
$$ \bar{y}_j | \theta_j \sim N(\theta_j, \sigma_j^2), \quad j = 1, \dots, J $$
ここで $\sigma_j^2 = \sigma^2/n_j$ は既知の分散です。$n_j$ はグループ $j$ のサンプルサイズで、サンプルが多いほど $\sigma_j^2$ は小さくなります。
第2段階（事前分布）: 各グループの真の平均 $\theta_j$ は、全体平均 $\mu$ の周りに正規分布でばらつきます。
$$ \theta_j | \mu, \tau \sim N(\mu, \tau^2) $$
$\mu$ は全体平均、$\tau^2$ はグループ間分散です。$\tau^2$ が小さければグループ間の差は小さく、大きければグループ間の差は大きいということを意味します。
このモデルの美しさは、事後分布が解析的に計算できることにあります。
事後分布の導出
$\mu$ と $\tau$ を固定したときの $\theta_j$ の事後分布を導出しましょう。ベイズの定理により、
$$ p(\theta_j | \bar{y}_j, \mu, \tau) \propto p(\bar{y}_j | \theta_j) \cdot p(\theta_j | \mu, \tau) $$
右辺は2つの正規分布の積です。尤度は $p(\bar{y}_j | \theta_j) \propto \exp\left(-\frac{(\bar{y}_j – \theta_j)^2}{2\sigma_j^2}\right)$ であり、事前分布は $p(\theta_j | \mu, \tau) \propto \exp\left(-\frac{(\theta_j – \mu)^2}{2\tau^2}\right)$ です。
指数部分を $\theta_j$ について平方完成します。まず指数部分の和を書き出すと、
$$ -\frac{(\bar{y}_j – \theta_j)^2}{2\sigma_j^2} – \frac{(\theta_j – \mu)^2}{2\tau^2} $$
$\theta_j$ の係数を整理するために展開すると、
$$ -\frac{1}{2}\left(\frac{1}{\sigma_j^2} + \frac{1}{\tau^2}\right)\theta_j^2 + \left(\frac{\bar{y}_j}{\sigma_j^2} + \frac{\mu}{\tau^2}\right)\theta_j + \text{const} $$
精度（分散の逆数）を導入すると見通しが良くなります。$\pi_j = 1/\sigma_j^2$（データの精度）、$\pi_0 = 1/\tau^2$（事前精度）とおくと、事後分布は正規分布になり、事後精度と事後平均は次のように計算されます。
事後精度: $\pi_j + \pi_0 = \frac{1}{\sigma_j^2} + \frac{1}{\tau^2}$
事後分散: $V_j = \frac{1}{\pi_j + \pi_0} = \frac{\sigma_j^2\tau^2}{\sigma_j^2 + \tau^2}$
事後平均: $\hat{\theta}_j = \frac{\pi_j \bar{y}_j + \pi_0 \mu}{\pi_j + \pi_0} = \frac{\bar{y}_j/\sigma_j^2 + \mu/\tau^2}{1/\sigma_j^2 + 1/\tau^2}$
この事後平均を、データの推定値 $\bar{y}_j$ と事前平均 $\mu$ の重み付き平均として書き直すと、
$$ \theta_j | \bar{y}_j, \mu, \tau \sim N\left(\hat{\theta}_j,\, V_j\right) $$
ここで事後平均と事後分散は
$$ \begin{equation} \hat{\theta}_j = (1 – B_j)\bar{y}_j + B_j \mu, \quad V_j = (1 – B_j)\sigma_j^2 \end{equation} $$
縮小係数 $B_j$ は
$$ \begin{equation} B_j = \frac{\sigma_j^2}{\sigma_j^2 + \tau^2} \end{equation} $$
この式が階層ベイズモデルの核心です。事後平均 $\hat{\theta}_j$ は、標本平均 $\bar{y}_j$ と全体平均 $\mu$ の加重平均であり、その重みが $B_j$ で制御されています。
縮小係数の解釈
$B_j$ はデータの分散 $\sigma_j^2$ とグループ間分散 $\tau^2$ の比で決まります。直感的に理解するために、極端な場合を考えましょう。
$\tau^2 \to \infty$ のとき（グループ間の差が非常に大きい場合）: $B_j = \sigma_j^2/(\sigma_j^2 + \tau^2) \to 0$ なので、$\hat{\theta}_j \to \bar{y}_j$（個別推定）に近づきます。グループ間の差が大きければ、他のグループの情報は参考にならないので、自分のデータだけを信頼するのが合理的です。
$\tau^2 \to 0$ のとき（グループ間の差がほとんどない場合）: $B_j \to 1$ なので、$\hat{\theta}_j \to \mu$（全体平均）に近づきます。グループ間の差がないのなら、すべてのデータをプールして推定するのが最善です。
$\sigma_j^2$ が小さいとき（データが多いグループ）: $B_j$ が小さくなり、個別推定 $\bar{y}_j$ を信頼します。データが十分にあるので、そのグループのデータから精度の高い推定ができます。
$\sigma_j^2$ が大きいとき（データが少ないグループ）: $B_j$ が大きくなり、全体平均 $\mu$ に縮小されます。データが少なく推定が不安定なので、他のグループの情報を借りて推定を安定化させます。
この適応的なメカニズムこそが「情報の借り入れ」（borrowing of strength）の正体です。データが豊富なグループはほとんど影響を受けませんが、データが乏しいグループは全体の傾向から大きな恩恵を受けるのです。
事後分散の解釈
事後分散 $V_j = (1 – B_j)\sigma_j^2$ も興味深い構造をしています。$0 < B_j < 1$ なので、常に $V_j < \sigma_j^2$ が成り立ちます。つまり、階層ベイズの推定量の分散は、個別推定の分散よりも小さくなります。
これは、事前情報（全体の傾向）を取り入れることで推定の不確実性が減少することを意味します。精度の向上はデータが少ないグループほど大きく、データの豊富なグループでは改善は限定的です。
次に、この理論的な結果をPythonで可視化し、縮小推定の効果を数値的に確認しましょう。
Pythonによる縮小推定の可視化
以下のコードでは、サンプルサイズが異なる10の学校のデータを生成し、個別推定と階層ベイズ推定を比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# --- 10校の模擬データ ---
J = 10 # 学校数
mu_true = 70 # 全体の真の平均
tau_true = 8 # グループ間の真の標準偏差
sigma = 15 # 個人の標準偏差
# 各学校の真のパラメータ
theta_true = np.random.normal(mu_true, tau_true, J)
# 各学校のサンプルサイズ（不均一）
n_j = np.array([5, 8, 10, 15, 20, 25, 30, 50, 80, 100])
# 各学校の標本平均
y_bar = np.array([np.random.normal(theta, sigma / np.sqrt(n))
for theta, n in zip(theta_true, n_j)])
sigma_j2 = sigma**2 / n_j # 各学校の分散
# --- 階層ベイズの推定（τを固定） ---
tau = tau_true # ここでは真の値を使用
mu = np.mean(y_bar) # 全体平均のプラグイン推定
B_j = sigma_j2 / (sigma_j2 + tau**2)
theta_hat = (1 - B_j) * y_bar + B_j * mu
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 縮小推定の可視化
ax = axes[0]
for j in range(J):
ax.plot([j, j], [y_bar[j], theta_hat[j]], "gray", linewidth=1, alpha=0.7)
ax.scatter(range(J), y_bar, s=80, color="red", marker="^", zorder=5,
label=r"$\bar{y}_j$ (sample mean)")
ax.scatter(range(J), theta_hat, s=80, color="blue", marker="o", zorder=5,
label=r"$\hat{\theta}_j$ (hierarchical)")
ax.scatter(range(J), theta_true, s=40, color="green", marker="D", zorder=5,
label=r"$\theta_j$ (true)")
ax.axhline(mu, color="purple", linewidth=1.5, linestyle="--",
label=rf"$\mu = {mu:.1f}$ (overall mean)")
ax.set_xlabel("School j", fontsize=12)
ax.set_ylabel("Score", fontsize=12)
ax.set_title("Shrinkage estimation", fontsize=13)
ax.legend(fontsize=9)
ax.set_xticks(range(J))
ax.set_xticklabels([f"n={n}" for n in n_j], fontsize=9, rotation=45)
ax.grid(True, alpha=0.3)
# (b) 縮小係数と推定精度
ax = axes[1]
mse_individual = (y_bar - theta_true)**2
mse_hierarchical = (theta_hat - theta_true)**2
x_pos = np.arange(J)
width = 0.35
ax.bar(x_pos - width/2, mse_individual, width, color="salmon", alpha=0.7,
label=rf"Individual MSE (avg={mse_individual.mean():.1f})")
ax.bar(x_pos + width/2, mse_hierarchical, width, color="steelblue", alpha=0.7,
label=rf"Hierarchical MSE (avg={mse_hierarchical.mean():.1f})")
ax.set_xlabel("School j", fontsize=12)
ax.set_ylabel("Squared error", fontsize=12)
ax.set_title("MSE comparison", fontsize=13)
ax.set_xticks(x_pos)
ax.set_xticklabels([f"n={n}" for n in n_j], fontsize=9, rotation=45)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("hierarchical_bayes.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"個別推定の平均MSE: {mse_individual.mean():.2f}")
print(f"階層ベイズの平均MSE: {mse_hierarchical.mean():.2f}")
print(f"改善率: {(1 - mse_hierarchical.mean()/mse_individual.mean())*100:.1f}%")
この可視化から、階層ベイズモデルの縮小推定の効果が明確に確認できます。
-
左図: 標本平均（赤の三角）が階層ベイズ推定値（青の丸）に比べて大きく散らばっている 。灰色の線は縮小の量を示しており、サンプルサイズが小さい学校（左側）ほど全体平均（紫の破線）に強く引き寄せられています。サンプルサイズが大きい学校（右側）ではほとんど縮小が起きず、標本平均に近い値を取っています。これはまさに縮小係数 $B_j$ の理論的予測と一致しています。
-
右図: ほとんどの学校で階層ベイズ推定（青）の方が個別推定（赤）よりも二乗誤差が小さい 。特にサンプルサイズが小さい学校で改善が顕著です。全体的な平均MSEも階層ベイズの方が大幅に小さく、情報の借り入れの効果が確認できます。サンプルサイズが大きい学校では差が小さいのは、$B_j$ が小さく縮小がほとんど起きないためです。
次に、縮小係数 $B_j$ の振る舞いをより詳しく可視化して、サンプルサイズとの関係を確認しましょう。
縮小係数の詳細分析
import numpy as np
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 縮小係数 B_j のサンプルサイズ依存性
ax = axes[0]
sigma = 15
n_range = np.linspace(1, 200, 500)
for tau in [2, 5, 8, 15, 30]:
sigma_j2 = sigma**2 / n_range
B = sigma_j2 / (sigma_j2 + tau**2)
ax.plot(n_range, B, linewidth=2, label=rf"$\tau = {tau}$")
ax.set_xlabel("Sample size $n_j$", fontsize=12)
ax.set_ylabel("Shrinkage factor $B_j$", fontsize=12)
ax.set_title(rf"Shrinkage factor vs sample size ($\sigma = {sigma}$)", fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 1)
ax.grid(True, alpha=0.3)
# (b) τの推定による影響
ax = axes[1]
tau_range = np.linspace(0.1, 30, 300)
n_vals = [5, 10, 30, 100]
for n_val in n_vals:
sigma_j2_val = sigma**2 / n_val
B_vals = sigma_j2_val / (sigma_j2_val + tau_range**2)
ax.plot(tau_range, B_vals, linewidth=2, label=rf"$n_j = {n_val}$")
ax.set_xlabel(r"Between-group SD $\tau$", fontsize=12)
ax.set_ylabel("Shrinkage factor $B_j$", fontsize=12)
ax.set_title(rf"Shrinkage factor vs $\tau$ ($\sigma = {sigma}$)", fontsize=13)
ax.legend(fontsize=10)
ax.set_ylim(0, 1)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("shrinkage_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
この分析から、縮小係数の振る舞いについて以下のことが明確になります。
-
左図: サンプルサイズが増えるほど縮小係数 $B_j$ は急速に減少する 。$\tau = 8$ の場合、$n_j = 5$ では $B_j \approx 0.78$（全体平均に強く縮小）ですが、$n_j = 100$ では $B_j \approx 0.03$（ほぼ個別推定）になります。$\tau$ が大きい（グループ間差が大きい）ほど、$B_j$ は全体的に小さく、個別推定を信頼する方向にシフトします。
-
右図: $\tau$ が大きくなるほど $B_j$ は減少する 。$\tau$ が小さい（グループ間差が小さい）場合、ほとんどのグループで縮小が強く効きます。これは直感に合っています。グループ間の差がほとんどないなら、全体平均が各グループの最良推定になるべきだからです。
経験ベイズ法との比較
経験ベイズの考え方
階層ベイズモデルでは $\mu$ と $\tau$ にも事前分布を置きますが、経験ベイズ法（empirical Bayes）ではこれらのハイパーパラメータをデータから推定し、推定値をプラグインします。
経験ベイズ法では、$\bar{y}_1, \dots, \bar{y}_J$ の周辺分布 $\bar{y}_j \sim N(\mu, \sigma_j^2 + \tau^2)$ を用いて $\mu$ と $\tau$ を推定します。具体的には、
$$ \hat{\mu} = \frac{\sum_j w_j \bar{y}_j}{\sum_j w_j}, \quad w_j = \frac{1}{\sigma_j^2 + \hat{\tau}^2} $$
$\hat{\tau}^2$ は周辺尤度を最大化する値として求められます。
階層ベイズとの違い
経験ベイズ法は計算が簡便ですが、ハイパーパラメータの推定における不確実性を無視するという欠点があります。$\hat{\mu}$ と $\hat{\tau}$ を固定値として扱うため、$\theta_j$ の事後分布が実際よりも狭く（不確実性が過小評価）なる傾向があります。
完全ベイズ法（fully Bayesian approach）では、$\mu$ と $\tau$ の事後分布上で積分を行い、ハイパーパラメータの不確実性も考慮します。
$$ p(\theta_j | \bm{y}) = \int p(\theta_j | \bar{y}_j, \mu, \tau) \cdot p(\mu, \tau | \bm{y})\, d\mu\, d\tau $$
この積分は一般に解析的には計算できないため、MCMCのような数値的手法が必要になります。
数値比較
両者の違いをコードで確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats, optimize
np.random.seed(42)
# データ生成
J = 10
mu_true = 70
tau_true = 8
sigma = 15
theta_true = np.random.normal(mu_true, tau_true, J)
n_j = np.array([5, 8, 10, 15, 20, 25, 30, 50, 80, 100])
y_bar = np.array([np.random.normal(theta, sigma / np.sqrt(n))
for theta, n in zip(theta_true, n_j)])
sigma_j2 = sigma**2 / n_j
# 経験ベイズ: τの周辺尤度推定
def neg_marginal_loglik(log_tau2, y_bar, sigma_j2):
tau2 = np.exp(log_tau2)
total_var = sigma_j2 + tau2
ll = -0.5 * np.sum(np.log(total_var) + (y_bar - np.mean(y_bar))**2 / total_var)
return -ll
result = optimize.minimize_scalar(lambda lt2: neg_marginal_loglik(lt2, y_bar, sigma_j2),
bounds=(-5, 10), method="bounded")
tau_eb = np.sqrt(np.exp(result.x))
mu_eb = np.average(y_bar, weights=1/(sigma_j2 + tau_eb**2))
# 経験ベイズ推定
B_eb = sigma_j2 / (sigma_j2 + tau_eb**2)
theta_eb = (1 - B_eb) * y_bar + B_eb * mu_eb
# 真のτを使った場合
B_true = sigma_j2 / (sigma_j2 + tau_true**2)
theta_true_shrink = (1 - B_true) * y_bar + B_true * mu_true
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 推定値の比較
ax = axes[0]
schools = np.arange(J)
ax.scatter(schools, y_bar, s=80, marker="^", color="red", zorder=5,
label="Sample mean")
ax.scatter(schools, theta_eb, s=80, marker="o", color="blue", zorder=5,
label=f"Empirical Bayes ($\\hat{{\\tau}}={tau_eb:.1f}$)")
ax.scatter(schools, theta_true_shrink, s=60, marker="s", color="green", zorder=5,
label=f"Oracle ($\\tau={tau_true}$)")
ax.scatter(schools, theta_true, s=40, marker="D", color="black", zorder=5,
label="True $\\theta_j$")
ax.set_xlabel("School j", fontsize=12)
ax.set_ylabel("Score", fontsize=12)
ax.set_title("Empirical Bayes vs Oracle", fontsize=13)
ax.legend(fontsize=9)
ax.set_xticks(schools)
ax.set_xticklabels([f"n={n}" for n in n_j], fontsize=9, rotation=45)
ax.grid(True, alpha=0.3)
# (b) MSE比較
ax = axes[1]
mse_ind = (y_bar - theta_true)**2
mse_eb = (theta_eb - theta_true)**2
mse_oracle = (theta_true_shrink - theta_true)**2
methods = ["Individual", f"Emp. Bayes\n(tau={tau_eb:.1f})", f"Oracle\n(tau={tau_true})"]
means = [mse_ind.mean(), mse_eb.mean(), mse_oracle.mean()]
colors = ["salmon", "steelblue", "lightgreen"]
ax.bar(methods, means, color=colors, edgecolor="gray", alpha=0.8)
ax.set_ylabel("Average MSE", fontsize=12)
ax.set_title("MSE comparison across methods", fontsize=13)
ax.grid(True, alpha=0.3, axis="y")
for i, v in enumerate(means):
ax.text(i, v + 1, f"{v:.1f}", ha="center", fontsize=11)
plt.tight_layout()
plt.savefig("empirical_bayes_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
この比較から、以下のことが確認できます。
-
左図: 経験ベイズ推定値（青）はオラクル推定値（緑、真の $\tau$ を使った場合）に近い 。データから $\tau$ を推定してもかなり良い推定が得られることがわかります。どちらも標本平均（赤）よりも真の値（黒）に近い傾向にあります。
-
右図: 経験ベイズ法は個別推定よりも大幅にMSEが小さい 。経験ベイズ法とオラクル法のMSEの差は比較的小さく、ハイパーパラメータの推定がうまくいっていることを示しています。実際にはオラクル法は利用できない（真の $\tau$ は未知）ので、経験ベイズ法は非常に実用的なアプローチです。
次に、ハイパーパラメータの不確実性も含めた完全ベイズ法をMCMCで実装してみましょう。
MCMCによる完全ベイズ推定
ギブスサンプリングの構造
正規-正規階層モデルでは、条件付き事後分布がすべて解析的に計算できるため、ギブスサンプリング（Gibbs sampling）が自然な選択です。
各パラメータの完全条件付き分布は次の通りです。
$\theta_j$ の条件付き分布（他のすべてのパラメータが固定された場合）:
$$ \theta_j | \bar{y}_j, \mu, \tau, \sigma \sim N\left((1-B_j)\bar{y}_j + B_j\mu,\; (1-B_j)\sigma_j^2\right) $$
$\mu$ の条件付き分布（$\theta_j$ と $\tau$ が固定された場合、無情報事前分布 $p(\mu) \propto 1$ を仮定）:
$$ \mu | \theta_1, \dots, \theta_J, \tau \sim N\left(\bar{\theta},\; \tau^2/J\right) $$
ここで $\bar{\theta} = \frac{1}{J}\sum_{j=1}^J \theta_j$ です。
$\tau^2$ の条件付き分布: 半コーシーなどの事前分布を使う場合はメトロポリス・ヘイスティングスが必要ですが、逆ガンマ事前分布 $\tau^2 \sim \text{Inv-Gamma}(a, b)$ を使えばギブスサンプリングが可能です。
$$ \tau^2 | \theta_1, \dots, \theta_J, \mu \sim \text{Inv-Gamma}\left(a + J/2,\; b + \frac{1}{2}\sum_{j=1}^J(\theta_j – \mu)^2\right) $$
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(42)
# データ（前のセクションと同じ）
J = 10
mu_true = 70
tau_true = 8
sigma = 15
theta_true = np.random.normal(mu_true, tau_true, J)
n_j = np.array([5, 8, 10, 15, 20, 25, 30, 50, 80, 100])
y_bar = np.array([np.random.normal(theta, sigma / np.sqrt(n))
for theta, n in zip(theta_true, n_j)])
sigma_j2 = sigma**2 / n_j
# --- ギブスサンプリング ---
n_iter = 10000
burn_in = 2000
# 保存用配列
theta_samples = np.zeros((n_iter, J))
mu_samples = np.zeros(n_iter)
tau2_samples = np.zeros(n_iter)
# 初期値
mu_current = np.mean(y_bar)
tau2_current = np.var(y_bar)
theta_current = y_bar.copy()
# 事前分布のハイパーパラメータ（弱い情報の事前分布）
a_prior = 1.0
b_prior = 1.0
for t in range(n_iter):
# Step 1: θ_j をサンプル
for j in range(J):
B_j = sigma_j2[j] / (sigma_j2[j] + tau2_current)
post_mean = (1 - B_j) * y_bar[j] + B_j * mu_current
post_var = (1 - B_j) * sigma_j2[j]
theta_current[j] = np.random.normal(post_mean, np.sqrt(post_var))
# Step 2: μ をサンプル
mu_post_mean = np.mean(theta_current)
mu_post_var = tau2_current / J
mu_current = np.random.normal(mu_post_mean, np.sqrt(mu_post_var))
# Step 3: τ² をサンプル（逆ガンマ事前分布）
a_post = a_prior + J / 2
b_post = b_prior + 0.5 * np.sum((theta_current - mu_current)**2)
tau2_current = 1.0 / np.random.gamma(a_post, 1.0 / b_post)
# 保存
theta_samples[t] = theta_current
mu_samples[t] = mu_current
tau2_samples[t] = tau2_current
# バーンインを除去
theta_post = theta_samples[burn_in:]
mu_post = mu_samples[burn_in:]
tau2_post = tau2_samples[burn_in:]
# 結果の可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
# (a) μのトレースプロットと事後分布
ax = axes[0, 0]
ax.plot(mu_samples, linewidth=0.3, alpha=0.7, color="steelblue")
ax.axhline(mu_true, color="red", linewidth=2, linestyle="--", label=f"True $\\mu = {mu_true}$")
ax.axvline(burn_in, color="gray", linewidth=1, linestyle=":", label="Burn-in")
ax.set_xlabel("Iteration", fontsize=12)
ax.set_ylabel(r"$\mu$", fontsize=12)
ax.set_title(r"Trace plot of $\mu$", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (b) τのトレースプロットと事後分布
ax = axes[0, 1]
ax.plot(np.sqrt(tau2_samples), linewidth=0.3, alpha=0.7, color="steelblue")
ax.axhline(tau_true, color="red", linewidth=2, linestyle="--", label=f"True $\\tau = {tau_true}$")
ax.axvline(burn_in, color="gray", linewidth=1, linestyle=":", label="Burn-in")
ax.set_xlabel("Iteration", fontsize=12)
ax.set_ylabel(r"$\tau$", fontsize=12)
ax.set_title(r"Trace plot of $\tau$", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (c) θ_j の事後分布
ax = axes[1, 0]
theta_post_means = theta_post.mean(axis=0)
theta_post_std = theta_post.std(axis=0)
ax.errorbar(range(J), theta_post_means, yerr=2*theta_post_std, fmt="o",
color="blue", capsize=5, label="Posterior mean ± 2SD")
ax.scatter(range(J), theta_true, s=80, color="green", marker="D", zorder=5,
label=r"True $\theta_j$")
ax.scatter(range(J), y_bar, s=60, color="red", marker="^", zorder=5,
label=r"Sample mean $\bar{y}_j$")
ax.axhline(mu_post.mean(), color="purple", linewidth=1, linestyle="--",
label=rf"Post. mean of $\mu = {mu_post.mean():.1f}$")
ax.set_xlabel("School j", fontsize=12)
ax.set_ylabel("Score", fontsize=12)
ax.set_title("Posterior estimates with uncertainty", fontsize=13)
ax.legend(fontsize=8)
ax.set_xticks(range(J))
ax.set_xticklabels([f"n={n}" for n in n_j], fontsize=9, rotation=45)
ax.grid(True, alpha=0.3)
# (d) ハイパーパラメータの事後分布
ax = axes[1, 1]
ax.hist(mu_post, bins=50, density=True, alpha=0.5, color="steelblue",
edgecolor="white", label=rf"$\mu$ (mean={mu_post.mean():.1f})")
ax.hist(np.sqrt(tau2_post), bins=50, density=True, alpha=0.5, color="salmon",
edgecolor="white", label=rf"$\tau$ (mean={np.sqrt(tau2_post).mean():.1f})")
ax.axvline(mu_true, color="blue", linewidth=2, linestyle="--")
ax.axvline(tau_true, color="red", linewidth=2, linestyle="--")
ax.set_xlabel("Value", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("Posterior distributions of hyperparameters", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("hierarchical_bayes_mcmc.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"μ の事後平均: {mu_post.mean():.2f} (真の値: {mu_true})")
print(f"τ の事後平均: {np.sqrt(tau2_post).mean():.2f} (真の値: {tau_true})")
このMCMCの結果から、完全ベイズ推定の特徴が確認できます。
-
上段のトレースプロット: $\mu$ と $\tau$ のサンプルがバーンイン後に安定した定常分布に達している 。真の値（赤の破線）の周辺でサンプルが分布しており、MCMCが正しく機能していることがわかります。
-
左下: 事後分布の平均（青の丸）は標本平均（赤の三角）よりも真の値（緑のひし形）に近い 。特にサンプルサイズが小さい学校（左側）で改善が顕著です。エラーバー（事後標準偏差の2倍）はサンプルサイズが小さい学校ほど広く、推定の不確実性を適切に反映しています。
-
右下: $\mu$ と $\tau$ の事後分布が真の値の近くにピークを持っている 。事後分布の幅はハイパーパラメータの推定に残る不確実性を表しています。経験ベイズ法ではこの不確実性が無視されるため、完全ベイズ法の方がより保守的で正直な推論を提供します。
3つ以上の階層
ここまでは2段階（データ → グループパラメータ → ハイパーパラメータ）の階層モデルを扱いましたが、実際の問題ではより多くの階層が必要になることがあります。
たとえば、教育データでは生徒 → クラス → 学校 → 学区という4段階の階層構造が自然です。各レベルで分散が定義され、「情報の借り入れ」は各レベルで起こります。
生徒 $i$ がクラス $c$、学校 $s$、学区 $d$ に属しているとき、
$$ y_{icsd} = \mu + \alpha_d + \beta_{s(d)} + \gamma_{c(s,d)} + \epsilon_{icsd} $$
$$ \alpha_d \sim N(0, \tau_d^2), \quad \beta_{s(d)} \sim N(0, \tau_s^2), \quad \gamma_{c(s,d)} \sim N(0, \tau_c^2) $$
このようなモデルはマルチレベルモデル（multilevel model）とも呼ばれ、社会科学やバイオインフォマティクスで広く使われています。階層が増えるほど解析的な計算は困難になりますが、MCMCやハミルトニアンモンテカルロ法（HMC）を使えばサンプリングが可能です。
階層ベイズモデルの設計指針
階層ベイズモデルを実際に構築する際には、以下の点に注意が必要です。
グループ数: グループ数 $J$ が小さすぎると（例えば $J < 5$）、ハイパーパラメータ $\tau$ の推定が不安定になります。$J$ が大きいほど、$\tau$ の推定精度が上がり、縮小推定の恩恵も大きくなります。
$\tau$ の事前分布: 完全ベイズ推定では、$\tau$（またはτ²）の事前分布の選択が重要です。逆ガンマ分布 $\text{Inv-Gamma}(\epsilon, \epsilon)$ はかつてよく使われましたが、$\epsilon \to 0$ の極限で不適切な事後分布になる場合があります。現在では、半コーシー分布 $\tau \sim \text{Half-Cauchy}(0, A)$ や半正規分布が推奨されています。
収縮の方向: 正規-正規モデルでは全体平均に向かって縮小しますが、問題に応じて適切な「基準点」を選ぶことが重要です。回帰モデルでは、係数をゼロに向かって縮小する正則化ベイズ（ラッソのベイズ版など）が広く使われています。
まとめ
本記事では、階層ベイズモデルの構造、理論、実装を解説しました。
- 階層ベイズモデルは、パラメータが共通の分布から生成されるという多段階の確率構造を持ち、個別推定とプール推定の間のバランスを自動的に実現する
- 縮小推定: 各グループの推定値は個別推定と全体平均の間に「縮小」され、データが少ないグループほど全体平均に近づく
- 縮小係数 $B_j = \sigma_j^2/(\sigma_j^2 + \tau^2)$ がデータ量に応じた適応的な重み付けを実現し、個別推定に比べて平均的な推定精度（MSE）が向上する（情報の借り入れの効果）
- 経験ベイズ法はハイパーパラメータをデータから推定する簡便な方法だが、ハイパーパラメータの不確実性を無視する。完全ベイズ法ではMCMCを用いてすべてのパラメータの事後分布を同時に推定する
- ギブスサンプリングにより、正規-正規階層モデルの完全条件付き分布からの逐次サンプリングが効率的に実行できる
次のステップとして、以下の記事も参考にしてください。
- ベイズノンパラメトリクス入門 — 無限次元の階層モデル
- 経験ベイズ法 — ハイパーパラメータをデータから推定する方法
- ベイズモデル比較 — 階層モデルの選択基準
- MCMC入門 — 階層モデルの計算手法