ガウス過程回帰の計算手法 - はじめよう実験計画

この記事について

本記事では、ガウス過程回帰（Gaussian Process Regression）を数式とRコードを用いて、実際に計算しながら理解します。

ガウス過程回帰とは

ガウス過程回帰は「データが多変量正規分布に従う」というシンプルな仮定から、予測値の平均と不確実性（分散）を同時に計算できる手法です。プロセス開発・材料開発・実験自動化などができるベイズ最適化の、内部モデルとしてよく使われます。私自身、ガウス過程回帰を半導体の材料開発で活用してきました。

本記事では以下の内容を解説します：

ガウス過程モデルの定式化
カーネル行列の計算方法と意味
ガウス過程回帰を計算するためのRコード
ガウス過程回帰の外挿における問題点

統計学の基礎知識（正規分布・行列演算）があれば読めるように書きました。

サンプルデータ

本記事では、下図のサンプルデータ（黒点）を例に計算を進めます。ここで、xは入力変数、yは観測値です。

このデータに対して、ガウス過程回帰モデルを適用し、下図のような予測分布を計算できるようにするのが、本記事の目標です。

図3 GP による予測（青線=予測平均、青帯=±1SD、黒線=真の関数、黒丸=訓練データ）。オレンジ縦線の外側が外挿領域。

Rコード

外部パッケージを使用せず、すべての行列計算を実装していますので、ガウス過程回帰の中身がちゃんとわかる内容になっています。Rコード全文は本記事の付録に含まれていますので、Rの実行環境がある場合はコピペで動かすことができます。

また、Rの実行環境がない場合もShinyliveを使って、こちらのリンクから試すこともできます。

記事の構成

1章ガウス過程回帰とは
2章新しい点での予測値
3章カーネルの計算方法と意味
4章予測分布を実際に計算
5章ガウス過程回帰を計算するRコード
6章外挿の危険性について

最終章の外挿の危険性については、こちらの記事でいかに克服するかを解説しましたので、合わせてお読みください。
www.doe-get-started.com

1. ガウス過程回帰とは

ガウス過程回帰（GP）はデータが「一つの大きな多変量正規分布」に従うというモデルです。

既知の観測値 $\mathbf{y}=[y_1,\dots,y_N]^T$ と予測値 $\mathbf{y_*}=[y_{1*},\dots,y_{M*}]^T$ をまとめると、ガウス過程は、

$\begin{bmatrix}y_1\\\vdots\\y_N\\y_{1*}\\\vdots\\y_{M*}\end{bmatrix}\sim N\!\left(\begin{bmatrix}0\\\vdots\\0\end{bmatrix},\begin{bmatrix}K&k_*\\k_*^T&k_{**}\end{bmatrix}\right)$

・・・(eq.1)

という形で、記述されます。

ここで、 $K$ 、 $k_*$ 、 $k_{**}$ はカーネルと言われる行列で、次節で説明するように、すべて入力 $[\mathbf{X}, \mathbf{X_*}]=[x_1,\dots,x_N,x_{1*},\dots,x_{M*}]^T$ から決まります（観測値 $\mathbf{y}$ には依存しません）。

各カーネル行列の意味：

$K$ : 訓練点どうしのカーネル行列（ $N\times N$ ）
$k_*$ : 訓練点と予測点のカーネル行列（ $N\times M$ ）
$k_{**}$ : 予測点どうしのカーネル行列（ $M\times M$ ）

2. 新しい点での予測値

今やりたいことは、既知のデータセット $D=\{(x_1,y_1),\dots,(x_N,y_N)\}$ が与えられた時、新しい点 $X_*$ における予測値 $y_*$ を計算することです。予測は条件付き分布(eq.2)として計算されます。この式が一番重要です。

$y_*|X^*,D\sim N\bigl(k_*^T(K+\sigma^2 I)^{-1}y,\;k_{**}-k_*^T(K+\sigma^2 I)^{-1}k_*\bigr)$

・・・(eq.2)

難解に思われるかもしれませんが、(eq.2)は(eq.1)に対して条件付き分布の"公式"を当てはめることで導出できるので、気になる方は付録を読んでください。

