機械学習と情報技術

ガウス過程回帰の理論と実装をわかりやすく解説

ガウス過程

ガウス過程回帰(Gaussian Process Regression, GPR)は、ベイズ推定の枠組みでカーネル法を用いた回帰手法です。予測値だけでなく予測の不確実性（分散）も同時に得られる点が大きな特徴で、ベイズ最適化やアクティブラーニングなどの応用でも広く利用されています。

本記事では、ガウス過程の定義からカーネル関数の役割、予測分布の導出まで丁寧に解説し、Pythonでスクラッチ実装します。

本記事の内容

  • ガウス過程の定義と直感的な理解
  • カーネル関数の役割
  • 予測分布の導出
  • Pythonでのスクラッチ実装と可視化

前提知識

この記事を読む前に、以下の概念を押さえておくと理解が深まります。

ガウス過程とは

ガウス過程(Gaussian Process, GP)とは、任意の有限個の点の関数値が多変量正規分布に従うような確率過程です。

大雑把に言うと、「関数そのもの」に対する確率分布を定義したものです。通常の正規分布は数値に対する分布ですが、ガウス過程は関数空間上の分布と考えることができます。

ガウス過程は平均関数 $m(\bm{x})$ とカーネル関数（共分散関数）$k(\bm{x}, \bm{x}’)$ によって特定されます。

$$ f(\bm{x}) \sim \mathcal{GP}(m(\bm{x}), k(\bm{x}, \bm{x}’)) $$

これは、任意の有限個の入力点 $\bm{x}_1, \dots, \bm{x}_n$ に対して、

$$ \begin{pmatrix} f(\bm{x}_1) \\ f(\bm{x}_2) \\ \vdots \\ f(\bm{x}_n) \end{pmatrix} \sim \mathcal{N}\left( \begin{pmatrix} m(\bm{x}_1) \\ m(\bm{x}_2) \\ \vdots \\ m(\bm{x}_n) \end{pmatrix}, \begin{pmatrix} k(\bm{x}_1, \bm{x}_1) & \cdots & k(\bm{x}_1, \bm{x}_n) \\ \vdots & \ddots & \vdots \\ k(\bm{x}_n, \bm{x}_1) & \cdots & k(\bm{x}_n, \bm{x}_n) \end{pmatrix} \right) $$

が成り立つことを意味します。

カーネル関数

カーネル関数は、2つの入力点がどれくらい「似ているか」を定量化する関数です。最もよく使われるのはRBF(Radial Basis Function)カーネルです。

$$ k(\bm{x}, \bm{x}’) = \sigma_f^2 \exp\left( -\frac{|\bm{x} – \bm{x}’|^2}{2\ell^2} \right) $$

ここで、

  • $\sigma_f^2$: シグナル分散（関数値の振幅を制御）
  • $\ell$: 長さスケール（入力空間での「滑らかさ」を制御）

$\ell$ が小さいと関数は急激に変化し、大きいと滑らかな関数になります。

予測分布の導出

訓練データ $\mathcal{D} = \{(\bm{x}_i, y_i)\}_{i=1}^{n}$ が与えられたとき、新しい入力点 $\bm{x}_*$ における関数値 $f_*$ の予測分布を求めます。

観測モデルとして $y = f(\bm{x}) + \epsilon$, $\epsilon \sim \mathcal{N}(0, \sigma_n^2)$ を仮定します。

訓練データの出力 $\bm{y}$ と予測値 $f_*$ の同時分布は、

$$ \begin{pmatrix} \bm{y} \\ f_* \end{pmatrix} \sim \mathcal{N}\left( \bm{0}, \begin{pmatrix} \bm{K} + \sigma_n^2 \bm{I} & \bm{k}_* \\ \bm{k}_*^T & k_{**} \end{pmatrix} \right) $$

ここで、

  • $\bm{K}$: 訓練データ間のカーネル行列 ($K_{ij} = k(\bm{x}_i, \bm{x}_j)$)
  • $\bm{k}_*$: 訓練データとテスト点間のカーネルベクトル ($k_{*i} = k(\bm{x}_i, \bm{x}_*)$)
  • $k_{**} = k(\bm{x}_*, \bm{x}_*)$

多変量正規分布の条件付き分布の公式から、予測分布は次のようになります。

$$ f_* | \bm{X}, \bm{y}, \bm{x}_* \sim \mathcal{N}(\bar{f}_*, \text{Var}(f_*)) $$

$$ \bar{f}_* = \bm{k}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{y} $$

