Gaussian-Process Regression

A Prior Over Functions

The previous two projects put priors on clusters. This one puts a prior on functions. A Gaussian process says that any finite collection of function values is jointly Gaussian with covariance given by a kernel, and with Gaussian noise everything that follows is closed form — posterior mean, posterior covariance, and a log marginal likelihood that trades fit against complexity automatically. There is no cross-validation anywhere in this project; the marginal likelihood picks the length-scale, the signal amplitude and the noise level by itself.

The motorcycle data, and the spline it turns out to be

On the motorcycle-crash data the GP finds the dip-and-rebound with no functional form imposed. The notebook is candid about what it gets wrong: with a single noise level the band is too wide over the flat opening and arguably too tight through the violent middle, because the data are heteroscedastic and a constant-noise GP cannot represent that. Set against loess and a smoothing spline, all three trace the same curve — unsurprising, since a smoothing spline is a GP posterior mean. What the GP adds is a principled uncertainty band and hyperparameters chosen by marginal likelihood rather than a hand-tuned span. R's mgcv makes the identity concrete: its bs="gp" basis is the same object approached from the frequentist side.

Building structure into the kernel

The Mauna Loa CO₂ series is where kernels earn their keep. Rather than one all-purpose smoother, the model is composed — a long-length-scale RBF for the rising trend, a periodic kernel locked to one year and damped by a slow RBF for the seasonal cycle, a short RBF for medium-term wiggle, and white noise. Sums and products of valid kernels are valid kernels, so structure can be built the way one would describe the physics. The decomposition then separates cleanly: a trend with a 56-year length-scale, and a seasonal component of 6.6 ppm peak-to-trough — matching, to a tenth of a ppm, what R's independent GAM recovers with a cyclic smooth.

The same method elsewhere

This GP is developed here on the motorcycle and CO₂ series, where the point is what the kernel can express. It is put to a very different test in Gaussian Processes & Splines — the Bayesian Kernel View, which runs the same machinery against gradient boosting and random forests on the 20,640-row California housing data. Because an exact GP costs O(n³), it is fitted there to a 1,500-point subsample and still reaches RMSE 0.611 with 90% intervals covering a measured 91.0% — calibrated, competitive with the ensembles, and beaten by them on accuracy. That scaling limit is itself the finding: it is why gradient boosting rather than the GP is the default on tabular data of this size. Its posterior mean is also, exactly, kernel ridge regression — verified to machine precision in Support Vector Machines & Kernel Methods.

The Band Grows, and Still Misses

Then the question that matters most, and the one a widening band invites: does it actually cover anything? Forecasting ten years past the data, it does not. The GP puts 2012 at 385.9 ppm with a 95% band of [383.3, 388.6]. The actual value is 393.8 — outside the interval, and wrong by 5.8 standard deviations. Describing that as honest extrapolation is precisely backwards.

The cause is specific and worth more than a shrug about forecasting being hard. Fit on the same window, a quadratic OLS lands within 1.4 ppm and a linear OLS is 11.5 ppm low. The GP sits with the linear model — because that is what its kernel assumes. An RBF trend term reverts toward the prior mean once it runs out of data, carrying the local slope forward and no curvature, while atmospheric CO₂ is accelerating. The posterior band is computed under the assumed kernel: it prices uncertainty about the function within that class and is silent about the class being wrong. That is why a 5.8σ miss is not a paradox. Marginal likelihood chooses hyperparameters, not assumptions.

The control that makes the diagnosis stick

The R engine supplies the control that makes the diagnosis stick. Its GAM forecasts 391.4 ppm with a band of [387.6, 395.2] — the truth falls inside, at 1.3 standard errors. A thin-plate spline retains some curvature past the boundary where the RBF kernel does not. Neither band, though, prices the risk that the trend model itself is wrong, which over a ten-year extrapolation of an accelerating series is the dominant risk. Both are honest about the function; neither is honest about the functional form, because no posterior band ever is.

Notebooks

Downloads

GP Module — Source Code

"""
gp.py -- GAUSSIAN-PROCESS REGRESSION (from scratch).

Backs the notebooks in  "Gaussian-Process Regression".

A Gaussian process is a prior over FUNCTIONS: any finite set of function values is jointly
Gaussian, with covariance given by a kernel k(x,x'). Regression with a GP prior and Gaussian
noise has everything in closed form. With y = f(X) + eps, eps ~ N(0, sigma_n^2):

    f ~ GP(0, k),   [f(X), f(X*)] jointly Gaussian
    posterior mean at X*:  m* = K(X*,X) [K(X,X)+sigma_n^2 I]^{-1} y
    posterior cov  at X*:  C* = K(X*,X*) - K(X*,X)[K(X,X)+sigma_n^2 I]^{-1} K(X,X*)

and the kernel HYPERPARAMETERS are learned by maximising the log marginal likelihood

    log p(y) = -1/2 y^T [K+sigma_n^2 I]^{-1} y - 1/2 log|K+sigma_n^2 I| - n/2 log 2pi,

which trades data fit against model complexity automatically -- no cross-validation needed.
Everything is computed stably through the Cholesky factor of K + sigma_n^2 I.

The kernel is where prior structure lives: an RBF/Matern kernel encodes smoothness and a
length-scale; a periodic kernel encodes a cycle; SUMS and PRODUCTS of kernels build structured
priors (the Mauna Loa CO2 model = long trend + periodic seasonal + short-term wiggle + noise).
This module is 1-D (a single input), which covers both notebook datasets.
"""

