

In [ ]:

    
# See https://github.com/yeahrmek/gp/blob/master/GaussianProcesses.ipynb

%matplotlib inline

import numpy as np
from matplotlib import pyplot

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel as C

import GPy



In [ ]:

    
def plot_with_sigma(x_train, y_train, x, y, sigma, label):
    pyplot.plot(x, y, label=label)
    pyplot.fill_between(x.ravel(), y.ravel() - 2 * sigma.ravel(), y.ravel() + 2 * sigma.ravel(), alpha=0.2)
    pyplot.plot(x_train, y_train, 'om', markersize=8, label='Data')
    pyplot.legend()



In [ ]:

    
import GPy
N = 50
noise_var = 0.05

np.random.seed(1)
x_train = np.linspace(0,10,50)[:,None]
k = GPy.kern.RBF(1)
y_train = np.random.multivariate_normal(np.zeros(N), k.K(x_train) + np.eye(N) * np.sqrt(noise_var)).reshape(-1,1)



In [ ]:

    
import GPy

kernel = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=0.1)
model = GPy.models.GPRegression(x_train, y_train, kernel)
model.optimize()
model.plot()

Bayesian Optimization

Suitable for optimization of "heavy" functions (hard to compute functions)
Uses approximation of the objective function
Takes into account the error of approximation



In [ ]:

    
def f(x):
    return (6 * x - 2)**2 * np.sin(12 * x - 4)  


def lower_confidence_bound(x, model, coefficient=2):
    prediction, std = model.predict(x, return_std=True)
    return prediction.ravel() - coefficient * std.ravel()

def get_new_point(model):
    x_grid = np.linspace(0, 1, 200).reshape(-1, 1)

    lcb = lower_confidence_bound(x_grid, model, 2)

    new_point_index = np.argmin(lcb)
    x_new = x_grid[new_point_index]
    lcb = lcb[new_point_index]
    return x_new, lcb


def plot(x_train, y_train, model, x_new, lcb):
    x_grid = np.linspace(0, 1, 100).reshape(-1, 1)
    y_grid = f(x_grid)

    prediction, std = model.predict(x_grid, return_std=True)
    prediction = prediction.ravel()
    std = std.ravel()

    pyplot.figure(figsize=(8, 6))
    pyplot.plot(x_train, y_train, 'or', markersize=8, label='Training set')
    pyplot.plot(x_grid, y_grid, '--b', linewidth=2, label='True function')
    pyplot.plot(x_grid, prediction, '-k', linewidth=2, label='Approximation')
    pyplot.fill_between(x_grid.ravel(), prediction - 2 * std, prediction + 2 * std, alpha=0.3)
    pyplot.plot(x_new, lcb, 'og', markersize=10, label='New point')
    pyplot.ylim([-15, 20])
    pyplot.legend(loc='best')

def optimization_step(x_train, y_train, model):
    model.fit(x_train, y_train)

    x_new, lcb = get_new_point(model)
    plot(x_train, y_train, model, x_new, lcb)
    pyplot.show()
    
    x_train = np.vstack([x_train, x_new.reshape(-1, 1)])
    y_train = np.vstack([y_train, f(x_new).reshape(-1, 1)])
    return x_train, y_train



In [ ]:

    
np.random.seed(239)

x_train = np.array([0.0, 0.58, 0.38, 0.95]).reshape(-1, 1)
y_train = f(x_train)

kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-1, 1e1))
model = GaussianProcessRegressor(kernel=kernel, normalize_y=True)
model.fit(x_train, y_train)



In [ ]:

    
x_grid = np.linspace(0, 1, 100).reshape(-1, 1)
y_grid = f(x_grid)
prediction, std = model.predict(x_grid, return_std=True)
prediction = prediction.ravel()
std = std.ravel()
pyplot.figure(figsize=(8, 6))
pyplot.plot(x_train, y_train, 'or', markersize=8, label='Training set')
pyplot.plot(x_grid, prediction, '-k', linewidth=2, label='Approximation')
pyplot.fill_between(x_grid.ravel(), prediction - 2 * std, prediction + 2 * std, alpha=0.3)
pyplot.ylim([-15, 20])
pyplot.legend(loc='best')



In [ ]:

    
x_train, y_train = optimization_step(x_train, y_train, model)



In [ ]:

    
x_train, y_train = optimization_step(x_train, y_train, model)



In [ ]:

    
x_train, y_train = optimization_step(x_train, y_train, model)



In [ ]:

    
x_train, y_train = optimization_step(x_train, y_train, model)



In [ ]:

    
x_train, y_train = optimization_step(x_train, y_train, model)



In [ ]:

    
x_train, y_train = optimization_step(x_train, y_train, model)



In [ ]:



In [ ]: