Intro to Gaussian Processes

import numpy as np
import seaborn as sns

# matplotlib animation
import matplotlib.pyplot as plt
from matplotlib import rc
from matplotlib.animation import FuncAnimation
from matplotlib import animation

from scipy.stats import multivariate_normal
from tqdm.notebook import tqdm

# Plotly
import plotly.graph_objects as go

plt.style.use("ggplot")
rc('animation', html='jshtml')

N = 100 # Number of points in mesh = N * N
frn = 30 # Number of frames to process in the animation
fps = 8 # Frames per second
mywriter = animation.PillowWriter(fps=fps) 

sns.set_style('darkgrid') # darkgrid, white grid, dark, white and ticks
plt.rc('axes', titlesize=18)     # fontsize of the axes title
plt.rc('axes', labelsize=14)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=13)    # fontsize of the tick labels
plt.rc('ytick', labelsize=13)    # fontsize of the tick labels
plt.rc('legend', fontsize=13)    # legend fontsize
plt.rc('font', size=13)          # controls default text sizes

Bi-variate Gaussian Distribution

We have all seen the single variable Gaussian Distribution. For Gaussian Processes, the next step is to check multi-variate Gaussian Distribution.

For $n$-dimensional Gaussian Distribution, the mean $\mu$ is described a 1D vector of size $n$ and a covariance matrix of size $(n \times n)$.

mean = [0.5, 0.3]
cov = [[1, 0.2], 
       [0.2, 1]]

start, end = -3, 3

distr = multivariate_normal(mean = mean, cov = cov, seed = 42)


x = np.linspace(start, end, num = N)
y = np.linspace(start, end, num = N)
pdf = np.zeros((len(x), len(y)))
X, Y = np.meshgrid(x, y)

Z = np.zeros((*X.shape, frn))

fig = go.Figure(data = [go.Surface(x = X, y = Y, z = Z)])
fig.update_layout(title='Gaussian Bivariate Distiru')
fig.show()

import plotly.graph_objects as go


import pandas as pd

# Read data from a csv
z_data = pd.read_csv('https://raw.githubusercontent.com/plotly/datasets/master/api_docs/mt_bruno_elevation.csv')

fig = go.Figure(data=[go.Surface(z=z_data.values)])

fig.update_layout(title='Mt Bruno Elevation', autosize=False,
                  width=500, height=500,
                  margin=dict(l=65, r=50, b=65, t=90))

fig.show()
# pyo.

What happens if we change Mean along $x_{0}$ dimension?

Better check with an animation.

mean_x_values = np.linspace(0, 1, num = frn)

for k in range(len(mean_x_values)):
    mean[0] = mean_x_values[k]
    distr = multivariate_normal(mean = mean, cov = cov, seed = 42)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i][j][k] = distr.pdf((X[i][j], Y[i][j]))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plot = [ax.plot_surface(x, y, Z[:, :, 0], color='0.75', rstride=1, cstride=1)]
ax.set_title("Mean along x_0 {:.4f}".format(mean_x_values[0]))
    
def change_plot(frame_number, Z, plot):
    plot[0].remove()
    plot[0] = ax.plot_surface(X, Y, Z[:, :, frame_number], cmap="coolwarm")
    ax.set_xlabel("x_0")
    ax.set_ylabel("x_1")
    ax.set_title("Mean along x_0 {:.3f}".format(mean_x_values[frame_number]))


ani = FuncAnimation(fig, change_plot, frn, fargs=(Z, plot), interval=1000 / fps)
plt.tight_layout()
display(ani)
mywriter = animation.PillowWriter(fps=fps) 
ani.save('./assets/2022-02-24-intro-gaussian-processes/1_mean_animation.gif',writer=mywriter)
plt.clf()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev = 90, azim = 0)
ax.set_title("Mean along x_0 {:.4f}".format(mean_x_values[0]))
plot = [ax.plot_surface(x, y, Z[:, :, 0], color='0.75', rstride=1, cstride=1)]

ani = FuncAnimation(fig, change_plot, frn, fargs=(Z, plot), interval=1000 / fps)
plt.tight_layout()
display(ani)
ani.save('task1_mean_animation_above.gif',writer=mywriter)
plt.clf()

Varying variance of x0 from 0.1 to 1

