階層モデルの実例と再パラメータ化

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：経験ベイズ　|　関連：ハミルトニアンモンテカルロとNUTS・収束診断

要点（BLUF）

階層モデルの事後は漏斗（funnel）型になりやすい：群間ばらつき $\tau$ が小さい領域で群効果 $\theta_j$ が狭く絞られ、サンプラーが首（neck）に入れません。
中心化パラメータ化（ $\theta_j$ を直接サンプル）は首で詰まります。非中心化（ $\theta_j=\mu+\tau z_j,\ z_j\sim\mathcal N(0,1)$ ）にすると幾何が真っ直ぐになり、よく混ざります。
Neal の漏斗で再現すると、中心化は分布の裾を取りこぼし（ $\mathrm{std}$ 過小）、非中心化が真の周辺 $\mathrm{std}=3$ を回収します。

1. 漏斗の幾何：なぜ首ができるのか

階層モデル $\theta_j\sim\mathcal N(\mu,\tau^2)$ では、群間ばらつき $\tau$ も推定対象でした（階層モデルの構造）。 $\tau$ と $\theta_j$ の同時事後を見ると、 $\tau$ が小さいほど $\theta_j$ は $\mu$ の周りに強く絞られ、 $\tau$ が大きいほど広く散らばる。この「 $\tau$ が下がると幅がすぼまる」形が、横から見ると漏斗になります。首の部分（小さい $\tau$ ・狭い $\theta$ ）は曲率が極端に強く、一定歩幅のサンプラーは——大きい歩幅では弾かれ、小さい歩幅では進めず——身動きが取れなくなります。

この病理を抽出したのが Neal の漏斗： $v\sim\mathcal N(0,3^2)$ 、 $x_i\sim\mathcal N(0,e^{v})$ 。 $v$ が $\log\tau^2$ 、 $x_i$ が群効果 $\theta_j$ の役回りです。 $v$ が負（ $\tau$ 小）だと $x_i$ は極端に狭く、首ができます。

2・3. 中心化の失敗と非中心化の処方

中心化：パラメータを $(v, x_1,\dots,x_D)$ のまま扱う。 $x_i$ の幅が $v$ に依存するので、首で詰まる。
非中心化： $x_i=e^{v/2}z_i,\ z_i\sim\mathcal N(0,1)$ と書き換え、 $(v, z_1,\dots,z_D)$ をサンプル。すると $v$ と $z$ が独立になり、漏斗が消えて素直な形に。サンプリング後に $x_i=e^{v/2}z_i$ で戻します。

\underbrace{x_i\sim\mathcal N(0,e^{v})}_{\text{中心化：幅が }v\text{ 依存}}\quad\Longleftrightarrow\quad \underbrace{z_i\sim\mathcal N(0,1),\ x_i=e^{v/2}z_i}_{\text{非中心化：幅が一定}}

両者を自作 MH（メトロポリスヘイスティングス）で回し、 $v$ の周辺（真は $\mathcal N(0,3^2)$ 、 $\mathrm{std}=3$ ）を回収できるか比べます。

import numpy as np

D = 9                                       # Neal の漏斗：x を 9 本
def mh(logpost, x0, step, N=80000, burn=15000, seed=0):
    rng = np.random.default_rng(seed)
    x = np.array(x0, float); lp = logpost(x); out = np.empty((N, len(x))); acc = 0
    for i in range(N):
        cand = x + rng.normal(0, step, len(x)); lc = logpost(cand)
        if np.log(rng.uniform()) < lc - lp: x, lp = cand, lc; acc += 1
        out[i] = x
    return out[burn:], acc/N

# 中心化：(v, x_1..x_D)。log p = -v²/18 - (D/2)v - ½Σx_i²/e^v
def logpost_centered(p):
    v, xs = p[0], p[1:]
    if v > 30: return -np.inf
    return -v**2/18 - 0.5*D*v - 0.5*np.sum(xs**2)/np.exp(v)

# 非中心化：(v, z_1..z_D)、x=e^{v/2}z → v と z は独立
def logpost_noncentered(p):
    v, z = p[0], p[1:]
    return -v**2/18 - 0.5*np.sum(z**2)

