Monte Carlo Integration

Learning Objectives

By the end of this section, you will be able to:

Where You'll Apply This: Computing expectations in high-dimensional spaces, training neural networks with stochastic gradient descent, uncertainty quantification, Bayesian inference, reinforcement learning, and any situation where exact integration is impossible.

The Big Picture

Many problems in statistics, machine learning, and physics require computing integrals that have no closed-form solution. Consider computing the expected value of a neural network's output over all possible inputs, or the posterior mean in a Bayesian model with thousands of parameters. These integrals are intractable - we cannot solve them analytically.

Monte Carlo integration provides a remarkably simple and powerful solution: instead of computing the integral exactly, we approximate it by drawing random samples. This idea - replacing integrals with random sampling - is one of the most important computational techniques in modern science and engineering.

Historical Development

The Fundamental Idea

The genius of Monte Carlo lies in its simplicity. To estimate an expected value$\mathbb{E}_p[f(X)]$, we:

1. Draw samples $X_1, X_2, \ldots, X_n$ from the distribution $p(x)$
2. Evaluate the function $f(X_i)$ at each sample point
3. Average the results: $\hat{I}_n = \frac{1}{n} \sum_{i=1}^{n} f(X_i)$

That's it! The sample average converges to the true expected value as the number of samples increases. This works for any distribution $p(x)$ we can sample from, and any function $f(x)$ we can evaluate.

Law of Large Numbers Connection

Why does Monte Carlo work? The answer is the Law of Large Numbers (LLN), which we covered in Chapter 10. The LLN guarantees that the sample average converges to the population mean:

Interactive: Basic Monte Carlo

Let's see Monte Carlo integration in action. We'll estimate the integral of a complex function using random samples. Watch how the estimate converges to the true value as you add more samples.

Convergence Analysis

Knowing that Monte Carlo converges isn't enough - we need to know how fastit converges and how accurate our estimates are. This is where the Central Limit Theorembecomes essential.

Variance and the CLT

By the CLT, the Monte Carlo estimator is approximately normally distributed:

This tells us several crucial things:

Interactive: Convergence Visualizer

Watch multiple independent Monte Carlo runs and observe how they all converge to the true value. The "confidence band" shrinks as $1/\sqrt{n}$, just as the theory predicts.

Why MC Beats Grids in High Dimensions

To understand Monte Carlo's power, compare it to deterministic grid-based integration (like the trapezoid rule):

Classic Example: Estimating π

The most famous demonstration of Monte Carlo integration is estimating π by throwing random "darts" at a square with an inscribed quarter circle. This example perfectly illustrates the core MC idea: probability ratios equal area ratios.

Interactive: π Estimator

Watch as random points accumulate and the estimate of π improves. Notice how the convergence follows the $1/\sqrt{n}$ law we derived above.

Variance Reduction Techniques

While basic Monte Carlo is powerful, we can often do much better by being clever abouthow we sample. Variance reduction techniques maintain the correct expected value while reducing the variance of our estimator - giving us the same accuracy with fewer samples.

Importance Sampling

The key observation is that when we estimate $\mathbb{E}_p[f(X)]$, not all regions of the input space contribute equally. If $f(x)$ varies greatly, we waste effort sampling where$f(x) \approx 0$ while missing regions where $f(x)$ is large.

The proposal distribution $q(x)$ should concentrate probability mass where $|f(x)| \cdot p(x)$ is large. The ideal choice (which achieves zero variance!) is:

Of course, we usually can't use the optimal proposal (computing it requires the integral we're trying to estimate!), but even approximate matches can dramatically reduce variance.

Interactive: Importance Sampling

Compare naive uniform sampling with importance sampling. Notice how importance sampling focuses samples where the integrand is large, achieving lower variance with the same sample count.

Monte Carlo in Bayesian Inference

One of the most important applications of Monte Carlo integration is in Bayesian inference. Given data, we want to compute quantities like:

• The posterior mean: $\mathbb{E}[\theta | \text{data}]$
• The posterior variance: $\text{Var}[\theta | \text{data}]$
• Posterior probabilities: $P(\theta > c | \text{data})$
• The predictive distribution: $p(y_{\text{new}} | \text{data})$

All of these are expectations with respect to the posterior distribution. Monte Carlo makes them trivial: just draw samples from the posterior and compute the sample average!

This is incredibly powerful: once you have posterior samples, you can estimate anyfunction of the parameters without recomputing anything. Need the mean? Average the samples. Need the variance? Compute sample variance. Need $P(\theta > 0.5)$? Count the fraction above 0.5.

Interactive: Bayesian Integration

See Monte Carlo in action for Bayesian inference. We compute the posterior mean for a Beta-Binomial model and compare the Monte Carlo estimate to the exact analytical solution.

Applications in AI/ML

Monte Carlo methods are ubiquitous in modern machine learning. Here are some key applications:

