Figure 6.10を再現するコード
わかりやすさ優先で冗長なコードにしてあります。
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from numpy.random import normal np.random.seed(20160317) num_data = 100 iter_num = 50 rate = 0.005 beta = 36.0 t0 = 1.0 eta1, eta2, eta3 = 1.0, 1.0, 1.0 def prep_data(): x1 = normal(0.0, 0.5, num_data) x2 = x1 + normal(0.0, 0.05, num_data) x3 = normal(0.0, 0.5, num_data) t = np.sin(2.0*np.pi*x1) + normal(0, 1.0/np.sqrt(beta), num_data) return t, x1, x2, x3 def calc_mean(t, xs): m = np.zeros(len(xs)) for i, x in enumerate(xs): k = np.zeros(num_data) for n in range(num_data): k[n] = t0 * np.exp(-0.5*eta1*(x1[n]-x)*(x1[n]-x)) c = matrix_c() c_inv = np.linalg.inv(c) m[i] = np.dot(k, np.dot(c_inv, t)) return m def matrix_c(): c = np.zeros(num_data*num_data).reshape(num_data,num_data) for n in range(num_data): for m in range(num_data): c[n][m] = t0 * np.exp(-0.5 * ( eta1*(x1[n]-x1[m])*(x1[n]-x1[m]) + eta2*(x2[n]-x2[m])*(x2[n]-x2[m]) + eta3*(x3[n]-x3[m])*(x3[n]-x3[m]) )) for n in range(num_data): c[n][n] += 1.0/beta return c def deru_c1(): c = matrix_c() for n in range(num_data): c[n][n] -= 1.0/beta for n in range(num_data): for m in range(num_data): c[n][m] *= -0.5 * (x1[n]-x1[m])*(x1[n]-x1[m]) return c def deru_c2(): c = matrix_c() for n in range(num_data): c[n][n] -= 1.0/beta for n in range(num_data): for m in range(num_data): c[n][m] *= -0.5 * (x2[n]-x2[m])*(x2[n]-x2[m]) return c def deru_c3(): c = matrix_c() for n in range(num_data): c[n][n] -= 1.0/beta for n in range(num_data): for m in range(num_data): c[n][m] *= -0.5 * (x3[n]-x3[m])*(x3[n]-x3[m]) return c # main if __name__ == '__main__': t, x1, x2, x3 = prep_data() eta1s, eta2s, eta3s = [eta1], [eta2], [eta3] for i in range(iter_num): c = matrix_c() c_inv = np.linalg.inv(c) d1, d2, d3 = deru_c1(), deru_c2(), deru_c3() tmp0 = np.dot(t, np.dot(c_inv, np.dot(d1, np.dot(c_inv, t)))) delta = -0.5 * np.trace(np.dot(c_inv, d1)) + 0.5 * tmp0 eta1 += rate * eta1 * delta tmp0 = np.dot(t, np.dot(c_inv, np.dot(d2, np.dot(c_inv, t)))) delta = -0.5 * np.trace(np.dot(c_inv, d2)) + 0.5 * tmp0 eta2 += rate * eta2 * delta tmp0 = np.dot(t, np.dot(c_inv, np.dot(d3, np.dot(c_inv, t)))) delta = -0.5 * np.trace(np.dot(c_inv, d3)) + 0.5 * tmp0 eta3 += rate * eta3 * delta eta1s.append(eta1) eta2s.append(eta2) eta3s.append(eta3) print eta1, eta2, eta3 fig = plt.figure() subplot = fig.add_subplot(2,1,1) subplot.set_yscale('log') subplot.plot(np.arange(0,iter_num+1), eta3s) subplot.plot(np.arange(0,iter_num+1), eta2s) subplot.plot(np.arange(0,iter_num+1), eta1s) subplot = fig.add_subplot(2,1,2) xs = np.linspace(-0.5,0.5,50) subplot.set_xlim(-0.5,0.5) subplot.scatter(x1,t) eta2, eta3 = 0, 0 ms = calc_mean(t, xs) subplot.plot(xs, ms) eta1 = 1.0 ms = calc_mean(t, xs) subplot.plot(xs, ms) plt.show()
実行結果はこちら。
詳細は、テキストを読んでいただくことにして・・・、ざっくり説明すると、次の通りです。
3つの特徴変数 を用いて、 の値を推定する問題です。実際には、 だけが と相関を持っており、 は、 に乱数を加えて擬似的な相関をもたせています。 はまったく無関係な乱数です。
このデータに対して、下記のカーネル関数を用いて、ガウス過程による回帰を実施します。
この時、ハイパーパラメーター の最適化を勾配降下法で実施したのが上のグラフで、(赤)と (緑)に対して、(青)は値が0に近くなります。
下のグラフは、 と置いて、 に対する の推定結果を図示しています。緑は、ハイパーパラメーターを最適化する前の結果で、青は、最適化した後の結果です。ハイパーパラメーターの最適化で大きく結果が変わることが分かります。