【初心者向け】ガウス過程回帰をわかりやすく説明〜Pythonで実装してみた〜

こんにちは！ぼりたそです！今回はベイズ最適化でも使用されるガウス過程回帰について初学者の方に向けてわかりやすく説明した記事になります。

この記事は以下の点をまとめています

Point

ガウス過程回帰とは？
　
ガウス過程回帰のメリット、デメリット
　
ガウス過程回帰のプロセス
　
Pythonで実行
　
オススメの書籍

それではガウス過程回帰についてご説明していきます。

ガウス過程回帰とは？
ガウス過程回帰のメリットとデメリット
ガウス過程回帰のプロセス
ガウス過程回帰をPythonで実行
オススメの書籍
終わりに

ガウス過程回帰とは？

ガウス過程回帰（Gaussian Process Regression）は、機械学習の一種であり、回帰に使用される手法です。

ガウス過程回帰の考え方として、データがガウス分布に従うと仮定しています。つまり、データの値は、確率的な変数として考えることができます。

この仮定に基づいて、与えられた入力データに対して、未知の出力値（予測値）とその予測の不確かさ（分散）を推定することができます。

この予測値の推測と同時に分散も推測できるのがガウス過程回帰の最大の特徴になります。ひとまずはこれだけ覚えていただければ問題ないと思います。

とはいえ、これだけだとあまりイメージが湧かないと思いますので、実際にガウス過程回帰をしたグラフを下に示します。

このグラフはy＝sin(x)に従う5点を入力データとしてガウス過程回帰を実行した結果になります。緑の破線が真のグラフであるのに対して、青線が予測のグラフになっています。

さらに、青色で塗りつぶされているのが予測の分散になります。真のグラフに対してデータが少ない領域は分散が大きく、予測値も外れていることがわかりますね。

このように入力データに対して、予測値とその分散を出力できるのがガウス過程回帰になります。

ガウス過程回帰のメリットとデメリット

次にガウス過程回帰のメリットとデメリットについてご説明します。

まずメリットは以下の通りです。

■メリット

不確実性の推定が可能:
ガウス過程回帰は、予測の不確実性を推定することができます。予測値のみならず、予測の信頼性や予測の範囲（信頼区間）も提供するため、モデルの予測結果がどれだけ信頼性があるかを理解することができます。
非線形な関数モデリングが可能:
ガウス過程回帰は非線形な関数のモデリングに適しています。カーネル関数を適切に選択することで、入力データの非線形な関係性を柔軟にモデル化することができます。
パラメータの最適化が容易:
ガウス過程回帰では、ハイパーパラメータの最適化が比較的容易です。最尤推定や交差検証などを用いて、カーネル関数のハイパーパラメータを最適な値に調整することができます。

次にデメリットについて以下に示します。

■デメリット

計算コスト:
ガウス過程回帰は、カーネル行列の計算や逆行列の計算などの計算コストが高いため、大規模なデータセットには適していません。データ数が増えると、計算時間が急激に増加することがあります。
データの依存性:
ガウス過程回帰は、データに強く依存するモデルです。入力データの分布やパターンに依存してモデルの性能が変わるため、十分なデータがない場合や入力データに偏りがある場合には適切な予測が困難になることがあります。

ガウス過程回帰は柔軟で強力なモデルであり、予測の不確実性を考慮に入れることができるため、多くの場面で利用されています。ただし、データの規模やモデルの選択に注意しながら利用する必要があると言えるでしょう。

ガウス過程回帰のプロセス

次にガウス過程回帰のプロセスについてご説明します。

ガウス過程回帰のプロセスとしては大きく5つに分けることができます。
具体的には以下の通りです。

入力データの準備
　
カーネル関数の選択
　
カーネル行列の計算
　
ハイパーパラメータの最適化
　
予測の実行

それでは詳細に解説していきます。

入力データの準備

