Monte Carlo integration · Probability & Statistics

The last two lessons produced samples. This one uses them — to compute integrals that have no closed form and that grid-based methods can't reach.

The idea

Any integral is an expectation in disguise:

$\int_a^b h(x)\,dx = (b-a)\int_a^b h(x)\cdot\frac{1}{b-a}\,dx = (b-a)\,E\big[h(U)\big]$

where U \sim \text{Uniform}(a,b). And by the Law of Large Numbers (§3), an expectation is estimated by an average:

$\hat I = \frac{b-a}{n}\sum_{i=1}^n h(U_i) \;\xrightarrow{a.s.}\; \int_a^b h(x)\,dx$

Draw points, evaluate, average, scale by the volume. That's the whole method.

More generally, for any density p you can sample from:

$\int h(x)p(x)\,dx = E_p[h(X)] \approx \frac{1}{n}\sum_{i=1}^n h(X_i)$

The error

The estimator is unbiased, and the CLT (§3) gives its distribution directly:

$\operatorname{SE}(\hat I) = \frac{\sigma}{\sqrt n}, \qquad \sigma^2 = \operatorname{Var}\big(h(U)\big)$

with \sigma estimable from the same sample. So you get a confidence interval for free:

$\hat I \pm 1.96\frac{\hat\sigma}{\sqrt n}$

That's a genuine advantage over deterministic quadrature, which gives you a number with no accompanying error estimate unless you do extra work.

Why it wins in high dimensions

The convergence rate O(n^{-1/2}) looks bad. The trapezoid rule is O(n^{-2}) and Simpson's is O(n^{-4}) in one dimension — far faster.

But in d dimensions, a grid method with n total points has only n^{1/d} points per axis, so its error becomes

$O\big(n^{-k/d}\big)$

for a method of order k. The rate degrades with dimension.

Monte Carlo's O(n^{-1/2}) does not depend on d at all. The variance \sigma^2 may grow, but the rate is untouched.

d	Simpson	Monte Carlo
1	n^{-4}	n^{-1/2}
4	n^{-1}	n^{-1/2}
8	n^{-1/2}	n^{-1/2}
20	n^{-1/5}	n^{-1/2}

The crossover is around d = 8, and beyond it Monte Carlo is the only practical option. That's why it dominates Bayesian inference, statistical physics, and financial derivative pricing — all high-dimensional integration problems.

There's also a robustness advantage: grid methods need smoothness, while Monte Carlo needs only finite variance.

Variance reduction

The rate is fixed at n^{-1/2}, so the only lever is \sigma. Halving \sigma is worth quadrupling n.

Antithetic variates. Use pairs (U, 1-U). If h is monotone the pair is negatively correlated, and

$\operatorname{Var}\left(\frac{h(U) + h(1-U)}{2}\right) = \frac{\sigma^2 + \operatorname{Cov}}{2} < \frac{\sigma^2}{2}$

Free, and often a large win.

Control variates. If g has a known integral and correlates with h:

$\hat I_{cv} = \overline{h} - c\left(\overline{g} - E[g]\right)$

The optimal c is \operatorname{Cov}(h,g)/\operatorname{Var}(g), and the variance is reduced by a factor 1 - \rho^2 — the same expression as §2's conditional variance, and for the same reason.

Stratified sampling. Split the domain and sample each piece — guarantees even coverage instead of leaving it to chance.

Importance sampling. Sample where h matters most. The last lesson of this course.

Worked example

Estimate \pi two ways.

Hit-or-miss. Throw darts uniformly at the unit square and count how many land inside the quarter circle x^2 + y^2 \le 1. That fraction estimates the area \pi/4, so

$\hat\pi = 4 \times \frac{\text{hits}}{n}$

As an integral. Note that

$\int_0^1 \frac{4}{1+x^2}\,dx = 4\arctan(1) = \pi$

so averaging h(x) = \frac{4}{1+x^2} over uniform draws estimates \pi directly.

Both are unbiased, but their variances differ enormously.

Hit-or-miss uses a Bernoulli indicator with p = \pi/4 \approx 0.785, so after the factor of 4 its variance is 16p(1-p) \approx 2.70.

The integral form uses h ranging only over [2, 4], giving

$\operatorname{Var}(h) = E[h^2] - \pi^2 = 16\int_0^1\frac{dx}{(1+x^2)^2} - \pi^2 \approx 10.28 - 9.87 = 0.41$

About 6.5 times smaller.

Same n, same rate, but the second estimator's error is roughly $\sqrt{6.5} \approx 2.6$ times smaller. How you set up the estimator matters more than how many samples you take — a 6.5× variance reduction is worth 6.5× the sample size, for no extra computation.

Doing it in Python

