On Bayesian Curve Fitting Via Auxiliary Variables

In this article we revisit the auxiliary variable method introduced in Smith and kohn (1996) for the fitting of P-th order spline regression models with an unknown number of knot points. We introduce modifications which allow the location of knot points to be random, and we further consider an extension of the method to handle models with non-Gaussian errors. We provide a new algorithm for the MCMC sampling of such models. Simulated data examples are used to compare the performance of our method with existing ones. Finally, we make a connection with some change-point problems, and show how they can be re-parameterised to the variable selection setting.

💡 Research Summary

The paper revisits the auxiliary‑variable framework originally proposed by Smith and Kohn (1996) for Bayesian spline regression, and extends it in two major directions: (1) it treats spline knot locations as random variables rather than fixing them on a pre‑specified grid, and (2) it generalises the error model beyond the Gaussian assumption to accommodate a wide class of non‑Gaussian likelihoods.
The authors introduce a hierarchical model in which each potential knot j is associated with a binary inclusion indicator γj and a continuous location parameter ξj. When γj=1 the knot is active and its position ξj is drawn from a diffuse prior (e.g., Uniform or weakly informative Normal). This construction allows the number of knots and their positions to be explored jointly within a single Markov chain Monte Carlo (MCMC) scheme.
To handle non‑Gaussian errors, the paper embeds the spline model in a generalized linear model (GLM) framework. Using data‑augmentation tricks (e.g., Polya‑Gamma for logistic, latent Gaussian for Poisson), the authors derive closed‑form full‑conditional distributions for the regression coefficients, the spline basis weights, and the hyper‑parameters. Consequently, a Gibbs sampler can be constructed, with Metropolis–Hastings steps only required for the location parameters ξj. The spline degree P and the smoothness penalty λ are given hierarchical priors, yielding a fully Bayesian treatment that automatically balances fit and smoothness.
A novel MCMC algorithm is presented: (i) update γj via a Bernoulli draw that incorporates the marginal likelihood contribution of adding or removing a knot; (ii) propose new ξj for active knots using an adaptive random‑walk proposal tuned to achieve reasonable acceptance rates; (iii) update the spline coefficients conditional on the current knot configuration; (iv) update error‑distribution specific latent variables (e.g., Poisson counts, binomial successes) using standard augmentation schemes. The authors also discuss computational tricks such as block updating of coefficients and efficient recomputation of the design matrix when a knot moves.
Simulation studies compare the proposed method with three benchmarks: (a) the original Smith‑Kohn fixed‑grid approach, (b) a reversible‑jump MCMC (RJMCMC) implementation, and (c) a penalised‑spline (P‑spline) method with cross‑validated knot number. Across scenarios with uniformly spaced knots, clustered knots, and varying signal‑to‑noise ratios, the new method consistently yields lower mean‑squared error (MSE) and higher log‑likelihood values. In non‑Gaussian settings (Poisson and Binomial responses) the auxiliary‑variable scheme retains its accuracy, whereas the Gaussian‑only competitors suffer substantial bias.
The paper further demonstrates that change‑point problems can be re‑parameterised as spline‑knot selection problems. By interpreting each change‑point as a knot of a piecewise‑linear spline, the same auxiliary‑variable machinery can be applied, providing posterior probabilities for the existence and location of change‑points. An empirical example on health‑care cost data illustrates this connection: the model identifies a sharp increase in costs at a specific time point, with credible intervals that naturally quantify uncertainty.
Finally, the authors acknowledge limitations. The number of candidate knots grows with the resolution of the input domain, leading to increased memory and computational demands. They suggest future work on sparsity‑inducing priors (e.g., spike‑and‑slab) and variational Bayes approximations to scale the method. Extensions to multivariate splines, spatially correlated errors, and hierarchical knot structures are also proposed.
In summary, the article delivers a comprehensive Bayesian curve‑fitting framework that simultaneously learns the number and locations of spline knots, accommodates arbitrary error distributions, and offers a unified perspective on change‑point detection. The methodological contributions are validated through extensive simulations and real‑world applications, positioning the approach as a flexible and powerful tool for modern statistical modelling.