ベイズ線形回帰の基礎:最小二乗法からベイズ推定まで

最小二乗法・最尤推定・MAP推定・ベイズ推定を比較解説。事後分布・予測分布をフルの行列計算で導出し、リッジ回帰との等価性を証明。sklearnのBayesianRidge/Ridgeとの数値一致、N=8とN=100での予測不確実性の収縮をPythonで実測検証します。

はじめに

回帰問題は、入力データとそれに対応する出力データの関係を学習し、未知の入力に対する出力を予測するタスクです。本記事では、最も基本的な回帰モデルである線形回帰を題材に、パラメータの推定方法として代表的な最小二乗法最尤推定MAP推定、そしてベイズ推定を比較し、その違いと特徴を解説します。

前半(1〜3節)で点推定の系譜(最小二乗法→最尤推定→MAP推定)を辿り、後半(4〜6節)でベイズ推定の事後分布・予測分布を行列計算を省略せずに導出します。特に4節では、「MAP推定=リッジ回帰」という主張を天下り的に受け入れるのではなく、係数の対応関係を明示的に証明します。最後に、すべての数式をPythonで実装し、scikit-learn との数値一致・過学習実験の実測値で検証します。

モデルとして、入力 \(x\) の非線形関数である基底関数 \(\phi(x)\) の線形結合を考えます。\(w\) をモデルのパラメータ(重み)、\(\epsilon\) を誤差とすると、モデルは以下のように表せます。

\[ y = \Phi w + \epsilon \]

ここで、\(\Phi\) は計画行列と呼ばれ、各データ点の基底関数ベクトルを並べたものです(\(N\) をデータ数、\(M\) を基底関数の数とすると \(\Phi \in \mathbb{R}^{N\times M}\) )。


1. 最小二乗法 (Least Squares Estimation)

最小二乗法は、モデルの予測値と実際の目標値との誤差の二乗和 \(S(w)\) を最小化するパラメータ \(\hat{w}\) を見つける手法です。

\[ S(w) = (y - \Phi w)^T (y - \Phi w) \]

\(S(w)\) を \(w\) で微分してゼロとおくことで、以下の解が得られます。

\[ \hat{w}\_{LS} = (\Phi^T \Phi)^{-1} \Phi^T y \]

これは解析的に解けるため、計算が非常に高速です。しかし、訓練データ数が少ない場合や、モデルの自由度が高い場合に**過学習(オーバーフィッティング)**を起こしやすいという欠点があります。また、\(\Phi^T\Phi\) が正則でない(特異、あるいは特異に近い)場合には逆行列が数値的に不安定になります。

2. 最尤推定 (Maximum Likelihood Estimation)

最尤推定は、観測データが得られる確率(尤度)を最大化するパラメータ \(\hat{w}\) を見つける手法です。

ここで、誤差項 \(\epsilon\) が平均 \(0\) 、分散 \(\sigma^2\) のガウス分布に従うと仮定します。すると、目標値 \(y\) の条件付き確率は、平均 \(\Phi w\) 、分散 \(\sigma^2\) のガウス分布となります。

\[ p(y | w, \sigma^2) = \mathcal{N}(y | \Phi w, \sigma^2 I) = \frac{1}{(2\pi\sigma^2)^{N/2}} \exp\lbrace-\frac{1}{2\sigma^2}(y-\Phi w)^T(y-\Phi w)\rbrace \]

この尤度の対数をとった対数尤度を最大化します。

\[ \ln p(y|w) = -\frac{N}{2}\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}(y-\Phi w)^T(y-\Phi w) \]

対数尤度を最大化することは、右辺第二項の二乗誤差項を最小化することと等価です。したがって、誤差にガウス分布を仮定した場合、最尤推定は最小二乗法と等価な解を与えます。

\[ \hat{w}\_{ML} = (\Phi^T\Phi)^{-1}\Phi^{T}y \]

以降では尤度の精度パラメータ(分散の逆数)を \(\beta = 1/\sigma^2\) と表記します。

