Monte Carlo Methods — Random Sampling for Simulation & Integration

1. The Core Idea: Sample → Estimate

Every Monte Carlo method follows the same three-step pattern:

1. Define a random variable X whose expected value E[X] equals the quantity you want.
2. Draw N independent samples x₁, x₂, …, x_N from the distribution of X.
3. Estimate the quantity as the sample mean: μ̂ = (x₁ + x₂ + … + xN) / N.

By the Law of Large Numbers, μ̂ converges to E[X] as N → ∞. By the Central Limit Theorem, the error is approximately Gaussian with standard deviation σ/√N, where σ² = Var(X).

2. Estimating π with Random Darts

The classic introduction: throw N random darts uniformly at the unit square [0,1]². Count the fraction that land inside the quarter-circle of radius 1 centred at the origin. Since the quarter-circle has area π/4, the fraction converges to π/4.

function estimatePi(N) {
  let inside = 0;
  for (let i = 0; i < N; i++) {
    const x = Math.random();
    const y = Math.random();
    if (x * x + y * y <= 1) inside++;
  }
  return 4 * inside / N;
}

console.log(estimatePi(1_000_000));  // ≈ 3.1415…  (error ~ 0.001)

With N = 10⁶ the typical error is about 0.002 — roughly 3 significant figures. To gain one more significant figure you need 100× more samples: N = 10⁸. This 1/√N convergence is characteristic of all Monte Carlo methods.

3. Monte Carlo Integration

The general goal is to compute the integral:

If we can sample points uniformly from Ω with volume |Ω|, the estimator is:

// Integrate f(x) = x² over [0, 1] — exact answer: 1/3
function mcIntegrate(f, a, b, N) {
  let sum = 0;
  for (let i = 0; i < N; i++) {
    const x = a + Math.random() * (b - a);
    sum += f(x);
  }
  return (b - a) * sum / N;
}

const result = mcIntegrate(x => x * x, 0, 1, 1e6);
console.log(result);   // ≈ 0.3333…

// Multi-dimensional: volume of 4D unit hypersphere (exact: π²/2 ≈ 4.935)
function hypersphereMC(d, N) {
  let inside = 0;
  for (let i = 0; i < N; i++) {
    let r2 = 0;
    for (let j = 0; j < d; j++) { const x = 2*Math.random()-1; r2 += x*x; }
    if (r2 <= 1) inside++;
  }
  return Math.pow(2, d) * inside / N;  // volume of enclosing hypercube = 2^d
}
console.log(hypersphereMC(4, 1e6));   // ≈ 4.93

4. Variance and Convergence Rate

The key insight: halving the error costs 4× more samples (O(1/ε²) in error ε). This is worse than deterministic methods for smooth 1D integrals (which converge as 1/N^k for some k ≥ 2), but better than deterministic methods in high dimensions (where the grid size explodes).

Reducing Variance without More Samples

Several techniques reduce σ(f) without increasing N:

• Control variates: subtract a correlated function g with known integral, then add it back analytically.
• Antithetic variates: pair each sample x with its mirror 1−x; their function values are negatively correlated.
• Stratified sampling: divide the domain into sub-regions and sample each proportionally.
• Quasi-Monte Carlo: replace pseudo-random sequences with low-discrepancy sequences (Halton, Sobol) for O(log(N)^d/N) convergence.

5. Importance Sampling

If f(x) is large only in a small region, most uniform samples contribute almost nothing to the estimator. Importance sampling concentrates samples where f is large by sampling from an alternative distribution p(x) and correcting with the likelihood ratio:

The ratio w(x) = f(x)/p(x) is called the importance weight. If p(x) closely mimics the shape of |f(x)|, the weights are nearly constant and variance plummets. The catch: we need a proposal p that we can sample from and that covers all regions where f is non-negligible (otherwise weights can explode — the notorious heavy-tail problem).

// Estimate ∫₀^∞ x·exp(−x) dx = 1 using Exponential(1) as proposal
// p(x) = exp(−x),  f(x)/p(x) = x·exp(-x)/exp(-x) = x
function sampleExponential() {
  return -Math.log(1 - Math.random());   // inverse-CDF method
}

