ガウシアンプロセスとベイズ最適化

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from kerasy.ML.linear import BayesianLinearRegression
from kerasy.utils.preprocessing import GaussianBaseTransformer
from kerasy.utils.data_generator import generateSin, generateGausian

N = 25; D = 1; M = 25
xmin=0; xmax=1
seed = 10

alpha = 1
beta = 25

x_test = np.linspace(xmin,xmax,1000)
y_test = np.sin(2*np.pi*x_test)
x_train,y_train = generateSin(size=N, xmin=xmin,xmax=xmax,seed=seed)
_,y_noise = generateGausian(size=N, x=x_train, loc=0, scale=1/beta, seed=seed)
y_train += y_noise

# transform
mu = np.linspace(xmin,xmax,M); sigma = 0.1
phi = GaussianBaseTransformer(mu=mu, sigma=sigma)
x_train_features = phi.transform(x_train)
x_test_features = phi.transform(x_test)

random_idx = np.random.RandomState(seed).choice(N, N)

fig = plt.figure(figsize=(16,8))

for i,n in enumerate([1,2,4,25]):
    ax = fig.add_subplot(2,2,i+1)
    model = BayesianLinearRegression(alpha=alpha, beta=beta)
    model.fit(x_train_features[random_idx[:n]],y_train[random_idx[:n]])
    y_pred,y_std = model.predict(x_test_features)
    ax.plot(x_test,y_test,color="black", label="$\sin(2\pi x)$")
    ax.plot(x_test,y_pred,color="red", label="predict($\mu$)")
    ax.fill_between(x_test, y_pred-y_std, y_pred+y_std, color="pink", alpha=0.5)
    ax.scatter(x_train[random_idx[:n]],y_train[random_idx[:n]],s=150,edgecolors='black',facecolor="white",label="train data")
    ax.set_ylim(-1.5,1.5), ax.legend(), ax.set_title(f"Bayesian Linear Regression (N={n})")
plt.tight_layout()
plt.show()

from kerasy.utils.preprocessing import GaussianBaseTransformer, SigmoidBaseTransformer, PolynomialBaseTransformer

X = np.linspace(-1,1,200)
M = 15
mu = np.linspace(-0.9,0.9,M);
sigma = 0.1
gausian_phi = GaussianBaseTransformer(mu=mu, sigma=sigma)

def mkSmootherMatrix(X, phi, alpha=1, beta=25):
    X_feature = phi.transform(X)
    N,M = X_feature.shape
    SN = np.linalg.inv(alpha*np.eye(M) + beta*X_feature.T.dot(X_feature))
    K = beta*np.dot(X_feature, np.dot(SN, X_feature.T))
    return K

K = mkSmootherMatrix(X, gausian_phi)

SP = [30, 120, 150]
colors = ["red","salmon","pink"]

fig = plt.figure(figsize=(12,6))

axKernel = plt.subplot2grid((3, 2), (0, 1), rowspan=3)
axKernel = sns.heatmap(K, vmin=np.min(K), vmax=np.max(K), cmap="jet", cbar=False, square=True, ax=axKernel)
axKernel.set_xticks([]), axKernel.set_yticks([])
axKernel.set_xlabel("$x$", fontsize=18), axKernel.set_ylabel("$x^{\prime}$", fontsize=18)

for i in range(3):
    axL = plt.subplot2grid((3, 2), (i, 0))
    axL.set_facecolor('lightcyan')
    axL.plot(K[SP[i]], color=colors[i], linewidth=4)
    axL.set_xticks([]), axL.set_yticks([])
    axKernel.axhline(SP[i], color=colors[i], linewidth=10, alpha=0.7)

plt.tight_layout()
plt.show()

Transformers = [
    (GaussianBaseTransformer, mu),
    (SigmoidBaseTransformer, mu),
    (PolynomialBaseTransformer, 10),
]

fig = plt.figure(figsize=(12,4))
for i,(phi_func, param) in enumerate(Transformers):
    ax = fig.add_subplot(1,3,i+1)
    phi = phi_func(param)
    K = mkSmootherMatrix(X, phi)
    ax.plot(K[100], color="blue")
    ax.scatter(100,0,color="red", marker="x", s=150, label="$x^{\prime}$")
    ax.set_xticks([]), ax.set_yticks([]), ax.set_title(f"$k(x,x')$: {phi.__name__}", fontsize=16)
    ax.set_xlabel("$x$", fontsize=18), ax.legend()
plt.tight_layout()
plt.show()

xmin = 0; xmax = 1
N=1000; M=25
seed = 1

alpha = 1
beta = 25

X = np.linspace(xmin,xmax,N)
Y = 5e1*X*(X-0.2)*(X-0.3)*(X-0.7)*(X-1)
_,noise = generateGausian(size=1000, x=X, loc=0, scale=1/beta, seed=seed)
Y_ob = Y+noise

mu = np.linspace(xmin,xmax,M)
sigma = 0.1
phi = GaussianBaseTransformer(mu=mu, sigma=sigma)

X_features = phi.transform(X)

fig = plt.figure(figsize=(18,8))

idx = np.random.RandomState(seed).choice(N, 2) # Initialization.
next_idx = idx
for i,_ in enumerate(range(12)):
    ax = fig.add_subplot(3,4,i+1)
    model = BayesianLinearRegression(alpha=alpha, beta=beta)
    model.fit(X_features[idx], Y_ob[idx])
    y_pred,y_std = model.predict(X_features)
    current_max = np.max(Y_ob[idx])
    max_improve_id = np.argmax(y_pred)
    
    ax.plot(X, Y, color="black", label="True distribution")
    ax.plot(X, y_pred,color="red", label="predict($\mu$)")
    ax.fill_between(X, y_pred-y_std, y_pred+y_std, color="pink", alpha=0.5)
    ax.scatter(X[idx],Y_ob[idx],s=50,edgecolors='black',facecolor="white",label="Observation")
    ax.scatter(X[next_idx],Y_ob[next_idx],s=150,edgecolors='black',facecolor="black",label="New Obsevation")
    ax.set_xlabel("Explanatory Variable"), ax.set_ylabel("Measured Value"), ax.set_title(f"N={len(idx)}, Max={current_max:.3f}")
    ax.set_ylim(-0.5,0.3), ax.set_xticks([]), ax.set_yticks([])
    ax.axhline(current_max, color="green", label="Current Maximum")
    ax.axvline(np.max(X[max_improve_id]), color="blue", label="Maximum Probability of Improvement")
    
    next_idx = max_improve_id if np.random.uniform()>0.5 else np.random.choice(N, 1)
    idx = np.r_[idx, next_idx]
ax.legend()
plt.tight_layout()    
plt.show()