目標変数の事前分布
関数y(x)の事前分布を定義したことで、回帰の目標変数tに対しても事前分布を与えることができる。入力変数{x_n}が与えられたとき、出力{t_n}が次のような回帰モデルに従うことを考える。
$$
t_n=y(x_n)+\epsilon_n,\ \ \ \ \ \epsilon_n\sim\mathcal{N}(0,\beta^{-1})
$$
つまり、y_n=y(x_n)が与えられれば、t_nは正規分布し、
p(t_n|y_n)=\mathcal{N}(y_n,\beta^{-1})\\
\Rightarrow p({\bf t}|{\bf y})=\mathcal{N}({\bf y},\beta^{-1}I_N)
である。今、ガウス過程の定義により、関数yの事前分布は
$$
p({\bf y})=\mathcal{N}({\bf 0},K)
$$
なのだから、出力tの周辺分布は
\begin{align}
p({\bf t})&=\int p({\bf y})p({\bf t}|{\bf y})d{\bf y}\\
&\propto \int\exp\{-\frac{\beta}{2}\ ^T({\bf t}-{\bf y})({\bf t}-{\bf y})-\frac{1}{2}\ ^T{\bf y}K^{-1}{\bf y}\}d{\bf y}
\end{align}
となる。指数関数の中身をyで平方完成する。
\begin{align}
& -\frac{\beta}{2}\ ^T({\bf t}-{\bf y})({\bf t}-{\bf y})-\frac{1}{2}\ ^T{\bf y}K^{-1}{\bf y}\\
&=–\frac{\beta}{2}[\ ^T{\bf y}\{I_N+(\beta K)^{-1}\}{\bf y}-\ ^T{\bf t}{\bf y}-\ ^T{\bf y}{\bf t}+\ ^T{\bf t}{\bf t}]\\
&=-\frac{\beta}{2}\ ^T[{\bf y}-\{I_N+(\beta K)^{-1}\}^{-1}{\bf t}]\{I_N+(\beta K)^{-1}\}[{\bf y}-I_N+(\beta K)^{-1}\}^{-1}{\bf t}]\\
&\ \ \ +\frac{\beta}{2}\ ^T{\bf t}\{I_N+(\beta K)^{-1}\}^{-1}{\bf t}-\frac{\beta}{2}\ ^T{\bf t}{\bf t}
\end{align}
右辺第一項は積分消去できるので、第二項をtで平方完成する。
\begin{align}
\frac{\beta}{2}\ ^T{\bf t}\{I_N+(\beta K)^{-1}\}^{-1}{\bf t}-\frac{\beta}{2}\ ^T{\bf t}{\bf t}&=\frac{\beta}{2}\ ^T{\bf t}\{I_N-(\beta K+I_N)^{-1}\}{\bf t}-\frac{\beta}{2}\ ^T{\bf t}{\bf t}\\
&= -\frac{1}{2}\ ^T{\bf t}(K+\beta^{-1}I_N)^{-1}{\bf t}
\end{align}
したがって、
$$
p({\bf t})=\mathcal{N}({\bf 0},C)
$$
ただし、Cは
C=\begin{pmatrix}
k({\bf x}_1,{\bf x}_1)+\beta^{-1} & k({\bf x}_1,{\bf x}_2) &\ldots&k({\bf x}_1,{\bf x}_N)\\
k({\bf x}_2,{\bf x}_1)&k({\bf x}_2,{\bf x}_2)+\beta^{-1}&&\vdots\\
\vdots&&&&\\
k({\bf x}_N,{\bf x}_1)&\ldots&&k({\bf x}_N,{\bf x}_N)+\beta^{-1}
\end{pmatrix}
である。k(x_n,x_m)を計算するカーネル関数として、
k(x_n,x_m)=\theta_0\exp\{-\frac{\theta_1}{2}\|{\bf x}_n-{\bf x}_m\|^2\}+\theta_2+\theta_3\ ^T{\bf x}_n{\bf x}_m
が広く用いられている。
以上より、ガウス過程の視点から目標データ点の集合上での事前分布をモデル化することができた。
推定
訓練集合として入力x_1,...,x_Nと対応する
{\bf t}=\begin{pmatrix}t_1\\ \vdots\\ t_N\end{pmatrix}
が与えられたとする。新たな入力x_{N+1}に対しt_{N+1}を予測することを考える。つまり、
p(t_{N+1}|{\bf x}_1,...,{\bf x}_N,{\bf t},{\bf x}_{N+1})
を求めたい。以後x_1,...,x_{N+1}は省略する。
{\bf t}_{N+1}=\begin{pmatrix}t_1\\ \vdots\\ t_N\\ t_{N+1}\end{pmatrix}
と置き、同時分布p(t_{N+1})を考える。今、
\begin{align}
p({\bf t})&=\mathcal{N}({\bf 0},C)\\
p({\bf t}_{N+1})&=\mathcal{N}({\bf 0},k({\bf x}_{N+1},{\bf x}_{N+1})+\beta^{-1})
\end{align}
より、
p({\bf t}_{N+1})=\mathcal{N}({\bf 0},C_{N+1})
ただし、
\begin{align}
C_{N+1}&=\begin{pmatrix}
k({\bf x}_1,{\bf x}_1)+\beta^{-1} & k({\bf x}_1,{\bf x}_2) &\ldots&k({\bf x}_1,{\bf x}_N) & k({\bf x}_1,{\bf x}_{N+1})\\
k({\bf x}_2,{\bf x}_1)&k({\bf x}_2,{\bf x}_2)+\beta^{-1}&&\vdots&\vdots\\
\vdots&&&&\\
k({\bf x}_N,{\bf x}_1)&\ldots&&k({\bf x}_N,{\bf x}_N)+\beta^{-1}&k({\bf x}_N,{\bf x}_{N+1})\\
k({\bf x}_{N+1},{\bf x}_1)&\ldots&&k({\bf x}_{N+1},{\bf x}_N&k({\bf x}_{N+1}),{\bf x}_{N+1})+\beta^{-1}
\end{pmatrix}\\
&\\
&=\begin{pmatrix} C_N &{\bf k}\\ ^T{\bf k}&k({\bf x}_{N+1},{\bf x}_{N+1})+\beta^{-1}\end{pmatrix}
\end{align}
(補題)を${\bf t}$とt_{N+1}の同時分布に対して適用することで、$p(t_{N+1}|{\bf t})$の平均と共分散は
\begin{align}
m({\bf x}_{N+1})&=0+\ ^T{\bf k}C_N^{-1}({\bf t}-{\bf 0})=\ ^T{\bf k}C_N^{-1}{\bf t}\\
\sigma^2({\bf x}_{N+1})&=k({\bf x}_{N+1},{\bf x}_{N+1})+\beta^{-1}-\ ^T{\bf k}C_N^{-1}{\bf k}
\end{align}
以上より、学習データが与えられたときの新しい入力x_{N+1}に対する出力t_{N+1}の条件付き分布が得られる。
補題
x_aとx_bの同時分布が
\begin{align}
p(\begin{pmatrix}{\bf x}_a\\ {\bf x}_b\end{pmatrix})&=\mathcal{N}({\bf \mu},\Sigma)\\
{\bf \mu}&=\begin{pmatrix}{\bf \mu}_a\\ {\bf \mu}_b\end{pmatrix}\\
\Sigma&=\begin{pmatrix}\Sigma_{aa} &\Sigma_{ab}\\ \Sigma_{ba}&\Sigma_{bb}\end{pmatrix}
\end{align}
で与えられるとき、条件付き分布p(x_a|x_b)の平均と共分散は
\begin{align}
{\bf \mu}_{a|b}&={\bf \mu}_a+\Sigma_{ab}\Sigma_{bb}^{-1}({\bf x}_b-{\bf \mu}_b)\\
\Sigma_{a|b}&= \Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba}
\end{align}
で与えられる。
超パラメータの学習
超パラメータは最尤推定する。
$$
p({\bf t}|{\bf \theta})=\mathcal{N}({\bf 0},C_N)
$$
より、対数尤度は
$$
\log p({\bf t}|{\bf \theta})=-\frac{1}{2}\log|C_N|-\frac{1}{2}\ ^T{\bf t}C_N^{-1}{\bf t}-\frac{N}{2}\log(2\pi)
$$
これを${\bf \theta}$について最大化すれば良い。
実装
import numpy as np
import random
import torch
import matplotlib.pyplot as plt
#---
# カーネル関数
#---
def Kernel(X, theta):
K_1 = torch.zeros((len(X), len(X)))
for i, x in enumerate(X):
K_1[i] = torch.sum((x - X)**2)
K_1 = theta[0] * torch.exp(-theta[1] * K_1 / 2)
K_2 = torch.ones((len(X), len(X))) * theta[2]
K_3 = theta[3] * (X @ X.T)
return K_1 + K_2 + K_3
#---
# 対数尤度関数
#---
def log_likelihood(C, t):
N = len(C)
det_C = torch.det(C)
inv_C = torch.inverse(C)
ln_p = -torch.log(det_C) / 2 - t.T @ inv_C @ t / 2 - N * torch.log(2*torch.tensor(np.pi)) / 2
return - ln_p
#---
# 勾配降下
#---
def grad_backward(theta, ln_p, lr = 0.001, EPOCH = 1000):
for epoch in range(EPOCH):
ln_p.backward(retain_graph=True)
with torch.no_grad():
theta - lr * theta.grad
theta.requires_grad = True
return theta
#---
# 回帰関数
#---
def gaussian_regression(train_X, train_t, theta, test_X, beta):
N = len(train_X)
N_ = len(test_X)
C = Kernel(train_X, theta) + torch.diag(torch.ones(N)) / beta
ln_p = log_likelihood(C, train_t)
theta = grad_backward(theta, ln_p)
m = np.zeros((N_))
sigma = np.zeros((N_))
for i, x in enumerate(test_X):
k_1 = torch.sum((x - train_X)**2, dim = 1)
k_1 = theta[0] * torch.exp(-theta[1] * k_1 / 2)
k_2 = torch.ones((len(train_X))) * theta[2]
k_3 = theta[3] * (x @ train_X.T)
k = k_1 + k_2 + k_3
m[i] = k.T @ torch.inverse(C) @ train_t
k_1 = theta[0] * torch.exp(-theta[1] * 0 / 2)
k_2 = theta[2]
k_3 = theta[3] * (x @ x.T)
k_ = k_1 + k_2 + k_3
sigma[i] = k_ + 1/beta - k.T @ torch.inverse(C) @ k
return m, sigma
#---
# データの準備
#---
x = np.linspace(0, 2*np.pi, 500)
t = np.sin(x)
train_idx = random.sample(range(500), k=20)
train_X = torch.from_numpy(x[train_idx].reshape(len(train_idx), -1))
train_t = torch.from_numpy(t[train_idx] + np.random.normal(0, 0.1, 20))
test_X = torch.from_numpy(x.reshape(len(x), -1))
#---
# 回帰
#---
theta = torch.tensor([0.7, 0.7, 0.7, 0.7], requires_grad = True)
beta = 0.5
m, sigma = gaussian_regression(train_X, train_t, theta, test_X, beta)
#---
# プロット
#---
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x, t,
linestyle = '-', color = 'orange', label = 'ground truth')
ax.plot(x[train_idx], train_t.numpy().ravel(),
linestyle = 'None', marker = '.',
color = 'blue', label = 'data')
ax.plot(x, m,
linestyle = '-', color = 'red', label = 'pred')
ax.fill_between(x, m - np.sqrt(sigma), m + np.sqrt(sigma),
color = 'pink', alpha = 0.5, label = 'pred std')
ax.legend(loc = 'upper right')
plt.show()
参考文献
パターン認識と機械学習 (C.M.ビショップ)