Outside conjugate models the posterior distribution is known only up to its normalising constant, since the marginal likelihood f(x‾)=∫Θf(x‾∣θ)π(θ)dθf\left(\underline{x}\right)=\int_{\Theta}f\left(\underline{x}\mid\theta\right)\pi\left(\theta\right)d\theta has no closed form, and posterior summaries such as E[g(θ)∣x‾]=∫Θg(θ)π(θ∣x‾)dθ\mathbb{E}\left[g\left(\theta\right)\mid\underline{x}\right]=\int_{\Theta}g\left(\theta\right)\pi\left(\theta\mid\underline{x}\right)d\theta are intractable integrals over a parameter space whose dimension is often large. Numerical integration on a grid is only feasible in one or two dimensions, since the number of grid points grows exponentially with the dimension. Computational inference approximates these integrals either by simulation from the posterior or by approximating the posterior with a tractable distribution. Below π~(θ)=f(x‾∣θ)π(θ)\tilde{\pi}\left(\theta\right)=f\left(\underline{x}\mid\theta\right)\pi\left(\theta\right) is the unnormalised posterior, which can be evaluated, and π(θ∣x‾)=π~(θ)f(x‾)\pi\left(\theta\mid\underline{x}\right)=\frac{\tilde{\pi}\left(\theta\right)}{f\left(\underline{x}\right)} is the target.

Let θ(1),…,θ(S)\theta^{(1)},\ldots,\theta^{(S)} be iid draws from the posterior distribution. The Monte Carlo estimate of a posterior expectation is the average of the function over the draws.

g^S:=1S∑s=1Sg(θ(s))≈E[g(θ)∣x‾]\hat{g}_S:=\frac{1}{S}\sum_{s=1}^S g\left(\theta^{(s)}\right)\approx\mathbb{E}\left[g\left(\theta\right)\mid\underline{x}\right]

The estimate is unbiased, converges almost surely to the posterior expectation by the strong law of large numbers and, if g(θ)g\left(\theta\right) has finite posterior variance σg2,\sigma_g^2, is asymptotically normal by the central limit theorem with Monte Carlo standard error σgS,\frac{\sigma_g}{\sqrt{S}}, estimated by replacing σg\sigma_g with the sample standard deviation of the g(θ(s)).g\left(\theta^{(s)}\right). The error decreases as S−12S^{-\frac{1}{2}} regardless of the dimension of Θ.\Theta. Posterior probabilities are estimated by averaging indicators, quantiles by sample quantiles of the draws and the posterior predictive distribution by drawing x~(s)∼f(x~∣θ(s))\tilde{x}^{(s)}\sim f\left(\tilde{x}\mid\theta^{(s)}\right) for each s.s. The difficulty is generating the draws.

Rejection sampling draws from the target using a proposal distribution qq that is easy to sample from and a constant MM with π~(θ)≤Mq(θ)  ∀ θ.\tilde{\pi}\left(\theta\right)\leq Mq\left(\theta\right)\,\,\forall\,\theta. A draw θ∼q\theta\sim q is accepted with probability π~(θ)Mq(θ)\frac{\tilde{\pi}\left(\theta\right)}{Mq\left(\theta\right)} and otherwise discarded, and the accepted draws are iid from the posterior. The probability that a draw is accepted is f(x‾)M,\frac{f\left(\underline{x}\right)}{M}, which is typically tiny in more than a few dimensions because no simple qq matches the posterior closely enough for MM to be small.

Importance sampling reweights draws from a proposal distribution rather than discarding them. For a proposal qq with q(θ)>0q\left(\theta\right)>0 wherever π~(θ)>0,\tilde{\pi}\left(\theta\right)>0,

