Integration: Simple, Complex, Intractable

Three Integrands, Three Levels of Difficulty

Randomness for a Deterministic Problem

The objectives are deterministic. We can induce randomness to achieve them.

A Brief History of Monte Carlo Methods

Origins.

• Named after Monaco’s Monte Carlo Casino
• Developed during Manhattan Project (1940s)
• Pioneers: Stanislaw Ulam, John von Neumann, Nicholas Metropolis

Statistical Applications

The Posterior Mean Is an Integral

Bayesian Computation. Want:

\[\hat{\theta} = \E[\theta|X] = \int \theta \, \pi(\theta|X) \, d\theta\]

Even if we know \(\pi(\theta|X)\), this integral may be impossible.

A Monte Carlo Estimate of an Intractable Integral

Implementation Considerations

Critical Numerical Issues.

Example: Computing Probabilities. Instead of \(\prod_{i=1}^n p_i\), compute \(\exp\left(\sum_{i=1}^n \log p_i\right)\).

The Log-Sum-Exp Trick

Problem. Compute \(y = \log\left(\sum_{i=1}^n e^{x_i}\right)\) when \(x_i\) are very large or small.

Log-Sum-Exp. For any constant \(c\):

\[y = \log\sum_i e^{x_i} = \log\left(e^{c}\sum_i e^{x_i - c}\right) = c + \log\sum_i e^{x_i - c}\]

Log-Sum-Exp: Choosing the Shift

• No overflow: \(\max_i e^{x_i - c} = e^0 = 1\)
• Minimal underflow impact

Application. Essential for stable computation in mixture models, HMMs, neural networks.

Random Variables: Foundation of Simulation

Random Variable. A variable whose outcomes are numerical values of a random phenomenon.

To simulate any distribution, we need to understand its CDF and relationships to other distributions.

Discrete and Continuous CDFs

Probability Mass and Density Functions

Computational Note. For simulation we usually need \(F^{-1}\).

Expected Value and Variance

Expected Value.

Variance.

\[\Var(X) = \E[(X - \E[X])^2] = \E[X^2] - (\E[X])^2\]

Computing a Variance Stably

The two expressions for \(\Var(X)\) are algebraically identical. Numerically they are not.

\(\E[X^2] - (\E[X])^2\) subtracts two large, nearly equal numbers whenever \((\E[X])^2 \gg \Var(X)\). Recall the catastrophic cancellations we talked about in the log-sum-exp trick? Compute sample variances from the centered form \(\frac{1}{n-1}\sum_i (x_i - \bar{x})^2\).

Two-pass algorithm. First pass computes \(\bar{x}\); second pass accumulates \(\sum_i(x_i - \bar{x})^2\). Stable, but requires storing all \(n\) values and reading them twice.

Key Discrete Distributions

Key Relationships.

• \(\Bin(n,p) \stackrel{d}{=} \sum_{i=1}^n \Ber(p)\) with the \(\Ber(p)\) independent
• \(\Bin(n,p) \approx \Pois(np)\) when \(n\) is large and \(p\) small with \(np\) fixed

Key Continuous Distributions

A Gallery of Densities and Mass Functions

The Weak Law of Large Numbers

Weak Law of Large Numbers (WLLN). If \(X_1, X_2, \ldots, X_n\) are i.i.d. with \(\E|X| < \infty\), then:

\[\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i \xrightarrow{P} \E[X] \text{ as } n \to \infty\]

Convergence in Probability: For any \(\epsilon > 0\):

\[\lim_{n \to \infty} P(|\bar{X}_n - \E[X]| > \epsilon) = 0\]

Essentially, this says that for large \(n\), large deviations of the sample mean from the true mean become highly unlikely.

The Strong Law of Large Numbers

Strong Law of Large Numbers (SLLN). Under the same conditions:

\[\bar{X}_n \xrightarrow{a.s.} \E[X] \text{ as } n \to \infty\]

