ベイズ線形回帰｜ベイズ統計

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

📎 前提：VAEとELBOの共通構造　|　数理：正規正規モデル・正則化（リッジ・Lasso）（統計）

要点（BLUF）

ベイズ線形回帰は、回帰係数 $w$ にガウス事前を置き、ガウス尤度と合わせると**事後も予測分布も閉形式（ガウス）**になります。多変量の正規–正規共役（正規正規モデル）の応用です。
事後予測分布は重み $w$ の不確実性を積分したもの。予測平均は $\phi(x_*)^\top m_N$ 、予測分散は雑音由来＋重み由来の和。
だから予測には不確実性の帯がつき、データから離れる（外挿）ほど帯が広がる——点推定の回帰にはない正直さです。

1. モデルと事後（閉形式）

入力 $x$ を特徴ベクトル $\phi(x)$ （ここでは $[1, x]$ ）に写し、 $y=\phi(x)^\top w+\varepsilon,\ \varepsilon\sim\mathcal N(0,\sigma^2)$ 。係数に事前 $w\sim\mathcal N(0,\alpha^{-1}I)$ を置きます。尤度も事前もガウスなので、事後は平方完成で再びガウス（正規正規モデルの多変量版）：

w\mid D\sim\mathcal N(m_N,S_N),\qquad S_N=(\alpha I+\beta\Phi^\top\Phi)^{-1},\qquad m_N=\beta\,S_N\Phi^\top y

ここで $\beta=1/\sigma^2$ （雑音精度）、 $\Phi$ は計画行列（各行が $\phi(x_i)$ ）。 $\alpha$ が事前精度で、大きいほど係数を $0$ へ強く引きます（正則化との関係のリッジ）。

2. コード：事後と事後予測

import numpy as np

rng = np.random.default_rng(1)
true_w = np.array([0.5, -1.2]); sigma = 0.4; beta = 1/sigma**2
alpha = 1.0                                   # 事前精度
x = np.sort(rng.uniform(-1, 1, 25))
X = np.c_[np.ones_like(x), x]                 # 計画行列 [1, x]
y = X @ true_w + rng.normal(0, sigma, len(x))

# 事後 w|D ~ N(m_N, S_N)（閉形式）
S_N = np.linalg.inv(alpha*np.eye(2) + beta*X.T@X)
m_N = beta * S_N @ X.T @ y
print(f"事後平均 w=[{m_N[0]:.3f}, {m_N[1]:.3f}] (真{true_w})")
print(f"事後SD   w=[{np.sqrt(S_N[0,0]):.3f}, {np.sqrt(S_N[1,1]):.3f}]")

# 事後予測 p(y*|x*)：平均=φ(x*)ᵀm_N, 分散=1/β + φ(x*)ᵀS_N φ(x*)
def pred(xs):
    Phi = np.c_[np.ones_like(xs), xs]
    mean = Phi @ m_N
    var = 1/beta + np.einsum('ij,jk,ik->i', Phi, S_N, Phi)
    return mean, np.sqrt(var)

for xs in [0.0, 1.0, 3.0]:                    # データ中心 / 端 / 外挿
    mean, sd = pred(np.array([xs]))
    print(f"x*={xs:>4}: 予測平均={mean[0]:+.3f}  予測SD={sd[0]:.3f}")

出力：

事後平均 w=[0.465, -1.113] (真[ 0.5 -1.2])
事後SD   w=[0.080, 0.141]
x*= 0.0: 予測平均=+0.465  予測SD=0.408
x*= 1.0: 予測平均=-0.648  予測SD=0.431
x*= 3.0: 予測平均=-2.874  予測SD=0.586

出力の意味：事後平均 $w=[0.47,-1.11]$ は真の $[0.5,-1.2]$ をよく回収。予測 SD はデータの中心 $x_*=0$ で $0.41$ 、端 $x_*=1$ で $0.43$ 、外挿 $x_*=3$ で $0.59$ とデータから離れるほど増えます。予測分散の式 $1/\beta+\phi(x_*)^\top S_N\phi(x_*)$ の第1項が「データ自体の雑音」、第2項が「係数の不確実性が予測に伝わるぶん」で、後者が外挿で膨らみます。点推定回帰なら一定幅の予測になるところを、ベイズは「知らない領域では自信がない」と帯を広げます。

3. 図：予測の不確実性帯

事後予測の平均と $\pm2$ SD の帯を、事後からサンプルした回帰直線とともに描きます。

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

rng = np.random.default_rng(1)
true_w = np.array([0.5, -1.2]); sigma = 0.4; beta = 1/sigma**2; alpha = 1.0
x = np.sort(rng.uniform(-1, 1, 25)); X = np.c_[np.ones_like(x), x]
y = X @ true_w + rng.normal(0, sigma, len(x))
S_N = np.linalg.inv(alpha*np.eye(2) + beta*X.T@X); m_N = beta*S_N@X.T@y

xs = np.linspace(-2.5, 2.5, 200); Phi = np.c_[np.ones_like(xs), xs]
mean = Phi@m_N
sd = np.sqrt(1/beta + np.einsum('ij,jk,ik->i', Phi, S_N, Phi))

plt.figure(figsize=(7.5, 4.5))
for w in rng.multivariate_normal(m_N, S_N, 30):   # 事後からの回帰直線サンプル
    plt.plot(xs, Phi@w, color="C1", alpha=0.15, lw=0.8)
plt.fill_between(xs, mean-2*sd, mean+2*sd, alpha=0.2, color="C0", label="予測 ±2SD")
plt.plot(xs, mean, color="C0", lw=2, label="予測平均")
plt.scatter(x, y, color="black", s=20, zorder=5, label="データ")
plt.axvspan(-1, 1, color="gray", alpha=0.08)      # データのある範囲
plt.xlabel("x"); plt.ylabel("y"); plt.legend(fontsize=9)
plt.title("ベイズ線形回帰：外挿で予測帯が広がる")
plt.tight_layout(); plt.show()

グラフの意味：データのある $[-1,1]$ （灰色帯）の中では予測帯が細く、事後からサンプルした直線（橙）もよく揃います。外側へ出ると直線が扇状に開き、予測帯（青）が広がる——「データのない所では傾きの不確実性が予測のばらつきになる」様子がそのまま見えます。

4. 正規–正規の応用、そして正則化へ

ベイズ線形回帰は、正規正規モデルの「分散既知の正規–正規」を多変量・回帰に拡張したものです。事後精度は $S_N^{-1}=\alpha I+\beta\Phi^\top\Phi$ ＝事前精度＋データ精度の足し算、事後平均は精度の加重平均——構図はまったく同じ。事前精度 $\alpha$ が係数を $0$ へ収縮させるので、事後平均（＝MAP）はリッジ回帰の解に一致します。この対応を次節正則化との関係で厳密に確かめます。非ガウス（ロジスティック・ポアソン）の回帰は共役が崩れるので、ラプラス近似や MCMC を使います（ベイズGLM）。

⚠️ よくある誤解

「ベイズ回帰の出力は1本の直線」ではない：出力は係数の事後分布＝直線の分布。予測も分布（帯）です。
「予測帯はどこでも一定幅」ではない：データから離れるほど広がります（§2-3）。一定幅は係数の不確実性を無視した点推定回帰。
「事前は結果に影響しない」ではない： $\alpha$ が係数を $0$ へ収縮させます（リッジ）。 $\alpha\to0$ で最小二乗に近づく。
「閉形式だから MCMC は不要」は線形・ガウスに限る：リンク関数が入る GLM（ベイズGLM）では共役が崩れ、近似や MCMC が要ります。