Computational Inference for Bayesian Models: MCMC, HMC, NUTS, and Variational Inference Without the Hand-Waving

Part 2 of 7. The machinery that turns a generative model into a posterior.

This is the second article of a seven-part series on Bayesian methods in industry. Part 1 covered the foundations of Bayesian thinking. Part 3 covers hierarchical models and GLMs in product systems. Part 4 covers Bayesian Marketing Mix Modeling end to end. Part 5 covers classical causal inference and debiased machine learning. Part 6 covers heterogeneous effects, uplift, and bandits. Part 7 covers Bayesian Deep Learning, Gaussian Processes, and Bayesian Optimization.

The Thesis

Most Bayesian disasters in industry are not modeling failures. They are inference failures. The model was correct, the priors were defensible, the likelihood was appropriate, and the practitioner concluded the model was wrong because the sampler produced nonsense. The nonsense was not the model. It was a divergence-riddled NUTS run on an unidentifiable parameterization, or a variational approximation that confidently underestimated the posterior variance, or a Gibbs sampler that mixed so poorly that the effective sample size was thirty after fifty thousand iterations.

Bayesian inference often reduces to approximating expectations under a posterior distribution that is analytically intractable. The choice of method (sampling, optimization, message passing, quadrature, or hybrids) determines what you can know about that distribution within a realistic compute budget.

The literature on Bayesian inference algorithms is enormous. Most of it is irrelevant for the practitioner who needs to ship a model on Monday. The set of algorithms most practitioners routinely rely on is relatively small: NUTS in its modern implementations (Nutpie, NumPyro, Stan, BlackJAX), variational inference in two flavors (ADVI and Pathfinder), conjugate updates where they are available, and INLA for the specific class of latent Gaussian models. Other families (Gibbs, SMC, EP, particle methods, ensemble samplers, annealed importance sampling, pseudo-marginal MCMC, stochastic-gradient MCMC) matter and are covered in later parts of this article when the problem calls for them. Knowing what each tool is, what assumptions it makes, when it fails, and how to diagnose the failure is the operational skill.

This article does not derive the Hamiltonian dynamics from first principles. It assumes that anyone serious about Bayesian inference has read the relevant chapters of Gelman or the Betancourt case studies and knows the math. What it covers is the operational vocabulary: when NUTS gives you divergences, what they mean, what to do about them, when to abandon NUTS for VI, how to tell that VI is lying to you, when conjugate updates beat anything fancier, and which combination of diagnostic statistics is sufficient to trust a posterior.

Part I. Why the Posterior Is Rarely Tractable

The posterior p(θ | y) ∝ p(y | θ) p(θ) is a probability density up to the normalizing constant p(y) = ∫ p(y | θ) p(θ) dθ. For most models in industry, this integral is high-dimensional and analytically intractable. The dimension grows with the number of parameters; a hierarchical model with random intercepts and slopes for a hundred regions has at least two hundred parameters before the global level, and the integration is over the joint distribution of all of them.

Closed-form posteriors exist for a small set of models with conjugate priors. Beta-Bernoulli, Gamma-Poisson, Normal-Normal with known variance, and a few others. For everything else, we approximate. The two dominant approaches are sample-based (MCMC and its modern variants) and density-based (variational inference and Laplace approximation).

The sample-based approach produces a set of draws {θ{(s)}}_{s=1}{S} such that any function of the posterior can be approximated by the corresponding Monte Carlo average. The posterior mean is approximated by the sample mean, the posterior variance by the sample variance, the posterior probability of an event by the proportion of samples in the event, the posterior predictive distribution by pushing samples through the likelihood. This is general-purpose: once you have samples, anything you want to compute about the posterior is a one-liner.

The density-based approach approximates the posterior with a tractable parametric family (typically Gaussian or factorized Gaussian) and finds the member of the family closest in some sense (typically KL divergence) to the true posterior. The result is a parametric description of the posterior, which is faster to evaluate but only as good as the assumption that the true posterior is in or near the family.

The tradeoff is the central engineering decision in computational Bayesian inference. Sample-based methods are consistent under appropriate conditions (ergodicity, sufficient mixing, correct integrator) and typically slow; in finite time they can still be biased by adaptation, slow mixing, multimodality, or numerical discretization. Density-based methods are fast and biased by construction, the bias being whatever the chosen family cannot represent. Neither is universally better. The practitioner’s job is to know which is appropriate for which problem.

Part II. The MCMC Family

Markov chain theory in one page

The core idea of MCMC is to construct a Markov chain whose stationary distribution is the posterior of interest. A Markov chain is a sequence of random variables where the distribution of the next state depends only on the current state, not on the history. Under conditions of irreducibility, aperiodicity, and recurrence, the chain has a unique stationary distribution, and ergodicity guarantees that long-run averages of functions of the chain converge to expectations under that stationary distribution.

To make this useful for sampling from a target distribution π(θ), you need a transition kernel q(θ’ | θ) that has π as its stationary distribution. The Metropolis-Hastings construction provides a generic recipe: propose θ’ from any proposal distribution g(θ’ | θ), accept with probability

and otherwise stay at θ. The detailed balance condition π(θ) q(θ’ | θ) = π(θ’) q(θ | θ’) is automatically satisfied by this construction, which guarantees that π is stationary. The catch is that detailed balance does not by itself produce a chain that mixes well. A poorly chosen proposal can produce a chain that takes millions of steps to traverse the posterior.

The practical question is therefore not whether MCMC is correct in theory (it is), but whether the specific MCMC algorithm you are running explores the posterior efficiently within your compute budget. Different algorithms make different proposals and have different properties.

Random-walk Metropolis

The simplest MCMC algorithm proposes θ’ from a symmetric distribution centered at θ, typically a Gaussian. The acceptance probability simplifies because the proposal densities cancel.

Random-walk Metropolis was the workhorse of Bayesian statistics from the early 1990s to roughly 2010. It is simple, requires only the ability to evaluate the unnormalized posterior, and works on essentially any model. It also scales catastrophically with dimension. The optimal acceptance rate for random-walk Metropolis in high dimensions is around 23 percent, and at that acceptance rate the effective sample size per second falls off as roughly d^{-1} where d is the dimension. A two-hundred-dimensional posterior, which is modest for a hierarchical model, requires hundreds of thousands of iterations to produce a few hundred effective samples.

The practical consequence is that random-walk Metropolis is no longer the right tool for any model with more than a few dimensions, except as a final-stage refinement or in places where gradient information is unavailable.

Gibbs sampling

Gibbs sampling exploits a special structure: if you can sample from the full conditional distributions p(θi | θ{-i}, y) exactly, you can construct a Markov chain that updates one parameter (or block of parameters) at a time. There is no rejection step; every proposal is accepted because it is drawn from the exact full conditional.

Gibbs is the natural fit for conjugate models, where the full conditionals are themselves tractable distributions. It is also the standard inside more complex algorithms; many Bayesian methods use Gibbs as a building block.

