GP Classification & Log-Gaussian Cox Processes

One Latent Function, Two Links

The previous project put a GP prior on a function observed with Gaussian noise, where everything stayed closed form. Push the same latent function through a link and that convenience disappears — but the modelling reach expands sharply. A logit link turns a GP into a classifier; a log link turns it into the intensity of a log-Gaussian Cox process, a prior over event rates in time or space. The posterior is no longer Gaussian, so it is handled here two ways: a Laplace approximation from scratch, and full HMC sampling in PyMC as the check on it.

Recovering a probability that is not monotone

The classification case is validated against a known truth first. On a simulated non-monotone probability curve the GP recovers $p(x)$ with RMSE 0.062 against logistic regression's 0.212 — logistic regression has log-odds linear in $x$, so it is structurally confined to a single monotone S and fits a nearly flat line through a wave. R's smooth binomial GAM reproduces the result independently at RMSE 0.084.

Where the GP does not win, stated plainly

The Pima diabetes data then supplies the honest negative. Fitted on glucose and BMI, the GP's accuracy is 0.776 against logistic regression's 0.773 — a dead heat, because the problem is close to linearly separable in those two features and the fitted length-scale of 5.5 standard deviations tells you so. The notebooks say this plainly rather than manufacturing a win. What the GP adds is not accuracy but calibration of ignorance: its uncertainty grows in the sparse corners of the feature space where any boundary is a guess, while logistic regression stays equally confident everywhere, including where it has seen almost nobody.

The same machinery as an event-rate model

Swapping the link for a log turns the same machinery into a Cox process. Applied to the classic coal-mining disaster record — 191 events from 1851 to 1962 — the fitted intensity falls from about 3.0 disasters per year around 1860 to 0.71 by 1950, with the sharpest decline near 1890, the period of the major mine-safety legislation. The Cox process adds what a kernel intensity estimate cannot: a credible band around the whole curve.

The same method elsewhere

The latent-GP-through-a-link construction used here for Pima is applied to a much larger classification problem in Gaussian Processes & Splines — the Bayesian Kernel View, where a GP classifier reaches AUC 0.742 on credit default from just 1,200 of 21,000 training rows — beating a logistic regression fitted to all of them.

A Factor of Two in the Exposure

That coal intensity is where the two engines have to be reconciled, because a Cox process is unusually easy to get wrong by exactly a factor of two.

The cause is a double correction on the exposure. The model is fitted with a $\log w$ offset for two-year bins, and the code then divided the prediction by $w$ a second time. An offset supplied through gam()'s offset argument is not carried into predict.gam(), so $\exp(\hat\eta)$ was already a per-year rate. Verified directly in R: $\exp(\hat\eta)=3.055$ at 1860 while fitted() gives 6.109, exactly $2\times$ that, confirming the offset lives in one and not the other. The observed raw rate is 3.154/yr over 1851–1890, which settles which figure was right. It was visible in the figure too — the fitted curve ran at half the height of the binned-rate bars it was supposed to pass through. Corrected, R now reports 3.1 and 0.71, and the claim of agreement is finally true.

Two further diagnostics

Two smaller corrections came with it. The R panel took its $y$-limits from the credible band alone, so the tallest binned-rate bars ran off the top of the frame. And the PyMC cross-check, which had concluded that sampling and Laplace "give the same intensity," rests on a maximum gap of 0.52 events/yr — about 16% of the peak intensity of 3.30, concentrated exactly where the rate is highest. Laplace is a good cheap stand-in here; good at the 10–20% level on the peak, not exact.

Notebooks

Downloads

GP-Link Module — Source Code

