GP Classification & Log-Gaussian Cox Processes¶
The latent Gaussian process, pushed through a link¶
Gaussian-process regression modelled a real-valued function. Send that latent function through a link and the same prior over functions models very different data:
- Classification — a binary or count outcome through a logit link: $$f\sim\mathcal{GP}(0,k),\qquad y_i\mid f\sim\text{Bernoulli}\big(\sigma(f(x_i))\big).$$ The GP becomes a flexible, nonlinear log-odds surface with calibrated uncertainty.
- Log-Gaussian Cox process — a point pattern whose intensity is a GP: $$\log\lambda(t)=f(t),\qquad \text{count in a bin}\sim\text{Poisson}(\lambda\cdot\text{width}).$$ The smooth, positive rate of events over time, with uncertainty.
Now the likelihood is non-Gaussian, so the posterior over $f$ has no closed form. The classic from-scratch fix is the Laplace approximation (Rasmussen & Williams, Alg. 3.1): Newton-iterate to the posterior mode of $f$ and approximate the posterior by the Gaussian there; the Laplace marginal likelihood then picks the kernel hyperparameters. We validate on a known non-monotone probability, fit the Pima diabetes data (a 2-D probability surface, against logistic regression), model the coal-mining disasters as a Cox process (against a kernel rate estimate), and cross-check the intensity in PyMC with gp.Latent.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.special import expit
import gplink as GL
rng = np.random.default_rng(3)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("A latent GP through a link: logit -> classification, log -> Cox-process intensity.")
A latent GP through a link: logit -> classification, log -> Cox-process intensity.
1. Does it work? — recovering a non-monotone probability¶
A binary outcome whose success probability rises and falls with $x$: $p(x)=\sigma(1.5\sin(1.2x))$. Logistic regression can only bend one way (its log-odds is linear in $x$), so it cannot capture the wave — a GP classifier can. We fit the latent GP by Laplace and squash the predictive latent through the link.
x=np.sort(rng.uniform(-3.2,3.2,220)); ptrue=lambda t: expit(1.5*np.sin(1.2*t))
y=(rng.random(len(x))<ptrue(x)).astype(float)
ker=lambda th,A,B: GL.k_rbf(A,B,np.exp(th[0]),np.exp(th[1]))
xg=np.linspace(-3.2,3.2,200)
res=GL.fit_predict(x,y,xg,GL.lik_bernoulli,None,ker,[0.,0.])
p_gp=GL.predictive_prob(res["mean"],res["var"]); sd=np.sqrt(res["var"])
p_lo=GL.predictive_prob(res["mean"]-1.96*sd, res["var"]*0); p_hi=GL.predictive_prob(res["mean"]+1.96*sd, res["var"]*0)
from sklearn.linear_model import LogisticRegression
lr=LogisticRegression().fit(x[:,None],y); p_lr=lr.predict_proba(xg[:,None])[:,1]
print(f"GP recovers p(x): RMSE {np.sqrt(np.mean((p_gp-ptrue(xg))**2)):.3f} vs logistic RMSE {np.sqrt(np.mean((p_lr-ptrue(xg))**2)):.3f}")
fig,ax=plt.subplots(figsize=(9,4.3))
ax.fill_between(xg,p_lo,p_hi,color=BLUE,alpha=.18,label="GP 95% band")
ax.plot(xg,ptrue(xg),color="k",lw=1.6,ls="--",label="true p(x)")
ax.plot(xg,p_gp,color=BLUE,lw=2.2,label="GP classifier")
ax.plot(xg,p_lr,color=RED,lw=1.8,ls=":",label="logistic regression")
ax.scatter(x,y*1.03-0.015,s=8,color=GREY,alpha=.4)
ax.set_xlabel("x"); ax.set_ylabel("P(y=1)"); ax.set_title("GP classifier recovers the wave; logistic cannot bend twice")
ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print("The GP tracks the rise-and-fall of the true probability with an honest band; logistic regression, whose")
print("log-odds is linear in x, is forced to a single monotone S -- it fits a nearly flat line and misses the wave.")
GP recovers p(x): RMSE 0.062 vs logistic RMSE 0.212
The GP tracks the rise-and-fall of the true probability with an honest band; logistic regression, whose log-odds is linear in x, is forced to a single monotone S -- it fits a nearly flat line and misses the wave.
2. The Pima data — a 2-D probability surface¶
532 Pima women, diabetes vs plasma glucose and BMI. A GP classifier learns a smooth 2-D log-odds surface; we map the predicted probability and the 0.5 decision boundary, and compare with logistic regression. On these two features the classes are close to linearly separable, so the accuracies nearly match — but the GP's real gift is the uncertainty: its probabilities are least confident where data is sparse, which the logistic surface never admits.
d=pd.read_csv("pima.csv"); X=d[["glu","bmi"]].to_numpy(); yb=d["diabetes"].to_numpy().astype(float)
mu,sdv=X.mean(0),X.std(0); Xs=(X-mu)/sdv
kg=lambda th,A,B: GL.k_rbf(A,B,np.exp(th[0]),np.exp(th[1]))
g1,g2=np.meshgrid(np.linspace(X[:,0].min(),X[:,0].max(),70), np.linspace(X[:,1].min(),X[:,1].max(),70))
grid=np.column_stack([g1.ravel(),g2.ravel()]); grids=(grid-mu)/sdv
rp=GL.fit_predict(Xs,yb,grids,GL.lik_bernoulli,None,kg,[0.5,0.])
P=GL.predictive_prob(rp["mean"],rp["var"]).reshape(g1.shape); U=np.sqrt(rp["var"]).reshape(g1.shape)
# in-sample accuracy for both
rp_in=GL.fit_predict(Xs,yb,Xs,GL.lik_bernoulli,None,kg,[0.5,0.]); pg=GL.predictive_prob(rp_in["mean"],rp_in["var"])
lr2=LogisticRegression().fit(Xs,yb); pl=lr2.predict_proba(Xs)[:,1]
print(f"accuracy GP {(( pg>.5)==yb).mean():.3f} logistic {((pl>.5)==yb).mean():.3f} (GP length-scale {np.exp(rp['theta'][0]):.1f} sd = nearly linear here)")
fig,ax=plt.subplots(1,2,figsize=(13,4.8))
c0=ax[0].contourf(g1,g2,P,levels=12,cmap="RdBu_r",vmin=0,vmax=1); ax[0].contour(g1,g2,P,levels=[0.5],colors="k",linewidths=2)
ax[0].scatter(X[yb==1,0],X[yb==1,1],s=10,color="darkred",alpha=.5,label="diabetes"); ax[0].scatter(X[yb==0,0],X[yb==0,1],s=10,color="navy",alpha=.4,label="none")
ax[0].set_xlabel("glucose"); ax[0].set_ylabel("BMI"); ax[0].set_title("GP P(diabetes) surface + 0.5 boundary"); ax[0].legend(frameon=False,fontsize=8); plt.colorbar(c0,ax=ax[0],fraction=.046)
c1=ax[1].contourf(g1,g2,U,levels=12,cmap="viridis"); ax[1].scatter(X[:,0],X[:,1],s=6,color="white",alpha=.35)
ax[1].set_xlabel("glucose"); ax[1].set_ylabel("BMI"); ax[1].set_title("GP predictive uncertainty (latent sd)"); plt.colorbar(c1,ax=ax[1],fraction=.046)
plt.tight_layout(); plt.show()
print("The boundary is smooth and, here, close to the logistic line -- Pima is nearly linearly separable in these")
print("two features, so the accuracies match. The second panel is what the GP adds: uncertainty that grows in the")
print("sparse corners of the feature space, where any single boundary is a guess. Logistic regression stays equally")
print("confident everywhere, including where it has seen almost no one.")
accuracy GP 0.776 logistic 0.773 (GP length-scale 5.5 sd = nearly linear here)
The boundary is smooth and, here, close to the logistic line -- Pima is nearly linearly separable in these two features, so the accuracies match. The second panel is what the GP adds: uncertainty that grows in the sparse corners of the feature space, where any single boundary is a guess. Logistic regression stays equally confident everywhere, including where it has seen almost no one.
3. Coal-mining disasters as a Cox process¶
191 fatal British coal-mining accidents, 1851–1962. Instead of a fixed rate we let the intensity be a GP: bin the events, model the counts as Poisson with $\lambda(t)=e^{f(t)}$, and infer the smooth rate. The famous feature — a sharp fall around 1890 (after mine-safety legislation) — appears as a drop in the intensity, now with a credible band. The frequentist counterpart is a kernel intensity estimate (a smoothed event histogram), which gives the same curve without the uncertainty.
This same disaster series is modelled a different way in the variable-selection arc (Reversible-Jump MCMC — Inference over Model Dimension), where a change-point process infers how many discrete jumps the rate undergoes. The two views are complementary: the Cox process assumes the rate changes smoothly, the change-point model assumes it is piecewise-constant with abrupt breaks — and both locate the late-Victorian fall near 1890.
ev=pd.read_csv("coal.csv",header=None)[0].to_numpy()
cx,cy,w=GL.bin_counts(ev,1851,1963,56)
kc=lambda th,A,B: GL.k_rbf(A,B,np.exp(th[0]),np.exp(th[1]))
cs=np.linspace(1851,1962,300)
rc=GL.fit_predict(cx,cy,cs,GL.lik_poisson,w,kc,[np.log(15),0.])
lam=np.exp(rc["mean"]+0.5*rc["var"]); sdf=np.sqrt(rc["var"])
lam_lo=np.exp(rc["mean"]-1.96*sdf); lam_hi=np.exp(rc["mean"]+1.96*sdf)
# frequentist kernel intensity: Gaussian-smoothed event rate
bw=6.0; kde=np.array([np.exp(-0.5*((t-ev)/bw)**2).sum()/(bw*np.sqrt(2*np.pi)) for t in cs])
fig,ax=plt.subplots(figsize=(9.5,4.3))
ax.bar(cx, cy/w, width=w*0.9, color=GREY, alpha=.35, label="binned rate (events/yr)")
ax.fill_between(cs, lam_lo, lam_hi, color=BLUE, alpha=.2, label="LGCP 95% band")
ax.plot(cs, lam, color=BLUE, lw=2.2, label="LGCP intensity (GP)")
ax.plot(cs, kde, color=RED, lw=1.8, ls="--", label="kernel intensity (frequentist)")
ax.set_xlabel("year"); ax.set_ylabel("disasters per year"); ax.set_title("Log-Gaussian Cox process: coal-mining disaster intensity")
ax.legend(frameon=False,fontsize=8); plt.tight_layout(); plt.show()
print(f"intensity: {lam[np.argmin(abs(cs-1860))]:.1f}/yr around 1860, {lam[np.argmin(abs(cs-1950))]:.2f}/yr by 1950 --")
print("a roughly four-fold decline, sharpest near 1890. The LGCP and the kernel estimate trace the same fall; the")
print("Cox process adds a credible band (wide early, where events are dense but the rate is high) around it.")
intensity: 3.0/yr around 1860, 0.71/yr by 1950 -- a roughly four-fold decline, sharpest near 1890. The LGCP and the kernel estimate trace the same fall; the Cox process adds a credible band (wide early, where events are dense but the rate is high) around it.
4. Cross-check in PyMC — gp.Latent for the Cox process¶
For a non-Gaussian likelihood the latent function must be sampled, which pm.gp.Latent does (non-centred, so it runs on Windows without scan). We put a GP on the log-intensity of the binned coal counts and let NUTS sample it, then compare the posterior intensity to the from-scratch Laplace fit. The same gp.Latent with a Bernoulli likelihood does classification.
import pymc as pm
cxs=(cx-cx.mean())/cx.std()
with pm.Model() as mod:
ell=pm.Gamma("ell",3,1); eta=pm.HalfNormal("eta",2)
cov=eta**2*pm.gp.cov.ExpQuad(1, ls=ell)
gpl=pm.gp.Latent(cov_func=cov)
f=gpl.prior("f", X=cxs[:,None])
pm.Poisson("y", mu=w*pm.math.exp(f), observed=cy)
idata=pm.sample(700, tune=1200, chains=4, target_accept=0.95, random_seed=4, progressbar=False)
fpost=idata.posterior["f"].stack(s=("chain","draw")).values # (nbins, S)
lam_pm=(np.exp(fpost)).mean(1); lam_pm_lo=np.percentile(np.exp(fpost),2.5,axis=1); lam_pm_hi=np.percentile(np.exp(fpost),97.5,axis=1)
lam_ls=np.exp(GL.fit_predict(cx,cy,cx,GL.lik_poisson,w,kc,[np.log(15),0.])["mean"])
print(f"intensity agreement at bin centres: max |PyMC - Laplace| = {np.abs(lam_pm-lam_ls).max():.2f} events/yr")
fig,ax=plt.subplots(figsize=(9,4.2))
ax.bar(cx,cy/w,width=w*0.9,color=GREY,alpha=.3,label="binned rate")
ax.fill_between(cx,lam_pm_lo,lam_pm_hi,color=GREEN,alpha=.2,label="PyMC 95%")
ax.plot(cx,lam_pm,color=GREEN,lw=2,label="PyMC gp.Latent"); ax.plot(cx,lam_ls,color=BLUE,lw=1.6,ls="--",label="from-scratch Laplace")
ax.set_xlabel("year"); ax.set_ylabel("disasters per year"); ax.set_title("LGCP intensity: PyMC vs from-scratch Laplace"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
gap = float(np.abs(lam_pm-lam_ls).max()); peak = float(lam_pm.max())
print("The sampled latent GP and the Laplace approximation trace the same fall at 1890. The absolute gap of")
print("%.2f events/yr is best read against the peak intensity of %.2f -- about %.0f%% -- and it sits where the"
% (gap, peak, 100*gap/peak))
print("intensity is highest and the counts are largest, not at the flat modern end. Laplace is a good, cheap")
print("stand-in for full sampling when the posterior over f is near-Gaussian, as it is here; 'good' at the")
print("10-20% level on the peak, not exact.")
g++ not available, if using conda: `conda install gxx`
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [ell, eta, f_rotated_]
Sampling 4 chains for 1_200 tune and 700 draw iterations (4_800 + 2_800 draws total) took 18 seconds.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
intensity agreement at bin centres: max |PyMC - Laplace| = 0.52 events/yr
The sampled latent GP and the Laplace approximation trace the same fall at 1890. The absolute gap of 0.52 events/yr is best read against the peak intensity of 3.30 -- about 16% -- and it sits where the intensity is highest and the counts are largest, not at the flat modern end. Laplace is a good, cheap stand-in for full sampling when the posterior over f is near-Gaussian, as it is here; 'good' at the 10-20% level on the peak, not exact.
5. Summary¶
One prior over functions, two models. Through a logit link the GP is a classifier — a nonlinear log-odds surface that recovered a non-monotone probability logistic regression could not touch, and on the Pima data drew a smooth boundary with honest uncertainty in the sparse corners of feature space. Through a log link it is a log-Gaussian Cox process — a smooth event intensity that reproduced the famous post-1890 fall in coal-mining disasters with a credible band.
The non-Gaussian likelihood cost us the closed form, so we used the Laplace approximation (Newton to the posterior mode, Gaussian curvature, marginal likelihood for the kernel) — and PyMC's sampled gp.Latent confirmed the intensity to within a fraction of an event per year, so the cheap approximation was faithful. The frequentist analogs — logistic regression and the kernel intensity estimate — trace the same curves but without the calibrated uncertainty a probability model provides.
These links tie the Gaussian-process arc back to the earlier models: GP classification is logistic regression with a nonparametric predictor, and the log-Gaussian Cox process is Poisson regression with a GP log-rate — the count and binary GLMs, freed of their linear form. Next the arc turns to nonparametric survival and hazards, and then to Polya trees.