ガウス過程分類(GPC)をわかりやすく解説

ガウス過程回帰(GPR)は、関数空間上にガウス過程の事前分布を置くことで、柔軟なノンパラメトリック回帰を実現する手法でした。では、この強力な枠組みを分類問題に拡張するにはどうすればよいでしょうか。

本記事では、ガウス過程を分類問題へと拡張したガウス過程分類(Gaussian Process Classification, GPC)について解説します。GPRとの違い、リンク関数の役割、ラプラス近似による事後分布の導出、そして scikit-learn を用いた Python 実装までを一通り扱います。GPCは予測の不確実性を自然に定量化できるベイズ的な分類手法であり、医療診断や異常検知のように「予測の確信度」が重要な場面で特に有用です。

本記事の内容

  • ガウス過程回帰から分類への拡張の考え方
  • シグモイド・プロビットリンク関数による尤度の構成
  • ラプラス近似による事後分布の導出(ヘッセ行列、ニュートン法)
  • 予測分布の計算方法
  • GPR と GPC の対比
  • scikit-learn による2次元分類問題の実装と決定境界の可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

ガウス過程の基本概念(カーネル関数、カーネル行列 $\bm{K}$、予測分布の導出など)を前提とします。

分類問題の定式化 — 回帰との違い

ガウス過程回帰の復習

ガウス過程回帰では、目的変数 $y$ は連続値であり、潜在関数 $f(\bm{x})$ にガウスノイズを加えたモデルを考えました。

$$ y = f(\bm{x}) + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma_n^2) $$

このとき、尤度関数 $p(y \mid f)$ はガウス分布となります。ガウス過程の事前分布もガウス分布であるため、事後分布 $p(\bm{f} \mid \bm{y})$ もまたガウス分布として解析的に求まるという嬉しい性質がありました。

分類問題への拡張

一方、分類問題では目的変数 $y$ は離散値(2値分類なら $y \in \{0, 1\}$)です。離散的なラベルに対してガウス分布の尤度をそのまま使うことはできません。

GPCでは、次のような2段階の構造を導入します。

  1. ガウス過程の事前分布に従う潜在関数 $f(\bm{x})$ を導入する
  2. 潜在関数 $f(\bm{x})$ をリンク関数で $[0, 1]$ の確率に変換する

つまり、クラス1に属する確率を次のようにモデル化します。

$$ p(y = 1 \mid \bm{x}) = \sigma(f(\bm{x})) $$

ここで $\sigma(\cdot)$ はリンク関数です。潜在関数 $f$ はガウス過程に従います。

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

$k(\bm{x}, \bm{x}’)$ はカーネル関数で、典型的にはRBF(ガウス)カーネルを使います。

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

$N$ 個の訓練データ $\bm{X} = \{\bm{x}_1, \dots, \bm{x}_N\}$ に対する潜在関数の値 $\bm{f} = (f_1, \dots, f_N)^T$ は、多変量ガウス分布に従います。

$$ p(\bm{f} \mid \bm{X}) = \mathcal{N}(\bm{f} \mid \bm{0}, \bm{K}) $$

ここで $\bm{K}$ はカーネル行列であり、$K_{ij} = k(\bm{x}_i, \bm{x}_j)$ です。

シグモイド(ロジスティック)関数

最もよく使われるリンク関数はロジスティックシグモイド関数です。

$$ \sigma(f) = \frac{1}{1 + \exp(-f)} $$

この関数は $f \in (-\infty, +\infty)$ を $(0, 1)$ に写像し、$f = 0$ のとき $\sigma(f) = 0.5$ となります。ロジスティック回帰と同じ変換であり、直感的にわかりやすい選択です。

プロビット関数

もう一つの代表的なリンク関数が、標準正規分布の累積分布関数(CDF)を用いたプロビット関数です。

$$ \Phi(f) = \int_{-\infty}^{f} \mathcal{N}(z \mid 0, 1) \, dz $$

プロビット関数はシグモイド関数と形状がよく似ていますが、裾の減衰がやや速いという違いがあります。プロビット関数を使う利点は、ガウス分布との親和性が高いことです。後述する予測分布の計算で、ガウス分布とプロビット関数の畳み込みが解析的に計算できるため、近似公式が簡潔になります。

