Free-knot splines — a curve with an unknown number of knots

Bayesian variable selection, Part 7 — DiMatteo, Genovese & Kass (2001)

The last unknown-dimension problem is nonparametric curve fitting. A cubic spline bends at its knots; more knots mean more flexibility. But how many knots, and where? Congdon (2005) Example 3.9 fits a cubic spline with the knot locations unknown but their number fixed at five. The reversible-jump extension (DiMatteo-Genovese-Kass 2001) lets the number of knots itself be inferred — the spline analogue of Parts 5 and 6.

The data is Congdon's 49-point curve: a response $Y$ that is nearly flat, then rises to a sharp peak near $X=895$, then falls back — the kind of localized feature that makes knot placement matter. We fit it with a from-scratch free-knot RJMCMC, and cross-check the recovered curve against a PyMC penalized spline and an R (mgcv) penalized regression spline.

The model and the reversible jump

A cubic spline in the truncated-power basis with $k$ knots $\xi_1<\dots<\xi_k$: $$ y_i = b_0 + b_1 x_i + b_2 x_i^2 + b_3 x_i^3 + \sum_{j=1}^{k} b_{3+j}\,(x_i-\xi_j)_+^3 + \varepsilon_i,\qquad \varepsilon_i\sim\mathcal N(0,\sigma^2), $$ with both the number of knots $k$ and their positions unknown. Exactly as in the change-point sampler, the regression coefficients and $\sigma^2$ are integrated out under a Zellner $g$-prior, leaving a closed-form marginal likelihood for a knot configuration (a function of the model's $R^2$ and its size). The reversible-jump moves then change dimension without a Jacobian:

• Birth — add a knot at a random empty grid position: accept with $\min\!\big(1,\ e^{\Delta\mathcal L}\,\lambda/(k+1)\big)$;
• Death — remove a random knot: $\min\!\big(1,\ e^{\Delta\mathcal L}\,k/\lambda\big)$;
• Move — relocate a knot (fixed dimension): $\min(1,e^{\Delta\mathcal L})$.

Hyperparameters (fixed). Knots live on a grid of 80 candidate interior positions; a $\text{Poisson}(\lambda=4)$ prior on their number supplies the parsimony penalty; $g=n$ (unit-information) in the marginal likelihood. $x$ is rescaled to $[0,1]$ for numerical stability. The fitted curve is the $g$-prior posterior mean, averaged over all sampled knot configurations (Bayesian model averaging over splines). Engine: freeknot.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import freeknot                                                     # from-scratch reversible-jump free-knot spline

d = pd.read_csv("spline_data.csv"); x = d["X"].values; y = d["Y"].values
res = freeknot.rj_spline(x, y, ncand=80, lam=4.0, niter=100000, burn=30000, seed=1)

kp = res["kpost"]
print("posterior on the number of knots k:")
for k in range(4, 9): print("  k=%d : %.3f" % (k, kp[k]))
print("\nE[k] = %.1f   mode = %d   |   fit RMSE = %.3f   |   knots concentrate near X = 895 (the peak)"
      % ((np.arange(len(kp))*kp).sum(), np.argmax(kp), np.sqrt(np.mean((res["fit"]-y)**2))))

# --- Fig: data, fitted curve, and where the sampler puts its knots ---
o = np.argsort(x)
fig, ax = plt.subplots(figsize=(9, 4.8))
ax.scatter(x, y, color='0.6', s=12, label='data')
ax.plot(x[o], res["fit"][o], color='firebrick', lw=2, label='free-knot fit (posterior mean)')
ax2 = ax.twinx()
ax2.fill_between(res["cand"], 0.0, res["loc"], color='steelblue', alpha=.25, step='mid')
ax2.plot(res["cand"], res["loc"], color='steelblue', lw=1.2, label='knot-location posterior')
ax2.set_ylim(bottom=0)
ax2.set_ylabel('knot-location posterior', color='steelblue')
ax2.tick_params(axis='y', colors='steelblue')
ax.set_xlabel('X'); ax.set_ylabel('Y')
ax.set_title('Free-knot spline: fitted curve and where the sampler places its knots')
ax.grid(alpha=.25)
h1, l1 = ax.get_legend_handles_labels(); h2, l2 = ax2.get_legend_handles_labels()
ax.legend(h1 + h2, l1 + l2, loc='upper left', framealpha=.9)
plt.tight_layout(); plt.show()

posterior on the number of knots k:
  k=4 : 0.002
  k=5 : 0.437
  k=6 : 0.380
  k=7 : 0.140
  k=8 : 0.031

E[k] = 5.8   mode = 5   |   fit RMSE = 0.014   |   knots concentrate near X = 895 (the peak)

Reading the free-knot fit

• Five knots, and it says where. The posterior on the number of knots peaks at $k=5$ (0.44), with $k=6$ almost as likely — $\mathbb E[k]\approx5.8$. Congdon's fixed choice of five knots turns out to be well-supported; the difference is that here it is inferred, with honest uncertainty, rather than assumed.
• Knots migrate to the curvature. The location posterior (blue, right axis) is concentrated around $X\approx895$ — the sharp peak — and sparse in the flat regions. A free-knot spline adapts: it spends its knots where the function bends and none where it is straight, which is exactly why it fits the peak so tightly (RMSE 0.014).
• This is selection for curves. The "variables" being selected are the knots; the posterior over their number and location is the curve-fitting version of the inclusion probabilities from Part 1.

import pymc as pm, patsy                                            # cross-check: a penalized (many-knot) spline
B = np.asarray(patsy.dmatrix("bs(x, df=16, degree=3, include_intercept=True) - 1", {"x": x}))
with pm.Model() as mod:
    sb  = pm.HalfNormal("sb", 1.0)
    b   = pm.GaussianRandomWalk("b", sigma=sb, shape=B.shape[1], init_dist=pm.Normal.dist(0, 5))  # RW smoothness
    sig = pm.HalfNormal("sig", 0.5)
    pm.Normal("y", pm.math.dot(B, b), sig, observed=y)
    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.95, random_seed=1, progressbar=False)
