信用区間と事後予測分布

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

📎 前提：点推定と損失関数　|　数理：事前分布・事後分布・ベイズ更新（統計）

要点（BLUF）

信用区間は「パラメータ $\theta$ がこの区間に入る確率が 95%」とそのまま読める区間。両裾を均等に切る等裾区間と、最も密度の高い範囲を取る HPD区間があります。
事後予測分布は、パラメータの不確実性を積分して次の観測を予測したもの。1点の推定値で予測するより幅が広くなり、それが正直な予測です。
ベータ–二項では、信用区間は $\mathrm{Beta}$ の分位点、事後予測は ベータ二項分布で閉じた形になります。

1. 信用区間：等裾とHPD

事後分布から、 $\theta$ が入る確率が $95\%$ になる区間を取ります。

P(\ell \le \theta \le u \mid D) = 0.95

取り方は一意ではなく、代表は2つ：

等裾区間（equal-tailed）：左右の裾をそれぞれ $2.5\%$ ずつ切る。分位点で簡単に求まる。
HPD区間（最高密度区間, highest posterior density）：区間内の密度が区間外より常に高くなるよう取る。同じ確率を最短の幅で覆える。

事後が対称なら両者は一致し、非対称なら HPD のほうが狭くなります。HPD は分位点だけでは出せず数値計算が要るので（arviz.hdi などは第4章で使います）、ここではまず等裾で実装します。

2. 事後予測分布：不確実性を積分して予測する

次に得るデータ $D_{\text{new}}$ を予測したいとき、「 $\theta$ を1点に決めて予測」では $\theta$ の不確実性を捨ててしまいます。ベイズでは $\theta$ を事後で積分して消します。

P(D_{\text{new}}\mid D) = \int P(D_{\text{new}}\mid\theta)\,P(\theta\mid D)\,d\theta

要するに：「ありうる $\theta$ すべてについて予測し、事後の重みで平均する」。だから事後予測分布は、データ自体のばらつき（二項のランダムさ）にパラメータの不確実性が上乗せされ、点推定での予測より広がります。ベータ事後 × 二項尤度の場合、この積分は閉じてベータ二項分布になります。

P(k'\mid D) = \int \mathrm{Binom}(k'\mid n', \theta)\,\mathrm{Beta}(\theta\mid a,b)\,d\theta = \mathrm{BetaBinomial}(k'\mid n', a, b)

flowchart LR
  PO["事後 P(θ|D)"] --> INT["θ で積分（重み付き平均）"]
  LK["観測モデル P(D_new|θ)"] --> INT
  INT --> PP["事後予測 P(D_new|D)"]

3. コードで信用区間と事後予測を求める

30球中20表 → 事後 $\mathrm{Beta}(21,11)$ で、信用区間と「次の10回の表数」の予測分布を出します。

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

# 30球中20表 → 事後 Beta(21,11)
post = stats.beta(21, 11)

# --- 95% 信用区間（等裾） ---
lo, hi = post.interval(0.95)
print(f"事後平均 = {post.mean():.3f}")
print(f"95%信用区間(等裾) = [{lo:.3f}, {hi:.3f}]")

# --- 事後予測分布：次に10回投げたときの表の数 k' ---
# ∫ Binom(k'|10,θ) Beta(θ|21,11) dθ = BetaBinomial(k'|10,21,11)
n_new = 10
k_new = np.arange(0, n_new + 1)
pred = stats.betabinom(n_new, 21, 11).pmf(k_new)
print(f"事後予測の平均 k' = {np.sum(k_new * pred):.2f}")
print(f"P(k' >= 8) = {pred[k_new >= 8].sum():.3f}")

fig, ax = plt.subplots(1, 2, figsize=(11, 4))
theta = np.linspace(0, 1, 300)
ax[0].plot(theta, post.pdf(theta), lw=2)
ax[0].axvspan(lo, hi, alpha=0.2, color="C0", label="95%信用区間")
ax[0].set_title("事後 Beta(21,11)"); ax[0].set_xlabel("θ"); ax[0].legend()
ax[1].bar(k_new, pred, color="C1")
ax[1].set_title("事後予測：次の10回の表の数 k'")
ax[1].set_xlabel("k'"); ax[1].set_ylabel("確率")
plt.tight_layout(); plt.show()

出力：

事後平均 = 0.656
95%信用区間(等裾) = [0.486, 0.808]
事後予測の平均 k' = 6.56
P(k' >= 8) = 0.308

出力の意味： $\theta$ の 95%信用区間は $[0.486, 0.808]$ ——「表の確率がこの範囲に入る確率が 95%」とそのまま読めます。右の事後予測分布は「次に10回投げたら表が何回出るか」で、平均 6.56 回、最頻は7回。そして 「8回以上表が出る確率は 0.308」 のように、未来の観測を確率で言えます。もし $\theta$ を点推定 $0.656$ に固定して二項分布で予測すると、この $\theta$ の不確実性ぶんだけ予測が過信（幅が狭く）になります。事後予測はそれを正しく広げています。

⚠️ よくある誤解

「信用区間＝信頼区間」ではない：信用区間は「 $\theta$ が入る確率」を直接表します。頻度論の信頼区間は「区間を作る手続きが 95% の確率で真値を覆う」という別物（頻度論とベイズの違い）。
「点推定で予測すれば十分」ではない： $\theta$ を1点に固定した予測は不確実性を捨て、予測区間が狭すぎる（過信）になります。意思決定に渡すなら事後予測を使います。
「等裾区間が唯一の信用区間」ではない：非対称な事後では HPD のほうが短く、極端な場合は等裾区間が密度の低い値を含むこともあります。目的に応じて選びます。

まとめ（Phase 1）

Phase 1 で、ベイズ推論の骨格——事後 ∝ 尤度 × 事前（事前・尤度・事後・周辺尤度）、データで信念を更新する逐次性（ベイズ更新と逐次推論）、損失で決まる点推定（点推定と損失関数）、そして区間と予測（本ノート）——が揃いました。次章では、ここで数値積分した事後を閉じた式で一気に求める共役事前分布に進みます。