Both estimators of \pi, with error bars:

import numpy as np

rng = np.random.default_rng(0)
n = 1_000_000

# 1. Hit-or-miss
x, y = rng.uniform(0, 1, n), rng.uniform(0, 1, n)
hits = (x**2 + y**2 <= 1)
est1 = 4 * hits.mean()
se1 = 4 * hits.std(ddof=1) / np.sqrt(n)

# 2. As an integral of 4/(1+x^2)
u = rng.uniform(0, 1, n)
h = 4 / (1 + u**2)
est2 = h.mean()
se2 = h.std(ddof=1) / np.sqrt(n)

print(f"{'method':>16} {'estimate':>12} {'std error':>12} {'95% interval':>26}")
for name, e, s in [("hit-or-miss", est1, se1), ("integral form", est2, se2)]:
    print(f"{name:>16} {e:>12.6f} {s:>12.6f} "
          f"{f'({e-1.96*s:.5f}, {e+1.96*s:.5f})':>26}")

print(f"\ntrue pi = {np.pi:.6f}")
print(f"variance ratio: {(se1/se2)**2:.1f}x -- same n, very different precision")

The \sqrt n rate, and why more samples help slowly:

import numpy as np

rng = np.random.default_rng(1)

print(f"{'n':>10} {'estimate':>12} {'error':>12} {'sigma/sqrt(n)':>15}")
for n in (100, 1_000, 10_000, 100_000, 1_000_000):
    u = rng.uniform(0, 1, n)
    h = 4 / (1 + u**2)
    est = h.mean()
    print(f"{n:>10} {est:>12.6f} {abs(est - np.pi):>12.6f} "
          f"{h.std(ddof=1)/np.sqrt(n):>15.6f}")

print("\nError tracks sigma/sqrt(n): 100x more samples for 10x more accuracy.")

Where Monte Carlo overtakes a grid — the dimension crossover:

import numpy as np

rng = np.random.default_rng(2)

# Integrate f(x) = exp(-sum(x^2)) over the unit cube in d dimensions.
def f(X):
    return np.exp(-np.sum(X**2, axis=-1))

# Exact value: (integral of exp(-x^2) from 0 to 1)^d
from scipy.integrate import quad
one_d = quad(lambda t: np.exp(-t**2), 0, 1)[0]

n = 10_000
print(f"{'d':>4} {'exact':>12} {'grid':>12} {'grid err':>11} "
      f"{'Monte Carlo':>13} {'MC err':>10}")
for d in (1, 2, 3, 5, 8):
    exact = one_d ** d

    # Grid: n^(1/d) points per axis (rounded down)
    per_axis = max(int(n ** (1 / d)), 2)
    axis = (np.arange(per_axis) + 0.5) / per_axis
    mesh = np.stack(np.meshgrid(*([axis] * d), indexing="ij"), axis=-1)
    grid_est = f(mesh).mean()

    mc_est = f(rng.uniform(0, 1, size=(n, d))).mean()

    print(f"{d:>4} {exact:>12.6f} {grid_est:>12.6f} {abs(grid_est-exact):>11.6f} "
          f"{mc_est:>13.6f} {abs(mc_est-exact):>10.6f}")

print("\nThe grid has only", int(n ** (1/8)), "points per axis at d = 8 --")
print("its accuracy collapses while Monte Carlo's is unchanged.")

Variance reduction, measured:

import numpy as np

rng = np.random.default_rng(3)
n, trials = 5_000, 2_000

def plain(rng):
    u = rng.uniform(0, 1, n)
    return (4 / (1 + u**2)).mean()

