More than 3 years have passed since last update.

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

Last updated at 2020-09-02Posted at 2020-08-30

ガウス過程回帰とは

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

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

真の関数の描画

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

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)

サンプリング

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

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)

先程の図にサンプル点をプロットする。
下半分はまだマシだが、上半分はサンプルが少なく、スカスカである。

ガウス過程回帰

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

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")