mean = [0.5, 0.3]
variance_x_values = np.linspace(0.1, 1, num = frn)

for k in range(len(variance_x_values)):
    cov[0][0] = variance_x_values[k]
    distr = multivariate_normal(mean = mean, cov = cov, seed = 42)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i][j][k] = distr.pdf((X[i][j], Y[i][j]))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("x_0")
ax.set_xlabel("x_1")
plot = [ax.plot_surface(x, y, Z[:, :, 0], color='0.75', rstride=1, cstride=1)]
ax.set_title("Variance along x_0 {:.3f}".format(variance_x_values[0]))
    
def change_plot(frame_number, Z, plot):
    plot[0].remove()
    plot[0] = ax.plot_surface(X, Y, Z[:, :, frame_number], cmap="coolwarm")
    ax.set_xlabel("x_0")
    ax.set_ylabel("x_1")
    ax.set_title("Variance along x_0 {:.3f}".format(variance_x_values[frame_number]))

ax.set_zlim(0, 0.5)

ani = FuncAnimation(fig, change_plot, frn, fargs=(Z, plot), interval=1000 / fps)
ani.save('task1_variance_animation.gif',writer=mywriter)
display(ani)
plt.clf()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev = 90, azim = 0)
ax.set_title("Variance along x_0 {:.3f}".format(variance_x_values[0]))
plot = [ax.plot_surface(x, y, Z[:, :, 0], color='0.75', rstride=1, cstride=1)]

ax.set_zlim(0, 0.5)
ani = FuncAnimation(fig, change_plot, frn, fargs=(Z, plot), interval=1000 / fps)
plt.tight_layout()
ani.save('task1_variance_animation_above.gif',writer=mywriter)

display(ani)
plt.clf()

Varying covariance from 0 to 1

mean = [0.5, 0.3]
cov = [[1, 0.2], 
       [0.2, 1]]
covariance_x_values = np.linspace(-0.99, 0.99, num = frn)

for k in range(len(covariance_x_values)):
    mean[0] = 0.5
    mean[1] = 0.3
    cov[0][0] = cov[1][1] = 1
    cov[0][1] = cov[1][0] = covariance_x_values[k]
    distr = multivariate_normal(mean = mean, cov = cov, seed = 42)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i][j][k] = distr.pdf((X[i][j], Y[i][j]))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("x_0")
ax.set_xlabel("x_1")
ax.set_title("Covariance {:.3f}".format(covariance_x_values[0]))
plot = [ax.plot_surface(x, y, Z[:, :, 0], color='0.75', rstride=1, cstride=1)]
    
def change_plot(frame_number, Z, plot):
    plot[0].remove()
    plot[0] = ax.plot_surface(X, Y, Z[:, :, frame_number], cmap="coolwarm")
    ax.set_xlabel("x_0")
    ax.set_ylabel("x_1")
    ax.set_title("Covariance {:.3f}".format(covariance_x_values[frame_number]))


ax.set_zlim(0, 1)

ani = FuncAnimation(fig, change_plot, frn, fargs=(Z, plot), interval=1000 / fps)
ani.save('task1_covariance_animation.gif',writer=mywriter)

display(ani)
plt.clf()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.view_init(elev = 90, azim = 0)
ax.set_title("Covariance {:.3f}".format(covariance_x_values[0]))
plot = [ax.plot_surface(x, y, Z[:, :, 0], color='0.75', rstride=1, cstride=1)]

ax.set_zlim(0, 1)
ani = FuncAnimation(fig, change_plot, frn, fargs=(Z, plot), interval=1000 / fps)
plt.tight_layout()
display(ani)
ani.save('task1_covariance_animation_above.gif',writer=mywriter)

plt.clf()

Task 2 and 3

def get_prob(X, mean, inv_cov, det_cov):
    k = inv_cov.shape[0]
    X = X.reshape((k, -1))
    mean_reshape = mean.reshape((k, -1))
    numerator = np.matmul(np.matmul((X - mean_reshape).T, inv_cov), X - mean_reshape)
    numerator = np.exp(-0.5 * numerator)
    denominator = np.sqrt(np.power(2 * np.pi, k) * det_cov)
    return numerator / denominator

