ガウス過程と回帰モデル（線形~線形基底~ガウス過程）

f:id:yumaloop:20190203151221g:plain

1. ガウス過程（GP, Gaussian Process）
2. さまざまな回帰モデルの比較
3. ガウス過程回帰をPythonで実装
4. 参考書籍

※ このポストの内容は、PRML下巻6章とGPML2章を参考にしています。

www.microsoft.com www.gaussianprocess.org

1. ガウス過程（GP, Gaussian Process）

1.1 ガウス分布の共役性

　任意の確率モデル $p(\cdot|w)$ に対して $w$ の事前分布 $p(w)$ としてガウス分布を仮定した場合、ベイズの定理によって求められる $w$ の事後分布 $p(w|X^{n})$ もガウス分布に従います。つまり、ガウス分布に従うパラメータ $w$ について、任意の確率モデル $p(x|w)$ とデータサンプル $X^{n}$ から、事後分布 $p(w|X^{n})$ を求める操作はガウス分布の中で閉じています。一般に、「事前分布としてある確率分布を仮定すると、事後分布も同じ確率分布に従う」という性質を共役性といい、ガウス分布は共役性をもちます。

$$ \begin{eqnarray} (Posterior~of~w) & = & \frac{(likelihood ~of~ X^{n})\times(Prior~of~w)}{(Evidence~of~the~model)} \\ \\ p(w|X^{n}) & = & \frac{p(X^{n}|w)p(w)}{\int p(X^{n}|w)p(w) ~dw} \end{eqnarray} $$ $$ p(w) \sim N(\mu, s) ~ \Rightarrow ~ p(w|X^{n}) \sim N(\mu', s') $$

1.2 ガウス分布の切断と周辺化

　それぞれがガウス分布に従う「互いに独立でない」確率変数列 $Y = \left\{ Y_1, \dots Y_N \right\}$ を考えます。このとき、確率変数列 $Y$ はNコの確率変数 $\left\{ Y_i \right\}_{i=1 \sim N}$ のペアであるため、確率変数列 $Y$ そのものを $\mathbb{R}^{N}$ 次元の確率変数として扱うことができます。同様に、確率変数列 $Y$ における「任意の $m$ コの部分列」は $\mathbb{R}^{m}$ 次元の確率変数として扱えます。

　平均が $\mu$ で、共分散行列が $\Sigma$ のN次元ガウス分布に従う確率変数 $X$ ：

$X \in \mathbb{R}^{N}, ~~~ X \sim N_{N} \bigl( \mu, \Sigma \bigr)$

に対して、 $X_a, X_b$ （ $a+b=N$ ）による分割を考えます。

$X = \left[ \begin{array}{c} X_a \\ X_b \end{array} \right] \sim$ $N_{a+b} \left( \left[ \begin{array}{c} \mu_a \\\ \mu_b \end{array} \right], ~ \left[ \begin{array}{c} \Sigma_{aa} \\ \Sigma_{ba} \end{array} \begin{array}{c} \Sigma_{ab} \\ \Sigma_{bb} \end{array} \right] \right) \\\$

$\mu = \left[ \begin{array}{c} \mu_a \\\ \mu_b \end{array} \right], ~~~ \Sigma = \left[ \begin{array}{c} \Sigma_{aa} \\ \Sigma_{ba} \end{array} \begin{array}{c} \Sigma_{ab} \\ \Sigma_{bb} \end{array} \right]$

このとき、 $X_a, X_b$ は以下のようなガウス分布に従います。

$X_{a} \sim N_a \left( \mu_a, \Sigma_{aa} \right) \\ X_{b} \sim N_b \left( \mu_b, \Sigma_{bb} \right)$

さらに、「 $X$ の共分散行列 $\Sigma$ は半正定値対称行列である」ことから、次式が成り立ちます。

$\Sigma_{ab} = \Sigma_{ba}^{-1}$

ここで、「 $X_a, X_b$ は独立ではない」ことに注意すると、条件つき確率 $p(X_a | X_b)$ は次のように与えられ、条件つき平均 $\mu_{a|b}$ と条件つき共分散行列 $\Sigma_{a|b}$ からなるガウス分布に従うことがわかります。

$X_{a} | X_{b} \sim N \left( \mu_{a|b}, \Sigma_{a|b} \right) \\$