一旦(eq.2)を認めていただくとして、ここで言いたいことは、 $X_*, D$ を決めると予測 $y_*$ の

平均 $k_*^T(K+\sigma^2 I)^{-1}y$
分散 $k_{**}-k_*^T(K+\sigma^2 I)^{-1}k_*$

は一意に定まるということです。

3. カーネル

次に、(eq.2)にカーネル $K$ 、 $k_*$ 、 $k_{**}$ をどのように計算するのかを説明します。
ガウス過程回帰ではRBFカーネルというカーネルがよく用いられます。

$k_{\mathrm{RBF}}(x_i,x_j)=\theta_1\exp\!\left(-\frac{(x_i-x_j)^2}{2\theta_{2}}\right)$

・・・(eq.3)

これはカーネル行列の各成分が、2つの入力点 $x_i, x_j$ とハイパーパラメータ $\theta_1, \theta_2$ を使って計算できます。

すべての $x_i, x_j$ の組に対して $K_{ij} = k_{\mathrm{RBF}}(x_i,x_j)$ を計算しカーネル行列 $K$ を構築していきます。

ここまでで、まだカーネルの意味はさっぱり分からないと思いますが、先にどんな計算をするかだけ紹介します。

3.1 カーネルの計算例

入力点 $(x_1, x_2, x_3, x_4, x_5)= (-1, -0.5, 0, 0.5, 1)$ として、RBFでカーネルを計算してみます。繰り返しになりますが、カーネルは入力点 $x$ だけで決まるもので、出力 $y$ は必要ありません。

そのままxに対して計算してもよいのですが、後々のために、xを次の式で標準化し、

$\displaystyle{ z^s = \frac{x-x_{\mathrm{mean}}}{\sigma_x} }$

$z^s = [-1.2649, -0.6325, 0.0000, +0.6325, +1.2649]$ に対してカーネルを計算してみましょう。RBFカーネルのハイパーパラメータは固定して、 $\theta_1=1.0, \quad \theta_2=0.4$ とします。

$\begin{aligned} k(z^s_{1}, z^s_{2}) &= \theta_1 \exp\!\left(-\frac{(z^s_{1}-z^s_{2})^2}{2\,\theta_2}\right) \\ &= 1.0 \times \exp\!\left(-\frac{(-1.2649 - (-0.6325))^2}{2 \times 0.4}\right) \\ &= \exp\!\left(-\frac{0.4000}{0.8}\right) \\ &= \exp(-0.5000) \\ &= 0.6065 \end{aligned}$

となります。すべての $(z^s_i, z^s_j)$ に対して同様に計算すると、

$K = \begin{bmatrix} 1.0000 & 0.6065 & 0.1353 & 0.0111 & 0.0003 \\ 0.6065 & 1.0000 & 0.6065 & 0.1353 & 0.0111 \\ 0.1353 & 0.6065 & 1.0000 & 0.6065 & 0.1353 \\ 0.0111 & 0.1353 & 0.6065 & 1.0000 & 0.6065 \\ 0.0003 & 0.0111 & 0.1353 & 0.6065 & 1.0000 \end{bmatrix}$

これが(eq.2)で用いるカーネル $K$ です。

以下に、カーネル行列を構築するRコードを示します。RBF以外も選べるようになっています。

# ---------- カーネル成分の計算(rbf以外も選択可能) ----------

kernel <- function(xi, xj, type = c("rbf", "linear", "quad", "rbf_linear"),

                   theta1 = 1, theta2 = 0.4) {

  type <- match.arg(type)

  switch(type,

         rbf = { return(theta1 * exp(-0.5*sum((xi - xj)^2) / theta2)) },

         linear    = { return(sum(xi*xj) + 1) },

         quad      = { return((sum(xi*xj) + 1)^2) },

         rbf_linear = { return(theta1 * exp(-0.5*sum((xi - xj)^2) / theta2) + sum(xi*xj) + 1) }

  )

}

# ---------- カーネル行列の構築 ----------

