Monte Carlo and Sampling Methods

The previous module closed with a family of approximate inference methods that are stochastic rather than deterministic: rather than optimizing over a restricted family of distributions, they generate random samples and use those samples' empirical behavior to approximate a quantity that cannot be computed in closed form. The general principle, Monte Carlo estimation, rests on the law of large numbers: if X1, ..., Xn are samples drawn from a distribution p, then the average of any function f applied to those samples converges, as n grows, to the true expectation of f under p. This lets an intractable integral or an exponentially large sum be replaced by a sample average whose accuracy is controllable simply by drawing more samples - at the cost of introducing random error that shrinks, but never fully disappears, as the sample count grows. The 1953 paper by Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller is a foundational point of origin for this style of sampling-based computation, applying it to simulate the equilibrium behavior of physical systems on early electronic computers - a use case that shares essentially the same underlying logic that later Bayesian applications would inherit.

Monte Carlo estimation as just described assumes samples can be drawn directly from the distribution of interest, p, which is often not true - p might be known only up to an unknown normalizing constant, or might simply be too complicated to sample from directly. Importance sampling addresses this by sampling instead from a different, easier-to-sample proposal distribution q, then correcting for the mismatch by reweighting each sample xi by the importance weight p(xi)/q(xi). The resulting weighted average remains an unbiased estimator of the expectation under p, provided q assigns nonzero probability everywhere p does. The catch is that if q is a poor match for p - assigning little probability mass to regions where p places substantial mass - the importance weights become highly unequal, a small number of samples end up dominating the weighted average, and the estimator's variance can be enormous even though it remains technically unbiased.

A particularly useful variant, self-normalized importance sampling, handles the common Bayesian situation where p (typically an unnormalized posterior) is known only up to its normalizing constant: rather than requiring that constant explicitly, the importance weights are normalized to sum to one across the drawn samples, and the resulting weighted average remains a consistent (though now slightly biased for finite samples) estimator of the target expectation, with the bias shrinking toward zero as more samples are drawn. This matters because, in Bayesian inference, the posterior is almost always known only as prior times likelihood, an unnormalized product whose normalizing constant (the marginal likelihood, obtained by integrating that product over all possible parameter values) is usually the single hardest quantity in the whole computation to obtain directly. Self-normalized importance sampling sidesteps that difficulty entirely, and this is exactly the situation exploited later in this domain by probabilistic programming systems that construct an unnormalized joint likelihood directly from a program's sample and observe statements without ever computing its normalizing constant in closed form.

Importance sampling can fail badly in high-dimensional spaces, where almost any fixed proposal distribution ends up placing negligible weight in the regions the target distribution actually cares about. Markov chain Monte Carlo (MCMC) takes a different approach: instead of drawing independent samples from a proposal and correcting them by reweighting, it constructs a Markov chain whose stationary distribution is exactly the target distribution, then simulates that chain and treats the (correlated) sequence of visited states as samples. The 1953 Metropolis algorithm did exactly this for the specific case of a symmetric proposal distribution: propose a candidate next state x' from the current state x using a proposal that is equally likely to propose x from x' as it is to propose x' from x, and accept the move with probability min(1, p(x')/p(x)) - always accepting a move to a more probable state, and sometimes accepting a move to a less probable one, which is exactly what lets the chain explore the whole distribution rather than getting stuck at a single mode.

