経緯

Pythonのガウス過程フレームワークGPyを用いたガウス過程回帰の計算過程を追った

題材にしたのは下記のコード（参考）

import GPy
import numpy as np

X = np.linspace(0.05,0.95,10)[:,None]
Y = -np.cos(np.pi*X) +np.sin(4*np.pi*X) + np.random.randn(10,1)*0.05

k = GPy.kern.rbf(input_dim=1, variance=1., lengthscale=1.)
m = GPy.models.GPRegression(X,Y,k)

やりたいこと

入力$\boldsymbol{x}$に対する出力$\boldsymbol{y}$が得られているとき（サンプル数$n$）、未知の入力$x^\ast$に対する出力$y^\ast$の確率分布$p(y^\ast|x^\ast, \boldsymbol{x}, \boldsymbol{y})$を求める

統計の前提知識

1次元ガウス分布（正規分布）

$$
\mathcal{N}(\mu,\sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}}\exp{\bigg[-\frac{(x-\mu)^2}{2\sigma^2}\bigg]
}$$

$x$: 確率変数
$\mu$: 平均
$\sigma$: 標準偏差
$\sigma^2$: 分散

N次元ガウス分布（多変量正規分布）

$$
\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{N/2} \sqrt{\Sigma}}\exp{\bigg[-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{T}{\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})}\bigg]
}$$

\begin{equation*}
\boldsymbol{\Sigma} =
    \begin{pmatrix}
    \sigma_{1}^{2} & \sigma_{1,2} & \cdots & \sigma_{1,N} \\
    \sigma_{2,1} & \sigma_{2}^{2} & \ddots & \vdots \\
    \vdots & \ddots & \ddots & \sigma_{N-1,N}\\
    \sigma_{N,1} & \cdots & \sigma_{N,N-1} & \sigma_{N}^{2}
    \end{pmatrix}
\end{equation*}

$\boldsymbol{x}=(x_1, x_2, \cdots, x_n)^{T}$: 確率変数ベクトル
$\boldsymbol{\mu}=(\mu_1, \mu_2, \cdots, \mu_n)^{T}$: 平均ベクトル
$\sigma$: 標準偏差
$\boldsymbol{\Sigma}$: 分散共分散行列（または共分散行列）
$\sigma_{i,j}$: $x_i$と$x_j$の共分散
$\sigma_{i}^{2}$: $x_i$の分散

導出

入出力関係のモデル化

$\boldsymbol{x}$と$\boldsymbol{y}$の関係を、$m$個の説明変数を用いた特徴ベクトル（または計画行列）$\boldsymbol{\Phi}$と、重みベクトル（または回帰係数ベクトル）$\boldsymbol{w}$を用いて下記のようにあらわす。

\boldsymbol{y}=\boldsymbol{\Phi}(\boldsymbol{x})\boldsymbol{w}

すなわち、

\begin{equation*}
    \begin{pmatrix}
    y_1 \\
    \vdots \\
    y_n
    \end{pmatrix}
    =
    \begin{pmatrix}
    \phi_{1}(x_1) & \cdots & \phi_{m}(x_1) \\
     & \vdots &  \\
    \phi_{1}(x_n) & \cdots & \phi_{m}(x_n) \\
    \end{pmatrix}
    \begin{pmatrix}
    w_1 \\
    \vdots \\
    w_n
    \end{pmatrix}
\end{equation*}

$\phi_i$: 基底関数

重回帰分析の場合は、任意の$i,j$に対して$x_i = \phi_j(x_i)$とする上で$\boldsymbol{w}を求める$

ここで、$\boldsymbol{w}$が、平均ベクトル（期待値）$\boldsymbol{0}$、分散共分散行列$\lambda^2 \boldsymbol{I}$となるガウス分布に従う（数式で書くと、$\boldsymbol{w} \sim \mathcal{N}(\boldsymbol{0},\lambda^2 \boldsymbol{I})$）と定義する。
このとき、ガウス分布に従う確率変数の線形変換はガウス分布に従うこと、$\boldsymbol{\Phi}$が定数であることから、$\boldsymbol{y}$もガウス分布に従うといえる。

出力がどのようなガウス分布に従うか

$\boldsymbol{y}$の平均（期待値）は、

\mathbb{E}[\boldsymbol{y}]=\boldsymbol{\Phi}\mathbb{E}[\boldsymbol{w}]=\boldsymbol{0}

$\boldsymbol{y}$の分散共分散行列$\boldsymbol{\Sigma_y}$は、