def draw_samples(mean, cov, a, b, n_points):
    points = []
    k = cov.shape[0]
    inv_cov = np.linalg.inv(cov)
    det_cov = np.linalg.det(cov)
    for i in tqdm(range(n_points)):
        while True:
            pt = (b - a) * np.random.random(size = k) + a
            pxyz = get_prob(np.array(pt), mean, inv_cov, det_cov)
            if pxyz >= np.random.random():
                points.append(pt)
                break
    points = np.array(points)
    return points

Task 3 (dim = 3)

mean = np.array([0.5, 0.3, 0.8])
cov = np.array([[0.8, 0.2, 0], 
       [0.2, 0.9, 0], 
       [0, 0, 0.6]])

points = draw_samples(mean, cov, -3, 3, int(1e4))

empirical_mean = np.mean(points, axis = 0)
print("Empirical Mean", empirical_mean)
print("Actual Mean", mean)

empirical_cov = np.zeros_like(cov)
for i in range(cov.shape[0]):
    for j in range(cov.shape[1]):
        empirical_cov[i][j] = np.mean((points[:, i] - empirical_mean[i]) * (points[:, j] - empirical_mean[j]))
display("Empirical Covariance matrix", empirical_cov)
display("Actual Covariance matrix", cov)

Task 3 (dim = 5)

mean = np.array([0.5, 0.3, 0.8, 0.1, 0.2])
cov = np.array(
    [  [0.8, 0.2, 0, 0, 0], 
       [0.2, 0.9, 0, 0, 0], 
       [0, 0, 0.6, 0, 0], 
       [0, 0, 0, 0.9, 0], 
       [0, 0, 0, 0, 0.8]])

points = draw_samples(mean, cov, -2, 2, int(1e4))

empirical_mean = np.mean(points, axis = 0)
print("Empirical Mean", empirical_mean)
print("Actual Mean", mean)

empirical_cov = np.zeros_like(cov)
for i in range(cov.shape[0]):
    for j in range(cov.shape[1]):
        empirical_cov[i][j] = np.mean((points[:, i] - empirical_mean[i]) * (points[:, j] - empirical_mean[j]))
display("Empirical Covariance matrix", empirical_cov)
display("Actual Covariance matrix", cov)

Box-Muller Method

n_points = int(1e4)
points = [(np.random.random(), np.random.random()) for i in range(n_points)]

R_func = lambda x, y : np.sqrt(-2 * np.log(x))
theta_func = lambda x, y : 2 * np.pi * y
gaussian_points = [ [R_func(x, y) * np.cos(theta_func(x, y)), (R_func(x, y) * np.sin(theta_func(x, y))) ] for (x, y) in points ]
gaussian_points = np.array(gaussian_points)

empirical_mean = np.mean(gaussian_points, axis = 0)
print("Empirical Mean", empirical_mean)
empirical_cov = np.zeros((2, 2))
for i in range(empirical_cov.shape[0]):
    for j in range(empirical_cov.shape[1]):
        empirical_cov[i][j] = np.mean((gaussian_points[:, i] - empirical_mean[i]) * (gaussian_points[:, j] - empirical_mean[j]))
display("Empirical Covariance matrix", empirical_cov)

def box_muller_polar_form_standard_normal(size):
    x = np.empty(size)
    for idx, _ in np.ndenumerate(x):
        while True:
            u, v = 2 * np.random.random() - 1, 2 * np.random.random() - 1
            s = u * u + v * v
            if s == 0 or s >= 1:
                continue
            r = np.sqrt((-2 * np.log(s)) / s)
            x[idx] = r * u
            break
    return x

Numpy Method

def numpy_method(mean, cov, size):
    mean = np.array(mean)
    cov = np.array(cov).astype(np.double)

    u, s, v = np.linalg.svd(cov)
    final_shape = [size, mean.shape[0]]
    # x = np.random.standard_normal(size = final_shape).reshape((-1, mean.shape[0]))
    x = box_muller_polar_form_standard_normal(size = final_shape).reshape((-1, mean.shape[0]))
    x = np.dot(x, np.sqrt(s)[:, None] * v)
    x += mean
    return x

samples = numpy_method(mean, cov, int(1e6))

empirical_mean = samples.mean(axis = 0)
print(empirical_mean)