The weakness of Gibbs is that it mixes poorly when parameters are highly correlated. The full conditionals only allow movement along the coordinate axes, so a posterior with a strong off-diagonal ridge requires many small Gibbs steps to traverse. This is one of the main reasons HMC became dominant: it can propose moves that are not aligned with the coordinate axes.

Metropolis-Adjusted Langevin Algorithm

MALA uses the gradient of the log posterior to construct a more informed proposal:

This biases the proposal toward higher-density regions of the posterior. The resulting chain mixes faster than random-walk Metropolis in high dimensions, and the gradient information is typically cheap to obtain through automatic differentiation in modern probabilistic programming frameworks.

MALA was the bridge from gradient-free MCMC to fully gradient-based methods. In practice it has been largely superseded by HMC and NUTS, which use the same gradient information more cleverly.

Part III. Hamiltonian Monte Carlo and NUTS

The physical intuition

Hamiltonian Monte Carlo treats the negative log posterior as a potential energy function. A fictitious momentum variable r is introduced, with its own Gaussian distribution, and the joint distribution p(θ, r) ∝ p(θ | y) N(r | 0, M) has the target posterior as its marginal in θ. The dynamics of the joint distribution under the Hamiltonian

conserve H along the trajectory. Integrating these dynamics for some time produces a proposed move in (θ, r) space, which is then accepted or rejected based on the small numerical error in H introduced by the integrator.

The crucial property of HMC is that the trajectory can move long distances through the posterior while keeping the acceptance probability high, because the dynamics are conservative. This is the source of HMC’s much better scaling with dimension compared to random-walk methods. Under favorable regularity conditions, asymptotic analyses (Beskos et al. 2013) suggest HMC scales closer to d^{1/4} rather than the d^{-1} of random-walk Metropolis. In practice this scaling argument is a useful intuition rather than a guarantee: geometry, curvature, conditioning, and the structure of the posterior dominate any clean asymptotic rate.

The leapfrog integrator is the standard numerical integrator for HMC. It alternates updates of position and momentum in a specific symplectic pattern that preserves volume in phase space. The two hyperparameters of HMC are the step size ε (controlling the granularity of the integration) and the number of leapfrog steps L (controlling the trajectory length). Tuning these by hand was the major operational obstacle to using HMC for two decades.

The No-U-Turn Sampler

Random-walk Metropolis vs HMC/NUTS on a correlated 2D posterior

NUTS (Hoffman and Gelman 2014) is the algorithm that made HMC usable without manual tuning. The key idea is to automatically choose the trajectory length by detecting when the trajectory starts to curve back on itself, at which point further integration is wasted. NUTS uses a recursive tree-doubling procedure that produces trajectories of variable length and balances integration error against trajectory exploration.

The hyperparameters that NUTS does still expose, and that occasionally need tuning, are:

• target_accept, the target acceptance probability during the adaptive warmup phase. The default is 0.8 in PyMC and Stan, and the recommendation when divergences appear is to increase it to 0.95 or 0.99. Higher target_accept means smaller step sizes, which is the right response when the geometry of the posterior is challenging.
• max_treedepth, the maximum depth of the tree-doubling. The default of 10 corresponds to up to 1024 leapfrog steps per iteration. When the model frequently hits max_treedepth, NUTS is failing to find a U-turn within its budget; the cause is usually a very long, thin posterior, and the fix is either reparameterization or, less ideally, a larger max_treedepth.
• step_size, which is adapted during warmup to achieve the target acceptance. You almost never set this manually after warmup.
• mass_matrix, which scales the momentum. The default is a diagonal mass matrix adapted from posterior variances during warmup. For posteriors with strong correlations, a dense mass matrix can help, though at higher computational cost per step.

The modern NUTS implementations have converged on essentially the same algorithm with minor variations.

Stan’s NUTS is the original and most mature implementation. The C++ backend and the long history of stress-testing make it the most trustworthy choice when reproducibility matters.

PyMC’s NUTS (currently defaulting to Nutpie in PyMC 5 and above) is a competitive implementation written in Rust. Nutpie is fast, supports JAX, and is the default for new PyMC installations.

NumPyro’s NUTS runs on JAX and scales naturally to GPUs. It is the right choice for very large models where the gradient computation can use accelerators.

BlackJAX is a JAX-native implementation that is the building block for many research papers and increasingly used in production JAX pipelines.

The choice between these is mostly a matter of which ecosystem you live in. The samplers should produce statistically similar posterior estimates when the model is well-behaved and the chains have converged, though the specific draws will differ because of adaptation differences, numerical precision, warmup behavior, mass-matrix support (diagonal versus dense), JAX-vs-C++ floating-point differences, and the exact handling of edge cases such as divergent transitions.

Part IV. Diagnostics

A converged posterior is not a guaranteed correct posterior. It is a posterior that the sampler successfully drew samples from. The diagnostics in this section check that the sampler converged. None of them check that the model is correct. That is a separate question covered by posterior predictive checks, which were discussed in the previous article.

R-hat

NUTS convergence diagnostics: good vs bad traceplots

R-hat (the potential scale reduction factor, originally Gelman and Rubin 1992, refined in Vehtari, Gelman, Simpson, Carpenter, Bürkner 2021) compares the variance within each chain to the variance between chains. If the chains have converged to the same stationary distribution, the two variances should be approximately equal, and R-hat should be close to 1. The modern recommendation is that R-hat should be below 1.01 for every parameter you care about.

The rank-normalized split-R-hat (the version in Vehtari et al. 2021) is the right diagnostic to use. It is robust to non-Gaussian posteriors and to chains that appear to have converged but are stuck in different modes. Older versions of R-hat can give false positives of convergence in those cases.

The crucial caveat is that R-hat checks whether the chains agree with each other. If all chains are stuck in the same local mode, R-hat will look good and the posterior will be wrong. This is rare in practice for NUTS on well-specified models, but it can happen for mixture models, models with label-switching pathologies, and bimodal posteriors. Looking at the trace plot and the pair plot is a second line of defense.

Effective sample size

The effective sample size (ESS) measures how many independent samples your correlated MCMC chain is worth. The relationship is

where n is the number of draws and ρ_k is the autocorrelation at lag k. A chain that produces near-independent samples has ESS close to n. A chain with strong autocorrelation has ESS much smaller than n.

The modern recommendation is to report two flavors: bulk-ESS and tail-ESS. Bulk-ESS measures how well the chain explores the center of the posterior. Tail-ESS measures how well it explores the tails. The minimum acceptable ESS depends on what you want to estimate. For posterior means and central credible intervals, bulk-ESS above 400 per chain is usually sufficient. For extreme quantiles (95th percentile and beyond), tail-ESS needs to be substantially larger.

If ESS is low, the chain is mixing poorly. The diagnostic does not tell you why; it just tells you that the sampler is not exploring efficiently. The investigation moves on to the trace plot, the pair plot, and the energy diagnostic.

