ガウス過程(1) 関数y(x)の事前分布のノンパラメトリックな定義 (Python)

PRMLでガウス過程を勉強していて詰まってしまったので[こちらのサイト](https://qiita.com/tyoshitake/items/7f1a75c43782f10d334b)を参考に実装してみる。

## ガウス過程

入力\{x_1,...,x_N\}に対して\{y_1,...,y_N\}を出力する関数をfとする。つまり、$$y_n=f(x_n)$$このとき、関数fは正規分布からランダムに生成される。

$$
f\sim \mathcal{N}(\boldsymbol{\mu},\Sigma)
$$

これは関数空間上の確率分布であって、無限次元。

今、入力xは1次元とする。このときfを特徴ベクトル
$$\boldsymbol{\phi}(x)=\begin{pmatrix}\phi_1(x)\\\ \vdots\\\ \phi_M(x)\end{pmatrix}$$
を用いて

$$
f(x)=\sum_iw_i\phi_i(x)=\ ^T\boldsymbol{w}\boldsymbol{\phi}(x)
$$

と、M次線形和として表す(Mは無限でも可)。今、入力\{x_1,...,x_N\}に対する出力\{y_1,...,y_N\}とは、

```math
\begin{align}
\boldsymbol{y}=\begin{pmatrix}y_1\\ \vdots\\ y_N\end{pmatrix}&=\begin{pmatrix}\phi_1(x_1) & \ldots&\phi_M(x_1)\\ \vdots&\ddots&\vdots\\ \phi_1(x_N)&\ldots&\phi_M(x_N)\end{pmatrix}\begin{pmatrix}w_1\\ \vdots\\ w_M\end{pmatrix}\\
&=\begin{pmatrix}^T\boldsymbol{\phi}(x_1)\\ \vdots\\ ^T\boldsymbol{\phi}(x_N)\end{pmatrix}\begin{pmatrix}w_1\\ \vdots\\ w_M\end{pmatrix}\\
&=\Phi\boldsymbol{w}
\end{align}
```

である。

今、$\boldsymbol{w}$の事前分布として

$$
p(\boldsymbol{w})=\mathcal{N}(\boldsymbol{0},\alpha^{-1}I)
$$

を仮定すれば、y_nはwの線形結合より$\boldsymbol{y}$は正規分布し、その平均と分散は

```math
\begin{align}
E[\boldsymbol{y}]&=\Phi E[\boldsymbol{w}]=\boldsymbol{0}\\
cov[\boldsymbol{y}]&=E[\boldsymbol{y}\ ^T\boldsymbol{y}]-E[\boldsymbol{y}]\ ^TE[\boldsymbol{y}]\\
&=E[\boldsymbol{y}\ ^T\boldsymbol{y}]\ \ \ \ \ \  (\ =E[\begin{pmatrix}y_1\\ \vdots\\ y_N\end{pmatrix}\begin{pmatrix}y_1& \ldots& y_N\end{pmatrix}]=\begin{pmatrix}
E[y(x_1)y(x_1)] & \ldots&E[y(x_1)y(x_N)]\\
\vdots & & \vdots\\
E[y(x_N)y(x_1)] & \ldots & E[y(x_N)y(x_N)]
\end{pmatrix}\ )\tag{1}\\
&=E[\Phi\boldsymbol{w}\ ^T\boldsymbol{w}\ ^T\Phi]\\
&=\Phi E[\boldsymbol{w}\ ^T\boldsymbol{w}]\ ^T\Phi\\
&=\frac{1}{\alpha}\Phi\ ^T\Phi\\
&=\frac{1}{\alpha}\begin{pmatrix}^T\boldsymbol{\phi}(x_1)\\ \vdots\\ ^T\boldsymbol{\phi}(x_N)\end{pmatrix}\begin{pmatrix}\boldsymbol{\phi}(x_1)&\ldots&\boldsymbol{\phi}(\boldsymbol{x}_N)\end{pmatrix}\\
&=\frac{1}{\alpha}\begin{pmatrix}^T\boldsymbol{\phi}(x_1)\boldsymbol{\phi}(x_1)&\ldots& ^T\boldsymbol{\phi}(x_1)\boldsymbol{\phi}(x_1)\\
\vdots&\ddots&\vdots\\
^T\boldsymbol{\phi}(x_N)\boldsymbol{\phi}(x_1)&\ldots&^T\boldsymbol{\phi}(x_N)\boldsymbol{\phi}(x_N)\end{pmatrix}\\
&=\frac{1}{\alpha}\begin{pmatrix}k(x_1,x_1) & \ldots & k(x_1,x_N)\\ \vdots&&\vdots\\ k(x_N,x_1) & \ldots&k(x_N,x_N)\end{pmatrix} \tag{2}
\end{align}
```

以上より、$$cov[\boldsymbol{y}]=K$$とおけば、入力$\\{x_1,...,x_N\\}$に対する出力$\\{y_1,...,y_N\\}$の同時事前分布は

$$
p(\boldsymbol{y})=\mathcal{N}(\boldsymbol{0},K)
$$

となる。

$k(x_n,x_m)=\ ^T\boldsymbol{\phi}(x_m)\boldsymbol{\phi}(x_n)$はカーネル関数であり、点x_m,x_nの類似度と解釈できる。(1)=(2)より、

$$
k(x_n,x_m)=\alpha E[y(x_n)y(x_m)]
$$

が成り立つ。すなわち、入力ベクトルx_n,x_mの特徴空間上の類似度が高い程、出力$y(x_n),y(x_m)$の相関も高くなることを意味する。今回は、ガウスカーネル

$$
k(x,x')=\exp(-\frac{\\|x-x'\\|^2}{\theta})
$$