empirical_cov = np.zeros_like(cov)
for i in range(cov.shape[0]):
    for j in range(cov.shape[1]):
        empirical_cov[i][j] = np.mean((samples[:, i] - empirical_mean[i]) * (samples[:, j] - empirical_mean[j]))
display("Empirical Covariance matrix", empirical_cov)
display("Actual Covariance matrix", cov)

Task 4

dim = 100
mean = np.zeros(dim)
points = np.linspace(-5, 5, dim)

Squared Exponential (RBF) Kernel

def square_exponential_kernel(t1, t2, variance, lengthscale):
    assert variance > 0
    assert lengthscale > 0
    return variance * np.exp( -1 * np.power(t1 - t2, 2) / (2 * np.power(lengthscale, 2)))

variance, lengthscale = 0.9, 1.22
cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(dim):
        cov[i][j] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

plt.imshow(cov, cmap = "coolwarm", interpolation='nearest')
plt.colorbar()

Periodic Kernel

def periodic_kernel(t1, t2, variance, lengthscale, periodicity):
    assert variance > 0
    assert lengthscale > 0
    assert periodicity > 0
    sine_value = np.sin(np.pi * np.abs(t1 - t2) / periodicity)
    exp_inner_term = -2 * np.power(sine_value, 2) / np.power(lengthscale, 2)
    return variance * np.exp(exp_inner_term)

variance, lengthscale, periodicity = 0.88, 0.93, 4
cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(dim):
        cov[i][j] = periodic_kernel(points[i], points[j], variance, lengthscale, periodicity)

plt.imshow(cov, cmap = "coolwarm", interpolation='nearest')
plt.colorbar()

Linear Kernel

def linear_kernel(t1, t2, variance, variance_b, offset):
    assert variance > 0
    assert variance_b > 0
    return variance_b + variance * (t1 - offset) * (t2 - offset)

variance, variance_b, offset = 0.3, 0.5, 1.0
cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(dim):
        cov[i][j] = linear_kernel(points[i], points[j], variance, variance_b, offset)

plt.imshow(cov, cmap = "coolwarm", interpolation='nearest')
plt.colorbar()

Task 5

dim = 100
mean = np.zeros(dim)
points = np.linspace(-5, 5, dim)

RBF Kernel

• Change in pathlength from 0 to 5

lengthscale_values = np.linspace(1e-3, 5, frn)
variance = 0.9
cov = np.zeros((dim, dim, frn))

for k in range(len(lengthscale_values)):
    lengthscale = lengthscale_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in lengthscale for RBF Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm")]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm")
    ax.title.set_text("Lengthscale of RBF Kernel {:.3f}".format(lengthscale_values[frame_number]))


plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('rbf_lengthscale_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

variance = 0.9
variance_values = np.linspace(1e-3, 5, frn)
lengthscale = 1.22
cov = np.zeros((dim, dim, frn))

for k in range(len(variance_values)):
    variance = variance_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in variance for RBF Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm", vmax = 5)]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap = "coolwarm", vmax = 5)
    ax.title.set_text("Variance in RBF Kernel {:.3f}".format(variance_values[frame_number]))


ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
plt.colorbar(plot[0])   
ani.save('rbf_variance_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

Periodic Kernel

lengthscale_values = np.linspace(1e-3, 5, frn)
variance, lengthscale, periodicity = 0.88, 0.93, 4
cov = np.zeros((dim, dim, frn))

for k in range(len(lengthscale_values)):
    lengthscale = lengthscale_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = periodic_kernel(points[i], points[j], variance, lengthscale, periodicity)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in lengthscale for Periodic Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm")]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm")
    ax.title.set_text("Lengthscale in Periodic Kernel {:.3f}".format(lengthscale_values[frame_number]))


plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('periodic_lengthscale_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

variance_values = np.linspace(1e-3, 5, frn)
variance, lengthscale, periodicity = 0.88, 0.93, 4
cov = np.zeros((dim, dim, frn))

for k in range(len(variance_values)):
    variance = variance_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = periodic_kernel(points[i], points[j], variance, lengthscale, periodicity)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in variance for Periodic Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm", vmax = 5)]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm", vmax = 5)
    ax.title.set_text("Variance in Periodic Kernel {:.3f}".format(variance_values[frame_number]))


plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('periodic_variance_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

periodicity_values = np.linspace(1e-3, 5, frn)
variance, lengthscale, periodicity = 0.88, 0.93, 4
cov = np.zeros((dim, dim, frn))

for k in range(len(periodicity_values)):
    periodicity = periodicity_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = periodic_kernel(points[i], points[j], variance, lengthscale, periodicity)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in periodicity for Periodic Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm")]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm")
    ax.title.set_text("Periodicity in Periodic Kernel {:.3f}".format(periodicity_values[frame_number]))


plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('periodic_periodicity_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

Linear Kernel

variance, variance_b, offset = 0.3, 0.5, 1.0
variance_values = np.linspace(1e-3, 5, frn)
cov = np.zeros((dim, dim, frn))

for k in range(len(variance_values)):
    variance = variance_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = linear_kernel(points[i], points[j], variance, variance_b, offset)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in variance for Linear Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm", vmax = 20)]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm", vmax = 20)
    ax.title.set_text("Variance in Linear Kernel {:.3f}".format(variance_values[frame_number]))

plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('linear_variance_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

variance, variance_b, offset = 0.3, 0.5, 1.0
variance_b_values = np.linspace(1e-3, 5, frn)
cov = np.zeros((dim, dim, frn))

for k in range(len(variance_b_values)):
    variance_b = variance_b_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = linear_kernel(points[i], points[j], variance, variance_b, offset)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in variance_b for Linear Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm", vmax = 7)]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm", vmax = 7)
    ax.title.set_text("variance_b in Linear Kernel {:.3f}".format(variance_b_values[frame_number]))


plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('linear_variance_b_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

variance, variance_b, offset = 0.3, 0.5, 1.0
offset_values = np.linspace(1e-3, 5, frn)
cov = np.zeros((dim, dim, frn))

for k in range(len(offset_values)):
    offset = offset_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = linear_kernel(points[i], points[j], variance, variance_b, offset)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.title.set_text("Change in offset for Linear Kernel")
plot = [ax.imshow(cov[:, :, 0], cmap = "coolwarm", vmax = 7)]
    
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax.imshow(cov[:, :, frame_number], cmap="coolwarm", vmax = 7)
    ax.title.set_text("offset in Linear Kernel {:.3f}".format(offset_values[frame_number]))


plt.colorbar(plot[0])
ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('linear_offset_parameter.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

Task 6

dim = 50
mean = np.zeros(dim)
points = np.linspace(0, dim, num = dim)

fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (50, 10))
variance, lengthscale = 0.9, 1.6

cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(i, dim):
        cov[i][j] = cov[j][i] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

print(np.linalg.det(cov))
ax1.imshow(cov, cmap = "coolwarm")

mv = np.random.multivariate_normal(np.zeros(dim), cov, size = 5)

xs = np.linspace(1, dim, num = dim)
for point in mv:
    ax2.plot(xs, point, "-")
    ax2.scatter(xs, point)

variance_values = np.linspace(1e-3, 10, num = frn)
n_samples = 5
variance, lengthscale = 0.9, 1.6
fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (20, 5))


cov = np.zeros((dim, dim, frn))
sampled_points = np.zeros((n_samples, dim, frn))

for k in range(frn):
    variance = variance_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = square_exponential_kernel(points[i], points[j], variance, lengthscale)
    sampled_points[:, :, k] = np.random.multivariate_normal(np.zeros(dim), cov[:, :, k], size = n_samples)

xs = np.linspace(1, dim - 1, num = dim)   
plot = [ax1.imshow(cov[:, :, 0], cmap = "coolwarm", vmax = 10)]
for i in range(1, n_samples + 1):
        plot.append(ax2.plot(xs, sampled_points[:, :, 0][i - 1])[0])
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax1.imshow(cov[:, :, frame_number], cmap="coolwarm", vmax = 10)
    for i in range(1, n_samples + 1):
        plot[i].remove()
        plot[i] = ax2.plot(xs, sampled_points[:, :, frame_number][i - 1])[0]
    ax1.title.set_text("Variance in RBF Kernel {:.3f}".format(variance_values[frame_number]))
    ax2.set_ylim(-14, 14)
    ax2.set_xticks(xs)




ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('rbf_pattern_variance.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

lengthscale_values = np.linspace(1e-3, 7, num = frn)
n_samples = 5
variance, lengthscale = 0.9, 1.6
fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (20, 5))

cov = np.zeros((dim, dim, frn))
sampled_points = np.zeros((n_samples, dim, frn))

for k in range(frn):
    lengthscale = lengthscale_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = square_exponential_kernel(points[i], points[j], variance, lengthscale)
    sampled_points[:, :, k] = np.random.multivariate_normal(np.zeros(dim), cov[:, :, k], size = n_samples)

xs = np.linspace(0, dim - 1, num = dim)   
plot = [ax1.imshow(cov[:, :, 0], cmap = "coolwarm")]
for i in range(1, n_samples + 1):
        plot.append(ax2.plot(xs, sampled_points[:, :, 0][i - 1])[0])
def change_plot(frame_number, cov, plot):
    plot[0].remove()
    plot[0] = ax1.imshow(cov[:, :, frame_number], cmap="coolwarm")
    for i in range(1, n_samples + 1):
        plot[i].remove()
        plot[i] = ax2.plot(xs, sampled_points[:, :, frame_number][i - 1])[0]
    ax1.title.set_text("Lengthscale in RBF Kernel {:.3f}".format(lengthscale_values[frame_number]))
    ax2.set_ylim(-4, 4)
    ax2.set_xticks(xs)
    

ani = FuncAnimation(fig, change_plot, frn, fargs=(cov, plot), interval=1000 / fps)
ani.save('rbf_pattern_lengthscale.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

Task 7

def add_train_points(known_points, mean, cov):
    dim = cov.shape[0]
    n_known_points = len(known_points)

    # Known and Unknown indices
    mean_reshape = mean.reshape((-1, 1))
    unknown_indices = [i for i in range(dim)]
    known_indices = np.array(list(known_points.keys()))
    for x in known_points.keys():
        unknown_indices.remove(x)
    unknown_indices = np.array(unknown_indices)

    a = mean_reshape[unknown_indices, :]
    b = mean_reshape[known_indices, :]
    A = cov[np.ix_(unknown_indices, unknown_indices)]
    B = cov[np.ix_(unknown_indices, known_indices)]
    C = cov[np.ix_(known_indices, known_indices)]

    y2 = np.array(list(known_points.values())).reshape(-1, 1)
    inv_C = np.linalg.inv(C)
    new_mean = np.zeros_like(mean_reshape)
    new_cov = np.zeros_like(cov)
    new_mean[known_indices] = y2
    new_mean[unknown_indices] = a + B @ inv_C @ (y2 - b)
    new_cov[unknown_indices[:, None], unknown_indices] = A - B @ inv_C @ B.T

    return new_mean.reshape(-1), new_cov

dim = 20
cov = np.zeros((dim, dim))
mean = np.zeros(dim)
points = np.linspace(0, dim - 1, num = dim)

fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (50, 10))
ax2.set_figwidth = 40
variance, lengthscale = 0.9, 1.6

cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(i, dim):
        cov[i][j] = cov[j][i] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

print(np.linalg.det(cov))
ax1.imshow(cov, cmap = "coolwarm")
mv = np.random.multivariate_normal(np.zeros(dim), cov, size = 5)

xs = np.linspace(0, dim - 1, num = dim)
ax2.set_xticks(xs)
for point in mv:
    ax2.scatter(xs, point)
    ax2.plot(point)

dim = 20
cov = np.zeros((dim, dim))
mean = np.zeros(dim)
points = np.linspace(0, dim - 1, num = dim)

fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (50, 10))
ax2.set_figwidth = 40
variance, lengthscale = 0.9, 1.6

cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(i, dim):
        cov[i][j] = cov[j][i] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

known_points = {5: 2,
                6: 3,
                10: 2,
                16:0.5}
new_mean, new_cov = add_train_points(known_points, mean, cov)
ax1.imshow(new_cov, cmap = "coolwarm")

mv = np.random.multivariate_normal(new_mean, new_cov, size = 5)

xs = np.linspace(0, dim - 1, num = dim)
ax2.set_xticks(xs)
for point in mv:
    ax2.scatter(xs, point)
    ax2.plot(point)

Task 8

dim = 20
mean = np.zeros(dim)
points = np.linspace(0, dim, num = dim)

fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (50, 10))
ax2.set_ylim(-10, 10)
variance, lengthscale = 0.6, 3

cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(i, dim):
        cov[i][j] = cov[j][i] = square_exponential_kernel(points[i], points[j], variance, lengthscale)

print(np.linalg.det(cov))
ax1.imshow(cov, cmap = "coolwarm")

mv = np.random.multivariate_normal(np.zeros(dim), cov, size = 2000)

ax2.plot(mean, color="black", alpha = 0.5)
curve1_ys = [mean[i] + np.sqrt(cov[i][i]) * 1.96 for i in range(dim)]
curve2_ys = [mean[i] -np.sqrt(cov[i][i]) * 1.96 for i in range(dim)]
xs = np.linspace(0, dim - 1, num = dim)
ax2.set_xticks(xs)
ax2.scatter(xs, mean, color = "red")
ax2.fill_between(xs, curve1_ys, curve2_ys, color = "gray", alpha = 0.5)