import numpy as np
from scipy.linalg import cholesky, cho_solve
from scipy.optimize import minimize


# --------------------------------------------------------------------------- #
#  kernel primitives (1-D)                                                      #
# --------------------------------------------------------------------------- #

def _D(xa, xb):
    xa = np.asarray(xa, float).ravel(); xb = np.asarray(xb, float).ravel()
    return xa[:, None] - xb[None, :]

def k_rbf(xa, xb, ell, eta):
    """squared-exponential (infinitely smooth), length-scale ell, variance eta."""
    return eta * np.exp(-0.5 * (_D(xa, xb) / ell) ** 2)

def k_matern52(xa, xb, ell, eta):
    """Matern-5/2 (twice differentiable -- rougher than RBF)."""
    r = np.abs(_D(xa, xb)) / ell
    return eta * (1 + np.sqrt(5) * r + 5.0 / 3.0 * r ** 2) * np.exp(-np.sqrt(5) * r)

def k_periodic(xa, xb, ell, period, eta):
    """exactly periodic kernel (Mackay), period `period`, smoothness ell."""
    return eta * np.exp(-2.0 * np.sin(np.pi * np.abs(_D(xa, xb)) / period) ** 2 / ell ** 2)

def k_linear(xa, xb, c, eta):
    xa = np.asarray(xa, float).ravel(); xb = np.asarray(xb, float).ravel()
    return eta * np.outer(xa - c, xb - c)


# --------------------------------------------------------------------------- #
#  fit (maximise log marginal likelihood) and predict                           #
# --------------------------------------------------------------------------- #

def _neg_lml(theta, X, y, build):
    kfun, nv = build(theta)
    n = len(X); K = kfun(X, X) + (nv + 1e-8) * np.eye(n)
    try:
        L = cholesky(K, lower=True)
    except np.linalg.LinAlgError:
        return 1e25
    alpha = cho_solve((L, True), y)
    return 0.5 * y @ alpha + np.log(np.diag(L)).sum() + 0.5 * n * np.log(2 * np.pi)

def fit_gp(X, y, build, theta0, bounds=None, restarts=1, rng=None):
    """maximise the log marginal likelihood over hyperparameters `theta` (usually log-scale).
    `build(theta)` returns (kfun, noise_var), where kfun(Xa,Xb) is the signal covariance.
    Returns (theta_hat, log_marginal_likelihood)."""
    best = None
    for r in range(restarts):
        t0 = np.asarray(theta0, float) if r == 0 else np.asarray(theta0, float) + 0.4 * rng.standard_normal(len(theta0))
        res = minimize(_neg_lml, t0, args=(X, y, build), method="L-BFGS-B", bounds=bounds)
        if best is None or res.fun < best.fun:
            best = res
    return best.x, -best.fun

def predict_gp(theta, X, y, Xs, build):
    """posterior predictive mean and covariance of the latent function at Xs."""
    kfun, nv = build(theta); n = len(X)
    K = kfun(X, X) + (nv + 1e-8) * np.eye(n); L = cholesky(K, lower=True)
    alpha = cho_solve((L, True), y)
    Ks = kfun(Xs, X); mean = Ks @ alpha
    v = cho_solve((L, True), Ks.T); cov = kfun(Xs, Xs) - Ks @ v
    return mean, cov

def sample_gp(mean, cov, n, rng):
    """draw n function samples from a Gaussian with the given mean and covariance."""
    m = np.asarray(mean, float); C = np.asarray(cov, float)
    L = cholesky(C + 1e-9 * np.eye(len(m)), lower=True)
    return m[:, None] + L @ rng.standard_normal((len(m), n))


# --------------------------------------------------------------------------- #
#  simulation for validation                                                    #
# --------------------------------------------------------------------------- #

def simulate_regression(f, x, noise_sd, rng):
    """y = f(x) + Normal(0, noise_sd) -- a known function to recover."""
    x = np.asarray(x, float)
    return f(x) + noise_sd * rng.standard_normal(len(x))

References

• Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. — the standard reference; the Mauna Loa composed kernel is its worked example
• Kimeldorf, G. & Wahba, G. (1970). A correspondence between Bayesian estimation on stochastic processes and smoothing by splines. Annals of Mathematical Statistics 41(2), 495–502. — why the smoothing spline is a GP posterior mean
• MacKay, D. J. C. (1998). Introduction to Gaussian processes. In Neural Networks and Machine Learning. — the periodic kernel used for the seasonal term
• Duvenaud, D. (2014). Automatic Model Construction with Gaussian Processes. PhD thesis, Cambridge. — kernel composition, and the extrapolation behaviour of the RBF that this project's forecast test exposes
• Silverman, B. W. (1985). Some aspects of the spline smoothing approach to non-parametric regression curve fitting. JRSS-B 47(1), 1–52. — the motorcycle data as a smoothing benchmark
• Wood, S. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.). CRC Press. — mgcv, its Gaussian-process basis and the additive CO₂ decomposition