$\begin{eqnarray} \left\{ \begin{array}{c} \mu_{a|b} & = & \mu_a + \Sigma_{ba}{\Sigma_{aa}}^{-1} (X_b - \mu_b) \\ \Sigma_{a|b} & = & \Sigma_{aa} - \Sigma_{ab} {\Sigma_{bb}}^{-1} \Sigma_{ba} \end{array} \right. \end{eqnarray}$

$X_a$ を手元にあるデータサンプル、 $X_b$ を新しいデータサンプルと捉えると、これらを結合した $X_{a+b}$ に対してガウス過程を導入することができます。

1.3 ガウス過程の定義

　任意の組が同時ガウス分布（joint Gaussian distribution）に従うような確率変数の集合 $\left\{ f(X_1), \cdots, f(X_n) \right\}$ を，ガウス過程（Gaussian Process, GP）と定義します。ガウス分布とガウス過程を比較すると以下のようになります。

ガウス分布（Gaussian distribution）
- 任意のガウス分布 $N()$ は、平均 $\mu$ と共分散行列 $\Sigma$ によって完全に記述される。

${\bf X} \sim N \left( \mu, \Sigma \right) as: \\\$

$\begin{eqnarray} {\bf \mu} & := & \mathbb{E} \left[ {\bf X} \right] \\ \Sigma & := & \mathbb{C o v} \left[ X_i, X_j \right] = \mathbb{E} \bigl[ X_i - \mu_i \bigr] \mathbb{E} \bigl[ X_j - \mu_j \bigr] \end{eqnarray}$

ガウス過程（Gaussian Process, GP）
- 任意のガウス過程 $\mathcal{GP}()$ は、平均関数 $m({\bf x})$ と共分散関数 $k({\bf x}, {\bf x}')$ によって完全に記述される。

$f({\bf x}) \sim \mathcal{GP} \left( m({\bf x}), ~ k({\bf x}, {\bf x}') \right) as: \\\$

$\begin{eqnarray} m({\bf x}) & := & \mathbb{E} \left[ f({\bf x}) \right] \\\ k({\bf x}, {\bf x}') & := & \mathbb{C o v} \left[ f({\bf x}), f({\bf x}') \right] = \mathbb{E} \bigl[ f({\bf x}) - m({\bf x}) \bigr] \mathbb{E} \bigl[ f({\bf x}') - m({\bf x}') \bigr] \end{eqnarray}$

2. さまざまな回帰モデルの比較

2.1 回帰問題とは？

　まず、回帰（Regression）に一般的な定義を与えてみます。

　確率変数 $X \in \mathbb{R}^{d}, Y \in \mathbb{R}$ の実現値からなるNコのデータサンプル：

$$ \mathcal{D} := \left\{ \left( \begin{array}{c} {\bf x}_1 \\ y_1 \end{array} \right), \dots, \left( \begin{array}{c} {\bf x}_N \\ y_N \end{array} \right) \right\} $$

を用いて、 $X$ の値から $Y$ の値を推定（予測）するモデル：

$$ f: X \in \mathbb{R}^{d} \to Y \in \mathbb{R} $$

を作りたい。

ここで、モデル $f$ による推定値 $f({\bf x}_i)$ と真値 $y_i$ との推定誤差 $\varepsilon_i$ を次のように定義し、

$$ \varepsilon_i := y_i - f({\bf x}_i) ~~~ (i=1, \dots, N) $$

さらに、各 $\varepsilon_i ~({\small i=1, \dots, N})$ はそれぞれ独立に同じパラメータをもつガウス分布に従うと仮定すると、

$$ y_i = f({\bf x}_i) + {\varepsilon}_i ~~~ (i=1, \dots N) \\ {\varepsilon}_i \sim N \left( 0, \frac{1}{\beta} \right) $$

という等式が得られます。これをまとめると、

$$ p( y_i | {\bf x}_i, f) = N \left( f({\bf x}_i), \frac{1}{\beta} \right) $$

となり、回帰問題の目的を、「 $f$ を定めることで，任意の $(x_i, y_i)$ に対する条件付き確率 $p({y_i | x_i, f})$ を定めること」とベイズ的に解釈することも可能です。すなわち、次式（ベイズの定理）で、 $p( f | y_i, {\bf x}_i )$ はデータサンプル $\mathcal{D}$ への「当てはまりの良さ（=尤度）」を表しますが、同時に、 $p( f | y_i, {\bf x}_i )$ によって、データサンプル $\mathcal{D}$ から $f$ の事後確率と予測分布が決定されます。

$$ p( f | y_i, {\bf x}_i ) = \frac{ p( y_i | {\bf x}_i, f ) p(f) } {p( y_i | {\bf x}_i )} $$

　一般化に、回帰問題では、データサンプル $\mathcal{D} = \left\{ {\bf X}, {\bf y} \right\}$ とモデル $f$ に対して

$$ p( f | {\bf y}, {\bf X} ) = \frac{ p( {\bf y} | {\bf X}, f ) p(f) } {p( {\bf y} | {\bf X} )}, ~~ i.e. ~~ p( f | \mathcal{D} ) = \frac{ p( \mathcal{D} | f ) p(f) } {p( \mathcal{D} )} $$

という構造が現れます。次章以降で、具体的な回帰モデルについて、解説したいと思います。

2.2 線形回帰モデル

2.2.1 線形回帰モデルの準備

　もっとも　素朴な回帰モデルとして「線形回帰モデル」を考えます。モデルパラメータとして、重みベクトル ${\bf w} \in \mathbb{R}^{d}$ を想定します。

$$ y_i = f({\bf x}_i) + \varepsilon_i, \\ f({\bf x}_i) = {\bf x}_i^{\rm T} {\bf w} = \left( \begin{array}{c} x_{i1}, \dots, x_{id} \end{array} \right) \cdot \left( \begin{array}{c} w_{1} \\ \vdots \\ w_{d} \end{array} \right), ~~~ \varepsilon_i \sim N \left( 0, \frac{1}{\beta} \right) $$

これらの関係式をまとめると、各 ${\bf x}_i$ に対する $y_i$ の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( y_i | x_i, {\bf w} ) = N \left( y_i | f({\bf x}_i), \frac{1}{\beta}{\bf I} \right) = N \left( {\bf x}_i^{\rm T} {\bf w}, \frac{1}{\beta} \right) $$

さらに、個々データではなくデータセット $\mathcal{D}$ 全体に対して、回帰モデル $f$ はどのように作用するのかをまとめて考えます。計画行列（design matrix） ${\bf X} \in \mathbb{R}^{d ~ \times ~ N}$ と、ベクトル ${\bf y}$ を以下のように定義します。

$$ {\bf X} := \left( {\bf x}_1, {\bf x}_2, \dots, {\bf x}_N \right)= \left( \begin{array}{cccc} x_{11} & x_{21} & \dots & x_{N1}\\ x_{12} & x_{22} & \dots & x_{N2}\\ \vdots & \vdots & \dots & \vdots \\ x_{1N} & x_{2N} & \dots & x_{NN} \end{array} \right), ~~~ {\bf y} := \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right), ~~~ {\boldsymbol \varepsilon} := \left( \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{array} \right) $$

すると、「線形回帰モデル」は以下のように記述され、

$$ {\bf y} = f({\bf X}) + {\boldsymbol \varepsilon}, \\ f({\bf X}) = {\bf X}^{\rm T}{\bf w}, ~~~ \boldsymbol \varepsilon \sim N \left( {\bf 0}, \frac{1}{\beta}{\bf I}\right) $$

これらの関係式をまとめると、データセット $\mathcal{D}$ における ${\bf y}$ の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( {\bf y} | {\bf X}, {\bf w} ) = N \left( {\bf y} | f({\bf X}), \frac{1}{\beta}{\bf I} \right) = N \left( {\bf X}^{\rm T} {\bf w}, \frac{1}{\beta}{\bf I} \right) $$

2.2.2 線形回帰モデルのベイズ的解釈

　線形回帰モデルにおいて、 ${\bf y}$ を推定する確率モデル $p( {\bf y} | {\bf X}, {\bf w} )$ の性能はパラメータ ${\bf w}$ の値によって左右されます。確率モデル $p( {\bf y} | {\bf X}, {\bf w} )$ をベイズの枠組み:

$$ p({\bf w} | {\bf y}, {\bf X} ) = \frac{ p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) } {p({\bf y} | {\bf X} )} = \frac{ p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) } { \int p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) ~ d{\bf w}} $$

