Nadaraya-WatsonモデルをPythonコードを交えて紹介

39659 ワード
テキストリンク

はじめに

PRML(Pattern Recognition and Machine Learning)でNadaraya-Watsonモデルについて学んだ内容をまとめて、人工的に生成したデータを使って学習しました。主に6.3節の内容です。pythonコードは必要に応じて読んでください。

表記について

  • 訓練データの説明変数: \mathbf{X} = \left\{\mathbf{x_1}, \mathbf{x_2}, \cdots \mathbf{x_N}\right\}, \mathbf{x_n} = (x_{n, 1}, x_{n, 2}, \cdots, x_{n_D})^T
  • 訓練データの目的変数: \mathbf{t} = \left\{t_1, t_2, \cdots t_N\right\}
  • 新たな入力データ: \mathbf{x_*}
  • \mathbf{x_*}に対する目的変数: t_*
  • \mathbf{x_*}に対する目的変数の予測: y(\mathbf{x}_*)

近いデータで予測する

「近い」データを使って、新たな入力データに対する出力を予測することを考えます。

import文
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm

np.random.seed(28)

データ

今回は以下のt = sin(x) + 0.1x^2となるデータを使います。

作図コード
line_x = np.linspace(0.0, 2*np.pi, num=500)
line_t = [np.sin(x) + 0.1*x**2 for x in line_x]

plt.plot(line_x, line_t)

func_vanilla
次のように45個のデータを生成しました。今回はこれを訓練データとして学習をします。

\begin{aligned} t_n = \mathcal{N}(t_n | sin(x_n) + 0.1x^2, 0.25) \end{aligned}
作図コード
beta = 0.5

data_x = np.random.rand(45) * 2 * np.pi
data_t = [np.sin(x) + 0.1*x**2 + np.random.normal(0, beta) for x in data_x]

plt.scatter(data_x, data_t)
plt.plot(line_x, line_t)

func_data

最近傍法

最近傍法では、訓練データの中でその説明変数と\mathbf{x_*}との距離（ユークリッド距離など）が最も近いものの目的変数でt_*を予測します。このアルゴリズムでは個々のデータに過度に適合してしまうと考えられます。

k-近傍法

k-近傍法は最近傍法と似ていますが、訓練データの中でその説明変数と\mathbf{x_*}との距離が最も近いものから順にk個選んで、その算術平均でt_*を予測します。適切なkの値を選べば、最近傍法より過適合しなくなり、よい予測ができそうです。

下の図で、左上のk=1となっているのが最近傍法です。kが大きくなるにつれて、予測分布の局所的な変化が小さくなっています。

学習と作図コード
def kNN(k, xs, ts, x_star):
    distance = [[abs(x - x_star), t] for (x, t) in zip(xs, ts)]
    distance.sort()
    t_star = 0
    for i in range(k):
        t_star += distance[i][1]
    t_star /= k
    
    return t_star

def RMSE(true_t, pre_t):
    sum = 0
    for (tt, pt) in zip(true_t, pre_t):
        sum += (tt - pt) ** 2
    rmse = np.sqrt(sum / len(true_t))
    
    return rmse

ks = list(range(1, 21))
kNN_RMSEs = []
nearest_neighbour_ts = [[] for _ in range(len(ks))]

for (i, k) in enumerate(ks):
    for x in line_x:
        nearest_neighbour_ts[i].append(kNN(k, data_x, data_t, x))
    kNN_RMSEs.append(RMSE(line_t, nearest_neighbour_ts[i]))

fig = plt.figure(figsize=(12, 9))

ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 2)
ax3 = fig.add_subplot(2, 2, 3)
ax4 = fig.add_subplot(2, 2, 4)

ax1.scatter(data_x, data_t)
ax1.plot(line_x, line_t)
ax1.plot(line_x, nearest_neighbour_ts[0], label="k=1")

ax2.scatter(data_x, data_t)
ax2.plot(line_x, line_t)
ax2.plot(line_x, nearest_neighbour_ts[2], label="k=3")

ax3.scatter(data_x, data_t)
ax3.plot(line_x, line_t)
ax3.plot(line_x, nearest_neighbour_ts[8], label="k=9")

ax4.scatter(data_x, data_t)
ax4.plot(line_x, line_t)
ax4.plot(line_x, nearest_neighbour_ts[19], label="k=20")

ax1.legend(loc = 'upper left')
ax2.legend(loc = 'upper left')
ax3.legend(loc = 'upper left')
ax4.legend(loc = 'upper left')

plt.show()

knn_prediction