尤度関数

2値分類($y \in \{0, 1\}$)の尤度関数は、リンク関数 $\sigma$ を用いて次のように書けます。

$$ p(y \mid f) = \sigma(f)^{y} \, (1 – \sigma(f))^{1 – y} $$

これはベルヌーイ分布そのものです。$y = 1$ のとき $p(y \mid f) = \sigma(f)$、$y = 0$ のとき $p(y \mid f) = 1 – \sigma(f)$ となります。

$N$ 個のデータが独立に観測されたとすると、尤度全体は次のようになります。

$$ p(\bm{y} \mid \bm{f}) = \prod_{i=1}^{N} \sigma(f_i)^{y_i} (1 – \sigma(f_i))^{1 – y_i} $$

GPCの事後分布 — なぜ近似が必要か

GPCでは、事後分布はベイズの定理より次のように書けます。

$$ p(\bm{f} \mid \bm{X}, \bm{y}) = \frac{p(\bm{y} \mid \bm{f}) \, p(\bm{f} \mid \bm{X})}{p(\bm{y} \mid \bm{X})} $$

ここで、

  • $p(\bm{f} \mid \bm{X}) = \mathcal{N}(\bm{f} \mid \bm{0}, \bm{K})$ はガウス過程の事前分布(ガウス分布)
  • $p(\bm{y} \mid \bm{f}) = \prod_{i=1}^{N} \sigma(f_i)^{y_i}(1 – \sigma(f_i))^{1 – y_i}$ はベルヌーイ尤度(非ガウス)

ガウス過程回帰では、事前分布(ガウス分布)と尤度(ガウス分布)の積がガウス分布となるため、事後分布を解析的に計算できました。しかし、GPCでは尤度がベルヌーイ分布(シグモイド関数を含む)であるため、事後分布 $p(\bm{f} \mid \bm{X}, \bm{y})$ はガウス分布にはなりません。したがって、何らかの近似手法が必要です。

代表的な近似手法を以下にまとめます。

近似手法 特徴
ラプラス近似 事後分布のモード(MAP推定値)の周りで2次のテイラー展開を行い、ガウス分布で近似する
期待値伝搬法(EP) 各尤度因子を順番にガウス分布で近似していく反復アルゴリズム
変分推論(VI) 変分下界を最大化して近似分布を求める

本記事では、最も基本的かつ理解しやすいラプラス近似を取り上げます。

ラプラス近似

ラプラス近似の基本的なアイデアは、事後分布のモード(MAP推定値)を中心として、対数事後分布を2次のテイラー展開し、ガウス分布で近似するというものです。

非正規化対数事後分布

事後分布の対数(正規化定数を除いたもの)を $\Psi(\bm{f})$ とおきます。

$$ \Psi(\bm{f}) = \log p(\bm{y} \mid \bm{f}) + \log p(\bm{f} \mid \bm{X}) $$

各項を具体的に書き下します。まず、対数尤度は次のようになります。

$$ \log p(\bm{y} \mid \bm{f}) = \sum_{i=1}^{N} \left[ y_i \log \sigma(f_i) + (1 – y_i) \log(1 – \sigma(f_i)) \right] $$

対数事前分布は次のとおりです。

$$ \log p(\bm{f} \mid \bm{X}) = -\frac{1}{2} \bm{f}^T \bm{K}^{-1} \bm{f} – \frac{N}{2} \log(2\pi) – \frac{1}{2} \log |\bm{K}| $$

定数項を除けば、非正規化対数事後分布は次のようになります。

$$ \Psi(\bm{f}) = \sum_{i=1}^{N} \left[ y_i \log \sigma(f_i) + (1 – y_i) \log(1 – \sigma(f_i)) \right] – \frac{1}{2} \bm{f}^T \bm{K}^{-1} \bm{f} + \text{const.} $$

MAPの解を求める — 勾配の計算

ラプラス近似では、まず $\Psi(\bm{f})$ を最大化する $\hat{\bm{f}}$(MAP推定値)を求めます。そのために勾配を計算しましょう。

まず、シグモイド関数の微分の性質を確認しておきます。