$$ \text{Var}(f_*) = k_{**} – \bm{k}_*^T (\bm{K} + \sigma_n^2 \bm{I})^{-1} \bm{k}_* $$

予測平均 $\bar{f}_*$ は訓練データの $y$ 値の重み付き線形結合であり、予測分散 $\text{Var}(f_*)$ はデータ点から離れるほど大きくなります。

Pythonでの実装

ガウス過程回帰をスクラッチ実装します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- カーネル関数（RBFカーネル）---
def rbf_kernel(X1, X2, sigma_f=1.0, length_scale=1.0):
    """RBFカーネル行列を計算"""
    sq_dist = np.sum(X1**2, axis=1).reshape(-1, 1) + \
              np.sum(X2**2, axis=1).reshape(1, -1) - \
              2 * X1 @ X2.T
    return sigma_f**2 * np.exp(-sq_dist / (2 * length_scale**2))

# --- ガウス過程回帰 ---
def gp_predict(X_train, y_train, X_test, sigma_f=1.0, length_scale=1.0, sigma_n=0.1):
    """
    ガウス過程回帰の予測分布を計算
    Returns: 予測平均, 予測分散
    """
    K = rbf_kernel(X_train, X_train, sigma_f, length_scale) + sigma_n**2 * np.eye(len(X_train))
    K_s = rbf_kernel(X_train, X_test, sigma_f, length_scale)
    K_ss = rbf_kernel(X_test, X_test, sigma_f, length_scale)

    # コレスキー分解で安定的に逆行列を計算
    L = np.linalg.cholesky(K)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_train))
    v = np.linalg.solve(L, K_s)

    # 予測平均と予測分散
    mu_pred = K_s.T @ alpha
    var_pred = np.diag(K_ss - v.T @ v)

    return mu_pred, var_pred

# --- テストデータ生成 ---
# 真の関数
def true_function(x):
    return np.sin(x) + 0.5 * np.sin(3 * x)

# 訓練データ
n_train = 10
X_train = np.sort(np.random.uniform(-5, 5, n_train)).reshape(-1, 1)
y_train = true_function(X_train.ravel()) + np.random.randn(n_train) * 0.2

# テストデータ
X_test = np.linspace(-6, 6, 200).reshape(-1, 1)
y_true = true_function(X_test.ravel())

# --- 予測 ---
mu_pred, var_pred = gp_predict(X_train, y_train, X_test,
                                sigma_f=1.0, length_scale=1.0, sigma_n=0.2)
std_pred = np.sqrt(var_pred)

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 予測結果
ax1 = axes[0]
ax1.plot(X_test, y_true, 'k--', label='True function', alpha=0.5)
ax1.plot(X_test, mu_pred, 'b-', label='GP mean')
ax1.fill_between(X_test.ravel(),
                 mu_pred - 2 * std_pred,
                 mu_pred + 2 * std_pred,
                 alpha=0.2, color='blue', label='95% confidence')
ax1.scatter(X_train, y_train, c='red', s=50, zorder=5, label='Training data')
ax1.set_xlabel('x')
ax1.set_ylabel('f(x)')
ax1.set_title('Gaussian Process Regression')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 長さスケールの影響
ax2 = axes[1]
for ls in [0.3, 1.0, 3.0]:
    mu, var = gp_predict(X_train, y_train, X_test,
                          sigma_f=1.0, length_scale=ls, sigma_n=0.2)
    ax2.plot(X_test, mu, label=f'l = {ls}')
ax2.scatter(X_train, y_train, c='red', s=50, zorder=5)
ax2.set_xlabel('x')
ax2.set_ylabel('f(x)')
ax2.set_title('Effect of Length Scale')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

このコードでは、RBFカーネルを用いたガウス過程回帰をスクラッチ実装しています。予測平均（青線）に加えて95%信頼区間（薄青の帯）が得られ、データ点から離れるほど不確実性が大きくなることがわかります。

まとめ

本記事では、ガウス過程回帰の理論と実装について解説しました。

  • ガウス過程は「関数空間上の確率分布」であり、平均関数とカーネル関数で特定される
  • RBFカーネルの長さスケール $\ell$ が関数の滑らかさを制御する
  • 予測分布は解析的に求まり、予測平均と予測分散を同時に得られる
  • データ点から離れると予測の不確実性が大きくなるという直感的な振る舞いをする

次のステップとして、以下の記事も参考にしてください。

関連記事

前の記事

確率空間とは？標本空間・事象・確率測度をわかりやすく解説

次の記事

【Python】ローレンツ方程式を実装して可視化する

関連タグ: カーネル法 ガウス過程 ベイズ推定 回帰 機械学習