Monte Carlo integration
Here is an idea that sounds like a practical joke: to measure the area of a shape, or the value of an integral, throw darts at it at random and count how many land inside. That is the whole method. It ignores the smoothness of the function, does not care how many dimensions you are integrating over, and returns an answer whose error you can read off the back of an envelope as σ/√N. The catch is the same back-of-envelope: halving the error costs four times the work, and no amount of cleverness removes that square-root wall. So the real subject of this part is not how to estimate an integral with randomness — that is four lines of arithmetic — but how to shrink the constant in front of 1/√N. Antithetic pairs, a control variate with a known mean, and stratification each buy extra accuracy for the same number of function evaluations, often a lot of it, and they are the difference between a simulation that finishes and one that does not.
The question
An integral you cannot do by hand, and a handful of darts
Plenty of integrals have no closed form. The normalising constant of a posterior, the expected loss of a policy over a high-dimensional action space, the average brightness of a scene seen from every possible camera pose — each is an integral or an expectation, and each is one you cannot evaluate with a clever substitution. Classical quadrature fixes a grid of points and weights, and it works beautifully in one dimension. In ten dimensions a grid with only ten points per axis already needs ten billion evaluations, and the cost explodes from there. Randomness sidesteps the grid entirely.
The move is to rewrite the integral as an expectation. For a function h on the unit interval, the average value of h over the interval is its integral, and the average value of h at a uniformly random point is the same number. So draw a point U from the uniform distribution on [0,1] and look at h(U). Its expectation is the integral. One random point is a terrible estimate; a million independent points average out. This is the Monte Carlo estimator, and the remarkable thing is that its accuracy depends on the variance of h, not on the dimension of the space you are sampling.
The price is visible in the error. Because the samples are independent, the variance of the average is the variance of one sample divided by N, so the typical error falls like 1/√N. That is slow. To gain one more correct digit you need a hundred times as many samples. Deterministic quadrature does far better in low dimensions — it can be exponentially accurate on smooth functions — but it cannot survive the curse of dimensionality, and Monte Carlo can. That trade — slower convergence, but indifferent to dimension — is why simulation is the default tool in statistics, physics, finance and machine learning.
The demo below is the oldest picture in the subject. Take the unit square, inscribe a quarter circle of radius one, and throw points uniformly at the square. The fraction that land inside the quarter circle estimates its area, π/4, so four times that fraction estimates π. As the points accumulate, the estimate wanders but its error band narrows like 1/√N. Watch both things at once: the estimate really does converge, and it really does keep moving.
The estimator and its error
Unbiased, and noisy by exactly the amount the variance says
Write the target as an expectation. If X is uniform on [0,1] and h is whatever you want to integrate, then $I=\int_0^1 h(x)\,dx=\mathbb{E}[h(X)]$. Draw N independent uniforms $x_1,\dots,x_N$, evaluate h at each, and average. The result $\hat I_N$ is unbiased: its expectation is the true integral, because expectation is linear and each sample has the right mean. Unbiasedness is the guarantee that the method points at the right answer.
The noise is the other half of the story. Independent samples make the variance of the average exactly σ^2/N, where $\sigma^2=\operatorname{Var}(h(X))$ is the variance of a single evaluation. Taking square roots gives the standard error σ/√N, which is the typical size of the gap between $\hat I_N$ and I. In practice you never know σ, so you estimate it from the same samples — the sample standard deviation of the h(x_i) divided by √N — and that estimated error is what turns the point estimate into an error bar.
That 1/√N is the whole economics of simulation. Four times the work halves the error; a hundred times the work removes one digit. It is called the square-root wall, and it cannot be demolished by better hardware, because more samples is exactly what hardware buys and the wall is made of samples. What can be moved is the constant σ: if two estimators are both unbiased for I, the one with smaller σ is the one you should run. The rest of this part is about making σ small without changing the answer.
Pi by darts is this estimator in disguise, and it is a good first test because the truth is known to ten digits. Sample a point (x_i,y_i) uniformly in the unit square and let h be the indicator that the point lies in the quarter disc $x^2+y^2\le 1$. The expectation of that indicator is the area of the quarter disc, π/4, so 4 times the sample mean estimates π. The per-sample variance is p(1-p) with p=π/4, giving a standard error near 1.64/√N — a constant times 1/√N, as promised.
Drag the number of darts and watch the running estimate on the right. The shaded band is centred on the true π and narrows like 1/√N; the curve should spend most of its time inside. Notice how slowly the band closes — at ten thousand darts it is still about a hundredth of a unit wide — and notice that a different seed redraws a completely different path inside the same band. The randomness is real, but its size is predictable, which is what makes a Monte Carlo answer reportable.
Blue darts land inside the quarter circle, grey darts outside. Each throw is a sample; the fraction inside estimates the area π/4.
The heavy line is $4\cdot(\text{inside}/N)$ as the darts accumulate; the dashed line is π; the shaded band is the $\pm 1/\sqrt{N}$-scaled error band.
One habit is worth forming here. Every random draw in this volume comes from a seeded stream, so the same seed always replays the same run. That is not a limitation; it is what makes a simulation a reproducible experiment rather than an anecdote. When you compare two estimators, you compare them on the same underlying randomness, which removes the luck of the draw from the comparison. The next section does exactly that, holding the target integral fixed and letting four estimators race.
Racing the variance reducers
Same answer, less noise, no extra samples
Take the target to be concrete: $I=\int_0^1 e^x\,dx=e-1\approx 1.71828$. The naive Monte Carlo estimate averages e^{x_i} over uniform points and carries a standard error of about 0.49/√N. Three classical tricks beat that constant while leaving the expectation untouched. None of them is exotic; all of them exploit a piece of structure in the integrand that the naive sampler refuses to look at.
Antithetic variates use the fact that a uniform sample can be reflected. For every point x, also evaluate at 1-x, and average the pair. Each point is still marginally uniform, so the estimate stays unbiased, but the two evaluations are negatively correlated whenever h is monotone — if e^x is unusually small at x, it is unusually large at 1-x. Averaging a negatively correlated pair cancels much of the noise, and since the estimator averages N/2 such pairs, its variance drops to $\tfrac12\sigma^2(1+\rho)$ per pair, with $\rho$ the correlation between h(X) and h(1-X). For a monotone integrand $\rho<0$, and the reflection is essentially free.
Control variates use a function whose mean you already know. Pick g with $\mathbb{E}[g(X)]=\mu_g$ known exactly — for the unit interval the natural choice is g(x)=x with $\mu_g=\tfrac12$. Estimate the integrand minus a multiple of how far g strays from its mean. The correction has zero mean, so the estimator is still unbiased, but subtracting the correlated part of g removes its contribution to the variance. The best coefficient is the regression slope $\beta^*=\operatorname{Cov}(h,g)/\operatorname{Var}(g)$, and the residual variance is the familiar $\sigma_h^2(1-\rho_{hg}^2)$: the fraction of the variance that g cannot explain is what remains. Here x and e^x are strongly related, so almost all of the noise is removed.
Stratification attacks the variance by dividing the problem instead of correlating samples. Split [0,1] into L slices, draw a proportionate number of points in each, and average the per-slice averages with weights w_l. Because the slice means are computed separately, none of the variation between slices leaks into the error; only the variation within a slice survives. Since e^x changes little across a narrow slice, the within-slice variances are tiny and the total variance, $\sum_l w_l^2\sigma_l^2/n_l$, is a small fraction of the naive $\sigma^2/N$. More strata means a stronger effect, and in the limit of infinitely many strata with one sample each the error can fall faster than 1/√N because the integrand's smoothness is finally being used.
The demo races all four estimators on the same integral. Each curve is the root-mean-square error of the estimator over many independent replications, plotted against the number of samples on log axes, so a straight line of slope -1/2 means exactly 1/√N convergence. Every method shows that slope: variance reduction moves the line down, it does not bend it. What changes is the vertical offset, and the offset is precisely the reduction in σ. The naive estimator sits on top; antithetic pairs pull it down sharply because the reflection cancels so much; the control variate and stratification each do better still, because they use more of the structure of e^x. The reference dashed line is the naive slope, drawn to make the constant visible.
RMS error against sample count for four unbiased estimators of $\int_0^1 e^x\,dx$, on log–log axes. Lower is better; the dashed guide has slope -1/2, the signature of $1/\sqrt{N}$.
Two warnings keep this honest. First, variance reduction must preserve unbiasedness, and the control variate only does so because $\mu_g$ is known exactly; estimating it from the same samples introduces a bias that is usually tiny but always there. Second, a badly chosen control variate is worse than none: if h and g are uncorrelated, the correction adds variance instead of removing it, and the optimal $\beta$ is simply zero. The same care applies to stratification, whose gains evaporate if the slices are chosen so that the integrand varies wildly inside them. Variance reduction is not free magic; it is a way of feeding known structure back into the sampler, and it rewards knowing your integrand.
Where this shows up
Every expectation you cannot compute directly
Policy gradients and baselines
Reinforcement learning from human feedback estimates an expected reward by sampling trajectories, and the estimator is a Monte Carlo average with all the noise that implies. The standard fix in RLHF — subtracting a baseline from the reward — is a control variate whose known mean is zero for a centred advantage, and it is the reason policy-gradient training converges at all.
Integrals past three dimensions
The Riemann sums of multiple integrals become hopeless once the dimension climbs, while the Monte Carlo error stays σ/√N regardless of dimension. That dimension-independence is why simulation replaces quadrature in statistical physics, rendering, and Bayesian computation.
Stochastic objectives in optimization
Robust fitting and pose estimation optimise expectations over noise and outliers, and stochastic methods in nonlinear optimization replace the exact objective with a sampled one. Variance reduction is what keeps those sampled gradients from drowning the descent direction in noise.
The samplers it depends on
Every estimator here assumes you can draw the right points, and in high dimensions that is its own problem. The inverse-CDF, rejection and importance-sampling machines of How to draw a sample are what turn a target distribution into the stream of x_i this part averages; importance sampling is itself a variance-reduction technique, aimed at a different representation of the same integral.
Cheat sheet
Every formula in one place
| Idea | Formula | Reading |
|---|---|---|
| Integral as expectation | $\int_0^1 h(x)\,dx=\mathbb{E}[h(U)]$, $U\sim\mathrm{Unif}(0,1)$ | The average value of the integrand is the integral. |
| Basic estimator | $\hat I_N=\frac1N\sum_i h(x_i)$ | Average the evaluations of independent samples; unbiased. |
| Variance and error | $\operatorname{Var}=\sigma^2/N$, error $\approx\sigma/\sqrt N$ | The square-root wall: 4× the samples halves the error. |
| Estimated error | $s/\sqrt N$, $s^2=\frac{1}{N-1}\sum_i(h(x_i)-\hat I_N)^2$ | Replace the unknown $\sigma$ with the sample spread. |
| Pi by darts | $\hat\pi_N=\frac4N\sum_i\mathbf 1\{x_i^2+y_i^2\le1\}$ | Area of the quarter disc is $\pi/4$; se $\approx1.64/\sqrt N$. |
| Antithetic variates | $\tfrac12[h(x)+h(1-x)]$, var $\tfrac12\sigma^2(1+\rho)$ | Reflected pairs are negatively correlated for monotone h. |
| Control variate | $h-\beta(g-\mu_g)$, $\beta^*=\operatorname{Cov}(h,g)/\operatorname{Var}(g)$ | Subtract the correlated part of a function with known mean $\mu_g$. |
| Residual variance | $\sigma_h^2(1-\rho_{hg}^2)/N$ | Only the part of h that g cannot explain remains. |
| Stratified estimator | $\sum_l w_l\,\bar h_l$, var $\sum_l w_l^2\sigma_l^2/n_l$ | Between-slice variation is removed; only within-slice noise survives. |
| What does not change | the expectation I | Variance reduction moves the constant, never the slope or the answer. |
Further reading
Where to go deeper
- J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods, 1964 — the classic treatment of the basic estimator, error bars and variance reduction, written before the field had a name for itself.
- Christian P. Robert and George Casella, Monte Carlo Statistical Methods, 2004 — the estimator in the context of inference, with careful chapters on antithetic and control variates.
- Art B. Owen, Monte Carlo Theory, Methods and Examples, 2013 — freely available and unusually readable, especially on stratification and the quasi-Monte Carlo limit.
- Paul Glasserman, Monte Carlo Methods in Financial Engineering, 2003 — variance reduction as a practical discipline, where every basis point of error is money.
- Sheldon M. Ross, Simulation, 2012 — an elementary account of antithetic, control and stratified sampling with worked examples you can reproduce by hand.