kernel_matrix <- function(X1, X2 = X1, type = "rbf", theta1 = 1, theta2 = 0.4,

                          noise = 0, jitter = 1e-5) {

  X1 <- as.matrix(X1); X2 <- as.matrix(X2)

  N1 <- nrow(X1); N2 <- nrow(X2)

  K  <- matrix(0, nrow = N1, ncol = N2)

  for (i in 1:N1) {

    xi <- X1[i, , drop = FALSE]

    for (j in 1:N2) {

      xj <- X2[j, , drop = FALSE]

      K[i, j] <- kernel(xi, xj, type = type, theta1 = theta1, theta2 = theta2)

    }

  }

  if (noise != 0 && identical(X1, X2)) { K <- K + noise * diag(N1) }

  return(K)

}

2.2 カーネルの意味

ここまでで、とりあえずカーネル行列は作れるようになりましたね。しかし、カーネルとは一体何を意味しているのでしょうか？本節はかなり抽象的になるので、飛ばしても問題ありません。

今、ある一つの入力 $x_1$ について無限次元の基底 $\Phi(x_1)=[\phi_0(x_1), \phi_1(x_1), ..., \phi_{\infty}(x_1)]$ を用意し、入力 $x_1$ に対する出力 $y_1$ を

$y_1= \begin{bmatrix}\phi_0(x_1), \phi_1(x_1), \cdots, \phi_{\infty}(x_1)\end{bmatrix} \begin{bmatrix} w_{01} \\ w_{11}\\ \vdots\\w_{1\infty} \end{bmatrix} = \Phi_1 \mathbf{w}_1, \quad\mathbf{w}_1\sim N(\mathbf{0},\mathbf{I})$

と表すことにしましょう。すると全データは、

$\mathbf{y}=[y_1, ..., y_N]^T = \begin{bmatrix}\phi_0(x_1) & \phi_1(x_1) & \cdots \\\phi_0(x_2) & \phi_1(x_2) & \cdots \\\vdots & \vdots & \ddots \\\phi_0(x_N) & \phi_1(x_N) & \cdots\end{bmatrix}\begin{bmatrix}w_1\\w_2\\\vdots\end{bmatrix}=\Phi\mathbf{w}$

のように書けますから、共分散行列は、

$\mathrm{Cov}(\mathbf{y})=E[\mathbf{y}\mathbf{y}^T]-E[\mathbf{y}]E[\mathbf{y}]^T=\Phi E[\mathbf{w}\mathbf{w}^T]\Phi^T=\Phi\Phi^{T}$

となります。これは $\mathbf{y}\sim N(\mathbf{0}, \Phi\Phi^{T})$ というガウス過程です。

以上を一行でまとめると、

$\Phi(x)=[\phi_0(x), \phi_1(x), ..., \phi_{\infty}(x)]$ のような無限次元の基底の線形結合モデル $y=\Phi w$ と、ガウス過程 $y\sim N(0, K), \quad K=\Phi\Phi^{T}$ が同じものだよ

ということです。

別の言い方をすれば、ガウス過程のカーネル $K$ を決めることは、基底関数の性質を決めることだと言えます。

そして、カーネル $K$ を決める際、基底 $\Phi$ を考える必要はない、というのがカーネルトリックと言われる、まさにトリックですね。

RBFカーネルの場合、Kの各成分は(eq.3)で与えられますが、これは背後で
$\phi_n(x)=\exp\!\left(-\frac{x^2}{2\theta_2}\right)\frac{x^n}{\theta_2^{n/2}\sqrt{n!}}$
という関数の無限次元の基底 $\Phi(x)=[\phi_0(x), \phi_1(x), ..., \phi_{\infty}(x)]$ を用いた線形結合モデル $y=\Phi w$ になっています。

4. 予測分布を実際に計算

前置きが長くなりましたので、何をしたいか思い出しておきましょう。我々は、ガウス過程回帰による新しい点の予測値 $y_*$ の平均と分散を、(eq.2)により計算したいのでした。

$y_*|X^*,D\sim N\bigl(k_*^T(K+\sigma^2 I)^{-1}y,\;k_{**}-k_*^T(K+\sigma^2 I)^{-1}k_*\bigr)$

・・・(eq.2 (再掲))

ここで、