3. MAP推定 (Maximum A Posteriori Estimation)

MAP推定は、過学習を抑制するための一般的な枠組みです。パラメータ \(w\) 自身も確率変数であると考え、その事前分布 \(p(w)\) を導入します。ベイズの定理を用いて、データが観測された後での \(w\) の事後分布 \(p(w|y)\) を考え、この事後確率が最大となる \(\hat{w}\) を求めます。

\[ p(w|y) = \frac{p(y|w)p(w)}{p(y)} \propto p(y|w)p(w) \]

\(w\) の事前分布として、平均 \(0\) 、共分散 \(\alpha^{-1}I\) のガウス分布を仮定します。これは「\(w\) の各要素は0に近い値をとるだろう」という事前知識をモデルに与えることに相当し、大きな値をとる重みにペナルティを課す正則化として機能します。

\[ p(w|\alpha) = \mathcal{N}(w|0, \alpha^{-1}I) \]

このとき、事後分布を最大化する \(\hat{w}_{MAP}\) は、以下のようになります(導出は4.2節で一般形から厳密に行います)。

\[ \hat{w}\_{MAP} = \left(\Phi^T\Phi + \frac{\alpha}{\beta}I\right)^{-1} \Phi^{T}y \]

この解は、リッジ回帰(L2正則化付き最小二乗法) と呼ばれる手法と同じ形をしており、正則化項 \(\frac{\alpha}{\beta}I\) のおかげで、\(\Phi^T\Phi\) が正則でない(逆行列を持たない)場合でも安定して解を求めることができ、過学習を抑制します。

注意: 比 \(\alpha/\beta\) は「事前分布の精度 ÷ 尤度の精度」であり、\(\beta/\alpha\) ではありません。分子分母を取り違えると正則化強度が大きく変わってしまう(\(\alpha=1,\beta=10\) なら \(0.1\) 倍と \(10\) 倍で100倍違う)ので注意してください。この対応関係は4.2節で明示的に証明します。

4. ベイズ推定 (Bayesian Estimation)

最小二乗法、最尤推定、MAP推定は、いずれも最適なパラメータ \(w\) を一つの値(点推定)として求めました。しかし、このアプローチでは「パラメータがどのくらい確からしいか」という不確実性を表現できません。

ベイズ推定では、\(w\) を点として求めるのではなく、事後分布 \(p(w|y)\) そのものを求めます。この分布は、データを見た後での、あり得る全ての \(w\) の値に対する確率分布を表します。ガウス事前分布とガウス尤度の組み合わせは共役であり、事後分布も厳密にガウス分布になります。以下、これを実際に行列計算で導出します。

4.1 事後分布の導出(完全な行列計算)

MAP推定では事前分布を \(\mathcal{N}(0, \alpha^{-1}I)\) という特殊形に限定しましたが、ここでは一般の平均 \(m_0\) ・共分散 \(S_0\) を持つガウス事前分布

\[ p(w) = \mathcal{N}(w \mid m_0, S_0) \]

から出発します(\(m_0=0, S_0=\alpha^{-1}I\) と置けば3節の設定に一致します)。尤度は

\[ p(y \mid w) = \mathcal{N}(y \mid \Phi w, \beta^{-1}I_N) \]

です。ベイズの定理より事後分布は \(p(w\mid y) \propto p(y\mid w)\, p(w)\) に比例するので、両者の対数(指数部)を足し合わせた

\[ E(w) = \frac{\beta}{2}(y-\Phi w)^T(y-\Phi w) + \frac{1}{2}(w-m_0)^TS_0^{-1}(w-m_0) \]

を \(w\) について整理します。まず各項を展開します。

\[ \begin{aligned} \frac{\beta}{2}(y-\Phi w)^T(y-\Phi w) &= \frac{\beta}{2}\left(y^Ty - 2w^T\Phi^Ty + w^T\Phi^T\Phi w\right) \\ \frac{1}{2}(w-m_0)^TS_0^{-1}(w-m_0) &= \frac{1}{2}\left(w^TS_0^{-1}w - 2w^TS_0^{-1}m_0 + m_0^TS_0^{-1}m_0\right) \end{aligned} \]

