Path plots, density plots, bands, histograms and tables

The final SDE lesson combines visual and tabular summaries while keeping individual paths, time-specific distributions and numerical uncertainty conceptually separate.

Choose a display from the biological question

Display	Question answered
Path plot	What possible histories and fluctuations look like
Mean with percentile band	How centre and time-specific spread change
Density at selected times	How the shape of the distribution evolves
Outcome histogram	How frequently final or trajectory-level outcomes occur
Table	What the exact reported summaries are

Lines and bands need precise labels

A mean path is not one realisation. A percentile band is pointwise unless a simultaneous construction is used. A few displayed sample paths illustrate variation but do not estimate frequencies.

Density versus histogram

A histogram shows bin counts or frequencies. A density display is scaled so its total area is approximately one, which permits comparison between samples of different sizes. Both depend on bin width, so avoid interpreting small bumps as biological discoveries.

Interactive Python laboratory

Interactive PythonComplete SDE results report

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 = 500
seed = 180

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)
low_I, high_I = np.quantile(
    I, [0.05, 0.95], axis=0
)

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],
    "median I": np.median(I[:, indices], axis=0),
    "SD I": I[:, indices].std(axis=0, ddof=1),
    "5th percentile": low_I[indices],
    "95th percentile": high_I[indices],
    "fraction at zero": (I[:, indices] == 0).mean(axis=0)
})

print("Time-specific ensemble summary")
print(summary.round(3).to_string(index=False))

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

# A. A readable subset of paths
for r in range(15):
    axes[0, 0].plot(times, I[r],
                    alpha=0.45, linewidth=0.8)
axes[0, 0].set_title("A. Fifteen possible paths")
axes[0, 0].set_xlabel("Time (days)")
axes[0, 0].set_ylabel("Infectious population")
axes[0, 0].grid(alpha=0.2)

# B. Ensemble mean and pointwise central band
axes[0, 1].fill_between(
    times, low_I, high_I,
    color="#b9ddc7", alpha=0.75,
    label="5th–95th percentiles"
)
axes[0, 1].plot(times, mean_I,
                color="#174c34", linewidth=2,
                label="Mean")
axes[0, 1].set_title("B. Centre and time-specific spread")
axes[0, 1].set_xlabel("Time (days)")
axes[0, 1].set_ylabel("Infectious population")
axes[0, 1].legend()
axes[0, 1].grid(alpha=0.2)

# C. Density estimates at selected times
density_days = [10, 40, 80]
for day in density_days:
    values = I[:, int(day / dt)]
    density, edges = np.histogram(
        values, bins=28, range=(0, N), density=True
    )
    centres = (edges[:-1] + edges[1:]) / 2
    axes[1, 0].plot(
        centres, density, label=f"Day {day}"
    )
axes[1, 0].set_title("C. Distribution changes through time")
axes[1, 0].set_xlabel("Infectious population")
axes[1, 0].set_ylabel("Density")
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.2)

# D. Final-time frequency histogram
axes[1, 1].hist(
    I[:, -1], bins=28,
    color="#79ae90", edgecolor="white"
)
axes[1, 1].axvline(
    np.median(I[:, -1]),
    color="#b42318", linestyle="--",
    label="Median"
)
axes[1, 1].set_title("D. Day-80 outcome frequencies")
axes[1, 1].set_xlabel("Infectious population")
axes[1, 1].set_ylabel("Number of paths")
axes[1, 1].legend()
axes[1, 1].grid(axis="y", alpha=0.2)

fig.suptitle("Stochastic SIS: paths and distributions")
fig.tight_layout()

Output

Run the code to see the result.

How to interpret the report

Panel A shows possible histories but only 15 of 500 paths.
Panel B uses every path and shows pointwise central spread around the mean.
Panel C shows that distribution shape can change, so mean and variance alone may be incomplete.
Panel D reports frequencies at one stated time rather than over the whole trajectory.
The printed table provides exact summaries supporting the visual claims.

Communication checklist

State model, initial conditions, parameter values, time step, number of paths and seed.
Put units and biological quantities on axes.
Say whether a band is percentile, SD or Monte Carlo error.
State whether summaries include all paths or only non-extinct paths.
Report numerical boundary treatment and correction frequency.
Use enough precision to reproduce results, but not false biological precision.

Biological interpretation

The ensemble describes intrinsic stochastic variation under fixed assumptions. Wide distributions indicate that the same biology permits substantially different prevalence levels. A mass at zero has direct extinction meaning only when the model and numerical boundary rule justify that interpretation.

SDE section completion

You can now derive epidemic diffusions from events, generate Brownian increments, program Euler–Maruyama, simulate one or many paths, calculate ensemble summaries, compare model classes, diagnose numerical validity, and communicate stochastic results accurately.