12 - Gaussian Process Regression Code Demo

1 R implementation

library(plgp)

Loading required package: mvtnorm

Loading required package: tgp

len <- 100
n_path <- 5
X <- seq(0, 5, length = len)
D <- distance(X)
eps <- sqrt(.Machine$double.eps) 
Sigma <- exp(-D/2) + diag(eps, len) 
y <- mvnfast::rmvn(n_path, mu = rep(0, len), sigma = Sigma)
matplot(X, t(y), type = "l", main = "Prior", lwd = 2,
        ylim = c(-3, 3), las = 1, lty = 1, xlab = "x", ylab = "f(x)")

n <- 5
X <- matrix(seq(0, 2*pi, length=n), ncol=1)
y <- sin(X)
D <- distance(X) 
Sigma <- exp(-D/2) + diag(eps, ncol(D))
XX <- matrix(seq(-0.5, 2*pi + 0.5, length=100), ncol=1)
DXX <- distance(XX)
SXX <- exp(-DXX/2) + diag(eps, ncol(DXX))
DX <- distance(XX, X)
SX <- exp(-DX/2) 
Si <- solve(Sigma)
mup <- SX %*% Si %*% y
Sigmap <- SXX - SX %*% Si %*% t(SX)
YY <- rmvnorm(100, mup, Sigmap)
q1 <- mup + qnorm(0.025, 0, sqrt(diag(Sigmap)))
q2 <- mup + qnorm(0.975, 0, sqrt(diag(Sigmap)))
matplot(XX, t(YY)[, 1:5], type="l", col=1:5, lty=1, xlab="x", ylab="f(x)",
        las = 1, main = "Posterior", lwd = 2)
points(X, y, pch=20, cex=2)

matplot(XX, t(YY), type="l", col=c(rep("grey", nrow(XX))), 
        lty=1, xlab="x", ylab="f(x)", las = 1, main = "Prediction with Uncertainty")
points(X, y, pch=20, cex=2)
lines(XX, t(YY)[, 1], col=3, lty = 1)
lines(XX, mup, lwd=2, col = "blue", lty=2)
lines(XX, q1, lwd=2, lty=2, col=2)
lines(XX, q2, lwd=2, lty=2, col=2)

2 Python implementation

import matplotlib.pyplot as plt
import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF

def plot_gpr_samples(gpr_model, n_samples, ax):
    """Plot samples drawn from the Gaussian process model.

    If the Gaussian process model is not trained then the drawn samples are
    drawn from the prior distribution. Otherwise, the samples are drawn from
    the posterior distribution. Be aware that a sample here corresponds to a
    function.

    Parameters
    ----------
    gpr_model : `GaussianProcessRegressor`
        A :class:`~sklearn.gaussian_process.GaussianProcessRegressor` model.
    n_samples : int
        The number of samples to draw from the Gaussian process distribution.
    ax : matplotlib axis
        The matplotlib axis where to plot the samples.
    """
    x = np.linspace(0, 2 * np.pi, 100)
    X = x.reshape(-1, 1)

    y_mean, y_std = gpr_model.predict(X, return_std=True)
    y_samples = gpr_model.sample_y(X, n_samples)

    for idx, single_prior in enumerate(y_samples.T):
        ax.plot(
            x,
            single_prior,
            linestyle="--",
            alpha=0.7,
            label=f"Sampled function #{idx + 1}",
        )
    ax.plot(x, y_mean, color="black", label="Mean")
    ax.fill_between(
        x,
        y_mean - y_std,
        y_mean + y_std,
        alpha=0.1,
        color="black",
        label=r"$\pm$ 1 std. dev.",
    )
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_ylim([-3, 3])

rng = np.random.RandomState(4)

X_train = np.linspace(0, 2 * np.pi, 5).reshape(-1, 1)
y_train = np.sin(X_train)

n_samples = 5

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF

kernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-1, 10.0))
gpr = GaussianProcessRegressor(kernel=kernel, random_state=0)

fig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))

# plot prior
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])
axs[0].set_title("Samples from prior distribution")

# plot posterior
gpr.fit(X_train, y_train)

GaussianProcessRegressor(kernel=1**2 * RBF(length_scale=1), random_state=0)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])
axs[1].scatter(X_train[:, 0], y_train, color="red", zorder=10, label="Observations")
axs[1].legend(bbox_to_anchor=(1.05, 1.5), loc="upper left")
axs[1].set_title("Samples from posterior distribution")

fig.suptitle("Radial Basis Function kernel", fontsize=18)
plt.tight_layout()
plt.show()