Divergences

A divergence in NUTS happens when the leapfrog integrator produces a large error in the conserved Hamiltonian during a trajectory. NUTS reports the number of divergent transitions per chain. The interpretation is that the posterior has a region of high curvature that the integrator cannot follow accurately.

Divergences are the most informative diagnostic when they appear, because they tell you where the trouble is. Plot the divergent samples against the rest of the posterior, typically as a pair plot with divergent points highlighted, and you will usually see them clustered in a specific region. That region is the geometric pathology.

The two canonical pathologies are funnels and ridges.

A funnel appears in hierarchical models when the group-level variance parameter τ is small. As τ shrinks, the conditional distribution of the group-level effects θ_g becomes very tight, and the geometry of the joint posterior becomes a narrow tunnel. The leapfrog integrator with a fixed step size cannot adapt locally to the changing scale, and it produces large errors at the tip of the funnel. The fix is reparameterization, which is the subject of the next section.

A ridge appears when two parameters are highly correlated. The posterior is concentrated along a line in the joint space, and the integrator must take very small steps perpendicular to the ridge while taking longer steps along it. The fix is either reparameterization to decorrelate the parameters or a dense mass matrix that captures the correlation.

When you have divergences, look at the pair plot first. The pattern of divergences usually identifies the pathology immediately.

Energy diagnostic

The energy diagnostic E-BFMI (Energy Bayesian Fraction of Missing Information) measures how well the momentum resampling at each iteration explores the energy distribution of the joint posterior. Low E-BFMI (below 0.3 is concerning) indicates that the momentum resampling is failing to give the trajectory enough energy to escape its current region of the posterior. The typical cause is a heavy-tailed posterior; the fix is reparameterization or a different parameterization of the heavy-tailed component.

Simulation-based calibration

Simulation-based calibration (Talts, Betancourt, Simpson, Vehtari, Gelman 2018) is one of the strongest end-to-end checks that the inference pipeline (model plus sampler) is jointly self-consistent. It validates the inference machinery on the model you wrote down; it does not, by itself, validate that the model is an adequate description of the world (that is the job of prior and posterior predictive checks). The idea is to draw parameter values from the prior, generate fake data conditional on those parameters, run the inference on the fake data, and check that the rank of the true parameter in the resulting posterior is uniformly distributed.

SBC is computationally expensive (you need to run inference hundreds of times) and is therefore typically used during model development, not during production runs. When it works, it provides strong evidence that the model and sampler are jointly correct. When it fails, it tells you something is wrong but does not always tell you what. The technique is the right thing to run when you are developing a new model class or when something has gone subtly wrong with an established model.

Part V. Reparameterization

Reparameterization is the most underused tool in the practitioner’s kit. It is also the single technique that most often turns a model that NUTS cannot fit into a model that NUTS fits in minutes.

The centered and non-centered parameterizations

Neal’s funnel: the centered parameterization (left) is pathological; the non-centered version (right) is NUTS-friendly

Consider a hierarchical model where the group-level effects θ_g have prior

This is the centered parameterization. It is intuitive: each group has its own effect drawn from a population distribution. The pathology is that when τ is small (groups are similar), the conditional posterior of θ_g given τ has scale proportional to τ, and the joint posterior has the funnel geometry described above.

The non-centered parameterization rewrites this as

Mathematically identical. Geometrically very different. The joint posterior of (z_g, τ) has no funnel; the z_g values are decoupled from τ by construction. NUTS samples this efficiently.

The non-centered parameterization is the default recommendation for hierarchical models with sparse or moderate within-group data. The centered parameterization is preferred only when there is so much within-group data that the conditional posterior of θ_g is well-determined regardless of τ; in that regime the non-centered version can actually be slower.

In PyMC the non-centered parameterization is sometimes available as a centered=False argument; otherwise you write it explicitly as a deterministic transformation of standard normal variables. In Stan it is written directly.

Other reparameterizations

Several other reparameterizations are part of the standard kit.

Log-transformation of positive parameters. A standard deviation or a positive scale parameter should be sampled on the log scale. NUTS works on the log scale natively for many distributions in PyMC and Stan; you do not need to do this manually.

Cholesky parameterization of covariance matrices. Sampling a full covariance matrix directly is hard. Sampling a Cholesky factor and reconstructing the matrix is much easier. The LKJ prior on correlation matrices is usually combined with Cholesky parameterization for efficiency.

Whitening of correlated parameters. If two parameters have a known correlation structure, sampling in the whitened (uncorrelated) coordinates and transforming back at the end is much more efficient than sampling in the correlated coordinates directly.

Concentration parameterization of Dirichlet. When some elements of a Dirichlet distribution are expected to be near zero, sampling on a stick-breaking representation or on log-scale concentrations is more efficient than the standard parameterization.

The general principle is that the parameterization of the model and the parameterization given to the sampler are different concerns. You should write the model in whatever form is most natural for thinking about the science, and reparameterize for the sampler. Modern probabilistic programming languages let you do both transparently.

Part VI. Variational Inference

When MCMC is too slow, the alternative is variational inference. The basic idea is to find a tractable distribution q(θ | φ) that approximates the true posterior p(θ | y) by minimizing the KL divergence

Minimizing this KL divergence is equivalent to maximizing the evidence lower bound (ELBO):

The first term is the expected log joint density, which can be approximated by Monte Carlo. The second term is the entropy of the variational distribution, which has a closed form for standard parametric families. Modern VI uses stochastic gradient descent on the ELBO with the reparameterization trick to compute unbiased gradients.

ADVI

Automatic Differentiation Variational Inference (Kucukelbir, Tran, Ranganath, Gelman, Blei 2017) is the algorithm that made VI broadly usable. ADVI transforms the parameter space to be unconstrained, fits a Gaussian (or factorized Gaussian) variational distribution in the transformed space, and uses automatic differentiation to compute the ELBO gradient.

ADVI is implemented in PyMC, Stan, and NumPyro. It is the default fast alternative to NUTS for any model expressible in those frameworks. It produces a result in seconds to minutes for problems that NUTS would take hours on.

The two flavors are mean-field (factorized Gaussian over all parameters) and full-rank (joint Gaussian with a full covariance matrix). Mean-field is faster but cannot capture correlations between parameters; full-rank can, but has many more variational parameters and is more expensive to optimize.

The crucial caveat about ADVI is that it systematically underestimates posterior variance. The KL divergence KL(q | p) penalizes putting mass where the true posterior has none, which produces concentrated approximations. The opposite KL divergence would have the opposite bias, but it is not what ADVI optimizes. The practical consequence is that ADVI credible intervals are too narrow, often by a factor of two or more for parameters whose posteriors have heavy tails.

For prototyping and exploration, this is acceptable. For reporting uncertainty to a stakeholder, it is not. The standard workflow is to use ADVI to iterate on the model and to use NUTS to produce the final result.

Pathfinder

