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

11900 ワード
PRML Python ガウス過程 機械学習 Python テキストリンク

PRMLでガウス過程を勉強していて詰まってしまったのでこちらのサイトを参考に実装してみる。

ガウス過程

入力{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}とは、

\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}$は正規分布し、その平均と分散は

\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$の事前分布を得ることができる。

実装

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]}')

続き：ガウス過程回帰

参考文献

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

ガウス過程の実装