Gaussian-Process Regression (R)¶
mgcv — GP smoothers, penalised splines, and the frequentist cousins¶
The R counterpart to gp_python.ipynb. Where the Python notebook builds the GP from scratch (kernel, Cholesky, marginal likelihood) and cross-checks it in PyMC, this notebook uses the standard R smoothing package mgcv — whose bs="gp" basis is a Gaussian process and whose penalised splines are its close cousin — alongside base R's loess and smooth.spline. On the motorcycle benchmark all trace the same curve (the smoothing spline is literally a GP posterior mean); a mgcv GAM then pulls the Mauna Loa CO₂ series into trend and seasonal parts.
options(repr.plot.width=12, repr.plot.height=4.4)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(mgcv))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"
cat("R engine: mgcv (bs='gp' Gaussian-process smoother, s() penalised spline) + base loess/smooth.spline\n")
R engine: mgcv (bs='gp' Gaussian-process smoother, s() penalised spline) + base loess/smooth.spline
1. The motorcycle data — a GP smoother and its cousins¶
133 accelerometer readings from a simulated crash. mgcv::gam(y ~ s(x, bs="gp")) fits an actual Gaussian-process smoother (with REML-chosen range), s(x) a penalised thin-plate spline, and base R's loess / smooth.spline the classical local and penalised smoothers. They should coincide — the cubic smoothing spline is a GP posterior mean.
d<-read.csv("mcycle.csv"); x<-d$times; y<-d$accel; xg<-seq(min(x),max(x),length=250)
g_gp <-gam(y~s(x, bs="gp", k=20)) # Gaussian-process smoother
g_spl<-gam(y~s(x, k=20)) # penalised thin-plate spline
lo<-loess(y~x, span=0.2); sp<-smooth.spline(x,y)
pg<-predict(g_gp, data.frame(x=xg), se.fit=TRUE)
par(mar=c(4,4,3,1))
ylim_m <- range(c(y, pg$fit - 1.96*pg$se.fit, pg$fit + 1.96*pg$se.fit))
plot(x,y,pch=19,col=GREY,cex=.6,ylim=ylim_m,xlab="time (ms)",ylab="acceleration (g)",main="Motorcycle: mgcv GP smoother vs spline / loess / smooth.spline")
polygon(c(xg,rev(xg)), c(pg$fit-1.96*pg$se.fit, rev(pg$fit+1.96*pg$se.fit)), col=adjustcolor(BLUE,.2), border=NA)
lines(xg, pg$fit, col=BLUE, lwd=2.6)
lines(xg, predict(g_spl,data.frame(x=xg)), col=ORANGE, lwd=1.8, lty=2)
lines(xg, predict(lo,xg), col=RED, lwd=1.6, lty=3); lines(sp, col=GREEN, lwd=1.6, lty=4)
legend("bottomright", c("mgcv GP smoother s(bs='gp')","mgcv penalised spline s()","loess","smooth.spline"),
col=c(BLUE,ORANGE,RED,GREEN), lwd=2, lty=1:4, bty="n", cex=.8)
cat(sprintf("mgcv GP effective df = %.1f; all four smoothers trace the same dip-and-rebound.\n", sum(g_gp$edf)))
cat("mgcv's bs='gp' is a Gaussian process and the smoothing spline is a GP posterior mean, so these are the\n")
cat("frequentist face of the from-scratch Bayesian GP (and PyMC gp.Marginal) in the Python notebook.\n")
mgcv GP effective df = 11.7; all four smoothers trace the same dip-and-rebound.
mgcv's bs='gp' is a Gaussian process and the smoothing spline is a GP posterior mean, so these are the
frequentist face of the from-scratch Bayesian GP (and PyMC gp.Marginal) in the Python notebook.
2. Mauna Loa CO₂ — trend and seasonal, decomposed by a GAM¶
The monthly CO₂ record is a rising trend with a yearly cycle. A mgcv GAM with a smooth trend term plus a cyclic seasonal term separates the two (the additive-model analogue of composing GP kernels) and predicts forward with confidence bands — a third graph zooms on the forecast region.
co<-read.csv("co2.csv"); co$month<-as.integer(format(as.Date(co$date),"%m"))
g<-gam(co2 ~ s(year, k=30) + s(month, bs="cc", k=12), data=co, knots=list(month=c(0.5,12.5)))
cat("GAM deviance explained:", round(summary(g)$dev.expl*100,2), "%\n")
xf<-data.frame(year=seq(1958,2012,length=600)); xf$month<-round((xf$year%%1)*12)+1
pf<-predict(g,newdata=xf,se.fit=TRUE)
par(mfrow=c(1,2), mar=c(4,4,3,1))
ylim_f <- range(c(co$co2, pf$fit - 1.96*pf$se.fit, pf$fit + 1.96*pf$se.fit))
plot(co$year,co$co2,pch=19,col=GREY,cex=.35,xlab="year",ylab="CO2 (ppm)",main="GAM fit + forecast",
xlim=c(1958,2012), ylim=ylim_f)
polygon(c(xf$year,rev(xf$year)), c(pf$fit-1.96*pf$se.fit, rev(pf$fit+1.96*pf$se.fit)), col=adjustcolor(BLUE,.25), border=NA)
lines(xf$year,pf$fit,col=BLUE,lwd=1.6); abline(v=max(co$year),lty=3)
mon<-1:12; se<-predict(g,newdata=data.frame(year=1985,month=mon),type="terms",terms="s(month)")
plot(mon, se, type="b", pch=19, col=GREEN, xlab="month", ylab="seasonal effect (ppm)", main="Seasonal cycle (cyclic smooth)")
par(mfrow=c(1,1))
# zoom on the forecast region
zi<-xf$year>=1996; di<-co$year>=1996
ylim_z <- range(c(co$co2[di], (pf$fit - 1.96*pf$se.fit)[zi], (pf$fit + 1.96*pf$se.fit)[zi]))
plot(co$year[di], co$co2[di], pch=19, col=GREY, cex=.5, xlim=c(1996,2012), ylim=ylim_z,
xlab="year", ylab="CO2 (ppm)",
main="Zoom on the forecast: annual cycle continues, band widens with lead time")
polygon(c(xf$year[zi],rev(xf$year[zi])), c((pf$fit-1.96*pf$se.fit)[zi], rev((pf$fit+1.96*pf$se.fit)[zi])), col=adjustcolor(BLUE,.25), border=NA)
lines(xf$year[zi], pf$fit[zi], col=BLUE, lwd=1.9); abline(v=max(co$year), lty=3); text(max(co$year)+0.2, min(co$co2[di]), "forecast", pos=4, cex=.8)
cat(sprintf("\nforecast 2012: %.1f ppm (actual ~394). Trend + a ~6 ppm yearly cycle, separated by the additive model;\n", pf$fit[which.min(abs(xf$year-2012))]))
i2012<-which.min(abs(xf$year-2012)); ACT<-393.8
lo<-pf$fit[i2012]-1.96*pf$se.fit[i2012]; hi<-pf$fit[i2012]+1.96*pf$se.fit[i2012]
cat("the band widens past the data. Whether it widens ENOUGH is a separate question, so we check:
")
cat(sprintf(" 95%% band at 2012 = [%.1f, %.1f]; actual %.1f -> inside? %s (error %+.1f ppm = %+.1f se)
",
lo, hi, ACT, ifelse(lo<=ACT & ACT<=hi,"YES","NO"), pf$fit[i2012]-ACT,
(pf$fit[i2012]-ACT)/pf$se.fit[i2012]))
cat("The GAM extrapolates better than the from-scratch GP in the Python notebook (which lands near 386 and
")
cat("misses by roughly 6 standard deviations), because a thin-plate spline on year keeps some curvature
")
cat("past the boundary while an RBF trend kernel reverts toward its mean and carries only the local slope.
")
cat("Neither band, however, prices the risk that the TREND MODEL itself is wrong -- which over a ten-year
")
cat("extrapolation of an accelerating series is the dominant risk. Both are honest about the function;
")
cat("neither is honest about the functional form, because no posterior band ever is.
")
GAM deviance explained: 99.97 %
forecast 2012: 391.4 ppm (actual ~394). Trend + a ~6 ppm yearly cycle, separated by the additive model;
the band widens past the data. Whether it widens ENOUGH is a separate question, so we check:
95% band at 2012 = [387.6, 395.2]; actual 393.8 -> inside? YES (error -2.4 ppm = -1.3 se)
The GAM extrapolates better than the from-scratch GP in the Python notebook (which lands near 386 and
misses by roughly 6 standard deviations), because a thin-plate spline on year keeps some curvature
past the boundary while an RBF trend kernel reverts toward its mean and carries only the local slope.
Neither band, however, 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.
3. Summary¶
The standard R smoothing package mgcv reproduces the Python result: its bs="gp" smoother is a Gaussian process, and it, the penalised spline, loess and smooth.spline all trace the same motorcycle curve — the cubic smoothing spline being a GP posterior mean, so the Bayesian, penalised and GP views coincide. A mgcv GAM split the Mauna Loa CO₂ series into a trend and a cyclic seasonal component (the additive echo of composing GP kernels) and forecast it with widening uncertainty.
These packages are the frequentist cross-check on the from-scratch Cholesky GP and PyMC gp.Marginal in gp_python.ipynb. The same Keeling curve is modelled a different way (SARIMA / structural state-space) in the time-series arc (Seasonal ARIMA (SARIMA)). Next the arc pushes the GP through a link function for classification and point-process intensity.