\(w\) について2次の項・1次の項・定数項をそれぞれまとめると、

\[ E(w) = \frac{1}{2}w^T\underbrace{(\beta\Phi^T\Phi + S_0^{-1})}_{\text{2次の係数}}w \;-\; w^T\underbrace{(\beta\Phi^Ty + S_0^{-1}m_0)}_{\text{1次の係数}} \;+\; \text{const} \]

となります。一方、事後分布がガウス分布 \(\mathcal{N}(w\mid m_N, S_N)\) であると仮定すると、その指数部は

\[ \frac{1}{2}(w-m_N)^TS_N^{-1}(w-m_N) = \frac{1}{2}w^TS_N^{-1}w - w^TS_N^{-1}m_N + \text{const} \]

という形になります(これを平方完成と呼びます)。この2つの式を係数比較する(2次の係数どうし、1次の係数どうしを等置する)ことで、

\[ \boxed{S_N^{-1} = S_0^{-1} + \beta\Phi^T\Phi} \] \[ S_N^{-1}m_N = S_0^{-1}m_0 + \beta\Phi^Ty \;\;\Longrightarrow\;\; \boxed{m_N = S_N\left(S_0^{-1}m_0 + \beta\Phi^Ty\right)} \]

が得られます。これが事後分布 \(p(w\mid y) = \mathcal{N}(w\mid m_N, S_N)\) の閉形式です。この式は、事前の「確信」(\(S_0^{-1}m_0\) )とデータからの「証拠」(\(\beta\Phi^Ty\) )を精度(分散の逆数)で重み付けして足し合わせたものが事後の精度重み付き平均になっている、と読むことができます。

4.2 正則化との等価性(リッジ回帰との対応の証明)

3節で天下り的に与えた MAP 推定の解を、4.1節の一般公式から厳密に導きます。\(m_0=0\) 、\(S_0=\alpha^{-1}I\) (ゼロ平均・等方性の事前分布)を代入すると、

\[ S_N^{-1} = \alpha I + \beta\Phi^T\Phi, \qquad m_N = S_N\left(\beta\Phi^Ty\right) = \beta\left(\alpha I + \beta\Phi^T\Phi\right)^{-1}\Phi^Ty \]

ここで \(\alpha I + \beta \Phi^T\Phi = \beta\left(\Phi^T\Phi + \frac{\alpha}{\beta}I\right)\) と \(\beta\) をくくり出すと、

\[ m_N = \beta \left[\beta\left(\Phi^T\Phi + \frac{\alpha}{\beta}I\right)\right]^{-1}\Phi^Ty = \beta \cdot \frac{1}{\beta}\left(\Phi^T\Phi + \frac{\alpha}{\beta}I\right)^{-1}\Phi^Ty = \left(\Phi^T\Phi + \frac{\alpha}{\beta}I\right)^{-1}\Phi^Ty \]

となり、\(\beta\) が完全にキャンセルします。これは正則化パラメータ \(\lambda = \alpha/\beta\) を持つリッジ回帰の解 \(w_{ridge} = (\Phi^T\Phi + \lambda I)^{-1}\Phi^Ty\) と厳密に一致します。

同じ結論は MAP の最適化問題を直接解いても得られます。事後分布を最大化することは \(E(w)\) を最小化することと同値なので、\(m_0=0,S_0=\alpha^{-1}I\) のもとで

\[ E(w) = \frac{\beta}{2}\|y-\Phi w\|^2 + \frac{\alpha}{2}\|w\|^2 \]

を \(w\) で微分してゼロとおくと、

\[ -\beta\Phi^T(y-\Phi w) + \alpha w = 0 \;\;\Longrightarrow\;\; (\beta\Phi^T\Phi + \alpha I)w = \beta\Phi^Ty \;\;\Longrightarrow\;\; w = \left(\Phi^T\Phi + \frac{\alpha}{\beta}I\right)^{-1}\Phi^Ty \]

