Monte Carlo Methods

Credible Intervals

In our discussion of Bayesian inference, we introduced posterior distributions for unknown parameters. The page on Markov chains developed models for sequential data and noted that Markov Chain Monte Carlo (MCMC) methods use the Markov property to sample from complex posterior distributions.

We begin with the problem of summarizing a posterior distribution. Since a full distribution is often complex (especially in high-dimensional parameter spaces) and difficult to interpret or communicate, we frequently summarize the uncertainty about a parameter using a credible interval. A credible interval is fundamentally different from a confidence interval in frequentist statistics. The former makes a direct probability statement about the parameter given the data, whereas the latter is a statement about the long-run frequency of the procedure. Nevertheless, they may sometimes numerically coincide.

Recall the credible interval. For a level \(\alpha \in (0,1)\), it is an interval \([l, u]\) with \(P(\theta \in [l, u] \mid \mathcal{D}) = 1 - \alpha\) under the posterior distribution \(p(\theta \mid \mathcal{D})\). For continuous posteriors, this condition does not uniquely determine the interval, so a common choice is the central credible interval, which places \(\alpha/2\) of the posterior probability in each tail.

If the cumulative distribution function (cdf) \(F\) of the posterior distribution is known, the central credible interval can be computed as: \[ \begin{align*} l &= F^{-1}\left( \frac{\alpha}{2} \right), \\\\ u &= F^{-1}\left( 1 - \frac{\alpha}{2} \right), \end{align*} \] where \(F^{-1}\) denotes the quantile function (the inverse of the cdf).

As an example, if the posterior distribution is approximately normal with posterior mean \(\mu\) and posterior standard deviation \(\sigma\), a rough \(95\%\) credible interval can be set as: \[ [l, u] = [\mu - 2\sigma, \mu + 2\sigma], \] because for a standard normal distribution \(\mathcal{N}(0,1)\), approximately \(95\%\) (more precisely \(95.45\%\)) of the probability mass lies within two standard deviations from the mean.

More precisely, when \(p(\theta \mid \mathcal{D}) = \mathcal{N}(\mu, \sigma^2)\), the central credible interval endpoints are: \[ \begin{align*} l &= \mu + \sigma \Phi^{-1}\left( \frac{\alpha}{2} \right), \\\\ u &= \mu + \sigma \Phi^{-1}\left( 1 - \frac{\alpha}{2} \right), \end{align*} \] where \(\Phi\) is the cdf of the standard normal distribution.

In particular, if the posterior distribution is the standard normal \(p(\theta \mid \mathcal{D}) = \mathcal{N}(0,1)\) and \(\alpha = 0.05\), then: \[ \begin{align*} l &= \Phi^{-1}(0.025) \approx -1.96, \\\\ u &= \Phi^{-1}(0.975) \approx 1.96, \end{align*} \] so the \(95\%\) credible interval is approximately \([-1.96, 1.96]\).

The computation above relies on knowing the quantile function \(F^{-1}\) in closed form. For many posterior distributions, especially those arising from non-conjugate models or high-dimensional problems, no such closed form is available. This motivates the use of Monte Carlo methods.

Monte Carlo Approximation

The Monte Carlo method is a general computational technique that uses random sampling to approximate quantities that are difficult or impossible to compute analytically. The fundamental idea is simple. To estimate an expectation, we draw independent samples and compute their average.

By choosing different functions \(f\), this estimator approximates many functionals of the posterior. Setting \(f(\theta) = \theta\) gives the posterior mean. For the posterior variance, the standard plug-in estimator is the sample variance \(\hat{\sigma}^2_S = \frac{1}{S} \sum_{s=1}^{S} (\theta^{(s)} - \hat{\theta}_S)^2\), where \(\hat{\theta}_S = \frac{1}{S} \sum_s \theta^{(s)}\) is the sample mean. Strictly speaking, this is a plug-in estimator rather than a direct Monte Carlo estimator, since \(f\) depends on the sample, but it inherits consistency from the joint convergence of \(\hat{\theta}_S\) and the second-moment estimator. Predictive distributions, posterior probabilities of events, and other integrable functionals follow the same template.

When the inverse cdf of the posterior distribution is not available in closed form, a practical alternative is to estimate the quantiles from samples. Suppose we have \(S\) independent samples \(\theta^{(1)}, \theta^{(2)}, \ldots, \theta^{(S)}\) from the posterior \(p(\theta \mid \mathcal{D})\). To approximate a \(100(1-\alpha)\%\) central credible interval:

• Sort the samples in increasing order and write \(\theta_{(1)} \leq \theta_{(2)} \leq \cdots \leq \theta_{(S)}\) for the order statistics.
• Set the lower endpoint \(l\) to be the \(\left(\frac{\alpha}{2}\right)\)-quantile, and the upper endpoint \(u\) to be the \(\left(1 - \frac{\alpha}{2}\right)\)-quantile.