"""
gplink.py -- GP CLASSIFICATION and LOG-GAUSSIAN COX PROCESSES (from scratch).

Backs the notebooks in  "GP Classification & Log-Gaussian Cox Processes".

A Gaussian-process regression puts a GP prior on a real-valued function. Push that latent
function through a LINK and it models other things:

  * GP CLASSIFICATION -- binary/count outcomes via a logit link:
        f ~ GP(0,k),   y_i | f ~ Bernoulli(sigmoid(f(x_i)))   (or Binomial)
    the GP is 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),   f ~ GP(0,k),   counts in a bin ~ Poisson(lambda * width)
    the smooth, positive intensity of events over time/space, with uncertainty.

Now the likelihood is non-Gaussian, so the posterior over f is not closed form. The classic
from-scratch solution is the LAPLACE APPROXIMATION (Rasmussen & Williams, Algorithm 3.1): Newton-
iterate to the posterior mode of f, then approximate the posterior by the Gaussian with that mode
and the curvature there. The mode step and the predictive variance are computed stably through the
Cholesky factor of B = I + W^{1/2} K W^{1/2}, where W is the (diagonal) negative Hessian of the
log-likelihood. The Laplace log marginal likelihood then selects the kernel hyperparameters, and
elliptical slice sampling gives an exact-posterior cross-check.

The kernels are D-dimensional (isotropic RBF over standardised inputs), covering the 2-D Pima
classification surface and the 1-D coal-mining intensity alike.
"""

import numpy as np
from scipy.linalg import cholesky, cho_solve, solve_triangular
from scipy.optimize import minimize
from scipy.special import expit


# --------------------------------------------------------------------------- #
#  kernels (D-dimensional)                                                      #
# --------------------------------------------------------------------------- #

def sqdist(Xa, Xb):
    Xa = np.atleast_2d(Xa); Xb = np.atleast_2d(Xb)
    return np.sum(Xa ** 2, 1)[:, None] + np.sum(Xb ** 2, 1)[None, :] - 2 * Xa @ Xb.T

def k_rbf(Xa, Xb, ell, eta):
    return eta * np.exp(-0.5 * sqdist(Xa, Xb) / ell ** 2)


# --------------------------------------------------------------------------- #
#  likelihoods: return (loglik, gradient d/df, W = -d2/df2 diagonal)            #
# --------------------------------------------------------------------------- #

def lik_bernoulli(f, y, extra=None):
    p = expit(f); ll = np.sum(y * np.log(p + 1e-12) + (1 - y) * np.log(1 - p + 1e-12))
    return ll, y - p, p * (1 - p)

def lik_binomial(f, y, n):
    p = expit(f); ll = np.sum(y * f - n * np.log1p(np.exp(f)))
    return ll, y - n * p, n * p * (1 - p)

def lik_poisson(f, y, e):
    lam = e * np.exp(f); ll = np.sum(y * (np.log(e + 1e-12) + f) - lam)
    return ll, y - lam, lam


# --------------------------------------------------------------------------- #
#  Laplace approximation (Rasmussen-Williams Algorithm 3.1)                      #
# --------------------------------------------------------------------------- #

def laplace_mode(K, y, likfun, extra, maxit=100, tol=1e-8):
    """find the posterior mode of the latent f and the Laplace log marginal likelihood."""
    n = len(y); f = np.zeros(n); obj_old = -np.inf; a = np.zeros(n); L = np.eye(n)
    for _ in range(maxit):
        _, g, W = likfun(f, y, extra); sW = np.sqrt(np.maximum(W, 1e-12))
        B = np.eye(n) + sW[:, None] * K * sW[None, :]
        L = cholesky(B, lower=True)
        b = W * f + g
        a = b - sW * cho_solve((L, True), sW * (K @ b))
        f = K @ a
        ll, _, _ = likfun(f, y, extra)
        obj = -0.5 * a @ f + ll - np.log(np.diag(L)).sum()
        if abs(obj - obj_old) < tol:
            break
        obj_old = obj
    return f, a, L, sW, obj