下の図はkに対する（真の分布に対する）平均二乗誤差の値を出力していますが、kが3から10くらいのところで誤差が小さくなっていることがわかります。kの値は大きすぎても小さすぎてもうまくいかないです。

作図コード
plt.plot(ks, kNN_RMSEs)
plt.xlabel("k")
plt.ylabel("RMSE")
plt.show()

knn_RMSE

Nadaraya-Watsonモデル

k-近傍法の考え方をもう少し広げてみて、\mathbf{x_*}と訓練データとの距離に応じて、重み付けを変える方法を考えてみます。ここで、距離が近いほど大きい値をとるような、近さを表現する非負の関数g(\cdot)を考えます。k(\cdot, \cdot)を次のように定義します。

\begin{aligned} k(\mathbf{x_*}, \mathbf{x_n}) = \frac{g(\mathbf{x_*} - \mathbf{x_n})}{\sum^N_{m=1}g(\mathbf{x_*}- \mathbf{x_m})} \end{aligned}

そして、この関数を用いて次のようにt_nに重み付けをして、以下のようなモデルでt_*を予測します。このモデルをNadaraya-Watsonモデルといいます。

\begin{aligned} y(\mathbf{x}_*) = \sum^N_{n=1}k(\mathbf{x_*}, \mathbf{x_n})t_n \end{aligned}

ここで定義より、k(\cdot, \cdot)は以下の2つの条件を満たしています。

\begin{aligned} k(\mathbf{x_*}, \mathbf{x_n}) &\geq 0 \\ \sum^N_{n=1}k(\mathbf{x_*}, \mathbf{x_n}) &= 1 \end{aligned}

次に、以下のように具体的なk(\cdot, \cdot)を定めて、先ほどのデータを用いて予測分布と平均二乗誤差を求めます。必ずしもg(\cdot)をガウス分布にする必要はありません。

\begin{aligned} k(\mathbf{x_*}, \mathbf{x_n}) = \frac{\mathcal{N}(\mathbf{x}_* | \mathbf{x}_n, \sigma^2)}{\sum^N_{m=1}\mathcal{N}(\mathbf{x}_* | \mathbf{x}_m, \sigma^2)} \\ \end{aligned}

k(\mathbf{x_*}, \mathbf{x_n})のパラメータ\sigmaの値によって、どれだけ個々のデータに適合するかが決まります。\sigma=0.36となっている右上の図を見てみると、予測分布がなめらかな曲線になっています。

学習と作図コード
def nadaraya_watson(sigma, xs, ts, x_star):
    prob_sum = 0
    tmp_t = 0
    for (x, t) in zip(xs, ts):
        prob = norm.pdf(x_star, x, sigma)
        prob_sum += prob
        tmp_t += prob * t
    pre_t = tmp_t / prob_sum

    return pre_t

sigmas = [(0.1*i) ** 2 for i in range(1, 21)]
nw_RMSEs = []
nadaraya_watson_ts = [[] for _ in range(len(sigmas))]

for (i, sigma) in enumerate(sigmas):
    for x in line_x:
        nadaraya_watson_ts[i].append(nadaraya_watson(sigma, data_x, data_t, x))
    nw_RMSEs.append(RMSE(line_t, nadaraya_watson_ts[i]))

fig = plt.figure(figsize=(12, 9))

ax1 = fig.add_subplot(2, 2, 1)
ax2 = fig.add_subplot(2, 2, 2)
ax3 = fig.add_subplot(2, 2, 3)
ax4 = fig.add_subplot(2, 2, 4)

ax1.scatter(data_x, data_t)
ax1.plot(line_x, line_t)
ax1.plot(line_x, nadaraya_watson_ts[0], label="$\sigma=0.01$")

ax2.scatter(data_x, data_t)
ax2.plot(line_x, line_t)
ax2.plot(line_x, nadaraya_watson_ts[5], label="$\sigma=0.36$")

ax3.scatter(data_x, data_t)
ax3.plot(line_x, line_t)
ax3.plot(line_x, nadaraya_watson_ts[9], label="$\sigma=1.0$")

ax4.scatter(data_x, data_t)
ax4.plot(line_x, line_t)
ax4.plot(line_x, nadaraya_watson_ts[19], label="$\sigma=4.0$")

ax1.legend(loc = 'upper left')
ax2.legend(loc = 'upper left')
ax3.legend(loc = 'upper left')
ax4.legend(loc = 'upper left')

plt.show()

nadaraya-watson_predition

作図コード
plt.plot(sigmas, nw_RMSEs)
plt.xlabel("sigma")
plt.ylabel("RMSE")
plt.show()

nadaraya-watson_RMSE

