最小二乗法などで回帰モデルの重みパラメータを決定することは統計学の基本ですが、今回はベイズ推論（ベイジアンモデリング）の枠組みで線形回帰を考えていきます。
この線形回帰問題をベイズの枠組みで行うことをベイズ線形回帰と呼びます。ベイズ線形回帰では、パラメータの点推定ではなく、パラメータの事後分布を求め、予測値の不確かさも含めた予測分布を得ることができます。
本記事の内容
- ベイズ線形回帰の定式化
- パラメータの事後分布の導出
- 予測分布の導出
- Python での実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
モデルの定式化
線形回帰モデル
入力 $\bm{x}$ に対する出力 $y$ を次のようにモデル化します。
$$ \begin{equation} y = \bm{w}^T \bm{\phi}(\bm{x}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \beta^{-1}) \end{equation} $$
ここで $\bm{w}$ はパラメータベクトル、$\bm{\phi}(\bm{x})$ は基底関数ベクトル、$\beta$ は観測ノイズの精度（分散の逆数）です。
尤度関数
$N$ 個のデータ $\{(\bm{x}_n, y_n)\}_{n=1}^{N}$ に対する尤度関数は、
$$ \begin{equation} p(\bm{y}|\bm{X}, \bm{w}) = \prod_{n=1}^{N} \mathcal{N}(y_n|\bm{w}^T\bm{\phi}(\bm{x}_n), \beta^{-1}) = \mathcal{N}(\bm{y}|\bm{\Phi}\bm{w}, \beta^{-1}\bm{I}) \end{equation} $$
ここで $\bm{\Phi}$ は計画行列で、第 $n$ 行が $\bm{\phi}(\bm{x}_n)^T$ です。
事前分布
パラメータ $\bm{w}$ の事前分布として多変量ガウス分布を設定します。
$$ \begin{equation} p(\bm{w}) = \mathcal{N}(\bm{w}|\bm{m}_0, \bm{S}_0) \end{equation} $$
簡単のため $\bm{m}_0 = \bm{0}$, $\bm{S}_0 = \alpha^{-1}\bm{I}$ とする場合が多いです（$\alpha$ は正則化パラメータ）。
パラメータの事後分布
ガウス分布の共役性により、事後分布もガウス分布になります。
$$ \begin{equation} p(\bm{w}|\bm{y}, \bm{X}) = \mathcal{N}(\bm{w}|\bm{m}_N, \bm{S}_N) \end{equation} $$
ここで、
$$ \begin{align} \bm{S}_N^{-1} &= \bm{S}_0^{-1} + \beta \bm{\Phi}^T \bm{\Phi} \\ \bm{m}_N &= \bm{S}_N (\bm{S}_0^{-1}\bm{m}_0 + \beta \bm{\Phi}^T \bm{y}) \end{align} $$
$\bm{m}_0 = \bm{0}$, $\bm{S}_0 = \alpha^{-1}\bm{I}$ の場合は、
$$ \begin{align} \bm{S}_N^{-1} &= \alpha\bm{I} + \beta \bm{\Phi}^T \bm{\Phi} \\ \bm{m}_N &= \beta \bm{S}_N \bm{\Phi}^T \bm{y} \end{align} $$
MAP推定と正則化
事後分布の最大値（MAP推定）は、
$$ \hat{\bm{w}}_{\text{MAP}} = \bm{m}_N = \beta \bm{S}_N \bm{\Phi}^T \bm{y} $$
これは L2 正則化（Ridge回帰）の解 $(\alpha\bm{I} + \beta\bm{\Phi}^T\bm{\Phi})^{-1}\beta\bm{\Phi}^T\bm{y}$ に一致します。
予測分布
新しい入力 $\bm{x}_*$ に対する予測分布は、
$$ \begin{align} p(y_*|\bm{x}_*, \bm{y}, \bm{X}) &= \int p(y_*|\bm{w}, \bm{x}_*) p(\bm{w}|\bm{y}, \bm{X}) \, d\bm{w} \\ &= \mathcal{N}(y_*|\bm{m}_N^T\bm{\phi}(\bm{x}_*), \sigma_*^2) \end{align} $$
ここで予測分散は、
$$ \begin{equation} \sigma_*^2 = \beta^{-1} + \bm{\phi}(\bm{x}_*)^T \bm{S}_N \bm{\phi}(\bm{x}_*) \end{equation} $$
第1項 $\beta^{-1}$ は観測ノイズに起因する不確かさ、第2項 $\bm{\phi}(\bm{x}_*)^T \bm{S}_N \bm{\phi}(\bm{x}_*)$ はパラメータの不確かさに起因する不確かさです。データが増えると第2項が減少し、予測の不確かさが小さくなります。
Python での実装
import numpy as np
import matplotlib.pyplot as plt
# 基底関数（多項式基底）
def polynomial_basis(x, degree):
return np.column_stack([x**i for i in range(degree + 1)])
# ベイズ線形回帰
class BayesianLinearRegression:
def __init__(self, alpha, beta, degree):
self.alpha = alpha # 事前分布の精度
self.beta = beta # 観測ノイズの精度
self.degree = degree
def fit(self, X, y):
Phi = polynomial_basis(X, self.degree)
M = Phi.shape[1]
# 事前分布
self.S0_inv = self.alpha * np.eye(M)
self.m0 = np.zeros(M)
# 事後分布
self.SN_inv = self.S0_inv + self.beta * Phi.T @ Phi
self.SN = np.linalg.inv(self.SN_inv)
self.mN = self.SN @ (self.S0_inv @ self.m0 + self.beta * Phi.T @ y)
def predict(self, X_new):
Phi_new = polynomial_basis(X_new, self.degree)
y_mean = Phi_new @ self.mN
y_var = 1.0 / self.beta + np.sum(Phi_new @ self.SN * Phi_new, axis=1)
return y_mean, y_var
# データ生成
np.random.seed(42)
true_w = np.array([1.0, -0.5, 0.3]) # y = 1 - 0.5x + 0.3x^2
N = 20
beta_true = 25 # 観測ノイズの精度
x_train = np.random.uniform(-3, 3, N)
y_train = polynomial_basis(x_train, 2) @ true_w + np.random.normal(0, 1/np.sqrt(beta_true), N)
# モデルの学習
model = BayesianLinearRegression(alpha=0.1, beta=beta_true, degree=2)
model.fit(x_train, y_train)
# 予測
x_test = np.linspace(-4, 4, 200)
y_mean, y_var = model.predict(x_test)
y_std = np.sqrt(y_var)
# 真の関数
y_true = polynomial_basis(x_test, 2) @ true_w
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左: 予測分布
ax1 = axes[0]
ax1.scatter(x_train, y_train, c='blue', s=40, zorder=5, label='Training data')
ax1.plot(x_test, y_true, 'g--', linewidth=1.5, label='True function')
ax1.plot(x_test, y_mean, 'r-', linewidth=2, label='Predictive mean')
ax1.fill_between(x_test, y_mean - 2*y_std, y_mean + 2*y_std,
alpha=0.2, color='red', label='95% CI')
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Bayesian Linear Regression: Prediction', fontsize=13)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-5, 8)
# 右: データ数による不確かさの変化
ax2 = axes[1]
n_values = [3, 5, 10, 20]
colors = ['purple', 'blue', 'green', 'red']
for n, color in zip(n_values, colors):
model_n = BayesianLinearRegression(alpha=0.1, beta=beta_true, degree=2)
model_n.fit(x_train[:n], y_train[:n])
y_m, y_v = model_n.predict(x_test)
y_s = np.sqrt(y_v)
ax2.fill_between(x_test, y_m - 2*y_s, y_m + 2*y_s, alpha=0.15, color=color)
ax2.plot(x_test, y_m, color=color, linewidth=1.5, label=f'N={n}')
ax2.plot(x_test, y_true, 'k--', linewidth=1.5, label='True')
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('y', fontsize=12)
ax2.set_title('Uncertainty Decreases with More Data', fontsize=13)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-5, 8)
plt.tight_layout()
plt.show()
左のグラフでは、予測平均が真の関数に近く、データが疎な領域では予測の不確かさ（信頼区間の幅）が大きくなっていることが確認できます。右のグラフでは、データ数が増えるにつれて不確かさが減少する様子が視覚的に確認できます。
通常の線形回帰との比較
| 観点 | 通常の線形回帰 | ベイズ線形回帰 |
|---|---|---|
| パラメータ推定 | 点推定 $\hat{\bm{w}}$ | 分布 $p(\bm{w} \mid \bm{y})$ |
| 予測 | 点予測 $\hat{y}$ | 予測分布 $p(y_* \mid \bm{x}_*)$ |
| 不確かさ | 信頼区間（漸近的） | 自然に得られる |
| 正則化 | 別途追加（Ridge, Lasso） | 事前分布として自然に組み込み |
| 過学習 | 起こりやすい | 事前分布が緩和 |
まとめ
本記事では、ベイズ線形回帰について解説しました。
- ベイズ線形回帰はパラメータの事後分布を解析的に求め、予測の不確かさを定量化できる
- 事後分布はガウス分布の共役性により閉じた形で得られる
- 予測分散は観測ノイズとパラメータの不確かさの和で構成される
- MAP推定はL2正則化（Ridge回帰）に対応する
次のステップとして、以下の記事も参考にしてください。