Prev
Unit 8 · Numerical Differentiation and Integration
Monte Carlo Integration
55m100 XP
Next
∫Numerical Differentiation and Integration

Monte Carlo Integration

Expectations as integrals, 1/√N convergence, the curse of dimensionality, uncertainty propagation

Los Angeles and Los Alamos, 1946
Stanisław Ulam, a Polish mathematician working at Los Alamos, is recovering from a serious illness. To pass the time he plays Canfield solitaire, and he begins to wonder: what is the chance that a randomly dealt game can be won? He tries to count the possibilities combinatorially and gets nowhere. Then he has the idea that changes scientific computing. It would be much simpler, he realises, to deal out a hundred games, play them, and count the wins. And the same trick could answer questions at Los Alamos that no formula could — how neutrons scatter, split and multiply through layers of material, each one following a random path. He describes the idea to John von Neumann, who sees at once how to do it on the new electronic computers; by 1948 ENIAC is running neutron-transport calculations by random sampling. Their colleague Nicholas Metropolis names the method after the casino in Monaco where Ulam's uncle liked to gamble. In 1949 Metropolis and Ulam publish The Monte Carlo Method.
Any integral can be written as an expected value. If X\mathbf X is uniformly distributed over a region Ω\Omega of volume VV, then
∫Ωf(x) dx=V E[f(X)]≈VN∑i=1Nf(Xi).\int_\Omega f(\mathbf x)\,d\mathbf x = V\,\mathbb E[f(\mathbf X)] \approx \frac{V}{N}\sum_{i=1}^{N}f(\mathbf X_i).
Draw NN random points, average the function, multiply by the volume. More generally, if X\mathbf X has density pp, E[g(X)]=∫g p\mathbb E[g(\mathbf X)] = \int g\,p: the mean, variance or failure probability of any quantity computed from random inputs is an integral, and sampling estimates it.
💡The 1/√N law
By the central limit theorem, the error of the sample mean is approximately normal with standard deviation σN,σ2=Var⁡f(X).\frac{\sigma}{\sqrt N}, \qquad \sigma^2 = \operatorname{Var} f(\mathbf X). Two consequences. It is slow: a hundred times more samples for one more digit — estimating π from random points, N=106N = 10^6 gives an error of about 10−310^{-3}, where Simpson's rule on a 1-D integral would give fourteen digits. It does not depend on the dimension: dd never appears. A product rule with mm points per axis needs mdm^d evaluations and its accuracy per evaluation degrades as dd grows — the curse of dimensionality. In ten dimensions, a grid with only 10 points per axis already needs 101010^{10} evaluations; Monte Carlo estimates the volume of the 10-dimensional unit ball to 2% with 10610^6. And the error bar comes free: estimate σ\sigma from the samples themselves.
Derivation
Where σ/N\sigma/\sqrt N comes from
Let Yi=f(Xi)Y_i = f(\mathbf X_i) for independent random points, each with mean μ\mu (the integral we want, divided by VV) and variance σ2\sigma^2. The estimate is their average, Yˉ=1N∑iYi\bar Y = \tfrac1N\sum_iY_i.
It is right on average: E[Yˉ]=1N∑iE[Yi]=μE[\bar Y] = \tfrac1N\sum_iE[Y_i] = \mu.
Its spread. Variances of independent quantities add, and scaling by 1N\tfrac1N scales a variance by 1N2\tfrac1{N^2}:
Var⁡(Yˉ)=1N2∑i=1NVar⁡(Yi)=Nσ2N2=σ2N⟹std(Yˉ)=σN.\operatorname{Var}(\bar Y) = \frac{1}{N^2}\sum_{i=1}^N\operatorname{Var}(Y_i) = \frac{N\sigma^2}{N^2} = \frac{\sigma^2}{N} \qquad\Longrightarrow\qquad \text{std}(\bar Y) = \frac{\sigma}{\sqrt N}.
Nothing in this calculation mentions the dimension of X\mathbf X. The central limit theorem adds that Yˉ\bar Y is approximately normally distributed, so about 95% of runs land within 1.96 σ/N1.96\,\sigma/\sqrt N of the truth.
For π. Each point contributes Y=4Y = 4 if it lands inside the quarter circle (probability p=π/4p = \pi/4) and 0 otherwise. Its variance is 16 p(1−p)16\,p(1 - p), so σ=4p(1−p)=40.7854×0.2146=1.642\sigma = 4\sqrt{p(1 - p)} = 4\sqrt{0.7854 \times 0.2146} = 1.642 — the constant printed by the lab.
📐Uncertainty propagation
A cantilever's tip deflection is δ=PL3/(3EI)\delta = PL^3/(3EI), nominally 10.00 mm. The load is uncertain by 10%, the length by 0.5%, the modulus by 5%. Sample 200,000 sets of inputs, compute δ\delta for each: mean 10.01 mm, standard deviation 1.13 mm, 99th percentile 12.76 mm, and the probability of exceeding the serviceability limit L/250=12L/250 = 12 mm is 4.3%. A first-order (linearised) analysis gives the same standard deviation, 1.13 mm — but only sampling gives the tail probability directly, and only sampling still works when the model is a nonlinear finite-element run with no formula at all.
Derivation
The first-order estimate, step by step
For small independent input errors, linearise the model: δy≈∑i∂y∂xi δxi\delta y \approx \sum_i \frac{\partial y}{\partial x_i}\,\delta x_i. Independent variances add (as in the derivation above), each multiplied by the square of its sensitivity:
σy2≈∑i(∂y∂xi)2σi2.\sigma_y^2 \approx \sum_i\Big(\frac{\partial y}{\partial x_i}\Big)^2\sigma_i^2.
For a product of powers like δ=PL3/(3EI)\delta = PL^3/(3EI) this is neatest in relative terms. Taking logarithms, ln⁡δ=ln⁡P+3ln⁡L−ln⁡E+const\ln\delta = \ln P + 3\ln L - \ln E + \text{const}, so a relative change dδ/δ=dP/P+3 dL/L−dE/E\mathrm d\delta/\delta = \mathrm dP/P + 3\,\mathrm dL/L - \mathrm dE/E. The relative standard deviations then combine as
(σδδ)2=(0.10)2+(3×0.005)2+(0.05)2=0.012725,σδδ=0.1128.\Big(\frac{\sigma_\delta}{\delta}\Big)^2 = (0.10)^2 + (3 \times 0.005)^2 + (0.05)^2 = 0.012725, \qquad \frac{\sigma_\delta}{\delta} = 0.1128.
On the nominal 10.00 mm that is 1.128 mm — the lab's first-order value. The length is the most precise input yet still contributes, because it enters cubed. Lesson 12-2 writes the same rule for vectors as JPJTJPJ^{\mathsf T}.
The mathematicians know a great deal about very little and the physicists very little about a great deal.
– Stanisław Ulam, On the Ergodic Behavior of Dynamical Systems, Los Alamos report LA-2055 (1955)
Variance reduction. Since the error is σ/N\sigma/\sqrt N, reducing σ\sigma is worth as much as increasing NN:
  • Antithetic variates: pair each sample X\mathbf X with its mirror image; errors partly cancel for monotone models.
  • Stratified and Latin hypercube sampling: force samples to cover every part of each input's range evenly. Standard in engineering reliability studies.
  • Importance sampling: sample more often where the integrand is large — near the failure region when estimating a 10−610^{-6} failure probability — and correct with weights.
  • Quasi-Monte Carlo: replace random points with low-discrepancy sequences (Sobol, Halton) that fill space more evenly, approaching 1/N1/N convergence for smooth integrands.
