正規正規モデル｜ベイズ統計

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

📎 前提：ガンマポアソンと指数　|　数理：共役事前分布（統計）・正規分布（標準正規・標準化）（統計）

要点（BLUF）

正規（分散既知）の平均 $\mu$ の共役は正規。精度 $\tau=1/\sigma^2$ で書くと、事後精度＝精度の足し算 $\tau_n=\tau_0+n\tau$ 、事後平均＝精度を重みにした加重平均になります。
「精度＝情報量」。データが増えるほど事後は標本平均 $\bar x$ へ収縮し鋭くなる。事前は弱い正則化として効き、これがリッジ回帰（ガウス事前の MAP）の素地です。
分散も未知なら逆ガンマ、平均・分散とも未知なら正規逆ガンマ（NIG）が共役。このとき事後予測は Student-t になり、分散の不確実性ぶん正規より裾が重くなります。

1. 分散既知：精度の言葉で見る

データの分散 $\sigma^2$ が既知で、平均 $\mu$ だけを推定します。式が一番すっきりするのは精度（precision） $\tau=1/\sigma^2$ で書くときです。精度とは「分散の逆数＝どれだけ情報が詰まって尖っているか」。

データ： $x_1,\dots,x_n\sim\mathcal N(\mu,\,1/\tau)$ 、標本平均 $\bar x$ 、データ精度 $\tau$ （既知）
事前： $\mu\sim\mathcal N(\mu_0,\,1/\tau_0)$

尤度も事前も $\mu$ の二次式 $\exp\{-\tfrac12(\cdots)\mu^2+(\cdots)\mu\}$ の形なので、掛けて指数の肩を平方完成すると再び正規のカーネルになります（正規 × 正規＝正規。完全な平方完成は共役事前分布）。結果は

\mu\mid x\sim\mathcal N(\mu_n,\,1/\tau_n),\qquad \tau_n=\tau_0+n\tau,\qquad \mu_n=\frac{\tau_0\mu_0+n\tau\,\bar x}{\tau_0+n\tau}

二つの式の読み方がこのトピックの核心です。

(i) 事後精度は精度の足し算： $\tau_n=\underbrace{\tau_0}_{\text{事前}}+\underbrace{n\tau}_{\text{データ }n\text{ 個ぶん}}$ 。情報（精度）は加算され、データが増えるほど事後が鋭くなります。

(ii) 事後平均は精度を重みにした加重平均：

\mu_n=\frac{\tau_0}{\tau_0+n\tau}\,\mu_0+\frac{n\tau}{\tau_0+n\tau}\,\bar x

「自信のある（精度の高い）側へ事後平均が寄る」。データが少なければ事前 $\mu_0$ 寄り、多ければ標本平均 $\bar x$ 寄りです。

import numpy as np
from scipy import stats
from scipy.integrate import trapezoid          # numpy 2.0+ で np.trapz は廃止

rng = np.random.default_rng(7)
sigma2 = 4.0                                     # 既知のデータ分散 → 精度 τ=1/σ²
tau = 1/sigma2
data = rng.normal(2.5, np.sqrt(sigma2), size=10)
n = len(data); xbar = data.mean()
mu0, tau0 = 0.0, 1.0                             # 事前 N(μ0, 1/τ0)

tau_n = tau0 + n*tau                             # 事後精度＝精度の足し算
mu_n  = (tau0*mu0 + n*tau*xbar) / tau_n          # 事後平均＝精度の加重平均
print(f"n={n}, x̄={xbar:.4f}")
print(f"事後 N(μ_n={mu_n:.4f}, 分散={1/tau_n:.4f}), 事後SD={np.sqrt(1/tau_n):.4f}")

# 閉じた式 vs 数値積分（μ のグリッド）で一致を確認
mu = np.linspace(-3, 6, 6000)
loglik = -0.5*tau*np.sum((data[:,None]-mu[None,:])**2, axis=0)   # 正規尤度の対数
like = np.exp(loglik - loglik.max())
pri  = stats.norm(mu0, np.sqrt(1/tau0)).pdf(mu)
post_num = (like*pri)/trapezoid(like*pri, mu)
post_closed = stats.norm(mu_n, np.sqrt(1/tau_n)).pdf(mu)
print(f"閉じた式 vs 数値積分 の最大差 = {np.max(np.abs(post_num-post_closed)):.2e}")

出力：

n=10, x̄=2.0953
事後 N(μ_n=1.4967, 分散=0.2857), 事後SD=0.5345
閉じた式 vs 数値積分 の最大差 = 5.55e-16

出力の意味：標本平均は $\bar x=2.10$ ですが、事後平均は $1.50$ ——事前 $\mu_0=0$ に引っぱられて、 $\bar x$ と $\mu_0$ の間に落ちています。事前精度 $\tau_0=1$ に対しデータ精度は $n\tau=10\times0.25=2.5$ なので、データのほうが少し強く、事後平均は $\bar x$ 寄りです。閉じた式とグリッド数値積分の差は $10^{-16}$ で完全一致。

