Gaussian Process from Scratch [WIP]
Implement Gaussian Process from Scratch by using numpy and scipy for learning GP purpose
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
import scipy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.gridspec as gridspec
import seaborn as sns
# Set matplotlib and seaborn plotting style
sns.set_style('darkgrid')
np.random.seed(42)
def exp_quadratic(xa, xb):
"""Exponentiated quadratic with σ=1"""
# L2 distance (Squared Euclidian)
sq_norm = -0.5 * scipy.spatial.distance.cdist(xa, xb, 'sqeuclidean')
return np.exp(sq_norm)
X = np.expand_dims(np.linspace(*xlim, 25), 1)
Σ = exp_quadratic(X, X)
plt.imshow(Σ, cmap=cm.YlGnBu);
zero = np.array([[0]])
Σ0 = exp_quadratic(X, zero)
plt.plot(X[:,0], Σ0[:,0]);
n_samples = 100
n_funcs = 8
X = np.expand_dims(np.linspace(-4,4, n_samples), 1)
Σ = exp_quadratic(X, X)
ys = np.random.multivariate_normal(mean=np.zeros(n_samples), cov=Σ, size=n_funcs)
for i in range(n_funcs):
plt.plot(X, ys[i], linestyle='-', marker='o', markersize=3)
plt.xlabel('$x$', fontsize=13)
plt.ylabel('$y = f(x)$', fontsize=13)
plt.title((
f'{n_funcs} different function realizations at {n_samples} points\n'
'sampled from a Gaussian process with exponentiated quadratic kernel'))
plt.xlim([-4, 4])
plt.show()
exponentiated_quadratic = exp_quadratic
A = np.array([[1,-2j],[2j,5]])
A, A.shape
L = np.linalg.cholesky(A)
np.dot(L, L.T.conj())
A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
np.linalg.cholesky(A) # an ndarray object is returned
np.linalg.cholesky(np.matrix(A))
def GP(X1, y1, X2, kernel_func):
cov11 = kernel_func(X1, X1)
cov12 = kernel_func(X1, X2)
solved = scipy.linalg.solve(cov11, cov12, assume_a='pos').T
mu2 = solved @ y1
cov22 = kernel_func(X2, X2)
cov2 = cov22 - (solved @ cov12)
return mu2, cov2
def GP2(X1, y1, X2, kernel_func):
K11 = kernel_func(X1, X1)
K12 = kernel_func(X1, X2)
K22 = kernel_func(X2, X2)
#L = np.linalg.cholesky(K11)
mu2 = K12.T.dot(np.linalg.inv(K11)).dot(y1)
cov2 = K22 - K12.T.dot(np.linalg.inv(K11)).dot(K12)
return mu2, cov2
ny = 10 # Number of functions
domain = (-6, 6)
domain[0]+2, domain[1]-2, (n1, 1)
X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
X1.shape
%%prun
f_sin = lambda x: (np.sin(x)).flatten()
n1 = 40 # Train points
n2 = 75 # Test points
ny = 5 # Number of functions
domain = (-6, 6)
X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
y1 = f_sin(X1)
X2 = np.linspace(domain[0], domain[1], n2).reshape(-1, 1)
# mu2, cov2 = GP(X1, y1, X2, exp_quadratic)
mu2, cov2 = GP2(X1, y1, X2, exp_quadratic)
sigma2 = np.sqrt(np.diag(cov2))
y2 = np.random.multivariate_normal(mean=mu2, cov=cov2, size=ny)
fig, (ax1, ax2) = plt.subplots(
nrows=2, ncols=1, figsize=(6, 6))
# Plot the distribution of the function (mean, covariance)
ax1.plot(X2, f_sin(X2), 'b--', label='$sin(x)$')
ax1.fill_between(X2.flat, mu2-2*sigma2, mu2+2*sigma2, color='red',
alpha=0.15, label='$2 \sigma_{2|1}$')
ax1.plot(X2, mu2, 'r-', lw=2, label='$\mu_{2|1}$')
ax1.plot(X1, y1, 'ko', linewidth=2, label='$(x_1, y_1)$')
# Plot some samples from this function
ax2.plot(X2, y2.T, '-')
ax2.set_xlabel('$x$', fontsize=13)
ax2.set_ylabel('$y$', fontsize=13)
ax2.set_title('5 different function realizations from posterior')
ax1.axis([domain[0], domain[1], -3, 3])
ax2.set_xlim([-6, 6])
plt.tight_layout()
plt.show()