Gaussian process regression from scratch

regression

Gaussian processes

ecology tutorial

Fitting a Gaussian process in base R along an ecological gradient: the RBF kernel, marginal-likelihood hyperparameters, the posterior mean and credible band, and honest length-scale limits.

Author

Tidy Ecology

Published

2026-08-10

Suppose you have measured a response along an environmental gradient: an abundance index across a moisture gradient, say, or leaf nitrogen across an elevation transect. The shape is clearly smooth but you have no mechanistic equation for it, the sampling is uneven, and there is a stretch of the gradient nobody visited. You want a flexible fit, and you want the uncertainty to grow honestly where the data thins out.

A Gaussian process (GP) does exactly this. Like the penalised splines behind a generalised additive model, a GP is a flexible smoother, but it is built from a covariance function rather than a basis, and it returns a principled predictive variance at every point. The same machinery, taken onto a map, is the engine behind kriging. Here we fit one from first principles in base R, show that its posterior mean is nothing more exotic than a penalised weighted average, and then look squarely at the part that is easy to gloss over: the length-scale that controls all of this is only weakly pinned down by a single sample.

A distribution over functions

A Gaussian process places a prior directly on the unknown function. The defining property is that the function values at any finite set of locations follow a multivariate normal distribution, with the covariance between two points set by a kernel. We use the squared-exponential (radial basis) kernel,

\[k(x, x') = \sigma_f^2 \, \exp\!\left(-\frac{(x - x')^2}{2\,\ell^2}\right),\]

which says two points are correlated when they are close and independent when they are far apart. The length-scale $\ell$ controls how quickly that correlation decays (how wiggly the function is allowed to be), and the signal variance $\sigma_f^2$ sets its vertical scale.

We simulate a smooth but non-trivial response on a standardised gradient, with a deliberate unsampled gap so that the uncertainty has somewhere to grow.

set.seed(4218)
f_true <- function(x) 1.2*exp(-((x-0.35)^2)/(2*0.09)) - 0.7*exp(-((x-0.78)^2)/(2*0.03)) + 0.25
sig_n_true <- 0.18
gap <- c(0.52, 0.68)                        # an unsampled stretch of the gradient
x  <- sort(runif(58, 0, 1)); x <- x[x < gap[1] | x > gap[2]]
n  <- length(x)
y  <- f_true(x) + rnorm(n, 0, sig_n_true)
ybar <- mean(y); yc <- y - ybar             # centre; the GP prior mean is 0 on the residual

Conditioning on the data

The observations are the function plus independent noise, $y_i = f(x_i) + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)$. Because everything is jointly normal, the posterior over the function at any set of new locations $x_\ast$ is again normal, with a closed form. Writing $K$ for the kernel matrix on the training points and $K_\ast$ for the cross-covariances to the new points, the posterior mean and variance are

\[\bar f_\ast = K_\ast (K + \sigma_n^2 I)^{-1} y, \qquad \operatorname{Var}(f_\ast) = \sigma_f^2 - K_\ast (K + \sigma_n^2 I)^{-1} K_\ast^{\top}.\]

We never invert the matrix directly. A Cholesky factor gives the same answer with far better numerical behaviour, and it also hands us the log determinant we will need for fitting.

D  <- as.matrix(dist(x))                     # n x n pairwise |x_i - x_j|
kern <- function(D, sf, ell) sf^2 * exp(-(D^2)/(2*ell^2))

nll <- function(par, D, yc) {                # negative log marginal likelihood
  ell <- exp(par[1]); sf <- exp(par[2]); sn <- exp(par[3])
  K   <- kern(D, sf, ell); Ky <- K + diag(sn^2, nrow(K))
  R   <- tryCatch(chol(Ky), error = function(e) NULL)
  if (is.null(R)) return(1e10)
  a   <- backsolve(R, forwardsolve(t(R), yc))
  0.5*sum(yc*a) + sum(log(diag(R))) + 0.5*length(yc)*log(2*pi)
}

