(導出)ガウス過程回帰を理解したい-Gaussian Process Regression

はじめに

最近、GPLVMに関係する論文を読んでたんですが、そもそもガウス過程について全く理解していなかったので勉強がてら記事にしてまとめてみました。

本記事ではガウス過程回帰の導出について書きました。

本記事の主な流れ

ガウス過程回帰を導出するに至って次の流れで求めていきま。

ガウス過程回帰への流れ

線形回帰モデルに対してベイズ推定の枠組みを導入
ベイズ推定の枠組みにおけるパラメータ$w$の事後分布を推定する
新規のデータに対して得られる予測分布を算出する
得られた予測分布に対してカーネルの枠組みを導入する
カーネル関数を用いた予測分布が得られる

線形回帰モデル

線形回帰モデルは基底関数における特徴ベクトルの数を増やすことで複雑な関数も表現できました。

$$f(x)=w_1\varphi_2(x)+w_2\varphi_2(x)+\dots+w_m\varphi_m(x)=\sum_{i}^{m}w_{i}\varphi_{i}(x)$$

上記の線形回帰のモデルの式からわかるように1次元の$x$に対して$N$個の特徴ベクトルで考えるとき、それに対応する基底関数の数も$N$個必要になります。

このとき、$x$の次元が2次元の場合は基底関数は$N^2$個必要になり、3次元の場合は$N^3$個必要になります。

このように$x$の次元が大きくなると必要な基底関数は指数関数的に増大します。

つまり、線形回帰のモデルは$x$の次元が低いときにのみ有効なのです。

: 線形回帰モデルのDualとカーネル関数について

この記事は古川研究室 Advent_calendar 19日目の記事です。本記事は古川研究室の学生 ...

続きを見る

ベイズ線形回帰

線形回帰では予測したい出力$y_n $y_{n}$ 、入力$x_n $x_{n}$ と重みベクトル$ $w$ によって決定しました。このとき、$x_n$と$y_n$は観測された変数ですが、重みベクトル$ $w$ は未知の変数です。線形回帰では未知の重みベクトル$w$を推定するのですが、ここにベイズ推定の枠組みを導入します。

ベイズ線形回帰では、推定対象となるパラメータ$w$に確率分布を仮定します。つまり、パラメータ$w$はある確率分布に従う確率変数であると考えます。

問題設定

Given

学習データセット

$$X=(x_1,y_1),\dots,(x_n,y_n)$$

パラメータ$w$の事前分布

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

Estimate

パラメータ$w$の事後分布$p(\mathbf{w}|X)$

$$p(\mathbf{w}|X)=\frac{p(X|\mathbf{w})p(\mathbf{w})}{p(X)}$$

新規の入力$x^{\ast}$に対する予測確率分布$p(y^{\ast}|x^{\ast},X)$

$$p(y^{\ast}|x^{\ast},X)=\int p(y^{\ast}|x^{\ast},\mathbf{w})p(\mathbf{w}|X)d \mathbf{w}$$

対数尤度を求める

ここから先は最尤推定やガウス分布などの知識が必要になります。

必要な知識の概要を簡単に書きました

必要知識の簡単な復習はこちら

確率密度関数

連続型確率変数$X$がある値$x$をとる確率密度を関数$f(x)$とすると、$f(x)$を「確率密度関数」と言います。

性質としては、連続型確率変数$A$の取りうる値の下限値を$ $A$$ ，上限値を$B$とおくとき、下限がないときは$A=-\infty$上限がないときは$B=-\infty$になります。$f(x)$が確率密度関数のとき$\int_A^Bf(x)dx=1$が成り立つ。

正規分布(ガウス分布)

確率関数$X$が平均$\mu$、分散$\sigma^2$ の正規分布に従うとき、その確率密度関数は以下のように表すことができる。

$$f(X, \mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma}}\exp(-\dfrac{(X-\mu)^ 2}{2\sigma^ 2})$$

主な性質は以下のものがある。

平均値と最頻値と中央値が一致
$x$ 軸が漸近線
平均値を中心にして左右対称
分散（標準偏差）が大きくなると左右に広がって平らになる。分散（標準偏差）が小さくなると、よりとんがった形になる。

最尤推定

