Gaussian-Process Regression

A prior over functions — learning smooth curves without choosing a form

Everything so far has clustered data or estimated densities. A Gaussian process does something different: it is a prior over functions. Instead of assuming $f$ is linear, quadratic, or has $k$ knots, we say any finite set of function values is jointly Gaussian, with covariance set by a kernel $k(x,x')$ that encodes smoothness and structure. Regression is then exact:

$$f\sim\mathcal{GP}(0,k),\quad y=f(X)+\varepsilon,\ \varepsilon\sim N(0,\sigma_n^2)\ \Rightarrow\ \begin{cases}\ \mathbb E[f(X_*)]=K_{*}[K+\sigma_n^2 I]^{-1}y\\[2pt] \ \operatorname{Cov}[f(X_*)]=K_{**}-K_{*}[K+\sigma_n^2 I]^{-1}K_{*}^\top.\end{cases}$$

The kernel hyperparameters (length-scale, variance, noise) are learned by maximising the log marginal likelihood, which balances fit against complexity with no cross-validation. Everything runs through the Cholesky factor of $K+\sigma_n^2 I$.

We build a GP from scratch, draw functions from the prior, fit the classic motorcycle-crash data (the heteroscedastic nonparametric-regression benchmark), place it beside its frequentist cousins — loess and the smoothing spline (which is a GP in disguise) — then compose a structured kernel for the Mauna Loa CO₂ series (trend + seasonal + noise) and forecast it with honest uncertainty. PyMC's gp.Marginal reproduces the fit.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import gp as G
rng = np.random.default_rng(2)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("A Gaussian process is a prior over functions; the kernel says what kind.")

A Gaussian process is a prior over functions; the kernel says what kind.

1. What is a GP? — functions drawn from a kernel

Before seeing data, a GP is fully described by its kernel. Sampling from the prior shows what each kernel believes: the RBF kernel gives very smooth wiggles whose typical width is the length-scale; Matérn-5/2 is rougher; the periodic kernel gives repeating functions. Shrinking the length-scale makes functions wigglier.

xs=np.linspace(0,10,200)
fig,ax=plt.subplots(1,3,figsize=(13,3.9))
for a,(name,Kf) in zip(ax,[("RBF, ell=1.5",lambda:G.k_rbf(xs,xs,1.5,1.0)),
                            ("Matern-5/2, ell=1.5",lambda:G.k_matern52(xs,xs,1.5,1.0)),
                            ("periodic, period=2",lambda:G.k_periodic(xs,xs,1.0,2.0,1.0))]):
    S=G.sample_gp(np.zeros(len(xs)), Kf(), 5, rng)
    a.plot(xs,S,lw=1.3); a.set_title(name); a.set_xlabel("x"); a.set_ylabel("f(x)"); a.set_ylim(-3,3)
plt.tight_layout(); plt.show()
# length-scale effect
fig,ax=plt.subplots(figsize=(9,3.4))
for ell,c in [(0.5,RED),(1.5,GREEN),(4.0,BLUE)]:
    ax.plot(xs, G.sample_gp(np.zeros(len(xs)),G.k_rbf(xs,xs,ell,1.0),1,rng)[:,0], color=c, lw=1.6, label=f"ell={ell}")
ax.set_title("RBF prior draws: shorter length-scale = wigglier functions"); ax.set_xlabel("x"); ax.legend(frameon=False)
plt.tight_layout(); plt.show()
print("Each kernel is a different belief about the unknown function. Learning the kernel's hyperparameters from")
print("data is how the GP decides how smooth, how wiggly, and how noisy the truth is.")

Each kernel is a different belief about the unknown function. Learning the kernel's hyperparameters from
data is how the GP decides how smooth, how wiggly, and how noisy the truth is.

2. The motorcycle data — GP regression with learned hyperparameters

133 accelerometer readings from a simulated motorcycle crash (time in ms vs head acceleration in g). It is the standard nonparametric-regression benchmark: flat at first, a sharp dip, a rebound, then a noisy tail. We fit a Matérn-5/2 GP, learning the length-scale, signal variance and noise by maximising the marginal likelihood, and read off the posterior mean and 95% band.

d=pd.read_csv("mcycle.csv"); xt=d.times.to_numpy(); yt=d.accel.to_numpy(); ym=yt.mean()
build=lambda t:((lambda A,B:G.k_matern52(A,B,np.exp(t[0]),np.exp(t[1]))), np.exp(t[2]))
th,lml=G.fit_gp(xt, yt-ym, build, [2.,7.,5.], restarts=4, rng=rng)
xg=np.linspace(xt.min(),xt.max(),250); m,C=G.predict_gp(th,xt,yt-ym,xg,build); m=m+ym; sd=np.sqrt(np.diag(C))
print(f"fitted: length-scale {np.exp(th[0]):.2f} ms, signal sd {np.sqrt(np.exp(th[1])):.1f} g, noise sd {np.sqrt(np.exp(th[2])):.1f} g   (log ML {lml:.1f})")
fig,ax=plt.subplots(figsize=(9,4.4))
ax.fill_between(xg, m-1.96*sd, m+1.96*sd, color=BLUE, alpha=.2, label="95% credible band")
ax.plot(xg, m, color=BLUE, lw=2, label="GP posterior mean")
ax.scatter(xt, yt, s=16, color=GREY, alpha=.7, label="data")
ax.set_xlabel("time (ms)"); ax.set_ylabel("acceleration (g)"); ax.set_title("Gaussian-process regression on the motorcycle data")
ax.legend(frameon=False); plt.tight_layout(); plt.show()
print("The GP finds the dip-and-rebound with no functional form imposed -- the marginal likelihood chose the")
print("length-scale. One caveat is visible: with a single NOISE level the band is too wide over the flat start")
print("and arguably too tight in the wild middle -- the data are heteroscedastic, which a constant-noise GP can't")
print("capture (an input-dependent noise GP would). The mean curve is nonetheless the benchmark fit.")

fitted: length-scale 6.55 ms, signal sd 45.7 g, noise sd 22.6 g   (log ML -622.7)

The GP finds the dip-and-rebound with no functional form imposed -- the marginal likelihood chose the
length-scale. One caveat is visible: with a single NOISE level the band is too wide over the flat start
and arguably too tight in the wild middle -- the data are heteroscedastic, which a constant-noise GP can't
capture (an input-dependent noise GP would). The mean curve is nonetheless the benchmark fit.

3. The frequentist counterparts — loess and the smoothing spline

The classical nonparametric-regression tools give a curve but not a probability model. Loess fits a local polynomial in a sliding window; the cubic smoothing spline minimises fit plus a roughness penalty $\lambda\int f''^2$. A beautiful fact ties them to today's model: the cubic smoothing spline is exactly the posterior mean of a GP with a particular (integrated-Wiener) kernel, and the penalty $\lambda$ plays the role of the noise-to-signal ratio. So these are the same curve seen without the uncertainty.

from statsmodels.nonparametric.smoothers_lowess import lowess
from scipy.interpolate import UnivariateSpline
lo = lowess(yt, xt, frac=0.2, return_sorted=True)
spl = UnivariateSpline(xt, yt, s=len(xt)*np.var(yt)*0.35)   # smoothing spline
fig,ax=plt.subplots(figsize=(9,4.4))
ax.scatter(xt, yt, s=15, color=GREY, alpha=.6, label="data")
ax.plot(xg, m, color=BLUE, lw=2.4, label="GP posterior mean")
ax.plot(lo[:,0], lo[:,1], color=RED, lw=1.8, ls="--", label="loess (frequentist)")
ax.plot(xg, spl(xg), color=GREEN, lw=1.8, ls=":", label="smoothing spline (frequentist)")
ax.set_xlabel("time (ms)"); ax.set_ylabel("acceleration (g)"); ax.set_title("GP vs loess vs smoothing spline")
ax.legend(frameon=False); plt.tight_layout(); plt.show()
print("All three trace essentially the same curve -- unsurprising, since the smoothing spline IS a GP posterior")
print("mean. The GP adds what the others lack: a principled uncertainty band and hyperparameters chosen by")
print("marginal likelihood rather than a hand-tuned span or penalty.")

All three trace essentially the same curve -- unsurprising, since the smoothing spline IS a GP posterior
mean. The GP adds what the others lack: a principled uncertainty band and hyperparameters chosen by
marginal likelihood rather than a hand-tuned span or penalty.

4. Mauna Loa CO₂ — building structure into the kernel

Kernels add and multiply, and that is how a GP encodes structure. The Mauna Loa atmospheric CO₂ record (monthly, 1958–2001) is a rising trend with a yearly seasonal cycle on top. We compose a kernel from pieces — a long RBF (trend) $+$ a periodic$\times$RBF (a seasonal cycle allowed to drift) $+$ a short RBF (small wiggles) $+$ noise — fit it, then decompose the posterior into its trend and seasonal parts and forecast past the data.

c=pd.read_csv("co2.csv"); xc=c.year.to_numpy(); yc=c.co2.to_numpy(); yc0=yc-yc.mean()
def buildc(t):
    p=np.exp(t)
    trend  = lambda A,B: G.k_rbf(A,B,p[0],p[1])
    season = lambda A,B: G.k_periodic(A,B,p[2],1.0,p[3])*G.k_rbf(A,B,p[4],1.0)
    short  = lambda A,B: G.k_rbf(A,B,p[5],p[6])
    return (lambda A,B: trend(A,B)+season(A,B)+short(A,B)), p[7], trend, season
build_full = lambda t: (buildc(t)[0], buildc(t)[1])
thc,lmlc = G.fit_gp(xc, yc0, build_full, np.log([60,3000, 1.5,50, 120, 1.0,1.0, 0.05]))
_,nv,trend,season = buildc(thc)
# forecast + component decomposition (posterior mean of a component = K_comp(Xf,X) alpha)
from scipy.linalg import cholesky, cho_solve
kfun,_=build_full(thc); K=kfun(xc,xc)+(nv+1e-8)*np.eye(len(xc)); L=cholesky(K,lower=True); alpha=cho_solve((L,True),yc0)
xf=np.linspace(1958,2012,600); mf,Cf=G.predict_gp(thc,xc,yc0,xf,build_full); mf=mf+yc.mean(); sdf=np.sqrt(np.clip(np.diag(Cf),0,None))
tr = trend(xf,xc)@alpha + yc.mean(); se = season(xf,xc)@alpha
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].scatter(xc,yc,s=6,color=GREY,alpha=.5,label="monthly CO2")
ax[0].fill_between(xf, mf-1.96*sdf, mf+1.96*sdf, color=BLUE, alpha=.2)
ax[0].plot(xf, mf, color=BLUE, lw=1.6, label="GP fit + forecast"); ax[0].axvline(xc.max(), color="k", ls=":", lw=1)
ax[0].set_xlabel("year"); ax[0].set_ylabel("CO2 (ppm)"); ax[0].set_title("Fit and forecast (band widens past the data)"); ax[0].legend(frameon=False)
ax[1].plot(xf, tr, color=RED, lw=2, label="trend component"); ax[1].plot(xf, yc.mean()+se, color=GREEN, lw=1, label="mean + seasonal")
ax[1].set_xlabel("year"); ax[1].set_ylabel("CO2 (ppm)"); ax[1].set_title("Kernel decomposition: trend vs seasonal")
ax[1].legend(frameon=False); plt.tight_layout(); plt.show()
# zoom into the forecast region so the continuing seasonal cycle and widening band are visible
zm = xf >= 1996; dm = xc >= 1996
fig,ax=plt.subplots(figsize=(9,4))
ax.fill_between(xf[zm], (mf-1.96*sdf)[zm], (mf+1.96*sdf)[zm], color=BLUE, alpha=.2, label="95% band")
ax.plot(xf[zm], mf[zm], color=BLUE, lw=1.8, label="GP mean + forecast")
ax.scatter(xc[dm], yc[dm], s=18, color=GREY, alpha=.8, label="observed (ends 2001)")
ax.axvline(xc.max(), color="k", ls=":", lw=1.3); ax.text(xc.max()+0.15, yc[dm].min(), "forecast →", fontsize=9)
ax.set_xlim(1996,2012); ax.set_xlabel("year"); ax.set_ylabel("CO2 (ppm)")
ax.set_title("Zoom on the forecast: the annual cycle continues and the band widens with lead time")
ax.legend(frameon=False, loc="upper left"); plt.tight_layout(); plt.show()
# Report the FITTED seasonal component, not the kernel's prior scale sqrt(eta): the periodic term
# is multiplied by a slow RBF decay and the posterior component is K_season(x*,X) alpha, so the
# prior variance is not the realised amplitude -- quoting it would misstate the cycle twofold.
win = (xf > xc.max() - 5) & (xf <= xc.max())
ptp = se[win].max() - se[win].min()
print(f"fitted trend length-scale {np.exp(thc[0]):.0f} yr; seasonal period pinned to 1.00 yr.")
print(f"fitted seasonal component over the last 5 years of data: {ptp:.1f} ppm peak-to-trough,")
print("which matches the '~6 ppm yearly cycle' the R notebook reports from an independent GAM.\n")