2. 収縮：データが増えると事後は標本平均へ

加重平均の重み $\dfrac{n\tau}{\tau_0+n\tau}$ は $n$ とともに 1 に近づきます。同じ $\bar x=2.5$ ・同じ事前 $\mathcal N(0,1)$ で、 $n$ を増やしたときの事後平均を見ます。

tau0, tau, mu0, xbar_fix = 1.0, 0.25, 0.0, 2.5
print(f"{'n':>5}{'事前重み':>10}{'データ重み':>10}{'事後平均':>10}")
for nn in [1, 5, 20, 100]:
    w_prior = tau0/(tau0 + nn*tau)
    w_data  = nn*tau/(tau0 + nn*tau)
    print(f"{nn:>5}{w_prior:>10.3f}{w_data:>10.3f}{w_prior*mu0 + w_data*xbar_fix:>10.3f}")

出力：

    n      事前重み     データ重み      事後平均
    1     0.800     0.200     0.500
    5     0.444     0.556     1.389
   20     0.167     0.833     2.083
  100     0.038     0.962     2.404

出力の意味： $n=1$ では事前重み $0.8$ で事後平均は $0.5$ （事前 $0$ にべったり）、 $n=100$ では事前重みが $0.04$ まで落ち、事後平均は $2.40$ とほぼ標本平均 $2.5$ 。 $n\to\infty$ で $\mu_n\to\bar x$ 、事前の影響は消えます。この「事前へ引き戻す」性質を**収縮（shrinkage）**と呼び、小標本の推定を安定させます（階層モデルでの収縮は収縮の数理で深掘り）。

正則化との接続：MAP（＝この正規では事後平均）は次を最小化する $\mu$ です。
$\arg\min_\mu\ \frac{\tau}{2}\sum_i (x_i-\mu)^2\ +\ \frac{\tau_0}{2}(\mu-\mu_0)^2$
第2項は「 $\mu$ を $\mu_0$ から離すと罰する」L2 ペナルティそのもの。ガウス事前の MAP ＝リッジ正則化で、事前精度 $\tau_0$ が正則化の強さに対応します（回帰での完全版は正則化との関係、機械学習側は正則化の理論・正則化（Ridge・Lasso・Elastic Net））。

3. 分散未知（平均既知）：逆ガンマが共役

平均 $\mu$ が既知で分散 $\sigma^2$ を推定する場合、共役事前は逆ガンマ分布 $\mathrm{Inv\text{-}Gamma}(a_0,b_0)$ です（精度 $\tau=1/\sigma^2$ で見ればガンマ）。観測の平方和 $SS=\sum_i(x_i-\mu)^2$ を使って、事後は

\sigma^2\mid x\sim\mathrm{Inv\text{-}Gamma}\!\Big(a_0+\tfrac{n}{2},\ b_0+\tfrac{1}{2}SS\Big)

import numpy as np
from scipy import stats

rng = np.random.default_rng(20)
mu_known = 2.5
data = rng.normal(mu_known, 2.0, size=12)
n = len(data); SS = np.sum((data - mu_known)**2)
a0, b0 = 3.0, 6.0                                # 事前 InvGamma(a0, b0)
a_n, b_n = a0 + n/2, b0 + SS/2                    # 事後 InvGamma
print(f"SS={SS:.3f}, n={n} → 事後 InvGamma(a={a_n:.1f}, b={b_n:.2f})")
print(f"σ² の事後平均 = b/(a-1) = {b_n/(a_n-1):.4f}")

出力：

SS=30.270, n=12 → 事後 InvGamma(a=9.0, b=21.13)
σ² の事後平均 = b/(a-1) = 2.6419

出力の意味：データは真の分散 $4$ （標準偏差 $2$ ）から生成しましたが、 $n=12$ の小標本なので、事後の $\sigma^2$ 平均は $2.64$ とやや小さめに出ています。 $a_0$ は「擬似的な観測数の半分」、 $b_0$ は「擬似的な平方和の半分」に相当し、ここでも擬似観測として読めます。

4. 両方未知：事後予測は Student-t

平均も分散も未知のとき、共役事前は正規逆ガンマ（NIG） $\;\mu\mid\sigma^2\sim\mathcal N(\mu_0,\sigma^2/\kappa_0),\ \sigma^2\sim\mathrm{Inv\text{-}Gamma}(\alpha_0,\beta_0)$ 。事後も NIG で、ハイパーパラメータは

\kappa_n=\kappa_0+n,\quad \mu_n=\frac{\kappa_0\mu_0+n\bar x}{\kappa_n},\quad \alpha_n=\alpha_0+\frac{n}{2},\quad \beta_n=\beta_0+\tfrac12\textstyle\sum(x_i-\bar x)^2+\frac{\kappa_0 n(\bar x-\mu_0)^2}{2\kappa_n}

と更新されます。重要なのは事後予測分布——次の観測 $x'$ について $\mu,\sigma^2$ を両方とも積分すると、正規ではなくStudent-t になります。

