Expected infectious population and its variance

Moment solutions become biologically useful when expected burden and variability are interpreted separately and together.

Two complementary summaries

\[\mu_I(t)=\mathbb E[I(t)],\qquad \sigma_I^2(t)=\operatorname{Var}(I(t)).\]

The mean describes the ensemble centre at each time. Variance describes path-to-path spread around that mean. Neither is one realised epidemic.

Standard deviation and coefficient of variation

\[\sigma_I(t)=\sqrt{\sigma_I^2(t)},\qquad \mathrm{CV}(t)=\frac{\sigma_I(t)}{\mu_I(t)}.\]

Standard deviation has units of people. CV is dimensionless and measures variation relative to the mean, but becomes unstable when the mean is close to zero.

Interactive Python laboratory

Interactive PythonInterpreting mean and variance

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

N, I0 = 100, 5
beta, gamma = 0.25, 0.15
times = np.linspace(0, 60, 601)

def rhs(t, y):
    m1, m2 = y
    m3 = 3 * m1 * m2 - 2 * m1 ** 3
    dm1 = (beta - gamma) * m1 - beta * m2 / N
    dm2 = (
        (2 * (beta - gamma) - beta / N) * m2
        - 2 * beta * m3 / N
        + (beta + gamma) * m1
    )
    return [dm1, dm2]

sol = solve_ivp(
    rhs, (0, 60), [I0, I0 ** 2],
    t_eval=times, rtol=1e-8, atol=1e-10
)
if not sol.success:
    raise RuntimeError(sol.message)

mean_I, second_I = sol.y
variance_I = np.maximum(second_I - mean_I ** 2, 0)
sd_I = np.sqrt(variance_I)
cv_I = np.divide(
    sd_I, mean_I,
    out=np.full_like(sd_I, np.nan),
    where=mean_I > 1e-10
)

peak_mean_index = np.argmax(mean_I)
peak_variance_index = np.argmax(variance_I)

print(f"Maximum mean in window: {mean_I[peak_mean_index]:.3f}")
print(f"Time of maximum mean in window: {times[peak_mean_index]:.3f}")
print(f"Peak variance: {variance_I[peak_variance_index]:.3f}")
print(f"Time of peak variance: {times[peak_variance_index]:.3f}")
print()

selected = [0, 10, 20, 40, 60]
idx = [int(t * 10) for t in selected]
table = pd.DataFrame({
    "time": selected,
    "expected I": mean_I[idx],
    "variance I": variance_I[idx],
    "SD I": sd_I[idx],
    "CV I": cv_I[idx]
})
print(table.round(4).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
plt.plot(times, mean_I, color="#174c34",
         linewidth=2, label="Expected I")
plt.fill_between(
    times,
    np.maximum(mean_I - sd_I, 0),
    mean_I + sd_I,
    color="#b9ddc7", alpha=0.65,
    label="Mean ± one SD"
)
plt.xlabel("Time")
plt.ylabel("Infectious population")
plt.title("Expected prevalence and stochastic spread")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

plt.figure(figsize=(9.5, 4.9))
plt.plot(times, variance_I, color="#b42318",
         label="Variance")
plt.axvline(
    times[peak_variance_index],
    color="black", linestyle="--",
    label="Variance peak"
)
plt.xlabel("Time")
plt.ylabel("Variance of I")
plt.title("Uncertainty can peak at a different time from the mean")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Mean ± SD is not automatically a probability interval

For a normal distribution, approximately 68% lies within one standard deviation of the mean. Epidemic distributions may be skewed, bounded or contain extinction mass, so this interpretation must not be assumed from moments alone.

Biological meaning

A high expected infectious population suggests high average burden. A high variance means scenarios with the same parameters can differ substantially. The timing of maximum uncertainty can be different from the timing of maximum expected burden.

What this lesson adds

You can now interpret expected prevalence, variance, standard deviation and relative variability, identify maxima within a stated observation window and avoid presenting a mean–SD display as an exact probability band.