ガウス過程回帰モデル

ガウス過程回帰 (GPR) モデルは、ノンパラメトリックなカーネルベースの確率モデルです。GPR モデルは、関数 fitrgp を使用して学習させることができます。

未知の分布から抽出された $$\{(x_i,y_i);i=1,2,\mathrm{...},n\}$$ という学習セットがあるとします。ここで、$$x_i\in {ℝ}^d$$ および $$y_i\in ℝ$$ です。GPR モデルは、与えられた新しい入力ベクトル $$x_{new}$$ と学習データから応答変数 $$y_{new}$$ の値を予測するという問題を扱います。線形回帰モデルは、次のような形式になります。

$$y=x^T\beta +\epsilon ,$$

ここで、$$\epsilon \sim N(0,{\sigma }^2\text{)}$$ です。誤差分散 σ² と係数 β は、データから推定されます。GPR モデルでは、ガウス過程 (GP) による潜在的変数 $$f\left(x_i\right),i=1,2,\mathrm{...},n$$ と明示的な基底関数 h を導入することにより応答を求めます。潜在的変数の共分散関数は応答の滑らかさを捕らえ、基底関数は入力 $$x$$ を p 次元の特徴空間に射影します。

GP は、どの有限個のサブセットも結合ガウス分布になる一連の確率変数です。$$\left\{f\left(x\right),x\in {ℝ}^d\right\}$$ が GP である場合、n 個の観測値 $$x_1,x_2,\mathrm{...},x_n$$ が与えられると、確率変数 $$f(x_1),f(x_2),\mathrm{...},f(x_n)$$ の同時分布はガウス分布になります。GP は、平均関数 $$m\left(x\right)$$ と共分散関数 $$k\left(x,x^{\prime }\right)$$ によって定義されます。つまり、$$\left\{f\left(x\right),x\in {ℝ}^d\right\}$$ がガウス過程である場合、$$E\left(f\left(x\right)\right)=m\left(x\right)$$ および $$Cov\left[f\left(x\right),f\left(x^{\prime }\right)\right]=E\left[\left\{f\left(x\right)-m\left(x\right)\right\}\left\{f\left(x^{\prime }\right)-m\left(x^{\prime }\right)\right\}\right]=k\left(x,x^{\prime }\right).$$ です。

次のようなモデルについて考えます。

$$h{(x)}^T\beta +f(x),$$

ここで、$$f\left(x\right)~GP\left(0,k\left(x,x^{\prime }\right)\right)$$、つまり f(x) は、共分散関数 $$k\left(x,x^{\prime }\right)$$ をもつゼロ平均 GP から派生しています。h(x) は、R^d における元の特徴ベクトル x を R^p における新しい特徴ベクトル h(x) に変換する一連の基底関数です。β は、基底関数の係数から構成される p 行 1 列のベクトルです。このモデルは、GPR モデルを表します。応答 y のインスタンスは、次のようにモデル化できます。

$$P\left(y_i|f\left(x_i\right),x_i\right) ~N\left(y_i|h{\left(x_i\right)}^T\beta +f\left(x_i\right),{\sigma }^2\right)$$

したがって、GPR モデルは確率モデルです。各観測値 $$x_i$$ について潜在的変数 f(x_i) を導入したので、GPR モデルはノンパラメトリックになります。ベクトル形式では、このモデルは次と等価です。

$$P(y|f,X)~N(y|H\beta +f,{\sigma }^{2}I),$$

ここで、次のようになります。

$$X=\left(\begin{array}{c}{x}_1^T\\ x_2^T\\ ⋮\\ x_n^T\end{array}\right),y=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right),H=\left(\begin{array}{c}h\left(x_1^T\right)\\ h\left(x_2^T\right)\\ ⋮\\ h\left(x_n^T\right)\end{array}\right),f=\left(\begin{array}{c}f\left(x_1\right)\\ f\left(x_2\right)\\ ⋮\\ f\left(x_n\right)\end{array}\right).$$

GPR モデルにおける潜在的変数 $$f\left(x_1\right),f\left(x_2\right),\mathrm{...},f\left(x_n\right)$$ の同時分布は次のようになり、線形回帰モデルに近くなっています。

$$P(f|X)~N\left(f|0,K\left(X,X\right)\right),$$

ここで、$$K\left(X,X\right)$$ は次のようになります。

$$K\left(X,X\right)=\left(\begin{array}{cccc}k\left(x_1,x_1\right)& k\left(x_1,x_2\right)& \cdots & k\left(x_1,x_n\right)\\ k\left(x_2,x_1\right)& k\left(x_2,x_2\right)& \cdots & k\left(x_2,x_n\right)\\ ⋮& ⋮& ⋮& ⋮\\ k\left(x_n,x_1\right)& k\left(x_n,x_2\right)& \cdots & k\left(x_n,x_n\right)\end{array}\right).$$