を用いる。

以上より、ガウス過程を用いることで入力ベクトル\{x_n\}だけを用いて関数$y$の事前分布を得ることができる。

## 実装

```python
import numpy as np
import matplotlib.pyplot as plt

def gaussian_kernel_matrix(x, theta):
    """
    カーネル行列を計算する。
    xはN個の入力ベクトル(np.array)で、thetaはスカラー。
    i行j列がx[i]とx[j]のカーネルであるような(N,N)配列を返す。
    """
    
    N = x.shape[0]
    x = x.reshape((N, 1))
    x_ = x.reshape((1, N))
    
    return np.exp( -(x-x_)**2 / theta )
    
"""    
 x = np.array([[1],
               [2],
               [3]])
 x_ = np.array([[1,2,3]])
 とすれば、ブロードキャスト演算により x-x_ は
[[ 0 -1 -2]
 [ 1  0 -1]
 [ 2  1  0]] 
 と、i行j列がx_i-x_jになる。さらに (x-x_)**2 は
[[0 1 4]
 [1 0 1]
 [4 1 0]]
 と、上の配列の各要素を2乗した形になる。さらに np.exp(-(x-x_)**2)は
[[1.         0.36787944 0.01831564]
 [0.36787944 1.         0.36787944]
 [0.01831564 0.36787944 1.        ]]
 と、各要素の指数を取る。これがガウスカーネル行列である。
"""

def gp_sampling(x, alpha, theta):
    """
    ガウス過程からのサンプリングを行う。
    xはN個の入力ベクトル(np.array)で、theta、alphaはスカラー。(N,)配列を返す。
    """
    
    N = x.shape[0]
    K = gaussian_kernel_matrix(x, theta) / alpha
    mean = np.zeros((N))
    
    return np.random.multivariate_normal(mean, K)

np.random.seed(2)

fig = plt.figure(figsize=(15, 5))

x = np.linspace(0, 10, 200, endpoint=False)
theta_list = [np.array([0.25]), np.array([0.5]), np.array([2.5])]

y_list = []
for theta in theta_list:
    for _ in range(3):
        y_list.append(gp_sampling(x, alpha = 1, theta = theta))

for i, y in enumerate(y_list):
    if i % 3 == 0:
        ax = fig.add_subplot(1,3,i//3+1)
    ax.plot(x, y, linestyle = '-')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(f'Gaussian Kernel, theta = {theta_list[i//3][0]}')
```

![スクリーンショット 2021-02-15 20.33.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/686685/8e8b8341-abf3-934b-aa85-24ff1eecbb33.png)

続き：[ガウス過程回帰](https://qiita.com/_R2_/items/fd5a8602470b4e6c6ea8)

## 参考文献

パターン認識と機械学習(下) C.M.ビショップ

[ガウス過程の実装](https://qiita.com/tyoshitake/items/7f1a75c43782f10d334b)