パラメータ$\theta$に従う分布の密度関数を$f(x;\theta)$とする。尤度関数を$L(\theta;x)=f(\theta;x)$とすると、$L(\theta;x)$を最大にするような推定量$\theta$を最尤推定量という。

「最尤」の意味は文字通り「一番尤もらしい」という意味になります。つまり、最尤法によるパラメータ推定とは、データから一番尤もらしいパラメータを推定する手法のことです。この手法によって推定されたパラメータを最尤推定量と呼びます。

$L(\theta;x)$が最大となるような$\theta$を考えるときに、微分を行いますが、そのまま微分をすると非常に複雑な計算を要するため尤度関数の対数を取った対数尤度$log L(\theta;x)$を考えます。なぜこんなことするのかというと、。対数を取ると、積ではなく、和の形になりますので、微分式を解くのにも計算が楽となるのと、ある関数の対数の最大化はその関数の最大化と同値なので、対数に置き換えても影響はないからです

入力とパラメータ$w$がわかっているときの観測される確率分布は次のガウス分布$N(y_n|f(x_n),\beta^{^-1})$に従います(尤度関数)

分散である$\beta^{^-1}$はノイズです。

$${\begin{eqnarray}p\left(y_n|x_n,\mathbf{w}\right)&=&N\left(y_n|f(x_n),\beta^{^-1}\right)\\&=&\sqrt\frac{\beta}{{2\pi}}\exp\left(-\frac{\beta}{2}(y_n-f(x_n))^{2}\right)\\&=&\sqrt\frac{\beta}{{2\pi}}\exp\left(-\frac{\beta}{2}\left(y_n-\mathbf{w}^{T}\boldsymbol{\varphi}\left(x_n\right)\right)^{2}\right)\end{eqnarray}}$$

この対数尤度を取ると

$$\ln{p(y_n|x_n,\mathbf{w})}=[-\frac{\beta}{2}(y_n-\mathbf{w}^{T}\boldsymbol{\varphi}(x_n))^2]+const$$

$c o n s t$ は$w$に関して定数であることを表しています。

これらは特定の$n$番目のデータにおける対数尤度でした。

これを全てのデータについて考えますと尤度関数と対数尤度はそれぞれ次のようになります。

$$p(X|\mathbf{w})=\prod_{n=1}^{N} p(y_n|x_n, \mathbf{w})$$

$${\begin{eqnarray}\ln{p(X|\mathbf{w})}&=&\sum_{n=1}^{N}[-\frac{\beta}{2}(y_n-\mathbf{w}^{T}\boldsymbol{\varphi}(x_n))^2]+const\\&=&-\frac{\beta}{2}||\mathbf{w}^{T}\Phi - y||^2+const \end{eqnarray}}$$

パラメータ$w$の事後分布を求める

ベイズの定理の復習はこちら

同時分布$P(A,B)$において、$B$がわかったときの$A$の確率分布(条件付き確率分布)、$A$がわかったときの$B$の確率分布(条件付き確率分布)は次のように表すことができます。

$$P(A|B)=\frac{P(A\cap B)}{P(B)}\,\,\,\,\,,\,\,\,\,\,P(B|A)=\frac{P(A\cap B)}{P(A)}$$

これらの式から次のことが導けます。

$$P(x|y)=\frac{P(y|x)P(x)}{P(y)}$$

P(x)は事前確率、P(x|y)は事後確率、P(y|x)は尤度（ゆうど）、P(y)は周辺尤度と呼ばれます。

$$事後分布=\frac{尤度 \times 事前分布}{周辺分布}$$

このとき、$A$の取りうる値が複数あるとすると、$A_i(i=1,2,\dots,n)$が互いに独立であるという条件の下では

$$P(B)=\sum_{i=1}^{n}P(A_i \cap B)=\sum_{i=1}^{n}P(B|A_i)P(A_i)$$

よって

$$P(A_i|B)=\frac{P(B|A_i)P(A_i)}{P(B)}=\frac{P(B|A_i)P(A_i)}{\sum_{i=1}^{n}P(A_i)P(B|A_i)}$$

これは離散型のベイズの定理ですが$\sum$を$\int$にすれば連続型になります