$K$ : 訓練点どうしのカーネル行列（ $N\times N$ ）
$k_*$ : 訓練点と予測点のカーネル行列（ $N\times M$ ）
$k_{**}$ : 予測点どうしのカーネル行列（ $M\times M$ ）

です。

また、 $\sigma^2$ はデータが得られるときのノイズを考慮した分散になります。

4.1 サンプルデータ

以下では、サンプルデータを使って、ガウス過程回帰の予測分布を計算していきます。5回の実験を行い、それぞれの入力xにおいて出力y（実験値）が、以下のように得られたとしましょう。

i	入力 $x$	観測値 $y$	標準化 $z^s$	標準化 $y^s$
1	−1.0	−8.200	−1.2649	−1.2785
2	−0.5	−1.125	−0.6325	−0.2212
3	0.0	−2.843	0.0000	−0.4779
4	+0.5	+8.697	+0.6325	+1.2463
5	+1.0	+5.250	+1.2649	+0.7313

これをplotすると図のようになります。

を行っておきます。生データと標準化後のデータは下図の通りです。

4.2 ガウス過程回帰による予測平均と分散

まず、カーネルを構築します。前述したカーネルを計算する関数kernel_matrixを利用しましょう。

# 訓練点どうしのカーネル行列（ノイズを含む: eq.2 における K+sigma2*I）

K <- kernel_matrix(X_train, X_train, type=type, theta1=theta1, theta2=theta2, noise=theta3)

# 訓練点と予測点のカーネル行列（ノイズを含まない）

k <- kernel_matrix(X_train, X_pred, type=type, theta1=theta1, theta2=theta2)

# 予測点どうしのカーネル行列（ノイズを含まない）

s <- kernel_matrix(X_pred, X_pred, type=type, theta1=theta1, theta2=theta2)

RBFカーネルのハイパーパラメータは固定して、theta1=1.0、theta2=0.4、sigma2=0.05 とすると

標準化後の入力 $z^s$ に対するK：

$K+\sigma^2 I = \begin{bmatrix} 1.0500 & 0.6065 & 0.1353 & 0.0111 & 0.0003 \\ 0.6065 & 1.0500 & 0.6065 & 0.1353 & 0.0111 \\ 0.1353 & 0.6065 & 1.0500 & 0.6065 & 0.1353 \\ 0.0111 & 0.1353 & 0.6065 & 1.0500 & 0.6065 \\ 0.0003 & 0.0111 & 0.1353 & 0.6065 & 1.0500 \end{bmatrix}$

4.2.1 内挿点 x*=0 における予測値

Step 1: 標準化（訓練データと同じ統計量を使用）
$x_{mean} = 0.0000, \quad \sigma_x = 0.7906$
$y_{mean} = 0.3560, \quad \sigma_y = 6.6926$
$z^s_* = (0.0 - 0.0) / 0.7906 = 0.0000$

Step 2: $k_*$ の計算（ $z^s_*=0$ と各訓練点 $z^s$ の間のカーネル値）

$k_* = \begin{bmatrix}0.1353 \\ 0.6065 \\ 1.0000 \\ 0.6065 \\ 0.1353\end{bmatrix},\quad k_{**} = k(z^s_*,z^s_*) = 1.0000$

Step 3: $(K+\sigma^2 I)^{-1}y$ の計算

$(K+\sigma^2 I)^{-1}y\quad(= \mathtt{yy}とおく)$ を求めます。

$\mathtt{yy} = \begin{bmatrix} 1.0500 & 0.6065 & 0.1353 & 0.0111 & 0.0003 \\ 0.6065 & 1.0500 & 0.6065 & 0.1353 & 0.0111 \\ 0.1353 & 0.6065 & 1.0500 & 0.6065 & 0.1353 \\ 0.0111 & 0.1353 & 0.6065 & 1.0500 & 0.6065 \\ 0.0003 & 0.0111 & 0.1353 & 0.6065 & 1.0500 \end{bmatrix}\begin{bmatrix}-1.2785\\-0.2212\\-0.4779\\+1.2463\\+0.7313\end{bmatrix} = \begin{bmatrix}-2.4666\\+2.9333\\-3.7414\\+3.5010\\-0.8739\end{bmatrix}$