によって拡張し、パラメータ ${\bf w}$ の事後分布 $p({\bf w} | {\bf y}, {\bf X} )$ を導出することで、N+1 番目の未知のデータサンプル $\left( \begin{array}{c} {\bf x}^{*} \\ y^{*} \end{array} \right)$ に対する「ばらつきを持った確率モデル」である予測分布を導出しましょう。

（１） ${\bf w}$ の事前分布

$$ p({\bf w}) = N \left( {\bf w} ~ | ~ {\bf 0}, ~ \frac{1}{\alpha} {\bf I} \right) $$

（２）尤度 $p({\bf y} | X, {\bf w})$

$$ \begin{eqnarray} p({\bf y} | X, {\bf w}) & = &\prod_{i=1}^{N} p(y_i | {\bf x}_i, {\bf w}) \\ & = & \prod_{i=1}^{N} \frac{1}{ \sqrt{ 2 \pi \beta^{-1}}} \exp \left( - \frac{{\left( y_i - {{\bf x}_i}^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & \frac{1}{ {\left( 2 \pi \beta^{-1} \right)}^{n/2}} \exp \left( - \frac{{\left( {\bf y} - X^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & N \left( {\bf y} ~ | ~ X^{\rm T}{\bf w}, ~ \frac{1}{\beta}{\bf I} \right) \end{eqnarray} $$

（３） ${\bf w}$ の事後分布

$$ \begin{eqnarray} p({\bf w} | {\bf X}, {\bf y}) & \propto & \exp \left( - \frac{1}{2} { ( {\bf y} - {\bf X}^{\rm T}{\bf w} ) }^{\rm T} { \left( \frac{1}{\beta} {\bf I} \right) }^{-1} ( {\bf y} - {\bf X}^{\rm T}{\bf w} ) \right) \exp \left( - \frac{1}{2} {\bf w}^{\rm T} { \left( \frac{1}{\alpha} {\bf I} \right) }^{-1} {\bf w} \right) \\ & \propto & \exp \left( { ( {\bf w} - {\hat{\bf w}}_{MAP} ) }^{\rm T} \left( \beta {\bf X}{\bf X}^{\rm T} + \alpha {\bf I} \right) ( {\bf w} - {\hat{\bf w}}_{MAP} ) \right) \\ & = & N \left( {\bf w} ~ | ~ {\hat{\bf w}}_{MAP}, ~ {\bf A}^{-1} \right) \end{eqnarray} $$

　ただし、

$$ \begin{eqnarray} {\hat{\bf w}}_{MAP} & := & \beta {\bf A }^{-1} {\bf X}{\bf y} \\ {\bf A} & := & \beta {\bf X}{\bf X}^{\rm T} + \alpha {\bf I} \end{eqnarray} $$

（４） $y^{*}$ の予測分布

$$ \begin{eqnarray} p^{*}(y^{*} | {\bf x}^{*}, {\bf X}, {\bf y}) & = & \int p(y^{*} | {\bf x}^{*}, {\bf w}) p( {\bf w} | X, {\bf y}) ~ d{\bf w} \\ & = & N \left( y^{*} ~|~ \beta {{\bf x}^{*}}^{\rm T} { \bf A }^{-1} {\bf X} {\bf y}, ~ { {\bf x}^{*}}^{\rm T} { \bf A }^{-1} {\bf x}^{*} \right) \end{eqnarray} $$

2.3 線形基底モデル

2.3.1 線形基底モデルの準備

　素朴な線形回帰モデルでは，入力ベクトル ${\bf x} \in \mathbb{R}^d$ のあるデータ空間（data space） $\subset \mathbb{R}^d$ の中で，出力 $y \in \mathbb{R}$ を推定します．しかし，この場合，データ空間を構成する基底ベクトルの線形結合によって $y$ を推定するため, モデルの表現力が制約されてします。そこで，入力ベクトル ${\bf x}$ を，基底関数（basis function） $\phi: \mathbb{R}^d \to \mathbb{R}^M$ によって特徴量 $\phi({\bf x}) \in \mathbb{R}^M$ へと変換し，高次元の特徴空間 $\subset \mathbb{R}^{M}$ （feature space）から $y \in \mathbb{R}$ の値を推定することを考えます．ここで，基底関数 $\phi$ による変換は線形でも非線形でもよいとします．

$$ {\bf x}_i \in \mathbb{R}^d \to \phi({\bf x}_i) \in \mathbb{R}^M \underset{\large \bf w}{\longrightarrow} f({\bf x}_i) \in \mathbb{R} \underset{ {\large \varepsilon_i} } {\longleftrightarrow} y_i \in \mathbb{R}, ~~~ (i=1, \dots, N) \\ $$

基底関数 $\phi: \mathbb{R}^d \to \mathbb{R}^M$ は、データを表現する基底を変える（空間を変える）関数と考えることができます。また、一般に、 $d \lt M$ であることから、 $\phi$ はデータサンプルを高次元空間に写像する（mapping）関数であることがわかります。

${\bf x}_i$ のデータ空間 $\subset \mathbb{R}^d$ における表現： $$ {\bf x}_i = \left( x_{i1}, x_{i2}, \dots, x_{id} \right) $$
${\bf x}_i$ の特徴空間 $\subset \mathbb{R}^M$ における表現： $$ \phi({\bf x}_i) = \left( \phi_{1}({\bf x}_i), \phi_{2}({\bf x}_i), \dots, \phi_{M}({\bf x}_i) \right) $$

実際に，「線形基底モデル」を構築しましょう．モデルパラメータとして，重みベクトル ${\bf w} \in \mathbb{R}^{M}$ を想定します．

$$ y_i = f({\bf x}_i) + \varepsilon_i \\ f({\bf x}_i) = { \phi ( {\bf x}_i ) }^{\rm T} {\bf w} = \left( \begin{array}{c} \phi_{1}({\bf x}_i), \dots, \phi_{M}({\bf x}_i) \end{array} \right) \cdot \left( \begin{array}{c} w_{1} \\ \vdots \\ w_{M} \end{array} \right), ~~~ \varepsilon_i \sim N \left( 0, \frac{1}{\beta} \right) $$

これらの関係式をまとめると、各 ${\bf x}_i$ に対する $y_i$ の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( y_i | x_i, {\bf w} ) = N \left( y_i | f( {\bf x}_i ), \frac{1}{\beta} \right) = N \left( y_i | {\phi ( {\bf x}_i ) }^{\rm T} {\bf w}, \frac{1}{\beta} \right) $$

さらに、個々データではなくデータセット $\mathcal{D}$ 全体に対して、回帰モデル $f$ はどのように作用するのかをまとめて考えます。計画行列（design matrix） $\Phi({\bf X}) \in \mathbb{R}^{M ~ \times ~ N}$ と、ベクトル ${\bf y}, {\boldsymbol \varepsilon}$ を以下のように定義します。

$$ \Phi({\bf X}) := \left( \phi({\bf x}_1), \phi({\bf x}_2), \dots, \phi({\bf x}_N) \right)= \left( \begin{array}{cccc} \phi_{1}({\bf x}_1) & \phi_{1}({\bf x}_2) & \dots & \phi_{1}({\bf x}_N) \\ \phi_{2}({\bf x}_1) & \phi_{2}({\bf x}_2) & \dots & \phi_{2}({\bf x}_N) \\ \vdots & \vdots & \dots & \vdots \\ \phi_{M}({\bf x}_1) & \phi_{M}({\bf x}_2) & \dots & \phi_{M}({\bf x}_N) \end{array} \right), ~~~ {\bf y} := \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right), ~~~ {\boldsymbol \varepsilon} := \left( \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{array} \right) $$

すると、「線形基底回帰モデル」は以下のように記述され、

$$ {\bf y} = f({\bf X}) + {\boldsymbol \varepsilon}, \\ f({\bf X}) = {\Phi({\bf X})}^{\rm T}{\bf w}, ~~~ \boldsymbol \varepsilon \sim N \left( {\bf 0}, \frac{1}{\beta}{\bf I}\right) $$

$$ p( {\bf y} | {\bf X}, {\bf w} ) = N \left( {\bf y} ~|~ f({\bf X}), \frac{1}{\beta}{\bf I} \right) = N \left( {\bf y} ~|~ {\Phi({\bf X})}^{\rm T}{\bf w}, \frac{1}{\beta}{\bf I} \right) $$

2.3.2 線形基底モデルのベイズ的解釈

　確率モデル $p( {\bf y} | {\bf X}, {\bf w} )$ をベイズの枠組み:

（１） ${\bf w}$ の事前分布

$$ p({\bf w}) = N \left( {\bf w} ~ | ~ {\bf 0}, ~ \frac{1}{\alpha} {\bf I} \right) $$

（２）尤度 $p({\bf y} | X, {\bf w})$

$$ \begin{eqnarray} p({\bf y} | \Phi({\bf X}), {\bf w}) & = &\prod_{i=1}^{N} p(y_i | \phi({\bf x}_i), {\bf w}) \\ & = & \prod_{i=1}^{N} \frac{1}{ \sqrt{ 2 \pi \beta^{-1}}} \exp \left( - \frac{{\left( y_i - {\phi({\bf x}_i)}^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & \frac{1}{ {\left( 2 \pi \beta^{-1} \right)}^{n/2}} \exp \left( - \frac{{\left( {\bf y} - {\Phi({\bf X})}^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & N \left( {\bf y} ~ | ~ {\Phi({\bf X})}^{\rm T}{\bf w}, ~ \frac{1}{\beta}{\bf I} \right) \end{eqnarray} $$

（３） ${\bf w}$ の事後分布 $p({\bf w} | \Phi({\bf X}), {\bf y})$

$$ \begin{eqnarray} p({\bf w} | \Phi({\bf X}), {\bf y}) & \propto & \exp \left( - \frac{1}{2} { ( {\bf y} - {\Phi({\bf X})}^{\rm T}{\bf w} ) }^{\rm T} { \left( \frac{1}{\beta} {\bf I} \right) }^{-1} ( {\bf y} - {\Phi({\bf X})}^{\rm T}{\bf w} ) \right) \exp \left( - \frac{1}{2} {\bf w}^{\rm T} { \left( \frac{1}{\alpha} {\bf I} \right) }^{-1} {\bf w} \right) \\ & \propto & \exp \left( { ( {\bf w} - {\hat{\bf w}}_{MAP} ) }^{\rm T} \left( \beta {\Phi({\bf X})}{\Phi({\bf X})}^{\rm T} + \alpha {\bf I} \right) ( {\bf w} - {\hat{\bf w}}_{MAP} ) \right) \\ & = & N \left( {\bf w} ~ | ~ {\hat{\bf w}}_{MAP}, ~ { A }^{-1} \right) \end{eqnarray} $$

　ただし、

$$ \begin{eqnarray} {\hat{\bf w}}_{MAP} & := & \beta { \bf A }^{-1} {\bf X}{\bf y} \\ A & := & \beta {\Phi({\bf X})}{\Phi({\bf X})}^{\rm T} + \alpha {\bf I} \end{eqnarray} $$

（４） $y^{*}$ の予測分布

$$ \begin{eqnarray} p^{*}(y^{*} | {\bf x}^{*}, \Phi({\bf X}), {\bf y}) & = & \int p(y^{*} | {\bf x}^{*}, {\bf w}) p( {\bf w} | \Phi({\bf X}), {\bf y}) ~ d{\bf w} \\ & = & N \left( y^{*} ~|~ \beta {{\bf x}^{*}}^{\rm T} { \bf A}^{-1} \Phi({\bf X}) {\bf y}, ~ { {\bf x}^{*}}^{\rm T} { \bf A}^{-1} {\bf x}^{*} \right) \end{eqnarray} $$

2.4 ガウス過程による線形基底モデルの表現

2.4.1 ガウス過程による回帰モデルの準備

　ガウス過程による回帰では、基底関数 $\phi()$ ではなく、回帰モデル $f({\bf X})$ の確率分布 $p(f({\bf X})) = p( {\Phi({\bf X})}^{\rm T}{\bf w} )$ の共分散行列 $K_{XX}$ を事前に定めます。 $K_{XX}$ をグラム行列（Gram Matrix）といいます。

$$ {\bf K}_{XX} = \left( \begin{array}{cccc} k({\bf x}_1, {\bf x}_1) & k({\bf x}_1, {\bf x}_2) & \dots & k({\bf x}_1, {\bf x}_N) \\ k({\bf x}_2, {\bf x}_1) & k({\bf x}_2, {\bf x}_2) & \dots & k({\bf x}_2, {\bf x}_N) \\ \vdots & \vdots & \ddots & \vdots \\ k({\bf x}_N, {\bf x}_1) & k({\bf x}_N, {\bf x}_2) & \dots & k({\bf x}_N, {\bf x}_N) \\ \end{array} \right) $$ ただし、 $$ \begin{eqnarray} k({\bf x}_i, {\bf x}_j) & = & C o v \left[ f( {\bf x}_i ), f( {\bf x}_j ) \right] \\ & = & { \phi( {\bf x}_i )}^{\rm T} \left( \frac{1}{\alpha} {\bf I} \right) \phi( {\bf x}_j ) \\ & = & \frac{1}{\alpha} \phi_{1}( {\bf x}_i ) \phi_{1}( {\bf x}_j ) + \cdots + \frac{1}{\alpha} \phi_{M}( {\bf x}_i ) \phi_{M}( {\bf x}_j ) \end{eqnarray} $$

　定義から、グラム行列 ${\bf K}_{XX}$ は $f({\bf x}_i)$ と $f({\bf x}_j)$ の分散からなる行列（=共分散行列）です。すなわち、グラム行列 ${\bf K}_{XX}$ は、基底関数（basis function） $\phi()$ と計画行列（design matrix） $\Phi({\bf X})$ および重みベクトル ${\bf w}$ の事前分布 $p(w) \sim N \left( {\bf 0}, \frac{1}{\alpha}{\bf I} \right)$ によっても記述されます。

$$ \begin{eqnarray} {\bf K}_{XX} & = & \left( \begin{array}{cccc} C o v \left[ f({\bf x}_1), f({\bf x}_1) \right] & C o v \left[ f({\bf x}_1), f({\bf x}_2) \right] & \dots & C o v \left[ f({\bf x}_1), f({\bf x}_N) \right] \\ C o v \left[ f({\bf x}_2), f({\bf x}_1) \right] & C o v \left[ f({\bf x}_2), f({\bf x}_2) \right] & \dots & C o v \left[ f({\bf x}_2), f({\bf x}_N) \right] \\ \vdots & \vdots & \ddots & \vdots \\ C o v \left[ f({\bf x}_N), f({\bf x}_1) \right] & C o v \left[ f({\bf x}_N), f({\bf x}_2) \right] & \dots & C o v \left[ f({\bf x}_N), f({\bf x}_N) \right] \\ \end{array} \right) \\ \\ & = & {\Phi({\bf X})}^{\rm T} \left( \frac{1}{\alpha}{\bf I} \right) \Phi({\bf X}) \\ & = & \frac{1}{\alpha} \Phi({\bf X}) { \Phi({\bf X}) }^{\rm T} \end{eqnarray} $$

すなわち、同一の線形基底モデルは，以下の2つの視点から記述することが可能です．ガウス過程を含む「カーネル法」は，後者の関数空間による視点（function space view）に根ざした多変量解析手法です．

基底関数 $\phi()$ と計画行列 $\Phi({\bf X})$ によって表現する（weight space view）

モデル設計者は，基底関数 $\phi()$ を事前に定める．

グラム行列 ${\bf K}_{XX}$ によって表現する（function space view）

モデル設計者は，カーネル関数 $k()$ を事前に定める．

以上より，モデルの出力値 $f({\bf X})$ と変数 ${\bf y}$ の確率分布は，次のように導出されます，

$$ \begin{eqnarray} p \bigl( f({\bf X}) \bigr) & = & N \left( f({\bf X}) ~|~ {\bf 0}, \frac{1}{\alpha}\Phi({\bf X}) {\Phi({\bf X})}^{\rm T} \right) \\ & = & N \Bigl( f({\bf X}) ~|~ {\bf 0}, {\bf K}_{XX} \Bigr) \\ \\ p \bigl( {\bf y} \bigr) & = & N \left( {\bf y} ~|~ {\bf 0}, \frac{1}{\alpha}\Phi({\bf X}) {\Phi({\bf X})}^{\rm T} + \frac{1}{\beta} {\bf I} \right) \\ & = & N \Bigl( {\bf y} ~|~ {\bf 0}, {\bf K}_{XX} + \frac{1}{\beta} {\bf I} \Bigr) \\ \end{eqnarray} $$

2.4.2 ガウス過程による回帰モデルによる予測

　ここで、N+1 番目の未知のデータ ${\bf x}^{*}$ に対する $y^{*}$ の予測値 $f( {\bf x}^{*} )$ を導出しましょう。

$$ y^{*} = f( {\bf x}^{*} ) + {\varepsilon}^{*} \\ f( {\bf x}^{*} ) = { \phi( {\bf x}^{*} ) }^{\rm T} {\bf w}, ~~~ {\varepsilon}^{*} \sim N \left( 0, \frac{1}{\beta} \right) $$

2.3.2 の結果から、 $f ( {\bf x}^{*} )$ の予測分布は以下のように記述されます。

$$ f ( {\bf x}^{*} ) \sim N \Bigl( \mathbb{E} \left[ f ( {\bf x}^{*} ) \right] , \mathbb{V} \left[ f ( {\bf x}^{*} ) \right] \Bigr) $$ $$ \begin{eqnarray} \left\{ \begin{array}{c} \mathbb{E} \left[ f ( {\bf x}^{*} ) \right] & = & { {\bf k}^{*} }^{\rm T} { \left( {\bf K}_{XX} + \frac{1}{\beta} {\bf I} \right) }^{-1} {\bf y} \\ \mathbb{V} \left[ f ( {\bf x}^{*} ) \right] & = & k ( {\bf x}^{*}, {\bf x}^{*} ) - { {\bf k}^{*} }^{\rm T} { \left( {\bf K}_{XX} + \frac{1}{\beta} {\bf I} \right) }^{-1} {\bf k}^{*} \end{array} \right. \end{eqnarray} $$

未知のデータ ${\bf x}^{*}$ をガウス過程 $\mathcal{GP} \left( {\bf 0}, {\bf K}_{XX} \right)$ に加えると、

$$ \begin{eqnarray} {{\bf K}_{XX}}^{*} & = & \left( \begin{array}{cccc} k( {\bf x}_1, {\bf x}_1 ) & \dots & k( {\bf x}_1, {\bf x}_N ) & k( {\bf x}_1, {\bf x}^{*} ) \\ \vdots & \ddots & \vdots & \vdots \\ k( {\bf x}_N, {\bf x}_1 ) & \dots & k( {\bf x}_N, {\bf x}_N ) & k( {\bf x}_N, {\bf x}^{*} ) \\ k( {\bf x}^{*}, {\bf x}_1 ) & \dots & k( {\bf x}^{*}, {\bf x}_N ) & k( {\bf x}^{*}, {\bf x}^{*} ) \\ \end{array} \right) \\ \\ & = & \left( \begin{array}{cc} {\bf K}_{XX} & {\bf k}^{*} \\ { {\bf k}^{*} }^{\rm T} & k( {\bf x}^{*}, {\bf x}^{*} ) \end{array} \right) \end{eqnarray} $$

ただし、 ${\bf k}^{*}$ は以下のように記述されるベクトルです。 $$ \begin{eqnarray} {\bf k}^{*} & := & \left( \begin{array}{c} k( {\bf x}_1, {\bf x}^{*} ) \\ \vdots \\ k( {\bf x}_N, {\bf x}^{*} ) \end{array} \right) \\ \\ & = & \left( \begin{array}{c} C o v \left[ f( {\bf x}_1 ), f( {\bf x}^{*} ) \right] \\ \vdots \\ C o v \left[ f( {\bf x}_N ), f( {\bf x}^{*} ) \right] \end{array} \right) \\ \\ & = & \left( \begin{array}{c} { \phi( {\bf x}_1 )}^{\rm T} \left( \frac{1}{\alpha} {\bf I} \right) \phi( {\bf x}^{*} ) \\ \vdots \\ { \phi( {\bf x}_N )}^{\rm T} \left( \frac{1}{\alpha} {\bf I} \right) \phi( {\bf x}^{*} ) \end{array} \right) \\ \\ & = & \left( \begin{array}{c} \frac{1}{\alpha} \phi_{1}( {\bf x}_1 ) \phi_{1}( {\bf x}^{*} ) + \cdots + \frac{1}{\alpha} \phi_{M}( {\bf x}_1 ) \phi_{M}( {\bf x}^{*} ) \\ \vdots \\ \frac{1}{\alpha} \phi_{1}( {\bf x}_N ) \phi_{1}( {\bf x}^{*} ) + \cdots + \frac{1}{\alpha} \phi_{M}( {\bf x}_N ) \phi_{M}( {\bf x}^{*} ) \end{array} \right) \\ \\ & = & { \phi( {\bf x}^{*} ) }^{\rm T} \frac{1}{\alpha} {\bf I} ~ \Phi( {\bf X} ) \end{eqnarray} $$

3. ガウス過程回帰をPythonで実装

　実際に、ある関数からデータサンプルを作成して、ガウス過程による回帰モデル $f$ によってこれを近似します。今回は、以下のような平均関数 $m(x)$ と共分散関数 $k(x, x')$ を用いたガウス過程 $\mathcal{GP}$ と、それに従う回帰モデル $f$ を考えます。ここで、共分散関数としてRBFカーネル（Radial basis function kernel, Gaussian kernel）を採用しています。

$$ f(x) \sim \mathcal{GP} \left( m(x), k(x, x') \right) $$

$\begin{eqnarray} \left\{ \begin{array}{c} m(x) & = & 0 \\ k(x, x') & = & \exp \left( - \frac{ { || x - x' || }^{2} }{2 {\sigma}^{2} } \right) = \exp \left( - \gamma { || x - x' || }^{2} \right) \end{array} \right. \end{eqnarray}$

この回帰モデルを用いて、関数：

$$ y = \sin(2 \pi x) + \cos(1.7 \times 2 \pi x) ~~~ (0 \leq x \leq 1) $$

を、生成されるデータサンプルから予測してみましょう。なお実装は、こちらのQuitta記事を参考にさせてもらいました。

3.1 RBFカーネルのためのクラス

class RBFKernel(object):
    
    def __init__(self, sigma):
        self.sigma = sigma
        
    def get_sigma(self):
        return np.copy(self.sigma)

    def __call__(self, x1, x2):
        return np.exp( (-1. / (self.sigma ** 2) ) *  (x1 - x2) ** 2)

    def derivatives(self, x1, x2):
        dif_sigma = np.exp( (-1. / (self.sigma ** 2) ) *  (x1 - x2) ** 2) * ( (x1 - x2) ** 2 ) / ( self.sigma ** 3)
        return dif_sigma

    def update_sigma(self, update):
        self.sigma += update

3.2 ガウス過程回帰を行うためのクラス

class GaussianProcessRegression(object):

    def __init__(self, kernel, x, y, beta=1.):
        self.kernel = kernel
        self.x = x
        self.y = y
        self.beta = beta
        
        # グラム行列を計算
        self.fit_kernel()

    def fit_kernel(self):
        ''' グラム行列を計算 '''
        
        # グラム行列 K
        self.gram = self.kernel(*np.meshgrid(self.x, self.x))
        # 共分散行列 K + 1/beta I
        self.covariance = self.gram + np.identity(len(self.x)) / self.beta
        # 精度行列 (共分散行列の逆行列)
        self.precision = np.linalg.inv(self.covariance)
    
    def predict_y(self, x_new):
        ''' y_new の予測値と分散を出力 '''
        
        # k*
        k = np.array([self.kernel(xi, x_new) for xi in self.x])

        # 平均と分散
        y_mean = k.dot(self.precision.dot(self.y))
        y_var = self.kernel(x_new, x_new) - k.dot(self.precision.dot(k))
        
        return y_mean, y_var

3.3 グラフの描画（100フレームのアニメーション）を作成する

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def func(x):
    ''' 回帰モデルによって近似を行う関数 '''
    return np.sin(2 * np.pi * x) + np.cos(1.7 * 2 * np.pi * x)


def create_data_sample(func, low=0, high=1., n=10, std=1.):
    ''' n 組のデータサンプル (x, y) を作成する '''
    x = np.random.uniform(low, high, n)
    y = func(x) + np.random.normal(scale=std, size=n)
    return x, y
    
def _update(frame, x_data, y_data):
    ''' フレームごとにグラフを描画する '''
        
    # 現在のグラフを消去する
    plt.cla()
    
    # データサンプルのframe数分の取得
    x = x_data[:frame]
    y = y_data[:frame]
    
    # 回帰モデルをつくる
    kernel = RBFKernel(sigma=0.5)
    regression = GaussianProcessRegression(kernel, x, y, beta=100)
    
    # 描画用にxのデータを生成
    x_test = np.linspace(0, 1, 100)
    
    # yの予測値の平均と分散を導出
    y_test=[]
    y_var=[]
    for x_new in x_test:
        y_new_mean, y_new_var = regression.predict_y(x_new)
        y_test.append(y_new_mean)
        y_var.append(y_new_var)
        
    # list を numpy.ndarrayに変換
    y_test = np.array(y_test)
    y_var = np.array(y_var)
    
    # グラフを描画
    plt.scatter(x, y, alpha=0.4, color="blue", label="data sample")
    plt.plot(x_test, func(x_test), color="blue", linestyle="dotted", label="actual function")
    plt.plot(x_test, y_test, color="red", label="predicted function")
    plt.fill_between(x_test, y_test - 3 * np.sqrt(y_var), y_test + 3 * np.sqrt(y_var), color="pink", alpha=0.5, label="95% confidence interval (2 * std)")
    plt.legend(loc="lower left")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.xlim([-0, 1])
    plt.ylim([-2.1, 2.1])
    plt.title("Gaussian Process Regression \n (RBF kernel, sigma=0.5, sample size = " + str(frame) + ")")
    
def main():
    # 描画領域
    fig = plt.figure(figsize=(8, 6))

    # データサンプルを作成
    x_data, y_data = create_toy_data(func, low=0, high=1, n=100, std=0.1)

    params = {
        'fig': fig,
        'func': _update,  # グラフを更新する関数
        'fargs': (x_data, y_data),  # 関数の引数 (フレーム番号を除く)
        'interval': 1,  # 更新間隔 (ミリ秒)
        'frames': np.arange(1, 101, 1),  # フレーム番号を生成するイテレータ, 1~100
        'repeat': False,  # 繰り返さない
    }

    # アニメーションを作成する
    ani = animation.FuncAnimation(**params)

    # アニメーションを保存する
    ani.save("./gpr_rbf.gif")
    # アニメーションを表示する
    # plt.show()