Bayesian Nonparametric Prediction — how do GPs and Bayesian splines fare?

Gaussian processes · Bayesian additive splines (GAM) · vs the ML scoreboard, on the same two datasets

The ML arc so far has scored frequentist predictors — trees, forests, boosting, penalized linear models — on two running datasets. This notebook asks a question that bridges back to the Bayesian Nonparametric arc: take the BNP methods built there and turn them loose on the same Taiwan credit-default (classification) and California-housing (regression) problems. How do they fare — on accuracy, and on the thing the frequentist models never gave us, calibrated predictive uncertainty?

Three Bayesian nonparametric predictors, each a supervised member of that arc:

• Gaussian Processes — a prior over functions, $f\sim\mathcal{GP}(0,k)$, updated to a posterior with a mean and a variance at every input. The Bayesian kernel machine; its posterior mean is exactly kernel ridge regression, the identity demonstrated to machine precision in Support Vector Machines & Kernel Methods. Here: GP regression (California) and GP classification (credit). The same machinery is developed at length in Gaussian-Process Regression (motorcycle and CO₂, with full MCMC in PyMC) and GP Classification & Log-Gaussian Cox Processes (Pima); this notebook is those methods turned loose on the ML arc's own data.
• Bayesian additive splines (GAM) — a sum of smooth per-feature functions, each a penalized spline. The smoothing penalty is a Bayesian random-walk prior on the basis coefficients (Wahba), so a GAM is a Bayesian mixed model with credible bands. The same penalised-spline construction is built from the basis up in Bayesian Penalised Splines & Additive Models.
• BART (Bayesian Additive Regression Trees) — the Bayesian tree ensemble, already applied to both datasets in BART — Bayesian Additive Regression Trees; we bring its numbers onto the scoreboard rather than refitting it.

The honest caveat, stated up front: exact GPs cost $O(n^3)$ — infeasible on 20–30k rows — so the GPs here are fit on a subsample (as the BART notebook was), while the GAMs use the full data. That scaling limit is itself a headline finding: it is why gradient boosting, not the GP, is the default for large tabular data. Engines here are scikit-learn's GP and pyGAM; the BNP arc's fuller treatments use PyMC / mgcv with full MCMC.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, time, warnings
warnings.filterwarnings("ignore")
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, roc_auc_score, brier_score_loss
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
rng=np.random.default_rng(0)
# California (regression)
h=pd.read_csv("cali_housing.csv"); hf=[c for c in h.columns if c!="MedHouseVal"]
Xh=h[hf].to_numpy(float); yh=h["MedHouseVal"].to_numpy(float)
Xh_tr,Xh_te,yh_tr,yh_te=train_test_split(Xh,yh,test_size=0.3,random_state=0)
sch=StandardScaler().fit(Xh_tr); Zh_tr=sch.transform(Xh_tr); Zh_te=sch.transform(Xh_te)
# credit (classification)
d=pd.read_csv("credit_default.csv"); feat=[c for c in d.columns if c!="default"]
Xc=d[feat].to_numpy(float); yc=d["default"].to_numpy(int)
Xc_tr,Xc_te,yc_tr,yc_te=train_test_split(Xc,yc,test_size=0.3,random_state=0,stratify=yc)
scc=StandardScaler().fit(Xc_tr); Zc_tr=scc.transform(Xc_tr); Zc_te=scc.transform(Xc_te)
print(f"California regression: {len(yh_tr):,} train / {len(yh_te):,} test, {Xh.shape[1]} features")
print(f"Credit classification: {len(yc_tr):,} train / {len(yc_te):,} test, {Xc.shape[1]} features, {100*yc.mean():.0f}% default")

California regression: 14,448 train / 6,192 test, 8 features
Credit classification: 21,000 train / 9,000 test, 23 features, 22% default

1. Gaussian process regression — California