More precisely, one common choice (corresponding to the lower-bound ceiling convention) is: \[ \begin{align*} l &= \theta_{\left(\left\lceil S \times \frac{\alpha}{2} \right\rceil\right)}, \\\\ u &= \theta_{\left(\left\lceil S \times \left(1 - \frac{\alpha}{2}\right) \right\rceil\right)}, \end{align*} \] where \(\lceil \cdot \rceil\) denotes the ceiling function. Other conventions exist. For example, numerical libraries (such as numpy's quantile) often use linear interpolation between adjacent order statistics. Provided the target quantiles are uniquely determined, all such conventions agree in the limit \(S \to \infty\).

The argument proceeds in two steps. First, for each fixed \(x\), the indicator variables \(\mathbb{1}\{X_s \leq x\}\) are iid with mean \(F(x)\), so the strong law of large numbers gives pointwise convergence \(\hat{F}_S(x) \to F(x)\) almost surely. Second, since both \(\hat{F}_S\) and \(F\) are monotone non-decreasing functions taking values in \([0, 1]\), a partitioning-by-rationals argument upgrades the countable pointwise convergence to uniform convergence.

Remark. The partitioning argument above goes through directly when \(F\) is continuous. The general case (where \(F\) may have atoms) requires an additional argument handling the jump points separately.

We return to the credible interval problem. Instead of computing the inverse cdf analytically, we use the empirical cumulative distribution function (ECDF) formed by the sorted samples to estimate the desired quantiles. By the Glivenko-Cantelli theorem above, the ECDF converges uniformly to the true cdf. Provided the target quantiles are uniquely determined (\(F\) is not flat at the levels \(\alpha/2\) and \(1 - \alpha/2\)), this implies convergence of the empirical quantiles to the true quantiles.

Highest Posterior Density (HPD)

As discussed above, central credible intervals can be suboptimal for asymmetric or multimodal posteriors because they allocate equal probability mass to each tail regardless of where the density is concentrated. The highest posterior density (HPD) region provides a more principled alternative by selecting the smallest region that contains the desired probability mass.

In one dimension, the HPD region corresponds to the interval (or union of intervals) containing the values of highest posterior density, with total mass equal to the target credibility level (for example, 95%). In this case, HPD regions are sometimes called highest density intervals (HDIs).

Unlike central credible intervals, HPD regions are not required to place probability \(\alpha/2\) in each tail. This makes them especially useful when the posterior distribution is skewed or multimodal. In such cases, HPD regions better capture the actual regions of highest posterior density. For unimodal distributions, the HPD region is in fact a single interval, the shortest one containing \(1-\alpha\) probability mass.

The demo above shows that the HPD interval is often narrower than the central credible interval and may exclude low-density regions that the central interval mistakenly includes. Estimating either the central credible interval or the HPD region by Monte Carlo requires samples from the posterior distribution. In the examples above, we assumed that such samples are readily available. In practice, however, directly sampling from the posterior is often impossible. The next section introduces the algorithmic framework that makes posterior sampling feasible.

Markov Chain Monte Carlo (MCMC)

The Monte Carlo methods described above assume that we can draw independent samples directly from the posterior distribution \(p(\theta \mid \mathcal{D})\). In practice, however, the full posterior is rarely available in closed form. By Bayes' theorem, stated for events on the linked page and used here in its density form, we have \[ p(\theta \mid \mathcal{D}) = \frac{p(\mathcal{D} \mid \theta) \, p(\theta)}{p(\mathcal{D})}, \] where computing the normalizing constant (the marginal likelihood or evidence) \[ p(\mathcal{D}) = \int p(\mathcal{D} \mid \theta) \, p(\theta) \, d\theta \] requires integrating over the entire parameter space, which is often intractable in high dimensions.

We can, however, easily evaluate the unnormalized posterior \[ \tilde{p}(\theta) := p(\mathcal{D} \mid \theta) \cdot p(\theta) \propto p(\theta \mid \mathcal{D}), \] since this only requires evaluating the likelihood and prior at a given point \(\theta\).

This is where Markov Chain Monte Carlo (MCMC) methods become essential. The key insight is that we can construct a Markov chain whose stationary distribution equals the target posterior \(p(\theta \mid \mathcal{D})\), using only the unnormalized density \(\tilde{p}(\theta)\). Although the samples produced by such a chain are not independent, under an ergodicity condition on the chain, which we do not develop here, time averages along the chain converge to expectations under the stationary distribution. This enables Monte Carlo estimation without computing \(p(\mathcal{D})\).

Three classic MCMC algorithms are foundational to Bayesian computation:

• Metropolis-Hastings:
  A general-purpose algorithm that proposes new samples from a proposal distribution and accepts or rejects them based on a ratio of unnormalized target densities, multiplied by a ratio of proposal densities when the proposal is not symmetric.
• Gibbs sampling:
  A special case of Metropolis-Hastings that is efficient when the full conditional distributions are available in closed form. It is commonly used in probabilistic graphical models.
• Hamiltonian Monte Carlo (HMC):
  A gradient-based method that uses the geometry of the posterior to propose distant yet high-probability samples, making it well-suited for high-dimensional problems. It is widely used in probabilistic programming frameworks such as Stan and PyMC.

MCMC occupies a specific place in the broader landscape of machine learning. In large-scale models such as transformers and large language models (LLMs), MCMC is not used for training due to its computational cost and poor scalability in very high-dimensional parameter spaces. These models are instead trained using gradient-based optimization (for example, stochastic gradient descent or Adam) with maximum likelihood or related objectives. However, MCMC remains widely used for Bayesian modeling in scientific applications, uncertainty quantification, simulation-based inference, and decision-making under uncertainty. These are settings where principled uncertainty estimates are essential.