\mathbb{E}[(\boldsymbol{y}-\mathbb{E}[\boldsymbol{y}])(\boldsymbol{y}-\mathbb{E}[\boldsymbol{y}])^{T}]=\mathbb{E}[\boldsymbol{y}\boldsymbol{y}^T]
=\boldsymbol{\Phi}\mathbb{E}[\boldsymbol{w}\boldsymbol{w}^T]\boldsymbol{\Phi}^{T}=\boldsymbol{\Phi}\lambda^2 \boldsymbol{I}\boldsymbol{\Phi}^{T}=\lambda^2\boldsymbol{\Phi}\boldsymbol{\Phi}^{T}

となるため、$\boldsymbol{y} \sim \mathcal{N}(\boldsymbol{0},\lambda^2\boldsymbol{\Phi}\boldsymbol{\Phi}^{T})$となる。したがって、$\boldsymbol{y}$がどのようなガウス分布に従うかを求める際、$\boldsymbol{w}$を求める必要がなく、$\boldsymbol{\Sigma_y}$を求めるのみでよいことがわかる。
ただし、一般的に$\boldsymbol{\Sigma_y} = \lambda^2\boldsymbol{\Phi}\boldsymbol{\Phi}^{T}$を直接求めることは計算コストが高い。

どのように出力の分散共分散行列を求めるか

分散共分散行列の特徴をとらえた別の関数でおきかえることで計算を軽くする手法がとられる（この手法をカーネルトリックといい、置き換える関数$k$をカーネル関数、分散共分散行列の要素をカーネル関数で置き換えた行列$\boldsymbol{K}$をカーネル行列またはグラム行列という）。
最終的にやりたいことは、得られているデータ$\boldsymbol{x}, \boldsymbol{y}$をもとに未知のデータ$x^{\ast}$に対する$y^{\ast}$を予測することであるため、任意の$x$に似た$x^{\ast}$に対しては、$y$に似た$y^{\ast}$を予測したい。
したがって、カーネル関数が再現すべき特徴は、$x_{i}$と$x_{j}$に相関が大きい場合、$y_{i}$と$y_{j}$に相関が大きくなるようにする。すなわち分散共分散行列の要素が大きくなるようにすることである。
このようなカーネル関数はいくつかあるが、題材にしたコードで用いているのは、RBF（Radial Basis Function、放射基底関数）という下記の関数である。

k(r)=\sigma^2 \exp{\bigg(-\frac{1}{2}r^2\bigg)}

ここで$r$はユークリッド空間上の距離であり、例えば$x_{i}$と$x_{j}$の距離$r_{i,j}$は$||x_i - x_j||$である。

以上により、すでに得られている$\boldsymbol{y}$が従う分布をモデル化できることがわかった。次に、未知のデータ$y^{\ast}$を含む場合どのように分布をモデル化するかを考える。

RBFカーネルは、二乗指数カーネル、べき乗二次カーネル、ガウスカーネルとも呼ばれるらしい

その他の関数や、分布のグラフについては下記の記事で紹介されています。

未知のデータに対する予測分布の求め方

未知データ$x^{\ast}, y^{\ast}$を含んだ出力ベクトルがどのような分布に従うかを書くと、下記のようになる。

\begin{equation*}
    \begin{pmatrix}
    y_1 \\
    \vdots \\
    y_n \\
    y^{\ast}
    \end{pmatrix}
    \sim
    \mathcal{N}
    \begin{pmatrix}
        \boldsymbol{0},
        \begin{pmatrix}
        k(r_{1,1}) & \cdots & k(r_{1,N}) & k(r_{1,\ast}) \\
        \vdots & \ddots &  \vdots & \vdots\\
        k(r_{N,1}) & \cdots & k(r_{N,N}) & k(r_{N,\ast}) \\
        k(r_{\ast,1}) & \cdots & k(r_{\ast,N}) & k(r_{\ast,\ast}) \\
        \end{pmatrix}
    \end{pmatrix}
\end{equation*}

ここで、$k_{\ast} = (k(r_{1,\ast}), k(r_{2,\ast}), \cdots, k(r_{N,\ast}))$とすると、最終的に得たい確率分布$p(y^\ast|x^\ast, \boldsymbol{x}, \boldsymbol{y})$は、下記のようにあらわすことができる（らしい）。

p(y^\ast|x^\ast, \boldsymbol{x}, \boldsymbol{y}) = \mathcal{N}(\boldsymbol{k}_{\ast}^{T}\boldsymbol{K}^{-1}\boldsymbol{y}, k(r_{\ast, \ast}) - \boldsymbol{k}_{\ast}^{T}\boldsymbol{K}^{-1}\boldsymbol{k}_{\ast})

以上のように、未知データに関する確率分布を得ている。

参考