# Does the widening band actually COVER the truth? Ten years past the data, check rather than assert.
i2012 = np.argmin(np.abs(xf-2012)); ACTUAL = 393.8      # Mauna Loa 2012 annual mean
lo, hi = mf[i2012]-1.96*sdf[i2012], mf[i2012]+1.96*sdf[i2012]
inside = lo <= ACTUAL <= hi
print(f"forecast for 2012: {mf[i2012]:.1f} ppm, 95% band [{lo:.1f}, {hi:.1f}]. Actual: {ACTUAL} ppm.")
print(f"Inside the band? {inside}.  Error {mf[i2012]-ACTUAL:+.1f} ppm = {(mf[i2012]-ACTUAL)/sdf[i2012]:+.1f} standard deviations.\n")
print("So the band grows -- and still misses, badly. This is the most useful thing in the notebook, so")
print("it is worth being precise about the cause rather than filing it under 'forecasting is hard'.")
q = np.polyfit(xc, yc, 2); l = np.polyfit(xc, yc, 1)
print(f"   quadratic OLS on the same window forecasts {np.polyval(q,2012):.1f} ppm  (error {np.polyval(q,2012)-ACTUAL:+.1f})")
print(f"   linear    OLS on the same window forecasts {np.polyval(l,2012):.1f} ppm  (error {np.polyval(l,2012)-ACTUAL:+.1f})")
print("   the GP                                     %.1f ppm  (error %+.1f)" % (mf[i2012], mf[i2012]-ACTUAL))
print("\nThe GP is sitting with the LINEAR model, and that is a direct consequence of its kernel. An RBF")
print("trend term reverts toward the prior mean once it runs out of data, so it carries the local slope")
print("forward and no curvature; atmospheric CO2, meanwhile, is accelerating. The posterior band is")
print("computed UNDER the assumed kernel -- it prices uncertainty about the function within that class,")
print("and is silent about the class itself being wrong. That is why a 5.8-sigma miss is not a paradox.")
print("Add a linear or quadratic kernel to the trend and the extrapolation changes completely; the band")
print("will still not cover the possibility that you picked the wrong kernel. Marginal likelihood chooses")
print("hyperparameters, not assumptions.")