となり、両辺を \(\beta\) で割った時点で正則化係数が \(\lambda=\alpha/\beta\) であることが確認できます。「ガウス事前分布のMAP推定=リッジ回帰」という等価性は、この2通りの経路(一般公式の特殊化/目的関数の直接最小化)のどちらからも同じ \(\lambda=\alpha/\beta\) に帰着することで証明されます。 5節でこの等価性を sklearn.linear_model.Ridge との数値一致として実際に検証します。

4.3 予測分布の導出

新しい入力 \(x_*\) (基底関数ベクトル \(\phi_* = \phi(x_*)\) )に対する予測は、\(w\) の事後分布 \(p(w\mid y)=\mathcal{N}(w\mid m_N, S_N)\) の下であり得るすべての \(w\) について周辺化することで得られます。これを予測分布と呼びます。

\[ p(y_* \mid x_*, y) = \int p(y_* \mid x_*, w)\, p(w \mid y)\, dw, \qquad p(y_* \mid x_*, w) = \mathcal{N}(y_* \mid \phi_*^Tw, \;\beta^{-1}) \]

この積分は、\(w\sim\mathcal{N}(m_N,S_N)\) の線形変換 \(y_* = \phi_*^Tw + \epsilon_*\) (\(\epsilon_*\sim\mathcal{N}(0,\beta^{-1})\) 、\(w\) と独立)の周辺分布を求める問題そのものです。一般に、\(x\sim\mathcal{N}(\mu,\Lambda^{-1})\) かつ \(y\mid x\sim\mathcal{N}(Ax+b, L^{-1})\) ならば、\(y\) の周辺分布は

\[ y \sim \mathcal{N}\left(A\mu+b,\; L^{-1} + A\Lambda^{-1}A^T\right) \]

というガウス周辺化の一般公式が成り立ちます(ビショップ『パターン認識と機械学習』式2.115)。ここで \(A=\phi_*^T\) (\(1\times M\) 行ベクトル)、\(b=0\) 、\(\Lambda^{-1}=S_N\) 、\(L^{-1}=\beta^{-1}\) を代入すると、

\[ p(y_* \mid x_*, y) = \mathcal{N}\left(y_* \;\middle|\; m_N^T\phi_*,\;\; \sigma_N^2(x_*)\right), \qquad \sigma_N^2(x_*) = \beta^{-1} + \phi_*^T S_N \phi_* \]

が得られます。この予測分散が2つの項の和になっている点が重要です。

  • \(\beta^{-1}\) : 観測ノイズそのものの分散。データがどれだけ増えても消えない**既約誤差(アレアトリック不確実性)**です。
  • \(\phi_*^T S_N \phi_*\) : パラメータ \(w\) の事後不確実性 \(S_N\) に由来する項(エピステミック不確実性)。\(S_N^{-1}=S_0^{-1}+\beta\Phi^T\Phi\) はデータが増えるほど(\(\Phi^T\Phi\) の蓄積により)大きくなるため、\(S_N\) 自体は単調に縮小していきます。したがって \(N\to\infty\) で \(\phi_*^TS_N\phi_*\to 0\) となり、予測分散は既約誤差の下限 \(\beta^{-1}\) に漸近します。

また、\(\phi_*^TS_N\phi_*\) の大きさは \(\phi_*\) の方向にも依存します。訓練データが密に分布する領域では \(\Phi^T\Phi\) がその方向に大きな固有値を持つため \(S_N\) はその方向で小さくなり、不確実性は抑えられます。逆に訓練データがまばら、あるいは全く存在しない領域(外挿領域)では対応する固有値が小さいままなので \(S_N\) は大きく残り、\(\phi_*\) (多項式基底では外挿点で \(\|\phi_*\|\) 自体も急増する)との積で不確実性が大きく増幅されます。これが「データのない領域で信頼区間が広がる」という直感の数学的な根拠です。5節ではこれを \(N=8\) と \(N=100\) の実データで定量的に確認します。

