Gaussian Processes & Splines — the Bayesian Kernel View

Python (scikit-learn, pyGAM) · R (kernlab, mgcv)

The Question

The machine-learning arc scores frequentist predictors — trees, forests, boosting, penalised linear models — on two running datasets. This example asks the obvious follow-up: take the Bayesian nonparametric methods built elsewhere in the collection and turn them loose on the same Taiwan credit-default and California-housing problems. Not a toy comparison — the identical train/test splits, the same scoreboard. Two questions: do they compete on accuracy, and what do they give that a point predictor cannot?

Gaussian processes — a mean and a variance everywhere

A Gaussian process puts a prior over functions, $f\sim\mathcal{GP}(0,k)$ , and returns a posterior with a mean and a variance at every input. Fitting maximises the marginal likelihood for the length-scale and noise. On California it reaches RMSE 0.611 from a 1,500-point subsample — and the error bars are the point: 90% intervals cover 91.0% of outcomes, and mean absolute error rises monotonically with the predicted standard deviation. The uncertainty is calibrated and informative; where the GP says it is unsure, it really is more wrong. On credit it reaches AUC 0.742 from just 1,200 of 21,000 training rows.

\bar f(x_*) = k(x_*,X)\,[K+\sigma^2 I]^{-1}y, \qquad \operatorname{Var}[f_*] = k(x_*,x_*) - k(x_*,X)\,[K+\sigma^2 I]^{-1}k(X,x_*)

Additive splines — nonlinearity you can read

A GAM is a sum of smooth per-feature functions, each a penalised spline — and the roughness penalty is exactly a Bayesian random-walk prior on the basis coefficients, which is why the smooths come with credible bands rather than confidence heuristics. It is additive, so no interactions unless asked, but it scales to the full data cheaply and every effect stays a readable curve. California RMSE 0.643, credit AUC 0.765 — within 0.01 of the tree ensembles while remaining interpretable.

g\big(E[y]\big) = \beta_0 + \sum_j f_j(x_j), \qquad f_j \text{ a penalised spline}

The Scoreboard

The honest verdict is competitive, not dominant. On California all the Bayesian methods clear the linear model comfortably; on credit the GAM comes within 0.01 AUC of the ensembles and the GP within 0.03 on a fraction of the data. Gradient boosting still edges raw accuracy at full scale. What none of the frequentist winners carry is a distribution.

method	credit AUC (higher better)	California RMSE (lower better)	uncertainty?
logistic / linear	0.715	0.737	no
random forest	0.775	0.523	no
XGBoost	0.774	0.494	no
Gaussian process (1,500-pt subsample)	0.742	0.611	yes — calibrated std
GAM (additive splines, full data)	0.765	0.643	yes — credible bands
BART (`pymc-bart`; 0.526 with `dbarts`)	0.760	0.665	yes — posterior intervals

BART, stated precisely

BART is the split verdict, and worth stating precisely rather than smoothing over. It is level with the forest on credit (0.760 against 0.775) but far off it on California (0.665 against 0.523) — and that gap is an implementation limit, not a limit of the method. The BART example diagnoses it: pymc-bart's in-sample RMSE is 0.641, meaning it underfits, while R's dbarts reaches 0.526 on the identical task. Carrying the R figure into the R scoreboard puts BART back alongside the forest.

Why the two languages are not pooled

A methodological point that the R companion makes explicit. Python and R use different 70/30 splits, so their numbers are not interchangeable — a scoreboard that mixes them is comparing across resamples, not across methods. Each notebook therefore refits its own logistic and linear baselines and carries ensemble figures only from notebooks on its own split: the R scoreboard uses R's forest (0.774 credit, 0.491 California) and R's XGBoost (0.776, 0.484), never Python's.

The R companion also shows where a modelling choice, not the language, moves the number. Its mgcv GAM reaches credit AUC 0.743 against pyGAM's 0.765, because the smooths are fitted only to the continuous features while the discrete PAY_* and demographic terms enter linearly; pyGAM splines every column. On California, where all features are continuous, mgcv comes out ahead at 0.607.

The cost

The standing caveat is scale. Exact GPs cost $O(n^3)$ , which is why every GP here is fit on a subsample — and that limit is itself the finding. It is precisely why gradient boosting, not the Gaussian process, is the industrial default for large tabular data. Sparse and inducing-point approximations close part of the gap.

The Bayesian Payoff

Accuracy is close to a tie; the difference is what the models report. A booster returns one number per case. Every Bayesian method here returns a distribution — GP predictive standard deviation, GAM credible band, BART posterior interval — and that distribution is calibrated (90% coverage measured at 0.910) and informative (error grows with the reported standard deviation). For credit decisions and portfolio construction that is not decoration: it is the difference between a point score and a risk-aware decision — size the position by the confidence, route the uncertain applications to manual review.

Where this sits

This example is the hinge between two sections. The methods come from Bayesian Nonparametrics — Gaussian-Process Regression develops the GP on the motorcycle and CO₂ series with full MCMC, GP Classification & Log-Gaussian Cox Processes the classification case, Bayesian Penalised Splines & Additive Models the spline construction from the basis up. The benchmark is the ML arc's. And the GP's posterior mean is kernel ridge regression: Support Vector Machines & Kernel Methods verifies that identity to machine precision, and resolves the apparent RMSE gap between the two — cross-validate kernel ridge's hyperparameters and it lands on this example's 0.611.

Notebooks

Python GP regression and classification with calibrated predictive intervals, pyGAM additive smooths with credible bands, the full scoreboard against the frequentist field, and a side-by-side of point prediction versus posterior interval (Python) Open → R Independent R implementation via kernlab::gausspr and mgcv — the same two tasks, R-computed baselines throughout, and the additive-structure choice that explains the credit-AUC difference from pyGAM (R) Open →

Downloads

cali_housing.csv California housing — the regression benchmark, shared across the whole machine-learning arc credit_default.csv Default of credit-card clients (UCI) — the classification benchmark, shared across the whole machine-learning arc

References

Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. — the GP, its marginal likelihood, and the kernel-ridge equivalence
Wahba, G. (1978). Improper priors, spline smoothing and the problem of guarding against model errors in regression. JRSS-B 40(3), 364–372. — the roughness penalty as a Bayesian prior
Wood, S. N. (2017). Generalized Additive Models: An Introduction with R, 2nd ed. CRC Press. — mgcv, REML smoothness selection and the credible bands used here
Hastie, T. & Tibshirani, R. (1990). Generalized Additive Models. Chapman & Hall. — the additive structure and its interpretability
Chipman, H., George, E. & McCulloch, R. (2010). BART: Bayesian additive regression trees. Annals of Applied Statistics 4(1), 266–298. — the Bayesian tree ensemble on the scoreboard
Nickisch, H. & Rasmussen, C. E. (2008). Approximations for binary Gaussian process classification. JMLR 9, 2035–2078. — the Laplace approximation behind the GP classifier