fitted trend length-scale 56 yr; seasonal period pinned to 1.00 yr.
fitted seasonal component over the last 5 years of data: 6.6 ppm peak-to-trough,
which matches the '~6 ppm yearly cycle' the R notebook reports from an independent GAM.

forecast for 2012: 385.9 ppm, 95% band [383.3, 388.6]. Actual: 393.8 ppm.
Inside the band? False.  Error -7.9 ppm = -5.8 standard deviations.

So the band grows -- and still misses, badly. This is the most useful thing in the notebook, so
it is worth being precise about the cause rather than filing it under 'forecasting is hard'.
   quadratic OLS on the same window forecasts 392.4 ppm  (error -1.4)
   linear    OLS on the same window forecasts 382.3 ppm  (error -11.5)
   the GP                                     385.9 ppm  (error -7.9)

The GP is sitting with the LINEAR model, and that is a direct consequence of its kernel. An RBF
trend term reverts toward the prior mean once it runs out of data, so it carries the local slope
forward and no curvature; atmospheric CO2, meanwhile, 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 itself being wrong. That is why a 5.8-sigma miss is not a paradox.
Add a linear or quadratic kernel to the trend and the extrapolation changes completely; the band
will still not cover the possibility that you picked the wrong kernel. Marginal likelihood chooses
hyperparameters, not assumptions.