E[g(θ)∣x‾]=∫g(θ)w(θ)q(θ)dθ∫w(θ)q(θ)dθ\mathbb{E}\left[g\left(\theta\right)\mid\underline{x}\right]=\frac{\int g\left(\theta\right)w\left(\theta\right)q\left(\theta\right)d\theta}{\int w\left(\theta\right)q\left(\theta\right)d\theta} w(θ):=π~(θ)q(θ)w\left(\theta\right):=\frac{\tilde{\pi}\left(\theta\right)}{q\left(\theta\right)}

so with draws θ(s)∼q\theta^{(s)}\sim q and importance weights ws=w(θ(s))w_s=w\left(\theta^{(s)}\right) the self-normalised importance sampling estimate is

g^S=∑s=1Swsg(θ(s))∑s=1Sws,\hat{g}_S=\frac{\sum_{s=1}^S w_sg\left(\theta^{(s)}\right)}{\sum_{s=1}^S w_s},

which is consistent and needs only the unnormalised posterior, while 1S∑s=1Sws\frac{1}{S}\sum_{s=1}^S w_s is an unbiased estimate of the marginal likelihood. Its variance is governed by the variability of the weights and is summarised by the effective sample size (∑sws)2∑sws2∈[1,S],\frac{\left(\sum_s w_s\right)^2}{\sum_s w_s^2}\in\left[1,S\right], the number of iid posterior draws that would give the same precision. The variance is infinite if qq has lighter tails than the posterior, and the effective sample size collapses in high dimensions, so importance sampling is used for low-dimensional problems or to correct draws from an approximate posterior.

Markov chain Monte Carlo (MCMC) constructs an irreducible, aperiodic Markov chain on Θ\Theta whose stationary distribution is the posterior and runs it to obtain dependent draws θ(1),θ(2),…\theta^{(1)},\theta^{(2)},\ldots By the ergodic theorem for Markov chains the sample path average 1S∑s=1Sg(θ(s))\frac{1}{S}\sum_{s=1}^S g\left(\theta^{(s)}\right) converges almost surely to the posterior expectation from any starting point, and the distribution of θ(s)\theta^{(s)} converges to the posterior. The chain is constructed to satisfy the detailed balance equations with respect to π(θ∣x‾),\pi\left(\theta\mid\underline{x}\right), which only requires ratios of the target and so needs only π~.\tilde{\pi}. The initial draws before the chain has approached its stationary distribution are the burn-in and are discarded. Because successive draws are correlated, the Monte Carlo standard error is larger than for iid draws: with autocorrelation ρk\rho_k at lag kk in the sequence g(θ(s))g\left(\theta^{(s)}\right) the effective sample size is S1+2∑k=1∞ρk.\frac{S}{1+2\sum_{k=1}^{\infty}\rho_k}. Convergence is assessed by running several chains from dispersed starting points and comparing the variance between chains with the variance within chains through the potential scale reduction factor R^,\hat{R}, which approaches one when the chains have mixed.

The Metropolis–Hastings algorithm generates the next state from the current state θ\theta by drawing a candidate θ′\theta' from a proposal distribution q(θ′∣θ)q\left(\theta'\mid\theta\right) and accepting it with the acceptance probability

α(θ,θ′):=min⁡{1,π~(θ′)q(θ∣θ′)π~(θ)q(θ′∣θ)},\alpha\left(\theta,\theta'\right):=\min\left\{1,\frac{\tilde{\pi}\left(\theta'\right)q\left(\theta\mid\theta'\right)}{\tilde{\pi}\left(\theta\right)q\left(\theta'\mid\theta\right)}\right\},

