カーネル回帰分析（２）実験結果 - 機械学習に詳しくなりたいブログ

先回カーネル回帰分析（１）で導出した結果を使って、回帰問題を解いてみたいと思います。せっかくなので回帰曲線だけではなく、予測分布 $p(t|\mathbf{x})$ も求めてみます。

まずは予測分布 $p(t|\mathbf{x})$ を求めるための準備。 $p(t|\mathbf{x})$ は条件付き確率、同時確率の式(1)より、

$\displaystyle p(t|\mathbf{x}) = \frac{p(\mathbf{x},t)}{p(\mathbf{x})} \tag{1}$

です。 $p(\mathbf{x},t)$ は、カーネル回帰分析（１）の式(1)より

$\displaystyle p(\mathbf{x},t) = \frac{1}{N}\sum_{n=1}^{N} f(\mathbf{x}-\mathbf{x}_n,t-t_{n}) \tag{2}$

です。 $p(\mathbf{x})$ は、条件付き確率、同時確率の式(10)より、

$\displaystyle p(\mathbf{x}) = \int p(\mathbf{x},t) dt \tag{3}$

ですから、

$\displaystyle p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N}\int f(\mathbf{x}-\mathbf{x}_n,t-t_{n})dt \tag{3}$

です。ここでカーネル回帰分析（１）の式(13)より、

$\displaystyle g(\mathbf{x}) = \int f(\mathbf{x}, t) dt \tag{4}$

でしたから、

$\displaystyle g(\mathbf{\mathbf{x}-\mathbf{x}_n}) = \int f(\mathbf{x}-\mathbf{x}_n, t) dt \tag{5}$

です。 $t$ は無限区間で積分しますから、 $-t_n$ の項は影響しませんので、

$\displaystyle p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N} g(\mathbf{\mathbf{x}-\mathbf{x}_n}) \tag{6}$

です。これで予測分布 $p(t|\mathbf{x})$ が求められるようになりました。

さて、肝心のカーネル関数はカーネル密度推定法と同じくガウスカーネルを使い、

$\displaystyle g(\mathbf{x}-\mathbf{x}_n) = \frac{1}{(2 \pi h^{2})^{\frac{D}{2}}} \exp \left( -\frac{ \| \mathbf{x}-\mathbf{x}_n \|^{2}}{2h^{2}} \right) \tag{7}$

とします。

ちょっと前置きが長かったですが、実験結果は以下のようになりました。

f:id:opabinia2:20190526230326p:plain まずは回帰曲線。ちゃんと近似できています。次に予測分布。

f:id:opabinia2:20190526230359p:plain

線形回帰をベイズ推定で解く（２）予測分布をプロットでは回帰曲線のモデルを多項式などと設定していましたが、ノンパラメトリックなカーネル回帰分析ではそういう設定がありませんので、同じ予測分布でもパラメトリック手法とは傾向が違って面白いですね。今回ガウスカーネルを使っていますので、予測分布の結果は、訓練データを中心とした正規分布が配置されていると思えばまあこうなるよねって感じですね。

上記の結果は、平滑化パラメータ $h$ を $0.1$ としていますが、 $0.05$ 、 $0.2$ としてみると、それぞれ次のような結果となりました。

f:id:opabinia2:20190526230843p:plain f:id:opabinia2:20190526230832p:plain

このあたりの結果はカーネル密度推定法の実験と同じですね。平滑化パラメータが小さければ訓練データに強く影響される結果となり、大きければ影響が弱くなだらかな分布となります。

今回のコードです。

# カーネル回帰　Nadaraya-Watsonモデル
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


# g(x-xn)（ガウスカーネル）
def g(x0):
    return np.exp(-(x0-x)**2/(2*(h**2))) / ((2*np.pi*h**2)**(D/2))


# カーネル k(x,x')
def kernel(x0):
    return g(x0)/np.sum(g(x0))


# p(x,t)
def pxt(x0, t0):
    # (x0-x)^2 + (t0-t)^2　は、 ||(x0,t0) - (x,t)||^2 を展開した形
    return np.average(np.exp(-((x0-x)**2 + (t0-t)**2)/(2*(h**2))) / ((2*np.pi*h**2)**(D/2)) )


# p(x)
def px(x0):
    return np.average(g(x0))


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

# ガウスカーネルのパラメータ
h = 0.1
D = 1

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

# 新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 1, 100)
y = np.empty(100)

# 予測値は訓練データとカーネル関数の線形和
for i in range(100):
    y[i] = t @ kernel(x2[i])

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

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