ガウス過程回帰とは pythonでの実装例

import numpy as np
from matplotlib import pyplot as plt
from itertools import product

def kernel(x, xd):    # カーネル関数
    K = np.zeros([len(x), len(xd)])
    for i, j in product(range(len(x)), range(len(xd))):
        K[i][j] = theta[0]*np.exp(-1*(x[i] - xd[j])**2 /theta[1])+(theta[2]*int(x[i]==xd[j]))
    return K

N = 20            # サンプル数
theta = np.array([1, 0.4, 0.1])        # 係数
X_smpl = np.random.uniform(-2*3.14, 2*3.14, N)    # サンプルX値
Y_smpl = np.sin(X_smpl)        # サンプルY値
X_est = np.linspace(-2*3.14, 2*3.14, 100)    # 推定したいX値
std = np.zeros([len(X_est)])        # 標準偏差

K = kernel(X_smpl, X_smpl)
W = np.dot(np.linalg.inv(K), Y_smpl)        # パラメータ
k = kernel(X_est, X_smpl)
mu = np.dot(k, W)        # 平均値
KK = np.dot(k, np.linalg.inv(K))
for i in range(len(X_est)):
    std[i] = np.sqrt(theta[0] + theta[2] - np.dot(KK[i], k[i].T))    # 標準偏差

plt.scatter(X_smpl, Y_smpl)
plt.plot(X_est, np.sin(X_est), c='blue')
plt.plot(X_est, mu, color='red')        # 平均値のプロット
plt.fill_between(X_est, mu + std, mu - std, alpha=.2, color='green')    # 回帰の信頼性(標準偏差)を範囲化
plt.show()