ガウス過程による事前分布からのサンプル

前提条件

まず初めに目的変数 $y$ は、ある重みのベクトル $w$ と、dataから得られた特徴ベクトル $\phi(x)$ の線形結合であらわされるとします。
要するに
$\displaystyle y = w^T \phi(x)$
となっているわけです。このとき、 $w$ の事前分布が
$p(w)= N(w|0,\alpha^{-1}I)$
で与えられているとすると、ガウス分布の線形和もガウス分布なので、yの分布もガウス分布となってくれます。

ガウス分布ということは、期待値と、共分散行列さえ求めてしまえばよいので、yの平均と共分散行列を求めると
$\displaystyle E[y]=\Phi E[w]=0\\ E[y y^T]=E[\Phi w (\Phi w)^T]=\Phi E[w w^T]\Phi = \Phi \alpha^{-1} \Phi^T=\alpha^{-1}\Phi\Phi^T$
となります。

これのイメージは、いくつかのxでのyの値がどういう分布になるのか知りたければ、今から知りたいと思うxの特徴ベクトル $\phi(x)$ をすべて計算して、それを掛け算して計算することのできるグラム行列 $K$ を計算すれば、yの同時分布がどういう多次元のガウス分布に従っているのか計算できますよ、ということです。

例えば $x=1,2,3$ でのyの値 $y(1),y(2),y(3)$ がどういう分布に従うのか計算したければ
1.まず各々の特徴ベクトルを計算して
2.グラム行列を作り
3.同時分布を得る

この同時分布が得られるというのが、ちょっと直感的にわかりづらいところだ個人的には思っていて、すべての点がこの場所にある確率がいくつ、という形で得られるので例えば $x\in[0,1]$ でどういう風なyの値になるのかを知りたいとすると、自分でxをある程度の数に区切って離散化して、そのグラム行列を計算してやる必要があります。なので厳密には連続的なyの分布を得ることはできません。（という理解でいます）

なんか要するに10次元とか、100次元とかの多次元ガウス分布を作ってそこからyを適当によってくると、それがガウス過程の事前分布から得られるあるサンプルのyの値ってことになるわけです。だから個人的にはPRML6章の図6.5のグラフはちょっとなめらかすぎて勘違いしそうだなぁ（というか僕が最初グラフの意味が分からなかったので…）とも思いました。はい。

ということを踏まえて図6.5を作るコードを書いてみました。
numpfspanっていうのが、離散化する数になっています。

# coding: UTF-8
import numpy as np
import matplotlib.pyplot as plt
class kernelfunction():
    def __init__(self,theta1,theta2,theta3,theta4):
        """
        kernel = Θ1*exp(-Θ2/2*||x1-x2||^2) + Θ3 + Θ4*<x1|x2>
        defaultはガウシアンカーネルを使用
        別途追加すれば別のカーネルに変更することもできるので
        """
        self.theta1 = theta1
        self.theta2 = theta2
        self.theta3 = theta3
        self.theta4 = theta4
        self.norm = "gaussian" # default
    def _kernel(self,x1,x2):
        if self.norm == "gaussian":
            return np.exp(-self.theta2/2.*np.inner(x1-x2,x1-x2))
        if self.norm == "ornstein process": #オルンシュタインウーレンベック過程
            return np.exp(-self.theta2*np.power(np.inner(x1-x2,x1-x2),.5))
    def kernel(self,x1,x2):
        """
        calculate kernel
        x1,x2: numpy.array, has same dimention
        """
        val = self.theta1 * self._kernel(x1,x2) + self.theta3 + self.theta4 * (np.inner(x1,x2))
        return val

kernel1 = kernelfunction(1,4.,0.,0.)#thata1=1,theta2=4 よくあるガウスカーネルになってます。
#kernel1.norm = "ornstein process"
numofspan = 100 #離散化する数 増やせばなめらかさは増しますが、計算コストも増えます。
gram_matrix = np.identity(numofspan)
x_range = np.linspace(-1,1,numofspan)
for i in range(numofspan):
    for k in range(numofspan):
        x1 = x_range[i]
        x2 = x_range[k]
        gram_matrix[i][k] = kernel1.kernel(np.array([x1]),np.array([x2]))
#plt.plot(range(len(gram_matrix[0])),np.abs(np.linalg.eig(gram_matrix)[0]))
#plt.show()
"""
#グラム行列の固有値がどうなってるか調べたかったので。
#本来はグラム行列が正定値である条件を満たしていないと、カーネルの条件に反してしまうので
#分布を作れないんですが、np.random.multivariate_normalは勝手に疑似逆行列を生成してくれるみたいで
#（正定値じゃないぞコラ）というエラーを吐きながらちゃんとサンプリングしてくれます。えらい！
"""

color = ["g","r","b","y"]
for i in range(10):
    y = np.random.multivariate_normal(np.zeros(numofspan),gram_matrix,1)
    plt.plot(x_range,y[0],color[i % len(color)])
plt.show()

(1.,4.,0.,0.)
f:id:dette:20150809050548p:plain
(1.,64.,0.,0.)
f:id:dette:20150809050551p:plain
(1.,4.,10.,0)
f:id:dette:20150809051244p:plain
どれもいい感じですね。傾向としては、ガウスカーネルのノルムにかかる係数が大きくなると、二つのxの値が離れた時に急速に0に共分散が近づいてしまうために、近いところのyの値にも相関がどんどんなくなっていることがわかります。逆に、近くないxにも大きな共分散を与えるような条件でやると、近い値どころか遠い値にも影響を及ぼしあって、あんまりyの値が変化しなくなっていることもわかります。