Monte Carlo integration

How do you find the average depth of a lake without surveying every square metre? Drop a measuring rod at many random spots, average the depths you read, and multiply by the lake's surface area — you get the volume, and the more random spots you sample, the better the estimate. Monte Carlo integration is exactly this idea turned into a way to compute integrals: estimate an integral by averaging the integrand at random points.

The core formula is one line. To approximate I = integral over a region D of f(x) dx, draw N points x_1, ..., x_N uniformly at random in D and form the average of f, then scale by the region's volume V: I is approximately (V / N) * sum of f(x_i). Why does this work? The integral equals V times the AVERAGE value of f over D, and the sample mean of f at random points estimates that average — the law of large numbers guarantees it converges to the true value as N grows. (More generally, if you sample points from a density p instead of uniformly, you average f(x_i)/p(x_i); this freedom is what importance sampling exploits.)

The headline property — and the reason this method is indispensable — is how its error behaves. The statistical error shrinks like O(1 / sqrt(N)), and crucially that rate is INDEPENDENT OF THE DIMENSION of the integral. Classical grid-based quadrature in d dimensions needs roughly k^d points for the same accuracy, which explodes (the curse of dimensionality) and becomes hopeless past a handful of dimensions; Monte Carlo just keeps plodding at 1/sqrt(N) whether the integral is over 3 variables or 300. The honest catch is that 1/sqrt(N) is SLOW: to cut the error in half you need four times as many samples, and to gain one decimal digit you need a hundred times as many. So Monte Carlo loses to grids in low dimensions but wins decisively in high ones — which is why it powers high-dimensional finance, physics, rendering, and Bayesian statistics.