5. Pythonによる数値検証

ここまでの導出が正しいことを、実際にコードを書いて検証します。基底関数は既存の実験結果(6節)と揃えて 9次の多項式 \(\phi_j(x)=x^j\ (j=0,\dots,9)\) を使います。

import numpy as np
from numpy.linalg import inv
from sklearn.linear_model import BayesianRidge, Ridge

DEGREE = 9

def design_matrix(x):
    x = np.asarray(x).reshape(-1)
    return np.vstack([x**j for j in range(DEGREE + 1)]).T  # shape (N, 10)

def true_function(x):
    return np.sin(2 * np.pi * x)

def make_data(n, noise_std, seed):
    rng = np.random.RandomState(seed)
    x = np.sort(rng.uniform(0, 1, size=n))
    y = true_function(x) + rng.normal(0, noise_std, size=n)
    return x, y

def posterior(Phi, y, alpha, beta, m0=None, S0=None):
    """4.1節の一般公式をそのまま実装:
       S_N^{-1} = S0^{-1} + beta * Phi^T Phi
       m_N = S_N (S0^{-1} m0 + beta * Phi^T y)"""
    M = Phi.shape[1]
    if S0 is None:
        S0 = np.eye(M) / alpha
    if m0 is None:
        m0 = np.zeros(M)
    S0_inv = inv(S0)
    S_N_inv = S0_inv + beta * Phi.T @ Phi
    S_N = inv(S_N_inv)
    m_N = S_N @ (S0_inv @ m0 + beta * Phi.T @ y)
    return m_N, S_N

def predictive(phi_star, m_N, S_N, beta):
    """4.3節の予測分布: 平均と分散(2項の和)"""
    mean = phi_star @ m_N
    var = 1.0 / beta + phi_star @ S_N @ phi_star
    return mean, var

(a) sklearnの BayesianRidge との一致確認

sklearn.linear_model.BayesianRidge は最尤(証拠最大化)でノイズ精度・重み精度を推定しますが、推定されたその2つの精度をそのまま4.1節の式に代入すればsklearn 内部の事後平均 coef_ ・事後共分散 sigma_ と一致するはずです(sklearnの命名は逆で、alpha_ がノイズ精度=本記事の \(\beta\) 、lambda_ が重み精度=本記事の \(\alpha\) に対応します)。

N_CHECK, NOISE_STD = 30, 0.2
x_chk, y_chk = make_data(N_CHECK, NOISE_STD, seed=1)   # y = sin(2πx) + noise
Phi_chk = design_matrix(x_chk)

br = BayesianRidge(fit_intercept=False, tol=1e-10, max_iter=1000)
br.fit(Phi_chk, y_chk)
beta_hat, alpha_hat = br.alpha_, br.lambda_             # sklearn命名の入れ替えに注意

m_N_scratch, S_N_scratch = posterior(Phi_chk, y_chk, alpha_hat, beta_hat)

print(np.max(np.abs(m_N_scratch - br.coef_)))    # 事後平均の最大絶対差
print(np.max(np.abs(S_N_scratch - br.sigma_)))   # 事後共分散の最大絶対差

実行結果:

sklearn推定 noise precision (beta)  : 38.033741
sklearn推定 weight precision (alpha): 0.008128
max |m_N(自前実装) - coef_(sklearn)| : 2.300e-11
max |S_N(自前実装) - sigma_(sklearn)|: 2.243e-10

自前実装の事後平均・事後共分散は、sklearnの結果と \(10^{-10}\) オーダーの数値誤差で一致しました。4.1節の導出が正しいことの直接的な裏付けです。

(b) リッジ回帰との等価性の数値検証

4.2節で証明した \(\lambda=\alpha/\beta\) の対応を、固定した \(\alpha=1.0,\ \beta=10.0\) (6節の実験と同じ値)で確認します。