predict_gp <- function(xs, x, yc, ybar, sf, ell, sn) {
  K  <- kern(as.matrix(dist(x)), sf, ell) + diag(sn^2, length(x))
  R  <- chol(K); a <- backsolve(R, forwardsolve(t(R), yc))
  Ds <- outer(xs, x, function(u, v) abs(u - v)); Ks <- kern(Ds, sf, ell)
  mu <- as.numeric(Ks %*% a) + ybar
  v  <- forwardsolve(t(R), t(Ks))
  var_f <- sf^2 - colSums(v^2); var_f[var_f < 0] <- 0
  list(mu = mu, sd_f = sqrt(var_f), sd_y = sqrt(var_f + sn^2))
}

Fitting the hyperparameters

The length-scale, signal variance and noise are not set by hand. They are chosen by maximising the marginal likelihood of the data, which after integrating the function out is a plain multivariate normal, $y \sim \mathcal{N}(0,\, K + \sigma_n^2 I)$. Its log density,

\[\log p(y) = -\tfrac12 y^{\top}(K + \sigma_n^2 I)^{-1} y - \tfrac12 \log\lvert K + \sigma_n^2 I\rvert - \tfrac{n}{2}\log 2\pi,\]

balances fit (the first term) against complexity (the log determinant), so it will not simply chase the noise. We optimise on the log scale to keep every parameter positive.

fit <- optim(c(log(0.15), log(sd(yc)), log(0.2)), nll, D = D, yc = yc,
             method = "BFGS", control = list(maxit = 500))
ell_hat <- exp(fit$par[1]); sf_hat <- exp(fit$par[2]); sn_hat <- exp(fit$par[3])
logML   <- -fit$value
xg  <- seq(-0.15, 1.15, length.out = 300)    # a fine grid, with a little extrapolation
pg  <- predict_gp(xg, x, yc, ybar, sf_hat, ell_hat, sn_hat)
inr <- xg >= min(x) & xg <= max(x)
rmse_gp <- sqrt(mean((pg$mu[inr] - f_true(xg[inr]))^2))

The fit recovers a length-scale of 0.24, a signal SD of 0.59, and a noise SD of 0.19, which is close to the true 0.18. The posterior mean tracks the truth well over the observed span (root mean squared error 0.070), and the band does the useful thing: it stays tight where points are dense and swells across the gap and at the edges, because the kernel has little to lean on there.

A green curve with a shaded band tracks scattered points along a gradient; the band is narrow where points are dense and clearly wider across a shaded unsampled gap and at the two ends. — Figure 1: Gaussian process fit along the gradient. Line: posterior mean; ribbon: 95% credible band for the function; dashed: the true response; points: observations. The band widens across the unsampled gap.

What the posterior mean actually is

It is tempting to treat a GP as a black box, but the posterior mean is transparent. Collect the weights $\alpha = (K + \sigma_n^2 I)^{-1} y$; then $\bar f_\ast = \sum_i \alpha_i\, k(x_\ast, x_i)$ is a weighted sum of kernels centred on the data. Pull the signal variance out of the kernel and the mean is exactly kernel ridge regression on the correlation kernel $\tilde K$, penalised by the noise-to-signal ratio $\lambda = \sigma_n^2 / \sigma_f^2$.

lambda_kr <- sn_hat^2 / sf_hat^2                 # noise-to-signal ratio = ridge penalty
Ktil   <- exp(-(D^2)/(2*ell_hat^2))              # correlation kernel (signal SD set to 1)
a_kr   <- solve(Ktil + diag(lambda_kr, n), yc)
Ds_g   <- outer(xg, x, function(u, v) abs(u - v)); Ktil_s <- exp(-(Ds_g^2)/(2*ell_hat^2))
mu_kr  <- as.numeric(Ktil_s %*% a_kr) + ybar
max_gp_vs_kr <- max(abs(pg$mu - mu_kr))          # agreement to machine precision