cen, ac = mh(logpost_centered,    [0.0]*(D+1), step=0.5)
non, an = mh(logpost_noncentered, [0.0]*(D+1), step=0.5)
print("真の周辺 v~N(0,3²)：std(v)=3.000 を回収したい")
print(f"中心化  : 受容率={ac:.2f}  std(v)={cen[:,0].std():.3f}  v最小={cen[:,0].min():.1f}")
print(f"非中心化: 受容率={an:.2f}  std(v)={non[:,0].std():.3f}  v最小={non[:,0].min():.1f}")

出力：

真の周辺 v~N(0,3²)：std(v)=3.000 を回収したい
中心化  : 受容率=0.36  std(v)=2.115  v最小=-3.5
非中心化: 受容率=0.47  std(v)=2.887  v最小=-10.3

出力の意味：中心化は $v$ の $\mathrm{std}$ を $2.12$ と過小評価し、 $v$ の最小も $-3.5$ までしか届きません——首（ $v$ が負の領域）に入れず、分布の裾を丸ごと取りこぼしている。一方、非中心化は $\mathrm{std}=2.89\approx3$ を回収し、 $v=-10.3$ まで深く首に入れています。同じサンプラー・同じ歩幅でも、書き方（パラメータ化）を変えるだけで結果が一変する。これは近似誤差ではなく、中心化では原理的に届かない領域があるためです。

flowchart LR
  C["中心化 θ_j ~ N(μ, τ²)<br/>幅が τ に依存 → 漏斗の首で詰まる"] --> P["再パラメータ化"]
  P --> N["非中心化 θ_j = μ + τ·z_j<br/>z_j ~ N(0,1)・幅一定 → 首が消える"]

4. 実務：発散・診断・PPL

発散（divergence）：HMC/NUTS（ハミルトニアンモンテカルロとNUTS）は首の強い曲率で軌道が破綻し「発散」を報告します。発散が多発したらまず非中心化を疑うのが定石。
診断：中心化の失敗は、 $\hat R$ の悪化・ $\tau$ 付近の低 ESS・トレースの張り付きに表れます（収束診断）。1本のチェーンだけ見ると「収束した」と誤認しがち。
PPL での扱い：PyMC・Stan・NumPyro では pm.Normal("theta", mu, tau) を非中心化形（mu + tau * z、z = pm.Normal("z",0,1)）に書き換えるのが定番。ライブラリによっては非中心化のヘルパもあります（API・既定は更新が速く要最新確認。PPL 概観は確率的プログラミング概観）。
使い分け：データが多くグループが情報豊富なら中心化が良いこともあります（首が浅くなる）。データが少なく $\tau$ が小さくなりがちな問題ほど非中心化が効きます。

⚠️ よくある誤解

「サンプラーが悪いだけ」ではない：漏斗は事後の幾何の問題。再パラメータ化（モデルの書き換え）が本質的な処方で、サンプラー変更だけでは限界があります。
「中心化と非中心化は別モデル」ではない：同じ事後の書き方違い。変数変換で完全に等価です（戻せば一致）。
「非中心化が常に最良」ではない：データが豊富で群効果が強い場合は中心化が有利なことも。問題に応じて選びます。
「 $\hat R\approx1$ なら漏斗は大丈夫」ではない：全チェーンが首を避ければ $\hat R$ は良く見えても裾を取りこぼします。発散の有無・ $\tau$ 付近の ESS も確認（収束診断）。

まとめ（Phase 5）

第5章では、グループ構造を持つデータをまとめて扱う階層ベイズを学びました——3つのプーリングと部分プーリングの収縮（階層モデルの構造）、収縮がなぜ推定を改善するかのジェームズ–スタイン（収縮の数理）、事前をデータから決める経験ベイズ（経験ベイズ）、そして実装の落とし穴と非中心化（本ノート）。階層モデルは第4章の MCMC を計算エンジンに、第2章の共役を条件付きの部品にして組み上がる、ベイズの総合力が問われる領域です。次章では、事後を近似でもっと速く求める変分推論へ進みます。