ALPHA_FIX, BETA_FIX = 1.0, 10.0
x15, y15 = make_data(15, 1.0 / np.sqrt(BETA_FIX), seed=2)
Phi15 = design_matrix(x15)

m_N_fix, S_N_fix = posterior(Phi15, y15, ALPHA_FIX, BETA_FIX)

lam = ALPHA_FIX / BETA_FIX  # = 0.1
ridge = Ridge(alpha=lam, fit_intercept=False, solver="cholesky").fit(Phi15, y15)
print(np.max(np.abs(m_N_fix - ridge.coef_)))

実行結果:

alpha=1.0, beta=10.0  =>  lambda = alpha/beta = 0.1
max |m_N - w_ridge|                         : 2.665e-15
(検算) 誤った比 beta/alpha = 10.0 を使うと    : max diff = 1.600e+00  (大きくずれる)

正しい比 \(\lambda=\alpha/\beta=0.1\) を使うと Ridge の係数と**倍精度の丸め誤差レベル(\(10^{-15}\) )**で一致し、誤って逆比 \(\beta/\alpha=10\) を使うと係数が最大 \(1.6\) もずれることが実測で確認できました。これは3節冒頭の注意で述べた「分子分母の取り違え」がいかに結果に影響するかを裏付けています。

下図は \(m_N\) (青)と Ridge の係数(赤)を10個の多項式係数ごとに並べたものです。棒がほぼ完全に重なっており、最大差はタイトルにある通り \(2.7\times10^{-15}\) です。

MAP推定の事後平均とsklearn Ridgeの係数比較(最大差2.7e-15)

(c) 過学習と予測不確実性の収縮:N=8 と N=100

4.3節で「予測分散はデータが増えるほど縮小する」と述べました。これを \(N=8\) (少数データ、過学習しやすい設定)と \(N=100\) (多数データ)で定量的に比較します。ノイズは \(\sigma=0.2\) (\(\beta=1/\sigma^2=25\) )、事前分布は \(\alpha=1.0\) に固定しています。

NOISE_STD, BETA_EXP, ALPHA_EXP = 0.2, 1.0 / 0.2**2, 1.0
x_interior, x_edge = 0.5, 1.1  # 0.5=訓練データ内, 1.1=訓練範囲[0,1]外(外挿)

for n in (8, 100):
    x_n, y_n = make_data(n, NOISE_STD, seed=42)
    Phi_n = design_matrix(x_n)
    m_N_n, S_N_n = posterior(Phi_n, y_n, ALPHA_EXP, BETA_EXP)
    _, var_i = predictive(design_matrix([x_interior])[0], m_N_n, S_N_n, BETA_EXP)
    _, var_e = predictive(design_matrix([x_edge])[0], m_N_n, S_N_n, BETA_EXP)
    print(n, np.sqrt(var_i), np.sqrt(var_e))

実行結果:

N=  8: 予測標準偏差 at x=0.5 (データ内) : 0.22929
N=  8: 予測標準偏差 at x=1.1 (外挿)    : 1.11805
N=100: 予測標準偏差 at x=0.5 (データ内) : 0.20285
N=100: 予測標準偏差 at x=1.1 (外挿)    : 0.52754

収縮率 at x=0.5: 0.22929 / 0.20285 = 1.13倍
収縮率 at x=1.1: 1.11805 / 0.52754 = 2.12倍

既約ノイズ下限 sqrt(1/beta) = 0.20000
N=100でのデータ内予測標準偏差 0.20285 はこの下限にほぼ到達している

訓練データの内側(\(x=0.5\) )では、\(N=8\) でもすでに近傍にデータ点があるため予測標準偏差の差は \(1.13\) 倍とマイルドです。一方、訓練範囲外への外挿点(\(x=1.1\) )では \(N=8\) で \(1.118\) 、\(N=100\) で \(0.528\) と、\(2.12\) 倍もの差が生じています。さらに \(N=100\) のデータ内予測標準偏差 \(0.20285\) は、既約ノイズの理論下限 \(\sqrt{1/\beta}=0.2\) にほぼ到達しており、4.3節で述べた「\(N\to\infty\) で予測分散はエピステミックな項が消えて既約誤差 \(\beta^{-1}\) に漸近する」という主張と整合します。