$$ \frac{d\sigma(f)}{df} = \sigma(f)(1 – \sigma(f)) $$

この性質を利用して、対数尤度の $f_i$ に関する微分を丁寧に計算します。

$$ \begin{align} \frac{\partial \log p(y_i \mid f_i)}{\partial f_i} &= \frac{\partial}{\partial f_i} \left[ y_i \log \sigma(f_i) + (1 – y_i) \log(1 – \sigma(f_i)) \right] \\ &= y_i \frac{\sigma(f_i)(1 – \sigma(f_i))}{\sigma(f_i)} + (1 – y_i) \frac{-\sigma(f_i)(1 – \sigma(f_i))}{1 – \sigma(f_i)} \\ &= y_i (1 – \sigma(f_i)) – (1 – y_i) \sigma(f_i) \\ &= y_i – y_i \sigma(f_i) – \sigma(f_i) + y_i \sigma(f_i) \\ &= y_i – \sigma(f_i) \end{align} $$

したがって、対数尤度の勾配ベクトルは次のようになります。

$$ \nabla_{\bm{f}} \log p(\bm{y} \mid \bm{f}) = \bm{y} – \bm{\pi} $$

ここで $\bm{\pi} = (\sigma(f_1), \sigma(f_2), \dots, \sigma(f_N))^T$ です。

事前分布の項の勾配は次のとおりです。

$$ \nabla_{\bm{f}} \log p(\bm{f} \mid \bm{X}) = -\bm{K}^{-1} \bm{f} $$

これらを合わせると、$\Psi(\bm{f})$ 全体の勾配は次のようになります。

$$ \nabla_{\bm{f}} \Psi(\bm{f}) = \bm{y} – \bm{\pi} – \bm{K}^{-1} \bm{f} $$

勾配を $\bm{0}$ とおくと、MAP条件が得られます。

$$ \hat{\bm{f}} = \bm{K}(\bm{y} – \hat{\bm{\pi}}) $$

ここで $\hat{\bm{\pi}} = \sigma(\hat{\bm{f}})$ です。この式は $\hat{\bm{f}}$ について非線形であるため、解析的には解けません。

ヘッセ行列の導出

ニュートン法にはヘッセ行列(2階微分)が必要です。$\Psi(\bm{f})$ の2階微分を計算しましょう。

対数尤度のヘッセ行列の各成分は次のようになります。

$$ \begin{align} \frac{\partial^2 \log p(y_i \mid f_i)}{\partial f_i^2} &= \frac{\partial}{\partial f_i} (y_i – \sigma(f_i)) \\ &= -\sigma(f_i)(1 – \sigma(f_i)) \\ &= -\pi_i(1 – \pi_i) \end{align} $$

異なる成分間の交差微分は $0$ なので、対数尤度のヘッセ行列は対角行列になります。対角行列 $\bm{W}$ を次のように定義します。

$$ \bm{W} = \text{diag}(\pi_1(1 – \pi_1), \, \pi_2(1 – \pi_2), \, \dots, \, \pi_N(1 – \pi_N)) $$

$\bm{W}$ の各対角成分はシグモイド関数の分散に相当し、常に正の値をとります。すると、$\Psi(\bm{f})$ 全体のヘッセ行列は次のように書けます。

$$ \nabla_{\bm{f}}^2 \Psi(\bm{f}) = -\bm{W} – \bm{K}^{-1} $$

$\bm{W}$ も $\bm{K}^{-1}$ も正定値行列であるため、ヘッセ行列は負定値です。これは $\Psi(\bm{f})$ が凹関数であることを意味し、ニュートン法がモードに確実に収束することを保証します。

ニュートン法の更新式

ニュートン法の一般的な更新式は次のとおりです。

$$ \bm{f}^{\text{new}} = \bm{f}^{\text{old}} – (\nabla_{\bm{f}}^2 \Psi)^{-1} \nabla_{\bm{f}} \Psi $$

勾配とヘッセ行列を代入すると、次の更新式が得られます。

$$ \begin{align} \bm{f}^{\text{new}} &= \bm{f}^{\text{old}} + (\bm{W} + \bm{K}^{-1})^{-1} \left[ (\bm{y} – \bm{\pi}) – \bm{K}^{-1} \bm{f}^{\text{old}} \right] \end{align} $$

