[PYTHON] Gaussian process regression using GPy

What is Gaussian process regression?

Based on the output `y``` corresponding to the input` `x``` that has already been sampled Create a regression model that returns the expected value and variance of the output prediction value y'``` for the new input `` x'```. It is used to predict the true shape of a function from a limited number of sample points. https://jp.mathworks.com/help/stats/gaussian-process-regression-models.html

To work with Gaussian process regression models, we use a python library called GPy. https://gpy.readthedocs.io/en/deploy/#

Drawing a true function

Let the inputs be two-dimensional, and the true function is the sum of the values passed through the cosine function.

temp.py

import numpy as np

#Function definition
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)

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

sampling

Sobol sequence, Latin hypercube sampling, etc. are available as methods for sampling without waste. We do not use them, and here we simply randomly determine the sample points. The number of sample points should be 20 at the beginning.

temp.py

import random
random.seed(1)

#Random sampling
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))

#drawing
ax1.scatter(Xa_rand, Ya_rand)

Plot the sample points in the previous figure. The lower half is still better, but the upper half has few samples and is squishy.

Gaussian process regression

Build a Gaussian process regression model.

GPy.Select a kernel function with kern. Here, it is a two-dimensional RBF kernel.

#### **`GPy.models.Build a regression model with GPRegression and model.Tune the model parameters with optimize.`**

temp.py

import GPy

#Training data
Input = np.stack([Xa_rand, Ya_rand], axis=1)
Output = Za_rand[:,None]

#Build Gaussian process regression model
kernel = GPy.kern.RBF(2)
model = GPy.models.GPRegression(Input, Output, kernel)
model.optimize(messages=True, max_iters=1e5)

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

Plot the response surface. Although the score was quite low at 20 points, the rough mountains and valleys were surprisingly reproduced. The error in the upper half is large.

Confidence interval drawing

Fix one of the two-dimensional inputs to 0 and look at the cross section of the response surface.

temp.py

# x2=0 cross section
model.plot(fixed_inputs=[(1, 0)])
plt.xlim(-1,1)
plt.ylim(-4,4)
plt.xlabel("x1")

# x1=0 cross section
model.plot(fixed_inputs=[(0, 0)])
plt.xlim(-1,1)
plt.ylim(-4,4)
plt.xlabel("x2")

x2 = 0 cross section

x1 = 0 cross section

The light blue band indicates the 2.5-97.5% confidence interval. The wider the confidence interval, the greater the variability in the regression results. It seems that the big part of x2 is still not confident.

When the number of samplings is increased

n_sample = 40

n_sample = 100

As the number of samplings increases, the width of the confidence interval becomes narrower and the variation in the regression results becomes smaller.

Summary

A Gaussian process regression model was constructed using GPy.