Step 4: 予測平均の計算

$\begin{aligned} \mu^s_* &= k_*^T(K+\sigma^2 I)^{-1}y = k_*^T \cdot \mathtt{yy} \\ &= 0.1353\times(-2.4666)+0.6065\times 2.9333+1.0000\times(-3.7414) \\ &\quad+0.6065\times 3.5010+0.1353\times(-0.8739) \\ &= -0.3338+1.7791-3.7414+2.1235-0.1183 \\ &= -0.2909 \quad\text{（標準化空間）} \end{aligned}$

Step 5: 予測分散 $\mathrm{var}^s_y$ の計算

$\begin{aligned} \mathrm{var}^s_{f*} &= k_{**}-k_*^T(K+\sigma^2 I)^{-1}k_* = 1.0000 - 0.9574 = 0.0426 \\ \mathrm{var}^s_{\mathrm{y}*} &= v_{\mathrm{gp}} + \sigma^2 = 0.0426 + 0.05 = 0.0926 \end{aligned}$

Step 6: 元スケールへの逆変換

$\begin{aligned} \mu_* &= -0.2909 \times 6.6926 + 0.3560 = -1.591 \\ \mathrm{var}_* &= \sqrt{0.0926} \times 6.6926 = 2.037 \end{aligned}$

よって予測値は-1.591 ± 2.037 となります。

4.2.2 外挿点 x*=2 における予測値

Step 1: 標準化
$z^s_* = (2.0 - 0.0) / 0.7906 = 2.5298$

Step 2: $k_*$ の計算（ $z^s_*=2.5298$ と各訓練点 $z^s$ の間のカーネル値）

$k_* = \begin{bmatrix}0.0000 \\ 0.0000 \\ 0.0003 \\ 0.0111 \\ 0.1353\end{bmatrix},\quad k_{**} = k(z^*,z^*) = 1.0000$

ここで、内挿点 $x*=0$ の時との違いは、 $k_*$ の成分がほぼゼロであることです。 $k_*$ は共分散ですから、外挿点 $x*=2$ と訓練データがほとんど互いに関係がないということを意味します。

Step 3: $\mathtt{yy}$ の計算（x*=0 と共通）

$\mathtt{yy} = \begin{bmatrix}-2.4666\\+2.9333\\-3.7414\\+3.5010\\-0.8739\end{bmatrix}$

Step 4: 予測平均の計算

$\begin{aligned} \mu^s_* &= k_*^T \cdot \mathtt{yy} \\ &= 0.0000\times(-2.4666)+0.0000\times 2.9333+0.0003\times(-3.7414) \\ &\quad+0.0111\times 3.5010+0.1353\times(-0.8739) \\ &= 0.0000+0.0000-0.0013+0.0389-0.1183 \\ &= -0.0806 \quad\text{（標準化空間）} \end{aligned}$

Step 5: 予測分散の計算

$\begin{aligned} \mathrm{var}^s_{f*} &= k_{**}-k_*^T(K+\sigma^2 I)^{-1}k_* = 1.0000 - 0.0261 = 0.9739 \\ \mathrm{var}^s_{\mathrm{y}*} &= \mathrm{var}^s_{f*} + \sigma^2 = 0.9739 + 0.05 = 1.0239 \end{aligned}$

Step 6: 元スケールへの逆変換

$\begin{aligned} \mu_* &= -0.0806 \times 6.6926 + 0.3560 = -0.184 \\ \mathrm{var}_* &= \sqrt{1.0239} \times 6.6926 = 6.772 \end{aligned}$

よって予測値は −0.184 ± 6.772 となります。

以下に、x*=0（内挿）と x*=2（外挿）の比較を示します：

	x*=0（内挿）	x*=2（外挿）
z*	0.0000	2.5298
$k_*$ 最大値	1.0000（中央点）	0.1353（最遠点のみ）
$\mu^s_*$ （標準化）	−0.2909	−0.0806（≈ 事前平均 0）
$\mathrm{var}^s_y$ （標準化）	0.0426（小）	0.9739（≈ 1）
$\mu_*$ （元スケール）	−1.591	−0.184（≈ 事前平均 0.356）
$\mathrm{var}_y$ （元スケール）	2.037	6.772（大）