The two predictions agree to machine precision (largest difference 2.4^{-15}), with penalty $\lambda =$ 0.10. This is worth holding onto for three reasons. It says a GP is a penalised smoother, the same family as the spline in a GAM, where the length-scale plays the role of the smoothing parameter; the geoadditive view (Kammann and Wand 2003) makes the link precise by writing both as a mixed model. It says the posterior mean is exactly the simple-kriging predictor used on maps in the kriging tutorial. And it says the length-scale and the signal-to-noise split are the levers that matter, which is where the honesty comes in.

The length-scale is only weakly identified

The marginal likelihood chose a single length-scale, but it did not choose it sharply. If we profile the marginal likelihood over $\ell$, re-optimising the variances at each value, the curve has a broad rather than a peaked maximum.

ell_grid <- exp(seq(log(0.03), log(0.60), length.out = 60))
prof <- sapply(ell_grid, function(le) {
  o <- optim(c(log(sd(yc)), log(0.2)),
             function(p) nll(c(log(le), p[1], p[2]), D, yc),
             method = "BFGS", control = list(maxit = 300))
  -o$value
})
ell_mle <- ell_grid[which.max(prof)]
within2 <- ell_grid[prof >= max(prof) - 2]       # length-scales within 2 log-ML of the best
ell_lo  <- min(within2); ell_hi <- max(within2); ratio_flat <- ell_hi / ell_lo

refit_at <- function(le) {
  o <- optim(c(log(sd(yc)), log(0.2)),
             function(p) nll(c(log(le), p[1], p[2]), D, yc),
             method = "BFGS", control = list(maxit = 300))
  list(ell = le, sf = exp(o$par[1]), sn = exp(o$par[2]))
}
plo <- refit_at(ell_lo); phi <- refit_at(ell_hi)
glo <- predict_gp(xg, x, yc, ybar, plo$sf, plo$ell, plo$sn)
ghi <- predict_gp(xg, x, yc, ybar, phi$sf, phi$ell, phi$sn)
dense <- inr & (xg < gap[1] | xg > gap[2])
maxdiff_dense <- max(abs(glo$mu[dense] - ghi$mu[dense]))
ig <- which.min(abs(xg - mean(gap)))
diff_gap    <- abs(glo$mu[ig] - ghi$mu[ig])
band_lo_gap <- 2*1.96*glo$sd_f[ig]; band_hi_gap <- 2*1.96*ghi$sd_f[ig]

Every length-scale between 0.14 and 0.34 sits within two log-likelihood units of the best value at 0.24: a factor of 2.4 in the quantity that governs the whole fit. Take the two ends of that plausible range and the consequence is clean. Where the points are dense the two fits are indistinguishable (largest difference 0.06), because the data pins the function down. In the gap they part company: their means differ by 0.10, and the shorter length-scale reports a 95% band nearly twice as wide as the longer one (0.47 against 0.24). The uncertainty you quote in the region with no data is set by a hyperparameter the data barely constrains.

Left: a broad-topped log-likelihood curve against length-scale with a shaded plateau. Right: two fitted curves overlapping where points are dense but separating across a shaded gap, each with its own band. — Figure 2: Weak identifiability of the length-scale. (A) The profile marginal log-likelihood is broad; the shaded band holds all length-scales within two units of the best. (B) Two length-scales from that band agree where data is dense but diverge across the gap, with different band widths.

This is not a defect of the code; it is intrinsic. For the Matern family, of which the squared-exponential is a limit, Zhang (2004) showed that under denser and denser sampling of a fixed region the variance and the range cannot both be estimated consistently: only a single combination of them, the one that matters for interpolation, is pinned down. A convergent fit and a tidy point estimate do not imply that the length-scale is identified.

Where this leaves you

A Gaussian process is a genuinely useful flexible smoother: it fits an unknown response without a functional form, and its band grows honestly where the data is sparse, which a bare spline fit does not advertise. But treat the length-scale and the signal-to-noise split as weakly identified quantities rather than settled numbers. The plug-in band conditions on a single hyperparameter estimate, so in data-poor regions it understates the true spread; propagating hyperparameter uncertainty, or at least checking how the fit moves across the plausible range, is the honest habit. And the plug-in variance should not be the only thing you trust: out-of-sample checks, which we take up in checking a predictive model, are what keep it honest. That is the same move that recurs across this site, from the E-value to honest quantile tails: name the modelling choice, and do not read a plug-in interval as the last word.