Pathfinder (Zhang, Carpenter, Gelman, Vehtari 2021) is a more recent variational method that uses quasi-Newton optimization to construct a sequence of Gaussian approximations along the optimization path. It is dramatically faster than ADVI for some problems and produces useful starting points for NUTS.

The most common use of Pathfinder in production is as a NUTS initialization. Rather than starting NUTS from random points (which wastes the warmup phase getting to the high-density region), Pathfinder gives you a Gaussian approximation that NUTS can refine. The combined runtime is often less than either method alone.

Normalizing flows and other modern VI

Mean-field and full-rank Gaussians are not the only variational families. Normalizing flows construct flexible distributions by transforming a simple base distribution through a sequence of invertible neural-network transformations. The result can approximate non-Gaussian posteriors much more accurately than ADVI.

Stein variational inference, neural-network-based VI, and amortized inference (where a neural network learns to produce variational parameters given the data) are all active research areas. In production, the standard choice is still ADVI or Pathfinder; the more exotic variants have not yet become standard tools.

Part VII. Other Approximations

Laplace approximation

The Laplace approximation finds the posterior mode and approximates the posterior as a Gaussian centered at the mode with covariance equal to the inverse of the negative log-posterior Hessian at the mode. It is essentially a second-order Taylor expansion of the log posterior.

For unimodal, approximately Gaussian posteriors, the Laplace approximation is fast and reasonably accurate. For posteriors that are skewed or heavy-tailed, it can be significantly biased. In modern practice, the Laplace approximation is most commonly used as a building block inside larger algorithms (notably INLA) rather than as a standalone method.

INLA

Integrated Nested Laplace Approximation (Rue, Martino, Chopin 2009) is a deterministic approximation for the class of latent Gaussian models. The class includes generalized linear mixed models, spatial models, time-series models, and many others. The defining feature is that the parameters can be split into a high-dimensional Gaussian latent field and a small number of non-Gaussian hyperparameters.

INLA exploits this structure by computing accurate Laplace approximations for the latent field conditional on the hyperparameters, then numerically integrating over the hyperparameters. The result is a deterministic, fast, and accurate approximation that is competitive with MCMC for the class of models it supports.

The R-INLA package is the standard implementation. INLA is the right choice when your model fits its structural assumptions and you need to run many instances at scale. The newer inlabru package extends INLA to a wider class of models.

Subsampling MCMC

For very large datasets, the cost of evaluating the likelihood at every iteration becomes the bottleneck. Subsampling MCMC methods (stochastic gradient MCMC, minibatch MH, and others) replace the full likelihood evaluation with an unbiased estimate based on a minibatch. The result trades bias for speed; the chain no longer samples exactly from the posterior, but if the bias is small relative to the variance, the effective rate of useful information per second can be much higher.

Subsampling MCMC is mostly relevant for problems with hundreds of thousands or millions of observations. For typical industrial problems (tens of thousands of weekly observations in MMM, tens of thousands of users in CLV) the standard NUTS is fast enough that subsampling is unnecessary.

Part VIII. Debugging in Practice

One useful debugging technique in PyMC is pm.model_to_graphviz(model), which produces a graphical representation of the model’s structure. This catches mistakes in how parameters are connected, which is a common source of subtle bugs. It is one tool among several: pair plots of divergent samples, prior predictive simulation, building the model up incrementally (add one component at a time and rerun), checking generated quantities against hand-computed values, and unit-testing individual likelihood pieces against scipy or analytic references are all equally important parts of a serious debugging workflow.

For inspecting intermediate computations, PyMC uses symbolic graphs (PyTensor), so a plain print statement inside the model definition will not print runtime values. The right tool is pytensor.printing.Print, which inserts a print node into the computational graph:

from pytensor.printing import Print my_var_printed = Print(‘my_var’)(my_var)

This will print the value of my_var every time the graph is evaluated, which is useful for catching NaN or infinity values that propagate through the model.

Common runtime errors and their fixes:

• Divergences after warmup: increase target_accept to 0.95 or 0.99. If that does not work, reparameterize. If that does not work, the model is misspecified.
• Chain failed to initialize: the initial point is outside the support of one of the priors. Pass init=0 or initvals={} with valid starting points.
• Bad initial energy: the log-posterior at the initialization is -inf or nan. Same fix as above, plus check whether any of your data have illegal values for the likelihood (negative counts in a Poisson, for example).
• Max treedepth exceeded: the posterior is very long and thin. Reparameterize to a more isotropic shape, or increase max_treedepth (less ideal).
• Sampler is very slow: the model has expensive likelihood evaluations or many parameters. Profile with pm.compile_pymc or switch to NumPyro for JAX acceleration.

Part IX. When to Use Which Tool

The decision tree below is the working reference.

When the model is conjugate (Beta-Bernoulli, Gamma-Poisson, Normal-Normal with known variance)

Use the closed-form update. There is no reason to fire up NUTS. The closed form is exact, fast, parallel, and trivial to ship to production. This is the right choice for production Thompson Sampling, for CLV models (Pareto/NBD, BG/NBD), and for any place where you need millisecond-latency posterior updates.

When the model is a standard GLM, GLMM, or regression with random effects

Use NUTS through PyMC, Stan, or NumPyro. Bambi or brms will write the code for you if the model fits the formula syntax. ADVI is a useful intermediate step for iteration, but the final fit should be NUTS.

When the model is hierarchical with sparse within-group data

Use NUTS with non-centered parameterization. This is the most common failure mode for novices: the centered parameterization works in toy examples and breaks on real data. Build the non-centered version from the start.

When the model is custom (state-space, MMM, neural with Bayesian heads)

Write it directly in PyMC, Stan, or NumPyro. Bambi and brms will not help here. Use NUTS for the final fit, ADVI or Pathfinder for iteration.

When you have millions of observations or need GPU acceleration

NumPyro on JAX is the right tool. The JAX compilation makes the gradient evaluations fast on GPU, and NUTS scales naturally. BlackJAX is the alternative if you want more control over the sampler internals.

When the model is a latent Gaussian model (spatial, temporal, GLMM)

INLA is competitive with NUTS for accuracy and dramatically faster. The R-INLA package is the standard.

When you need a fast iteration loop on model structure

ADVI or Pathfinder for the initial fits, NUTS for the final result. Switching between them in PyMC is changing one function call.

Part X. Operational Assumptions

Every inference algorithm makes assumptions, and most failures of inference are violations of those assumptions rather than bugs in the algorithm.

MCMC assumes ergodicity. The chain must be able to reach every region of the posterior eventually. This breaks for multimodal posteriors with widely separated modes, where the sampler can get stuck in one mode for impossibly long times.

HMC and NUTS assume differentiability. The log-posterior must be differentiable with respect to the parameters. This excludes models with discrete latent variables unless they are marginalized out, and models with non-smooth likelihoods.

ADVI mean-field assumes posterior factorization. Parameters are assumed to be approximately independent in the posterior. This breaks when there are strong correlations (which is most of the time in interesting models).