外挿点では $k_*$ が小さくなるため、 $\mu_*=k_*^T\cdot\mathtt{yy}\approx 0$ となり予測平均が標準化空間でゼロ（元スケールで $\bar{y}=0.356$ ）に収束します。また $\mathrm{var}_y\approx k_{**}=1$ となり予測分散が事前分散まで膨らみます。

5. ガウス過程回帰で新しい点を予測するRコード

以下に、新しい点の予測平均・分散を計算する関数 GP_pred を示します。

# ---------- コレスキー分解により Vx=b を x について解く ----------

chol_solve <- function(L, b) {

  y <- forwardsolve(L, b, upper.tri = FALSE, transpose = FALSE)

  x <- backsolve(t(L), y, upper.tri = TRUE, transpose = FALSE)

  x

}

# ---------- ガウス過程回帰の予測分布(mu, var_y)の計算 ----------

GP_pred <- function(X_train, Y_train, X_pred,

                    theta1 = 1, theta2 = 0.4, theta3 = 0,

                    type = c("rbf", "linear", "quad", "rbf_linear")) {

  X_train <- as.matrix(X_train)

  Y_train <- as.matrix(Y_train)

  X_pred  <- as.matrix(X_pred)

  N <- nrow(X_train); M <- nrow(X_pred)

  if (N != length(Y_train)) stop("length(x_train) must equal length(y_train)")

  # カーネル行列の構築（theta2 = 0.4 をそのまま渡す）

  K <- kernel_matrix(X_train, X_train, type=type, theta1=theta1, theta2=theta2, noise=theta3)

  k <- kernel_matrix(X_train, X_pred,  type=type, theta1=theta1, theta2=theta2)

  s <- kernel_matrix(X_pred,  X_pred,  type=type, theta1=theta1, theta2=theta2)

  # コレスキー分解

  L   <- t(chol(K))

  yy  <- chol_solve(L, Y_train)   # yy = K^{-1} y

  v   <- chol_solve(L, k)         # v  = K^{-1} k_*

  mu    <- t(k) %*% yy            # 予測平均

  var_f <- s - t(k) %*% v         # 予測分散（関数）

  if (abs(theta3) > .Machine$double.eps^0.5) {

    var_y <- var_f + theta3 * diag(M)

    return(list(mu = mu, var_f = var_f, var_y = var_y))

  } else {

    return(list(mu = mu, var_f = var_f))

  }

}

これまでの関数（kernel_matrixおよびGP_pred）とその他必要な処理すべてを含んだコードを示します。以下のコードをすべてコピペすると、サンプルデータを使ったガウス過程回帰のプロットが得られます。

# -----------------------------------------------------------------------

# Kernel処理関数

# -----------------------------------------------------------------------

chol_solve <- function(L, b) {

  y <- forwardsolve(L, b, upper.tri = FALSE, transpose = FALSE)

  x <- backsolve(t(L), y, upper.tri = TRUE, transpose = FALSE)

  x

}

kernel <- function(xi, xj, type = c("rbf", "linear", "quad", "rbf_linear"),

                   theta1 = 1, theta2 = 0.4) {

  type <- match.arg(type)

  switch(type,

    rbf = { d2 <- sum((xi - xj)^2); return(theta1 * exp(-0.5*d2 / theta2)) },

    linear = { return(sum(xi*xj) + 1) },

    quad = { return((sum(xi*xj) + 1)^2) },

    rbf_linear = { d2 <- sum((xi - xj)^2); return(theta1 * exp(-0.5*d2 / theta2) + sum(xi*xj) + 1) }

  )

}