まず、回帰分析をしたいデータセットを用意します。各データポイントには、説明変数とその対応する目的変数が必要になります。

一見当たり前だと思うかもしれませんが、実際には使用できないデータなどが多くあり、入力データを準備するプロセスが一番時間がかかってしまうケースも多いですね。

カーネル関数の選択

ガウス過程回帰では、データポイント間の類似性を計算するためにカーネル関数を選択します。

最も一般的なカーネル関数としては、ガウスカーネル（RBFカーネル）が挙げられます。他にも線形カーネル、多項式カーネル、シグモイドカーネルなどがあります。

選択したカーネル関数によって、モデルの性能が大きく変わることがあるので注意が必要ですね。

カーネルについて以下の記事をまとめていますので、詳細が気になる方は参照いただければと思います。

合わせて読みたい

カーネル行列の計算

選択したカーネル関数を使用して、データセット内のすべてのデータ間の類似性を表すカーネル行列を計算します。

カーネル行列は、データのペアごとにカーネル関数を適用した結果を要素として持つ対称行列です。

カーネル行列について

カーネル行列は、ガウス過程回帰においてデータセット内のすべてのデータポイントのペアに対するカーネル関数の値を格納した行列です。

カーネル行列はガウス過程回帰の学習および予測の際に使用されます。

訓練データの間の類似性を表すカーネル行列を計算し、新しいデータポイントと訓練データの間の類似性を示すカーネルベクトルを使って予測を行います。

カーネル行列を計算する際には、カーネル関数と訓練データの特徴ベクトルを用いて行列の要素を計算します。

カーネルのハイパーパラメータの最適化

ガウス過程回帰にはいくつかのハイパーパラメータ（予め設定しておく値）が存在します。主にカーネル関数のハイパーパラメータを調整する必要があります。

パラメータ最適化の一般的な方法としては、最尤法や交差検証などの方法挙げられます。

ハイパーパラメータ最適化について以下の記事をまとめているので、詳細が気になる方は参照ください。

合わせて読みたい

予測の実行

新しい入力データポイントに対して、予測値と予測の不確かさを推定します。予測の実行には、学習フェーズで得られたカーネル行列とハイパーパラメータが使用されます。

ガウス過程回帰をPythonで実行

それでは実際にPythonを使用してガウス過程回帰を実行してみましょう。

今回はy=sin(x)に従う5つのデータからガウス過程回帰を実行するコードをご紹介します。

このコードではカーネル関数はRBF、ハイパーパラメータの最適化は最尤推定から実行しています。

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

# 真の関数
def true_function(X):
    return np.sin(X).ravel()

# 一元のデータを作成
np.random.seed(1)
X = np.sort(5 * np.random.rand(5, 1), axis=0)
y = true_function(X)

# ガウス過程回帰モデルの構築（最尤推定を使用）
kernel = C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, optimizer=None)

# モデルの学習
gp.fit(X, y)

# 予測のための新しいデータ点を作成
x_pred = np.linspace(0, 5, 100)[:, np.newaxis]

# 予測値と予測の不確かさ（分散）を計算
y_pred, std_dev = gp.predict(x_pred, return_std=True)
lower_bound = y_pred - 1.96 * std_dev
upper_bound = y_pred + 1.96 * std_dev

# プロット
plt.figure(figsize=(10,6))
plt.scatter(X, y, c='r', label='data')
plt.plot(x_pred, y_pred, 'b', label='prediction')
plt.fill_between(x_pred[:, 0], lower_bound, upper_bound, alpha=0.2, color='blue',label='uncertainty')

# 真の関数をプロット
true_y = true_function(x_pred)
plt.plot(x_pred, true_y, 'g', label='true function', linestyle='dashed')

plt.xlabel('X')
plt.ylabel('y')
plt.title('Gaussian Process Regression with Maximum Likelihood Estimation')
plt.legend()
plt.show()

実際に得られたグラフがこちらになります。