ADVI full-rank assumes Gaussianity. The posterior is approximated as a Gaussian. This breaks for multimodal or strongly non-Gaussian posteriors.

Laplace approximation assumes unimodality. Same caveat as ADVI full-rank, but more aggressive because it only uses the mode.

INLA assumes latent Gaussian structure. The model must have the right structure for the algorithm to apply.

Simulation-based calibration assumes exact samples. Autocorrelated MCMC samples need to be thinned before SBC; otherwise SBC will produce false positives of calibration error.

For moderate-dimensional differentiable models in classical industrial Bayesian work, NUTS remains a strong default reference method. It makes few assumptions, scales well beyond a few dimensions, and its diagnostics are the most informative when something goes wrong. The picture is broader at the frontier: SMC and particle MCMC, deep variational models, amortized inference, simulation-based inference, diffusion-based posterior samplers, and probabilistic transformers are growing fast and are the right choice in their respective regimes. For very high-dimensional VI systems, large-scale amortized inference, SVI in deep probabilistic models, latent-variable architectures, and simulation-based inference, NUTS is either too expensive or simply not applicable, and the appropriate method comes from the relevant section above.

Part XI. FAQ for Computational Diagnostics

Q: NUTS gives me divergences after warmup. What do I try first?

Increase target_accept to 0.95. If that does not reduce divergences below a handful, increase to 0.99. If that still leaves divergences, reparameterize. The single most likely fix is the non-centered parameterization for any hierarchical structure. After that, look at the pair plot of divergent samples; the pattern of divergences will tell you where the geometric pathology is.

Q: R-hat is 1.00 everywhere. Do I trust the posterior?

Only if bulk-ESS and tail-ESS are also healthy (above 400 per chain), and you have run at least four chains. R-hat detects gross non-convergence; it can miss subtler failures where all chains are stuck in the same local mode.

Q: My chains are sampling but the trace plot looks like a hairball.

Good. That is what a converged chain should look like. A pretty smooth trace usually indicates a problem with mixing.

Q: My chains seem to converge but I see one chain stuck in a different region.

Either a true multimodal posterior or a label-switching problem in a mixture model. For label-switching, impose an ordering constraint on the components (using pm.Potential with a penalty for unordered labels, or with the ordered transform). For true multimodality, NUTS will struggle; consider tempered transitions or parallel tempering.

Q: ADVI converges in 30 seconds but says something different from NUTS.

In most practical cases, trust the NUTS fit unless there is strong evidence the chains failed to converge. ADVI has just told you that the mean-field approximation is poor for your model. Use ADVI for iteration on structure, NUTS for results.

Q: My Stan model compiles but fails to initialize.

The initial parameter values are outside the support of the priors or the likelihood. Specify init=0 to start at the constrained-space zero, or pass explicit initial values that are valid. Watch for log-transformed parameters where a value of zero is constrained-space zero but unconstrained-space infinity.

Q: NumPyro on GPU is slower than PyMC on CPU.

The model is small enough that JAX compilation overhead dominates. JAX on GPU pays off when the gradient evaluation is expensive, typically for models with many observations or many parameters. For small models, CPU NUTS is fast enough.

Q: My model has tens of thousands of parameters and NUTS is taking forever.

Three options. First, profile to see where the time is being spent. Second, switch to NumPyro on GPU; the gradient evaluation scales much better with hardware. Third, accept that NUTS is overkill and use ADVI with the understanding that the variance estimates will be too narrow.

Q: My posterior is bimodal. NUTS only samples one mode.

NUTS cannot easily jump between widely separated modes. Three options: reparameterize so the modes are closer in the sampling space; use parallel tempering; or accept that the bimodality is a model problem and respecify the model so it is unimodal.

Q: What is a divergence telling me about my model versus my sampler?

A divergence is a sampler-level event but it is informative about the model. Divergences happen where the posterior geometry is hard for the integrator. The geometry is a property of the model, so the divergence is telling you that your model produces a posterior with that hard geometry. Reparameterization fixes the geometry without changing the model. If reparameterization is not enough, the geometry induced by the model specification is likely the problem: heavy tails, multimodality, weak identifiability, pathological priors, or numerical instability in the likelihood are all candidates. Sampler settings (target_accept, max_treedepth, warmup length, mass-matrix choice) can also still matter and should be checked before declaring the model unsalvageable.

Part XII. Advanced Variational Inference

ADVI with mean-field or full-rank Gaussian is the entry-level VI. The frontier has moved well past it. Modern variational inference uses more flexible variational families, more careful optimization, and methods that combine the speed of VI with calibration closer to MCMC.

Normalizing Flows

A normalizing flow (Rezende and Mohamed 2015) constructs a flexible distribution by transforming a simple base distribution through a sequence of invertible, differentiable mappings. The resulting density has a closed form via the change-of-variables formula:

The flow can represent multimodal, skewed, and otherwise non-Gaussian distributions while remaining tractable for VI.

Several flow families are operationally important. Planar and radial flows are the original constructions, simple but limited. Inverse Autoregressive Flow (IAF, Kingma et al. 2016) uses autoregressive transformations and is the dominant family for high-dimensional variational posteriors. Real NVP and Glow use coupling layers and are the dominant family for generative modeling. Neural Spline Flows use rational quadratic transformations and have stronger expressiveness for the same parameter count.

For variational inference of complex posteriors, IAF or Neural Spline Flow as the variational family produces calibration much closer to MCMC than mean-field VI, at modest computational cost. The implementations are in Pyro and NumPyro through the AutoNormalizingFlow guide.

Stein Variational Gradient Descent

SVGD (Liu and Wang 2016) is a non-parametric variational method that maintains a set of particles and updates them along a deterministic flow that minimizes the KL divergence to the posterior. The update rule combines a gradient component (move particles toward higher posterior density) and a repulsion component (keep particles apart to maintain diversity):

The first term inside the sum pushes toward the mode; the second term pushes particles apart. The kernel k controls the smoothing.

SVGD is more flexible than parametric VI (no Gaussianity assumption) and cheaper than MCMC (deterministic, no rejections). Empirically it produces good results for moderate-dimensional posteriors. The libraries pyro.infer.SVGD and numpyro.infer.SteinVI implement it.

Pathfinder

Pathfinder (Zhang, Carpenter, Gelman, Vehtari 2021), already introduced in Part VII, deserves more detail. The method constructs a sequence of Gaussian approximations along the trajectory of an L-BFGS optimization, picks the best by importance sampling, and uses Pareto smoothing to validate the choice. The output is a single Gaussian approximation that is typically much better than ADVI mean-field and converges in seconds.

Pathfinder is the workhorse modern VI tool when you need a fast approximation that you can trust. Its combination with NUTS (Pathfinder for initialization, NUTS for refinement) is the recommended workflow in PyMC and Stan for hard problems.

Amortized variational inference

Amortized inference replaces the per-observation variational parameters with a neural network that predicts them from the data. The procedure is:

1. Define a variational family q(z | x; phi) parameterized by a neural network with weights phi.

2. Train phi to maximize the ELBO across the dataset.

3. At inference, for each new x, evaluate the network to get the variational parameters.

The benefit is dramatic at inference: no per-instance optimization, just a forward pass through the network. The cost is that the network must be expressive enough to capture the conditional structure.

Variational autoencoders (Kingma and Welling 2014) are the canonical amortized inference architecture. They are used as building blocks in many production systems, particularly for representation learning, generative modeling, and any setting where inference at scale is critical.

Part XIII. Sequential Monte Carlo and Particle Filters

When the model has temporal structure (state-space, time-series with latent dynamics, sequential decision making), the natural inference paradigm is sequential. Sequential Monte Carlo (SMC) methods, also called particle filters in the state-space context, are the standard tools.

The basic particle filter

The state-space model is

with latent states z and observations y. The goal is the filtering distribution p(zt | y{1:t}) at each time t.

The bootstrap particle filter (Gordon, Salmond, Smith 1993) maintains a set of N particles representing the current filtering distribution. At each step, the particles are propagated forward through the dynamics, weighted by the new observation likelihood, and resampled in proportion to the weights. The result is an approximation of the filtering distribution that adapts to non-linear dynamics and non-Gaussian observations.

Particle filters are the right tool for any state-space model where the Kalman filter assumptions (linear Gaussian) do not hold. Applications include object tracking, financial volatility modeling, biological state estimation, and robotics SLAM.

SMC samplers

SMC samplers (Del Moral, Doucet, Jasra 2006) generalize particle filters from sequential observation to sequential tempering. The idea is to start with particles drawn from the prior and transition them through a sequence of intermediate distributions that interpolate from the prior to the posterior. At each step, the particles are reweighted, resampled, and moved via MCMC kernels.

SMC samplers are competitive with NUTS for difficult posteriors, particularly multimodal ones. The tempering schedule helps the particles cross between modes. They are also naturally parallelizable across particles.

The Pyro library implements SMC samplers through pyro.infer.SMC. The Stan ecosystem provides similar functionality through external packages.

Auxiliary particle filters

The auxiliary particle filter (Pitt and Shephard 1999) improves on the bootstrap filter by proposing particles using a one-step-ahead prediction that incorporates the next observation. This reduces the variance of importance weights and improves the effective sample size.

For state-space models in production (stochastic volatility, hidden Markov models with continuous states, target tracking), auxiliary particle filters are typically the right tool.

Pseudo-marginal MCMC

Pseudo-marginal MCMC (Andrieu and Roberts 2009) is a remarkable result that extends MCMC to settings where the likelihood is intractable but can be estimated unbiasedly. The standard MH acceptance ratio is replaced by an estimate of the likelihood ratio. Remarkably, as long as the likelihood estimator is unbiased, the resulting chain has the correct posterior as its stationary distribution, regardless of how noisy the estimator is.

The implications are large. Pseudo-marginal MCMC enables Bayesian inference in models where the likelihood involves an intractable integral (latent variable models, doubly intractable models with normalizing constants, some PDE-constrained models). The standard pattern uses particle filters to estimate the likelihood within an MH outer loop, producing the “particle MCMC” framework (Andrieu, Doucet, Holenstein 2010).

Part XIV. Likelihood-Free and Simulation-Based Inference

The Bayesian models covered in this article have explicit likelihood functions p(y | theta). Many scientific and industrial models do not. A complex simulator (epidemiological model, ecological model, financial market model, physics simulator) can generate data from parameters but does not provide a tractable likelihood. Likelihood-free inference covers the methods for Bayesian inference in this regime.

Approximate Bayesian Computation

ABC is the historical entry point. The basic ABC rejection sampler is:

1. Sample theta from the prior.

2. Simulate data y_sim ~ simulator(theta).

3. Accept theta if distance(y_sim, y_obs) < epsilon.

4. Repeat.

The accepted thetas approximate the posterior conditional on the distance being below epsilon. For epsilon = 0, this is exact Bayesian inference. For epsilon > 0, the procedure is approximate, with the bias decreasing as epsilon decreases.

The standard refinements of ABC (sequential Monte Carlo ABC, ABC with regression adjustment, semi-automatic summary statistics) make the basic procedure more efficient. The pyabc Python library implements modern ABC.

ABC is the right tool for simple simulator-based problems with low-dimensional parameters and informative summary statistics. It scales poorly with the dimension of the parameter and the dimension of the data.

Simulation-Based Inference with Neural Networks

Modern simulation-based inference (SBI, Cranmer, Brehmer, Louppe 2020) replaces the rejection step with a neural network that learns the posterior, the likelihood, or the likelihood ratio from simulations. The standard variants are:

Neural Posterior Estimation (NPE): train a normalizing flow on simulated (theta, y) pairs to learn p(theta | y) directly. At inference, condition the flow on observed y to get the posterior.

Neural Likelihood Estimation (NLE): train a flow to learn p(y | theta). At inference, use MCMC or VI with the learned likelihood.

Neural Ratio Estimation (NRE): train a classifier to distinguish (theta, y) pairs from theta and y drawn independently. The classifier’s output recovers the likelihood ratio, which is sufficient for inference.

The sbi library implements all three. SBI scales to higher-dimensional problems than ABC and is the modern standard for likelihood-free Bayesian inference.

The applications are dominantly in scientific domains (cosmology, neuroscience, epidemiology, particle physics) but also in some industrial settings where the data-generating process is described by a complex simulator (financial risk models, demand simulators).

Part XV. Riemannian HMC and Geometric MCMC

Standard HMC uses a Euclidean metric: the kinetic energy is r^T r / 2. Riemannian HMC (Girolami and Calderhead 2011) uses a position-dependent metric, with the local mass matrix determined by the local curvature of the log-posterior. The result is dramatically better performance on posteriors with strong local correlations or pathological geometry.

The cost is the computation and inversion of the metric at each step, which is expensive. RHMC is therefore most useful for low-dimensional problems with hard geometry rather than for general-purpose inference.

The successor methods (Lagrangian Monte Carlo, Magnetic HMC) explore further geometric variants. None has displaced NUTS in production, but for specific difficult posteriors (Bayesian inverse problems, models with deep ridges) they are the appropriate tools.

Part XVI. Probabilistic Programming Language Internals

The user-facing API of PyMC, Stan, and NumPyro hides a substantial amount of machinery. For senior practitioners, particularly those working on infrastructure or developing custom inference algorithms, the internals matter.

Effect handlers in Pyro

Pyro is built on the concept of effect handlers (Pyro’s poutine module). An effect handler is a context manager that intercepts calls to pyro.sample and pyro.param and modifies their behavior. Effect handlers are how Pyro implements:

· Conditioning (the condition handler observes the value of a sample site).

· Tracing (the trace handler records every sample and its log probability).

