Monte Carlo integration approximates integrals by random sampling. The core idea is to rewrite an integral as an expectation and then estimate that expectation by an average of random samples.
A geometric intuition is measuring the area of an irregular pond inside a known field area \(A\). If \(N\) random throws are made and \(N_{\text{pond}}\) hits land in the pond, then
\[ A_{\text{pond}}\approx \frac{N_{\text{pond}}}{N}A. \]
From Integrals to Expectations
Let \(f:[a,b]\to[0,\infty)\) and let \(p\) be a density on \([a,b]\) with \(\int_a^b p(x)\,dx=1\). Define
\[ I(f):=\int_a^b f(x)\,dx. \]
Then, whenever \(p(x)>0\) on the support of \(f\):
\[ \begin{aligned} I(f) &=\int_a^b f(x)\,dx\\ &=\int_a^b \frac{f(x)}{p(x)}p(x)\,dx\\ &=\mathbb{E}[g(X)],\quad g(x):=\frac{f(x)}{p(x)},\; X\sim p. \end{aligned} \]
So the Monte Carlo estimator with samples \(X_i\sim p\) is
\[ I_N(f):=\frac{1}{N}\sum_{i=1}^N g(X_i) =\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)}. \]
By the law of large numbers, \(I_N(f)\to I(f)\) as \(N\to\infty\).
Uniform Sampling on \([a,b]\)
If \(X_i\sim\mathrm{unif}(a,b)\), then \(p(x)=\frac{1}{b-a}\), so
\[ I_N(f) =\frac{1}{N}\sum_{i=1}^N\frac{f(X_i)}{p(X_i)} =\frac{b-a}{N}\sum_{i=1}^N f(X_i). \]
This is the standard one-dimensional Monte Carlo integration formula.
Error and Convergence Rate
The estimator is unbiased:
\[ \mathbb{E}[I_N(f)] = I(f). \]
Its variance is
\[ \mathrm{Var}(I_N(f)) = \frac{\mathrm{Var}(g(X))}{N}, \]
so the standard error scales as
\[ \sqrt{\mathrm{Var}(I_N(f))} = \frac{\sigma_g}{\sqrt{N}}. \]
Hence reducing error by a factor of \(10\) usually requires roughly \(100\) times more samples.
Variance Reduction
Importance Sampling
Choose \(p\) so that samples are concentrated where \(f\) contributes most. Intuitively, this reduces \(\mathrm{Var}(g)\). A common principle is to choose \(p\) approximately proportional to \(|f|\) on the integration interval.
Stratified Sampling
Split \([a,b]\) into sub-intervals (strata), sample each stratum separately, and combine weighted estimates. This enforces coverage of the whole domain and often lowers variance.
Rejection Sampling
When direct sampling from the desired \(p\) is difficult, use proposal samples from an easier density and accept/reject points with a correction rule. Accepted points follow the target distribution.
Examples
Constant Function
For \(f(x)=k\), the exact integral is \(k(b-a)\), and Monte Carlo recovers the same value in expectation:
\[ \int_a^b f(x)\,dx =\int_a^b k\,dx \approx \frac{b-a}{N}\sum_{i=1}^N k =k(b-a). \]
Approximating \(\pi\)
Consider the unit circle \(x^2+y^2=1\) restricted to the first quadrant, where \(x\in[0,1]\) and \(y\in[0,1]\). Its area formula \(A=\pi r^2\) with \(r=1\) gives the full circle area \(\pi\), so the quarter-circle confined to the unit square \([0,1]^2\) has area
\[ A_{\text{circle}}=\frac{\pi r^2}{4}=\frac{\pi}{4}. \]
The bounding square has area \(A_{\text{square}}=1\), so the ratio of the two areas is exactly \(\pi/4\), and solving for \(\pi\) gives
\[ \pi = 4\cdot\frac{A_{\text{circle}}}{A_{\text{square}}}. \]
This ratio is what Monte Carlo integration estimates. With uniform random points \((x_i,y_i)\in[0,1]^2\), define the indicator
\[ h(x_i,y_i)= \begin{cases} 1, & x_i^2+y_i^2\le 1,\\ 0, & \text{otherwise}, \end{cases} \]
so that \(A_{\text{circle}}\) is the two-dimensional integral of \(h\) over the unit square, estimated the same way as \(I_N(f)\) above:
\[ A_{\text{circle}}=\int_0^1\int_0^1 h(x,y)\,dx\,dy\approx\frac{1}{N}\sum_{i=1}^N h(x_i,y_i). \]
Multiplying by \(4\) turns this into an estimator for \(\pi\) itself:
\[ \pi \approx \frac{4}{N}\sum_{i=1}^N h(x_i,y_i). \]
Implementation
A JavaScript implementation for uniform Monte Carlo integration:
function integrate(f, a, b, N) {
let sum = 0;
for (let i = 0; i < N; i++) {
sum += f(a + Math.random() * (b - a));
}
return (b - a) * sum / N;
}