Reverse-Engineering Hamiltonian Monte Carlo: The MCMC Engine Behind Modern Bayesian Inference

Hamiltonian Monte Carlo (HMC), a variant of Markov Chain Monte Carlo (MCMC), is the computational engine behind many modern Bayesian models, with roots in physics and applications in epidemiology, including COVID-19 intervention modeling in Europe (Flaxman et al., 2020). The article explains how HMC generates posterior distributions, outlines its requirements, and demonstrates sampling on models of varying complexity, building on previous tutorials using PyMC for Bayesian regression on wildfire hectares burnt.

The past few articles we’ve written covered the building blocks of Bayesian Inference and Statistical Machine Learning from the ground up: Starting with Bayes’ Theorem https://pub.towardsai.net/why-most-introductory-examples-of-bayesian-statistics-misrepresent-it-d2e12ac69278 , before learning how to set up models using PyMC https://medium.com/towards-artificial-intelligence/playing-with-fire-and-priors-learning-the-limits-of-bayesian-linear-regression-with-pymc-1897962c25c5 , then reasoning about their causal structures and the data generating processes https://pub.towardsai.net/rethinking-predictors-why-causal-reasoning-matters-in-data-science-part-1-f1d4c1e08068 , and using them to inform our modelling process for real data https://medium.com/towards-artificial-intelligence/putting-dags-to-the-test-what-regression-reveals-about-wildfire-drivers-part-2-c03d4f8a9b13 . With each Bayesian model we built, we’ve had to formally notate our prior knowledge about our parameters’ probability distributions before fitting them to our data. How the model generates a posterior distribution from our priors and likelihoods/data used to feel completely abstract. But now ‘used to’ here is the key phrase that’s a thing of the past. It’s finally clicked. With our previous article exploring the earliest variants of the Markov Chain Monte Carlo https://medium.com/towards-artificial-intelligence/explaining-markov-chain-monte-carlo-using-wildfire-forensics-a334fecaefb3 algorithm, we’ve been able to peel back some of that abstraction. However, the mechanisms that modern probabilistic frameworks like The Hamiltonian Monte Carlo HMC , for some context, is a variant of the Markov Chain Monte Carlo MCMC family of sampling algorithms that’s quietly running as the computational engine behind a significant share of today’s Bayesian models. Although it’s often assumed that MCMC is exclusive to Bayesian Statistics, its roots actually originate from the field of physics where it was originally developed to simulate the behaviour of the hydrogen bomb. MCMC algorithms, particularly HMC, have since had a positive impact on various diverse fields outside of Bayesian Statistics such as in epidemiology where it played a key role in modelling the effects of interventions in Europe during the spread of COVID-19 Flaxman et al., 2020 . In fact, we’ve already utilized HMC without knowing it in all our previous Bayesian Regression articles whenever we generated posterior distributions on total hectares burnt https://medium.com/towards-artificial-intelligence/playing-with-fire-and-priors-learning-the-limits-of-bayesian-linear-regression-with-pymc-1897962c25c5 in a wildfire. In this article, we’ll take a close look as to how HMC conceptually generates a posterior distribution, cover the requirements necessary to get the algorithm running, and then apply these concepts by sampling models of varying complexity. Let’s get started Before we hit the ground running on Hamiltonian Monte Carlo, let’s take a step back and revisit the first iteration and the conceptual foundation of Markov Chain Monte Carlo — the Metropolis-Hastings algorithm. To quickly recap, the goal of the Metropolis-Hastings algorithm is to obtain samples from multi-dimensional probability distributions where sampling them directly is too difficult and complex. You might then ask: How then do MCMC-class algorithms avoid having to compute a highly complex, multi-dimensional distribution directly? Metropolis-Hastings, and other variants, accomplish this by constructing a random walk through parameter space where, over many draws, the visiting frequencies match the shape of the distribution they’re trying to sample from. Reviewing the basics of Metropolis-Hastings gives us the architectural foundation that HMC inherited and highlights how it resolved its past limitations. If you want a more in-depth refresher on this topic, feel free to navigate here https://medium.com/towards-artificial-intelligence/explaining-markov-chain-monte-carlo-using-wildfire-forensics-a334fecaefb3 . Here’s the pseudocode for how the Metropolis-Hastings MCMC algorithm works: 1.Proposal Step : From your current position, propose a direction to move to and sample. The proposal needs to be generated randomly, such as by rolling a die. 2. Acceptance Ratio p move - Note that the earliest iteration of this step in 1953 assumed asymmetric proposalwhere the probability of proposing a move in one direction equaled the probability of proposing the reverse i.e. a fair die . However, Professor W.K. Hastings 1970 improved this step by adding acorrection factor,to account forq,asymmetric proposalswhich are instances where some proposals were structurally more likely than others i.e. a loaded die . This matters in practice whenever a parameter has a natural constraint, such as with a standard deviation where its parameter values must be strictly positive. 3. Accept/Reject Step: If the acceptance ratio is greater than or equal to 1 pmove ≥ 1 , that means the proposed position is more probable than the current one so we can move there with certainty. And if - Ifp, then we can move to the proposed position. But ifmove up, then we must stay put. The lower thatmove≤upis, the more likely the draw fails which therefore means we must stay put.move For all innovations that Metropolis-Hastings brought to the table in terms of sampling from complex probability distributions, it wasn’t able to keep up as data and raw computational power scaled over time. Metropolis-Hastings, along with older variants of MCMC like Gibbs Sampling, hit a scaling wall when sampling from models with hundreds, let alone, thousands of parameters because it could not navigate a model’s sampling space efficiently. This was in part due to the fact that these algorithms had no sense of a posterior’s overall shape. Older variants would essentially explore a sample space blindfolded because their proposal mechanism in step 1 was random without any sense of where the high-probability regions were. On the other hand, Hamiltonian Monte Carlo was able to overcome this limitation and sampled high-dimensional distributions more efficiently by taking advantage of their gradient information to guide its proposal mechanism. To understand how it works, we need to think of the algorithm as a simulation of a “Hamiltonian particle” that behaves similarly to a marble moving through an abstract landscape with varying peaks and valleys. However unlike a marble which will roll continuously before losing energy and eventually stopping at the lowest point due to friction, our Hamiltonian particle doesn’t have any friction to deal with. In this simulation of Hamiltonian mechanics, the particle’s trajectory also behaves like a ball being tossed in the air. During its flight upwards, gravity will slow its momentum before hitting a momentary “pause” where the ball/particle hits zero velocity at the peak of its arc before coming back down to Earth. This is exactly what happens to a Hamiltonian particle moving into a low-probability region of the landscape characterized by upward-facing slopes. The particle slows down the way a ball does climbing toward the top of its arc, and speeds back up moving into a high-probability region, the way a ball accelerates falling back toward the ground. The downward-sloping valleys in the landscape are the high probability regions. However unlike the ball toss, a Hamiltonian particle doesn’t wait for its momentum to permanently settle to zero before we take a sample. Instead, we let the particle travel for a fixed stretch of simulated time, however far the landscape’s shape allows it to travel in that timespan, and wherever it ends up becomes the “sample” we collect for our distribution. A steep uphill section might mean the particle barely gets anywhere while a smooth downhill stretch might carry it a longer way. Either way, the sampling point isn’t chosen by the particle running out of steam, but by us. Now that we’ve covered our simulation of the Hamiltonian particle, let’s zoom out to see how it relates to the generation of a posterior distribution from a Bayesian model and discover how each individual detail we covered plays a role in this process. To build our intuition of how Hamiltonian Monte Carlo samples a Bayesian model that’s been fitted to some data, let’s first generate a dataset with two columns of 50 random values: python import arviz as azimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport pymc as pmimport requestsimport sqlite3import warningsfrom scipy import statsfrom sklearn.preprocessing import StandardScalerwarnings.simplefilter action="ignore", category=FutureWarning pd.set option "mode.chained assignment", None %config Inline.figure format = 'retina'az.style.use "arviz-darkgrid" az.rcParams "stats.hdi prob" = 0.89 sets default credible interval used by arviz print "Generating our test data... \n\n" np.random.seed 42 test datareal = stats.multivariate normal 0, 0 , np.identity 2 x, y = real.rvs 50 .Tprint f"x len = {len x } : {x}, \n\n \y len = {len y } : {y}" Going back to our simulation, the geometric landscape that influences where the Hamiltonian particle is more likely to travel is calculated based on the model we’re sampling from. Notice how Figure 2 illustrates the movement of the particle in a two-dimensional landscape? To put it simply, each parameter in a model actually represents a dimension of that landscape. Although it’s not possible for us to visualize how this process works with a hundred parameter model, we can be confident that the principles for the Hamiltonian particle’s motion remain the same regardless of its landscape’s dimension size. To compute the landscape for our Hamiltonian particle to traverse through, we need to write what’s called a negative log-posterior function . The purpose of the negative log-posterior function is to compute the negative log-probability of the data and its parameters so that the algorithm knows what the “elevation” is for a given set of values at that position. For some that might have been a mouthful so let’s take a step back. Think of every possible combination of parameter values as a point in a landscape, where the number of dimensions equals the number of parameters in your model. With our synthetic dataset, our model has exactly two parameters, μx and Now if we take our generic model from earlier, the negative log-probability function will look like the following where essentially each line in the model is added together. Notice how the first two terms on the left are summed across every data point as per our rule above since these represent the log-likelihood of our data, given the current parameters. On the other hand, the last two appear only once each since there’s just a single μx and And here’s the code form of that equation: The Negative Log-Probabilitydef calc U q, x, y, a=0, b=0.5, k=0, d=0.5 : mu y, mu x = q U = np.sum stats.norm.logpdf y, loc=mu y, scale=1 likelihood + np.sum stats.norm.logpdf x, loc=mu x, scale=1 likelihood + stats.norm.logpdf mu y, loc=a, scale=b prior + stats.norm.logpdf mu x, loc=k, scale=d prior return -U Great, so hopefully we understand how the “landscape” for our Hamiltonian particle is generated based on a simple model. Now, what does this calc U function look like if we’re having to sample from a significantly more complicated one? Let’s test that by bringing back a revamped version of our old Multiple Linear Regression https://medium.com/towards-artificial-intelligence/putting-dags-to-the-test-what-regression-reveals-about-wildfire-drivers-part-2-c03d4f8a9b13 that we used to model wildfire size but instead, we’ll upgrade it to a Hierarchical Multi-Level Regression model DATA EXTRACTION db link = "https://raw.githubusercontent.com/vanislekahuna/wps-labs/main/data/historical bc wildfires/bc wildfires.db"response = requests.get db link with open "bc wildfires.db", "wb" as f: f.write response.content conn = sqlite3.connect 'bc wildfires.db' cursor = conn.cursor join query = """WITH weather AS SELECT wind speed ms, wind direction deg, wind direction, temperature c AS temp c ign date, humidity dewpoint temperature 2m AS humidity dewpoint temp k ign date, soil temperature level 1 AS soil temp lvl1 ign date, fire label AS weather fire label FROM weather data ,terraclimate AS SELECT monthly mean temp c, monthly total precip mm, monthly mean humidity vpd kPa, monthly mean soil moisture mm, fire label AS terraclimate fire label FROM terraclimate bc weather SELECT FROM historical bc wildfiresLEFT JOIN weather ON historical bc wildfires.FIRELABEL = weather.weather fire labelLEFT JOIN terraclimate ON historical bc wildfires.FIRELABEL = terraclimate.terraclimate fire labelWHERE historical bc wildfires.SIZE HA IS NOT NULLAND historical bc wildfires.SIZE HA 0AND historical bc wildfires.FIRELABEL NOT IN '1951-R00037', '1951-R00050', '1951-R00060', '1951-R00067', '1951-R00069', '1951-R00070', '1956-R00107', '1958-V00283', '1956-R00160', '1968-R00088', '1985-V70083', '1985-V70088', '1987-V90016', '1990-V50012' ;"""df = pd.read sql query join query, conn wildfire df = df "FRCNTR", "FIRELABEL", "IGN DATE", "FIRE CAUSE", "SIZE HA", "wind speed ms", "temp c ign date", "humidity dewpoint temp k ign date", "soil temp lvl1 ign date", "monthly mean temp c", "monthly total precip mm", "monthly mean humidity vpd kPa", "monthly mean soil moisture mm" .copy print f"wildfire df shape: {wildfire df.shape} \n" DATA TRANSFORMATION centres = wildfire df "FRCNTR" .unique wildfire df "size log1ptransformed" = np.log1p wildfire df "SIZE HA" wildfire df "wind speed log1ptransformed" = np.log1p wildfire df "wind speed ms" wildfire df "soil temp c lvl1 ign date" = wildfire df "soil temp lvl1 ign date" - 273.15one hot encodding = pd.get dummies wildfire df "FIRE CAUSE" , dtype=int wildfire df = pd.concat wildfire df, one hot encodding , axis=1 large fire df = wildfire df wildfire df 'SIZE HA' = 100 .dropna print f"\n large fire df shape: {large fire df.shape} \n" DATA LOADING/PRE-PROCESSING continuous cols = "wind speed log1ptransformed", "temp c ign date", "monthly total precip mm", "monthly mean humidity vpd kPa" X cont std = StandardScaler .fit transform large fire df continuous cols X mat = np.column stack X cont std :, 0 , large fire df "Lightning" .values, large fire df "Person" .values, X cont std :, 1 , X cont std :, 2 , X cont std :, 3 , zone labels = sorted large fire df "FRCNTR" .astype str .unique zone to idx = {z: i for i, z in enumerate zone labels }zone idx = large fire df "FRCNTR" .astype str .map zone to idx .valuesprint f"Value count of generated zone labels: \n{pd.Series zone idx .value counts .sort index }" y = large fire df "size log1ptransformed" .valuespy HIERARCHICAL MULTI-LEVEL BAYESIAN REGRESSION MODELLING coords = { "zone": zone labels, "predictor": "wind speed", "lightning", "person", "temp", "precip", "vpd" , "obs id": np.arange len y ,}with pm.Model coords=coords as hier model: zone idx = pm.Data "zone idx", zone idx, dims="obs id" X = pm.Data "X mat", X mat, dims= "obs id", "predictor" Population-level grand mean across zones mu alpha = pm.Normal "mu alpha", mu=4.0, sigma=0.5 log 100+1 ≈ 4.6, our size floor mu beta = pm.Normal "mu beta", mu=0.0, sigma=1.0, dims="predictor" Between-zone SD — magnitude of cross-zone heterogeneity per effect sigma alpha = pm.HalfNormal "sigma alpha", sigma=0.5 sigma beta = pm.HalfNormal "sigma beta", sigma=0.5, dims="predictor" Non-centered parameterization z alpha = pm.Normal "z alpha", 0.0, 1.0, dims="zone" z beta = pm.Normal "z beta", 0.0, 1.0, dims= "zone", "predictor" alpha zone = pm.Deterministic "alpha zone", mu alpha + sigma alpha z alpha, dims="zone" beta zone = pm.Deterministic "beta zone", mu beta + sigma beta z beta, dims= "zone", "predictor" mu = alpha zone zone idx + beta zone zone idx X .sum axis=-1 sigma y = pm.HalfNormal "sigma y", sigma=0.5 y obs = pm.Normal "y obs", mu=mu, sigma=sigma y, observed=y, dims="obs id" idata = pm.sample draws=2000, tune=2000, chains=4, target accept=0.95, higher than default — hierarchical funnels need this random seed=42, idata kwargs={"log likelihood": True}, needed for LOO below n div = idata.sample stats "diverging" .sum .item print f"Divergences: {n div}" az.summary idata, var names= "mu alpha", "mu beta", "sigma alpha", "sigma beta", "sigma y" ZONE-LEVEL HETEROGENEITY: Population-level average effect across all zones az.plot forest idata, var names= "mu beta" , combined=True, hdi prob=0.89 plt.axvline 0, color='black', linestyle='--', linewidth=1 plt.title "Zone-level Heterogeneity: Population-level average effect 89% HDI " ZONE-LEVEL HETEROGENEITY: Zone-specific coefficients, partially pooledaz.plot forest idata, var names= "beta zone" , combined=True, r hat=True, ess=True plt.gcf .suptitle "Zone-level Heterogeneity: Zone-specific Coefficients, partially pooled 89% HDI " PREDICTIVE PREFORMANCE: Bayesian R² Gelman et al. 2019 formulation y pred samples = idata.posterior predictive "y obs" .values.reshape -1, len y bayes r2 = az.r2 score y, y pred samples print bayes r2 LOO — The natural comparison point against a non-hierarchical complete-pooling version of this same model to formally test whether the hierarchy earns its added complexityloo = az.loo idata, pointwise=True print loo high k = loo.pareto k.values 0.7 .sum print f"Observations with unreliable LOO estimates: {high k}" With the following equation + code, the important piece here is not to get stuck on the technical complexity of the model itself but instead on the position of each parameter, whether that’s the prior or the likelihood. The take-home point here is that regardless of our model’s complexity, the elevation at any given point in our Hamiltonian particle’s landscape is essentially just a sum of log-probability contributions from every one of its parameters : python def unpack params q, J, P : """ Slice a flat parameter vector q into its named pieces. This mirrors what PyTensor does internally inside PyMC -- every model, no matter how many groups or predictors it has, ultimately reduces to one long vector that the sampler walks through. """ idx = 0 mu alpha = q idx ; idx += 1 mu beta = q idx:idx+P ; idx += P sigma alpha = q idx ; idx += 1 sigma beta = q idx:idx+P ; idx += P z alpha = q idx:idx+J ; idx += J z beta = q idx:idx+J P .reshape J, P ; idx += J P sigma y = q idx ; idx += 1 return mu alpha, mu beta, sigma alpha, sigma beta, z alpha, z beta, sigma ydef hierarchical calc U q, X, y, zone idx, J, P, mu alpha prior= 4.0, 0.5 , mu beta prior= 0.0, 1.0 , sigma alpha prior=0.5, sigma beta prior=0.5, sigma y prior=0.5 : """Negative log-posterior potential energy at a single point q.""" mu alpha, mu beta, sigma alpha, sigma beta, z alpha, z beta, sigma y = \ unpack params q, J, P HalfNormal parameters have zero density outside their support. Real NUTS avoids this entirely by sampling on a log-transformed, constrained scale internally. We stay on the natural scale here. if sigma alpha <= 0 or np.any sigma beta <= 0 or sigma y <= 0: return np.inf Non-centered transform matches the pm.Deterministic lines alpha zone = mu alpha + sigma alpha z alpha beta zone = mu beta + sigma beta z beta Linear predictor mu = alpha zone zone idx + np.sum beta zone zone idx X, axis=-1 log lik = np.sum stats.norm.logpdf y, loc=mu, scale=sigma y log prior = stats.norm.logpdf mu alpha, mu alpha prior + np.sum stats.norm.logpdf mu beta, mu beta prior one sum per predictor - 6 predictors 6 parameters + stats.halfnorm.logpdf sigma alpha, scale=sigma alpha prior + np.sum stats.halfnorm.logpdf sigma beta, scale=sigma beta prior one sum per predictor - 6 predictors 6 parameters + np.sum stats.norm.logpdf z alpha, 0.0, 1.0 one sum per zone - 6 zones 6 parameters + np.sum stats.norm.logpdf z beta, 0.0, 1.0 one sum per zone per predictor - 6 zones x 6 predictors 36 parameters + stats.halfnorm.logpdf sigma y, scale=sigma y prior return - log lik + log prior The second function that the Hamiltonian Monte Carlo algorithm needs is the gradient of the negative log-probability which is a separate slope measurement for each individual parameter in the model. Going back to our Hamiltonian mechanics simulation, the element there that represents the gradient of the negative log probability is the force acting on the particle at its current location. If we’re reusing the ball toss analogy, the gradient here is functioning as the gravity in the toss. For contrast, while the negative log-probability function is calculating the shape of the landscape and determining where the peaks and valleys are, the gradients add value by acting as the “gravity” force in the simulation. However, unlike real-life gravity where its strength remains constant, the gravity of the gradient in the Hamiltonian mechanics simulation varies by location. The pull is strongest in the steepest terrains and essentially non-existent on flat surfaces. Now coming back to our simple model, since we only have two parameters in the model, μx and And the associated Python function for this specific, two-parameter model gradient would be: Gradient of the Negative Log-Probabilitydef calc U gradient q, x, y, a=0, b=0.5, k=0, d=0.5 : mu y, mu x = q G1 = np.sum y - mu y + a - mu y / b 2 dU/dmuy G2 = np.sum x - mu x + k - mu x / d 2 dU/dmux return np.array -G1, -G2 Ok, now that we’ve covered the most complex parts of HMC with the functions, let’s quickly cover the two settings we need to tune in order to get the algorithm started. The first is the leapfrog steps L which is simply the count of how many steps our Hamiltonian particle takes to complete one trajectory. If you remember earlier in the ball toss analogy, we mentioned how we specify the span of time that the ball spends in the air? Together, these two settings determine how long our particle’s trajectory runs before we stop and record a sample. The second setting we need to tune is the step size ϵ which is the distance the simulation covers in a single leapfrog set. For example with our simple model notation from earlier, we specified that: “The number of leapfrog steps to L = 11 and the step size to An issue that may arise from these two settings is that if we specify a poor combination of the leapfrog steps and the step size then we can run into an issue known as the U-turn problem. The U-turn problem arises when a Hamiltonian particle’s simulated trajectory travels far enough to loop back to its own starting point or the general area it started from. As a consequence of this issue, we’ll end up with consecutive samples that are highly autocorrelated and will require far more computation for the particle to explore its landscape. Luckily, Hoffman and Gelman 2014 solved this problem over a decade ago now with their proposal of the No-U-Turn Sampler NUTS . In essence, what the No-U-Turn Sampler NUTS does is automatically detect when a trajectory starts to double back on itself and stops the simulation right before that happens, instead of relying on a fixed, manually-tuned number of leapfrog steps L or step size The last requirement for HMC is the starting point which we’ll arbitrarily assign to q = −0.1, 0.2 . The starting point will kick off our “Hamiltonian particle” simulation which allows us to sample from a posterior distribution based on which areas the particle visits most frequently. Once we plug in our two functions and two settings, with the flexibility to plug in whichever model’s functions we’re sampling from, then our entire Hamiltonian Monte Carlo sampler will look like the following: python def Hamiltonian Monte Carlo U, grad U, epsilon, L, current q, args : q = current q.copy p = np.random.normal loc=0, scale=1, size=len q random flick - p is momentum current p = p.copy Make a half step for momentum at the beginning p -= epsilon grad U q, args / 2 initialize bookkeeping - saves trajectory qtraj = np.full L + 1, len q , np.nan ptraj = qtraj.copy qtraj 0, : = current q ptraj 0, : = p Code 9.9 starts here Alternate full steps for position and momentum for i in range L : q += epsilon p Full step for the position qtraj i + 1, : = q Make a full step for the momentum, except at the end of trajectory if i = L - 1: p -= epsilon grad U q, args ptraj i + 1, : = p Make a half step for momentum at the end p -= epsilon grad U q, args / 2 ptraj L, : = p Negate momentum at end of trajectory to make the proposal symmetric p = -1 Evaluate potential and kinetic energies at start and end of trajectory current U = U current q, args current K = np.sum current p 2 / 2 proposed U = U q, args proposed K = np.sum p 2 / 2 Accept or reject the state at end of trajectory, returning either the position at the end of the trajectory or the initial position accept = False if np.random.uniform < np.exp current U - proposed U + current K - proposed K : new q = q accept accept = True else: new q = current q reject return dict q=new q, traj=qtraj, ptraj=ptraj, accept=accept And often one of the results we like to see immediately after running the sampler is a trace plot to evaluate its preformance. With a healthy trace plot, we generally want to see two things. The first thing to look for within each individual chain is that the particle should fluctuate rapidly and densely around a stable value, not drift steadily upward or downward over time. The second thing to look for across all four chains is that all four chains should converge on that same stable value. Chains that settle into visibly different regions signal a problem, even if each one looks stable on its own. The trace plot below is an example of a high preforming chain that was generated for the complex Hierarchical Bayesian Regression from earlier. az.plot trace idata, var names= "mu alpha", "mu beta", "sigma beta" , compact=True plt.gcf .suptitle "Convergence Diagnostic Trace Plot" plt.tight layout We set out to reverse-engineer the machine behind modern Bayesian inference and by now, we’ve succeeded in taking it apart Let’s recap what we accomplished. What began as an abstract black box that somehow turned priors and likelihoods into a posterior distribution is a black box no longer. We then traced the lineage from Metropolis-Hastings’ blind, random walk through parameter space, to the moment gradient information transformed that walk into something far more deliberate: A Hamiltonian particle, gliding through a landscape shaped by our own model’s mathematics, accelerating into high-probability valleys and decelerating up low-probability slopes. Lastly, we built the two functions that define that landscape and the forces acting on it, tuned the two settings that govern how far our particle travels, and watched all five pieces come together into a single, working sampler. What makes this worth understanding isn’t just the mechanics themselves: It’s the foundation with which we’re using Bayesian Inference to drive real-world impact. Every time pymc.sample https://www.pymc.io/projects/docs/en/stable/api/generated/pymc.sample.html pymc.sample runs, this sophisticated, and highly robust, simulation of a “Hamiltonian particle” is what's quietly happening underneath it all. Hamiltonian Monte Carlo isn't a niche statistical curiosity but rather, the quiet computational engine behind a significant share of the inference done in Bayesian statistics today. For example, we noted in the intro how it played a direct role in modeling the effects of interventions across Europe during COVID-19 which helped shape real public health decisions during a genuine crisis Flaxman et al., 2020 . We even saw its reach firsthand, too. The same skeleton that generated a posterior distribution for a toy two-parameter model was able to scale, without needing to be reinvented, to a complex Hierarchical Bayesian Regression model spanning dozens of parameters across BC’s Fire Centre zones. With our HMC quietly working in the background, our Hierarchical Bayesian Regression model was able to estimate how wind speed, lightning, precipitation, and a handful of other atmospheric variables could shape the size of a BC wildfire, with the uncertainty in that estimate quantified honestly rather than papered over. Our toy example of HMC scaled to a complex machine running on real data, producing a real answer to a real question We came into this article treating the posterior as something that simply appeared. Now we can now proudly say we’re leaving it with a working understanding of exactly how it gets there though one leapfrog step, one gradient, and one accepted proposal at a time. Thank you for reading through to the very end. This was an amazing journey I’m glad we were able to travel through together. Reverse-Engineering Hamiltonian Monte Carlo: The MCMC Engine Behind Modern Bayesian Inference https://pub.towardsai.net/reverse-engineering-hamiltonian-monte-carlo-the-mcmc-engine-behind-modern-bayesian-inference-e1d6b54a8c79 was originally published in Towards AI https://pub.towardsai.net on Medium, where people are continuing the conversation by highlighting and responding to this story.