The composed kernel separates cleanly into a smooth rising trend and a steady ≈6 ppm yearly seasonal swing, and the forecast continues both while the credible band widens past the last observation — the zoom makes the continuing sawtooth and the growing uncertainty explicit. The point forecast is a touch conservative (the RBF trend mean-reverts slightly under extrapolation), a known and visible property, not a hidden failure. Structure lived entirely in the kernel; nothing about "trend plus seasonality" was hard-coded into a functional form.

The same Keeling curve is modelled a different way in the time-series arc (Seasonal ARIMA (SARIMA)), where a SARIMA / structural state-space model captures the trend and seasonality through differencing rather than a kernel. The two are worth reading together: a frequentist Box–Jenkins decomposition versus a Bayesian-nonparametric one on identical data.

5. Cross-check in PyMC — gp.Marginal

PyMC's gp.Marginal implements the same conjugate GP regression (the latent function is marginalised, so it runs cleanly on Windows without scan). We fit the motorcycle data with a Matérn-5/2 kernel and compare the MAP hyperparameters and predictive mean to the from-scratch result. For large series the Hilbert-space approximation gp.HSGP gives the same fit far faster.

import pymc as pm
with pm.Model() as mod:
    ell = pm.Gamma("ell", 2, 0.2); eta = pm.HalfNormal("eta", 60); sig = pm.HalfNormal("sig", 30)
    cov = eta**2 * pm.gp.cov.Matern52(1, ls=ell)
    gpm = pm.gp.Marginal(cov_func=cov)
    gpm.marginal_likelihood("y", X=xt[:,None], y=yt-ym, sigma=sig)
    mp = pm.find_MAP(progressbar=False)
    mu_pm, var_pm = gpm.predict(xg[:,None], point=mp, diag=True, pred_noise=False)