Almost Sure Convergence:

\[P\!\left(\lim_{n \to \infty} \bar{X}_n = \E[X]\right) = 1\]

We can think of this as saying that the sample mean converges to the true mean for “almost every” sequence of random draws.

Both results are asymptotic: they guarantee convergence as \(n \to \infty\), not for any fixed \(n\).

WLLN vs SLLN: Why the Distinction Matters

Both laws say \(\bar{X}_n \to \E[X]\), but neither says how fast. The CLT does.

The Central Limit Theorem

Central Limit Theorem (CLT). If \(X_1, \ldots, X_n\) are i.i.d. with \(\E[X] = \mu\) and \(\Var(X) = \sigma^2 < \infty\), then:

\[\sqrt{n}(\bar{X}_n - \mu) \xrightarrow{d} \Normal(0, \sigma^2)\]

From CLT to Error Bounds

\[ \begin{aligned} \text{From CLT: } & \bar{X}_n \approx \Normal\left(\mu, \frac{\sigma^2}{n}\right) \text{ for large } n\\ \text{Standard Error: } & \text{SE}(\bar{X}_n) = \frac{\sigma}{\sqrt{n}} \\ \text{95% CI: } & \bar{X}_n \pm 1.96 \cdot \frac{\sigma}{\sqrt{n}} \\ \text{Error: } & |\bar{X}_n - \mu| = O_P(n^{-1/2}) \end{aligned} \]

Three Asymptotic Results, Three Roles

Monte Carlo Integration

Definition. Monte Carlo integration is a method to approximate an integral using random sampling:

\[\theta = \int_a^b h(x)\,dx \approx \frac{b-a}{n}\sum_{i=1}^n h(U_i), \quad U_i \sim \Uniform(a,b)\]

Replace a deterministic sum with a random average, then apply the LLN.

Why use it.

• Analytical integration is hard or impossible
• High-dimensional integrals make grid methods infeasible
• Monte Carlo error rate: \(O_P(n^{-1/2})\) regardless of dimension

From Integrals to Expectations

The Transformation. For any integral \(\theta = \int_0^1 h(x)dx\):

\[ \begin{aligned} \theta &= \int_0^1 h(x)dx \\ &= \int_0^1 h(x) \cdot 1 \, dx \\ &= \int_0^1 h(x) \cdot f_U(x) \, dx \quad \text{where } f_U(x) = 1 \text{ for } x \in [0,1] \\ &= \E[h(U)] \quad \text{where } U \sim \Uniform(0,1) \end{aligned} \]

Any Integral Is an Expectation

• We’ve transformed an integral into an expectation
• Now we can apply SLLN to approximate it
• Sample averages \(\xrightarrow{a.s.}\) True integral

Works for any distribution. This is the foundation of Monte Carlo methods.

Why SLLN Provides Stronger Assurance

Monte Carlo Setting. We run one simulation and increase \(n\), we observe one sequence \(\hat{\theta}_1, \hat{\theta}_2, \ldots\), and we want this sequence to converge to the truth.

SLLN gives stronger assurance, CLT gives practical error bounds.

One Sample Path, Converging

Monte Carlo Error Falls at Rate \(n^{-1/2}\)

Shrinking bands come from CLT, not SLLN. SLLN ensures eventual convergence, CLT quantifies error for finite \(n\).

What You Actually Do in Practice

The theory describes a sequence \(\hat{\theta}_1, \hat{\theta}_2, \dots\) converging to \(\theta\). In practice you never watch that sequence. You fix \(n\), draw one sample \(X_1, \dots, X_n\), compute one number \(\hat{\theta}_n\), and report it. That’s the recipe for Monte Carlo estimation: one sample, one estimate, one confidence interval.

The roles of the two results.

