Mean, variance, covariance and probability visualisations

The final lesson matches each moment or probability to a display with correct units, scale and biological interpretation.

Each quantity needs its own interpretation

Quantity	Units	Suitable display
Mean	Original population units	Line through time
Variance	Squared population units	Separate line or panel
Covariance	Product of two variables’ units	Signed line, matrix or heatmap
Probability	Dimensionless, between 0 and 1	Probability curve or table

Putting all four on one vertical axis can conceal scale and meaning. Separate panels preserve honest labels.

Scenario

Use 3,000 stochastic SIR epidemics and observe \(S(t)\) and \(I(t)\) daily. Report expected infectious prevalence, its variance, susceptible–infectious covariance and \(\Pr(I(t)\ge30)\).

Interactive Python laboratory

Interactive PythonMoment and probability report

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

N, I0 = 200, 3
beta, gamma = 0.25, 0.10
times = np.arange(0, 81)
M = 3000
threshold = 30
rng = np.random.default_rng(320)

def sir_path():
    S, I, R = N - I0, I0, 0
    time = 0.0
    event_times = [0.0]
    S_values = [S]
    I_values = [I]

while I > 0 and time < times[-1]:
        b = beta * S * I / N
        d = gamma * I
        wait = rng.exponential(1 / (b + d))
        if time + wait > times[-1]:
            break
        time += wait

if rng.random() < b / (b + d):
            S -= 1
            I += 1
        else:
            I -= 1
            R += 1

event_times.append(time)
        S_values.append(S)
        I_values.append(I)

pos = np.searchsorted(
        event_times, times, side="right"
    ) - 1
    return (
        np.array(S_values)[pos],
        np.array(I_values)[pos]
    )

S_paths = np.zeros((M, len(times)))
I_paths = np.zeros((M, len(times)))

for run in range(M):
    S_paths[run], I_paths[run] = sir_path()

mean_I = I_paths.mean(axis=0)
variance_I = I_paths.var(axis=0, ddof=1)
covariance_SI = np.array([
    np.cov(S_paths[:, n], I_paths[:, n], ddof=1)[0, 1]
    for n in range(len(times))
])
probability_threshold = (
    I_paths >= threshold
).mean(axis=0)

selected = [0, 10, 20, 40, 60, 80]
table = pd.DataFrame({
    "day": selected,
    "mean I": mean_I[selected],
    "variance I": variance_I[selected],
    "Cov(S,I)": covariance_SI[selected],
    "P(I>=30)": probability_threshold[selected]
})

print("Moment and probability summary")
print(table.round(3).to_string(index=False))

fig, axes = plt.subplots(2, 2, figsize=(11, 8.2))

axes[0, 0].plot(times, mean_I, color="#174c34")
axes[0, 0].set_title("A. Expected infectious population")
axes[0, 0].set_ylabel("People")
axes[0, 0].grid(alpha=0.2)

axes[0, 1].plot(times, variance_I, color="#b42318")
axes[0, 1].set_title("B. Variance of infectious population")
axes[0, 1].set_ylabel("People²")
axes[0, 1].grid(alpha=0.2)

axes[1, 0].plot(times, covariance_SI, color="#2563eb")
axes[1, 0].axhline(0, color="black", linewidth=0.8)
axes[1, 0].set_title("C. Susceptible–infectious covariance")
axes[1, 0].set_xlabel("Day")
axes[1, 0].set_ylabel("People²")
axes[1, 0].grid(alpha=0.2)

axes[1, 1].plot(
    times, probability_threshold, color="#d97706"
)
axes[1, 1].set_title("D. Probability I is at least 30")
axes[1, 1].set_xlabel("Day")
axes[1, 1].set_ylabel("Probability")
axes[1, 1].set_ylim(0, 1)
axes[1, 1].grid(alpha=0.2)

fig.suptitle("SIR ensemble: moments and a threshold probability")
fig.tight_layout()

Output

Run the code to see the result.

Read the panels together

The mean shows expected prevalence burden.
Variance shows when paths differ most strongly.
Negative covariance shows that epidemics with fewer susceptible people tend to have more infectious people at that time.
The threshold probability directly answers a stated risk question and is not interchangeable with the mean.

Presentation rules

State whether moments are analytical, closed approximations or Monte Carlo estimates.
Give sample size and seed for simulated moments.
Use squared units for variance and covariance.
Keep probability axes between 0 and 1.
Define thresholds and observation times.
Accompany plots with a compact exact-value table.

Part 5 completion

You can now define and estimate moments, derive first- and second-moment equations, identify unclosed nonlinear hierarchies, apply and test closures, solve retained ODEs, validate against simulation, approximate probabilities and communicate mean, variance, covariance and risk accurately.