Bayesian Penalised Splines & Additive Models

Too Many Bumps, Held Together by a Penalty

A Gaussian process places its prior on the function directly. A penalised spline gets to the same place by a different and much cheaper route: lay down a generous B-spline basis — deliberately more bumps than the data can support — and then penalise the roughness of the coefficients. The second-difference penalty ties each coefficient to its neighbours, so the fit cannot wiggle freely between knots. Read as a prior, that penalty is a random walk on the coefficients, which makes the whole construction Bayesian without changing a line of the algebra.

Inferring the smoothness instead of choosing it

Everything then turns on how hard the tie is pulled. On the airquality data — 111 complete days, log ozone against temperature — the two extremes are drawn explicitly: a tiny penalty chases noise, a huge one collapses to a straight line. The Bayesian version does not select between them. It puts a prior on the smoothing variance and integrates over it, returning both the curve and a credible band, with the smoothness inferred rather than tuned.

Against cross-validation, measured rather than asserted

The frequentist counterpart tunes the same penalty by generalised cross-validation — which estimates, in closed form and without actually leaving anything out, how well each candidate $\lambda$ would predict a held-out point, and picks the minimiser. The comparison with the Bayesian answer is worth making quantitative rather than qualitative. GCV selects $\lambda=38$; the posterior implies an effective penalty of $\sigma^2/\mathbb{E}[\tau^2]=29$. That is a factor of 1.3 apart — and the two fitted curves differ by at most 0.069 log-ppb, 1.3% of the response range. The reason the gap in $\lambda$ matters so little is structural: the penalty enters only as a ratio inside a matrix inverse, so a factor of a few in $\lambda$ moves the curve far less than it moves $\lambda$. One approach tunes the penalty, the other averages over it, and here they agree.

One smooth per covariate

Stacking one penalised spline per covariate gives an additive model, each term with its own inferred smoothness, all estimated jointly by backfitting inside the Gibbs sampler. The three partial effects on ozone are physically legible: temperature raises it strongly and nonlinearly, wind lowers it (calm days trap pollution), and solar radiation raises it then flattens. R's mgcv recovers the same three shapes with REML, at 70.7% deviance explained and roughly 2 effective degrees of freedom per term — semiparametric regression with no functional forms assumed anywhere.

The same method elsewhere

The penalised spline built here from the basis up reappears as a predictor in Gaussian Processes & Splines — the Bayesian Kernel View, where a sum of these smooths — an additive model — comes within 0.01 AUC of the tree ensembles on credit default while every effect stays a readable curve with a credible band. The same roughness penalty, the same random-walk prior, scored against machine-learning baselines instead of simulated curves.

Where NUTS Struggles and Gibbs Does Not

The PyMC cross-check produced the more interesting caveat, and it is about the sampler rather than the answer. Encoding the penalty as a random-walk Potential reproduces the from-scratch fit closely — the two curves differ by 0.051 log-ppb, 1.0% of the response range — but 6.2% of NUTS draws diverge. That is the classic funnel geometry of a hierarchical variance, and it is precisely what the conjugate Gibbs sampler never feels, because it draws $\beta$ and $\tau^2$ from their exact full conditionals rather than exploring the joint density with a gradient. The agreement between the two is the evidence that the divergences have not distorted the posterior mean here; a non-centred reparameterisation would be the fix if they had.

Notebooks

Downloads

P-Spline Module — Source Code