連続型ベイズの定理の分母（周辺尤度）について着目すると、積分することから周辺尤度は定数になります。つまり、比例の記号を用いて

$$事後分布 \propto 尤度 \times事前分布$$

が成り立ちます。

共役事前分布についてはこちら

ベイズ推定は式の形から、事前分布が複雑であると、その事後分布の計算も複雑なります。

そこで、事後分布の形が簡単になるような、尤度関数と事前分布を選択したいところです。

一般的には、事前分布としてある確率分布を選び、尤度関数にかけたら、その事後分布も同じ確率分布の形で表せられるような確率分布が存在していることが知られています。

ある確率分布を選んで事前分布としてベイズ推定を行った時、その事後分布もその確率分布の形で表されるとき、その事前分布（確率分布）のことを共役事前分布と言います。

したがって、共役事前分布とは、ベイズ推定において「事後分布」と「事前分布」が同じ形になるように選んだ「事前分布」のことです。

共役事前分布は尤度関数の形に応じて決められるので、代表的なものをいくつかあげると、次の表のようになります。

事前分布	尤度	事後分布
ガンマ分布	ポアソン分布	ガンマ分布
逆ガンマ分布	正規分布	逆ガンマ分布
ベータ分布	二項分布	ベータ分布
正規分布	正規分布	正規分布

パラメータ$w$の事後分布を求めます。

ベイズの定理より$p(\mathbf{w}|X)=\frac{p(X|\mathbf{w})p({\mathbf{w}})}{p({\mathbf{X}})}$

対数をとって

$${\begin{eqnarray} \ln p(\mathbf{w}|X)&=&\ln{p(X|\mathbf{w})}+\ln p({\mathbf{w}})-\ln p({\mathbf{X}})\\&=&-\frac{\beta}{2}||y-\mathbf{w}^{T}\Phi||^2-\frac{\alpha}{2}||\mathbf{w}||^2+const\\&=&-\frac{\beta}{2}(\mathbf{w}^{T}\Phi - y)^T(\mathbf{w}^{T}\Phi - y)-\frac{\alpha}{2}\mathbf{w}^T\mathbf{w}+const\\&=&-\frac{1}{2}(\beta\mathbf{w}^T\Phi^T\Phi\mathbf{w}+\alpha\mathbf{w}^T\mathbf{w})+(\beta\mathbf{y}^T\Phi)\mathbf{w}+const\\&=&-\frac{1}{2}\mathbf{w}^T(\beta\Phi^T\Phi+\alpha\mathbf{I})\mathbf{w}+(\beta\mathbf{y}^T\Phi)\mathbf{w}+const\end{eqnarray}}$$

constはパラメータ$w$と無関係な値になっています。

パラメータ$w$の事後分布が求まりました。ここで、尤度（関数）と事前分布の関係が共役な関係になるように事前分布を選ぶと事後分布が事前分布と同じ分布になるので，事後分布は多変量ガウス分布$N(\mathbf{w}|\mathbf{\mu},S)$になります。

つまり$p(\mathbf{w}|X)=N(\mathbf{w}|\mathbf{\mu},S)$です。

$n$変数の多変量正規分布は${f(\vec{x}) = \frac{1}{\sqrt{(2\pi)^n |\sum|}}exp \left \{-\frac{1}{2}{}^t (\vec{x}-\vec{\mu}) {\sum}^{-1} (\vec{x}-\vec{\mu}) \right \} }$の形で表すことができるので次の式が導けます。

$${\begin{eqnarray}\ln p(\mathbf{w}|X)&=&-\frac{1}{2}(\mathbf{w}-\mathbf{\mu})^TS^{-1}(\mathbf{w}-\mathbf{\mu})+const\\&=&-\frac{1}{2}\mathbf{w}^T S^{-1}\mathbf{w}+\mathbf{\mu} S^{-1}\mathbf{w}+const\end{eqnarray}}$$

2つの$p(\mathbf{w}|X)$の式を比較してみます。

$$\ln p(\mathbf{w}|X)=-\frac{1}{2}\mathbf{w}^T(\beta\Phi^T\Phi+\alpha\mathbf{I})\mathbf{w}+(\beta\mathbf{y}^T\Phi)\mathbf{w}+const$$

