機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

ガウス過程

引き続きカーネル法を勉強しています。今回の目的は、線形回帰における出力\mathbf{y}の事前分布はカーネル関数を使って表すことができる、ということの確認です。計算自体は難しくないのですが、自分の中ですっきり理解したと思えるまで大変でした。

まずは線形回帰で何度もやってきた通り、回帰式のモデルを


y = \mathbf{w}^{T}\mathbf{x} \tag{1}

と設定します。線形回帰をMAP推定で解くで計算したのと同じく、係数\mathbf{w}の事前分布を


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

とします。平均が0で、各パラメータには相関がなく、全て等しい分散\alpha^{-1}であるということです。線形回帰をMAP推定で解くの式(5)で \mathbf m_0=0 \mathbf S_0 = \alpha^{-1}\mathbf Iとすれば同じことです。

ここで最小二乗法の解の導出の式(4)と同じように計画行列\mathbf{X}を用いて


\mathbf{y} = \mathbf{X}\mathbf{w} \tag{3}

とします。\mathbf{y} は訓練データx_1, \cdots , x_Nに対する値y_1, \cdots , y_Nをまとめたベクトルになっています。冒頭で書いた通り、今回はこの\mathbf{y}の分布を求めることが目標です。

さて、正規分布に従う独立な確率変数の和は、やはり正規分布に従うという性質があります。正規分布に従うと仮定した係数\mathbf{w}と入力\mathbf{x}の線形和がyですから、この性質よりyは正規分布に従うことがわかります。(参考:再生性 - Wikipedia

今回の目的は\mathbf{y}の分布を求めることでしたから、正規分布に従うということは平均と共分散がわかれば良いということになります。まずは平均を求めますと、


\begin{eqnarray*}
E[\mathbf{y}] &=& E[ \mathbf{X}\mathbf{w}] \tag{4} \\
&=& \mathbf{X} E[\mathbf{w}] \tag{5} \\
&=& \mathbf{0} \tag{6}
\end{eqnarray*}

です。\mathbf{w}は式(2)で平均\mathbf{0}としていますから式(5)は\mathbf{0}と計算できます。

次に\mathbf{y}の共分散を求めます。多変量正規分布の式(2)の通り、


\mathrm{Cov}(y_a,y_b) = E[y_a y_b] - E[y_a]E[y_b]  \tag{7}

です。ここで式(6)より E [ \mathbf{y} ] = \mathbf{0}でしたから、


\mathrm{Cov}(y_a,y_b) = E[y_a y_b]   \tag{8}

です。共分散行列は多変量正規分布の式(4)ですので、


\begin{eqnarray*}
\mathrm{Cov} [ \mathbf{y}  ] &=& E[  \mathbf{y}\mathbf{y}^{T} ] \tag{9} \\
&=& E[ \mathbf{X}\mathbf{w} \mathbf{w}^{T}\mathbf{X}^{T} ] \tag{10} \\
&=&\mathbf{X} E[ \mathbf{w} \mathbf{w}^{T} ] \mathbf{X}^{T} \tag{11} \\
\displaystyle &=&  \mathbf{X} \frac{1}{\alpha} \mathbf{I} \mathbf{X}^{T} \tag{12} \\
\displaystyle &=&  \frac{1}{\alpha} \mathbf{X}  \mathbf{X}^{T} \tag{13} \\
\end{eqnarray*}

と計算できます。式(9)から式(10)は転置行列の定理を使っています。ここでカーネル法による正則化最小二乗法(1)の式(10)より、


\mathbf{K}= \mathbf{X}\mathbf{X}^{T}   \tag{14}

でした。以上をまとめると、求めたい分布は


p(\mathbf{y}) = N(\mathbf{y}|\mathbf{0},\alpha^{-1}\mathbf{K}) \tag{15}

となることがわかります。カーネル関数さえ定まればいきなりyの事前分布が求まってしまうっていうことなんですが、なんか驚きです。

さて、ここでようやくタイトルのガウス過程の話です。ガウス過程とは、確率過程において任意の個数の入力を選んだとき、それに対する出力が常に多次元正規分布に従うことをいいます。ちょっと説明が正確でないかもしれません、、、Wikiを参考にしました:ガウス過程 - Wikipedia p(\mathbf{y})は式(15)の通り多次元正規分布ですから、これを満たしています。つまりガウス過程となっています。

実験結果です。分布は式(15)の通り、多次元正規分布になることは明らかですが、どのような回帰曲線が得られるかいくつかサンプリングしてみました。

f:id:opabinia2:20190603234940p:plain f:id:opabinia2:20190603234950p:plain

上のグラフは、5次の多項式モデルを考え、\mathbf{K}を計算しています。下のグラフはガウスカーネル\mathbf{K}を計算しています。情報のない事前分布なので、なんとなく0近辺をうろうろしているだけという感じでしょうか。

続きは次回:ガウス過程による回帰(1)

今回のコードです。

# ガウス過程からのサンプリング
import matplotlib.pyplot as plt
import numpy as np


# カーネル
def kernel(x1, x2):
    # 5次の多項式モデルを考えたときのカーネル関数(内積)
    ret = 0
    for i in range(5):
        ret += x1**i * x2**i
    return ret
    # ガウスカーネル
    # sigma = 0.3
    # return np.exp(-(x1 - x2) ** 2 / (2 * sigma ** 2))



# グラム行列
def gram_matrix(x):
    gram = np.empty([N, N])
    for i in range(N):
        for j in range(N):
            gram[i, j] = kernel(x[i], x[j])
    return gram


# x
N = 100
x = np.linspace(-1, 1, N)

# α
ALPHA = 1.0

# グラム行列K
K = gram_matrix(x) / ALPHA

for i in range(5):
    y = np.random.multivariate_normal(np.zeros(N), K)
    plt.plot(x, y)


# 結果の表示
plt.xlim(-1, 1)
plt.show()