"""
pspline.py -- BAYESIAN PENALISED SPLINES (P-splines) and ADDITIVE MODELS (from scratch).

Backs the notebooks in  "Bayesian Penalised Splines & Additive Models".

Two routes to a flexible curve. The variable-selection arc's free-knot notebook chooses the
NUMBER and LOCATION of a few knots by trans-dimensional MCMC. The PENALISED-spline route
(Eilers & Marx 1996; Bayesian version: Lang & Brezger 2004) does the opposite: lay down MANY
equally spaced B-spline basis functions -- far more than needed -- and control smoothness not by
the knots but by a PENALTY on the coefficients. As a Bayesian prior this penalty is a random walk:

    y_i = sum_j B_j(x_i) beta_j + eps_i,   eps ~ N(0, sigma^2)
    beta_j - 2 beta_{j-1} + beta_{j-2} ~ N(0, tau^2)          (2nd-order random-walk prior)

i.e. beta ~ N(0, tau^2 (D'D)^{-1}) with D the 2nd-difference matrix. The smoothing parameter is
the ratio sigma^2/tau^2, and -- crucially -- it is INFERRED from the data through tau^2's own
prior, so there is nothing to cross-validate. The full conditionals are all conjugate (a Gaussian
Markov random field):

    beta   ~ N( (B'B/sigma^2 + D'D/tau^2)^{-1} B'y/sigma^2,  (B'B/sigma^2 + D'D/tau^2)^{-1} )
    sigma^2 ~ InvGamma(a + n/2,      b + 1/2 ||y - B beta||^2)
    tau^2   ~ InvGamma(a + rank(D)/2, b + 1/2 beta' D'D beta)

ADDITIVE MODELS (GAMs) stack one penalised spline per covariate, y = beta0 + sum_k f_k(x_k),
each with its own smoothing variance, and are sampled by the same Gibbs sweep applied to partial
residuals (backfitting inside Gibbs). A P-spline is also exactly a Gaussian process with a
particular kernel, so this connects straight back to the GP-regression notebook.
"""

import numpy as np
from scipy.interpolate import BSpline
from scipy.linalg import cholesky, cho_solve, solve_triangular


# --------------------------------------------------------------------------- #
#  B-spline design matrix (equally spaced knots) + difference penalty           #
# --------------------------------------------------------------------------- #

def make_knots(xl, xr, ndx, deg):
    dx = (xr - xl) / ndx
    return xl + dx * np.arange(-deg, ndx + deg + 1)

def bspline_design(x, knots, deg):
    nb = len(knots) - deg - 1; B = np.empty((len(x), nb))
    for j in range(nb):
        c = np.zeros(nb); c[j] = 1.0
        B[:, j] = BSpline(knots, c, deg, extrapolate=True)(x)
    return B

def diff_penalty(nb, order=2):
    return np.diff(np.eye(nb), order, axis=0)          # (nb-order, nb); K = D'D


def _draw_beta(BtB, Bty, K, sig2, tau2, rng):
    nb = BtB.shape[0]
    Prec = BtB / sig2 + K / tau2 + 1e-8 * np.eye(nb)
    L = cholesky(Prec, lower=True)
    mean = cho_solve((L, True), Bty / sig2)
    return mean + solve_triangular(L.T, rng.standard_normal(nb))


# --------------------------------------------------------------------------- #
#  single penalised spline                                                      #
# --------------------------------------------------------------------------- #

def pspline_gibbs(x, y, rng, ndx=25, deg=3, penord=2, draws=2500, burn=1000,
                  grid=None, a0=1e-3, b0=1e-3):
    x = np.asarray(x, float); y = np.asarray(y, float); n = len(y)
    xl, xr = x.min(), x.max(); knots = make_knots(xl, xr, ndx, deg)
    B = bspline_design(x, knots, deg); nb = B.shape[1]
    K = diff_penalty(nb, penord); K = K.T @ K
    if grid is None: grid = np.linspace(xl, xr, 200)
    Bg = bspline_design(grid, knots, deg)
    BtB = B.T @ B; Bty = B.T @ y; rankD = nb - penord
    tau2 = 1.0; sig2 = float(np.var(y)); FIT = np.empty((draws, len(grid)))
    TAU = np.empty(draws)
    for it in range(draws + burn):
        beta = _draw_beta(BtB, Bty, K, sig2, tau2, rng)
        resid = y - B @ beta
        sig2 = 1.0 / rng.gamma(a0 + n / 2, 1.0 / (b0 + 0.5 * resid @ resid))
        tau2 = 1.0 / rng.gamma(a0 + rankD / 2, 1.0 / (b0 + 0.5 * beta @ K @ beta))
        if it >= burn:
            FIT[it - burn] = Bg @ beta; TAU[it - burn] = tau2
    return dict(grid=grid, fit=FIT, tau2=TAU, sig2=sig2)


