Reverse-Engineering Hamiltonian Monte Carlo: The MCMC Engine Behind Modern Bayesian Inference
Every time a Bayesian model gets fit in PyMC, the Hamiltonian Monte Carlo algorithm (HMC) is quietly doing the work behind the scenes of simulating a physical particle moving across a probability landscape to generate a posterior distribution, all without ever needing to compute that distribution directly. This article takes that machine apart. Building on our earlier article covering the foundational Markov Chain Monte Carlo (MCMC) and it’s original variant - the Metropolis-Hastings algorithm, we trace the specific limitations that made HMC necessary then reconstruct it piece by piece. Our tour of HMC starts with the negative log-posterior function that defines the landscape, the gradient function that acts as gravity pulling the particle toward high-probability regions, before concluding with the leapfrog integration settings that govern how far the particle travels before a sample gets recorded. Along the way, we confront the U-turn problem that plagues poorly-tuned HMC and see how the No-U-Turn Sampler (NUTS) solves it automatically. Rather than applying this innovation on toy problems, we do the opposite by utilizing them in a complex Hierarchical Multi-Level Bayesian Regression model designed to estimate wildfire size across British Columbia’s Fire Centre zones to demonstrate how the same skeleton powering a two-parameter demonstration scales, largely unchanged, to a model with dozens of parameters. By the end of it, we’ll find that the black box behind pymc.sample() is a black box to us no longer.
