Monte Carlo integration

Deep dive

How it works

In: a sampler for X, an integrand f, sample budget N (or target error) · Out: estimate of E[f(X)] with a standard error · State: count, running mean and M₂ (Welford); one RNG stream per worker · Cost: N evaluations of f; error σ/√N

1. Generator: a statistically strong, seedable PRNG (PCG/xoshiro); independent streams per thread.
2. Sampler: transform uniforms into the target distribution (inverse CDF, Box-Muller or ziggurat for normals, rejection sampling, MCMC when direct sampling is impossible).
3. Estimator: accumulate f(Xᵢ) with Welford's online mean/variance: δ = x − mean; mean += δ/n; M₂ += δ·(x − mean); variance = M₂/(n − 1).
4. Error bar: report mean ± z·√(variance/n); stop when the half-width reaches the target (sequential stopping slightly biases coverage, so use a margin).

Variance reduction:

• Antithetic variates: pair U with 1 − U (or Z with −Z); for monotone f the pair is negatively correlated, so the pair average varies less.
• Control variates: find g correlated with f and with known E[g] (e.g. the geometric-average Asian option as control for the arithmetic one). Estimate β from the same samples (small bias, usually negligible) or from a pilot run.
• Importance sampling: draw from q concentrated where p·|f| is large, average f·p/q. Zero variance in theory when q ∝ p·|f|.
• Stratified sampling: split the domain into strata and sample each proportionally; never increases variance.

Code

/// SplitMix64-based uniform generator (see prng-pcg-xoshiro for stronger options).
pub struct Rng(u64);
impl Rng {
    pub fn new(seed: u64) -> Self { Rng(seed) }
    pub fn next_u64(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9e3779b97f4a7c15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
        z ^ (z >> 31)
    }
    /// Uniform in (0, 1): never exactly 0, safe for ln().
    pub fn uniform(&mut self) -> f64 { ((self.next_u64() >> 11) as f64 + 0.5) * (1.0 / (1u64 << 53) as f64) }
    /// Standard normal pair by Box-Muller.
    pub fn normal_pair(&mut self) -> (f64, f64) {
        let (u1, u2) = (self.uniform(), self.uniform());
        let r = (-2.0 * u1.ln()).sqrt();
        let t = 2.0 * std::f64::consts::PI * u2;
        (r * t.cos(), r * t.sin())
    }
}

/// Welford's streaming mean and variance.
#[derive(Default, Clone, Copy)]
pub struct Stats { pub n: u64, mean: f64, m2: f64 }
impl Stats {
    pub fn push(&mut self, x: f64) {
        self.n += 1;
        let d = x - self.mean;
        self.mean += d / self.n as f64;
        self.m2 += d * (x - self.mean);
    }
    /// Combine two partial results (parallel reduction, Chan et al.).
    pub fn merge(a: Stats, b: Stats) -> Stats {
        if a.n == 0 { return b; }
        if b.n == 0 { return a; }
        let n = a.n + b.n;
        let d = b.mean - a.mean;
        Stats { n, mean: a.mean + d * b.n as f64 / n as f64, m2: a.m2 + b.m2 + d * d * (a.n as f64 * b.n as f64) / n as f64 }
    }
    pub fn mean(&self) -> f64 { self.mean }
    pub fn variance(&self) -> f64 { if self.n > 1 { self.m2 / (self.n - 1) as f64 } else { f64::NAN } }
    /// Standard error of the mean; mean +- 1.96 * stderr is an approximate 95 % interval.
    pub fn stderr(&self) -> f64 { (self.variance() / self.n as f64).sqrt() }
}

/// Plain Monte Carlo of E[f(U)], U ~ Uniform(0, 1).
pub fn plain(f: impl Fn(f64) -> f64, n: u64, rng: &mut Rng) -> Stats {
    let mut s = Stats::default();
    for _ in 0..n { s.push(f(rng.uniform())); }
    s
}

/// Antithetic pairs: each sample is (f(U) + f(1 - U)) / 2; n pairs = 2n evaluations.
pub fn antithetic(f: impl Fn(f64) -> f64, n: u64, rng: &mut Rng) -> Stats {
    let mut s = Stats::default();
    for _ in 0..n { let u = rng.uniform(); s.push(0.5 * (f(u) + f(1.0 - u))); }
    s
}

/// Control variate g with known mean `g_mean`; beta estimated from the same samples.
pub fn control_variate(f: impl Fn(f64) -> f64, g: impl Fn(f64) -> f64, g_mean: f64, n: u64, rng: &mut Rng) -> (f64, f64) {
    let (mut sf, mut sg, mut sfg, mut sgg, mut sff) = (0.0, 0.0, 0.0, 0.0, 0.0);
    for _ in 0..n {
        let u = rng.uniform();
        let (a, b) = (f(u), g(u));
        sf += a; sg += b; sfg += a * b; sgg += b * b; sff += a * a;
    }
    let nf = n as f64;
    let (mf, mg) = (sf / nf, sg / nf);
    let cov = sfg / nf - mf * mg;   // fine here: values O(1); use Welford-style updates for large means
    let var_g = sgg / nf - mg * mg;
    let var_f = sff / nf - mf * mf;
    let beta = cov / var_g;
    let est = mf - beta * (mg - g_mean);
    let var_cv = var_f - cov * cov / var_g; // = var_f (1 - rho^2)
    (est, (var_cv / nf).sqrt())
}

