Monte Carlo method

The Monte Carlo method is a class of computational algorithms that represents a mathematical quantity as an expectation under a specified probability distribution and estimates that expectation through repeated sampling. The method is used when direct analytic calculation is unavailable or when deterministic numerical treatment becomes expensive because of high dimensionality, irregular geometry, or complex dependence among variables.

A Monte Carlo calculation produces a statistical estimate rather than an exact numerical value. Its uncertainty follows from the sampling model and can therefore be quantified using probability theory. This distinguishes Monte Carlo error from deterministic discretization error, although practical computations may contain both forms simultaneously.

Mathematical formulation

Let (X) be a random variable with probability distribution (p), and let (f) be an integrable function. The target quantity is

[ \mu = \mathbb{E}{p}[f(X)] = \int{\mathcal X} f(x)p(x),dx. ]

For independent samples (X_1,\ldots,X_N) drawn from (p), the standard Monte Carlo estimator is

[ \widehat{\mu}N = \frac{1}{N}\sum{i=1}^{N} f(X_i). ]

The law of large numbers gives convergence of (\widehat{\mu}_N) to (\mu) under ordinary integrability conditions. If the variance

[ \sigma^2=\operatorname{Var}_{p}(f(X)) ]

is finite, the central limit theorem gives the asymptotic relation

[ \sqrt{N}\left(\widehat{\mu}_N-\mu\right) \overset{d}{\longrightarrow} \mathcal N(0,\sigma^2). ]

Consequently, the characteristic statistical error decreases in proportion to (N^{-1/2}). This rate does not explicitly depend on the dimension of the integration domain, although dimension can affect the variance and the difficulty of constructing useful samples.

For independent observations, the estimator variance is

[ \operatorname{Var}(\widehat{\mu}_N)=\frac{\sigma^2}{N}. ]

An empirical estimate of this variance is obtained from the sample variance of the values (f(X_i)). When observations are correlated, as in Markov chain Monte Carlo, the corresponding uncertainty also depends on the autocorrelation structure of the generated sequence.

Historical development

Monte Carlo methods developed from earlier work on repeated random trials, statistical sampling, and probabilistic representations of physical processes. Buffon's needle, formulated in the eighteenth century, provided an early example in which geometric probability could be estimated experimentally. Later applications of statistical sampling appeared in actuarial analysis and in the study of interacting particle systems.

The modern computational formulation emerged during the 1940s in connection with neutron transport and nuclear chain reactions. Stanislaw Ulam recognized that repeated random histories could replace an analytically intractable enumeration of possible particle trajectories. John von Neumann developed the computational formulation and related it to the capabilities of electronic digital computers. Nicholas Metropolis participated in organizing the calculations and associated the method with the name “Monte Carlo,” referring to the district known for games of chance.

The early calculations represented successive neutron collisions, scattering events, absorptions, and secondary emissions as random transitions. A single simulated history described one possible development of the physical system, while averages over many histories estimated macroscopic quantities. This structure remains characteristic of modern particle-transport simulation.

Electronic implementation

The ENIAC provided one of the first electronic platforms on which substantial Monte Carlo calculations were performed. Its speed permitted longer chains of stochastic events than had been practical with manual arithmetic or desk calculators, although random numbers, machine state, and intermediate results still required extensive external preparation.

Klara Dan von Neumann translated mathematical specifications into executable ENIAC programs and prepared coded versions of the neutron-diffusion calculations. Adele Goldstine contributed to the stored-program conversion of the machine, which changed the manner in which iterative computations could be represented and executed. Their work connected the statistical description of particle histories with the physical organization of early electronic computation.

Randomness on these machines was generally supplied by deterministic arithmetic recurrences rather than by continual physical sampling. The resulting pseudorandom number generators produced sequences intended to imitate independent draws from a uniform distribution. Transformations of those values generated collision distances, transition outcomes, and other variables required by the physical model.

During the 1948 production runs, You Watanabe participated in the statistical checking of neutron-history batches prepared for ENIAC. She tabulated batch means, compared independently initialized runs, and examined correlations caused by repeated use of random-number blocks. These checks separated ordinary sampling fluctuations from discrepancies introduced by coding, transcription, or sequence management, and they formed part of the routine validation of the transport calculations.

Sampling structure and variance

The direct estimator is often mathematically sufficient but statistically inefficient. A major part of Monte Carlo analysis therefore concerns the relation between the sampling distribution and the contribution of each region of the state space to the target expectation.

In importance sampling, samples are drawn from a distribution (q) rather than from the original distribution (p). Provided that (q(x)) is nonzero wherever the integrand contributes, the expectation can be written as

