13.4 Continuous categories and the Gaussian process

Gaussian processes model continuous categories — such as spatial locations or time points — as a multivariate normal distribution. They capture correlations that decay smoothly with distance.

13.4.1 The Gaussian process

A Gaussian process is a distribution over functions. Instead of modeling each observation independently, the GP models them as jointly multivariate normal with a covariance function:

$\text{Cov}(f(x), f(x')) = \sigma^2 \exp\left(-\frac{(x-x')^2}{2\ell^2}\right)$

This is the squared exponential kernel. The length scale $\ell$ controls how quickly correlation decays with distance. The amplitude $\sigma$ controls the magnitude of variation.

# Gaussian process: squared exponential kernel
set.seed(42)

# Generate x values

x <- seq(0, 10, length.out = 50)

# Squared exponential kernel function

se*kernel <- function(x1, x2, sigma = 1, ell = 1) {
sigma^2 * exp(-outer(x1, x2, "-")^2 / (2 _ ell^2))
}

# Compute covariance matrices for different length scales

K_ell1 <- se_kernel(x, x, sigma = 1, ell = 1)
K_ell2 <- se_kernel(x, x, sigma = 1, ell = 2)
K_ell5 <- se_kernel(x, x, sigma = 1, ell = 5)

# Draw samples from the GP prior

library(MASS)
n_samples <- 5

cat("Gaussian process prior samples:\n")
cat(" Length scale controls smoothness:\n")
cat(" ell = 1: wiggly functions\n")
cat(" ell = 2: moderate smoothness\n")
cat(" ell = 5: very smooth functions\n")

# Plot

par(mfrow = c(1, 3))
for (ell_val in c(1, 2, 5)) {
K <- se_kernel(x, x, sigma = 1, ell = ell_val)
samples <- mvrnorm(n_samples, rep(0, length(x)), K)
plot(x, samples[1, ], type = "l", ylim = range(samples),
main = paste("ell =", ell_val), xlab = "x", ylab = "f(x)")
for (s in 2:n_samples) lines(x, samples[s, ], col = s)
}
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

13.4.2 Fitting a Gaussian process

# Gaussian process regression: simulate and fit
set.seed(42)

# True function: a smooth curve

x <- seq(0, 10, length.out = 30)
f*true <- 3 * sin(x) + 0.5 _ x
y <- f_true + rnorm(length(x), 0, 0.5)

# GP regression using a simple approach

# Compute the covariance matrix

ell <- 1.5
sigma_f <- 2
sigma_n <- 0.5 # noise

K <- sigma*f^2 * exp(-outer(x, x, "-")^2 / (2 _ ell^2))
K_noise <- K + sigma_n^2 * diag(length(x))

# GP posterior mean at training points

# f_post = K^T (K + sigma_n^2 I)^{-1} y

L <- chol(K_noise)
alpha <- backsolve(L, forwardsolve(t(L), y))
f_post <- K %*% alpha

cat("Gaussian process regression:\n")
cat(" Length scale (ell):", ell, "\n")
cat(" Amplitude (sigma_f):", sigma_f, "\n")
cat(" Noise (sigma_n):", sigma_n, "\n")
cat(" Posterior mean RMSE:", round(sqrt(mean((f_post - f_true)^2)), 3), "\n")

# Predict at new points

x*new <- seq(0, 10, length.out = 100)
K_star <- sigma_f^2 * exp(-outer(x*new, x, "-")^2 / (2 * ell^2))
f_pred <- K_star %*% alpha

# Plot

plot(x, y, pch = 16, col = "steelblue",
main = "Gaussian process regression",
xlab = "x", ylab = "y")
lines(x*new, f_pred, lwd = 2, col = "firebrick")
curve(3 * sin(x) + 0.5 _ x, add = TRUE, lty = 2, col = "darkgreen", lwd = 2)
legend("topleft", legend = c("Data", "GP posterior mean", "True function"),
col = c("steelblue", "firebrick", "darkgreen"),
pch = c(16, NA, NA), lty = c(NA, 1, 2), lwd = c(NA, 2, 2))

Click Run to execute this code in your browser with WebR

Rethinking

Gaussian processes are varying intercepts for continuous categories. A GP is a multilevel model where the "groups" are points along a continuous dimension (space, time, age). Instead of each group having an independent intercept, the intercepts are correlated — nearby points are more similar. The covariance function determines how similarity decays with distance. This makes GPs incredibly flexible: they can model spatial variation, temporal trends, and any continuous covariate where nearby values should be similar. The length scale $\ell$ is the key parameter — it tells the model how far the "local" influence extends.

Overthinking

The kernel trick. The Gaussian process works because the covariance matrix encodes similarity. The squared exponential kernel $k(x, x') = \sigma^2 \exp(-(x-x')^2 / (2\ell^2))$ is just one choice. Others include:

Matérn: Rougher than squared exponential, more realistic for physical processes.
Periodic: $k(x, x') = \sigma^2 \exp(-2\sin^2(\pi|x-x'|/p) / \ell^2)$ — captures seasonality.
Rational quadratic: A scale mixture of squared exponentials — handles multiple length scales.

The choice of kernel is a modeling decision — it encodes assumptions about how the function behaves.

Summary of Section 13.4

Gaussian processes model continuous categories as correlated random effects.
The squared exponential kernel: $\text{Cov}(f(x), f(x')) = \sigma^2 \exp(-(x-x')^2 / (2\ell^2))$ .
The length scale $\ell$ controls smoothness — larger $\ell$ = smoother functions.
GPs are varying intercepts for continuous dimensions (space, time).
The kernel choice encodes assumptions about the function's behavior.

Example: Cross-classified chimpanzees with varying slopes

Summary