def antithetic(rng):
    u = rng.uniform(0, 1, n // 2)
    h = 4 / (1 + u**2)
    h_anti = 4 / (1 + (1 - u)**2)
    return ((h + h_anti) / 2).mean()

def control_variate(rng):
    u = rng.uniform(0, 1, n)
    h = 4 / (1 + u**2)
    g = u                                   # E[g] = 0.5, and g correlates with h
    c = np.cov(h, g)[0, 1] / np.var(g)
    return (h - c * (g - 0.5)).mean()

results = {}
for name, method in [("plain", plain), ("antithetic", antithetic),
                     ("control variate", control_variate)]:
    results[name] = np.array([method(rng) for _ in range(trials)])

# Plain Monte Carlo is the baseline; the ratio says how many times more
# samples plain would need to match this method's variance.
base_var = results["plain"].var()
for name, ests in results.items():
    print(f"{name:>17}: mean {ests.mean():.6f}  sd {ests.std():.6f}  "
          f"variance ratio vs plain {base_var / ests.var():6.1f}x")

print("\nAntithetic and control variates both cut the spread substantially,")
print("at the same computational cost -- same n, visibly smaller error bars.")

A high-dimensional integral with no realistic alternative:

import numpy as np

rng = np.random.default_rng(4)

# E[max(X_1, ..., X_20)] for independent standard normals -- a 20-dimensional
# integral with no closed form. Grid methods are hopeless here.
d, n = 20, 500_000
X = rng.standard_normal((n, d))
maxima = X.max(axis=1)

est = maxima.mean()
se = maxima.std(ddof=1) / np.sqrt(n)
print(f"E[max of {d} standard normals] = {est:.5f} +/- {1.96*se:.5f}")
print(f"\n(A grid with just 3 points per axis would need 3^{d} = "
      f"{3**d:,} evaluations.)")

Your turn

1. How many samples to halve the Monte Carlo error?

2. Why doesn't the convergence rate depend on dimension?

3. You estimate \int_0^1 e^x dx with n = 10{,}000 and get 1.72 with \hat\sigma = 0.5. Give a 95% interval.

Solutions

1. Four times as many.

The error scales as \sigma/\sqrt n, so

$\frac{\sigma}{\sqrt{4n}} = \frac{\sigma}{2\sqrt n}$

To reduce the error tenfold you'd need 100× the samples; for three more decimal places, a million times.

This is the same \sqrt n tax as §3's standard error, and it's why variance reduction matters more than brute force. Cutting \sigma in half is equivalent to quadrupling n — but usually far cheaper.

2. Because the estimator is an average of iid values, and the CLT gives \sigma/\sqrt n for any average regardless of where the samples came from.

h(X_1), \dots, h(X_n) are iid real numbers whether X_i lives in \mathbb R or \mathbb R^{1000}. The dimension of the input never enters the variance calculation for the output.

Contrast with a grid: covering d dimensions with m points per axis costs m^d evaluations, so n total points gives only n^{1/d} per axis. The spacing — and hence the error — degrades exponentially. Monte Carlo doesn't try to cover the space systematically, so it pays no such penalty.

The caveat: while the rate is dimension-free, \sigma^2 itself may grow with d. Monte Carlo escapes the exponential blow-up, not all difficulty.

3. Standard error:

$\operatorname{SE} = \frac{\hat\sigma}{\sqrt n} = \frac{0.5}{\sqrt{10000}} = \frac{0.5}{100} = 0.005$

95% interval:

$1.72 \pm 1.96(0.005) = 1.72 \pm 0.0098 = (1.710,\; 1.730)$

The true value is e - 1 \approx 1.71828, which falls inside. ✓

Note how routine this is — the interval came from the same sample that produced the estimate, using nothing but the CLT. Deterministic quadrature gives you a number with no error bar unless you separately analyse the method's truncation error, which usually requires bounding derivatives you don't have.

Check yourself in code

Estimate \pi both ways and confirm the integral form has much lower variance.

Print exactly this:

hit-or-miss 3.141368
integral form 3.141768
hit-or-miss SE 0.001642
integral SE 0.000643
integral form is better: True

Use default_rng(0) with 1000000 samples for each method. Round the estimates to 6 significant decimal places as shown and the standard errors to 6 decimal places. Use ddof=1 for the standard deviations.

import numpy as np

rng = np.random.default_rng(0)
n = 1_000_000

x, y = rng.uniform(0, 1, n), rng.uniform(0, 1, n)
hits = (x**2 + y**2 <= 1)
est1 = 4 * hits.mean()
print("hit-or-miss", round(est1, 6))

# Estimate pi again as the mean of 4/(1+u^2) over uniform draws, then print
# both standard errors and whether the integral form has the smaller one.

import numpy as np

rng = np.random.default_rng(0)
n = 1_000_000

x, y = rng.uniform(0, 1, n), rng.uniform(0, 1, n)
hits = (x**2 + y**2 <= 1)
est1 = 4 * hits.mean()
print("hit-or-miss", round(est1, 6))

u = rng.uniform(0, 1, n)
h = 4 / (1 + u**2)
est2 = h.mean()
print("integral form", round(est2, 6))

se1 = 4 * hits.std(ddof=1) / np.sqrt(n)
se2 = h.std(ddof=1) / np.sqrt(n)
print("hit-or-miss SE", round(se1, 6))
print("integral SE", round(se2, 6))
print("integral form is better:", bool(se2 < se1))

Monte Carlo integration rewrites an integral as an expectation and estimates it by averaging. Convergence is O(n^{-1/2}) — slow in one dimension, but independent of dimension, which is why it takes over above about d = 8. The CLT supplies a free error bar, and since the rate is fixed, the real leverage is in reducing \sigma rather than increasing n.

Next: the variance-reduction method that makes rare events computable.

66. Monte Carlo integration