Monte Carlo method

The Monte Carlo method is a class of computational algorithms that estimate numerical quantities by repeated random sampling. Rather than solving a problem analytically, a Monte Carlo procedure draws many random inputs from a probability distribution, evaluates a deterministic function on each, and aggregates the results to approximate expectations, integrals, or the distribution of an output. The approach is especially useful for high-dimensional integration, for systems with many coupled degrees of freedom, and for problems whose deterministic solution is intractable, where the estimation error decreases at a rate that is largely independent of the dimension of the problem.

Because the estimate is itself a random variable, a Monte Carlo computation is normally reported together with a statistical measure of its uncertainty, such as a standard error or confidence interval, and much of the theory of the method concerns reducing that uncertainty for a given amount of computation.

Principle

Most Monte Carlo estimators target an expectation of the form $\mu =\mathbb{𝔼}[f(X)]$, where $X$ is a random variable with a known distribution and $f$ is a function of interest. Drawing independent samples $X_1,\dots ,X_n$ and forming the sample mean ${\hat{\mu }}_n=\frac{1}{n}\sum _{i=1}^{n}f(X_i)$ yields an unbiased estimator whose consistency follows from the Law of large numbers and whose fluctuations are described, for large $n$, by the Central limit theorem. A defining feature is that the standard error of ${\hat{\mu }}_n$ shrinks in proportion to $1/\sqrt{n}$ regardless of the dimension of $X$, which contrasts with deterministic quadrature rules whose cost grows exponentially with dimension. This convergence rate is the principal reason the method is applied to multidimensional integrals arising in statistical physics, Bayesian inference, and operations research.

Numerical integration can be recast in this form: an integral $\underset{\mathrm{\Omega }}{\int }g(x)dx$ equals the expectation of $g(X)/p(X)$ when $X$ is sampled from a density $p$ supported on $\mathrm{\Omega }$, so any integral can in principle be estimated by sampling. The quality of the estimate depends heavily on the choice of sampling distribution, motivating the variance-reduction techniques described below.

History

The modern method was developed in the 1940s at the Los Alamos Laboratory, where Stanislaw Ulam, John von Neumann, and Nicholas Metropolis applied statistical sampling on early electronic computers to problems in neutron diffusion that were difficult to treat by deterministic means. Ulam recounted that the idea occurred to him while considering the probability of winning a game of solitaire, and the code name "Monte Carlo", after the casino, was suggested by Metropolis in reference to an uncle of Ulam's who borrowed money to gamble there. The first published exposition aimed at a statistical audience was the 1949 paper by Metropolis and Ulam, which framed the technique as a way of studying integro-differential equations through random models.

A landmark extension came in 1953, when Metropolis and co-authors introduced an algorithm for sampling from the Boltzmann distribution of a many-particle system by simulating a Markov chain that accepts or rejects proposed moves according to an energy criterion. Hastings generalized this construction in 1970 to arbitrary proposal distributions, giving the Metropolis–Hastings algorithm its general statistical form. The method entered mainstream statistics in the following decades, catalyzed by the Gibbs-sampling framework of Geman and Geman and by the demonstration by Gelfand and Smith that Markov chain sampling could compute the marginal posterior densities central to Bayesian inference.

Variance reduction

Because the error of a crude Monte Carlo estimate scales as $\sigma /\sqrt{n}$, reducing the variance ${\sigma }^2$ of the per-sample estimator delivers the same accuracy with fewer samples, and a family of variance-reduction techniques has been developed for this purpose. Importance sampling draws from an alternative distribution chosen to concentrate samples where the integrand is large and reweights them accordingly, which can dramatically lower variance and is central to estimating the probabilities of rare events. Control variates exploit a correlated auxiliary quantity with known expectation to subtract off part of the estimator's fluctuation, while antithetic variates and stratified sampling reduce variance by inducing negative correlation or by partitioning the sample space. Related deterministic constructions known as quasi-Monte Carlo methods replace pseudorandom points with low-discrepancy sequences and can achieve faster convergence for sufficiently smooth integrands.

Markov chain Monte Carlo

When direct sampling from a target distribution is infeasible—as is typical for high-dimensional Bayesian posteriors—Markov chain Monte Carlo (MCMC) constructs a Markov chain whose stationary distribution is the desired target, then treats successive states of the chain as (correlated) samples. The Metropolis–Hastings algorithm builds such a chain by proposing a candidate state and accepting it with a probability that preserves the target as the invariant distribution. Gibbs sampling is a special case in which each coordinate is updated in turn from its full conditional distribution, an approach introduced in an image-restoration context and later shown to be broadly applicable to statistical models. Because MCMC samples are serially correlated, assessing convergence to the stationary distribution and estimating the effective sample size are essential parts of a sound analysis.

Workforce management and contact centers

In workforce management, analytical models such as the Erlang formulas provide closed-form staffing estimates only under restrictive assumptions—typically stationary arrivals, exponential service times, and a single homogeneous skill. When those assumptions fail, discrete-event simulation, a form of Monte Carlo experiment in which random arrival and service times drive an event-scheduled model of the system over simulated time, is the standard alternative for evaluating staffing plans. A contact center can be represented as a queueing system whose randomness—interarrival times, handle times, abandonment, and agent availability—is sampled repeatedly so that performance measures such as service level, occupancy, and average speed of answer are obtained as statistical estimates with confidence intervals rather than single-point analytical values. This makes simulation well suited to features that resist closed-form treatment, including time-varying (non-stationary) call volumes, multi-skill routing, customer patience and abandonment, and shrinkage, at the cost of requiring careful input modeling, replication, and output analysis to control sampling error. Rare-event scenarios relevant to service-level risk, such as the probability of a severe understaffing event during irregular operations, can in turn be estimated more efficiently with importance sampling than with naive replication.

See also

• Stochastic optimization
• Mathematical optimization