dim = 20
mean = np.zeros(dim)
points = np.linspace(0, dim - 1, num = dim)
variance, lengthscale = 0.6, 3

fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (50, 10))
cov = np.zeros((dim, dim))
for i in range(dim):
    for j in range(i, dim):
        cov[i][j] = cov[j][i] = square_exponential_kernel(points[i], points[j], variance, lengthscale)


print(np.linalg.det(cov))
taken_points = np.random.choice(np.linspace(0, dim - 1, dim), size = int(0.2 * dim))

known_points = dict()

for pt in taken_points:
    known_points[int(pt)] = np.random.random() * 5

print(known_points)
new_mean, new_cov = add_train_points(known_points, mean, cov)
# new_cov = get_near_psd(new_cov)
ax1.imshow(new_cov, cmap = "coolwarm")
print(np.linalg.det(new_cov))

mv = np.random.multivariate_normal(new_mean, new_cov, size = 20000)

ax2.plot(new_mean, color="black", alpha = 0.5)
curve1_ys = [new_mean[i] + np.sqrt(new_cov[i][i]) * 1.96 for i in range(dim)]
curve2_ys = [new_mean[i] -np.sqrt(new_cov[i][i]) * 1.96 for i in range(dim)]
xs = np.linspace(0, dim - 1, num = dim)
ax2.set_xticks(xs)
scatter_colors = np.array([(255, 0, 0) for i in range(dim)])
scatter_colors[list(map(int, known_points.keys())), :] = [0, 0, 255]
ax2.scatter(xs, new_mean, c = scatter_colors / 255)
ax2.fill_between(xs, curve1_ys, curve2_ys, color = "gray", alpha = 0.5)

Bonus Task

dim = 20
lengthscale_values = np.linspace(1e-3, 5, num = frn)
variance, lengthscale = 0.6, 3
fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (20, 5))

known_points = {0:-2, 1:3, 9:-1, 14:-1}

cov = np.zeros((dim, dim, frn))
new_cov = np.zeros((dim, dim, frn))
new_mean = np.zeros((dim, frn))
mean = np.zeros(dim)
curve1s = np.zeros((dim, frn))
curve2s = np.zeros((dim, frn))

scatter_colors = np.array([(255, 0, 0) for i in range(dim)])
scatter_colors[list(map(int, known_points.keys())), :] = [0, 0, 255]

for k in range(frn):
    lengthscale = lengthscale_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = square_exponential_kernel(points[i], points[j], variance, lengthscale)
    new_mean[:, k], new_cov[:, :, k] = add_train_points(known_points, mean, cov[:, :, k])
    new_cov[:, :, k] = get_near_psd(new_cov[:, :, k])
    curve1s[:, k] = [new_mean[:, k][i] + np.sqrt(new_cov[:, :, k][i][i]) * 1.96 for i in range(dim)]
    curve2s[:, k] = [new_mean[:, k][i] -np.sqrt(new_cov[:, :, k][i][i]) * 1.96 for i in range(dim)]

