GM

import numpy as np
from matplotlib.pyplot import figure, show, plot, scatter, gca


def k(xs, ys, sigma=1, l=1):
    dx = np.expand_dims(xs, 1) - np.expand_dims(ys, 0)
    return (sigma ** 2) * np.exp(-((dx / l) ** 2) / 2)


def m(x):
    """The mean function. As discussed, we can let the mean always be zero."""
    return np.zeros_like(x)


def f(x):
    # coefs[i] is the coefficient of x^i
    coefs = [6, -2.5, -2.4, -0.1, 0.2, 0.03]

    total = 0

    for i, coef in enumerate(coefs):
        total += coef * (x ** i)
    return total


x_obs = np.array([-4, -1.5, 0, 1.5, 2.5, 2.7])
y_obs = f(x_obs)

x_s = np.linspace(-8, 7, 80)

K = k(x_obs, x_obs)
K_s = k(x_obs, x_s)
K_ss = k(x_s, x_s)

K_sTKinv = np.matmul(K_s.T, np.linalg.pinv(K))

mu_s = m(x_s) + np.matmul(K_sTKinv, y_obs - m(x_obs))
Sigma_s = K_ss - np.matmul(K_sTKinv, K_s)

p = figure(figsize=(8, 6))
axes = gca()
axes.set_ylim([-7, 8])

# real line
y_true = f(x_s)
plot(x_s, y_true, linewidth=3, color='black', linestyle="-", alpha=0.4)

scatter(x_obs, y_obs, color='red')

stds = np.sqrt(Sigma_s.diagonal())
err_xs = np.concatenate((x_s, np.flip(x_s, 0)))
err_ys = np.concatenate((mu_s + 2 * stds, np.flip(mu_s - 2 * stds, 0)))
# patch(err_xs, err_ys, alpha=0.2, line_width=0, color='grey')

for color in ['yellow', 'green']:
    y_s = np.random.multivariate_normal(mu_s, Sigma_s)
    plot(x_s, y_s, linewidth=1, color=color)

plot(x_s, mu_s, linewidth=3, color='blue', alpha=0.4)
show(p)