x'\mid D\sim t_{2\alpha_n}\!\left(\mu_n,\ \frac{\beta_n(\kappa_n+1)}{\alpha_n\kappa_n}\right)

「分散がまだ不確かだから、正規より裾が重い」。これがベイズで予測区間が現実的な広さになる理由です。モンテカルロで確かめます。

import numpy as np
from scipy import stats

rng = np.random.default_rng(42)
data = rng.normal(2.5, 2.0, size=15)
n = len(data); xb = data.mean(); ss = np.sum((data - xb)**2)
mu0, kappa0, alpha0, beta0 = 0.0, 1.0, 2.0, 2.0      # NIG 事前

kappa_n = kappa0 + n
mu_n    = (kappa0*mu0 + n*xb) / kappa_n
alpha_n = alpha0 + n/2
beta_n  = beta0 + 0.5*ss + 0.5*kappa0*n*(xb - mu0)**2 / kappa_n
df    = 2*alpha_n
scale = np.sqrt(beta_n*(kappa_n + 1)/(alpha_n*kappa_n))
print(f"n={n}, x̄={xb:.4f}")
print(f"事後予測 = t(df={df:.0f}, loc={mu_n:.4f}, scale={scale:.4f})")

# モンテカルロ：σ²→μ→x' と階層的にサンプルし、Student-t と比べる
tdist = stats.t(df=df, loc=mu_n, scale=scale)
s2   = stats.invgamma(a=alpha_n, scale=beta_n).rvs(400000, random_state=rng)
mu_s = rng.normal(mu_n, np.sqrt(s2/kappa_n))
xnew = rng.normal(mu_s, np.sqrt(s2))
for q in [0.05, 0.5, 0.95]:
    print(f"  分位{q}: MC={np.quantile(xnew, q):+.2f}  t理論={tdist.ppf(q):+.2f}")
print(f"  予測SD: MC={xnew.std():.3f}  t理論={tdist.std():.3f}  正規近似={np.sqrt(beta_n/(alpha_n-1)):.3f}")

出力：

n=15, x̄=2.4950
事後予測 = t(df=19, loc=2.3390, scale=1.7885)
  分位0.05: MC=-0.74  t理論=-0.75
  分位0.5: MC=+2.34  t理論=+2.34
  分位0.95: MC=+5.44  t理論=+5.43
  予測SD: MC=1.888  t理論=1.891  正規近似=1.834

出力の意味：階層的にサンプルした予測の分位点が、自由度 19 の Student-t とほぼ一致します（残差はモンテカルロ誤差）。予測 SD は $1.89$ で、分散を点に固定した「正規近似」 $1.83$ より広い。差は小さく見えますが、 $t$ は裾が重く、95%予測区間の幅や外れ値の確率で効いてきます。自由度 $2\alpha_n$ はデータが増えると大きくなり、 $t$ は正規に近づきます——「分散の不確実性が減れば予測も正規に戻る」わけです。

裾の重さを重ねて描くと違いが見えます。

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

# §4 で求めた事後予測 t(df=19, loc=2.339, scale=1.789) と正規近似を重ねる
df, mu_n, scale, sd_norm = 19, 2.339, 1.789, 1.834
x = np.linspace(mu_n - 8, mu_n + 8, 400)
plt.figure(figsize=(7, 4))
plt.plot(x, stats.t(df=df, loc=mu_n, scale=scale).pdf(x), lw=2,
         label=f"事後予測 t(df={df})")
plt.plot(x, stats.norm(mu_n, sd_norm).pdf(x), ls="--", lw=2,
         label="正規近似（分散を点に固定）")
plt.xlabel("x'（次の観測）"); plt.ylabel("予測密度")
plt.title("事後予測 Student-t は正規より裾が重い")
plt.legend(); plt.tight_layout(); plt.show()

グラフの意味：中心付近はほぼ重なりますが、両裾で Student-t（実線）が正規（破線）より上に出ます。この「裾の厚み」が、分散をまだ確定できていない不確実性の正直な表れです。点に固定した正規で予測区間を作ると、この裾を取りこぼして外れ値の確率を過小評価します。

⚠️ よくある誤解

「正規の共役はいつも正規」ではない：未知母数で変わります。平均未知・分散既知→正規、分散未知・平均既知→逆ガンマ、両方未知→正規逆ガンマ。分散を推定する問題で「正規だ」と短絡すると誤りです。
「精度と分散の混同」：更新が足し算になるのは精度 $\tau=1/\sigma^2$ で書いたとき。分散のまま「足す」と間違えます。
「事後予測は正規」ではない：分散が未知なら予測は Student-t。正規で予測区間を作ると、分散の不確実性ぶん過信（狭すぎ）になります。
「無情報のつもりで $\tau_0$ を大きく」：事前精度 $\tau_0$ を大きくするのは「 $\mu_0$ を強く信じる」こと。無情報にしたいなら $\tau_0$ を小さく（事前分散を大きく）します。