setting θ(s+1)=θ′\theta^{(s+1)}=\theta' if accepted and θ(s+1)=θ\theta^{(s+1)}=\theta otherwise. The resulting chain satisfies detailed balance with respect to the posterior for any proposal that makes it irreducible. With a symmetric proposal such as θ′∼N(θ,τ2I)\theta'\sim\mathrm{N}\left(\theta,\tau^2I\right) the algorithm is random walk Metropolis and the acceptance probability reduces to min⁡{1,π~(θ′)π~(θ)},\min\left\{1,\frac{\tilde{\pi}\left(\theta'\right)}{\tilde{\pi}\left(\theta\right)}\right\}, so moves to higher posterior density are always accepted and moves to lower density are accepted with probability equal to the density ratio. The independence sampler proposes from a fixed distribution q(θ′)q\left(\theta'\right) that approximates the posterior. The proposal scale τ\tau trades off the size of moves against the acceptance rate; too small and the chain moves slowly, too large and almost every candidate is rejected, with an acceptance rate near 0.230.23 optimal for random walk proposals in high dimensions. Random walk Metropolis explores a dd-dimensional posterior in a number of steps that grows like d,d, which is slow for large d.d.

Let θ=(θ1,…,θd)\theta=\left(\theta_1,\ldots,\theta_d\right) and θ−j\theta_{-j} denote all components except θj.\theta_j. The full conditional distribution of θj\theta_j is its posterior distribution given the other components, π(θj∣θ−j,x‾)∝π~(θ)\pi\left(\theta_j\mid\theta_{-j},\underline{x}\right)\propto\tilde{\pi}\left(\theta\right) viewed as a function of θj\theta_j alone. The Gibbs sampler updates one component at a time by drawing it from its full conditional distribution with the other components fixed at their current values, cycling through j=1,…,dj=1,\ldots,d to produce each new state. It is the Metropolis–Hastings algorithm with the full conditional as proposal, for which every candidate is accepted, and it requires no tuning but does require the full conditionals to be standard distributions, which is the case when each parameter has a conjugate prior given the others. Strongly correlated components make the sampler move slowly since each update is restricted to a coordinate direction, and components that are not conjugate can be updated with a Metropolis–Hastings step inside the cycle. Introducing latent variables whose full conditionals, together with those of θ,\theta, are tractable is data augmentation and is how mixture and censored models are sampled.

Hamiltonian Monte Carlo (HMC) uses the gradient of the log posterior to propose distant states with high acceptance probability. The parameter is augmented with an auxiliary momentum vector p‾∈Rd\underline{p}\in\mathbb{R}^d with distribution Nd(0‾,M)\mathrm{N}_d\left(\underline{0},M\right) for a positive definite mass matrix M,M, and the Hamiltonian is the negative log joint density of (θ,p‾),\left(\theta,\underline{p}\right),

H(θ,p‾):=−log⁡eπ~(θ)+12p‾⊺M−1p‾,H\left(\theta,\underline{p}\right):=-\log_e\tilde{\pi}\left(\theta\right)+\frac{1}{2}\underline{p}^\intercal M^{-1}\underline{p},

the sum of a potential energy and a kinetic energy. Each iteration draws a fresh momentum, then simulates Hamilton’s equations dθdt=M−1p‾,\frac{d\theta}{dt}=M^{-1}\underline{p}, dp‾dt=∇θlog⁡eπ~(θ)\frac{d\underline{p}}{dt}=\nabla_{\theta}\log_e\tilde{\pi}\left(\theta\right) for a fixed time using the leapfrog integrator with step size ε\varepsilon and LL steps, which alternates half-step momentum updates with full-step position updates, and accepts the end point (θ′,p‾′)\left(\theta',\underline{p}'\right) with probability min⁡{1,exp⁡(H(θ,p‾)−H(θ′,p‾′))}.\min\left\{1,\exp\left(H\left(\theta,\underline{p}\right)-H\left(\theta',\underline{p}'\right)\right)\right\}. Exact Hamiltonian dynamics conserve HH and would always be accepted, and the leapfrog integrator is reversible and volume preserving so that the Metropolis correction exactly removes its discretisation error. The trajectory follows the geometry of the posterior rather than diffusing, so HMC explores a dd-dimensional posterior in a number of steps that grows like d14d^{\frac{1}{4}} and far outperforms random walk methods for large d.d. It requires a differentiable log posterior and tuning of ε\varepsilon and L;L; the No-U-Turn sampler (NUTS) chooses LL adaptively by extending the trajectory until it starts to double back, and is the default sampler in modern probabilistic programming software.

Analytic approximations replace the posterior by a tractable distribution rather than sampling from it. The Laplace approximation fits a normal distribution at the posterior mode using the second-order Taylor expansion of the log posterior,

π(θ∣x‾)≈Nd(θ^MAP,H−1)\pi\left(\theta\mid\underline{x}\right)\approx\mathrm{N}_d\left(\hat{\theta}_{\mathrm{MAP}},H^{-1}\right) H:=−∇θ2log⁡eπ~(θ^MAP)H:=-\nabla^2_{\theta}\log_e\tilde{\pi}\left(\hat{\theta}_{\mathrm{MAP}}\right)

where HH is the negative Hessian matrix of the log posterior at the mode, and gives the approximation of the marginal likelihood f(x‾)≈π~(θ^MAP)(2π)d2∣H∣−12f\left(\underline{x}\right)\approx\tilde{\pi}\left(\hat{\theta}_{\mathrm{MAP}}\right)\left(2\pi\right)^{\frac{d}{2}}\lvert H\rvert^{-\frac{1}{2}} by integrating the normal density. It is accurate when the posterior is close to normal, which by the Bernstein–von Mises theorem holds for large nn in regular models, and poor for skewed, multimodal or bounded posteriors.

Variational inference chooses a family Q\mathcal{Q} of tractable distributions and finds the member closest to the posterior in Kullback–Leibler divergence,

q∗:=arg min⁡q∈QKL(q ∥ π(⋅∣x‾)),q^*:=\operatorname*{arg\,min}_{q\in\mathcal{Q}}\mathrm{KL}\left(q\,\|\,\pi\left(\cdot\mid\underline{x}\right)\right),

which turns integration into optimisation. The divergence involves the unknown marginal likelihood, but expanding it gives

log⁡ef(x‾)=ELBO(q)+KL(q ∥ π(⋅∣x‾))\log_e f\left(\underline{x}\right)=\mathrm{ELBO}\left(q\right)+\mathrm{KL}\left(q\,\|\,\pi\left(\cdot\mid\underline{x}\right)\right) ELBO(q):=Eq[log⁡eπ~(θ)]−Eq[log⁡eq(θ)]\mathrm{ELBO}\left(q\right):=\mathbb{E}_q\left[\log_e\tilde{\pi}\left(\theta\right)\right]-\mathbb{E}_q\left[\log_e q\left(\theta\right)\right]

where the evidence lower bound (ELBO) depends only on the unnormalised posterior and, since the divergence is non-negative, is a lower bound on the log marginal likelihood that is tight when qq equals the posterior. Minimising the divergence is therefore equivalent to maximising the ELBO over Q.\mathcal{Q}. The mean-field family consists of distributions that factorise over components, q(θ)=∏j=1dqj(θj),q\left(\theta\right)=\prod_{j=1}^d q_j\left(\theta_j\right), for which coordinate ascent variational inference maximises the ELBO one factor at a time with the closed-form update

log⁡eqj∗(θj)=Eq−j[log⁡eπ~(θ)]+constant,\log_e q_j^*\left(\theta_j\right)=\mathbb{E}_{q_{-j}}\left[\log_e\tilde{\pi}\left(\theta\right)\right]+\mathrm{constant},

the expectation of the log joint density over the other factors, mirroring the Gibbs sampler with expectations in place of draws. Alternatively, Q\mathcal{Q} is a parametric family such as multivariate normal distributions and the ELBO is maximised by gradient methods using Monte Carlo estimates of its gradient. Variational inference is much faster than MCMC and scales to large data sets, but the fitted q∗q^* underestimates the posterior variance, because minimising KL(q ∥ π)\mathrm{KL}\left(q\,\|\,\pi\right) penalises placing mass of qq where the posterior is small more than it penalises missing mass of the posterior, and the mean-field family cannot represent posterior correlations.