GPyを使ってガウス過程回帰

## ガウス過程回帰とは
すでにサンプリング済みの入力```x```に対応する出力```y```をもとにして、
新しい入力```x'```に対する出力の予測値```y'```の期待値と分散を返す回帰モデルを作成する。
限られたサンプル点から関数の真の形を予測するときに用いられる。
https://jp.mathworks.com/help/stats/gaussian-process-regression-models.html

ガウス過程回帰モデルを扱うには、GPyというpythonのライブラリを使う。
https://gpy.readthedocs.io/en/deploy/#

## 真の関数の描画
入力を２次元とし、真の関数はそれらを余弦関数に通した値の足し合わせであるとする。

```python:temp.py
import numpy as np

# 関数の定義
def func(x):
    fx = np.sum(np.cos(2 * np.pi * x))
    return fx

xa = np.linspace(-1, 1, 101)
ya = np.linspace(-1, 1, 101)
Xa, Ya = np.meshgrid(xa, ya)

Za = np.zeros([101, 101])
for i in range(len(Xa)):
    for j in range(len(Ya)):
        x = np.array([Xa[i,j], Ya[i,j]])
        Za[i,j] = func(x)

# 描画
import matplotlib.pyplot as plt
fig1 = plt.figure(figsize=(8,8))
ax1 = fig1.add_subplot(111)
ax1.contour(Xa, Ya, Za, cmap="jet", levels=10, alpha=1)
plt.xlim(-1,1)
plt.ylim(-1,1)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/0bbef8e7-79b6-24f2-efd3-4b24a0be0d51.png)

## サンプリング
無駄なくサンプリングする方法として、Sobol sequenceやLatin hypercube sampling等がある。
それらは使わず、ここでは単純にランダムにサンプル点を決定する。
サンプル点数は最初は少なめに、20点とする。

```python:temp.py
import random
random.seed(1)

# ランダムにサンプリング
n_sample = 20
Xa_rand = [random.random()* 2 - 1 for i in range(n_sample)]
Ya_rand = [random.random()* 2 - 1 for i in range(n_sample)]

xlist = np.stack([Xa_rand, Ya_rand], axis=1)
Za_rand = []
for x in xlist:
    Za_rand = np.append(Za_rand, func(x))

# 描画
ax1.scatter(Xa_rand, Ya_rand)
```
先程の図にサンプル点をプロットする。
下半分はまだマシだが、上半分はサンプルが少なく、スカスカである。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/7046e7ba-1cf6-f360-960a-4393f78c0b9a.png)

## ガウス過程回帰
ガウス過程回帰モデルを構築する。
```GPy.kern```でカーネル関数を選択する。ここでは2次元のRBFカーネルとした。
```GPy.models.GPRegression```で回帰モデルを構築し、```model.optimize```でモデルのパラメータをチューニングさせる。

```python:temp.py
import GPy

#　学習用データ
Input = np.stack([Xa_rand, Ya_rand], axis=1)
Output = Za_rand[:,None]

# ガウス過程回帰モデルを構築
kernel = GPy.kern.RBF(2)
model = GPy.models.GPRegression(Input, Output, kernel)
model.optimize(messages=True, max_iters=1e5)

# 描画
model.plot(levels=10)
plt.gcf().set_size_inches(8, 8, forward=True)
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.xlabel("x1")
plt.ylabel("x2")
```

応答曲面をプロットする。
20点とかなり点数が少なかったが、意外と大まかな山谷は再現できている。
上半分の誤差は大きい。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/c787e0bc-cee9-61bb-3562-af8f372cc78d.png)


## 信頼区間の描画
応答曲面の断面を見てみる。

x2=0断面
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/723cbe89-5219-46b5-2ad5-1495a2487231.png)

x1=0断面
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/4c3086db-d84f-d840-dc94-35abffac6f4f.png)

水色の帯は2.5~97.5%の信頼区間を示す。
信頼区間の幅が大きいほど、回帰結果にばらつきが大きい。
x2の大きいところはやはり自信がないらしい。

## サンプリング数を増やした場合
```n_sample = 40```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/5ed4f9f2-0e44-26ab-18ed-5b6381d6f69a.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/21ded3a6-750c-f0f8-3104-e51d735033fe.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/295890c5-08a9-3c8a-e0ae-e1753792dcd9.png)

```n_sample = 100```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/38b79ff2-1089-bbb1-30a9-29d1a8c405f0.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/c13858d7-3cf8-aad3-5080-faef75742f0e.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/578196/3672f742-2a38-6b7e-78f3-ca72cfa00ae2.png)

サンプリングの数が多くなるほど、信頼区間の幅は狭くなり、回帰結果のばらつきが小さくなる。

## まとめ
GPyを用いて、ガウス過程回帰モデルを構築した。