pm_fit = (idata.posterior["b"].values.reshape(-1, B.shape[1]) @ B.T).mean(0)
pm_rmse = np.sqrt(np.mean((pm_fit-y)**2)); fk_rmse = np.sqrt(np.mean((res["fit"]-y)**2))
print("PyMC penalized spline (16 knots + random-walk smoothness): fit RMSE = %.3f" % pm_rmse)
print("-> recovers the same curve, but oversmooths the sharp peak that free knots capture.")

# --- Fig: posterior over the number of knots, and free-knot vs penalized spline ---
o = np.argsort(x)
kp = res["kpost"]; ks = np.arange(3, 10)
fig, (axA, axB) = plt.subplots(1, 2, figsize=(12.4, 4.3))

axA.bar(ks, kp[ks], color='firebrick', alpha=.85)
axA.set_xlabel('number of knots k'); axA.set_ylabel('posterior P(k)')
axA.set_title('Posterior over the number of knots')
axA.grid(alpha=.25, axis='y')

axB.scatter(x, y, color='0.6', s=12, label='data')
axB.plot(x[o], res["fit"][o], color='firebrick', lw=2, label=f'free-knot RJMCMC (RMSE {fk_rmse:.3f})')
axB.plot(x[o], pm_fit[o], color='steelblue', lw=2, ls='--', label=f'PyMC penalized (RMSE {pm_rmse:.3f})')
axB.set_xlabel('X'); axB.set_ylabel('Y')
axB.set_title('Free-knot vs penalized spline')
axB.grid(alpha=.25); axB.legend(loc='upper left', framealpha=.9)

plt.tight_layout(); plt.show()

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\threadpoolctl.py:1226: RuntimeWarning: 
Found Intel OpenMP ('libiomp') and LLVM OpenMP ('libomp') loaded at
the same time. Both libraries are known to be incompatible and this
can cause random crashes or deadlocks on Linux when loaded in the
same Python program.
Using threadpoolctl may cause crashes or deadlocks. For more
information and possible workarounds, please see
    https://github.com/joblib/threadpoolctl/blob/master/multiple_openmp.md

  warnings.warn(msg, RuntimeWarning)
Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [sb, b, sig]

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 6 seconds.

PyMC penalized spline (16 knots + random-walk smoothness): fit RMSE = 0.067
-> recovers the same curve, but oversmooths the sharp peak that free knots capture.

Results — two philosophies of flexibility

• Same curve, different route. The PyMC penalized spline (many fixed knots, coefficients tied by a random-walk smoothness prior) and R's mgcv (next notebook) recover the same overall shape as the free-knot RJMCMC — a reassuring agreement across implementations.
• But adaptivity shows at the peak. The penalized spline spreads a single smoothness everywhere, so it both undershoots the sharp peak and wiggles across the flat baseline (RMSE 0.067); the free-knot sampler, by placing knots only where the curve bends, captures the peak cleanly and leaves the baseline flat (RMSE 0.014). Selecting where the flexibility goes beats applying it uniformly when the function has a localized feature — the classic argument for free-knot over penalized splines.
• The unknown-dimension trilogy, complete. Change-points (Part 5), mixture components (Part 6) and now spline knots — three "how many?" problems, all solved by the same move: put a conjugate prior on the within-piece parameters, integrate them out, and let a reversible-jump (or collapsed) sampler roam over the discrete structure.

References

• DiMatteo, I., Genovese, C. R. & Kass, R. E. (2001). Bayesian curve-fitting with free-knot splines. Biometrika 88, 1055–1071.
• Green, P. J. (1995). Reversible jump MCMC computation and Bayesian model determination. Biometrika 82, 711–732.
• Congdon, P. (2005). Bayesian Models for Categorical Data, Example 3.9. Wiley.

Data: spline_data.csv — Congdon's 49-point curve ($Y$ vs $X$), a flat baseline with a sharp peak near $X=895$. Next: freeknot_R.ipynb — the R cross-check (mgcv penalized regression spline).