# --------------------------------------------------------------------------- #
#  additive model (GAM) -- one penalised spline per covariate                   #
# --------------------------------------------------------------------------- #

def gam_gibbs(Xcols, y, rng, ndx=15, deg=3, penord=2, draws=2500, burn=1000,
              grids=None, a0=1e-3, b0=1e-3):
    """Xcols: list of 1-D covariate arrays. Returns per-term partial-effect draws on `grids`
    (centred so each smooth sums to zero; the intercept carries the level)."""
    y = np.asarray(y, float); n = len(y); p = len(Xcols)
    knots = [make_knots(x.min(), x.max(), ndx, deg) for x in Xcols]
    Braw = [bspline_design(np.asarray(x, float), knots[k], deg) for k, x in enumerate(Xcols)]
    cmean = [b.mean(0) for b in Braw]
    B = [Braw[k] - cmean[k] for k in range(p)]                  # centred bases (identifiability)
    K = [diff_penalty(b.shape[1], penord) for b in B]; K = [d.T @ d for d in K]
    BtB = [B[k].T @ B[k] for k in range(p)]
    if grids is None:
        grids = [np.linspace(x.min(), x.max(), 200) for x in Xcols]
    Bg = [bspline_design(grids[k], knots[k], deg) - cmean[k] for k in range(p)]
    beta = [np.zeros(b.shape[1]) for b in B]; b0i = float(y.mean())
    tau2 = [1.0] * p; sig2 = float(np.var(y)); rankD = [b.shape[1] - penord for b in B]
    PART = [np.empty((draws, len(grids[k]))) for k in range(p)]
    B0 = np.empty(draws)
    for it in range(draws + burn):
        fitsum = sum(B[k] @ beta[k] for k in range(p))
        b0i = (y - fitsum).mean() + rng.standard_normal() * np.sqrt(sig2 / n)
        for k in range(p):
            r = y - b0i - sum(B[j] @ beta[j] for j in range(p) if j != k)
            beta[k] = _draw_beta(BtB[k], B[k].T @ r, K[k], sig2, tau2[k], rng)
            tau2[k] = 1.0 / rng.gamma(a0 + rankD[k] / 2, 1.0 / (b0 + 0.5 * beta[k] @ K[k] @ beta[k]))
        resid = y - b0i - sum(B[k] @ beta[k] for k in range(p))
        sig2 = 1.0 / rng.gamma(a0 + n / 2, 1.0 / (b0 + 0.5 * resid @ resid))
        if it >= burn:
            B0[it - burn] = b0i
            for k in range(p):
                PART[k][it - burn] = Bg[k] @ beta[k]
    return dict(partial=PART, grids=grids, intercept=B0, sig2=sig2)


# --------------------------------------------------------------------------- #
#  simulation for validation                                                    #
# --------------------------------------------------------------------------- #

def simulate_additive(n, funcs, ranges, noise_sd, rng, intercept=0.0):
    """y = intercept + sum_k f_k(x_k) + noise; funcs = list of callables, ranges = list of (lo,hi)."""
    Xcols = [rng.uniform(lo, hi, n) for (lo, hi) in ranges]
    mu = intercept + sum(f(x) for f, x in zip(funcs, Xcols))
    y = mu + noise_sd * rng.standard_normal(n)
    return Xcols, y

References

• Eilers, P. H. C. & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science 11(2), 89–121. — P-splines, the basis-and-penalty construction used here
• Lang, S. & Brezger, A. (2004). Bayesian P-splines. JCGS 13(1), 183–212. — the random-walk prior reading of the penalty, and the Gibbs sampler
• Ruppert, D., Wand, M. P. & Carroll, R. J. (2003). Semiparametric Regression. Cambridge University Press. — penalised splines as mixed models, and additive extensions
• Hastie, T. & Tibshirani, R. (1990). Generalized Additive Models. Chapman & Hall. — backfitting and the additive structure
• Wood, S. N. (2011). Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. JRSS-B 73(1), 3–36. — the REML selection mgcv uses in the R engine
• Betancourt, M. & Girolami, M. (2015). Hamiltonian Monte Carlo for hierarchical models. In Current Trends in Bayesian Methodology. — the funnel geometry behind the divergences reported here