Hastings's 1970 generalization removed the requirement that the proposal be symmetric. The Metropolis-Hastings acceptance rule, accept with probability min(1, [p(x')q(x|x')] / [p(x)q(x'|x)]), adds a correction ratio q(x|x')/q(x'|x) that exactly compensates for any asymmetry in how the proposal favors certain moves over others, while preserving the same essential guarantee via a property called detailed balance: the target distribution remains the chain's unique stationary distribution regardless of what proposal is used, as long as the acceptance rule correctly compensates for it. Crucially, both versions of the rule only ever require the target distribution's value up to an unknown multiplicative constant, since any constant cancels in the ratio p(x')/p(x) - precisely the situation Bayesian posterior inference is normally in, since the posterior's normalizing constant (the marginal likelihood) is usually the hardest part of the computation to obtain directly. Detailed balance itself is a local, pairwise condition - that the flow of probability from state x to state x' matches the flow back from x' to x under the target distribution - and it is straightforward to verify from the acceptance rule directly, which is part of why Metropolis-Hastings is comparatively easy to prove correct even for complicated, high-dimensional target distributions where almost nothing else about the distribution is easy to reason about analytically.

Gibbs sampling, applied by Geman and Geman to Bayesian image restoration, is a particularly convenient special case: rather than proposing a move for all variables at once and computing an acceptance ratio, it cycles through the variables one at a time, resampling each one directly from its full conditional distribution given the current values of all the others. Because sampling exactly from the full conditional is equivalent to a Metropolis-Hastings proposal that is always accepted, Gibbs sampling requires no explicit accept/reject step at all, provided those full conditional distributions are themselves tractable to sample from. In practice, using any MCMC method responsibly requires attention to burn-in (discarding early samples before the chain has approached its stationary distribution), mixing (how quickly the chain moves through the full support of the distribution rather than getting stuck near one region), and autocorrelation between consecutive samples, all of which affect how many effectively independent samples a given run of the chain is actually worth.

The state-space filtering problem from the previous module - maintaining a belief distribution over a hidden state as new observations arrive one at a time - has a closed-form exact solution, the Kalman filter, only under linear-Gaussian assumptions. When the dynamics or observation model are meaningfully nonlinear or non-Gaussian, as they routinely are in robotics and tracking applications with sensors that produce multimodal or heavy-tailed noise, that exactness is generally lost, and neither the discrete forward algorithm from the hidden Markov model setting nor the Gaussian Kalman recursion applies directly to recover it. Particle filtering, a sequential (online) form of Monte Carlo estimation, fills this gap by representing the filtered belief distribution not as a table of discrete state probabilities or a Gaussian mean and covariance, but as a weighted set of sample points, or particles, each representing one concrete hypothesis about the hidden state's current value, together forming an empirical approximation to the true, possibly complicated posterior shape.

Gordon, Salmond, and Smith's bootstrap filter, introduced as a general solution to nonlinear, non-Gaussian Bayesian state estimation, defines the now-standard sequential Monte Carlo recursion in three steps repeated at every time step. A prediction step propagates each particle forward through the (possibly nonlinear) dynamics model, sampling a new state for each particle from the transition distribution. An update (reweighting) step adjusts each particle's importance weight in proportion to how likely the new observation is given that particle's predicted state, using the observation model's likelihood - particles whose predicted state fits the new evidence well gain relative weight, and those that fit it poorly lose relative weight. A resampling step then duplicates high-weight particles and discards low-weight ones, drawing a fresh set of equally-weighted particles from the current weighted set, which the paper showed dramatically improves performance relative to simply carrying unequal weights forward indefinitely. This predict-reweight-resample cycle repeats at every new time step, so the particle set is continually refreshed to track wherever the true posterior currently has most of its probability mass, without ever requiring the dynamics or observation model to be linear or the noise to be Gaussian - the only real requirement is that both can be sampled from and evaluated pointwise.

Resampling exists to solve a specific failure mode called weight degeneracy: without it, after enough time steps nearly all of the particle set's total weight collapses onto a vanishingly small number of particles, so that the set as a whole stops usefully representing the posterior even though it nominally still contains many particles. Resampling is not free, however - repeatedly duplicating the same high-weight particles over many time steps can lead to particle impoverishment, in which the particle set loses diversity and starts to represent only a narrow region of the true posterior, an especially acute risk when the state space has many dimensions. The effective sample size, a diagnostic estimating how many independent samples the current weighted particle set is actually worth, is one standard tool for monitoring this tradeoff in practice. Particle filters are widely used for robot localization and target tracking precisely because they degrade gracefully in the nonlinear, non-Gaussian regimes where the Kalman filter's exact assumptions fail outright.