Monte Carlo Integration Simulation

Area as an average

If U is uniform on [a, b], its density is 1/(b − a), so the expected value of f(U) is the integral divided by the interval length. Turned around:

Replace the expectation with a sample mean over n independent uniform points and you have the plain (or crude) Monte Carlo estimator:

In the upper panel the dashed line is that sample mean of the heights. The rectangle under it, width b − a, is the estimate; the shaded region under the curve is the truth. With few points the two areas visibly disagree; press Play and watch the dashed line settle.

An error bar from the same samples

The estimator is an average of independent terms, so its variance is σ²/n, where σ² is the variance of (b − a)·f(U). The lab estimates σ² from the sample itself and reports

where s is the sample standard deviation of the f(x_i). By the central limit theorem the interval covers the true integral in about 95% of runs — so roughly one seed in twenty shows “no” in the coverage readout, and that is the method working as advertised, not failing. Switch the integrand to e^−x² with seed 42 to see one of the misses.

The constant σ depends on the integrand. For sin x on [0, π] it is √(π²/2 − 4) ≈ 0.967, so 1,000 samples give a standard error near 0.03. Flatter integrands have smaller σ and converge faster; the error still shrinks only like 1/√n, so each extra correct digit costs 100 times more samples. How accurate a Monte Carlo estimate is covers the reasoning in general.

Why not just use the midpoint rule?

In one dimension you should. The midpoint rule evaluates f at the centers of n equal subintervals; for a smooth integrand its error falls like 1/n², against 1/√n for Monte Carlo. The last readout shows the midpoint error with the same number of function evaluations, and it is smaller by orders of magnitude. The quarter circle √(1 − x²) is instructive: its slope is infinite at x = 1, and the midpoint error there falls only like n^−1.5 — still far ahead of random sampling.

The picture changes with dimension. A grid with n points in d dimensions has only n^1/d points per axis, so the product midpoint rule's error behaves like n^−2/d. The Monte Carlo error is σ/√n in every dimension, because the argument above never used the dimension: an average of independent samples has variance σ²/n whether each sample is a number or a point in a 100-dimensional cube. For d > 4 the rate of Monte Carlo is better, and for the high-dimensional integrals of physics, statistics and finance — such as pricing an option by simulation — it is often the only practical choice. Keep in mind that σ can also grow with dimension; the rate is dimension-free, the constant is not.

Algorithm

1. Draw x₁, …, x_n independently and uniformly from [a, b].
2. Evaluate y_i = f(x_i) and keep a running mean and sum of squared deviations (Welford's update, which avoids cancellation).
3. Return Î = (b − a)·ȳ and SE = (b − a)·s/√n.

Example result

With sin x on [0, π], 1,000 samples and seed 42, the estimate is 1.95578 against the exact value 2.00000, an error of 0.04422 with a standard error of 0.03048. The 95% interval runs from 1.8960 to 2.0155 (covers the exact value: yes). The midpoint rule with the same 1,000 evaluations is off by only 8.2e-7.

Implementation

The page runs simulation.js from the package below; the Python file is a line-by-line port with the same random-number generator, so a given seed produces the same result in both languages. The excerpts show the core of each.

// Plain Monte Carlo: I = ∫ f(x) dx over [a, b] = (b − a) · E[f(U)], U ~ Uniform(a, b).
// Draw x uniformly, keep a running mean of f(x) and its sum of squared deviations (Welford).
export function sample(st, k, rng, onSample) {
  const { a, b, f } = st.g;
  const end = Math.min(st.n, st.i + k);
  while (st.i < end) {
    const x = a + (b - a) * rng.uniform();
    const y = f(x);
    st.i++;
    const d = y - st.mean;
    st.mean += d / st.i;
    st.m2 += d * (y - st.mean);
    if (onSample) onSample(x, y);
  }
}

// Estimate, standard error from the sample variance, and a normal-theory 95% interval.
export function summarize(st) {
  const { a, b, exact } = st.g;
  const w = b - a;
  const estimate = w * st.mean;
  const se = st.i > 1 ? w * Math.sqrt(st.m2 / (st.i - 1) / st.i) : null;
  return {
    points: st.i,
    estimate,
    exact,
    error: Math.abs(estimate - exact),
    std_error: se,
    ci_low: se === null ? null : estimate - 1.96 * se,
    ci_high: se === null ? null : estimate + 1.96 * se,
    covers: se === null ? null : Math.abs(estimate - exact) <= 1.96 * se,
  };
}

def sample(st, k, rng):
    """Draw up to k points x ~ U(a, b); keep a running mean of f(x) and its sum of squares (Welford)."""
    a, b, f = st["g"]["a"], st["g"]["b"], st["g"]["f"]
    end = min(st["n"], st["i"] + k)
    while st["i"] < end:
        x = a + (b - a) * rng.uniform()
        y = f(x)
        st["i"] += 1
        d = y - st["mean"]
        st["mean"] += d / st["i"]
        st["m2"] += d * (y - st["mean"])


def summarize(st):
    """I = (b - a) * mean f(U); SE from the sample variance; normal-theory 95% interval."""
    g = st["g"]
    w = g["b"] - g["a"]
    estimate = w * st["mean"]
    se = w * math.sqrt(st["m2"] / (st["i"] - 1) / st["i"]) if st["i"] > 1 else None
    return {
        "points": st["i"],
        "estimate": estimate,
        "exact": g["exact"],
        "error": abs(estimate - g["exact"]),
        "std_error": se,
        "ci_low": None if se is None else estimate - 1.96 * se,
        "ci_high": None if se is None else estimate + 1.96 * se,
        "covers": None if se is None else abs(estimate - g["exact"]) <= 1.96 * se,
    }

Run the full model from a terminal:

node simulation.js --integrand sin --n 1000 --seed 42
python simulation.py --integrand sin --n 1000 --seed 42

Download the model

Assumptions & limitations

• Uniform sampling only. This is the plain estimator. Importance sampling, stratification and control variates reduce σ, often dramatically, and are the first things to try when plain Monte Carlo is too slow.
• Finite variance. The 1/√n rate and the interval need f(U) to have finite variance. For integrands with strong singularities the estimate can converge slowly or erratically.
• Normal approximation. The 95% interval uses the central limit theorem and the sample standard deviation; for small n or skewed integrands its real coverage can differ from 95%.
• Pseudo-random points. Samples come from a seeded generator (xoshiro128**); see pseudo-random numbers and seeds. Quasi-random (low-discrepancy) points are a different method with different error behavior.
• One-dimensional demo. The integrands are one-dimensional so the result can be drawn and checked; the advantage of Monte Carlo only appears in higher dimensions.

Related simulations

• OS-001Monte Carlo πEstimate π from random points in a unit square and watch the error shrink like 1/√n.
• OS-003Buffon's needleDrop needles on a lined floor; the share that cross a line is 2ℓ/(πt), which gives an estimate of π.
• OS-005Metropolis–HastingsA Markov chain that samples any density known up to a constant; tune the step size and watch it mix.
• OS-011Central limit theoremAverage n draws from almost any distribution and the means turn normal — at a speed you can measure.

Background: How Accurate Is a Monte Carlo Estimate? · Pseudo-Random Numbers and Seeds in Simulation

References

1. Metropolis, N. & Ulam, S. (1949). The Monte Carlo Method. Journal of the American Statistical Association, 44(247), 335–341. doi.org/10.1080/01621459.1949.10483310
2. Owen, A. B. (2013). Monte Carlo theory, methods and examples. Stanford University. artowen.su.domains/mc
3. Robert, C. P. & Casella, G. (2004). Monte Carlo Statistical Methods, 2nd ed. Springer.