$$\ln p(\mathbf{w}|X)=-\frac{1}{2}\mathbf{w}^T S^{-1}\mathbf{w}+\mathbf{\mu} S^{-1}\mathbf{w}+const$$

したがって次のことが成り立ちます。

$$S^{-1}=\beta\Phi^T\Phi+\alpha\mathbf{I}$$

$$\mathbf{\mu}S^{-1}=\beta\mathbf{y}^T\Phi$$

Sについて解くと

$$\mathbf{S}=(\beta\Phi^{T}\Phi+\alpha \mathbf{I})^{-1}$$

$\mu$について解くと

$${\begin{eqnarray} \mathbf{\mu}^T\mathbf{S}^{-1}&=&\beta\mathbf{y}^T\Phi\\\mathbf{\mu}^T&=&\beta\mathbf{y}^T\Phi S\\ \mathbf{\mu}^T&=&\beta\mathbf{y}^T\Phi (\beta\Phi^{T}\Phi+\alpha \mathbf{I})^{-1}\\ \mathbf{\mu}^T&=&\mathbf y^T\Phi(\Phi^T\Phi+\lambda \mathbf I)^{-1}\\ \end{eqnarray}}$$

ゆえに

$${\begin{eqnarray} \mathbf{\mathbf{\mu}} &=&(y^T\Phi(\Phi^T\Phi+\lambda \mathbf I)^{-1})^T\\ &=&((\Phi^{T}\Phi+\lambda \mathbf{I})^{-1})^T(\mathbf y^T\Phi)^T\\ &=&(\Phi^{T}\Phi+\lambda \mathbf{I})^{-1}\Phi^{T}\mathbf{y} \end{eqnarray} }$$

このとき$(\Phi^{T}\Phi+\lambda \mathbf{I})^{-1}$は対称行列、$\lambda$は$\lambda=\frac{\alpha}{\beta}$です。

まとめますと

事後分布が多変量ガウス分布$N(\mathbf{w}|\mathbf{\mu},S)$のとき，その$\mu$と$S$は，

$${\begin{eqnarray} \mathbf{S}&=&(\beta\Phi^{T}\Phi+\alpha \mathbf{I})^{-1}\\ \mathbf{\mu}&=&(\Phi^{T}\Phi+\lambda \mathbf{I})^{-1}\Phi^{T}\mathbf{y} \end{eqnarray} }$$

新規の入力$x^{\ast}$に対する予測確率分布$p(y^{\ast}|x^{\ast},X)$を求める

推定された事後分布から、新規の入力$x^{\ast}$に対する予測確率分布$p(y^{\ast}|x^{\ast},X)$を考えていきます。

$${\begin{eqnarray}p(y^{\ast}|x^{\ast},X)&=&\int p(y^{\ast}|x^{\ast},\mathbf{w})p(\mathbf{w}|X)d \mathbf{w}\\&=&\int N\left(y_n|f(x_n),\beta^{^-1}\right)N(\mathbf{w}|\mathbf{\mu},S)d \mathbf{w}\end{eqnarray} }$$

この予測分布の式は、ガウス分布同士の畳み込み積分になっています。ここで，ガウス分布同士の畳み込み積分を求める公式を使って計算していくと次の式が得られます。非常に難しい計算で膨大になるので省略して結果だけのせます。(ここはいつかLatexで書いて証明したい)

ガウス過程回帰の式

$$p(y^{\ast}|\mathbf{x}^{\ast},X,Y)=\mathcal{N}(\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{Y},\mathbf{k}_{\ast,\ast}-\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{k}_{\ast})$$

式の左辺の$Y$は今まで省略していました。

$\mathbf{k}_{\ast}$は$X$と$\mathbf{x}^{\ast}$のカーネル、$\mathbf{K}$は$X$と$X$のカーネル、$\mathbf{K}_{\ast,\ast}$は$\mathbf{x}^{\ast}$と$\mathbf{x}^{\ast}$のカーネルです．

終わりに

長い記事になりましたが最終的にガウス過程回帰の式を導出することができました。

実装やガウスの性質については別記事に書きたいと思います。

ここまで読んでくださってありがとうございました。

見にくいところや間違っているところの編集リクエストもお待ちしています。