Free-Knot Splines — an Unknown Number of Knots

Python · R · Part 7 of the variable-selection arc — unknown model dimension

Overview

Part 7 of the variable-selection arc: nonparametric curve fitting, the third unknown-dimension problem. A cubic spline bends at its knots — more knots mean more flexibility — but how many knots, and where? Congdon (2005) 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 the change-point (Part 5) and mixture-component (Part 6) samplers. The data is Congdon's 49-point curve: a response nearly flat, then rising to a sharp peak near $X=895$ , then falling — the kind of localized feature that makes knot placement matter. Fit with a from-scratch free-knot RJMCMC, cross-checked against a PyMC penalized spline and R's mgcv.

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 k,\ \{\xi_j\}\ \text{unknown}

Free-knot cubic splines by reversible jump

A cubic spline in the truncated-power basis with $k$ knots $\xi_1<\dots<\xi_k$ adds a term $b_{3+j}(x_i-\xi_j)_+^3$ per knot, 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 with no Jacobian: birth (add a knot at a random empty grid position), death (remove a random knot), and move (relocate one at fixed dimension), accepted on the change in marginal log-likelihood times the Poisson prior ratio. 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. The fitted curve is the $g$ -prior posterior mean averaged over all sampled knot configurations — Bayesian model averaging over splines. Engine: a self-contained freeknot.py.

Results

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]\approx 5.8$ ) — Congdon's fixed choice of five is well-supported, but here it is inferred, with honest uncertainty, rather than assumed. Knots migrate to the curvature: the location posterior is concentrated around $X\approx 895$ — the sharp peak — and sparse in the flat regions. A free-knot spline adapts, spending its knots where the function bends and none where it is straight, which is why it fits the peak so tightly (RMSE 0.014). The "variables" being selected are the knots; the posterior over their number and location is the curve-fitting version of Part 1's inclusion probabilities. Two philosophies of flexibility: the PyMC penalized spline (many fixed knots, coefficients tied by a single random-walk smoothness prior) spreads that smoothness everywhere, so it undershoots the sharp peak while wiggling across the flat baseline (RMSE 0.067); R's mgcv, whose penalty is tuned by REML, traces the same shape far more closely (RMSE 0.025) but spends ~19 effective degrees of freedom to do it. The free-knot sampler fits the peak tightest of all (RMSE 0.014) with only ~5 knots, by placing them where the curve bends rather than smoothing uniformly — the classic argument for free-knot over penalized splines when the function has a localized feature. This completes the unknown-dimension trilogy — change-points, mixture components, and now spline knots — three "how many?" problems, all solved by the same move: conjugate prior on the within-piece parameters, integrate out, and let a reversible-jump (or collapsed) sampler roam over the discrete structure.

Notebooks

Python Free-knot spline from scratch — reversible-jump birth/death/move over the number and location of knots with coefficients integrated out, plus a PyMC penalized-spline cross-check, on Congdon’s sharp-peak curve (Python) Open → R R cross-check — an mgcv penalized regression spline recovering the same curve (R) Open →

Downloads

freeknot.py The from-scratch reversible-jump free-knot spline engine — g-prior marginal likelihoods and birth/death/move over the knot configuration (new self-contained module) spline_data.csv Congdon's 49-point curve (Y vs X) — a flat baseline with a sharp peak near X = 895 (Congdon Example 3.9)

References

DiMatteo, I., Genovese, C. R. & Kass, R. E. (2001). Bayesian curve-fitting with free-knot splines. Biometrika 88(4), 1055–1071. — the free-knot reversible-jump spline reproduced here
Green, P. J. (1995). Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82(4), 711–732. — the reversible-jump machinery for moving between dimensions
Denison, D. G. T., Mallick, B. K. & Smith, A. F. M. (1998). Automatic Bayesian curve fitting. Journal of the Royal Statistical Society B 60(2), 333–350. — an earlier free-knot curve-fitting sampler
Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g-prior distributions. In Bayesian Inference and Decision Techniques, 233–243. Elsevier. — the g-prior giving the closed-form marginal likelihood
Wood, S. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.). Chapman & Hall/CRC. — the R mgcv penalized-spline cross-check
Congdon, P. (2005). Bayesian Models for Categorical Data, Example 3.9. Wiley. — the 49-point sharp-peak curve and the fixed-knot spline it extends