通常、共分散関数 $$k\left(x,x^{\prime }\right)$$ は一連のカーネル パラメーターまたはハイパーパラメーター $$\theta$$ によりパラメーター表現されます。多くの場合、$$\theta$$ に依存することを明示的に示すため、$$k\left(x,x^{\prime }\right)$$ は $$k\left(x,x^{\prime }|\theta \right)$$ と記述されます。

fitrgp は、GPR モデルに学習させるときに、カーネル関数の基底関数係数 $$\beta$$、ノイズ分散 $${\sigma }^2$$ およびハイパーパラメーター $$\theta$$ をデータから推定します。基底関数、カーネル (共分散) 関数およびパラメーターの初期値を指定できます。

GPR モデルは確率的なので、学習済みのモデルを使用して予測区間を計算できます (predict と resubPredict を参照してください)。

学習済みの GPR モデルを使用して回帰誤差の計算もできます (loss と resubLoss を参照してください)。

GPR モデルの予測区間の比較

この例では、GPR モデルをノイズのないデータ セットとノイズを含むデータ セットに当てはめます。この例では、2 つの GPR 近似モデルの予測応答と予測区間を比較します。

関数 $\mathit{g}\left(\mathit{x}\right)=\mathit{x}\cdot \sin \left(\mathit{x}\right)$ から 2 つの観測データ セットを生成します。

rng('default') % For reproducibility
x_observed = linspace(0,10,21)';
y_observed1 = x_observed.*sin(x_observed);
y_observed2 = y_observed1 + 0.5*randn(size(x_observed));

y_observed1 の値にはノイズがなく、y_observed2 の値にはランダム ノイズが含まれています。

GPR モデルを観測されたデータ セットに当てはめます。

gprMdl1 = fitrgp(x_observed,y_observed1);
gprMdl2 = fitrgp(x_observed,y_observed2);

当てはめたモデルを使用して、予測応答と 95% の予測区間を計算します。

x = linspace(0,10)';
[ypred1,~,yint1] = predict(gprMdl1,x);
[ypred2,~,yint2] = predict(gprMdl2,x);

Figure のサイズを変更して、1 つの Figure に 2 つのプロットを表示します。

fig = figure;
fig.Position(3) = fig.Position(3)*2;

1 行 2 列のタイル表示チャート レイアウトを作成します。

tiledlayout(1,2,'TileSpacing','compact')

各タイルに、観測されたデータ点の散布図と $\mathit{x}\cdot \sin \left(\mathit{x}\right)$ の関数プロットを描画します。次に、GP 予測応答のプロットと予測区間のパッチを追加します。

nexttile
hold on
scatter(x_observed,y_observed1,'r') % Observed data points
fplot(@(x) x.*sin(x),[0,10],'--r')  % Function plot of x*sin(x)
plot(x,ypred1,'g')                  % GPR predictions
patch([x;flipud(x)],[yint1(:,1);flipud(yint1(:,2))],'k','FaceAlpha',0.1); % Prediction intervals
hold off
title('GPR Fit of Noise-Free Observations')
legend({'Noise-free observations','g(x) = x*sin(x)','GPR predictions','95% prediction intervals'},'Location','best')

nexttile
hold on
scatter(x_observed,y_observed2,'xr') % Observed data points
fplot(@(x) x.*sin(x),[0,10],'--r')   % Function plot of x*sin(x)
plot(x,ypred2,'g')                   % GPR predictions
patch([x;flipud(x)],[yint2(:,1);flipud(yint2(:,2))],'k','FaceAlpha',0.1); % Prediction intervals
hold off
title('GPR Fit of Noisy Observations')
legend({'Noisy observations','g(x) = x*sin(x)','GPR predictions','95% prediction intervals'},'Location','best')

観測値にノイズがない場合、GPR 近似の予測応答は観測値と交差します。予測応答の標準偏差はほぼゼロです。したがって、予測区間は非常に狭くなります。観測値にノイズが含まれる場合、予測応答は観測値と交差せず、予測間隔が広くなります。

参照

[1] Rasmussen, C. E. and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press. Cambridge, Massachusetts, 2006.

参考

fitrgp | RegressionGP | predict

関連するトピック

• 厳密 GPR 法
• GPR モデルのデータ サブセット近似
• GPR モデルの回帰変数サブセット近似
• GPR モデルの完全独立条件近似
• GPR モデルのブロック座標降下近似