xs = np.linspace(0, dim - 1, num = dim)
plot = [ax1.imshow(cov[:, :, 0], cmap = "coolwarm"), 
        ax2.plot(new_mean[:, 0], color="black", alpha = 0.5)[0],
        ax2.scatter(xs, new_mean[:, 0], c = scatter_colors / 255), 
        ax2.fill_between(xs, curve1s[:, 0], curve2s[:, 0], color = "gray", alpha = 0.5)]
                         
def change_plot(frame_number, new_mean, new_cov, plot):
    plot[0].remove()
    plot[0] = ax1.imshow(new_cov[:, :, frame_number], cmap="coolwarm")
    plot[1].remove()
    plot[1] = ax2.plot(new_mean[:, frame_number], color="black", alpha = 0.5)[0]
    plot[2].remove()
    plot[2] = ax2.scatter(xs, new_mean[:, frame_number], c = scatter_colors / 255)
    plot[3].remove()
    plot[3] = ax2.fill_between(xs, curve1s[:, frame_number], curve2s[:, frame_number], color = "gray", alpha = 0.5)
    ax1.title.set_text("lengthscale in RBF Kernel {:.3f}".format(lengthscale_values[frame_number]))
    ax2.set_ylim(-4, 12)
    ax2.set_xticks(xs)

ani = FuncAnimation(fig, change_plot, frn, fargs=(new_mean, new_cov, plot), interval=1000 / fps)
ani.save('bonus_pathlength_pattern.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()

dim = 20
variance_values = np.linspace(1e-3, 5, num = frn)
variance, lengthscale = 0.6, 3
fig, (ax1, ax2) = plt.subplots(nrows = 1, ncols = 2, figsize = (20, 5))

known_points = {0:-2, 1:3, 9:-1, 14:-1}

cov = np.zeros((dim, dim, frn))
new_cov = np.zeros((dim, dim, frn))
new_mean = np.zeros((dim, frn))
mean = np.zeros(dim)
curve1s = np.zeros((dim, frn))
curve2s = np.zeros((dim, frn))

scatter_colors = np.array([(255, 0, 0) for i in range(dim)])
scatter_colors[list(map(int, known_points.keys())), :] = [0, 0, 255]

for k in range(frn):
    variance = variance_values[k]
    for i in range(dim):
        for j in range(dim):
            cov[i][j][k] = square_exponential_kernel(points[i], points[j], variance, lengthscale)
    new_mean[:, k], new_cov[:, :, k] = add_train_points(known_points, mean, cov[:, :, k])
    new_cov[:, :, k] = get_near_psd(new_cov[:, :, k])
    curve1s[:, k] = [new_mean[:, k][i] + np.sqrt(new_cov[:, :, k][i][i]) * 1.96 for i in range(dim)]
    curve2s[:, k] = [new_mean[:, k][i] -np.sqrt(new_cov[:, :, k][i][i]) * 1.96 for i in range(dim)]

xs = np.linspace(0, dim - 1, num = dim)
plot = [ax1.imshow(cov[:, :, 0], cmap = "coolwarm"), 
        ax2.plot(new_mean[:, 0], color="black", alpha = 0.5)[0],
        ax2.scatter(xs, new_mean[:, 0], c = scatter_colors / 255), 
        ax2.fill_between(xs, curve1s[:, 0], curve2s[:, 0], color = "gray", alpha = 0.5)]
                         
def change_plot(frame_number, new_mean, new_cov, plot):
    plot[0].remove()
    plot[0] = ax1.imshow(new_cov[:, :, frame_number], cmap="coolwarm", vmax = 4)
    plot[1].remove()
    plot[1] = ax2.plot(new_mean[:, frame_number], color="black", alpha = 0.5)[0]
    plot[2].remove()
    plot[2] = ax2.scatter(xs, new_mean[:, frame_number], c = scatter_colors / 255)
    plot[3].remove()
    plot[3] = ax2.fill_between(xs, curve1s[:, frame_number], curve2s[:, frame_number], color = "gray", alpha = 0.5)
    ax1.title.set_text("variance in RBF Kernel {:.3f}".format(variance_values[frame_number]))
    ax2.set_ylim(-4, 12)
    ax2.set_xticks(xs)

ani = FuncAnimation(fig, change_plot, frn, fargs=(new_mean, new_cov, plot), interval=1000 / fps)
ani.save('bonus_variance_pattern.gif',writer=mywriter)
plt.tight_layout()
display(ani)
plt.clf()