下図は \(N=8\) (左)と \(N=100\) (右)それぞれについて、予測平均(青実線)と \(\pm2\sigma\) の信頼区間(青帯)、真の関数(灰破線)、訓練データ(黒点)を重ねたものです。灰色の帯は訓練範囲 \([0,1]\) の外(外挿領域)を示します。\(N=8\) では外挿領域で信頼区間が大きく広がり、真の関数から予測平均が大きく外れているのに対し、\(N=100\) では同じ外挿領域でも信頼区間の広がりが明確に小さいことが視覚的にも確認できます。

N=8とN=100でのベイズ予測分布の比較。データが増えると外挿領域の不確実性が大きく縮小する

6. 実験結果(既存の可視化)

  • データ: \(y = \sin(2\pi x)\) からノイズを加えて生成した15点の訓練データ
  • 基底関数: 9次の多項式 (\(f_j(x) = x^j, j=0, ..., 9\) )
  • ハイパーパラメータ: \(\alpha=1.0, \beta=10.0\) (5節(b)の数値検証と同じ設定)

最小二乗推定 / 最尤推定

訓練データに強く適合しようとするため、データのない領域で予測が大きく振動し、過学習を起こしています。 最小二乗推定による回帰結果(過学習の例)

MAP推定

事前分布(正則化)の効果により、過学習が抑制され、より滑らかな予測曲線が得られています。 MAP推定による回帰結果(正則化で過学習を抑制)

ベイズ推定

MAP推定と同様に滑らかな予測(平均値、実線)が得られると同時に、訓練データが少ない領域では予測の不確実性が増大し、信頼区間(青い影の領域)が広がっていることがわかります。 ベイズ推定による回帰結果(信頼区間付き)

モデルの評価(決定係数 \(R^2\) )

決定係数は、モデルの当てはまりの良さを示す指標(1に近いほど良い)です。MAP推定とベイズ推定(の平均)が、最小二乗法よりも高いスコアを示しています。 各推定手法の決定係数R2の比較

ハイパーパラメータ \(\beta\) の影響

ベイズ推定において、尤度の精度パラメータ \(\beta\) (ノイズの逆分散)を大きくすると、モデルは訓練データにより強く適合しようとします。値を大きくしすぎると、信頼区間が狭まり、過学習に近づいていく様子がわかります。

  • \(\beta=50\) ベイズ推定 beta=50のときの予測結果
  • \(\beta=100\) ベイズ推定 beta=100のときの予測結果
  • \(\beta=1000\) ベイズ推定 beta=1000のときの予測結果

まとめ

  • 最小二乗法・最尤推定: シンプルで高速だが、過学習しやすい。
  • MAP推定: 事前分布(正則化)を導入することで、過学習を抑制できる。ガウス事前分布のMAP推定は、正則化係数 \(\lambda=\alpha/\beta\) のリッジ回帰と厳密に等価である(4.2節で証明、5節(b)で \(10^{-15}\) 精度の数値一致を確認)。
  • ベイズ推定: パラメータの不確実性を事後共分散 \(S_N\) として保持し、予測分布の分散 \(\sigma_N^2(x_*)=\beta^{-1}+\phi_*^TS_N\phi_*\) を通じて「既約ノイズ」と「パラメータ不確実性」を分離して定量化できる。データが増えるほど後者は縮小し、特に訓練データの外側(外挿領域)でその効果が顕著(実測で2.12倍の収縮)。

線形回帰のパラメータ空間で行った今回のベイズ推定は、基底関数を無限次元に飛ばし、関数そのものに直接ガウス過程の事前分布を置く関数空間の視点として一般化できます。この立場は ガウス過程回帰の基礎とPython実装 で扱っています。

参考文献