And all of it depends on the quality of the random numbers themselves — a subject for Unit 12.
💡Where you will use this
Geotechnical and Structural reliability: probabilities of failure with uncertain soil strength, loads and material properties. Orbital Mechanics: collision probability and dispersion of re-entry footprints. Rocket Propulsion and Guidance: Monte Carlo dispersion analysis of trajectories — every launch-vehicle flight is cleared with thousands of simulated flights. GNSS and Navigation: particle filters represent probability distributions by samples.
📚Further reading
  • Art Owen, Monte Carlo Theory, Methods and Examples (2013) — a free online book on sampling, variance reduction and quasi-Monte Carlo, by one of the field's leading researchers.
  • Nicholas Metropolis and Stanisław Ulam, 'The Monte Carlo method', Journal of the American Statistical Association 44(247), 335–341 (1949) — the paper that named it.
  • Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 13, 'Random Numbers and Stochastic Simulation'.
↓But there is a problem
Everything in this unit computed integrals of given functions. The most important integrals in engineering are not of functions we know, but of laws we know: Newton's second law, the rate equations of a combustion chamber, the attitude dynamics of a spacecraft. A differential equation tells us how the state changes at every instant; integrating it forward in time gives the trajectory. In 1768 Leonhard Euler proposed the simplest possible way to do it: take the rate of change now, and step forward a little. Most of what comes after — Runge, Kutta, Adams, Gear and every ODE solver in SciPy — is a story of making that step smarter.