はじめに
ガウス過程について一通り勉強した僕は一点だけもやもやしていた点があった。
それは予測分布の導出方法である。
ガウス過程の説明において、予測分布は以下のように導出されることが多い。
「$\mathbf{y}$と$\mathbf{y}_{*}$の同時分布はガウス過程に従うので、以下のように表せます。」
\left(
\begin{matrix}
\mathbf{y} \\
\mathbf{y}_{*}
\end{matrix}
\right)\sim\mathcal{N}\left(\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),\left(
\begin{matrix}
\mathbf{K} &\mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)\right)
「ガウス分布の公式によって、予測分布は以下のようになります」
p(\mathbf{y}_{*}|\mathbf{y})=\mathcal{N}\left(\mathbf{k}_{*}^{T}\mathbf{K}^{-1}\mathbf{y},\mathbf{k}_{**}-\mathbf{k}_{*}^{T}\mathbf{K}^{-1}\mathbf{k}_{*}\right)
だが、よくあるベイズ的には
- 事前分布と確率モデル（尤度関数）を決める。
- ベイズの定理に従って事後分布を求める。
- 事後分布で確率モデルの期待値を取ることで予測分布を求める。
という手続きを踏んで初めて、予測分布に到達する。
上記のガウス過程の説明だと、いきなり予測分布に飛ぶので少し怖い。
そこで、この記事では、ガウス過程回帰における事前分布・事後分布・予測分布の導出についてきちんと整理したいと思う。
基本的なガウス過程の考え方に関しては他の記事や文献に譲る。
全体の同時分布
まず、無心になってこの最強の同時分布を計算する。
p(\mathbf{f},\mathbf{f}_{*},\mathbf{y},\mathbf{y}_{*}|\mathbf{X},\mathbf{X}_{*})
これを計算しておくと後々便利なのである。
$\mathbf f$は入力$\bf X$（データが観測された点）における、ノイズが乗る前のデータベクトル（観測できない）。
$\mathbf f_ *$は入力$\mathbf X_ *$（基本的にデータがない点）における、ノイズが乗る前のデータベクトル（観測できない）。
$\mathbf y$は入力$\bf X$（データが観測された点）において観測されたデータ。
$\mathbf y_ *$は入力$\mathbf X_ *$（基本的にデータがない点）における推定値。
以降は$p(\mathbf f,\mathbf f _ *,\mathbf y,\mathbf y _*)$という省略した形で表す。
\begin{align}
p(\mathbf{f},\mathbf{f} _ *,\mathbf{y},\mathbf{y} _*)&= p(\mathbf{y},\mathbf{y} _*|\mathbf{f},\mathbf{f}_{*})p(\mathbf{f},\mathbf{f}_{*})\\
&=\mathcal{N}\left( \left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*}
\end{matrix}
\right),\sigma^2\mathbf{I}\right)\mathcal{N}\left( \left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),\mathbf{K}'\right) \\
&\propto\exp\left(-\frac{1}{2}\left(
\begin{matrix}
\mathbf{y}-\mathbf{f} \\
\mathbf{y}_{*}-\mathbf{f}_{*}
\end{matrix}
\right)^{T}(\sigma^2\mathbf{I})^{-1}\left(
\begin{matrix}
\mathbf{y}-\mathbf{f} \\
\mathbf{y}_{*}-\mathbf{f}_{*}
\end{matrix}
\right)\right)\\
&\times \exp\left(-\frac{1}{2}\left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*}
\end{matrix}
\right)^{T}(\mathbf{K}')^{-1}\left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*}
\end{matrix}
\right)\right)\\
&=\exp\left(-\frac{1}{2}\left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*} \\
\mathbf{y} \\
\mathbf{y}_{*}
\end{matrix}
\right)^{T}\left(
\begin{matrix}
(\mathbf{K}')^{-1}+(\sigma^2\mathbf{I})^{-1}&-(\sigma^2\mathbf{I})^{-1}\\
-(\sigma^2\mathbf{I})^{-1}&(\sigma^2\mathbf{I})^{-1}
\end{matrix}
\right)\left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*} \\
\mathbf{y} \\
\mathbf{y}_{*}
\end{matrix}
\right)\right)
\end{align}
ここで
\mathbf{K}'=\left(
\begin{matrix}
\mathbf{K}&\mathbf{k}_{*}\\
\mathbf{k}_{*}^{T}&\mathbf{k}_{**}
\end{matrix}
\right)
上記の式の真ん中の行列の逆行列が、同時分布の分散共分散行列になる。
PRMLの2.3.3章を参考にして計算すると分散共分散行列$\Sigma$は
\begin{align}
\Sigma&=\left(
\begin{matrix}
\mathbf{K}'&\mathbf{K}'\\
\mathbf{K}'&\mathbf{K}+\sigma^2\mathbf{I}'
\end{matrix}
\right)\\
&=\left(
\begin{matrix}
\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)&\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)\\
\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)&\left(
\begin{matrix}
\mathbf{K}+\sigma^2\mathbf{I} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}+\sigma^2\mathbf{I}
\end{matrix}
\right)
\end{matrix}
\right)
\end{align}
と表せる。つまり、
\begin{align}
\left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*} \\
\mathbf{y} \\
\mathbf{y}_{*}
\end{matrix}
\right)\sim
\mathcal{N}\left(
\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0} \\
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),\left(
\begin{matrix}
\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)&\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)\\
\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)&\left(
\begin{matrix}
\mathbf{K}+\sigma^2\mathbf{I} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}+\sigma^2\mathbf{I}
\end{matrix}
\right)
\end{matrix}
\right)\right)
\end{align}
実は、これを周辺化することで3つの分布を計算できてしまう。
事前分布
ガウス過程回帰における事前分布は、いわゆる「関数を生成する確率分布」である。
だが実際は$N$次元の実数ベクトル$\mathbf{f}$を生成しているだけの、ごく普通の多変量ガウス分布である。
ここに「関数」という概念を間に無理やり入れ込んで
入力$\mathbf{X}$が与えられる $\to$ 事前分布の形が決定する $\to$ $\mathbf{f}$がサンプリングされる
という流れを、
入力$\mathbf{X}$が与えられる $\to$ 事前分布の形が決定する $\to$ $\mathbf{X}$を$\mathbf{f}$に写す関数$f$が生成される $\to$ $\mathbf{X}$が$\mathbf{f}$に写される
とわざわざ解釈しているという話である（たぶん）。
ちなみに$\mathbf{f}=(f_1,f_2,\cdots,f_N)=(f(\mathbf x_1),f(\mathbf x_2),\cdots,f(\mathbf x_N))$と書ける。
さて、その事前分布は
p(\mathbf{f}|\mathbf{X})
p(\mathbf{f}_{*}|\mathbf{X}_{*})
p(\mathbf{f},\mathbf{f}_{*}|\mathbf{X},\mathbf{X}_{*})
という3パターンの形式で書ける。
形は違うが、それは入力が違うだけでどれも本質的には同じことを表している。
どれも平均が$\mathbf{0}$で、要素がカーネルの分散共分散行列の正規分布なので
\mathbf{f}\sim\mathcal{N}\left(\mathbf{0},\mathbf{K}\right)
\mathbf{f}_{*}\sim\mathcal{N}\left(\mathbf{0},\mathbf{k}_{**}\right)
\begin{align}
\left(
\begin{matrix}
\mathbf{f} \\
\mathbf{f}_{*}
\end{matrix}
\right)\sim
\mathcal{N}\left(
\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),
\left(
\begin{matrix}
\mathbf{K} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)\right)
\end{align}
となる。
ちなみに先ほどの同時分布を周辺化してみても、同じ結果が得られることがわかる。
事後分布
ガウス過程回帰における事後分布もまた、関数を生成する確率分布である。
ただしデータを使っているので、事前分布よりもデータに合わせた尤もらしい関数が出てくる。
事後分布は数式的には、
p(\mathbf{f}|\mathbf{y},\mathbf{X})
p(\mathbf{f}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})
p(\mathbf{f},\mathbf{f}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})
こちらについても先ほどの同時分布を周辺化すると、
\begin{align}
\left(
\begin{matrix}
\mathbf{y} \\
\mathbf{f}
\end{matrix}
\right)\sim
\mathcal{N}\left(
\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),
\left(
\begin{matrix}
\mathbf{K}+\sigma^2\mathbf{I} & \mathbf{K}\\
\mathbf{K} & \mathbf{K}
\end{matrix}
\right)\right)
\end{align}
\begin{align}
\left(
\begin{matrix}
\mathbf{y} \\
\mathbf{f}_{*}
\end{matrix}
\right)\sim
\mathcal{N}\left(
\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),
\left(
\begin{matrix}
\mathbf{K}+\sigma^2\mathbf{I} & \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)\right)
\end{align}
\begin{align}
\left(
\begin{matrix}
\mathbf{y} \\
\mathbf{f}\\
\mathbf{f}_{*}
\end{matrix}
\right)\sim
\mathcal{N}\left(
\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),
\left(
\begin{matrix}
\mathbf{K}+\sigma^2\mathbf{I} & \left(
\begin{matrix}\mathbf{K}&\mathbf{k}_{*} \end{matrix}
\right)\\
\left(
\begin{matrix}\mathbf{K}\\
\mathbf{k}_{*}^{T} \end{matrix}
\right)& \left(
\begin{matrix}
\mathbf{K}& \mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}
\end{matrix}
\right)
\end{matrix}
\right)\right)
\end{align}
が出てくる。
例の公式を使うことで、
\begin{align}
p(\mathbf{f}|\mathbf{y},\mathbf{X})=
\mathcal{N}\left(\mathbf{K}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y},\mathbf{K}-\mathbf{K}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{K}\right)
\end{align}
\begin{align}
p(\mathbf{f}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})=
\mathcal{N}\left(\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y},\mathbf{k}_{**}-\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{k}_{*}\right)
\end{align}
\begin{align}
p(\mathbf{f},\mathbf{f}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})=
\mathcal{N}\left(\left(
\begin{matrix}\mathbf{K}\\
\mathbf{k}_{*}^{T} \end{matrix}
\right)(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y},\mathbf{K}'-\left(
\begin{matrix}\mathbf{K}\\
\mathbf{k}_{*}^{T} \end{matrix}
\right)(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\left(
\begin{matrix}\mathbf{K}&\mathbf{k}_{*} \end{matrix}
\right)\right)
\end{align}
となって、事後分布の式が出てきた。
こちらもいずれも内容は同じで、入力が異なる3つのパターンをいちいち書いているだけである。
ベイズ的な計算法の場合、例えば一つ目の式ならば
\begin{align}
p(\mathbf{f}|\mathbf{y},\mathbf{X})=\frac{p(\mathbf{y}|\mathbf{f},\mathbf{X})p(\mathbf{f}|\mathbf{X})}{\int p(\mathbf{y}|\mathbf{f},\mathbf{X})p(\mathbf{f}|\mathbf{X})d\mathbf{f}}
\end{align}
というベイズの定理の式に従って事後分布を計算するのが王道だが、この場合は事前分布も尤度関数もガウス分布なので、どっちにしろ実際の計算は同時分布を計算して周辺化するという流れになる。
詳しくはやはりPRML2.3.3章を参考にしていただきたい。
（実は最初は共役事前分布の公式から事後分布を直接求めようとしたが結果がうまく合わず、次に愚直に指数の中身を計算しようとしたが墜落し、結局この方法にたどり着いた）
予測分布
とりあえず同様に周辺化すると
\left(
\begin{matrix}
\mathbf{y} \\
\mathbf{y}_{*}
\end{matrix}
\right)\sim\mathcal{N}\left(\left(
\begin{matrix}
\mathbf{0} \\
\mathbf{0}
\end{matrix}
\right),\left(
\begin{matrix}
\mathbf{K}+\sigma^2\mathbf{I} &\mathbf{k}_{*}\\
\mathbf{k}_{*}^{T} & \mathbf{k}_{**}+\sigma^2\mathbf{I}
\end{matrix}
\right)\right)
\begin{align}
p(\mathbf{y}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})=
\mathcal{N}\left(\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y},\mathbf{k}_{**}+\sigma^2\mathbf{I}-\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{k}_{*}\right)
\end{align}
があっさり出てくる。
一般的なベイズの流れに照らし合わせて考えると、
\begin{align}
p(\mathbf{y}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})=\int
\end{align}p(\mathbf{y}_{*}|\mathbf{f}_{*},\mathbf{X}_{*})p(\mathbf{f}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})d\mathbf{f}_{*}
というように、確率モデルを事後分布で期待値を取ることで予測分布は導き出される。
今回の場合は、
p(\mathbf{y}_{*}|\mathbf{f}_{*},\mathbf{X}_{*})=\mathcal{N}\left(\mathbf{f}_{*},\sigma^2\mathbf{I}\right)\\
p(\mathbf{f}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})=
\mathcal{N}\left(\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y},\mathbf{k}_{**}-\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{k}_{*}\right)
と、どちらもガウス分布なので、共役事前分布の性質を使えるので、
予測分布は事後分布の分散共分散行列に$\sigma^2\mathbf{I}$を足したものになり、結局
\begin{align}
p(\mathbf{y}_{*}|\mathbf{y},\mathbf{X},\mathbf{X}_{*})=
\mathcal{N}\left(\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{y},\mathbf{k}_{**}+\sigma^2\mathbf{I}-\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma^2\mathbf{I})^{-1}\mathbf{k}_{*}\right)
\end{align}
が導き出される。
参考文献
ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)
パターン認識と機械学習
ガウス過程の基礎と教師なし学習
ガウス過程回帰の基礎 - J-Stage
Conjugate prior