これを整理しましょう。$(\bm{W} + \bm{K}^{-1})^{-1}$ でまとめると、

$$ \begin{align} \bm{f}^{\text{new}} &= (\bm{W} + \bm{K}^{-1})^{-1} \left[ (\bm{W} + \bm{K}^{-1}) \bm{f}^{\text{old}} + \bm{y} – \bm{\pi} – \bm{K}^{-1} \bm{f}^{\text{old}} \right] \\ &= (\bm{K}^{-1} + \bm{W})^{-1} (\bm{W} \bm{f}^{\text{old}} + \bm{y} – \bm{\pi}) \end{align} $$

擬似目的変数 $\bm{b} = \bm{W} \bm{f}^{\text{old}} + \bm{y} – \bm{\pi}$ を定義すれば、更新式は次のように簡潔に書けます。

$$ \bm{f}^{\text{new}} = (\bm{K}^{-1} + \bm{W})^{-1} \bm{b} $$

この反復を $\bm{f}$ が収束するまで繰り返すことで、MAP推定値 $\hat{\bm{f}}$ が得られます。

ガウス近似の事後分布

MAP推定値 $\hat{\bm{f}}$ が求まったら、$\Psi(\bm{f})$ を $\hat{\bm{f}}$ の周りで2次のテイラー展開します。$\hat{\bm{f}}$ はモードであるため1次の項は $0$ になり、

$$ \Psi(\bm{f}) \approx \Psi(\hat{\bm{f}}) + \frac{1}{2} (\bm{f} – \hat{\bm{f}})^T \nabla_{\bm{f}}^2 \Psi(\hat{\bm{f}}) \, (\bm{f} – \hat{\bm{f}}) $$

これは対数ガウス分布の形をしています。したがって、近似事後分布は次のガウス分布となります。

$$ q(\bm{f} \mid \bm{X}, \bm{y}) = \mathcal{N}(\bm{f} \mid \hat{\bm{f}}, \, (\bm{K}^{-1} + \bm{W})^{-1}) $$

つまり、ラプラス近似の結果、事後分布は次のパラメータを持つガウス分布で近似されます。

$$ \begin{align} \text{平均:} \quad & \hat{\bm{f}} = \bm{K}(\bm{y} – \hat{\bm{\pi}}) \\ \text{共分散:} \quad & \bm{\Sigma} = (\bm{K}^{-1} + \bm{W})^{-1} \end{align} $$

予測分布の計算

新しいテスト入力 $\bm{x}_*$ に対する予測は、2つのステップで行います。まず潜在関数 $f_*$ の予測分布を求め、次にリンク関数を通してクラス確率に変換します。

潜在関数の予測分布

テスト入力 $\bm{x}_*$ に対する潜在関数 $f_*$ の予測分布は、近似事後分布を用いて次のように計算されます。

$$ q(f_* \mid \bm{X}, \bm{y}, \bm{x}_*) = \int p(f_* \mid \bm{X}, \bm{x}_*, \bm{f}) \, q(\bm{f} \mid \bm{X}, \bm{y}) \, d\bm{f} $$

ガウス過程の性質と近似事後分布がどちらもガウス分布であることから、この積分は解析的に計算でき、結果は次のガウス分布になります。

$$ q(f_* \mid \bm{X}, \bm{y}, \bm{x}_*) = \mathcal{N}(f_* \mid \mu_{f_*}, \sigma_{f_*}^2) $$

ここで、平均と分散は次のとおりです。

$$ \mu_{f_*} = \bm{k}_*^T (\bm{y} – \hat{\bm{\pi}}) $$

$$ \sigma_{f_*}^2 = k(\bm{x}_*, \bm{x}_*) – \bm{k}_*^T (\bm{K} + \bm{W}^{-1})^{-1} \bm{k}_* $$

$\bm{k}_* = (k(\bm{x}_*, \bm{x}_1), \dots, k(\bm{x}_*, \bm{x}_N))^T$ は、訓練データとテスト入力間のカーネルベクトルです。

