Mean paths, variance and uncertainty bands

An ensemble becomes scientifically useful when its centre, spread and simulation error are calculated and labelled correctly.

Time-specific ensemble summaries

For \(M\) simulated paths evaluated at time \(t_n\), the sample mean and variance are

\[\overline X(t_n)=\frac{1}{M}\sum_{r=1}^{M}X_r(t_n),\]

\[s^2(t_n)=\frac{1}{M-1}\sum_{r=1}^{M}[X_r(t_n)-\overline X(t_n)]^2.\]

The mean may be fractional even when the biological variable represents people.

Three different intervals

Display	Interpretation
Mean ± one SD	Scale of path-to-path spread; not generally a fixed probability interval.
5th–95th percentiles	Central 90% of simulated outcomes at each time.
Mean ± 1.96 SE	Monte Carlo uncertainty in the estimated mean, where \(\mathrm{SE}=s/\sqrt M\).

These answer different questions and must not share the vague label “confidence band.”

Interactive Python laboratory

Interactive PythonMean, variance and bands

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

N = 500
I0 = 10
beta = 0.25
gamma = 0.15
T = 80
dt = 0.02
M = 1000
seed = 94

steps = int(round(T / dt))
times = np.linspace(0, T, steps + 1)
rng = np.random.default_rng(seed)

I = np.zeros((M, steps + 1))
I[:, 0] = I0

for n in range(steps):
    current = I[:, n]
    infection_rate = beta * (N - current) * current / N
    recovery_rate = gamma * current
    drift = infection_rate - recovery_rate
    diffusion = np.sqrt(np.maximum(
        infection_rate + recovery_rate, 0
    ))
    dW = np.sqrt(dt) * rng.normal(0, 1, size=M)
    I[:, n + 1] = np.clip(
        current + drift * dt + diffusion * dW,
        0, N
    )

mean_I = I.mean(axis=0)
variance_I = I.var(axis=0, ddof=1)
sd_I = np.sqrt(variance_I)
lower_quantile = np.quantile(I, 0.05, axis=0)
upper_quantile = np.quantile(I, 0.95, axis=0)

standard_error = sd_I / np.sqrt(M)
mean_ci_lower = mean_I - 1.96 * standard_error
mean_ci_upper = mean_I + 1.96 * standard_error

selected_days = [0, 10, 20, 40, 60, 80]
indices = [int(day / dt) for day in selected_days]

summary = pd.DataFrame({
    "day": selected_days,
    "mean I": mean_I[indices],
    "variance I": variance_I[indices],
    "SD I": sd_I[indices],
    "5th percentile": lower_quantile[indices],
    "95th percentile": upper_quantile[indices],
    "MC SE of mean": standard_error[indices]
})

print(summary.round(3).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
plt.fill_between(
    times, lower_quantile, upper_quantile,
    color="#b9ddc7", alpha=0.75,
    label="5th–95th percentiles"
)
plt.plot(times, mean_I, color="#174c34",
         linewidth=2, label="Mean path")
plt.xlabel("Time (days)")
plt.ylabel("Infectious population")
plt.title("Path variability: central 90% simulation band")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

plt.figure(figsize=(9.5, 5.2))
plt.fill_between(
    times, mean_ci_lower, mean_ci_upper,
    color="#f2c7c2", alpha=0.8,
    label="Approximate 95% Monte Carlo interval for mean"
)
plt.plot(times, mean_I, color="#b42318",
         linewidth=1.8, label="Estimated mean")
plt.xlabel("Time (days)")
plt.ylabel("Mean infectious population")
plt.title("Precision of the estimated mean")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Pointwise does not mean simultaneous

The 5th and 95th percentiles are calculated separately at each time. Approximately 90% of path values lie within them at any selected time, but this does not mean 90% of complete trajectories remain inside the band for the entire interval.

Variance is not parameter uncertainty

The variance above measures intrinsic path-to-path stochasticity with fixed parameters. If \(\beta\) and \(\gamma\) are uncertain, parameter uncertainty requires another layer of sampling or statistical inference.

What this lesson adds

You can now calculate time-dependent mean and sample variance, construct percentile bands, quantify Monte Carlo error in the mean, and state clearly which kind of uncertainty each display represents.