ガウス過程による回帰（３）実験結果 - 機械学習に詳しくなりたいブログ

概要

今回は、以下の記事を経て導出した結果

$\displaystyle m(\mathbf{x}_{N+1}) = \sum_{n=1}^{N} a_{n} k(\mathbf{x}_n,\mathbf{x}_{N+1}) \tag{1}$

を使って、ガウス過程による回帰を実験してみました。

www.iwanttobeacat.com

ガウス過程による回帰の実験結果

早速結果です。カーネル関数はガウスカーネルを使いました。また、予測分布 $p(t_{N+1}| \mathbf{t})$ も求めてみました。カーネルを用いた回帰でも、カーネル回帰分析（２）実験結果とはずいぶん傾向が違いますね。今回の訓練データの場合では、ガウスカーネルでなくても、多項式モデルから $\mathbf{K}$ を計算してもちゃんと予測分布が求められました。事前分布のパラメータ $\alpha$ の定数倍を微調整する必要がありましたが。ちなみに、線形回帰をベイズ推定で解く（２）予測分布をプロットで求めた予測分布は、線形モデルのパラメータ $\mathbf{w}$ の事前分布を設定し、訓練データから事後分布を求め、 $\mathbf{w}$ について積分することで予測分布を求めていました。ガウス過程による回帰も、出発点はパラメータ $\mathbf{w}$ の事前分布を設定することでしたが、考え方がいろいろあって面白いですね。

ガウスカーネルの分散を小さくしてみたら面白い結果になりました。ガウスカーネルの分散が小さいということは、隣の訓練データの影響が小さくなるのでいびつな回帰曲線になってしまうんですね。

サンプリング

ガウス過程で実験してみたようにサンプリングしてみました。

左がガウス過程の記事でサンプリングした結果。 $y$ の期待値は0でしたから、0付近をうろうろした結果となっています。一方で平均値が式(1)で求められた右のサンプリングはどれも訓練データである正弦波を再現するようなサンプリングとなりました。

今回のコードです。

# ガウス過程による回帰
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


def k(x0, x1):
    sigma = 0.3
    return np.exp(-(x0 - x1) ** 2 / (2 * sigma ** 2))
    # ret = 0
    # for i in range(5):
    #  ret += x0**i * x1**i
    # return ret/0.01


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


# p(t_N+1|t)
def pt(x1, y1):
    k_t = np.vectorize(k)(x.T, np.repeat(x1, N))
    mu = k_t @ a
    v = k(x1, x1) + 1.0/BETA - k_t @ Cn_inv @ k_t
    return np.exp(-(y1 - mu) ** 2 / (2 * v)) / np.sqrt(2*np.pi*v)


# ランダムシードを固定
np.random.seed(0)
# 訓練データ数
N = 10

# 訓練データ
x = np.linspace(0, 1, N)
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)

# 係数aを求める
BETA = 5.0
Cn_inv = np.linalg.inv(gram_matrix(x) + 1.0/BETA*np.eye(N))
a = Cn_inv @ t


# 新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 1, 100)
y = np.empty(100)
for i in range(x2.size):
    y[i] = np.vectorize(k)(x.T, np.repeat(x2[i], N)) @ a


# 結果の表示
plt.xlim(0.0, 1.0)
plt.ylim(-1.5, 1.5)
plt.scatter(x, t)
plt.plot(x2, y)
plt.show()

# 予測分布 p(t_N+1|t)
X, Y = np.meshgrid(np.linspace(-0, 1, 100), np.linspace(-1.5, 1.5, 100))
Z = np.vectorize(pt)(X, Y)
plt.contourf(X, Y, Z, cmap=cm.coolwarm)
plt.scatter(x, t)
plt.xlim(0, 1)
plt.ylim(-1.5, 1.5)
plt.show()