• SLLN certifies this recipe with explicit guarantees: the estimator targets the right quantity, and you can get arbitrarily close by increasing \(n\). It says nothing about the \(n\) you chose.
• CLT enters at the \(n\) you chose: \(\hat{\theta}_n \approx \Normal(\theta,\, \sigma^2/n)\), so \(\hat{\theta}_n \pm 1.96\,\hat{\sigma}/\sqrt{n}\) is a 95% confidence interval. This is what you actually report.

SLLN tells you the method works. CLT tells you how well it worked at the \(n\) you used.

General Monte Carlo Principle

Monte Carlo Approximation. For any function \(g\) and random variable \(X\) with \(\E|g(X)| < \infty\):

\[\E[g(X)] \approx \frac{1}{n}\sum_{i=1}^n g(X_i)\]

where \(X_1, \ldots, X_n \iid X\).

Applications.

• Integration: \(\int_a^b h(x)dx = (b-a)\E[h(U)]\), \(U \sim \Uniform(a,b)\)
• Probability: \(P[X \in A] = \E[\mathbf{1}_A(X)]\)

If in addition \(\Var[g(X)] < \infty\), the CLT gives \(\left|\frac{1}{n}\sum_i g(X_i) - \E[g(X)]\right| = O_P(n^{-1/2})\), regardless of dimension.

Halving the Error Costs Four Times the Samples

• Error decreases as \(1/\sqrt{n}\) (not \(1/n\))
• To halve the error, need four times as many samples
• Dimension-independent: Same rate for \(\int_0^1\) or \(\int_{[0,1]^{100}}\)

The Curse of Dimensionality in Numerical Integration

Traditional Quadrature Methods. A product rule with \(m\) points per axis costs \(N = m^d\) evaluations:

Exponential Explosion

For \(m = 10\) points per dimension:

• 3D: \(10^3 = 1{,}000\) evaluations
• 10D: \(10^{10} = 10\) billion evaluations
• 100D: \(10^{100}\) evaluations, more than the number of atoms in the observable universe

Traditional methods become computationally impossible.

Quadrature Error Against the Monte Carlo Rate

Product midpoint rule on \(\prod_i e^{x_i}\) over \([0,1]^d\) (exact value \((e-1)^d\)): relative error \(O(N^{-2/d})\) in the total evaluations \(N\). Plain Monte Carlo: \(O_P(N^{-1/2})\) in every \(d\), needing only finite variance. The exponents coincide at \(d = 4\); quadrature wins below it, Monte Carlo above it, and the gap widens with \(d\). Constants keep the midpoint rule ahead at \(d = 4\) here; at \(d = 8\) the rate takes over.

Monte Carlo: Breaking the Curse

Dimension-Independent Convergence. For any dimension \(d\), provided \(\Var[h(\mathbf{U})] = \sigma^2 < \infty\), the Monte Carlo estimator for \(\theta = \int_{[a,b]^d} h(\mathbf{x})d\mathbf{x}\):

\[\hat{\theta}_n = \frac{(b-a)^d}{n}\sum_{i=1}^n h(\mathbf{U}_i), \quad \mathbf{U}_i \sim \Uniform([a,b]^d)\]

has error \(|\hat{\theta}_n - \theta| = O_P(n^{-1/2})\) regardless of \(d\).

The Rate Does Not Depend on the Dimension

Punchlines, and What Comes Next

• An integral you cannot do becomes an expectation you can average, as long as you can sample from a density covering the integrand.
• The WLLN tells you that you targeted the right object with the right statistic.
• The SLLN says that your statistic converges almost surely to the target.
• When the variance is finite, the CLT budgets the approximation cost (at \(n^{-1/2}\)) and hands you a confidence interval with it.
• Furthermore, that rate is the same in one dimension and in a million!

The Monte Carlo method, together with the LLNs and the CLT, hands us a universal recipe for approximating integrals and expectations. Without the last two, the statistician in us would be left unsatisfied, not ever knowing how close we are to the truth. With them, we can quantify the uncertainty in our approximation and report it.