【python】ガウス過程回帰の仕組みと実務での応用|ノンパラメトリック機械学習
こんにちは、青の統計学です。
今回はガウス過程回帰について解説いたします。
製造業の現場など、n=20やそこらぐらいのデータセットで予測を行う必要がある時によく使われます。
ガウス過程は少数データとの相性がよく、予測値の分散まで出せるのがメリットです。
また、ガウス過程モデルはハイパーパラメータチューニングの一つ「ベイズ最適化」に応用されます。
【kaggle】ベイズ最適化とXGBでtitanicの予測問題を解く|python
ガウス過程回帰(gaussian processing regression)
観測データに基づいて連続的な関数を推定するためのベイズ推定手法です。
ガウス過程回帰の基本的な考え方は、観測データが多変量正規分布（Multivariate Normal Distribution）から生成されるという仮定に基づいています。
ベイズ線型回帰とは異なり、基底関数を明示的に与えなくて良いというメリットがあり、関数の分布を直接仮定するという特徴があります。
ガウス過程とは(gaussian processing)
ガウス過程回帰では、ガウス過程は関数の不確実性をモデル化するために使用されます。
まずはガウス過程とは何か理解する必要があります。
まず入力\(\boldsymbol{x}_1,\boldsymbol{x}_2 ,…,\boldsymbol{x}_N \in X \)に対する出力ベクトルを考えます。
当然以下のようなN次元ベクトル空間になりますね。
$$\boldsymbol{f}=(f(\boldsymbol{x}_1 ),f(\boldsymbol{x}_2 ),..,f(\boldsymbol{x}_N )) $$
この同時分布がガウス分布に従うときに、fはガウス過程に従うと言います。
$$f〜N(m,\boldsymbol{K} )$$
ガウス分布についてまずおさらいしたいという方は以下のコンテンツをご覧ください。
平均関数\(m(\boldsymbol{x} )\)について
平均関数は\(m(x)\) は、ガウス過程における任意の点 \(x\)における関数の期待値を表します。
つまり、関数の出力がどの程度の値になるかの「平均的な」予測を提供します。
多くの場合、平均関数はゼロ関数（すなわち、すべての x に対して\(m(x)=0\)）として設定されます。
これは、ガウス過程が主に関数の形状をモデル化するため、関数が0の周辺にあることを仮定しています。
カーネル関数\(k(\boldsymbol{x},\boldsymbol{x}’)\)について
\(K\)はカーネル関数で定義される\(n×n\)のグラム行列です。
カーネル関数は、カーネルSVMやニューラルネットワークなどにも使用される出力の共分散(言い換えれば類似度)を表現する関数です。
一つカーネル関数をご紹介します。
-EQカーネル(exponential squared kernel)-
$$k(x, x’) = \exp\left(-\frac{||x – x’||^2}{2l^2}\right)$$
EQカーネルは、ガウス過程回帰で最も一般的に使用されるカーネルの一つであり、ガウスカーネルや指数二次元カーネルと呼ばれたりします。
\(l\)はスケールパラメータ（または長さスケール）です。このパラメータは、関数の滑らかさを制御します。\(l\) が大きいほど、関数は滑らかになります。
こちらのコンテンツでもカーネルに関して取り上げております。
【python】カーネルSVMとは？kernel関数を利用した非線形データの判別問題に挑戦|機械学習
ちなみにカーネル関数が基底ベクトルの内積で与えられる場合、ガウス過程回帰の予測分布はベイズ線型回帰の予測分布と一致します。
カーネル関数のスケールパラメータ\(l\)と不確実性を表す振幅(下画像の青線内)が、ガウス過程回帰のハイパーパラメータになります。
このハイパーパラメータのチューニングは、周辺尤度を最大化する過程で探索します。
平均関数とカーネル関数についてご紹介しました。
カーネルに関してはデータ点間の類似度で共分散を定義し、これを回帰に使うことになります。
ガウス過程は、簡単にいえば「無限次元のガウス分布」として理解することができます。
つまり、それは関数の集合であり、それぞれの関数はある平均と共分散構造を持つということです。
では、実際に回帰にガウス過程を使う意義や仕組みに触れていきます。
数学的背景
ガウス過程回帰は、どんな入力に対しても同時分布がガウス分布になるという仮定をおきます。
つまり、教師データとして与えられた\(f(\boldsymbol{x}_1),…,f(\boldsymbol{x}_N)\)と未知のデータ\(f(\boldsymbol{x}_1),…,f(\boldsymbol{x}_T)\)の同時分布を考えます。
上記の同時分布は、ガウス過程に従うため以下のような関係が数式として出せます。
$$p(\boldsymbol{f}_N,\boldsymbol{f}_T)〜N(\begin{bmatrix}\boldsymbol{m}_N\\\boldsymbol{m}_T\\ \end{bmatrix},\begin{bmatrix}\boldsymbol{K}_N&& \boldsymbol{K}_{NT} \\ \boldsymbol{K}_{TN}&& \boldsymbol{K}_T\\ \end{bmatrix})$$
Kはグラム行列ですね。
尤度に関しては、出力が分散\(σ_{v}^2\)のガウスノイズを持つとみなせるし、\(f(x)\)と出力\(y\)の差分がガウス分布に従うことを考えると、
$$p(\boldsymbol{y}_N,\boldsymbol{y}_T|\boldsymbol{f}_N,\boldsymbol{f}_T)〜N(\begin{bmatrix}\boldsymbol{f}_N\\ \boldsymbol{f}_T\\ \end{bmatrix},\begin{bmatrix}σ_{v}^2 \boldsymbol{I}_N&& \boldsymbol{0} \\ \boldsymbol{0}&& σ_{v}^2 \boldsymbol{I}_T\\ \end{bmatrix})$$
Iは単位行列ですね。
上の二つの式から、同時周辺分布は以下のように表すことができます。
$$p(\boldsymbol{y}_N,\boldsymbol{y}_T)〜N(\begin{bmatrix}\boldsymbol{m}_N\\ \boldsymbol{m}_T\\ \end{bmatrix},\begin{bmatrix} \boldsymbol{K}_N+σ_{v}^2 \boldsymbol{I}_N&& \boldsymbol{K}_{NT} \\ \boldsymbol{K}_{TN}&& \boldsymbol{K}_{T}+σ_{v}^2 \boldsymbol{I}_T\\ \end{bmatrix})$$
ここで、肝心の予測分布を求めるために多変量ガウス分布の条件付き分布の公式を使います。
$$\begin{align} \mu_{X|Y} &= \mu_X + \Sigma_{XY} \Sigma_{YY}^{-1} (y – \mu_Y) \ \Sigma_{X|Y} &= \Sigma_{XX} – \Sigma_{XY} \Sigma_{YY}^{-1} \Sigma_{YX} \end{align}$$
この公式に尤度と同時周辺分布を与えてみます。
$$p(\boldsymbol{y}_T|\boldsymbol{y}_N))=N(\boldsymbol{μ},\boldsymbol{Σ})$$
$$\boldsymbol{μ}=\boldsymbol{m}_T+\boldsymbol{K}_{TN}(\boldsymbol{K}_N+σ_{v}^2\boldsymbol{I}_{N})^{-1}(\boldsymbol{y}_{N}-\boldsymbol{m}_{N})$$
$$\boldsymbol{Σ}=\boldsymbol{K}_T+σ_{v}^2\boldsymbol{I}_{N}\boldsymbol{K}_{TN}(\boldsymbol{K}_N+σ_{v}^2\boldsymbol{I}_{N})^{-1}\boldsymbol{K}_{NT})$$
普通は平均関数が0なのでもう少し条件付き平均はシンプルになります。
CODE|python
# 必要なライブラリをインポート
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
from sklearn.metrics import mean_squared_error
import numpy as np
# 合成データを生成
def f(x):
return x * np.sin(x)
X = np.atleast_2d(np.linspace(0, 10, 100)).T
y = f(X).ravel()
# データを訓練データとテストデータに分割
X_train = X[::5]
y_train = y[::5]
X_test = np.concatenate((X[1::5], X[2::5], X[3::5], X[4::5]))
y_test = np.concatenate((y[1::5], y[2::5], y[3::5], y[4::5]))
# ガウス過程回帰モデルを定義
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
# モデルを訓練データに適合
gp.fit(X_train, y_train)
# テストデータに対する予測を行う
y_pred = gp.predict(X_test)
# 予測の精度を計算
mse = mean_squared_error(y_test, y_pred)
print('Mean Squared Error: ', mse)
Mean Squared Error: 1.1393927523823373e-08
1次元の合成データに対するガウス過程回帰を行い、テストデータに対する予測の精度（平均二乗誤差）を計算します。
基本は多次元での回帰になるため、予測曲線とその振幅の描画は難しいです。
