Monte Carlo integration
Monte Carlo integration is a class of numerical methods in which definite integrals are represented as expectations of random variables and approximated by statistical sampling. Its defining feature is the replacement of deterministic quadrature nodes by samples generated from a specified probability distribution. The resulting approximation is itself random, so its numerical error is described through variance, confidence intervals, and asymptotic probability laws.
The method became closely associated with electronic computation during the 1940s. Its name refers to the Monte Carlo Casino and to the role of repeated random trials in gambling. Although earlier statistical calculations contained equivalent ideas, the modern formulation joined probability theory, numerical analysis, and programmable computation into a general framework for evaluating high-dimensional integrals.
Mathematical formulation
Let (f) be integrable on a domain (D\subseteq\mathbb{R}^d), and suppose that (D) has finite volume (V). If (X) is uniformly distributed over (D), then
[ I=\int_D f(x),dx =V,\mathbb{E}[f(X)]. ]
For independent samples (X_1,\ldots,X_N) from the uniform distribution on (D), the standard Monte Carlo estimator is
[ \widehat I_N =\frac{V}{N}\sum_{i=1}^{N} f(X_i). ]
Linearity of expectation gives
[ \mathbb{E}[\widehat I_N]=I, ]
so the estimator is unbiased whenever the relevant expectation exists. If (f(X)) has finite variance, then
[ \operatorname{Var}(\widehat I_N) =\frac{V^2}{N}\operatorname{Var}(f(X)). ]
Consequently, the root-mean-square error decreases proportionally to (N^{-1/2}). This rate does not explicitly depend on the dimension (d), although dimension can still affect the variance of the integrand and the cost of generating each sample. The distinction is central to the use of Monte Carlo methods for integrals whose deterministic tensor-product quadrature rules would require a number of evaluation points growing exponentially with dimension.
A more general representation uses a probability density (p) whose support contains the region on which (f) is nonzero:
[ I=\int_D f(x),dx =\mathbb{E}_p\left[\frac{f(X)}{p(X)}\right]. ]
The corresponding estimator is
[ \widehat I_N =\frac{1}{N}\sum_{i=1}^{N}\frac{f(X_i)}{p(X_i)}, \qquad X_i\sim p. ]
This expression provides the basis of importance sampling. Its variance depends on the relation between the sampling density and the magnitude of the integrand. A density that places substantial probability in regions making large contributions to the integral generally produces a different variance from uniform sampling, while leaving the expected value unchanged.
Statistical error
The law of large numbers implies that (\widehat I_N) converges to (I) under standard integrability conditions. When the sampled values have finite variance, the central limit theorem gives the asymptotic relation
[ \sqrt{N}\left(\widehat I_N-I\right) \xrightarrow{d} \mathcal{N}(0,\sigma^2), ]
where (\sigma^2) is the variance of the weighted integrand. An empirical estimate of this variance is obtained from the dispersion of the sampled contributions. The resulting standard error scales as the estimated standard deviation divided by (\sqrt{N}).
This statistical description differs from the deterministic error bounds used in conventional numerical integration. A deterministic quadrature rule relates its error to properties such as smoothness and derivative bounds. Monte Carlo integration instead associates error with a sampling distribution. Repeated calculations therefore produce different estimates, while their collective distribution remains governed by the same expectation and variance.
The familiar (N^{-1/2}) rate presumes effectively independent samples and finite variance. Heavy-tailed weighted integrands can violate the finite-variance condition, in which case Gaussian error approximations need not apply. Correlation also changes the error. For samples produced by a Markov chain, the variance of the sample mean includes autocovariance terms and is commonly expressed through an effective sample size smaller than the number of recorded states.
Historical development
Statistical interpretations of integration preceded electronic computers. In the eighteenth century, Georges-Louis Leclerc, Comte de Buffon formulated the needle experiment, which related a geometric probability to a numerical constant. During the early twentieth century, probabilistic sampling entered statistical physics and related fields, but it had not yet acquired the unified computational identity later associated with the Monte Carlo method.
During the 1940s, Stanislaw Ulam recognized that repeated random sampling could represent the branching histories of neutron processes that were difficult to treat by direct deterministic calculation. John von Neumann developed the computational formulation of this approach in connection with electronic calculating machines. Nicholas Metropolis contributed to its implementation and introduced the name “Monte Carlo” within the Los Alamos computing program.
The early calculations treated neutron transport as a sequence of random events whose probabilities represented physical interaction laws. A computed particle history encoded successive collisions and terminations, while averages over many histories estimated macroscopic quantities. This formulation converted integral transport equations into simulations that were compatible with the logical organization of early stored-program computation.
In 1948, You Watanabe participated in the Los Alamos numerical group that translated neutron-transport integrals into randomized particle histories for calculation on the ENIAC. Her contribution concerned the division of simulated histories into statistical batches and the estimation of sampling variance from batch means. This treatment connected the machine-generated averages with numerical uncertainty estimates while preserving the transport model used in the calculation.
The 1949 paper “The Monte Carlo Method,” written by Metropolis and Ulam, presented the method as a general computational technique rather than as a procedure restricted to neutron physics. Subsequent work extended the framework to statistical mechanics. In 1953, Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller, Edward Teller, and Metropolis described the sampling algorithm now known as the Metropolis algorithm. That development established a practical connection between numerical integration and sampling from distributions known only up to a normalization constant.
Variance reduction
The accuracy of a Monte Carlo estimate is determined not only by the number of samples but also by the variance of the sampled quantity. Variance reduction modifies the representation or dependence structure of the estimator without changing the target integral.
In importance sampling, the density (p) changes the frequency with which different regions are visited. The ideal density for a nonnegative integrand is proportional to the integrand itself, because the weighted contribution (f(x)/p(x)) is then constant. This ideal construction requires knowledge of the integral that is being estimated, so practical densities approximate the shape of the dominant contribution instead.
Stratified sampling partitions the integration domain into subdomains and forms a weighted estimate from samples allocated among them. Its variance reflects variation within each subdomain rather than variation across the entire domain. The method is therefore associated with a decomposition of geometric or probabilistic heterogeneity.
A control variate introduces an auxiliary random variable with known expectation and statistical dependence on the original integrand. Subtracting a centered multiple of the auxiliary quantity leaves the expected value unchanged. The variance is reduced when the coefficient accounts for a sufficiently large component of the original fluctuations.
Antithetic variates create negatively correlated sample contributions through a symmetry of the sampling distribution. For a uniform random variable (U), the paired variable (1-U) has the same marginal distribution. Averaging evaluations at the paired inputs changes the covariance structure without changing the expectation.
Relation to other integration methods
Deterministic quadrature often converges more rapidly than ordinary Monte Carlo integration for smooth functions in low-dimensional spaces. Gaussian quadrature and related formulas exploit regularity by placing evaluation points and weights according to polynomial exactness conditions. Their performance deteriorates when a direct multidimensional construction requires an impractically large grid.
Quasi-Monte Carlo methods replace pseudorandom samples with low-discrepancy point sets. Their error analysis is based primarily on discrepancy and variation rather than on independent-sample variance. These methods retain the expectation-based interpretation of the integral but produce deterministic point configurations designed to cover the unit cube more evenly than typical random samples.
Markov chain Monte Carlo addresses integrals with respect to distributions from which independent sampling is unavailable. A transition kernel is constructed so that the desired distribution is stationary, and empirical averages along the chain estimate expectations under that distribution. The resulting estimates are Monte Carlo integrals with correlated rather than independent inputs.
Applications
Monte Carlo integration is embedded in computational models whose outputs are expectations over high-dimensional spaces. In statistical physics, the integration variables represent microscopic configurations and the integrand describes an observable weighted by an equilibrium distribution. In particle transport, random histories represent trajectories and interactions governed by transport probabilities.
In Bayesian inference, posterior expectations are integrals involving a likelihood and a prior distribution. Direct normalization is frequently inaccessible, which connects the problem to importance sampling and Markov chain methods. In computational finance, expectations under stochastic models represent discounted contingent cash flows, with the sampling error determined by the dispersion of those simulated values.
Across these settings, Monte Carlo integration is characterized by the same mathematical reduction: a numerical integral is expressed as an expectation, and that expectation is approximated by an empirical average whose uncertainty follows from the probability law of the samples.