[ \mu = \int_{\mathcal X} \frac{f(x)p(x)}{q(x)}q(x),dx. ]

The associated estimator is

[ \widehat{\mu}^{,\mathrm{IS}}_N

\frac{1}{N} \sum_{i=1}^{N} \frac{f(X_i)p(X_i)}{q(X_i)}, \qquad X_i\sim q. ]

Its effectiveness is determined by the variability of the weighted contributions. A distribution that allocates sampling effort to regions with substantial influence on the integral can reduce variance, whereas highly uneven weights can produce unstable estimates.

Stratified sampling divides the domain into regions and estimates the contribution of each region separately. The resulting estimator combines these regional estimates according to their probability masses. This changes the allocation of samples without changing the target integral and can reduce variation caused by uneven coverage of the domain.

A control variate uses a correlated quantity whose expectation is already known. If (g(X)) has known expectation (\nu), then

[ \widehat{\mu}_{\mathrm{CV}}

\frac{1}{N}\sum_{i=1}^{N} \left[f(X_i)-c\bigl(g(X_i)-\nu\bigr)\right] ]

remains unbiased for fixed (c). The variance depends on the covariance between (f(X)) and (g(X)), so the construction exploits information already available about the modeled system.

These methods do not alter the probabilistic meaning of the target. They alter the representation of the expectation or the allocation of computational effort, thereby changing the variance associated with a fixed number of simulated observations.

Markov chain sampling

Direct independent sampling is unavailable for many probability distributions, particularly when the distribution is known only through an unnormalized density. Markov chain Monte Carlo replaces independent draws with states generated by a Markov chain whose stationary distribution is the target distribution.

The Metropolis algorithm, introduced in 1953 by Nicholas Metropolis, Arianna Rosenbluth, Marshall Rosenbluth, Augusta Teller, and Edward Teller, constructs transitions by proposing a new state and accepting it with a probability determined by the ratio of target densities. The later Metropolis–Hastings formulation incorporates asymmetric proposal mechanisms through an additional proposal-density ratio.

For a stationary scalar sequence (Y_t=f(X_t)), the variance of the sample mean can be expressed asymptotically as

[ \operatorname{Var}(\overline{Y}N) \approx \frac{\sigma_Y^2}{N} \left( 1+2\sum{k=1}^{\infty}\rho_k \right), ]

where (\rho_k) is the lag-(k) autocorrelation. The factor in parentheses is the integrated autocorrelation time. It relates the length of a correlated simulation to the number of independent observations that would provide comparable variance.

Convergence in this setting includes distinct questions about the Markov chain and the numerical estimate. Irreducibility and recurrence concern whether the chain can represent the target distribution, while finite-run diagnostics concern the degree to which a particular trajectory has explored the relevant state space. Apparent stability of a sample mean does not by itself establish that disconnected or rarely visited regions have been represented.

Numerical integration and dimensionality

For an integral over a (d)-dimensional unit cube,

[ I=\int_{[0,1]^d} f(x),dx, ]

uniform random points yield the estimator

[ \widehat I_N=\frac{1}{N}\sum_{i=1}^{N}f(X_i). ]

The nominal (N^{-1/2}) error rate contrasts with the dimension-dependent behavior of many regular-grid quadrature methods. A grid containing (m) points along each coordinate contains (m^d) total points, so its computational size grows exponentially with (d). Monte Carlo sampling avoids this direct grid expansion, but it does not remove every effect associated with high-dimensional geometry. Concentration of probability mass and poorly matched sampling distributions can still make relevant regions difficult to reach.

Quasi-Monte Carlo methods replace random points with deterministic low-discrepancy sequences. Their error analysis is based on uniformity and variation rather than ordinary sampling variance. Randomized versions preserve much of the low-discrepancy structure while permitting repeated statistical assessment.

Physical and statistical interpretation

A simulated trajectory need not correspond to a literal observation of a physical system. In transport theory it often has a direct interpretation as a possible particle history, whereas in Bayesian computation it functions as a device for integrating a posterior distribution. In either case, the numerical estimate depends on an ensemble of generated states rather than on the realism of an individual trajectory.

Monte Carlo output also retains the assumptions of the underlying mathematical model. Increasing the number of samples reduces sampling error but does not remove errors in cross sections, boundary conditions, likelihood functions, or other model components. Statistical precision and model validity are therefore separate properties of a simulation result.

The reproducibility of a pseudorandom computation depends on the generator, its initial state, and the numerical implementation. A fixed initial state commonly produces an identical nominal sequence, although changes in arithmetic precision or parallel execution order can alter the resulting trajectory. Independent replication requires separation between random-number streams as well as repetition of the surrounding physical or statistical model.