ガウス過程｜ベイズ統計

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：ベイズ線形回帰　|　関連：正規正規モデル

要点（BLUF）

ガウス過程（GP）は、関数そのものに事前を置くノンパラメトリックな手法。任意の有限個の点での関数値が多変量正規に従うと定める。カーネルが滑らかさ・相関を決めます。
GP 回帰の事後も閉形式（多変量正規の条件付き）。予測平均はデータを補間し、予測帯はデータから離れるほど事前の不確実性へ広がる。
ベイズ線形回帰（ベイズ線形回帰）の特徴を無限次元にした極限。カーネルトリックで明示的な特徴ベクトルを使わず、内積（カーネル）だけで計算します。

1. 関数への事前：ガウス過程

これまでは有限個のパラメータ（係数 $w$ など）に事前を置きました。GP は発想が違い、関数 $f$ 全体に事前を置きます。

定義：関数 $f$ がガウス過程 $f\sim\mathcal{GP}(m(x),k(x,x'))$ に従うとは、任意の有限個の入力 $x_1,\dots,x_n$ について、関数値のベクトル $[f(x_1),\dots,f(x_n)]$ が多変量正規 $\mathcal N(\mathbf m,K)$ に従うこと（ $K_{ij}=k(x_i,x_j)$ ）。

平均関数 $m(x)$ は通常 $0$ 。鍵は共分散関数（カーネル） $k(x,x')$ で、「近い入力の関数値ほど強く相関する」という滑らかさの仮定を表します。

2. RBF カーネルと滑らかさ

最もよく使う RBF カーネル（ガウスカーネル）：

k(x,x')=\sigma_f^2\exp\!\left(-\frac{(x-x')^2}{2\ell^2}\right)

$\ell$ （length scale）：相関が効く距離。小さいと細かくうねる関数、大きいとなだらかな関数。
$\sigma_f$ （振幅）：関数の縦の変動幅。

$x=x'$ で $k=\sigma_f^2$ （最大）、離れると指数的に $0$ へ。だから近い点は強く連動し、遠い点は独立——これが「滑らかな関数」の事前になります。

3. GP 回帰の事後（閉形式）

観測 $y_i=f(x_i)+\varepsilon_i,\ \varepsilon\sim\mathcal N(0,\sigma_n^2)$ 。新しい点 $X_*$ での予測は、訓練点と $X_*$ の関数値の同時正規分布を条件付けるだけ——多変量正規の条件付き公式（正規正規モデルの多変量版）で閉形式です。

\bar f_*=K_{*}\,(K+\sigma_n^2 I)^{-1}y,\qquad \mathrm{Cov}(f_*)=K_{**}-K_{*}\,(K+\sigma_n^2 I)^{-1}K_{*}^\top

ここで $K=k(X,X)$ 、 $K_*=k(X_*,X)$ 、 $K_{**}=k(X_*,X_*)$ 。実装します。

import numpy as np

def rbf(a, b, ell=1.0, sf=1.0):                  # RBFカーネル
    d2 = (a[:,None]-b[None,:])**2
    return sf**2*np.exp(-0.5*d2/ell**2)

rng = np.random.default_rng(0)
f = lambda x: np.sin(2*x) + 0.3*x
Xtr = np.array([-3, -2, -1.2, 0.5, 1.5, 2.5, 3.5])    # 訓練入力
sigma_n = 0.1
ytr = f(Xtr) + rng.normal(0, sigma_n, len(Xtr))
Xte = np.linspace(-5, 6, 200)

# 事後：mean=K*(K+σ²I)⁻¹y, var=diag(K**) - 行ごとの K*(K+σ²I)⁻¹K*ᵀ
K  = rbf(Xtr, Xtr) + sigma_n**2*np.eye(len(Xtr))
Ks = rbf(Xte, Xtr); Kss = rbf(Xte, Xte)
L = np.linalg.cholesky(K)
alpha = np.linalg.solve(L.T, np.linalg.solve(L, ytr))
mean = Ks @ alpha
v = np.linalg.solve(L, Ks.T)
sd = np.sqrt(np.maximum(np.diag(Kss) - np.sum(v**2, axis=0), 0))

def at(x0):
    i = np.argmin(np.abs(Xte - x0)); return mean[i], sd[i]
for x0 in [0.5, 0.0, 5.5]:
    m, s = at(x0)
    tag = "訓練点上" if min(abs(Xtr-x0))<0.1 else ("データ内" if -3<=x0<=3.5 else "外挿")
    print(f"x={x0:>4}: 予測平均={m:+.3f} 予測SD={s:.3f} （{tag}）")
print(f"訓練点の平均SD={np.mean([at(x)[1] for x in Xtr]):.3f}（≈雑音{sigma_n}）, 外挿x=5.5でSD={at(5.5)[1]:.3f}")

出力：

x= 0.5: 予測平均=+0.955 予測SD=0.101 （訓練点上）
x= 0.0: 予測平均=+0.045 予測SD=0.261 （データ内）
x= 5.5: 予測平均=+0.451 予測SD=0.985 （外挿）
訓練点の平均SD=0.100（≈雑音0.1）, 外挿x=5.5でSD=0.985

出力の意味：訓練点の上では予測 SD が雑音レベル $0.10$ まで縮み（データを信じて補間）、訓練点の間（ $x=0$ ）では $0.26$ とやや広く、データの全くない外挿（ $x=5.5$ ）では $0.99$ ＝事前の振幅 $\sigma_f=1$ へ戻ります。GP は「知っている所では自信を持ち、知らない所では正直に分からないと言う」——ベイズ線形回帰（ベイズ線形回帰）の予測帯が外挿で広がったのと同じ性質を、関数の形を仮定せずに実現します。

4. 図：補間と外挿の不確実性

事後平均・ $\pm2$ SD 帯と、事後からサンプルした関数を描きます。

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

def rbf(a, b, ell=1.0, sf=1.0):
    return sf**2*np.exp(-0.5*(a[:,None]-b[None,:])**2/ell**2)
rng = np.random.default_rng(0)
f = lambda x: np.sin(2*x) + 0.3*x
Xtr = np.array([-3,-2,-1.2,0.5,1.5,2.5,3.5]); sn = 0.1
ytr = f(Xtr) + rng.normal(0, sn, len(Xtr)); Xte = np.linspace(-5, 6, 200)
K = rbf(Xtr,Xtr)+sn**2*np.eye(len(Xtr)); Ks = rbf(Xte,Xtr); Kss = rbf(Xte,Xte)
L = np.linalg.cholesky(K); alpha = np.linalg.solve(L.T, np.linalg.solve(L, ytr))
mean = Ks@alpha; V = np.linalg.solve(L, Ks.T)
cov = Kss - V.T@V; sd = np.sqrt(np.maximum(np.diag(cov), 0))

plt.figure(figsize=(8, 4.5))
for s in rng.multivariate_normal(mean, cov + 1e-9*np.eye(len(Xte)), 8):
    plt.plot(Xte, s, color="C1", alpha=0.25, lw=0.8)
plt.fill_between(Xte, mean-2*sd, mean+2*sd, alpha=0.2, color="C0", label="予測 ±2SD")
plt.plot(Xte, mean, color="C0", lw=2, label="予測平均")
plt.scatter(Xtr, ytr, color="black", zorder=5, label="データ")
plt.xlabel("x"); plt.ylabel("f(x)"); plt.legend(fontsize=9)
plt.title("ガウス過程回帰：データを補間し外挿で帯が広がる")
plt.tight_layout(); plt.show()

グラフの意味：事後からサンプルした関数（橙）はデータ点を通り、データのある $[-3,3.5]$ では互いに揃って予測帯（青）が細く、外側へ出ると関数がばらけて帯が事前の幅へ広がります。関数の明示式を仮定せず、カーネルの滑らかさだけでこの推論ができるのが GP の強みです。

5. ベイズ線形回帰との関係：無限次元の極限

GP はベイズ線形回帰（ベイズ線形回帰）と地続きです。基底関数 $\phi(x)$ を使った線形回帰 $f(x)=\phi(x)^\top w$ 、 $w\sim\mathcal N(0,\Sigma_p)$ では、 $f$ は GP になり、そのカーネルは $k(x,x')=\phi(x)^\top\Sigma_p\phi(x')$ 。基底を無限個にした極限が RBF カーネル GP に対応します。GP は明示的な $\phi$ を持たず、カーネル（内積）だけで計算するので、無限次元の特徴を扱えます（カーネルトリック）。

ハイパーパラメータ（ $\ell,\sigma_f,\sigma_n$ ）は、周辺尤度 $p(y\mid X)=\mathcal N(y\mid 0,K+\sigma_n^2I)$ を最大化して学習します（経験ベイズ＝経験ベイズの GP 版）。計算コストは $K$ の逆行列で $O(n^3)$ ——大規模データでは疎近似・誘導点などの手法が要ります（この領域は実装が速く動くので要最新確認）。

⚠️ よくある誤解

「GP はパラメータがない」ではない：関数を無限次元のパラメータと見るノンパラメトリック。代わりにカーネルのハイパーパラメータ（ $\ell$ 等）を学習します。
「予測帯はどこでも一定」ではない：データの近くで細く、外挿で事前へ広がります（§3）。これが GP の不確実性定量化の核心。
「GP は重くない」ではない： $O(n^3)$ で大規模データに弱い。疎近似が要ります。
「カーネルは何でもよい」ではない：カーネルが事前の仮定（滑らかさ・周期性など）を決めます。問題に合うカーネル選択が肝心。