/// P(Z > a) for standard normal Z by importance sampling from N(a, 1): weight = exp(-a z + a^2 / 2).
pub fn tail_prob_is(a: f64, n: u64, rng: &mut Rng) -> Stats {
    let mut s = Stats::default();
    for _ in 0..n / 2 {
        let (z1, z2) = rng.normal_pair();
        for z in [z1 + a, z2 + a] {
            s.push(if z > a { (-a * z + 0.5 * a * a).exp() } else { 0.0 });
        }
    }
    s
}

Complexity and constants

• Error after N samples ≈ σ/√N; 10⁶ samples of an O(1)-variance integrand give about 3 correct digits.
• Measured with the code above, μ = ∫₀¹ eᵘ du = e − 1, 10⁶ evaluations each: plain standard error 4.9·10⁻⁴ (σ² = 0.242); antithetic 8.9·10⁻⁵ (variance 31× smaller at equal cost); control variate g(u) = u 6.3·10⁻⁵ (61× smaller). All three intervals contain the true value.
• Rare event P(Z > 4) = 3.17·10⁻⁵ with 10⁶ samples: plain MC sees 35 hits, estimate 3.5·10⁻⁵ ± 5.9·10⁻⁶ (17 % relative error); importance sampling from N(4, 1) gives 3.173·10⁻⁵ ± 6.7·10⁻⁸ (0.2 %), a variance reduction near 8,000×.
• Quasi-Monte Carlo error is O((log N)^d / N) for integrands of bounded variation; in practice randomized Sobol points often beat plain MC substantially for smooth integrands, even in tens of dimensions.

Efficiency

• Cost: N evaluations of f plus the sampler; error σ/√N, so 10⁶ samples of an O(1)-variance integrand give about 3 digits and halving the error costs 4× the samples. Measured with the code above (∫₀¹ eᵘ du, 10⁶ evaluations): plain standard error 4.9·10⁻⁴, antithetic 8.9·10⁻⁵, control variate 6.3·10⁻⁵. Memory is O(1) per stream. Bound: cost of f (compute, or memory latency in path tracing).
• Levers:
  • antithetic pairs — variance 31× smaller at equal cost (measured above) — monotone integrands
  • control variate with known mean — variance 61× smaller (measured above) — a correlated cheap g exists
  • importance sampling — about 8,000× variance reduction for P(Z > 4) (measured above) — rare events, q covering the tails of p·|f|
  • quasi-Monte Carlo (randomized Sobol) — error near O(1/N) instead of O(1/√N) (page above) — smooth integrands
  • multilevel Monte Carlo — big savings by mixing coarse and fine simulations (Giles 2008) — SDE path simulation
  • independent stream per thread, merged Stats — near-linear scaling — all parallel runs
• Naive → standard → state of the art:
  • plain MC, fixed N · error σ/√N (4.9·10⁻⁴ at N = 10⁶ above)
  • plain MC with antithetic or control variates, Welford error bars · 6.3·10⁻⁵ at the same N
  • randomized quasi-MC or multilevel MC with importance sampling where it fits · error decays faster than 1/√N

Implementation notes

• Parallel: one generator stream per thread (PCG increments, xoshiro jump()), local Stats, then Stats::merge; results are reproducible if the thread-to-stream mapping is fixed.
• Report the error bar with every estimate; check that it shrinks like 1/√N (if not, the variance may be infinite, e.g. heavy-tailed importance weights).
• Effective sample size for weighted samples: (Σw)²/Σw²; a small ESS means importance sampling is failing.
• Common random numbers: compare two systems (A/B, parameter sensitivity) with the same random stream; the difference has far lower variance than independent runs.
• Normals: Box-Muller is simple; ziggurat is faster; inverse CDF (Acklam/Wichura) is needed for quasi-Monte Carlo and stratification.
• Rendering: multiple importance sampling (Veach's balance heuristic) combines light and BSDF sampling robustly.
• MCMC (Metropolis-Hastings, HMC) when direct sampling is impossible; samples are correlated, so the error uses the effective sample size from autocorrelation, not N.

Pitfalls

• Naive variance formula E[x²] − E[x]² on data with a large mean: negative or wildly wrong variances (floating-point error); measured on 10⁹ + (i mod 10), it returns 2048 instead of 8.25, Welford returns 8.25.
• Importance distribution with lighter tails than the target: finite-looking estimates with infinite variance, error bars that lie.
• Reusing the same seed for "independent" replicas, or seeding parallel streams with consecutive integers in a weak generator.
• Using rand() % n or low-quality generators (old LCGs): correlations show up in high-dimensional sampling.
• Stopping as soon as the interval looks good (optional stopping): undercoverage; set N in advance or use sequential methods with corrected thresholds.

Variants and alternatives

In production

Finance: derivative pricing and risk engines (QuantLib Monte Carlo framework). Physics: Geant4 and MCNP/OpenMC particle transport. Graphics: path tracers (pbrt, Cycles, Arnold, RenderMan). Statistics: Stan and PyMC (MCMC), SciPy stats.qmc (Sobol, Halton). Earth and health: climate and epidemiology ensembles.

Implementations

Sources

• Paul Glasserman, Monte Carlo Methods in Financial Engineering, Springer, 2003. book
• Art B. Owen, Monte Carlo Theory, Methods and Examples, online book, 2013. book
• B. P. Welford, Note on a Method for Calculating Corrected Sums of Squares and Products, Technometrics 4(3), 1962. paper
• Tony F. Chan, Gene H. Golub, Randall J. LeVeque, Algorithms for Computing the Sample Variance: Analysis and Recommendations, The American Statistician 37(3), 1983. paper
• Eric Veach, Robust Monte Carlo Methods for Light Transport Simulation, PhD thesis, Stanford University, 1997. thesis
• Michael B. Giles, Multilevel Monte Carlo Path Simulation, Operations Research 56(3), 2008. paper