function importanceSamplingEstimate(N) {
  let sum = 0;
  for (let i = 0; i < N; i++) {
    const x = sampleExponential();   // x ~ Exp(1)
    sum += x;                         // weight f(x)/p(x) = x
  }
  return sum / N;   // converges to 1 with zero variance for this integrand!
}
console.log(importanceSamplingEstimate(100));   // ≈ 1.000 (exact!)

6. Markov Chain Monte Carlo (MCMC)

When the target distribution p(x) is complex (e.g. high-dimensional posterior in Bayesian inference) we cannot sample from it directly. MCMC constructs a Markov chain whose stationary distribution is p(x), then uses the chain's trajectory as correlated samples.

The Metropolis–Hastings Algorithm

1. Start at any state x₀.
2. Propose a new state x′ from a symmetric proposal q(x′|xᵢ) (e.g. Gaussian step).
3. Compute acceptance ratio: α = min(1, p(x′)/p(xᵢ)).
4. Accept x′ with probability α; otherwise stay at xᵢ.
5. Repeat, building up a trajectory (x₀, x₁, x₂, …).

// Metropolis–Hastings: sample from p(x) ∝ exp(−x²/2) (standard Normal)
function metropolis(logP, start, stepSize, N) {
  const samples = [start];
  let current = start;
  let logPcur  = logP(current);

  for (let i = 1; i < N; i++) {
    const proposal = current + (Math.random() * 2 - 1) * stepSize;
    const logPprop  = logP(proposal);
    const logAlpha  = logPprop - logPcur;       // log acceptance ratio

    if (Math.log(Math.random()) < logAlpha) {   // accept
      current = proposal;
      logPcur = logPprop;
    }
    samples.push(current);
  }
  return samples;
}

// Target: standard normal   log p(x) = −x²/2  (constant dropped)
const draws = metropolis(x => -x * x / 2, 0, 1, 50000);
// After warm-up, draws follow N(0,1) — compute mean and variance to verify
const mean = draws.reduce((a, b) => a + b) / draws.length;   // ≈ 0

MCMC requires a burn-in period (early samples are autocorrelated and biased by the starting point) and has effective sample size smaller than N due to autocorrelation. Practical diagnostics include the Gelman–Rubin R̂ statistic and trace-plot inspection.

7. JavaScript: Practical Tips

Pseudo-Random Number Generator Quality

Math.random() in V8 uses the xorshift128+ algorithm — fast, 128-bit state, period 2¹²⁸ − 1. It passes all standard statistical tests and is fine for simulations. For cryptographic applications use crypto.getRandomValues() instead.

Seeding for Reproducibility

// Simple LCG PRNG (not crypto-safe — for reproducible simulations only)
function makePRNG(seed) {
  let s = seed & 0xffffffff;
  return function() {
    s = (Math.imul(1664525, s) + 1013904223) & 0xffffffff;
    return (s >>> 0) / 4294967296;
  };
}
const rng = makePRNG(42);
console.log(rng(), rng());   // same sequence every time

Vectorised Batches with TypedArrays

// Fill a Float64Array with N random points — faster than looping Math.random()
function mcPiBatch(N) {
  const data = new Float64Array(N * 2);
  for (let i = 0; i < data.length; i++) data[i] = Math.random();
  let inside = 0;
  for (let i = 0; i < N; i++) {
    const x = data[2*i], y = data[2*i+1];
    if (x*x + y*y <= 1) inside++;
  }
  return 4 * inside / N;
}

8. Applications across Disciplines

Particle Physics & High-Energy Simulations

The GEANT4 toolkit simulates particle interactions in detectors using Monte Carlo sampling of cross-sections, decay modes, and scattering angles. A single LHC collision can generate millions of secondary particles, all tracked probabilistically through the detector material.

Financial Option Pricing (Black–Scholes MC)

The price of a European option is the discounted expected payoff under the risk-neutral measure. Since the stock path follows geometric Brownian motion, we simulate N independent paths, compute the payoff for each, and average:

Bayesian Inference

MCMC (particularly Hamiltonian Monte Carlo / NUTS, as in Stan and PyMC) allows scientists to infer posterior distributions of parameters given data, even with hundreds of correlated parameters. This powers everything from climate model calibration to genomic association studies.

Computer Graphics: Path Tracing

Rendering a physically based image means evaluating the rendering integral — the incoming radiance integrated over all directions and all light paths. Path tracing estimates this integral by randomly sampling light bounces. Modern renderers (Cycles, Arnold, V-Ray) are all Monte Carlo path tracers.