# PRML6.4.2 Gaussian processes for regressionのメモ

#### 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.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)

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に近くなります。