ベイズGLM｜ベイズ統計

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

📎 前提：正則化との関係　|　数理：指数型分布族と共役事前

要点（BLUF）

**一般化線形モデル（GLM）**は、リンク関数で線形予測子 $x^\top w$ を確率や率に繋ぎます（ロジスティック＝ロジットリンク、ポアソン＝対数リンク）。ガウス共役が崩れ、事後は閉形式になりません。
ラプラス近似：MAP を求め、事後を MAP 周りのガウス（共分散＝負の対数事後のヘッセ行列の逆）で近似します。手早く事後平均・SD が出ます。
ベイズロジスティック回帰をラプラス近似で解き、フル MCMC とほぼ一致することを確かめます。非対称・多峰では近似が崩れるので MCMC を使います。

1. GLM とリンク関数：なぜ共役が崩れるか

線形回帰（ベイズ線形回帰）は出力が連続でガウスでした。GLM は出力が二値・カウントなどのとき、指数型分布族（指数型分布族と共役事前）の尤度にリンク関数を挟みます。

ロジスティック回帰： $y\in\{0,1\}$ 、 $p=\sigma(x^\top w)$ （ $\sigma$ はロジスティック関数）。ロジット $\log\frac{p}{1-p}=x^\top w$ 。
ポアソン回帰： $y\in\{0,1,2,\dots\}$ 、 $\lambda=\exp(x^\top w)$ 。対数リンク $\log\lambda=x^\top w$ 。

これらの尤度にガウス事前 $w\sim\mathcal N(0,\alpha^{-1}I)$ を掛けても、 $\sigma$ や $\exp$ が非線形に絡むため、事後が標準分布になりません（共役でない）。事後を解析的に書けないので近似が要ります。

2. ラプラス近似：MAP 周りをガウスで

最も手軽な近似がラプラス近似。対数事後 $\log p(w\mid D)$ を MAP $w_{\mathrm{MAP}}$ の周りで2次までテイラー展開すると、対数が2次関数＝ガウスのカーネルになります。

p(w\mid D)\approx\mathcal N\big(w_{\mathrm{MAP}},\ H^{-1}\big),\qquad H=-\nabla^2\log p(w\mid D)\big|_{w_{\mathrm{MAP}}}

$H$ は負の対数事後のヘッセ行列（曲率）。ロジスティック回帰では $H=\Phi^\top S\Phi+\alpha I$ （ $S=\mathrm{diag}(p_i(1-p_i))$ ）。MAP はニュートン法（＝IRLS）で速く求まります。事後の山を1つのガウスで近似するので、単峰でほぼ対称な事後によく当たります。

3. コード：ベイズロジスティック回帰（ラプラス vs MCMC）

ラプラス近似が、計算の重いフル MCMC（メトロポリスヘイスティングス）とどれだけ一致するか確かめます。

import numpy as np

rng = np.random.default_rng(0)
n = 200
true_w = np.array([0.4, 1.5, -1.0])                      # 切片, x1, x2
X = np.c_[np.ones(n), rng.standard_normal((n, 2))]
y = (rng.uniform(size=n) < 1/(1+np.exp(-X@true_w))).astype(float)
alpha = 1.0                                              # 事前 N(0, 1/α)

sig = lambda z: 1/(1+np.exp(-z))
def neg_logpost(w):
    pz = sig(X@w)
    return -(y@np.log(pz+1e-12) + (1-y)@np.log(1-pz+1e-12)) + 0.5*alpha*w@w

# MAP（ニュートン法）→ ラプラス近似 事後≈N(w_map, H⁻¹)
w = np.zeros(3)
for _ in range(50):
    pz = sig(X@w)
    grad = X.T@(pz - y) + alpha*w
    H = X.T@(X*(pz*(1-pz))[:,None]) + alpha*np.eye(3)
    w = w - np.linalg.solve(H, grad)
w_map = w
pz = sig(X@w_map)
H = X.T@(X*(pz*(1-pz))[:,None]) + alpha*np.eye(3)
sd_laplace = np.sqrt(np.diag(np.linalg.inv(H)))
print(f"ラプラス: 事後平均={np.round(w_map,3)}  事後SD={np.round(sd_laplace,3)}")

# MH で検証
def mh(N=120000, burn=20000, step=0.15, seed=1):
    rng = np.random.default_rng(seed); w = w_map.copy(); lp = -neg_logpost(w)
    out = np.empty((N,3))
    for i in range(N):
        cand = w + rng.normal(0, step, 3); lc = -neg_logpost(cand)
        if np.log(rng.uniform()) < lc-lp: w, lp = cand, lc
        out[i] = w
    return out[burn:]
s = mh()
print(f"MH      : 事後平均={np.round(s.mean(0),3)}  事後SD={np.round(s.std(0),3)}  (真{true_w})")

出力：

ラプラス: 事後平均=[ 0.381  1.652 -1.104]  事後SD=[0.18  0.246 0.219]
MH      : 事後平均=[ 0.385  1.688 -1.13 ]  事後SD=[0.182 0.251 0.223]  (真[ 0.4  1.5 -1. ])

出力の意味：ラプラス近似の事後平均 $[0.38,1.65,-1.10]$ ・SD $[0.18,0.25,0.22]$ が、フル MCMC の結果とほぼ一致します（差は小数第2位）。どちらも真の係数 $[0.4,1.5,-1.0]$ を回収。ロジスティック回帰の事後はこの規模のデータでほぼ単峰・準対称なので、MAP とヘッセ行列だけで MCMC 並みの事後が得られる——ニュートン法数十回 vs MH 十数万回で、ラプラスが圧倒的に軽い。係数の SD（不確実性）まで出るのが点推定の最尤ロジスティックとの違いです。

4. ポアソン回帰も同じ枠組み

カウントデータのポアソン回帰 $\lambda=\exp(x^\top w)$ も同様です。対数事後の MAP をニュートン法で求め、ヘッセ行列でラプラス近似——リンク関数と尤度が変わるだけで、手順は共通。GLM はすべて指数型分布族（指数型分布族と共役事前）の尤度＋リンク関数で、自然パラメータ $\eta=x^\top w$ を線形予測子に結ぶ構造を共有します。だから「指数型分布族の自然パラメータを線形にする」のが GLM、と一言で言えます。

精度がもっと要る・事後が非対称や多峰なら、ラプラスでなく MCMC（第4章）や変分推論（第6章）を使います。実務では PyMC などで pm.Bernoulli／pm.Poisson を宣言して NUTS で解くのが定番（要最新確認、確率的プログラミング概観）。

⚠️ よくある誤解

「GLM もガウス共役で閉形式」ではない：リンク関数の非線形で共役が崩れます。近似（ラプラス・VI）か MCMC が必要。
「ラプラス近似はいつでも正確」ではない：単峰・準対称な事後に有効。多峰や強い歪み・少数データ（分離）では外します。そのとき MCMC へ。
「MAP だけで十分」ではない：ラプラスは MAP に加え**ヘッセ行列で不確実性（SD）**も出します。最尤ロジスティックの点推定との差はここ。
「ロジスティック回帰のベイズ化は重い」ではない：ラプラス近似ならニュートン法数十回で済み、MCMC より遥かに軽い（§3）。