kernel_matrix <- function(X1, X2 = X1, type = "rbf", theta1 = 1, theta2 = 1, noise = 0, jitter = 1e-5) {

  X1 <- as.matrix(X1); X2 <- as.matrix(X2)

  N1 <- nrow(X1); N2 <- nrow(X2)

  K <- matrix(0, nrow = N1, ncol = N2)

  for (i in 1:N1) {

    xi <- X1[i, , drop = FALSE]

    for (j in 1:N2) { K[i, j] <- kernel(xi, X2[j, , drop = FALSE], type = type, theta1 = theta1, theta2 = theta2) }

  }

  if (noise != 0 && identical(X1, X2)) K <- K + noise * diag(N1)

  return(K)

}

# -----------------------------------------------------------------------

# Gaussian Process処理関数

# -----------------------------------------------------------------------

GP_pred <- function(X_train, Y_train, X_pred, theta1 = 1, theta2 = 0.4, theta3 = 0, type = c("rbf", "linear", "quad", "rbf_linear")) {

  X_train <- as.matrix(X_train); Y_train <- as.matrix(Y_train); X_pred <- as.matrix(X_pred)

  N <- nrow(X_train); M <- nrow(X_pred)

  if (N != length(Y_train)) { stop("length(x_train) must equal length(y_train)") }

  K <- kernel_matrix(X_train, X_train, type = type, theta1 = theta1, theta2 = theta2, noise = theta3)

  k <- kernel_matrix(X_train, X_pred,  type = type, theta1 = theta1, theta2 = theta2)

  s <- kernel_matrix(X_pred,  X_pred,  type = type, theta1 = theta1, theta2 = theta2)

  L <- t(chol(K)); yy <- chol_solve(L, Y_train); v <- chol_solve(L, k)

  mu <- t(k) %*% yy; var_f <- s - t(k) %*% v

  if (abs(theta3) > .Machine$double.eps^0.5) {

    var_y <- var_f + theta3 * diag(M)

    return(list(mu = mu, var_f = var_f, var_y = var_y))

  } else {

    return(list(mu = mu, var_f = var_f))

  }

}

# -----------------------------------------------------------------------

# Standardize処理関数

# -----------------------------------------------------------------------

standardize_fit <- function(X, Z, y, standardize_X = TRUE, intercept_col = 1L) {

  X <- as.matrix(X); Z <- as.matrix(Z); y <- as.numeric(y)

  n <- length(y); p <- ncol(X); d <- ncol(Z)

  y_mean <- mean(y); y_sd <- sd(y); if (y_sd == 0) y_sd <- 1

  y_s <- (y - y_mean) / y_sd

  Z_mean <- colMeans(Z); Z_sd <- apply(Z, 2, sd); Z_sd[Z_sd == 0] <- 1

  Z_s <- sweep(sweep(Z, 2, Z_mean, "-"), 2, Z_sd, "/")

  X_s <- X; X_mean <- rep(0, p); X_sd <- rep(1, p)

  if (standardize_X) {

    for (j in seq_len(p)) {

      if (!is.null(intercept_col) && j == intercept_col) { X_mean[j] <- 0; X_sd[j] <- 1

      } else {

        X_mean[j] <- mean(X[, j]); X_sd[j] <- sd(X[, j]); if (X_sd[j] == 0) X_sd[j] <- 1

        X_s[, j] <- (X[, j] - X_mean[j]) / X_sd[j]

      }

    }

  }

  list(X = X_s, Z = Z_s, y = y_s, y_mean = y_mean, y_sd = y_sd,

       Z_mean = Z_mean, Z_sd = Z_sd, X_mean = X_mean, X_sd = X_sd,

       standardize_X = standardize_X, intercept_col = intercept_col)

}

standardize_apply <- function(Xnew, Znew, std) {

  Xnew <- as.matrix(Xnew); Znew <- as.matrix(Znew)

  Z_s <- sweep(sweep(Znew, 2, std$Z_mean, "-"), 2, std$Z_sd, "/")

  X_s <- Xnew

  if (isTRUE(std$standardize_X)) {

    p <- ncol(Xnew)

    for (j in seq_len(p)) {

      if (!is.null(std$intercept_col) && j == std$intercept_col) next

      X_s[, j] <- (Xnew[, j] - std$X_mean[j]) / std$X_sd[j]

    }

  }

  list(X = X_s, Z = Z_s)

}

destandardize_y    <- function(y_scaled, std) { std$y_mean + std$y_sd * y_scaled }

