ガウス過程(2) ガウス過程回帰

## 目標変数の事前分布

[関数y(x)の事前分布を定義した](https://qiita.com/_R2_/items/bdbe25fd81f9355515bb)ことで、回帰の目標変数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は正規分布し、

```math
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の周辺分布は

```math
\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で平方完成する。

```math
\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で平方完成する。

```math
\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は

```math
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)を計算するカーネル関数として、

```math
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と対応する

```math
{\bf t}=\begin{pmatrix}t_1\\ \vdots\\ t_N\end{pmatrix}
```

が与えられたとする。新たな入力x_{N+1}に対しt_{N+1}を予測することを考える。つまり、

```math
p(t_{N+1}|{\bf x}_1,...,{\bf x}_N,{\bf t},{\bf x}_{N+1})
```

を求めたい。以後x_1,...,x_{N+1}は省略する。

```math
{\bf t}_{N+1}=\begin{pmatrix}t_1\\ \vdots\\ t_N\\ t_{N+1}\end{pmatrix}
```

と置き、同時分布p(t_{N+1})を考える。今、

```math
\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}
```

より、

```math
p({\bf t}_{N+1})=\mathcal{N}({\bf 0},C_{N+1})
```
ただし、

```math
\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})$の平均と共分散は

```math
\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の同時分布が

```math
\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)の平均と共分散は

```math
\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}$について最大化すれば良い。

## 実装

```python
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()
```

![スクリーンショット 2021-02-16 20.36.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/686685/7222c11a-86a1-268b-9035-ce4037700271.png)

## 参考文献

パターン認識と機械学習　(C.M.ビショップ)