Python Implementation

Let's implement Monte Carlo integration from scratch. Click on code lines to see detailed explanations of each component.

np.random.uniform(0, 1, 10000) - draw 10K samples

For [0, 1]: each x has equal probability of being anywhere in the interval

For f(x)=x² on [0,1]: estimate = 1 × mean(x_i²) ≈ 1/3

ESS = n means all weights are equal (perfect proposal). ESS = 1 means one sample dominates.

1import numpy as np
2from scipy import stats
3
4# ================================================
5# Monte Carlo Integration: Core Implementation
6# ================================================
7
8def monte_carlo_integrate(f, a, b, n_samples=10000, seed=None):
9    """
10    Estimate integral of f(x) from a to b using Monte Carlo.
11
12    Returns: (estimate, standard_error)
13    """
14    if seed is not None:
15        np.random.seed(seed)
16
17    # Draw uniform samples from [a, b]
18    x_samples = np.random.uniform(a, b, n_samples)
19
20    # Evaluate function at sample points
21    f_values = f(x_samples)
22
23    # Monte Carlo estimate: (b-a) * mean(f)
24    estimate = (b - a) * np.mean(f_values)
25
26    # Standard error estimate
27    std_error = (b - a) * np.std(f_values) / np.sqrt(n_samples)
28
29    return estimate, std_error
30
31
32def importance_sampling(f, p_pdf, q_pdf, q_sampler, n_samples=10000):
33    """
34    Estimate E_p[f(X)] using importance sampling with proposal q.
35
36    Args:
37        f: Function to integrate
38        p_pdf: Target density p(x)
39        q_pdf: Proposal density q(x)
40        q_sampler: Function that draws samples from q
41
42    Returns: (estimate, effective_sample_size)
43    """
44    # Sample from proposal distribution
45    x_samples = q_sampler(n_samples)
46
47    # Compute importance weights: w(x) = p(x) / q(x)
48    weights = p_pdf(x_samples) / q_pdf(x_samples)
49
50    # Normalize weights (self-normalized importance sampling)
51    normalized_weights = weights / np.sum(weights)
52
53    # Weighted estimate
54    estimate = np.sum(normalized_weights * f(x_samples))
55
56    # Effective sample size (diagnostic for weight quality)
57    ess = 1 / np.sum(normalized_weights ** 2)
58
59    return estimate, ess
60
61
62# ================================================
63# Example: Posterior Mean via Monte Carlo
64# ================================================
65
66def bayesian_mc_example():
67    """
68    Compute E[theta | data] for Beta-Binomial model.
69    """
70    # Prior: Beta(2, 2), Data: 15 successes in 20 trials
71    prior_a, prior_b = 2, 2
72    successes, trials = 15, 20
73
74    # Posterior: Beta(prior_a + successes, prior_b + failures)
75    post_a = prior_a + successes
76    post_b = prior_b + (trials - successes)
77
78    # Exact posterior mean (for comparison)
79    exact_mean = post_a / (post_a + post_b)
80
81    # Monte Carlo estimate
82    n_samples = 10000
83    theta_samples = np.random.beta(post_a, post_b, n_samples)
84    mc_mean = np.mean(theta_samples)
85    mc_std = np.std(theta_samples) / np.sqrt(n_samples)
86
87    print(f"Posterior: Beta({post_a}, {post_b})")
88    print(f"Exact Mean:     {exact_mean:.6f}")
89    print(f"MC Estimate:    {mc_mean:.6f} ± {mc_std:.6f}")
90    print(f"MC Error:       {abs(mc_mean - exact_mean):.6f}")
91
92    # Also estimate posterior probability P(theta > 0.7 | data)
93    prob_gt_07 = np.mean(theta_samples > 0.7)
94    exact_prob = 1 - stats.beta.cdf(0.7, post_a, post_b)
95    print(f"\nP(θ > 0.7 | data):")
96    print(f"  MC Estimate:  {prob_gt_07:.4f}")
97    print(f"  Exact:        {exact_prob:.4f}")
98
99
100# Run the example
101if __name__ == "__main__":
102    # Basic integration example
103    f = lambda x: np.exp(-x**2) * (1 + np.sin(3*x)**2)
104    est, se = monte_carlo_integrate(f, -3, 3, n_samples=100000)
105    print(f"Integral estimate: {est:.4f} ± {se:.4f}")
106
107    print("\n" + "="*50 + "\n")
108    bayesian_mc_example()

Common Pitfalls

Knowledge Check

Test your understanding of Monte Carlo integration with this interactive quiz.

Summary

Looking Ahead: In the next section, we'll tackle a fundamental problem: what if we can't directly sample from our target distribution? Markov Chain Monte Carlo(MCMC) provides the answer, using carefully constructed random walks to generate samples from complex, high-dimensional distributions.