Principles of Monte Carlo Integration

Bayesian inference relies on understanding the posterior distribution $p(\theta | D)$. Often, the entire distribution itself is not the primary focus; instead, summaries derived from it are of interest. These summaries typically take the form of expectations. For instance, the posterior mean, a common point estimate for parameters $\theta$, is calculated as:

$$E[\theta | D] = \int \theta p(\theta | D) d\theta$$

Similarly, we might want to compute the posterior variance, credible intervals, or the probability of a specific hypothesis, all of which involve integrating some function $f(\theta)$ with respect to the posterior distribution:

$$E[f(\theta) | D] = \int f(\theta) p(\theta | D) d\theta$$

The core difficulty arises because the posterior $p(\theta | D)$ is frequently high-dimensional and lacks a standard analytical form. Direct computation of these integrals is usually impossible.

This is where the idea of Monte Carlo integration comes into play. It's a numerical method that approximates the value of a definite integral using random sampling. The intuition is grounded in the Law of Large Numbers (LLN).

The Monte Carlo Estimator

Suppose we want to estimate the expectation $E[f(\theta)] = \int f(\theta) p(\theta) d\theta$, where $\theta$ is drawn from the distribution $p(\theta)$. If we could somehow generate $S$ independent samples $\theta^{(1)}, \theta^{(2)}, ..., \theta^{(S)}$ directly from the distribution $p(\theta)$, the LLN tells us that the sample average of $f(\theta)$ will converge to the true expectation as the number of samples increases.

We define the Monte Carlo estimator of $E[f(\theta)]$ as:

$$\hat{E}_S[f(\theta)] = \frac{1}{S} \sum_{s=1}^{S} f(\theta^{(s)})$$

where $\theta^{(s)} \sim p(\theta)$ for $s=1, ..., S$.

As $S \to \infty$, this estimator converges to the true expectation:

$$\lim_{S \to \infty} \frac{1}{S} \sum_{s=1}^{S} f(\theta^{(s)}) = E[f(\theta)]$$

Furthermore, the Central Limit Theorem (CLT) often applies. If the variance of $f(\theta)$, denoted $\text{Var}[f(\theta)]$, is finite, then for large $S$, the distribution of the estimator $\hat{E}_S[f(\theta)]$ is approximately Normal:

$$\hat{E}_S[f(\theta)] \approx N\left(E[f(\theta)], \frac{\text{Var}[f(\theta)]}{S}\right)$$

This is immensely useful because it tells us how the accuracy of our estimate improves with more samples. The standard error of the estimate decreases proportionally to $1/\sqrt{S}$. We can estimate $\text{Var}[f(\theta)]$ using the sample variance of the computed $f(\theta^{(s)})$ values, giving us a way to quantify the uncertainty in our Monte Carlo estimate.

An Illustrative Example: Estimating a Beta Distribution Mean

Let's consider a simple example. Suppose our parameter $\theta$ follows a Beta distribution, $\theta \sim \text{Beta}(\alpha, \beta)$, perhaps representing the unknown probability of success in a series of Bernoulli trials. The probability density function is $p(\theta | \alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1} (1-\theta)^{\beta-1}$ for $\theta \in [0, 1]$.

We want to estimate the mean $E[\theta]$. Analytically, we know $E[\theta] = \frac{\alpha}{\alpha+\beta}$. Let's choose $\alpha=3$ and $\beta=7$. The true mean is $E[\theta] = 3 / (3+7) = 0.3$.

Now, let's pretend we don't know the analytical formula and use Monte Carlo integration. Here, $f(\theta) = \theta$. We can easily draw samples $\theta^{(s)}$ directly from a $\text{Beta}(3, 7)$ distribution using standard statistical software libraries.

We draw $S$ samples and compute the Monte Carlo estimate:

$$\hat{E}_S[\theta] = \frac{1}{S} \sum_{s=1}^{S} \theta^{(s)}$$

The plot below shows how the estimate $\hat{E}_S[\theta]$ converges towards the true mean (0.3) as the number of samples $S$ increases.

Convergence of the Monte Carlo estimate for the mean of a Beta(3, 7) distribution. The blue line shows the estimate based on $S$ samples (log scale), approaching the true mean (red dashed line) as $S$ increases.

This simple example illustrates the core principle: if you can sample directly from the target distribution $p(\theta)$, you can approximate expectations using sample averages.

The Catch: Sampling from the Posterior

The power of Monte Carlo integration relies entirely on our ability to draw independent samples $\theta^{(s)}$ from the target distribution $p(\theta)$. In Bayesian inference, our target distribution is the posterior $p(\theta | D)$.

Unfortunately, for most interesting models, the posterior distribution $p(\theta | D) \propto p(D | \theta) p(\theta)$ is a complex, high-dimensional function known only up to a normalizing constant (the evidence $p(D)$). Directly generating independent samples from such distributions is typically intractable. We cannot simply call a function like numpy.random.beta to get samples from an arbitrary posterior.

Connecting to Markov Chain Monte Carlo

This significant limitation is precisely why we need Markov Chain Monte Carlo (MCMC) methods. MCMC techniques provide an ingenious workaround. Instead of generating independent samples directly from $p(\theta | D)$, they construct a Markov chain whose states are possible values of $\theta$, and whose equilibrium (or stationary) distribution is exactly the target posterior $p(\theta | D)$.

By running this Markov chain for a sufficiently long time, the states visited by the chain can be treated as samples from the posterior distribution. While these samples are typically correlated (unlike the independent samples assumed in basic Monte Carlo), they can still be used within the Monte Carlo integration framework. The Law of Large Numbers still holds under certain conditions (ergodicity) for Markov chains.

So, MCMC methods allow us to generate the necessary (though correlated) samples $\theta^{(1)}, \theta^{(2)}, ..., \theta^{(S)}$ from distributions we cannot sample from directly. Once we have these samples, we use the standard Monte Carlo estimator $\frac{1}{S} \sum_{s=1}^{S} f(\theta^{(s)})$ to approximate posterior expectations $E[f(\theta) | D]$.

The subsequent sections of this chapter will detail how to build and use these Markov chains effectively, starting with the foundational Metropolis-Hastings algorithm.