· Replay (the replay handler reuses recorded values at sample sites).

· Block (the block handler hides sample sites from outer handlers).

· Plates (the plate handler annotates conditional independence).

Custom inference algorithms in Pyro are typically composed by stacking effect handlers. Writing a new sampler often comes down to defining a new handler. This compositional design is one of Pyro’s distinctive features and a marker of probabilistic programming sophistication.

Plates and conditional independence

A plate in Pyro and NumPyro annotates a batch dimension over which observations are conditionally independent given the parameters. The plate context manager tells the inference engine that the sample statements inside the plate are independent across the plate dimension. This enables vectorization, parallelization, and efficient gradient computation.

Plates can be nested for multi-dimensional structures (subjects within groups within experiments). The plates correspond to the natural conditional independence structure of hierarchical models.

AutoGuide families

Pyro and NumPyro provide a family of “guides” that automatically construct variational distributions. The hierarchy is:

· AutoNormal: mean-field Normal.

· AutoMultivariateNormal: full-rank Normal (joint).

· AutoLowRankMultivariateNormal: low-rank + diagonal Gaussian.

· AutoNormalizingFlow: normalizing flow.

· AutoLaplaceApproximation: Laplace at the MAP.

· AutoDelta: point estimate (MAP).

The right choice depends on the posterior. AutoNormal is the entry point. AutoLowRankMultivariateNormal is a good compromise between expressiveness and parameter count. AutoNormalizingFlow is the most expressive but most expensive.

NumPyro JAX primitives

NumPyro is implemented on JAX, which gives access to functional transformations: jit for compilation, vmap for vectorization, pmap for multi-device parallelism, grad for automatic differentiation. The full power of JAX is available within probabilistic programs.

The patterns that matter for production:

· lax.scan for efficient loops over time-series.

· lax.cond for differentiable branching.

· vmap over plates for vectorized inference.

· pmap over data shards for distributed inference.

· jit on the entire inference loop for compiled performance.

For large-scale Bayesian inference on GPU clusters, NumPyro on JAX is the most mature platform in 2026.

Stan transformed parameters and generated quantities

Stan’s structure is more rigid than Pyro’s but has its own depth. The four blocks (data, parameters, transformed parameters, model, generated quantities) correspond to a specific computational pipeline:

· transformed parameters is computed once per leapfrog step and is autodiffed.

· generated quantities is computed only at saved iterations and is not autodiffed.

The distinction matters operationally. Quantities used in the model (priors, likelihoods) belong in transformed parameters; post-inference quantities (predictions, derived statistics, posterior checks) belong in generated quantities. Mixing them up either slows down inference dramatically or produces incorrect autodiffed gradients.

Stan’s user-defined functions (defined in the functions block) provide modularity. Custom likelihoods (target += … syntax) are how non-standard models are implemented.

Part XVII. Expectation Propagation and Variational Message Passing

Beyond MCMC and variational inference, a third paradigm of approximate inference exists. Message passing methods exploit the conditional independence structure of probabilistic graphical models to compute (or approximate) marginals by passing messages along the graph. The two dominant message-passing methods for Bayesian inference are Variational Message Passing (VMP) and Expectation Propagation (EP).

Probabilistic graphical models and factor graphs

Factor graph: variables (circles) and factors (squares) of a probabilistic model

Before discussing the algorithms, the underlying structure deserves coverage. A probabilistic graphical model decomposes a joint distribution into a product of factors:

where each factor f_k depends on a subset S_k of the variables. The structure is naturally represented as a factor graph: a bipartite graph with variables on one side and factors on the other, with edges connecting each factor to the variables it involves.

Directed graphical models (Bayesian networks) and undirected graphical models (Markov random fields) are both special cases of factor graphs. The factor graph is the most general representation and is the natural data structure for message-passing algorithms.

For models with low treewidth (the graph has limited connectivity), exact inference is tractable via the junction tree algorithm. The procedure constructs a tree of cliques and passes messages between them to compute exact marginals. For high-treewidth or loopy graphs, exact inference is intractable and approximations are needed.

Belief propagation and the sum-product algorithm

The sum-product algorithm is the canonical message-passing procedure. Each factor sends a message to each connected variable, summarizing the factor’s belief about that variable. Each variable sends a message to each connected factor, combining incoming messages from other factors. The procedure iterates until convergence, at which point the variable marginals can be read off from the products of incoming messages.

On tree-structured graphs, sum-product converges in one forward and one backward pass and produces exact marginals. On loopy graphs, the procedure is loopy belief propagation, which is approximate but often empirically effective.

The max-product algorithm is the analogue for MAP inference: replace sum with max in the message updates, and the procedure computes the most likely joint configuration. The Viterbi algorithm for hidden Markov models is a special case.

Variational Message Passing

VMP (Winn and Bishop 2005) reformulates variational inference as message passing on a factor graph. The variational distribution factorizes across nodes in the graph, and the optimal updates have a message-passing form: each node updates its variational parameters based on messages from its neighbors.

The result is the same as mean-field VI in the final objective, but the implementation is much more modular. Adding a new model is a matter of specifying the factors; the message-passing engine handles the inference automatically. This is the design behind Infer.NET (Microsoft Research), one of the most influential probabilistic programming systems for industrial use.

VMP excels on models that factorize naturally as graphical models: large hierarchical models, factor analysis, latent variable models, and most models with conjugate exponential family structure. For these, VMP is often much faster than NUTS and produces comparable accuracy.

Expectation Propagation

EP (Minka 2001) is an alternative message-passing method that, instead of minimizing reverse KL divergence (as VI does), minimizes forward KL divergence between approximate and true conditionals at each factor. The result is an approximation that tends to be mode-covering rather than mode-seeking: EP produces broader posteriors than VI, often closer to the true posterior in calibration.

The procedure iteratively removes a factor’s contribution to the approximate posterior, computes the true conditional with that factor reinstated, projects back to the approximating family, and updates the factor’s message. EP is more delicate in practice than its clean derivation suggests: convergence is not guaranteed, oscillations and outright divergence can occur for non-Gaussian or skewed factors, damping and double-loop variants (Heskes and Zoeter 2002) are often necessary to stabilize the updates, and there are no general optimization-style guarantees. Practitioners typically tune the damping factor empirically and fall back to VI or MCMC when EP fails to settle.

When it does settle, EP is the standard inference engine for Gaussian Process classification (where VI struggles with the heavy tails of the posterior). It is also the foundation of several long-running production systems at Microsoft Research (AdPredictor, the Bing ad ranking system).

When to use message passing

Message passing is the right tool when:

· The model has natural factor graph structure with low or moderate treewidth.

· The factors are tractable (exponential family with conjugate structure).

· Modularity and extensibility matter (adding a new factor type should not require rewriting the inference).

· Calibration matters more than mode-seeking (use EP rather than VI).

· The model fits the Infer.NET paradigm.

The libraries are Infer.NET (.NET ecosystem), Stan (limited support), Pyro (via custom guides), and dimod for graphical model inference broadly. For specifically GP classification, EP is built into GPyTorch and GPflow.

Part XVIII. Slice Sampling

Slice sampling (Neal 2003) is a general-purpose MCMC algorithm that requires only the ability to evaluate the unnormalized posterior density and does not require gradients or hand-tuned step sizes. The procedure introduces an auxiliary variable that uniformly samples below the density and steps along level sets to explore.

The basic algorithm:

1. Sample an auxiliary y uniformly between 0 and the current density value f(x).

2. Find the level set {x : f(x) > y} (the “slice”).

3. Sample x uniformly from the slice.

4. Return the new x as the next state.

The level set is typically found by stepping out from the current point until the density falls below y, then shrinking the bracket until a sample is accepted.

Why slice sampling matters in practice

Slice sampling has several practical advantages:

· No step size tuning: the step-out and shrink procedure adapts automatically.

· No gradient required: works on models where HMC and NUTS cannot apply.

· Robust to scale differences: handles parameters at very different scales without preconditioning.

· Simple to implement: substantially less code than a tuned HMC sampler, with fewer moving parts.

The cost is that slice sampling mixes slower than HMC on smooth gradient-friendly posteriors. For high-dimensional well-conditioned problems, NUTS is much faster.

The right use of slice sampling in modern Bayesian practice is as a fallback when NUTS does not apply: models with discrete latent variables that cannot be marginalized, models with non-differentiable likelihoods, models where automatic differentiation is broken or unstable. The PyMC and Stan ecosystems both include slice samplers as alternative kernels.

Part XIX. Parallel Tempering for Multimodal Posteriors

NUTS cannot easily jump between widely separated modes of the posterior. When the model has a multimodal posterior with deep wells separated by low-density regions, standard MCMC samplers get stuck in one mode and produce biased results.

Parallel tempering (Geyer 1991, Earl and Deem 2005) addresses this by running multiple chains in parallel at different “temperatures.” A chain at temperature T samples from a tempered posterior p(theta | y)^{1/T}, which flattens the posterior and makes it easier to traverse modes. Periodically, neighboring chains in temperature swap states with a Metropolis acceptance criterion that maintains correct stationarity.

The high-temperature chains explore the full state space freely; the low-temperature chains target the actual posterior; the swap moves let modes be discovered at high temperature and refined at low temperature.

When parallel tempering is needed

Parallel tempering is the right tool when:

· The posterior is genuinely multimodal with separated modes (not label-switching, which has its own fixes).

· The modes correspond to substantively different model interpretations.

· Standard NUTS gets stuck in one mode and the chains do not agree.

· The compute budget allows multiple parallel chains at different temperatures (typically 4 to 16).

Multimodal posteriors arise in mixture models (without identifiability constraints), in physics-based models with multiple local equilibria, in models with discrete latent structure, and in some economic structural models. For these, parallel tempering or its variants (simulated tempering, Wang-Landau, stochastic approximation Monte Carlo) is the appropriate approach.

The implementations are in emcee (Goodman-Weare ensemble samplers with tempering), custom Stan code, and the dedicated parallel tempering libraries in physics packages.

Part XX. Streaming Bayesian Inference

The MCMC and VI methods covered so far assume a fixed dataset. In production, data arrive continuously, and the model must update its posterior incrementally without re-running the full inference each time.

The conjugate path

For conjugate models, streaming Bayesian inference is exact and trivial. The posterior after observing data D is the prior for the next observation, and the conjugate update is one closed-form computation per new datum. Beta-Bernoulli for click-through rates, Gamma-Poisson for event counts, Normal-Normal for continuous metrics: each supports millisecond-latency updates and is the natural fit for streaming bandits, real-time A/B tests, and online dashboards.

Variational filtering

For non-conjugate models, the streaming analogue of VI is sequential variational inference. The posterior after observing D_{1:t-1} serves as the prior for observation t. The variational distribution is updated by minimizing the KL divergence to the posterior given the new datum.

The procedure is straightforward in principle and computationally light: each update is a few gradient steps on the ELBO conditional on the new observation. The challenges are:

· Catastrophic forgetting: the variational parameters can drift if the model is non-conjugate, losing information from old data.

· Variance underestimation accumulating: each variational update slightly underestimates the variance, and the bias accumulates over many updates.

The standard fixes are to maintain a buffer of recent observations and re-process them periodically, and to apply variance correction factors that compensate for accumulated bias.

Particle filters for streaming

Sequential Monte Carlo methods (Part XIII) are the natural framework for streaming inference on state-space models. The particle filter maintains a set of particles representing the current filtering distribution and updates them as each new observation arrives. The cost is constant per observation, making it well-suited to high-throughput streaming.

For models that fit the state-space form (target tracking, online demand forecasting, real-time risk monitoring), particle filters are the production tool.

Streaming in production

The patterns in production:

· Online learning systems (bandits, recommenders, fraud detection) use conjugate updates where possible.

· Real-time monitoring uses Bayesian Online Change Point Detection (covered in Article 6).

· State-space models in production use particle filters or extended Kalman filters with Bayesian uncertainty.

· Periodic batch reinference is the most common pattern overall: run the full Bayesian inference on a recent window of data on a schedule (daily, hourly), serve the result until the next refresh.

The libraries that support streaming Bayesian inference include river (online learning with Bayesian variants), tensorflow-probability (sequential Monte Carlo), pyro.contrib.tracking (particle filtering), and custom production code at Meta, Netflix, Amazon, and others.

Closing

Computational inference for Bayesian models is a small set of tools used well. NUTS for the general case, conjugate updates where they apply, variational inference for iteration, INLA for latent Gaussian models. The skill is not memorizing dozens of algorithms; it is knowing which tool to reach for and how to diagnose when it is failing.

In the next article we move from general inference to specific industrial models. Hierarchical models and Bayesian GLMs are the workhorses of every production system: pricing, churn, customer lifetime value, demand forecasting, retention. The structure of these models, and the patterns that make them robust at scale, are the operational content of the next part.

References

· Hoffman, Gelman. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, JMLR (2014).

· Betancourt. A Conceptual Introduction to Hamiltonian Monte Carlo, arXiv 1701.02434.

· Vehtari, Gelman, Simpson, Carpenter, Bürkner. Rank-Normalization, Folding, and Localization: An Improved R-hat for Assessing Convergence of MCMC, Bayesian Analysis (2021).

· Talts, Betancourt, Simpson, Vehtari, Gelman. Validating Bayesian Inference Algorithms with Simulation-Based Calibration, arXiv 1804.06788.

· Kucukelbir, Tran, Ranganath, Gelman, Blei. Automatic Differentiation Variational Inference, JMLR (2017).

· Zhang, Carpenter, Gelman, Vehtari. Pathfinder: Parallel quasi-Newton variational inference, arXiv 2108.03782.

· Rue, Martino, Chopin. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations, JRSS-B (2009).