理論

モデル

目的変数t_nが、y(\mathbf{z}_n)にガウス分布に従う誤差が加えられて生成されていると考えます。また、\mathbf{z}_nは直接観測できず、何らかの分布g(\mathbf{\xi}_n)に従う\mathbf{\xi}_nを用いて、\mathbf{x}_n = \mathbf{z}_n + \mathbf{\xi}_nの形で観測されると考えます。出力の誤差がガウス分布に従うため、観測値\{\mathbf{x}_n, t_n\} \; (n = 1, \cdots N)と対応する2乗和誤差関数を考えます。

\begin{align} E = \frac{1}{2}\sum^N_{n=1}\int\{y(\mathbf{x}_n - \mathbf{\xi}_n) - t_n\}^2g(\mathbf{\xi}_n)d\mathbf{\xi}_n \end{align}

これを変分法を用いてEをy(\mathbf{z})に関して最小化することで、y(\mathbf{x})の最適解がNadaraya-Watsonモデルで与えられることを示します。

解の導出

まず、y(\mathbf{x})の微小変化を考えます。\epsilonを十分小さくとり、\eta(\mathbf{x})を任意の関数として、

\begin{aligned} y(\mathbf{x}) \rightarrow y(\mathbf{x}) + \epsilon\eta(\mathbf{x}) \end{aligned}

とします。これを(1)の右辺のyに代入して、\tilde{E}とおくと、

\begin{aligned} \tilde{E} &= \frac{1}{2}\sum^N_{n=1}\int\{y(\mathbf{x}_n - \mathbf{\xi}_n) + \epsilon\eta(\mathbf{x}_n - \mathbf{\xi}_n) - t_n\}^2g(\mathbf{\xi}_n)d\mathbf{\xi}_n \\ &= E + \epsilon \cdot \frac{1}{2}\sum^N_{n=1}\int\{y(\mathbf{x}_n - \mathbf{\xi}_n) - t_n\}\eta(\mathbf{x}_n - \mathbf{\xi}_n)g(\mathbf{\xi}_n)d\mathbf{\xi}_n + O(\epsilon^2) \end{aligned}

となります。\epsilonの係数が0になるとき、Eは最小となるので、\epsilonの係数を0とおいて、

\begin{align} \frac{1}{2}\sum^N_{n=1}\int\{y(\mathbf{x}_n - \mathbf{\xi}_n) - t_n\}\eta(\mathbf{x}_n - \mathbf{\xi}_n)g(\mathbf{\xi}_n)d\mathbf{\xi}_n = 0 \end{align}

とします。これは全ての\eta(\mathbf{x})に対して成り立つ必要があるので、

\begin{aligned} \eta(\mathbf{x}) = \delta(\mathbf{x} - \mathbf{z}) \end{aligned}

を選べば良いです。また、こうすることで、\mathbf{\xi}_n上の積分を評価できます。

\begin{aligned} (2) &\Leftrightarrow \sum^N_{n=1}\int\{y(\mathbf{x}_n - \mathbf{\xi}_n) - t_n\}\delta(\mathbf{x}_n + \mathbf{\xi}_n - \mathbf{z})g(\mathbf{\xi}_n)d\mathbf{\xi}_n = 0 \\ &\Leftrightarrow \sum^N_{n=1}\{y(\mathbf{z}) - t_n\}g(\mathbf{z} - \mathbf{x}_n) = 0 \\ &\Leftrightarrow y(\mathbf{z})\sum^N_{n=1}g(\mathbf{z} - \mathbf{x}_n) = \sum^N_{n=1}t_ng(\mathbf{z} - \mathbf{x}_n) \\ &\Leftrightarrow y(\mathbf{z}) = \sum^N_{n=1}t_n\frac{g(\mathbf{z} - \mathbf{x}_n)}{\sum^N_{m=1}g(\mathbf{z} - \mathbf{x}_m)} \\ \end{aligned}

今、\mathbf{z}を\mathbf{x}_*とすれば、

\begin{aligned} y(\mathbf{x}_*) &= k(\mathbf{x}_*, \mathbf{x}_n)t_n \\ where \quad k(\mathbf{x}_*, \mathbf{x}_n) &= \frac{g(\mathbf{x}_* - \mathbf{x}_n)}{\sum^N_{m=1}g(\mathbf{x}_* - \mathbf{x}_m)} \end{aligned}

となります。これで、確かにy(\mathbf{x})の最適解がNadaraya-Watsonモデルで与えられることが示せました。

Nuxt.jsによるVuexのストア分割の仕組み
コンポーネントに型安全にReact Hook Formを適用する