A Gaussian process is one flexible smoother among several. Local regression reaches a similar fit from a different direction, weighting nearby points rather than placing a prior over functions, and it makes the bias hidden in every smoother easy to see. And when a calibrated interval matters more than a model-based one, conformal prediction wraps any of these fits in a band with a coverage guarantee that does not depend on the length-scale being right. The same idea carried onto a map is kriging: there the distance is measured on the ground, but everything above carries over unchanged.

References

Rasmussen CE, Williams CKI 2006. Gaussian Processes for Machine Learning. MIT Press. ISBN 978-0-262-18253-9.

Zhang H 2004. Journal of the American Statistical Association 99(465):250-261 (10.1198/016214504000000241).

Kammann EE, Wand MP 2003. Journal of the Royal Statistical Society Series C 52(1):1-18 (10.1111/1467-9876.00385).

Dormann CF, McPherson JM, Araujo MB, and others 2007. Ecography 30(5):609-628 (10.1111/j.2007.0906-7590.05171.x).

--- title: "Gaussian process regression from scratch" description: "Fitting a Gaussian process in base R along an ecological gradient: the RBF kernel, marginal-likelihood hyperparameters, the posterior mean and credible band, and honest length-scale limits." date: "2026-08-10 09:00" categories: [R, regression, Gaussian processes, ecology tutorial] image: thumbnail.png image-alt: "A Gaussian process fit to points along a gradient, with a shaded credible band that widens across an unsampled gap." --- Suppose you have measured a response along an environmental gradient: an abundance index across a moisture gradient, say, or leaf nitrogen across an elevation transect. The shape is clearly smooth but you have no mechanistic equation for it, the sampling is uneven, and there is a stretch of the gradient nobody visited. You want a flexible fit, and you want the uncertainty to grow honestly where the data thins out. A Gaussian process (GP) does exactly this. Like the penalised splines behind a [generalised additive model](../penalised-regression-splines/), a GP is a flexible smoother, but it is built from a covariance function rather than a basis, and it returns a principled predictive variance at every point. The same machinery, taken onto a map, is the engine behind [kriging](../kriging-spatial-interpolation/). Here we fit one from first principles in base R, show that its posterior mean is nothing more exotic than a penalised weighted average, and then look squarely at the part that is easy to gloss over: the length-scale that controls all of this is only weakly pinned down by a single sample. ## A distribution over functions A Gaussian process places a prior directly on the unknown function. The defining property is that the function values at any finite set of locations follow a multivariate normal distribution, with the covariance between two points set by a kernel. We use the squared-exponential (radial basis) kernel, $$k(x, x') = \sigma_f^2 \, \exp\!\left(-\frac{(x - x')^2}{2\,\ell^2}\right),$$ which says two points are correlated when they are close and independent when they are far apart. The length-scale $\ell$ controls how quickly that correlation decays (how wiggly the function is allowed to be), and the signal variance $\sigma_f^2$ sets its vertical scale. We simulate a smooth but non-trivial response on a standardised gradient, with a deliberate unsampled gap so that the uncertainty has somewhere to grow. ```{r} #| label: dgp set.seed(4218) f_true <- function(x) 1.2*exp(-((x-0.35)^2)/(2*0.09)) - 0.7*exp(-((x-0.78)^2)/(2*0.03)) + 0.25 sig_n_true <- 0.18 gap <- c(0.52, 0.68) # an unsampled stretch of the gradient x <- sort(runif(58, 0, 1)); x <- x[x < gap[1] | x > gap[2]] n <- length(x) y <- f_true(x) + rnorm(n, 0, sig_n_true) ybar <- mean(y); yc <- y - ybar # centre; the GP prior mean is 0 on the residual ``` ## Conditioning on the data The observations are the function plus independent noise, $y_i = f(x_i) + \varepsilon_i$ with $\varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)$. Because everything is jointly normal, the posterior over the function at any set of new locations $x_\ast$ is again normal, with a closed form. Writing $K$ for the kernel matrix on the training points and $K_\ast$ for the cross-covariances to the new points, the posterior mean and variance are $$\bar f_\ast = K_\ast (K + \sigma_n^2 I)^{-1} y, \qquad \operatorname{Var}(f_\ast) = \sigma_f^2 - K_\ast (K + \sigma_n^2 I)^{-1} K_\ast^{\top}.$$ We never invert the matrix directly. A Cholesky factor gives the same answer with far better numerical behaviour, and it also hands us the log determinant we will need for fitting. ```{r} #| label: gp-machinery D <- as.matrix(dist(x)) # n x n pairwise |x_i - x_j| kern <- function(D, sf, ell) sf^2 * exp(-(D^2)/(2*ell^2)) nll <- function(par, D, yc) { # negative log marginal likelihood ell <- exp(par[1]); sf <- exp(par[2]); sn <- exp(par[3]) K <- kern(D, sf, ell); Ky <- K + diag(sn^2, nrow(K)) R <- tryCatch(chol(Ky), error = function(e) NULL) if (is.null(R)) return(1e10) a <- backsolve(R, forwardsolve(t(R), yc)) 0.5*sum(yc*a) + sum(log(diag(R))) + 0.5*length(yc)*log(2*pi) } predict_gp <- function(xs, x, yc, ybar, sf, ell, sn) { K <- kern(as.matrix(dist(x)), sf, ell) + diag(sn^2, length(x)) R <- chol(K); a <- backsolve(R, forwardsolve(t(R), yc)) Ds <- outer(xs, x, function(u, v) abs(u - v)); Ks <- kern(Ds, sf, ell) mu <- as.numeric(Ks %*% a) + ybar v <- forwardsolve(t(R), t(Ks)) var_f <- sf^2 - colSums(v^2); var_f[var_f < 0] <- 0 list(mu = mu, sd_f = sqrt(var_f), sd_y = sqrt(var_f + sn^2)) } ``` ## Fitting the hyperparameters The length-scale, signal variance and noise are not set by hand. They are chosen by maximising the marginal likelihood of the data, which after integrating the function out is a plain multivariate normal, $y \sim \mathcal{N}(0,\, K + \sigma_n^2 I)$. Its log density, $$\log p(y) = -\tfrac12 y^{\top}(K + \sigma_n^2 I)^{-1} y - \tfrac12 \log\lvert K + \sigma_n^2 I\rvert - \tfrac{n}{2}\log 2\pi,$$ balances fit (the first term) against complexity (the log determinant), so it will not simply chase the noise. We optimise on the log scale to keep every parameter positive. ```{r} #| label: fit fit <- optim(c(log(0.15), log(sd(yc)), log(0.2)), nll, D = D, yc = yc, method = "BFGS", control = list(maxit = 500)) ell_hat <- exp(fit$par[1]); sf_hat <- exp(fit$par[2]); sn_hat <- exp(fit$par[3]) logML <- -fit$value xg <- seq(-0.15, 1.15, length.out = 300) # a fine grid, with a little extrapolation pg <- predict_gp(xg, x, yc, ybar, sf_hat, ell_hat, sn_hat) inr <- xg >= min(x) & xg <= max(x) rmse_gp <- sqrt(mean((pg$mu[inr] - f_true(xg[inr]))^2)) ``` The fit recovers a length-scale of `r sprintf("%.2f", ell_hat)`, a signal SD of `r sprintf("%.2f", sf_hat)`, and a noise SD of `r sprintf("%.2f", sn_hat)`, which is close to the true `r sig_n_true`. The posterior mean tracks the truth well over the observed span (root mean squared error `r sprintf("%.3f", rmse_gp)`), and the band does the useful thing: it stays tight where points are dense and swells across the gap and at the edges, because the kernel has little to lean on there. ```{r} #| label: fig-fit #| echo: false #| fig-cap: "Gaussian process fit along the gradient. Line: posterior mean; ribbon: 95% credible band for the function; dashed: the true response; points: observations. The band widens across the unsampled gap." #| fig-alt: "A green curve with a shaded band tracks scattered points along a gradient; the band is narrow where points are dense and clearly wider across a shaded unsampled gap and at the two ends." suppressPackageStartupMessages(library(ggplot2)) theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#dad9ca", linewidth = 0.3), axis.title = element_text(colour = "#2c3a31"), axis.text = element_text(colour = "#5d6b61"), plot.title = element_text(colour = "#16241d", face = "bold", size = 13), plot.subtitle = element_text(colour = "#46604a", size = 10), strip.text = element_text(colour = "#275139", face = "bold"), legend.title = element_text(colour = "#2c3a31"), legend.text = element_text(colour = "#46604a"), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA)) } forest <- "#275139"; sage <- "#93a87f"; ink <- "#16241d"; terra <- "#b5534e"; gold <- "#c9b458" dfit <- data.frame(x = xg, mu = pg$mu, lo = pg$mu - 1.96*pg$sd_f, hi = pg$mu + 1.96*pg$sd_f, truth = f_true(xg)) dobs <- data.frame(x = x, y = y) ggplot() + annotate("rect", xmin = gap[1], xmax = gap[2], ymin = -Inf, ymax = Inf, fill = "#dad9ca", alpha = 0.35) + annotate("text", x = mean(gap), y = 1.85, label = "unsampled gap", colour = "#5d6b61", size = 3) + geom_ribbon(data = dfit, aes(x, ymin = lo, ymax = hi), fill = sage, alpha = 0.35) + geom_line(data = dfit, aes(x, truth), colour = terra, linetype = "dashed", linewidth = 0.7) + geom_line(data = dfit, aes(x, mu), colour = forest, linewidth = 1) + geom_point(data = dobs, aes(x, y), colour = ink, size = 1.5, alpha = 0.8) + labs(x = "environmental gradient (standardised)", y = "response") + coord_cartesian(ylim = c(-0.7, 1.95)) + theme_te() ``` ## What the posterior mean actually is It is tempting to treat a GP as a black box, but the posterior mean is transparent. Collect the weights $\alpha = (K + \sigma_n^2 I)^{-1} y$; then $\bar f_\ast = \sum_i \alpha_i\, k(x_\ast, x_i)$ is a weighted sum of kernels centred on the data. Pull the signal variance out of the kernel and the mean is exactly kernel ridge regression on the correlation kernel $\tilde K$, penalised by the noise-to-signal ratio $\lambda = \sigma_n^2 / \sigma_f^2$. ```{r} #| label: equivalence lambda_kr <- sn_hat^2 / sf_hat^2 # noise-to-signal ratio = ridge penalty Ktil <- exp(-(D^2)/(2*ell_hat^2)) # correlation kernel (signal SD set to 1) a_kr <- solve(Ktil + diag(lambda_kr, n), yc) Ds_g <- outer(xg, x, function(u, v) abs(u - v)); Ktil_s <- exp(-(Ds_g^2)/(2*ell_hat^2)) mu_kr <- as.numeric(Ktil_s %*% a_kr) + ybar max_gp_vs_kr <- max(abs(pg$mu - mu_kr)) # agreement to machine precision ``` The two predictions agree to machine precision (largest difference `r signif(max_gp_vs_kr, 2)`), with penalty $\lambda =$ `r sprintf("%.2f", lambda_kr)`. This is worth holding onto for three reasons. It says a GP is a penalised smoother, the same family as the spline in a GAM, where the length-scale plays the role of the smoothing parameter; the [geoadditive](../penalised-regression-splines/) view (Kammann and Wand 2003) makes the link precise by writing both as a mixed model. It says the posterior mean is exactly the simple-kriging predictor used on maps in [the kriging tutorial](../kriging-spatial-interpolation/). And it says the length-scale and the signal-to-noise split are the levers that matter, which is where the honesty comes in. ## The length-scale is only weakly identified The marginal likelihood chose a single length-scale, but it did not choose it sharply. If we profile the marginal likelihood over $\ell$, re-optimising the variances at each value, the curve has a broad rather than a peaked maximum. ```{r} #| label: profile ell_grid <- exp(seq(log(0.03), log(0.60), length.out = 60)) prof <- sapply(ell_grid, function(le) { o <- optim(c(log(sd(yc)), log(0.2)), function(p) nll(c(log(le), p[1], p[2]), D, yc), method = "BFGS", control = list(maxit = 300)) -o$value }) ell_mle <- ell_grid[which.max(prof)] within2 <- ell_grid[prof >= max(prof) - 2] # length-scales within 2 log-ML of the best ell_lo <- min(within2); ell_hi <- max(within2); ratio_flat <- ell_hi / ell_lo refit_at <- function(le) { o <- optim(c(log(sd(yc)), log(0.2)), function(p) nll(c(log(le), p[1], p[2]), D, yc), method = "BFGS", control = list(maxit = 300)) list(ell = le, sf = exp(o$par[1]), sn = exp(o$par[2])) } plo <- refit_at(ell_lo); phi <- refit_at(ell_hi) glo <- predict_gp(xg, x, yc, ybar, plo$sf, plo$ell, plo$sn) ghi <- predict_gp(xg, x, yc, ybar, phi$sf, phi$ell, phi$sn) dense <- inr & (xg < gap[1] | xg > gap[2]) maxdiff_dense <- max(abs(glo$mu[dense] - ghi$mu[dense])) ig <- which.min(abs(xg - mean(gap))) diff_gap <- abs(glo$mu[ig] - ghi$mu[ig]) band_lo_gap <- 2*1.96*glo$sd_f[ig]; band_hi_gap <- 2*1.96*ghi$sd_f[ig] ``` Every length-scale between `r sprintf("%.2f", ell_lo)` and `r sprintf("%.2f", ell_hi)` sits within two log-likelihood units of the best value at `r sprintf("%.2f", ell_mle)`: a factor of `r sprintf("%.1f", ratio_flat)` in the quantity that governs the whole fit. Take the two ends of that plausible range and the consequence is clean. Where the points are dense the two fits are indistinguishable (largest difference `r sprintf("%.2f", maxdiff_dense)`), because the data pins the function down. In the gap they part company: their means differ by `r sprintf("%.2f", diff_gap)`, and the shorter length-scale reports a 95% band nearly twice as wide as the longer one (`r sprintf("%.2f", band_lo_gap)` against `r sprintf("%.2f", band_hi_gap)`). The uncertainty you quote in the region with no data is set by a hyperparameter the data barely constrains. ```{r} #| label: fig-identifiability #| echo: false #| fig-cap: "Weak identifiability of the length-scale. (A) The profile marginal log-likelihood is broad; the shaded band holds all length-scales within two units of the best. (B) Two length-scales from that band agree where data is dense but diverge across the gap, with different band widths." #| fig-alt: "Left: a broad-topped log-likelihood curve against length-scale with a shaded plateau. Right: two fitted curves overlapping where points are dense but separating across a shaded gap, each with its own band." dprof <- data.frame(ell = ell_grid, lml = prof) pA <- ggplot(dprof, aes(ell, lml)) + annotate("rect", xmin = ell_lo, xmax = ell_hi, ymin = -Inf, ymax = Inf, fill = sage, alpha = 0.30) + geom_line(colour = forest, linewidth = 1) + geom_vline(xintercept = ell_mle, colour = ink, linetype = "dotted") + labs(x = "length-scale", y = "profile log marginal likelihood", subtitle = "(A) the length-scale is weakly identified") + theme_te() mk <- function(g, lab) data.frame(x = xg, mu = g$mu, lo = g$mu - 1.96*g$sd_f, hi = g$mu + 1.96*g$sd_f, which = lab) d2 <- rbind(mk(glo, sprintf("short (%.2f)", plo$ell)), mk(ghi, sprintf("long (%.2f)", phi$ell))) pB <- ggplot() + annotate("rect", xmin = gap[1], xmax = gap[2], ymin = -Inf, ymax = Inf, fill = "#dad9ca", alpha = 0.35) + geom_ribbon(data = d2, aes(x, ymin = lo, ymax = hi, fill = which), alpha = 0.20) + geom_line(data = data.frame(x = xg, truth = f_true(xg)), aes(x, truth), colour = terra, linetype = "dashed", linewidth = 0.6) + geom_line(data = d2, aes(x, mu, colour = which), linewidth = 0.9) + geom_point(data = dobs, aes(x, y), colour = ink, size = 1, alpha = 0.7) + scale_colour_manual(values = c(forest, gold), name = "length-scale") + scale_fill_manual(values = c(forest, gold), name = "length-scale") + coord_cartesian(ylim = c(-0.7, 1.95)) + labs(x = "environmental gradient", y = "response", subtitle = "(B) agree where dense, disagree in the gap") + theme_te() + theme(legend.position = c(0.80, 0.86), legend.background = element_rect(fill = "white", colour = "#dad9ca"), legend.key.size = unit(0.8, "lines")) grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2, widths = grid::unit(c(1, 1.15), "null")))) print(pA, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` This is not a defect of the code; it is intrinsic. For the Matern family, of which the squared-exponential is a limit, Zhang (2004) showed that under denser and denser sampling of a fixed region the variance and the range cannot both be estimated consistently: only a single combination of them, the one that matters for interpolation, is pinned down. A convergent fit and a tidy point estimate do not imply that the length-scale is identified. ## Where this leaves you A Gaussian process is a genuinely useful flexible smoother: it fits an unknown response without a functional form, and its band grows honestly where the data is sparse, which a bare spline fit does not advertise. But treat the length-scale and the signal-to-noise split as weakly identified quantities rather than settled numbers. The plug-in band conditions on a single hyperparameter estimate, so in data-poor regions it understates the true spread; propagating hyperparameter uncertainty, or at least checking how the fit moves across the plausible range, is the honest habit. And the plug-in variance should not be the only thing you trust: out-of-sample checks, which we take up in [checking a predictive model](../checking-predictive-calibration/), are what keep it honest. That is the same move that recurs across this site, from the [E-value](../sensitivity-to-unmeasured-confounding/) to [honest quantile tails](../checking-a-quantile-regression/): name the modelling choice, and do not read a plug-in interval as the last word. A Gaussian process is one flexible smoother among several. [Local regression](../local-regression-from-scratch/) reaches a similar fit from a different direction, weighting nearby points rather than placing a prior over functions, and it makes the bias hidden in every smoother easy to see. And when a calibrated interval matters more than a model-based one, [conformal prediction](../conformal-prediction-intervals/) wraps any of these fits in a band with a coverage guarantee that does not depend on the length-scale being right. The same idea carried onto a map is [kriging](../kriging-spatial-interpolation/): there the distance is measured on the ground, but everything above carries over unchanged. ## References Rasmussen CE, Williams CKI 2006. Gaussian Processes for Machine Learning. MIT Press. ISBN 978-0-262-18253-9. Zhang H 2004. Journal of the American Statistical Association 99(465):250-261 (10.1198/016214504000000241). Kammann EE, Wand MP 2003. Journal of the Royal Statistical Society Series C 52(1):1-18 (10.1111/1467-9876.00385). Dormann CF, McPherson JM, Araujo MB, and others 2007. Ecography 30(5):609-628 (10.1111/j.2007.0906-7590.05171.x). ## Related tutorials - [Local regression from scratch](../local-regression-from-scratch/) - [Conformal prediction intervals from scratch](../conformal-prediction-intervals/) - [Checking a predictive model: calibration and coverage](../checking-predictive-calibration/) - [Penalised regression splines from scratch](../penalised-regression-splines/)