Estimating moments from Monte Carlo simulations

Monte Carlo moments are estimates whose precision depends on the number of independent simulations.

Monte Carlo estimators

For independent values \(X_1,\ldots,X_M\), estimate the first two raw moments and variance by

\[ \widehat m_1=\frac1M\sum_rX_r, \qquad \widehat m_2=\frac1M\sum_rX_r^2, \qquad s^2=\frac1{M-1}\sum_r(X_r-\widehat m_1)^2. \]

Two different kinds of variation

Quantity	Meaning
Sample variance \(s^2\)	Estimated intrinsic variation between epidemic outcomes.
Monte Carlo SE of the mean \(s/\sqrt M\)	Estimated numerical sampling error in \(\widehat m_1\).

More simulations reduce Monte Carlo error, but they do not remove the biological stochastic variance represented by the model.

Interactive Python laboratory

Interactive PythonMonte Carlo moment estimates

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def sis_value(rng, N=100, I0=5, beta=0.25,
              gamma=0.15, T=20):
    I = I0
    time = 0.0
    while I > 0 and time < T:
        b = beta * (N - I) * I / N
        d = gamma * I
        wait = rng.exponential(1 / (b + d))
        if time + wait > T:
            break
        time += wait
        I += 1 if rng.random() < b / (b + d) else -1
    return I

rng = np.random.default_rng(240)
maximum_size = 10000
sample = np.array([
    sis_value(rng) for _ in range(maximum_size)
])

sizes = [50, 100, 500, 1000, 5000, 10000]
rows = []

for M in sizes:
    values = sample[:M]
    mean = values.mean()
    second = np.mean(values ** 2)
    variance = values.var(ddof=1)
    se_mean = np.sqrt(variance / M)

rows.append({
        "M": M,
        "first moment": mean,
        "second moment": second,
        "variance": variance,
        "MC SE of mean": se_mean,
        "95% lower": mean - 1.96 * se_mean,
        "95% upper": mean + 1.96 * se_mean
    })

results = pd.DataFrame(rows)
print(results.round(4).to_string(index=False))

reference_mean = sample.mean()
absolute_error = np.abs(
    results["first moment"] - reference_mean
)

plt.figure(figsize=(9, 5))
plt.plot(results["M"], results["MC SE of mean"],
         "o-", label="Estimated MC SE")
plt.plot(results["M"], absolute_error,
         "s--", label="Absolute difference from M=10000 estimate")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Number of independent simulations")
plt.ylabel("Magnitude (log scale)")
plt.title("Monte Carlo precision improves with sample size")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Interpret the interval correctly

The displayed 95% interval describes Monte Carlo uncertainty in the estimated mean under a normal approximation. It is not the range containing 95% of epidemic outcomes and does not include parameter uncertainty.

Reproducibility checklist

State the seed, sample size, model and observation time.
Use independent simulations.
Report ddof=1 for sample variance.
Increase \(M\) and check that conclusions stabilise.
Do not treat the largest run as the exact truth.

What this lesson adds

You can now estimate raw moments and variance from Monte Carlo samples, calculate the standard error of a mean estimate, and distinguish biological stochasticity from finite-simulation error.