Monte Carlo simulation

The formula

The whole method is one estimator. You want some quantity Q that depends on random inputs — an average, a probability, an integral. Draw N samples x⁽¹⁾, x⁽²⁾, …, x⁽ᴺ⁾ from the input distributions, push each through the model f, and average:

Q̂ = (1/N) · Σ f(x⁽ⁱ⁾)

Reading: the sample mean of the model outputs is your estimate of the true expected value. The law of large numbers guarantees it converges; everything else is bookkeeping.

The estimate comes with its own error bar for free:

SE(Q̂) = s / √N        (s = standard deviation of the f(x⁽ⁱ⁾) outputs)

Reading: the standard error shrinks like 1/√N — to halve the error you need four times the runs, and the model's dimensionality never appears. Ten inputs or ten thousand, the rate is the same. That dimension-independence is the entire reason the method exists.

For estimating a probability — the case engineers care about most — count hits:

p̂ = k/N,     relative SE ≈ 1/√(p·N)

Reading: if k of your N runs land in the failure region, k/N estimates the failure probability, and the relative error depends on how many failures you actually observed, not on N alone. Rare events need enormous N (see below).

Where you meet it

Tolerance stack-up at the drawing review. Worst-case says the five-part stack is ±0.050 in and the assembly won't close. Draw each part uniformly across its ±0.010 band, run 200,000 virtual assemblies, and the stack comes out with σ ≈ 0.0129, so ±3σ ≈ ±0.039 — and now you can also read off the exact fraction that violates the gap requirement, something the RSS formula can't tell you when the requirement is one-sided or the model is nonlinear.
RF link budget margin. The spreadsheet version adds nominal gains and losses and quotes one margin number. The Monte Carlo version draws antenna gain, pointing loss, rain fade, and noise figure from their distributions and reports the probability the link closes — which is what the customer actually asked for. A 3 dB nominal margin can hide a 20% outage probability if two loss terms are wide and skewed.
Cost and schedule risk at the program review. Every line item gets a low/likely/high triangle distribution, the simulation rolls the whole WBS ten thousand times, and the board gets an S-curve: "70% confidence of finishing under $4.2M." That chart in the quarterly review is a Monte Carlo output, whether the analyst said so or not.
Measurement uncertainty when the GUM math breaks down. When the measurement equation is nonlinear or an input distribution is asymmetric, the linearized GUM propagation misleads. GUM Supplement 1 (JCGM 101:2008) makes Monte Carlo the official alternative: propagate the full distributions, not just the variances, with runs on the order of 10⁶ for a trustworthy 95% coverage interval ^[7].

How it works

Three steps, no cleverness required: characterize each input as a distribution, sample all inputs jointly and run the deterministic model once per sample, then treat the pile of outputs as data — histogram it, take percentiles, count exceedances. The model can be a closed-form equation, a circuit simulator, or a 6-DOF trajectory code; Monte Carlo doesn't care, which is why it's the tool of last resort that also happens to work first.

The classic demonstration is estimating π: throw random points in a unit square, count the fraction landing inside the quarter circle, multiply by 4. With 10⁶ points the standard error is 4·√(p(1−p)/N) ≈ 0.0016, and a typical run lands within about 0.001 of π. Slow — but the same machinery prices an option, sizes a heat shield margin, or clears a stack-up, unchanged.

What actually bites people:

Rare events eat samples. To estimate a 10⁻⁴ failure probability to 10% relative error you need about 10⁶ runs, because relative SE ≈ 1/√(p·N) and you need roughly 100 observed failures before the estimate means anything. Ten thousand runs with zero failures does not demonstrate p < 10⁻⁴; it demonstrates you didn't run enough. Importance sampling and subset methods exist precisely to push samples into the tail.
Correlated inputs. Sampling each input independently when they aren't — resistors from the same reel, machined features cut in the same fixture setup, rain fade and antenna wetting loss — is the most common silent error. Independence assumptions usually make the output look tighter than reality. If two inputs share a cause, they need a joint distribution or a correlation structure, not two separate draws.
The input distributions are the analysis. The simulation propagates whatever you feed it with complete fidelity, including fiction. A ±1% resistor is not Gaussian centered at nominal — vendors sort, and the distribution can be truncated, bimodal, or biased to one side of the band. Garbage distributions in, confident-looking garbage out, now with a professional histogram.
Convergence is honest but slow. 1/√N never accelerates. For smooth low-dimensional integrals, quadrature or a linearized propagation beats Monte Carlo badly; the method wins when the model is nonlinear, discontinuous, high-dimensional, or a black box. Quasi-random (Sobol) sequences and Latin hypercube sampling buy a real speedup for the same run count and cost nothing to adopt.
Report the seed and the error bar. A Monte Carlo result without N and a standard error is a random number. Fix the RNG seed so the analysis is repeatable, and check stability by doubling N once — if the answer moves more than the claimed error bar, the claimed error bar was wrong.

The mistake people make at review boards: quoting a percentile deep in the tail of the output histogram — a "99.87% reliable" from 5,000 runs — as if it were measured. The middle of a Monte Carlo histogram converges fast; the tails converge last, and they also lean hardest on the distribution assumptions you were least sure about. The method quantifies the uncertainty you modeled. It says nothing about the uncertainty you didn't.

History

The idea has an ancestor: in 1777, Georges-Louis Leclerc, Comte de Buffon, showed that a needle of length l dropped on a floor ruled with lines spaced d apart (needle no longer than the spacing, l ≤ d) crosses a line with probability 2l/(π·d) — meaning you could, in principle, estimate π by throwing needles and counting ^[6]^[8]. Physical random sampling as a computation. It stayed a curiosity for 170 years because nobody could throw needles fast enough.

Enrico Fermi got there next, quietly. In 1930s Rome he used statistical sampling to work neutron-diffusion problems by hand, kept the trick to himself, and never published ^[1]^[2]. In 1947 at Los Alamos, with the electronic computers occupied elsewhere, he had L.D.P. King build the FERMIAC — a hand-cranked brass trolley that rolled across scale drawings of a reactor, tracing individual sampled neutron histories through the geometry ^[2]^[5]. An analog Monte Carlo machine, sitting in a museum case now.

The method proper starts with a sick man and a deck of cards. In 1946 Stanislaw Ulam, recovering from emergency surgery for encephalitis, was playing Canfield solitaire and wondered what fraction of deals could be won. The combinatorics were hopeless; he realized that dealing a hundred hands and counting wins would answer the question — and that the same trick would work on the neutron multiplication problems stalling the weapons program ^[1]^[2]^[3]. John von Neumann saw it immediately, worked out how to make a deterministic computer produce the needed random numbers (his middle-square generator), and in the spring of 1948 the team ran the first automated statistical sampling of a fission core on the ENIAC ^[1]^[2]. Nicholas Metropolis supplied the code name: Ulam had an uncle who borrowed money from relatives because he "just had to go to Monte Carlo" ^[2]. The name stuck, and secrecy required one anyway. Metropolis and Ulam published the open literature version, "The Monte Carlo Method," in the Journal of the American Statistical Association in September 1949 ^[3]^[4] — and the method has scaled with every computer built since, because embarrassingly parallel was its birthright.

Related tools

/tools/link-budget — the canonical candidate for Monte Carlo treatment: replace nominal gains and losses with distributions and get outage probability instead of one margin number
/tools/voltage-divider — resistor-tolerance stack-up on the output ratio is a two-input Monte Carlo you can sanity-check by hand
/tools/resistor-e-series — the tolerance bands you'd sample from; real vendor distributions inside those bands are the gotcha noted above
/tools/bearing-life-l10 — Weibull life models feed fleet-level Monte Carlo reliability simulations
/tools/tsiolkovsky-delta-v — draw Isp and dry mass from their uncertainty bands and the propellant margin question becomes a Monte Carlo run

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →