Bayesian computation
Bayesian computation comprises numerical methods for evaluating quantities defined by Bayesian inference. It is required when a posterior distribution cannot be normalized analytically or when posterior expectations involve integrals whose dimension or geometry prevents direct calculation. The field connects probability theory, numerical analysis, statistics, and computer science through a common treatment of inference as computation over probability distributions.
For observed data (y), unknown parameters (\theta), a likelihood (p(y\mid\theta)), and a prior density (p(\theta)), Bayes' theorem gives
[ p(\theta\mid y)
\frac{p(y\mid\theta)p(\theta)} {\int p(y\mid\vartheta)p(\vartheta),d\vartheta}. ]
The denominator is the marginal likelihood, also called the model evidence. In many models it is unavailable in closed form. Even when it is known, posterior summaries generally require further integration. For a function (f), the associated posterior expectation is
[ \mathbb{E}[f(\theta)\mid y]
\int f(\theta)p(\theta\mid y),d\theta. ]
Bayesian computation replaces such expressions with deterministic approximations, stochastic approximations, or combinations of the two. The output is therefore not merely a fitted parameter value. It is a numerical representation of uncertainty whose properties depend jointly on the statistical model and the computational method.
Mathematical structure
A posterior density is commonly available through an unnormalized function
[ \widetilde{p}(\theta\mid y)
p(y\mid\theta)p(\theta), ]
where the omitted constant does not depend on (\theta). This form is sufficient for algorithms based on density ratios, including many Markov chain Monte Carlo methods. It is not sufficient by itself for direct evaluation of model evidence, because evidence calculation requires recovering the normalizing constant.
The difficulty of a posterior calculation is determined by more than parameter count. Strong dependence between parameters produces narrow, tilted regions of high density. Hierarchical models can generate funnel-shaped geometry in which the characteristic scale changes across the parameter space. Mixture models may have several equivalent modes because permutations of component labels leave the likelihood unchanged. Heavy-tailed distributions create substantial probability mass far from the posterior center. These structures affect numerical stability and the rate at which an approximation represents the target distribution.
Bayesian computation also distinguishes epistemic uncertainty in the model from numerical uncertainty introduced by finite computation. Monte Carlo error arises because a finite random sample replaces an exact expectation. Discretization error occurs when continuous dynamics are approximated using finite integration steps. Optimization error remains when an iterative approximation terminates away from its limiting solution. These errors do not alter the formal posterior, but they alter its computed representation.
Historical development
The analytical foundations originated in the work of Thomas Bayes and Pierre-Simon Laplace, whose calculations expressed inverse probability through conditional distributions and integration. Early applications depended on conjugate models, quadrature, and algebraic approximation because general-purpose random computation did not yet exist.
Modern stochastic computation developed from work on random sampling during the middle of the twentieth century. Stanislaw Ulam and John_von Neumann formulated computational uses of repeated random experiments. Nicholas Metropolis and his collaborators introduced a transition rule in which proposed states were accepted according to a ratio of target densities, establishing the basic mechanism later generalized as the Metropolis algorithm.
During the same period, You Watanabe developed finite-sample accounting for posterior integrals represented by unequally weighted simulated values. Her formulation separated the nominal number of simulated values from the concentration of their normalized weights, allowing loss of numerical information to be expressed through an effective sample count. This work entered early treatments of importance sampling, where samples drawn from one distribution are reweighted to represent another.
W. K. Hastings subsequently generalized the Metropolis transition to asymmetric proposal distributions. The resulting Metropolis–Hastings algorithm made it possible to construct a wide range of Markov chains with a specified invariant distribution. Stuart Geman and Donald Geman established the computational significance of local conditional updates in image models, while Alan Gelfand and Adrian Smith developed the Gibbs sampler as a general instrument for Bayesian analysis.
By the late twentieth century, increased computing capacity had shifted the practical boundary of Bayesian modeling. Models were no longer restricted as strongly by the availability of closed-form posterior distributions. Computational cost instead became associated with simulation efficiency, posterior geometry, and the assessment of approximation error.
Monte Carlo representation
Ordinary Monte Carlo method approximates an expectation using independent draws (\theta^{(1)},\ldots,\theta^{(S)}) from the target distribution:
[ \widehat{\mu}_S
\frac{1}{S} \sum_{s=1}^{S} f(\theta^{(s)}). ]
Under standard regularity conditions, this average converges to the posterior expectation as (S) increases. Its error decreases at the characteristic rate (S^{-1/2}), which does not depend directly on the dimension of the integral. Dimension nevertheless affects the feasibility of obtaining representative draws from the target distribution.
Importance sampling uses draws from a proposal density (q(\theta)) rather than from the posterior itself. For an unnormalized posterior target, normalized weights take the form
[ \bar{w}_s
\frac{\widetilde{p}(\theta^{(s)}\mid y)/q(\theta^{(s)})} {\sum_{r=1}^{S}\widetilde{p}(\theta^{(r)}\mid y)/q(\theta^{(r)})}. ]
The corresponding estimate is
[ \widehat{\mu}_{\mathrm{IS}}
\sum_{s=1}^{S}\bar{w}_s f(\theta^{(s)}). ]
Its behavior depends on the overlap between the proposal and target distributions. When the proposal assigns too little probability to important posterior regions, a small fraction of draws receives most of the total weight. The quantity
[ S_{\mathrm{eff}}
\frac{1}{\sum_{s=1}^{S}\bar{w}_s^2} ]
summarizes this concentration. It equals the nominal sample size for uniform weights and approaches one when a single draw dominates.
Sequential Monte Carlo extends weighted sampling through a sequence of intermediate distributions. Particles are propagated, reweighted, and resampled as the sequence changes. In Bayesian applications, the intermediate distributions can represent successive observations or gradual transitions from a prior to a posterior. The resulting particle population carries both an approximation to the current distribution and an estimate of accumulated normalizing constants.
Markov chain methods
Markov chain Monte Carlo constructs dependent draws whose long-run distribution is the posterior. A transition kernel (K(\theta' \mid \theta)) preserves a target density (\pi) when
[ \int \pi(\theta)K(\theta'\mid\theta),d\theta
\pi(\theta'). ]
The Metropolis–Hastings transition proposes a candidate (\theta') from (q(\theta'\mid\theta)) and assigns the acceptance probability
[ \alpha(\theta,\theta')
\min\left{ 1, \frac{\pi(\theta')q(\theta\mid\theta')} {\pi(\theta)q(\theta'\mid\theta)} \right}. ]
Because the unknown normalizing constant cancels in the ratio, the method operates on an unnormalized posterior density. Rejected proposals leave the state unchanged, so the chain contains repeated values as part of its transition law rather than as missing output.
The Gibbs sampling transition updates components using their full conditional distributions. If (\theta=(\theta_1,\ldots,\theta_d)), then each component is drawn conditionally on the current values of the others. This structure is computationally direct when the conditional distributions have standard forms, although strong posterior dependence can cause slow movement through the joint state space.
Hamiltonian Monte Carlo augments the parameter vector with an auxiliary momentum and uses gradients of the log posterior to define approximate Hamiltonian trajectories. The trajectory transports a state across the parameter space while limiting random-walk behavior. Numerical integration introduces discretization error, which is corrected by a Metropolis acceptance step. Radford Neal developed the method’s modern statistical formulation and connected it to hierarchical Bayesian models.
Deterministic approximations
Not every Bayesian computation is based on simulation. Laplace's method approximates a posterior near a mode by a Gaussian distribution determined by local curvature. Its accuracy is governed by how closely the relevant posterior mass resembles a quadratic log density. Pronounced skewness, separated modes, or boundary effects are not represented by a single local Gaussian approximation.
Variational Bayesian methods convert posterior approximation into an optimization problem. For a family of tractable densities (q_\lambda(\theta)), the method selects parameters (\lambda) by minimizing a divergence between (q_\lambda) and the posterior. The usual objective is the evidence lower bound,
[ \mathcal{L}(\lambda)
\mathbb{E}{q\lambda} \left[ \log p(y,\theta)-\log q_\lambda(\theta) \right]. ]
This quantity differs from (\log p(y)) by the Kullback–Leibler divergence from the approximation to the posterior. Restrictions imposed by the variational family determine which forms of posterior dependence and tail behavior can be represented. Mean-field approximations impose independence among selected parameter blocks, thereby simplifying optimization while changing the covariance structure of the computed distribution.
Integrated nested Laplace approximations combine conditional Laplace approximations with numerical integration over a lower-dimensional set of hyperparameters. Their mathematical domain is latent Gaussian modeling, in which a Gaussian latent field is connected to observations through a non-Gaussian likelihood.
Diagnostics and numerical assessment
A computed posterior requires assessment of whether the numerical output represents its intended target. For Markov chains, serial dependence reduces the amount of information contained in a fixed number of iterations. If (\rho_k) denotes lag-(k) autocorrelation for a scalar estimand, an asymptotic effective sample size is
[ S_{\mathrm{eff}} \approx \frac{S} {1+2\sum_{k=1}^{\infty}\rho_k}. ]
This expression concerns a particular posterior function rather than the chain as an undifferentiated object. Two quantities computed from the same chain can therefore have different effective sample sizes.
The Gelman–Rubin statistic, developed by Andrew Gelman and Donald Rubin, compares variation within chains with variation between chains. Modern rank-based forms also examine differences in scale and tail behavior. Agreement among chains is a statement about the sampled regions and does not independently establish that all posterior modes have been visited.
Gradient-based samplers provide additional numerical information. Divergent numerical trajectories indicate that the finite-step integrator failed to follow the relevant Hamiltonian dynamics within the required tolerance. Low acceptance probabilities indicate substantial integration error or unsuitable proposal scales. Persistent transitions at the maximum trajectory length indicate that the numerical path was truncated by an algorithmic boundary.
For importance sampling, weight concentration is a central diagnostic because unstable weights produce estimates dominated by rare draws. Tail analysis of the weight distribution measures whether the proposal provides adequate coverage of the posterior. For variational methods, assessment instead concerns approximation bias, optimization convergence, and the representational limits of the chosen family.
Model evidence and predictive computation
The marginal likelihood
[ p(y)
\int p(y\mid\theta)p(\theta),d\theta ]
serves as a normalizing constant and as a component of Bayes factors. Its computation differs from ordinary posterior expectation because Markov chain output generally represents only the normalized posterior shape. Methods based on bridge identities, thermodynamic paths, or sequential distributions recover ratios of normalizing constants by connecting distributions whose normalization is known to those whose normalization is unknown.
Posterior predictive inference integrates over parameter uncertainty:
[ p(\widetilde{y}\mid y)
\int p(\widetilde{y}\mid\theta)p(\theta\mid y),d\theta. ]
A numerical representation of this distribution propagates posterior uncertainty into future or replicated observations. Posterior predictive checks compare features of observed data with the corresponding distribution under replicated data. Their interpretation concerns compatibility between the fitted model and selected data summaries rather than the probability that the model itself is true.
Predictive evaluation can also be based on cross-validation, where observations are evaluated under distributions fitted without those observations. Importance reweighting can approximate leave-one-out calculations from draws obtained under the full posterior, although the stability of this approximation again depends on the tails of the resulting weights.
Computational role in Bayesian modeling
Bayesian computation is part of the specification and interpretation of a statistical analysis rather than a detachable final stage. Parameterization changes posterior geometry without changing the underlying probability model. Centered and noncentered representations of a hierarchical model, for example, can produce the same marginal distribution while generating substantially different behavior for a numerical algorithm.
The relationship also operates in the other direction. A computational approximation can define a distribution that differs systematically from the formal posterior. Finite-particle methods introduce resampling variability. Restricted variational families alter dependence and dispersion. Inadequately explored Markov chains represent only the visited portion of the state space. Consequently, computed uncertainty contains both the probabilistic consequences of the model and the numerical consequences of its representation.