事後予測チェック｜ベイズ統計

← ベイズ統計一覧

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

📎 前提：階層回帰と混合効果モデル　|　関連：信用区間と事後予測分布・事前の選び方と感度分析

要点（BLUF）

事後予測チェック（PPC）は、事後予測分布から複製データを何組も生成し、観測データと検定統計量で比べて、モデルが観測の特徴を再現できるか診断します。
ベイズ p 値は、複製の統計量が観測値以上になる割合。 $0.5$ 付近なら適合、 $0$ や $1$ 付近なら不適合のサイン。
過分散データにポアソンを当てると、PPC が分散の不適合を検出（ベイズ p 値 $0.000$ ）。事前予測（事前の選び方と感度分析）が事前の点検なら、PPC は事後の点検です。

1. 事後予測チェックの考え方

モデルが妥当なら、そのモデルから生成した「ありえたデータ」は、実際の観測データと似た特徴を持つはずです。PPC はこれを使います。

事後 $p(\theta\mid D)$ からパラメータをサンプル。
そのパラメータで複製データ $D^{\text{rep}}$ （観測と同じサイズ）を生成（信用区間と事後予測分布の事後予測分布から引く）。
これを何組も繰り返し、複製データの検定統計量 $T(D^{\text{rep}})$ の分布を、観測の $T(D)$ と比べる。

観測 $T(D)$ が複製の分布の中心付近にあれば「モデルはこの特徴を再現できている」、裾に外れていれば「再現できていない＝モデルの欠陥」。

2. 検定統計量とベイズ p 値

ずれを数値化するのが**事後予測 p 値（ベイズ p 値）**です。

p_B=\Pr\big(T(D^{\text{rep}})\ge T(D)\,\big|\,D\big)\approx\frac{1}{M}\sum_{m}\mathbb 1\big[T(D^{\text{rep}}_m)\ge T(D)\big]

$p_B\approx0.5$ なら観測は複製の真ん中＝適合。 $p_B$ が $0$ や $1$ に近いと、観測がモデルの予測の外＝不適合。頻度論の p 値と違い、「有意・棄却」の道具ではなく、モデルのどの側面が合っていないかを探す診断に使います。

3. コード：ポアソンの過分散を PPC で検出

過分散（分散 ≫ 平均）のカウントデータにポアソン回帰（分散＝平均を仮定）を当て、PPC が見抜くか確かめます。対照として、真にポアソンなデータも調べます。

import numpy as np
from scipy import stats

def ppc_poisson(data, label, seed=0):
    rng = np.random.default_rng(seed)
    n = len(data); S = data.sum()
    post = stats.gamma(a=1+S, scale=1/(1+n))         # ポアソン事後 λ~Gamma(1+Σx,1+n)
    T_var, T_zero, T_max = [], [], []
    for _ in range(4000):                            # 事後予測から複製データを生成
        lam = post.rvs(random_state=rng)
        rep = rng.poisson(lam, n)
        T_var.append(rep.var()); T_zero.append(np.mean(rep==0)); T_max.append(rep.max())
    T_var, T_zero, T_max = map(np.array, (T_var, T_zero, T_max))
    pv = lambda rep, obs: np.mean(rep >= obs)        # ベイズ p 値
    print(f"[{label}] 観測 平均={data.mean():.2f} 分散={data.var():.2f} "
          f"ゼロ率={np.mean(data==0):.2f} 最大={data.max()}")
    print(f"   ベイズp値: 分散={pv(T_var,data.var()):.3f} "
          f"ゼロ率={pv(T_zero,np.mean(data==0)):.3f} 最大={pv(T_max,data.max()):.3f}\n")

rng = np.random.default_rng(7)
over = stats.nbinom(n=1.5, p=1/3).rvs(120, random_state=rng)   # 過分散（平均3,分散9）
ppc_poisson(over, "ポアソンを過分散データに")
ppc_poisson(rng.poisson(3.0, 120), "ポアソンをポアソンデータに（対照）")

出力：

[ポアソンを過分散データに] 観測 平均=3.07 分散=7.50 ゼロ率=0.16 最大=13
   ベイズp値: 分散=0.000 ゼロ率=0.000 最大=0.003

[ポアソンをポアソンデータに（対照）] 観測 平均=2.77 分散=2.44 ゼロ率=0.06 最大=7
   ベイズp値: 分散=0.752 ゼロ率=0.626 最大=0.927

出力の意味：過分散データにポアソンを当てると、分散のベイズ p 値が $0.000$ ——観測の分散 $7.5$ は、ポアソン複製が作る分散（平均≒3 付近）よりはるかに大きく、複製の分布の右端を完全に外れています。ゼロ率・最大値も $0$ 付近で、「ポアソンでは過分散・多いゼロ・大きい外れ値を再現できない」と明確に検出。対照（真にポアソン）では p 値が $0.5$ 前後で適合。PPC は「平均は合うがばらつきの構造が合わない」というモデルの欠陥をピンポイントで指します（処方は負の二項回帰など）。

4. 図：複製の分散分布と観測のずれ

過分散ケースで、複製データの分散のヒストグラムに観測の分散を重ねます。

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

rng = np.random.default_rng(7)
data = stats.nbinom(n=1.5, p=1/3).rvs(120, random_state=rng)
post = stats.gamma(a=1+data.sum(), scale=1/(1+len(data)))
T_var = np.array([rng.poisson(post.rvs(random_state=rng), len(data)).var() for _ in range(4000)])

plt.figure(figsize=(7, 4))
plt.hist(T_var, bins=40, alpha=0.7, label="複製データの分散（ポアソン）")
plt.axvline(data.var(), color="crimson", lw=2, label=f"観測の分散 {data.var():.1f}")
plt.xlabel("分散の検定統計量"); plt.ylabel("頻度")
plt.title("事後予測チェック：観測の分散が複製の外＝不適合")
plt.legend(); plt.tight_layout(); plt.show()

グラフの意味：複製データの分散は $2$ 〜 $4$ あたりに固まり（ポアソンは分散＝平均だから）、観測の分散 $7.5$ （赤線）はその塊から大きく右に外れます。この「赤線が山の外」がベイズ p 値 $0.000$ の正体で、モデル不適合を一目で示します。

⚠️ よくある誤解

「ベイズ p 値は仮説検定の p 値」ではない：棄却の道具ではなく、どの側面が合わないかを探す診断。 $0.5$ から外れるほど不適合の手がかり。
「同じデータでフィットしてチェックは二重使用で無効」ではない：PPC は二重使用ぶん**保守的（不適合を見逃しやすい）**になりますが、それでも検出できれば確かな欠陥。厳密には交差検証（交差検証とLOO）を併用。
「平均が合えばモデルは正しい」ではない：分散・ゼロ率・最大など複数の統計量で見ます。平均だけ合って分散が外れる例が§3。
「PPC が通れば真のモデル」ではない：選んだ統計量で矛盾が出なかっただけ。別の統計量では破れるかもしれません。