def fit_predict(X, y, Xs, likfun, extra, kernel, theta0, bounds=None):
    """fit kernel hyperparameters by maximising the Laplace log marginal likelihood, then return
    the predictive latent mean and variance at Xs. kernel(theta, Xa, Xb) -> covariance."""
    X = np.asarray(X, float); Xs = np.asarray(Xs, float)
    if X.ndim == 1: X = X[:, None]
    if Xs.ndim == 1: Xs = Xs[:, None]
    def neg_obj(theta):
        K = kernel(theta, X, X) + 1e-6 * np.eye(len(y))
        return -laplace_mode(K, y, likfun, extra)[4]
    res = minimize(neg_obj, theta0, method="L-BFGS-B", bounds=bounds)
    th = res.x; K = kernel(th, X, X) + 1e-6 * np.eye(len(y))
    f, a, L, sW, obj = laplace_mode(K, y, likfun, extra)
    Ks = kernel(th, Xs, X); fs_mean = Ks @ a
    v = solve_triangular(L, sW[:, None] * Ks.T, lower=True)
    fs_var = np.maximum(np.diag(kernel(th, Xs, Xs)) - np.sum(v ** 2, 0), 0.0)
    return dict(theta=th, f=f, mean=fs_mean, var=fs_var, logml=obj)


# --------------------------------------------------------------------------- #
#  elliptical slice sampling -- exact-posterior cross-check                      #
# --------------------------------------------------------------------------- #

def ess_sample(K, y, likfun, extra, draws=400, burn=200, rng=None):
    """draw latent-function samples f | y by elliptical slice sampling (Murray et al. 2010)."""
    n = len(y); Lk = cholesky(K + 1e-6 * np.eye(n), lower=True)
    f = Lk @ rng.standard_normal(n); out = []
    for it in range(draws + burn):
        nu = Lk @ rng.standard_normal(n)
        ll0 = likfun(f, y, extra)[0] + np.log(rng.random())
        th = rng.uniform(0, 2 * np.pi); lo, hi = th - 2 * np.pi, th
        while True:
            fp = f * np.cos(th) + nu * np.sin(th)
            if likfun(fp, y, extra)[0] > ll0:
                f = fp; break
            if th < 0: lo = th
            else: hi = th
            th = rng.uniform(lo, hi)
        if it >= burn: out.append(f.copy())
    return np.array(out)


# --------------------------------------------------------------------------- #
#  helpers: link transforms + LGCP binning                                      #
# --------------------------------------------------------------------------- #

def predictive_prob(mean, var):
    """MacKay probit approximation to E[sigmoid(f)] under f ~ N(mean, var)."""
    return expit(mean / np.sqrt(1 + np.pi * var / 8))

def bin_counts(events, lo, hi, nbins):
    """bin event times into counts; returns bin centres, counts, and bin width (exposure)."""
    edges = np.linspace(lo, hi, nbins + 1); width = edges[1] - edges[0]
    counts, _ = np.histogram(events, bins=edges)
    centres = 0.5 * (edges[:-1] + edges[1:])
    return centres, counts.astype(float), width

References

• Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning, ch. 3. MIT Press. — GP classification and the Laplace approximation implemented here
• Møller, J., Syversveen, A. R. & Waagepetersen, R. P. (1998). Log Gaussian Cox processes. Scandinavian Journal of Statistics 25(3), 451–482. — the intensity model
• Murray, I., Adams, R. P. & MacKay, D. J. C. (2010). Elliptical slice sampling. AISTATS. — the sampler provided alongside Laplace in the module
• Jarrett, R. G. (1979). A note on the intervals between coal-mining disasters. Biometrika 66(1), 191–193. — the dataset, and the published event rate the intensity is checked against
• Smith, J. W., Everhart, J. E. et al. (1988). Using the ADAP learning algorithm to forecast the onset of diabetes mellitus. Proc. Symp. Computer Applications in Medical Care. — the Pima data
• Wood, S. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.). CRC Press. — mgcv, and the offset semantics at the centre of the correction