平均 $\mu_{f_*}$ の式を GPR の予測平均 $\bm{k}_*^T \bm{K}^{-1} \bm{y}$ と比較すると、構造が似ていることがわかります。回帰では重みが $\bm{K}^{-1} \bm{y}$ であるのに対し、分類では $\bm{y} – \hat{\bm{\pi}}$ に変わっています。分散の式も GPR と似ていますが、$\bm{W}^{-1}$ が加わっている点が異なります。

クラス確率への変換

最終的に求めたいのは、テスト入力 $\bm{x}_*$ がクラス1に属する確率です。潜在関数の予測分布をリンク関数で変換して得られます。

$$ \bar{\pi}_* = p(y_* = 1 \mid \bm{X}, \bm{y}, \bm{x}_*) = \int \sigma(f_*) \, q(f_* \mid \bm{X}, \bm{y}, \bm{x}_*) \, df_* $$

この積分は、シグモイド関数とガウス分布の畳み込みであり、解析的には計算できません。しかし、プロビット近似を用いると、次の近似式が得られます。

$$ \bar{\pi}_* \approx \sigma\!\left( \frac{\mu_{f_*}}{\sqrt{1 + \frac{\pi}{8} \sigma_{f_*}^2}} \right) $$

ここで分母の $\frac{\pi}{8}$ は、シグモイド関数をプロビット関数(ガウス分布のCDF)で近似したときに現れるスケーリング定数です。

この近似式は、GPCの予測が持つ重要な性質を端的に示しています。

  • $\mu_{f_*}$ が大きいほど $\bar{\pi}_*$ は1に近づく(確信度の高い予測)
  • $\sigma_{f_*}^2$ が大きい(不確実性が高い)ほど、分母が大きくなり $\bar{\pi}_*$ は0.5に近づく(「わからない」という予測)

つまり、GPCはデータが少ない領域では自然に「自信のない」予測を返すという、ベイズ的な手法ならではの望ましい性質を持っています。

GPR と GPC の対比

ガウス過程回帰と分類の違いを表にまとめます。

項目 GPR(回帰) GPC(分類)
目的変数 連続値 $y \in \mathbb{R}$ 離散値 $y \in \{0, 1\}$
尤度関数 ガウス分布 $\mathcal{N}(y \mid f, \sigma_n^2)$ ベルヌーイ分布 $\text{Bern}(y \mid \sigma(f))$
事後分布 解析的にガウス分布 非ガウス(近似が必要)
近似手法 不要 ラプラス近似 / EP / VI
予測の出力 平均と分散(連続値) クラス確率 $p(y=1 \mid \bm{x})$
潜在関数の意味 直接的に目的変数の予測値 リンク関数を通す前の対数オッズに対応
予測平均の重み $\bm{K}^{-1} \bm{y}$ $\bm{y} – \hat{\bm{\pi}}$
予測分散の補正 $\bm{k}_*^T \bm{K}^{-1} \bm{k}_*$ を引く $\bm{k}_*^T (\bm{K} + \bm{W}^{-1})^{-1} \bm{k}_*$ を引く
計算コスト $\mathcal{O}(N^3)$ $\mathcal{O}(N^3 \times \text{反復回数})$
カーネル関数 RBF, Matern 等 共通(RBF, Matern 等)

Pythonでの実装

scikit-learn の GaussianProcessClassifier を使って、2次元の分類問題を解き、決定境界を可視化してみましょう。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
from sklearn.datasets import make_moons

# --- データ生成 ---
# 半月形の2クラス分類データ
X_train, y_train = make_moons(n_samples=100, noise=0.3, random_state=42)

# --- カーネルの定義 ---
# ConstantKernel * RBFカーネル
kernel = C(1.0, (1e-3, 1e3)) * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e2))

# --- GPCモデルの学習 ---
gpc = GaussianProcessClassifier(kernel=kernel, n_restarts_optimizer=10, random_state=42)
gpc.fit(X_train, y_train)

print(f"学習後のカーネル: {gpc.kernel_}")
print(f"対数周辺尤度: {gpc.log_marginal_likelihood_value_:.3f}")

# --- 決定境界の可視化 ---
# 予測用のグリッド作成
x_min, x_max = X_train[:, 0].min() - 1, X_train[:, 0].max() + 1
y_min, y_max = X_train[:, 1].min() - 1, X_train[:, 1].max() + 1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, 200),
                      np.linspace(y_min, y_max, 200))