mu_pm = mu_pm + ym
print(f"PyMC MAP: length-scale {float(mp['ell']):.2f} ms, signal sd {float(mp['eta']):.1f} g, noise sd {float(mp['sig']):.1f} g")
print(f"from-scratch: length-scale {np.exp(th[0]):.2f} ms, signal sd {np.sqrt(np.exp(th[1])):.1f} g, noise sd {np.sqrt(np.exp(th[2])):.1f} g")
print(f"predictive-mean agreement: max |PyMC - from-scratch| = {np.abs(mu_pm-m).max():.2f} g")
fig,ax=plt.subplots(figsize=(9,4.2))
ax.scatter(xt,yt,s=14,color=GREY,alpha=.6,label="data")
ax.plot(xg, m, color=BLUE, lw=2.4, label="from-scratch GP")
ax.plot(xg, mu_pm, color=RED, lw=1.5, ls="--", label="PyMC gp.Marginal")
ax.set_xlabel("time (ms)"); ax.set_ylabel("acceleration (g)"); ax.set_title("From-scratch vs PyMC GP"); ax.legend(frameon=False)
plt.tight_layout(); plt.show()
print("Same hyperparameters, same curve: the closed-form GP and PyMC's gp.Marginal are one model.")

g++ not available, if using conda: `conda install gxx`

PyMC MAP: length-scale 6.27 ms, signal sd 43.1 g, noise sd 22.5 g
from-scratch: length-scale 6.55 ms, signal sd 45.7 g, noise sd 22.6 g
predictive-mean agreement: max |PyMC - from-scratch| = 0.35 g

Same hyperparameters, same curve: the closed-form GP and PyMC's gp.Marginal are one model.

6. Summary

A Gaussian process is a prior over functions: pick a kernel, and regression is exact — a posterior mean and a calibrated uncertainty band, with the smoothness/noise learned from the data by marginal likelihood rather than by cross-validation. On the motorcycle benchmark the GP found the dip-and-rebound with no functional form imposed (and honestly exposed its one weakness, the constant-noise assumption on heteroscedastic data). On Mauna Loa CO₂, kernel addition and multiplication built a structured prior — trend + drifting seasonal + wiggle + noise — that decomposed cleanly and forecast with a band that widens beyond the data.

Its frequentist counterparts, loess and the cubic smoothing spline, trace essentially the same curve; indeed the smoothing spline is a GP posterior mean, with the roughness penalty standing in for the noise-to-signal ratio — the Bayesian and penalised views are two faces of one estimator. The from-scratch Cholesky GP and PyMC's gp.Marginal agree exactly.

This is the workhorse of Bayesian nonparametric regression. Next we let the GP output something other than a real number — a latent GP through a link gives GP classification and log-Gaussian Cox processes for point patterns, connecting this arc back to the logit and counts models.