ガウス過程 〜実践編〜

はじめに

前回、ガウス過程の理論を解説した。今回はPyMC3を用いて実際に計算を行ってみる。ガウス過程の計算では、パラメータを含むカーネルをあらかじめ「主観的」に設定し、観測データを用いてこのパラメータを推定する。パラメータを推定する手法として、最尤推定法、事後確率最大化法、ベイズ推定法がある。今回はベイズ推定を行う。

ベイズ推定による定式化

回帰問題を考える。観測値を(X,Y)、パラメータを\thetaとして、これらの同時確率分布p(X,Y,\theta)から次式を得る。

(1)   \begin{equation*} p(\theta|X,Y)=\frac{p(Y|X,\theta)\;p(\theta)}{p(Y|X)} \end{equation*}

適当な尤度p(Y|X,\theta)と事前確率p(\theta)を設定し、事後確率p(\theta|X,Y)を求めることがベイズ推定の目的である。いま

(2)   \begin{eqnarray*} X&=&(x_1,\cdots,x_N), \;\;x_n \in \mathbb{R}^M\\ Y&=&(y_1,\cdots,y_N), \;\;\;y_n \in \mathbb{R} \end{eqnarray*}

とおき、パラメータ\thetaをカーネルに関わる部分\theta_Kとそうでない部分\theta_Lの2つに分離する。さらに、XとYの間に次式を仮定する。

(3)   \begin{eqnarray*} y_n&=&f(x_n)+\epsilon_n \\ \epsilon_n&\sim&\mathcal{N}(\epsilon_n|0,\sigma^2) \end{eqnarray*}

ここで、\mathcal{N}(\epsilon_n|0,\sigma^2)は平均0、分散\sigma^2の正規分布を表す。各観測値は独立同分布から生成されると仮定すると、尤度は

(4)   \begin{eqnarray*} p(Y|X,\theta_L) &=&\prod_{n=1}^{N} p(y_n|x_n,\theta_L) \\ &=&\prod_{n=1}^{N} \mathcal{N}(y_n|f(x_n),\sigma^2) \end{eqnarray*}

と書くことができる。ここで、\theta_L=(f(\cdot),\sigma)とした。\sigmaの事前分布としてp(\sigma)を、関数f(\cdot)の事前分布として次のガウス過程を仮定する。

(5)   \begin{equation*} p(f|\theta_K)=\mathcal{N}(f|0_N,K_{NN}^{\theta_K}) \end{equation*}

ここで、f=\left(f(x_1),\cdots,f(x_N)\right)とした。0_NはN次元のゼロベクトル、K_{NN}^{\theta_K}はパラメータ\theta_Kを持つN\times Nの共分散行列、その成分(K_{NN}^{\theta_K})_{n,m}=k_{NN}^{\theta_K}(x_n,x_m)がカーネル関数となる。さらに、\theta_Kの事前分布としてp(\theta_K)も導入する。ここまでの結果をまとめると、式(1)は次のようになる。

(6)   \begin{equation*} p(f,\sigma,\theta_K|X,Y)=\frac{p(Y|X,\sigma,f)\;p(\sigma)\;p(f|\theta_K)\;p(\theta_K)}{p(Y|X)} \end{equation*}

fの事前分布p(f|\theta_K)がガウス過程であることに注意する。上の式の右辺をベイズ推定で評価し、事後確率p(f|X,Y),p(\sigma|X,Y),p(\theta_K|X,Y)を求める。これらを用いると、未知の値x_{\star}に対するf_\star=f(x_\star)の確率分布を次式により計算することができる。

(7)   \begin{equation*} p(f_\star|X,Y)=\int df\;d\theta_K\;p(f_\star|f,\theta_K)\;p(f|X,Y)\;p(\theta_K|X,Y) \end{equation*}

右辺の積分内の最初の項p(f_\star|f,\theta_K)はp(f(x_\star)|f(x_1),\cdots,f(x_N),\theta_K)を表し、式(5)から解析的に求めることができる条件付き確率分布である。上式からf_{\star}の分布が分かれば、未知の値x_{\star}に対するy_\starの値は次式から決定することができる。

(8)   \begin{eqnarray*} f_\star&\sim&p(f_\star|X,Y) \\ \sigma&\sim&p(\sigma|X,Y) \\ \epsilon_\star&\sim&\mathcal{N}(\epsilon_\star|0,\sigma^2) \\ y_\star&=&f_\star+\epsilon_\star \end{eqnarray*}

以上で定式化できたので、実際に計算を行う。

ソースの場所

今回のソースはここのsample_pymc3_latent.ipynbである。

データの生成

人工的に生成した以下のデータを使う。

計算式は以下の通り。

(9)   \begin{eqnarray*} f(x)&=&\cos{\left(\frac{x^2}{2}\right)}+\frac{x^2}{10}\\ y&=&f(x)+\epsilon\\ \epsilon&\sim&\mathcal{N}(\epsilon|0,\sigma^2)\\ \sigma&=&0.2 \end{eqnarray*}

PyMC3による計算

ベイズ推定のコードは以下の通り。

  • 2行目から4行目まででカーネルを定義する。\ellと\etaはカーネルパラメータ\theta_Kに相当する。このパラメータに対しても事前分布（ガンマ分布、半コーシ分布）を設定していることに注意する。
  • 6行目で式(5)のガウス過程p(f|\theta_K)を定義する。
  • 10行目から12行目までで式(4)の尤度を定義する。上の定式化では尤度としてガウス分布を用いたが、コードではスチューデントのt分布を用いた。外れ値に強いためである。この分布のパラメータ\sigma,\nuに対しても事前分布を定義する。パラメータ\sigma,\nuが\theta_Lに相当する。
  • 14行目でMCMCを実行する。

上の計算のあと次のコードで未知データx_\starに対するf_\star=f(x_\star)の分布を計算する。

  • 1行目と2行目で未知データx_\starを作成する。
  • 5行目で条件付き確率p(f_\star|f,\theta_K)が設定される。
  • 6行目で、この条件付き確率とtrace内に格納されている各種事後確率を用いて、式(7)で与えたp(f_\star|X,Y)が計算される。samplesは生成する曲線f_{\star}の数。結果は以下の通り

    生成された1000本の曲線を赤で示した。赤の濃さは曲線の密度を表し、曲線の揺れ幅が予測の不確かさを表す。データの密度が高い場所では不確かさは小さく、データの密度が低い場所では不確かさは大きくなる。

    まとめ

    ガウス過程による回帰の計算は、PyMC3を用いれば容易に実行することができることを示した。また、式とコードを明確に関連付けた。理論式と対応づけることで理解を深めることができる。
    ガウス過程で重要となるのは、カーネルとして何を設定すべきかという点である。上のコードではPyMC3のチュートリアルで紹介されているカーネルをそのまま用いた。チュートリアルを読むと、どのカーネルを選択すべきかのおおよその目安は分かる。カーネルの選択についての解説は機会を改めたい。

    参照文献

    1. Pythonによるベイズ統計モデリング　良書ですが、理論式との対応付けが明確でないので何を計算しているのか分からなくなることがある。
    2. PyMC3
Kumada Seiya

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

最近の記事

  • 関連記事
  • おすすめ記事
  • 特集記事

アーカイブ

カテゴリー

PAGE TOP