X_grid = np.c_[xx.ravel(), yy.ravel()]

# 各グリッド点でのクラス1の予測確率
Z_prob = gpc.predict_proba(X_grid)[:, 1].reshape(xx.shape)

# プロット
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左図: 決定境界 + クラス確率
ax = axes[0]
contour = ax.contourf(xx, yy, Z_prob, levels=np.linspace(0, 1, 21),
                       cmap="RdBu_r", alpha=0.8)
ax.contour(xx, yy, Z_prob, levels=[0.5], colors="k", linewidths=2)
ax.scatter(X_train[y_train == 0, 0], X_train[y_train == 0, 1],
           c="blue", edgecolors="k", label="Class 0", s=40)
ax.scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1],
           c="red", edgecolors="k", label="Class 1", s=40)
ax.set_xlabel("$x_1$")
ax.set_ylabel("$x_2$")
ax.set_title("GPC: Decision Boundary and Class Probability")
ax.legend(loc="upper left")
fig.colorbar(contour, ax=ax, label="$p(y=1 | x)$")

# 右図: 予測の不確実性(エントロピー)
entropy = -Z_prob * np.log2(np.clip(Z_prob, 1e-10, 1)) \
          - (1 - Z_prob) * np.log2(np.clip(1 - Z_prob, 1e-10, 1))
ax2 = axes[1]
contour2 = ax2.contourf(xx, yy, entropy, levels=20, cmap="YlOrRd", alpha=0.8)
ax2.contour(xx, yy, Z_prob, levels=[0.5], colors="k", linewidths=2)
ax2.scatter(X_train[y_train == 0, 0], X_train[y_train == 0, 1],
            c="blue", edgecolors="k", label="Class 0", s=40)
ax2.scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1],
            c="red", edgecolors="k", label="Class 1", s=40)
ax2.set_xlabel("$x_1$")
ax2.set_ylabel("$x_2$")
ax2.set_title("Prediction Uncertainty (Entropy)")
ax2.legend(loc="upper left")
fig.colorbar(contour2, ax=ax2, label="Entropy (bits)")

plt.tight_layout()
plt.show()

このコードを実行すると、2つの図が得られます。

左図は、GPCが学習した決定境界とクラス確率を示しています。色は予測確率 $p(y=1 \mid \bm{x})$ を表し、黒線が決定境界($p = 0.5$)に対応します。GPCではSVMのような硬い境界ではなく、確率的に滑らかな境界が得られることがわかります。

右図は、予測のエントロピー(不確実性)を可視化しています。決定境界付近やデータが存在しない領域でエントロピーが高く(不確実性が大きく)なっていることが確認できます。これはGPCの大きな利点であり、「この予測にどれだけ自信があるか」を定量的に評価できます。

まとめ

本記事では、ガウス過程分類(GPC)について、理論の導出から Python での実装まで解説しました。

  • ガウス過程回帰から分類への拡張では、リンク関数(シグモイド/プロビット)を導入して、潜在関数の出力をクラス確率に変換する
  • 尤度関数がベルヌーイ分布(非ガウス)になるため、事後分布 $p(\bm{f} \mid \bm{X}, \bm{y})$ は解析的に計算できず、近似手法が必要
  • ラプラス近似では、非正規化対数事後分布 $\Psi(\bm{f})$ のMAP推定値 $\hat{\bm{f}}$ をニュートン法で求め、ヘッセ行列 $-\bm{W} – \bm{K}^{-1}$ を用いてガウス分布で近似する
  • 予測分布は、潜在関数の予測分布(ガウス分布)をリンク関数で変換して得られ、プロビット近似により効率的に計算できる
  • GPCは分類の不確実性を定量化できる点が大きな利点であり、データが少ない領域では予測が自然に $0.5$ に近づく

次のステップとして、ガウス過程の重要な応用であるベイズ最適化について学ぶことをおすすめします。ベイズ最適化では、ガウス過程をブラックボックス関数のサロゲートモデルとして用い、獲得関数を最大化することで効率的に最適解を探索します。