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.

Read More

Explaining Markov Chain Monte Carlo using Wildfire Forensics

A wildfire has burned overnight in a remote stretch of forest, and Prof. Markov has been called in to find out where it started. He has a burn map divided into a 5-by-5 grid (25 total cells), some historical records on lightning strikes and logging roads, and a single die with eight sides. That’s it. No coordinates, no shortcuts, no way to survey the entire area all at once. However, what he does have is a procedure: At every cell, ask whether the next one over is more or less likely than where he’s standing, and let that comparison do the rest after running for a hundred thousand iterations. By the end, the places he visited most often are the answers to his original question of where the fire started, down to a full probability distribution rather than a single guess.

Read More

Putting DAGs to the Test: What Regression Reveals about Wildfire Drivers (Part 2)

This article puts the hypothesized causal DAG from Part 1 to the test by fitting various Multiple Linear Regression models to ~4,000 large BC wildfires (>100 hectares) and examining whether the atmospheric relationships we drew actually hold up in the data. Before the main analysis, we take a pedagogical detour through Bayesian Multiple Linear Regression using PyMC to demonstrate what happens when weakly informative priors meet a large dataset and why that finding justifies switching to simpler Frequentist methods for the rest of the analysis. The honest result: Our final model explains roughly 11% of variance in fire size, which sounds underwhelming until you consider what that 11% actually represents. Two of our hypothesized mediators don’t survive contact with the data, and our revised DAG is far better for it.

Read More

Rethinking Predictors: Why Causal Reasoning Matters in Data Science (Part 1)

This article makes the case for causal inference as a framework for thinking more deliberately about the models we build and the predictors we choose. Using multiple linear regression as our primary analytical lens, we introduce the core ideas behind causal reasoning, such as the data generating processes, mediators, confounding variables, and spurious correlations, and apply them to a real dataset of large BC wildfires (>100 hectares). The analysis culminates in a hypothesized causal DAG mapping the atmospheric conditions we believe drive wildfire size, which we’ll formally test with regression in Part 2. Spoiler: not every arrow in our DAG will survive contact with the data, and that’s precisely the point.

Read More

Playing With Fire and Priors: Learning the Limits of Bayesian Linear Regression with PyMC

This notebook demonstrates the complete workflow for building a Bayesian Simple Linear Regression model using PyMC to predict wildfire sizes solely from wind speed data. In it we formalize the model using statistical notation, implement prior predictive simulations to validate our assumptions, and then generate posterior distributions through Markov Chain Monte Carlo (MCMC) sampling. Unfortunately, the analysis revealed that wind speed alone is a weak (quite honestly, a terrible) predictor of fire size. Nonetheless, the value-add in our work was that it provided practical insights about the importance of model diagnostics and the pitfalls of violating assumptions of linearity and heteroscedasticity when working with real-world data!

Read More

An Intro to PyMC and the Language for Describing Statistical Models

This notebook introduces the language of statistical model notation used to describe Bayesian models, demonstrating how mathematical notation translates directly into executable PyMC code. Using a real-world example of Vancouver Island Coastal Wolves (a recently discovered subspecies known for displaying distinguished behaviour from other wolf populations such the ability to swim great distances or the significant prevelance of marine organisms in their diet), we build a Bayesian model to estimate gender ratios from limited sample data to emphasize how posterior distributions reveal the full range of plausible outcomes rather than single point estimates.

Read More

Why Most Introductory Examples of Bayesian Statistics Misrepresent It

This notebook introduces the basic idea behind Bayes’ Theorem and highlights some key differences between Bayesian Statistics and the widely taught traditional branch of statistics, commonly known as Frequentist Statistics. Furthermore we challenge the traditional medical testing example used in countless textbooks to introduce Bayesian Inference, arguing that using fixed constants (also known as point estimates) misrepresents the true nature of Bayesian Statistics. Rather than inputting single values in Bayes’ Theorem, we utilize a more faithful Bayesian approach by considering probability distributions of possible outcomes, thus revealing how disease prevalence uncertainty affects diagnostic accuracy.

Read More

How to Augment Wildfire Datasets with Historical Weather Data using Python and Google Earth Engine

This notebook demonstrates how to augment wildfire datasets with historical weather data using Python and Google Earth Engine’s ERA5 dataset. Here we transform basic fire records (coordinates and timestamps) into enriched datasets containing comprehensive environmental context including temperature, wind patterns, humidity, and soil conditions which are all critical information for fire risk modeling and analysis. The workflow includes automatic retry logic for handling API timeouts and incremental batch saving to prevent data loss during long processing runs.

Read More

You're up and running!

Next you can update your site name, avatar and other options using the _config.yml file in the root of your repository (shown below).

Read More