Python
PRML
MachineLearning

ガウス過程回帰 (PRML 第 6 章の図版再現)

More than 1 year has passed since last update.

Machine Learning Advent Calendar 2016の 2 日目です。

初日?は非常に気合入った内容でしたね。本日はあっさりした内容になっております。

今年から機械学習の勉強をはじめました。理解を深めるためにこの記事では PRML 第 6 章のガウス過程による回帰の、python での実装を試みています。

この記事の対象読者は、自分と同じように PRML を勉強している人となります。PRML 片手に読むとよいです。

この記事のゴールは、PRML の図 6.8 のような、右の方に行くに従って不確かさが大きくなる様子を出来る限り再現することです。

実装にあたっては @naoya_t さんの記事(カーネル法とガウス過程(PRML下巻 第6章))を特に参考にしました。(というか結果的にほとんど python に訳しただけになってしまいました)

楽しかったです ✌ ('ω' ✌ )三 ✌ ('ω') ✌ 三( ✌ 'ω') ✌


結果

とりあえず結果です。

完成図

こちらはオリジナルのグラフです

Screen Shot 2016-12-02 at 6.06.31 PM.png

終点は元のグラフより小さめですが、だいたい望みどおりのグラフが書けました。

ここで記事を終えてもいいのですが、一応以下で簡単に説明します。また、PRML を参照しやすいように、コードに対応する式の番号を書いています。


ガウス過程について

ガウス過程回帰では、予測時に訓練データ点の全部または一部を使います。これは、ベイズ線形回帰が訓練時にしか訓練データを使わないのと異なります。ベイズ線形回帰に比べて複雑な予測をすることができますが、過学習しやすくもあります。

ガウス過程回帰の計算量は訓練点の個数をNとすると$O(N^3)$かかります。一方でベイズ線形回帰では、基底関数の個数をMとすると$O(M^3)$です。訓練点は多くなりがちなため、ガウス過程回帰では計算に多くの時間を必要とします。


実装と解説

実装を簡単にするために入力は1次元としています。

行列計算は numpy を使い、pylab でグラフを描画しました。

python のバージョンは 2.7 です。


入力

同じ図版を再現するためには少なくとも同じ入力が必要です。入力に使われているのは PRML 上巻付録 A の人工データ(正弦関数)の右から 3 つの点を取り除いたものです。残念ながら公式からすでに消えているため、webarchive から探してきました。

https://web.archive.org/web/20110206181056/http://research.microsoft.com/en-us/um/people/cmbishop/prml/webdatasets/curvefitting.txt

コード中では、PRML の図と同じく最後の3点を取り除いて使っています。

  x_train = [

0.000000,
0.111111,
0.222222,
0.333333,
0.444444,
0.555556,
0.666667,
]

y_train = [
0.349486,
0.830839,
1.007332,
0.971507,
0.133066,
0.166823,
-0.848307,
]


カーネル関数

カーネル関数としてはガウスカーネル(6.23)や、

k(x, x') = exp(-||x-x'||^2/2\sigma^2)

def gaussian_kernel(sigma):

return lambda x1, x2: math.exp(-(x1 - x2) ** 2 / (2 * sigma ** 2))

2次形式の指数を取ったものに、定数と線形の項を加えたもの(6.63)などがあります。

k(x, x') = \theta_0exp \bigl\{- \frac{\theta_1}{2}||x-x'||^2 \bigl\} + \theta_2 + \theta_3 x x'

def regression_kernel(t0, t1, t2, t3):

return lambda x1, x2: t0 * math.exp(- (t1 / 2) * (x1 - x2) ** 2) + t2 + t3 * x1 * x2

ガウス過程回帰では後者が広く使われているようです。今回の実装では後者を使いました。


平均と共分散

ガウス過程によって、未知の点 $x_{N+1}$ の平均と共分散は、次の式(6.66, 6.67)で得られます。

m(x_{N+1}) = \mathbf{k^TC^{-1}_Nt}

def mean(xv, tv, x, kernel, beta):

kv = k(xv, x, kernel)
cn = Cn(xv, x, kernel, beta)
return kv.T.dot(numpy.linalg.inv(cn)).dot(tv)

\sigma^2(x_{N+1}) = c - \mathbf{k^TC^{-1}_Nk}

def s2(xv, x, kernel, beta):

c = kernel(x, x) + 1 / beta
cn = Cn(xv, x, kernel, beta)
kv = k(xv, x, kernel)
return c - kv.T.dot(numpy.linalg.inv(cn)).dot(kv)

上記コード中の kernel はカーネル関数です。

ここでの $\mathbf{C_N}$は$N \times N$の共分散行列で、要素として $k(x_n, x_m)$ (カーネル関数に $x_n, x_m$ を与えたもの)を持つ行列に、単位行列を$\beta^{-1}$倍した行列を足しあわせたものになります。$C(x_n, x_m)$は次のようになります(6.62)。

C(x_n, x_m) = k(x_n, x_m) + \beta ^{-1} \delta_{nm}

$\mathbf{C_N}$をコードにすると次のようになります。

def Cn(xv, x, kernel, beta):

N = len(xv)
I = numpy.identity(N)
return numpy.array(
[kernel(xi, xj) for xi in xv for xj in xv]).reshape(N, N) + I / beta

また $\mathbf{k}$ は要素 $k(x_n, x_{N+1}) (n=1, \cdots , N)$ を持つベクトルです。

def k(xv, x, kernel):

N = len(xv)
return numpy.array([kernel(xi, x) for xi in xv]).reshape(N, 1)


描画

グラフは分散の範囲を PRML と同じく 95% 区間とします。また、図 6.8 を目測した結果、 x軸 0.0 〜 1.0、y軸 -1.4 〜 1.4 の区間を表示します。

最初から与えられた点は青いドットで表示し、もとの関数は緑の線で、与えられた点から予想される線は赤で表示しています。

PRML のようなグラフにするためにパラメータを手探りで色々試し、$(\theta_0, \theta_1, \theta_2, \theta_3, \beta) = (0.5, 12, 0, 0, 25)$ とすることで結果のグラフが得られました。


コード

上記をすべて記述したコードは github でご覧いただけます。関数と PRML の該当箇所の対応がわかりやすいように、コメントで対応する式の番号を記述してあります。

github


終わりに

とりあえず動くものは実装できましたが、記事を書く中で自分がいかに理解できていないかがわかりました。勉強のため、また何か実装して記事にしたいと思います。

明日は @ytakashina さんによる「Petuum について」です!