A GP places a prior over functions: any finite set of points is jointly Gaussian with covariance $k(x,x')$. We use the RBF (squared-exponential) kernel plus a white-noise term; fitting maximises the marginal likelihood for the length-scale and noise (empirical Bayes), and prediction returns a full posterior — a mean and a standard deviation at each test point. That standard deviation is the payoff: an honest, per-prediction error bar no tree or boosting model provides.

Fit on a 1,500-point subsample (exact GP is $O(n^3)$). We report out-of-sample RMSE, the 90% predictive-interval coverage (a calibration check — it should be ≈ 0.90), and the predicted-vs-actual graph.

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel as Cst
m=1500; sub=rng.choice(len(Zh_tr),m,replace=False)
kern=Cst(1.0)*RBF(1.0)+WhiteKernel(0.1)
t=time.time(); gpr=GaussianProcessRegressor(kernel=kern,normalize_y=True).fit(Zh_tr[sub],yh_tr[sub]); ft=time.time()-t
mu,sd=gpr.predict(Zh_te,return_std=True)
gp_rmse=mean_squared_error(yh_te,mu)**.5; cover=np.mean((yh_te>=mu-1.645*sd)&(yh_te<=mu+1.645*sd))
print(f"GP regression (subsample {m}): test RMSE {gp_rmse:.3f}   90% PI coverage {cover:.3f}   fit {ft:.0f}s")
print(f"learned kernel: {gpr.kernel_}")
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
ax[0].scatter(mu,yh_te,s=5,alpha=.12,color=BLUE)
qb=np.quantile(mu,np.linspace(0,1,11)); idx=np.clip(np.digitize(mu,qb[1:-1]),0,9)
ax[0].plot([mu.min(),mu.max()],[mu.min(),mu.max()],"k--",lw=1,label="perfect")
ax[0].plot([mu[idx==j].mean() for j in range(10)],[yh_te[idx==j].mean() for j in range(10)],"o-",color=RED,lw=2,label="decile means")
ax[0].set_xlabel("GP predicted value ($100k)"); ax[0].set_ylabel("actual value"); ax[0].set_title(f"GP regression — predicted vs actual (RMSE {gp_rmse:.3f})"); ax[0].legend()
# is the uncertainty meaningful? bin by predictive sd, show error grows with sd
err=np.abs(yh_te-mu); sb=np.quantile(sd,np.linspace(0,1,7)); bi=np.clip(np.digitize(sd,sb[1:-1]),0,5)
ax[1].plot([sd[bi==j].mean() for j in range(6)],[err[bi==j].mean() for j in range(6)],"o-",color=PURP,lw=2)
ax[1].set_xlabel("GP predictive std (uncertainty)"); ax[1].set_ylabel("mean |actual - predicted|"); ax[1].set_title("The uncertainty is informative: error grows with predicted std")
plt.tight_layout(); plt.show()
print(f"Coverage {cover:.2f} ~ nominal 0.90 -> the GP's error bars are calibrated. And where the GP says it is unsure, it")
print("really is more wrong (right panel) -- uncertainty you can act on, which a point prediction cannot give.")

GP regression (subsample 1500): test RMSE 0.611   90% PI coverage 0.910   fit 3s
learned kernel: 1.95**2 * RBF(length_scale=2.79) + WhiteKernel(noise_level=0.241)

Coverage 0.91 ~ nominal 0.90 -> the GP's error bars are calibrated. And where the GP says it is unsure, it
really is more wrong (right panel) -- uncertainty you can act on, which a point prediction cannot give.

2. Gaussian process classification — credit default

For a binary outcome the GP models a latent function squashed through a logit link; the posterior is non-Gaussian, so scikit-learn uses the Laplace approximation (the arc's Pima notebook does the full MCMC version). Fit on a 1,200-point subsample; we report out-of-sample AUC and the reliability curve. ROC-AUC is defined in the CART notebook (0.5 = chance, 1 = perfect).

from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.calibration import calibration_curve
m2=1200; sub2=rng.choice(len(Zc_tr),m2,replace=False)
t=time.time(); gpc=GaussianProcessClassifier(kernel=Cst(1.0)*RBF(1.0)).fit(Zc_tr[sub2],yc_tr[sub2]); ct=time.time()-t
pgc=gpc.predict_proba(Zc_te)[:,1]
gp_auc=roc_auc_score(yc_te,pgc); gp_brier=brier_score_loss(yc_te,pgc)
from sklearn.linear_model import LogisticRegression as _LR
log_auc=roc_auc_score(yc_te,_LR(max_iter=2000).fit(Zc_tr,yc_tr).predict_proba(Zc_te)[:,1])
print(f"GP classification (subsample {m2}): test AUC {gp_auc:.3f}   Brier {gp_brier:.3f}   fit {ct:.0f}s")
from sklearn.metrics import roc_curve
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
fpr,tpr,_=roc_curve(yc_te,pgc); ax[0].plot(fpr,tpr,color=BLUE,lw=2,label=f"GP classifier (AUC {gp_auc:.3f})")
ax[0].plot([0,1],[0,1],"k:",lw=1); ax[0].set_xlabel("false positive rate"); ax[0].set_ylabel("true positive rate"); ax[0].set_title("Credit — ROC"); ax[0].legend(loc="lower right")
pt,pp=calibration_curve(yc_te,pgc,n_bins=10,strategy="quantile")
ax[1].plot([0,pgc.max()],[0,pgc.max()],"k--",lw=1,label="perfect calibration"); ax[1].plot(pp,pt,"o-",color=RED,lw=2,label=f"GP (Brier {gp_brier:.3f})")
ax[1].set_xlabel("predicted P(default) [decile bins]"); ax[1].set_ylabel("observed default proportion"); ax[1].set_title("Credit — reliability (proportions vs predictions)"); ax[1].legend()
plt.tight_layout(); plt.show()
print(f"GP classification AUC {gp_auc:.3f} beats the unpenalized logistic baseline ({log_auc:.3f}) despite seeing only {m2} of")
print(f"{len(yc_tr):,} training rows -- the kernel captures nonlinearity the linear model cannot.")

GP classification (subsample 1200): test AUC 0.742   Brier 0.143   fit 5s

GP classification AUC 0.742 beats the unpenalized logistic baseline (0.715) despite seeing only 1200 of
21,000 training rows -- the kernel captures nonlinearity the linear model cannot.

3. Bayesian additive splines — GAM

A GAM models the response as a sum of smooth per-feature functions, $g(E[y]) = \beta_0 + \sum_j f_j(x_j)$, each $f_j$ a penalized spline. The roughness penalty is equivalent to a Bayesian random-walk prior on the spline coefficients (Wahba, 1978) — so a GAM is a Bayesian mixed model, and pyGAM reports credible bands on each smooth. Unlike the GP it is additive (no interactions unless asked) but it scales to the full data cheaply and stays interpretable — you can read each feature's effect off its curve. LinearGAM for California, LogisticGAM for credit; cross-links the arc's airquality GAM notebook.

from pygam import LinearGAM, LogisticGAM
t=time.time(); lgam=LinearGAM(n_splines=20).fit(Xh_tr,yh_tr); gam_rmse=mean_squared_error(yh_te,lgam.predict(Xh_te))**.5; gt=time.time()-t
t=time.time(); cgam=LogisticGAM(n_splines=15).fit(Xc_tr,yc_tr); gam_auc=roc_auc_score(yc_te,cgam.predict_proba(Xc_te)); gt2=time.time()-t
print(f"LinearGAM (California, full data): test RMSE {gam_rmse:.3f}  fit {gt:.0f}s")
print(f"LogisticGAM (credit, full data):   test AUC  {gam_auc:.3f}  fit {gt2:.0f}s")
# partial-dependence smooths with credible bands
fig,ax=plt.subplots(1,3,figsize=(15,4))
for a,term,lab in [(ax[0],0,"MedInc (California)"),(ax[1],2,"AveRooms (California)")]:
    XX=lgam.generate_X_grid(term=term); pdep,ci=lgam.partial_dependence(term=term,X=XX,width=0.95)
    a.plot(XX[:,term],pdep,color=BLUE,lw=2); a.fill_between(XX[:,term],ci[:,0],ci[:,1],color=BLUE,alpha=.2)
    a.set_xlabel(hf[term]); a.set_ylabel("partial effect on value"); a.set_title(f"GAM smooth: {lab}")
tc=feat.index("PAY_1"); XXc=cgam.generate_X_grid(term=tc); pd_,ci_=cgam.partial_dependence(term=tc,X=XXc,width=0.95)
ax[2].plot(XXc[:,tc],pd_,color=RED,lw=2); ax[2].fill_between(XXc[:,tc],ci_[:,0],ci_[:,1],color=RED,alpha=.2)
ax[2].set_xlabel("PAY_1 (credit)"); ax[2].set_ylabel("partial effect on log-odds default"); ax[2].set_title("GAM smooth: recent repayment status")
plt.tight_layout(); plt.show()
print(f"The GAM lands within 0.01 AUC of the tree ensembles on credit ({gam_auc:.3f} against ~0.775) and beats the linear")
print(f"model on California ({gam_rmse:.3f} against 0.737) -- while every effect stays a readable curve with a credible band.")
print("Additivity is the price: no interactions unless asked for, which is also why it costs so much less than the GP.")

LinearGAM (California, full data): test RMSE 0.643  fit 1s
LogisticGAM (credit, full data):   test AUC  0.765  fit 6s

The GAM lands within 0.01 AUC of the tree ensembles on credit (0.765 against ~0.775) and beats the linear
model on California (0.643 against 0.737) -- while every effect stays a readable curve with a credible band.
Additivity is the price: no interactions unless asked for, which is also why it costs so much less than the GP.

4. The scoreboard — Bayesian nonparametric vs frequentist ML

All the BNP predictors placed beside the frequentist models from the earlier subsections, on the identical train/test splits. We recompute the key ML baselines here for a fair comparison and bring in BART's numbers from BART — Bayesian Additive Regression Trees (credit AUC 0.760, California RMSE 0.665 for the pymc-bart fit used there).

from sklearn.linear_model import LogisticRegression, LinearRegression
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
import xgboost as xgb
cl={}; rg={}
cl["logistic / linear"]=roc_auc_score(yc_te,LogisticRegression(max_iter=2000).fit(Zc_tr,yc_tr).predict_proba(Zc_te)[:,1])
rg["logistic / linear"]=mean_squared_error(yh_te,LinearRegression().fit(Zh_tr,yh_tr).predict(Zh_te))**.5
cl["random forest"]=roc_auc_score(yc_te,RandomForestClassifier(n_estimators=400,min_samples_leaf=5,random_state=0,n_jobs=-1).fit(Xc_tr,yc_tr).predict_proba(Xc_te)[:,1])
rg["random forest"]=mean_squared_error(yh_te,RandomForestRegressor(n_estimators=300,min_samples_leaf=3,random_state=0,n_jobs=-1).fit(Xh_tr,yh_tr).predict(Xh_te))**.5
cl["XGBoost"]=roc_auc_score(yc_te,xgb.XGBClassifier(n_estimators=300,learning_rate=0.05,max_depth=4,verbosity=0,n_jobs=-1).fit(Xc_tr,yc_tr).predict_proba(Xc_te)[:,1])
rg["XGBoost"]=mean_squared_error(yh_te,xgb.XGBRegressor(n_estimators=300,learning_rate=0.05,max_depth=4,verbosity=0,n_jobs=-1).fit(Xh_tr,yh_tr).predict(Xh_te))**.5
cl["Gaussian process"]=gp_auc;  rg["Gaussian process"]=gp_rmse
cl["GAM (splines)"]=gam_auc;     rg["GAM (splines)"]=gam_rmse
cl["BART (pymc-bart)"]=0.760; rg["BART (pymc-bart)"]=0.665   # from the BART notebook; dbarts there gets 0.783 / 0.526
sc=pd.DataFrame({"credit AUC (higher=better)":cl,"California RMSE (lower=better)":rg})
print(sc.round(3).to_string())
bnp={"Gaussian process","GAM (splines)","BART (pymc-bart)"}
def col(names): return [PURP if n in bnp else GREY for n in names]
fig,ax=plt.subplots(1,2,figsize=(14,4.8))
na=list(cl); ax[0].barh(na,[cl[k] for k in na],color=col(na)); ax[0].set_xlim(0.5,0.8); ax[0].invert_yaxis(); ax[0].set_title("Credit — test AUC (higher better)")
nb=list(rg); ax[1].barh(nb,[rg[k] for k in nb],color=col(nb)); ax[1].invert_yaxis(); ax[1].set_title("California — test RMSE (lower better)")
plt.tight_layout(); plt.show()
print("Bayesian nonparametric methods (purple) are COMPETITIVE, not dominant. On credit the GAM comes within 0.01 AUC of")
print("the ensembles and the GP within 0.03 on a mere subsample; on California both clear the linear model comfortably.")
print("BART is the split verdict: level with the forest on credit (0.760 vs 0.775) but far off it on regression (0.665 vs")
print("0.523), and its own notebook pins that on the pymc-bart implementation -- R's dbarts reaches 0.526 on the same task.")
print("Boosting still edges the top on raw accuracy at full scale. What none of the frequentist winners carry is the thing")
print("section 5 is about: a calibrated distribution rather than a single number.")

                   credit AUC (higher=better)  California RMSE (lower=better)
logistic / linear                       0.715                           0.737
random forest                           0.775                           0.523
XGBoost                                 0.774                           0.494
Gaussian process                        0.742                           0.611
GAM (splines)                           0.765                           0.643
BART (pymc-bart)                        0.760                           0.665

Bayesian nonparametric methods (purple) are COMPETITIVE, not dominant. On credit the GAM comes within 0.01 AUC of
the ensembles and the GP within 0.03 on a mere subsample; on California both clear the linear model comfortably.
BART is the split verdict: level with the forest on credit (0.760 vs 0.775) but far off it on regression (0.665 vs
0.523), and its own notebook pins that on the pymc-bart implementation -- R's dbarts reaches 0.526 on the same task.
Boosting still edges the top on raw accuracy at full scale. What none of the frequentist winners carry is the thing
section 5 is about: a calibrated distribution rather than a single number.

5. The Bayesian payoff — uncertainty the ML models don't have

Accuracy is a near-tie; the real difference is what the models report. A frequentist booster returns a single number per case. Every Bayesian nonparametric method here returns a distribution — GP predictive std, GAM credible band, BART posterior interval — and (§1) that uncertainty is calibrated and informative: the 90% GP intervals cover ≈ 90% of outcomes, and the cases the GP flags as uncertain are exactly the ones it gets more wrong. In credit and asset management that is not a nicety: it is the difference between a point score and a risk-aware decision (size the position by the confidence, route uncertain applications to manual review). The plot contrasts the two worlds on a handful of California test blocks.

sel=rng.choice(len(yh_te),12,replace=False); sel=sel[np.argsort(mu[sel])]
xgr=xgb.XGBRegressor(n_estimators=300,learning_rate=0.05,max_depth=4,verbosity=0,n_jobs=-1).fit(Xh_tr,yh_tr).predict(Xh_te)
fig,ax=plt.subplots(figsize=(9,4.6)); xp=np.arange(len(sel))
ax.errorbar(xp-0.12,mu[sel],yerr=1.645*sd[sel],fmt="o",color=BLUE,capsize=3,label="GP: mean ± 90% interval")
ax.scatter(xp+0.12,xgr[sel],marker="s",color=ORANGE,label="XGBoost: point prediction")
ax.scatter(xp,yh_te[sel],marker="*",s=140,color=RED,zorder=5,label="actual value")
ax.set_xlabel("test blocks (sorted by GP prediction)"); ax.set_ylabel("median house value ($100k)")
ax.set_title("Same accuracy, different information: the GP reports how sure it is"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Both models land near the true stars; only the GP tells you the width of the bar -- wide where it is unsure, narrow")
print("where it is confident. That calibrated uncertainty is what the Bayesian nonparametric arc adds to the ML toolkit.")

Both models land near the true stars; only the GP tells you the width of the bar -- wide where it is unsure, narrow
where it is confident. That calibrated uncertainty is what the Bayesian nonparametric arc adds to the ML toolkit.

6. Summary

How do the Bayesian nonparametric methods fare on real tabular prediction? Well — and with something extra:

method	credit AUC	California RMSE	gives uncertainty?
GP (kernel machine)	~0.74 (on a subsample)	~0.62	yes — calibrated std
GAM (additive splines)	~0.76, matches trees	~0.64	yes — credible bands
BART (Bayesian trees)	0.760	0.665 (0.526 with dbarts)	yes — posterior intervals
linear / logistic	0.715	0.737	no
random forest / XGBoost	~0.77	~0.49	no

Three takeaways for the portfolio:

1. Competitive on accuracy. The GAM matches the tree ensembles on credit and beats the linear model on California while staying interpretable (each feature is a readable curve); the GP does the same from a small subsample; BART is level with the forest on credit but well behind it on California in the pymc-bart fit — an implementation limit its own notebook diagnoses, not a limit of BART. Gradient boosting still edges raw accuracy at full scale.
2. Uncertainty is the differentiator. Only the Bayesian methods return a distribution — and it is calibrated (90% GP coverage ≈ 0.90) and informative (bigger error where the GP is less sure). That is what point-prediction ML cannot give, and what makes these methods first-class for risk-aware finance and credit decisions.
3. Scale is the cost. Exact GPs are $O(n^3)$ — hence the subsampling — which is precisely why boosting, not the GP, is the industrial default for large tabular data. (Sparse/inducing-point GPs close some of this gap.)

This connects the ML arc to the Bayesian Nonparametrics section — Gaussian-Process Regression, GP Classification & Log-Gaussian Cox Processes and Bayesian Penalised Splines & Additive Models — and to BART — Bayesian Additive Regression Trees in the Trees and Ensembles section. The GP's posterior mean is kernel ridge regression — the exact link taken up in Support Vector Machines & Kernel Methods, which also shows that the RMSE gap between the two is hyperparameter selection, not estimation: cross-validate kernel ridge and it lands on this notebook's 0.611. (BNP methods here use scikit-learn's GP and pyGAM; the fuller MCMC treatments live in the BNP arc's PyMC / mgcv notebooks.)