PyMC Sampling Algorithms for Bayesian Inference
PyMC is a leading Python library for probabilistic programming and Bayesian statistical modeling that provides an array of sampling algorithms to approximate complex posterior distributions. For continuous parameters, it predominantly uses gradient-based Markov Chain Monte Carlo (MCMC) algorithms like the No-U-Turn Sampler (NUTS), while also offering classical algorithms like Metropolis-Hastings and Slice sampling for non-differentiable or discrete spaces. Additionally, PyMC incorporates advanced population methods like Sequential Monte Carlo (SMC) to handle multimodal geometries and facilitate model comparison.
No-U-Turn Sampler (NUTS)
NUTS is the default MCMC engine in PyMC for continuous random variables. It is an extension of Hamiltonian Monte Carlo (HMC) that simulates physical dynamics along the posterior surface using gradient information provided by PyMC's computational backend (PyTensor).
Unlike standard HMC, which requires manual tuning of trajectory lengths and step sizes, NUTS automatically tunes the step size during the warm-up phase and determines when the trajectory begins to turn back on itself. This prevents wasted computation, eliminates random walk behavior, and allows efficient exploration of high-dimensional parameter spaces with complex correlations.
Metropolis-Hastings
The pm.Metropolis step method implements the classical
random-walk Metropolis-Hastings algorithm. It proposes transitions based
on a Gaussian distribution centered at the current state and accepts or
rejects them according to the Metropolis acceptance probability.
Because it does not rely on gradient calculations, Metropolis is computationally lightweight per step, making it a reliable fallback for models with non-differentiable likelihoods or when gradients are unavailable. However, it suffers from slow mixing and random-walk behavior in high dimensions.
Binary and Categorical Samplers
For discrete parameters where gradient-based methods cannot operate, PyMC provides specialized discrete step methods:
- BinaryMetropolis: Tailored specifically for binary variables (Bernoulli or deterministic boolean switches), proposing flips between 0 and 1.
- CategoricalGibbsMetropolis: Designed for discrete variables with categorical outcomes, updating states by cycling through potential categories based on conditional probabilities.
Slice Sampling
The pm.Slice algorithm draws samples uniformly from the
region under the probability density function. It adaptively defines an
interval around the current state and shrinks it until an acceptable
point is sampled. Slice sampling is more robust to step-size tuning than
Metropolis-Hastings and mixes better across flat or constrained targets,
though it requires multiple likelihood evaluations per step.
Sequential Monte Carlo (SMC)
Sequential Monte Carlo (pm.sample_smc) is a
particle-based sampling algorithm fundamentally different from standard
MCMC chains. It maintains a population of particles that gradually
transitions from the prior distribution to the posterior distribution
through a series of intermediate stages.
SMC is particularly effective for:
- Multimodal distributions: Overcoming energy barriers that trap single MCMC chains.
- Marginal likelihood estimation: Directly computing model evidence for Bayesian model selection.
- Approximate Bayesian Computation (SMC-ABC): Inferring parameters for simulator-based models where the likelihood function is intractable.
Differential Evolution Metropolis (DEMetropolis)
PyMC includes pm.DEMetropolis (and its variations, such
as DEMetropolisZ), which is a population MCMC method. It
runs multiple chains in parallel and uses vector differences between
pairs of chains to generate proposals. This approach allows the sampler
to efficiently learn the orientation and scale of correlated parameters
without requiring analytical gradients.
Compound Step Methods
In models containing both continuous and discrete parameters, PyMC uses a compound step method. It automatically inspects the model graph and assigns the optimal sampler to each variable subset—typically assigning NUTS to all continuous variables and Metropolis, BinaryMetropolis, or CategoricalGibbsMetropolis to discrete ones—executing them sequentially within each iteration.