destandardize_y_sd <- function(sd_scaled, std) { std$y_sd * sd_scaled }
# -----------------------------------------------------------------------

# サンプルデータを使ったGaussian Processの計算

# -----------------------------------------------------------------------

theta1 <- 1.0; theta2 <- 0.4; theta3 <- 0.05

x1 <- c(-1, -0.5, 0, 0.5, 1)

y  <- c(-8.200343, -1.124568, -2.842514, 8.697183, 5.250363)

X  <- cbind(1, x1); Z <- cbind(x1)

std   <- standardize_fit(X, Z, y, standardize_X = TRUE, intercept_col = 1L)

X_s   <- std$X; Z_s <- std$Z; y_s <- std$y

m     <- 100; x1n <- seq(-2, 2, length.out = m)

y_true_line <- 1 + 3*x1n - 0.5*exp(-2*(x1n))

new_s <- standardize_apply(cbind(1, x1n), cbind(x1n), std)

res   <- GP_pred(X_train = Z_s, Y_train = y_s, X_pred = new_s$Z,

                 theta1 = theta1, theta2 = theta2, theta3 = theta3, type = "rbf")

mu_pred <- destandardize_y(res$mu, std)

sd_pred <- destandardize_y_sd(sqrt(diag(res$var_y)), std)

y_lo <- mu_pred - sd_pred; y_hi <- mu_pred + sd_pred

plot(x1n, mu_pred, type = "l", col = "steelblue", lwd = 2.5,

     ylim = c(-20, 20), xlab = "x", ylab = "y", main = "GP Prediction")

polygon(c(x1n, rev(x1n)), c(y_lo, rev(y_hi)), col = rgb(0.27, 0.51, 0.71, 0.2), border = NA)

lines(x1n, y_true_line, col = "black", lwd = 2)

points(x1, y, pch = 16, cex = 1.3)

abline(v = c(-1, 1), col = "orange", lty = 2, lwd = 1.5)

legend("topleft",

       legend = c("GP mean", "+-1 SD", "True f(x)", "Training data"),

       lty = c(1, NA, 1, NA), pch = c(NA, 15, NA, 16),

       col = c("steelblue", rgb(0.27, 0.51, 0.71, 0.3), "black", "black"),

       pt.cex = c(NA, 2, NA, 1.3), bty = "n")

上記のコードによって得られる予測グラフは下図の通りです。

ガウス過程回帰の計算方法については、本章で終わりです。長文を読んでいただきありがとうございました。

6. ガウス過程回帰による外挿の危険性について

せっかくガウス過程回帰を理解できたところではありますが、ひとつ注意すべき点があります。それは外挿に対する予測の不確かさです。

図3を見てみると、オレンジ色の外挿域ではガウス過程回帰の予測はゼロに収束していっています。
これは、ガウス過程の平均が $\mu=k_*^T\cdot \mathbf{yy}$ で計算されるところで、外挿域では $k_*^T\approx 0$ だからです。

ガウス過程回帰は与えらえたデータ付近の局所的な構造にしかフィットしないと言うこともできます。

つまり、人間が物理的な前提知識や、大局的なデータのトレンドを知っていても、予測モデルに反映することができないわけで、ガウス過程回帰は本質的には外挿に弱いモデルです。

ここで、RBFカーネルと線形カーネルを組み合わせたカーネルを使えば、線形項により大域的な傾向をとらえられるのでは？と反論があるかもしれません。実際、図3に対してカーネルを調整すれば、外挿域でもフィットさせることはできます。

しかし、カーネルとして線形項を入れると、基底は明示的に表現していない（カーネルトリック）ので、例えば特定の変数の回帰係数がいくらか、という解釈ができません。

私は、このガウス過程の「外挿に弱い」「解釈が難しい」という課題を解決するのがセミパラメトリックベイズだと考えていて、その解説は別の記事にまとめましたので、根気のある方は本記事と合わせて読んでいただけると嬉しいです。
www.doe-get-started.com

付録

条件付き正規分布の公式

$\begin{pmatrix}a\\b\end{pmatrix}\sim N\!\left(\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix},\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\right)$