Histograms, sample paths, probability plots and tables

The same stochastic results can be communicated well or badly. The correct display depends on the question: individual histories, outcome variability, event probability through time or precise numerical estimates.

The purpose and boundary of this lesson

This lesson introduces no new epidemic transitions, probability formulas or simulation algorithms. It reuses the DTMC and Monte Carlo methods from Lessons 1–13 and focuses on selecting, labelling and interpreting outputs.

A display should answer one clearly stated question. Trying to make one graph communicate trajectories, distributions, thresholds, estimates and uncertainty simultaneously usually reduces clarity.

Start with the question, not the graph type

What can individual epidemics look like?

Show several sample paths against time.

Use: step-path plot

How variable is one outcome?

Collect one value per trajectory, such as final size.

Use: histogram or empirical distribution

How does an event probability change with time?

Estimate the proportion extinct by each common time.

Use: probability curve

What exact values should be reported?

Present definitions, simulation count, estimates and uncertainty.

Use: numerical table

1. Sample-path plots

A sample-path plot displays a small selection of realised integer-valued trajectories:

\[ I_0^{(m)},I_1^{(m)},\ldots,I_K^{(m)}. \]

Use a step plot because the DTMC state changes at fixed time points. Display only enough paths to reveal variation without creating an unreadable mass of lines.

Label the time unit and compartment count.
State how many paths are shown and how many were simulated.
Use transparency so overlapping paths remain visible.
Do not describe one selected path as the predicted epidemic.

Example caption: “Fifteen of 600 simulated SIR DTMC trajectories. Each grey step line is one realisation from the same initial condition and parameter set; variation arises from random infection and recovery events.”

2. Histograms of trajectory outcomes

A histogram should use one comparable outcome per independent trajectory, such as:

completed final epidemic size;
peak infectious count;
peak time; or
extinction time.

The horizontal axis shows the outcome value or interval. The vertical axis must say whether bars represent counts, relative frequencies or probability density.

Histogram appearance depends on bin width. Very wide bins hide structure; very narrow bins create a noisy display. Check that the biological conclusion is not an artefact of one arbitrary bin choice.

If a major-outbreak threshold is relevant, mark it with a vertical line and identify it in the legend or caption.

3. Probability plots

Here “probability plot” means a probability shown against time or a threshold—not a statistical Q–Q plot.

Because extinction is absorbing in the closed SIR model:

\[ \widehat P(\text{extinct by }t_n) = \frac{1}{M} \sum_{m=1}^{M} \mathbf 1\{I_n^{(m)}=0\}. \]

The resulting curve must be non-decreasing. If it falls, either the event is not absorbing, the event definition changed, or the code is incorrect.

Keep the vertical probability axis between 0 and 1.
State the denominator \(M\).
State the initial condition and parameter scenario.
Do not call an estimated probability exact.

4. Numerical tables

Graphs reveal patterns; tables preserve exact reported values. A probability table should normally contain:

Field	Why it is needed
Event	Defines precisely what was counted.
Threshold or time	Prevents ambiguous terms such as “major” or “early.”
Count	Shows the number of simulations satisfying the event.
Total simulations	Provides the denominator.
Probability estimate	Reports count divided by total simulations.
Monte Carlo uncertainty	Shows the numerical precision of the simulation estimate.

Round values consistently, but retain enough digits to reflect the Monte Carlo precision. Reporting many more decimal places than the standard error supports gives false precision.

One scenario, four complementary outputs

The runnable example uses one shared dataset from 600 SIR simulations:

\[ N=100,\quad(S_0,I_0,R_0)=(98,2,0), \quad\beta=0.3,\quad\gamma=0.1, \quad\Delta t=0.02. \]

The major-outbreak threshold is 20 people ever infected. All outputs use the same simulations, ensuring that differences between displays arise from the questions they answer—not from different model settings.

Interactive Python laboratory

Click Run code to generate a communication-ready summary table and three deliberately separate figures.

Interactive PythonCommunicating DTMC results

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

def simulate_one_sir(
    N, initial_I, beta, gamma, dt, max_days, rng
):
    max_steps = int(max_days / dt)

S = N - initial_I
    I = initial_I
    R = 0

I_path = np.empty(max_steps + 1, dtype=int)
    I_path[0] = I

peak_I = I
    extinction_time = None

for step in range(max_steps):
        p_infection = beta * S * I / N * dt
        p_recovery = gamma * I * dt
        p_no_change = 1 - p_infection - p_recovery

if p_no_change < -1e-12:
            raise ValueError("dt gives invalid probabilities.")

u = rng.random()

if u < p_infection:
            S -= 1
            I += 1
        elif u < p_infection + p_recovery:
            I -= 1
            R += 1

I_path[step + 1] = I
        peak_I = max(peak_I, I)

if I == 0:
            extinction_time = (step + 1) * dt
            I_path[step + 2:] = 0
            break

return {
        "I_path": I_path,
        "peak_I": peak_I,
        "final_size": N - S,
        "extinction_time": extinction_time,
        "completed": extinction_time is not None
    }

def wilson_interval(successes, total, z=1.96):
    p_hat = successes / total
    denominator = 1 + z**2 / total
    centre = (
        p_hat + z**2 / (2 * total)
    ) / denominator
    half_width = (
        z / denominator
        * np.sqrt(
            p_hat * (1 - p_hat) / total
            + z**2 / (4 * total**2)
        )
    )
    return centre - half_width, centre + half_width

# Shared model and simulation settings
N = 100
initial_I = 2
beta = 0.30
gamma = 0.10
dt = 0.02
max_days = 160

number_of_simulations = 600
major_threshold = 20
seed = 41

max_steps = int(max_days / dt)
time_values = np.arange(max_steps + 1) * dt
rng = np.random.default_rng(seed)

all_I_paths = np.empty(
    (number_of_simulations, max_steps + 1),
    dtype=int
)
rows = []

for simulation in range(number_of_simulations):
    result = simulate_one_sir(
        N, initial_I, beta, gamma, dt, max_days, rng
    )

all_I_paths[simulation] = result["I_path"]

rows.append({
        "simulation": simulation,
        "peak_I": result["peak_I"],
        "final_size": result["final_size"],
        "extinction_time": result["extinction_time"],
        "completed": result["completed"]
    })

outcomes = pd.DataFrame(rows)

if not outcomes["completed"].all():
    raise ValueError(
        "Increase max_days before reporting completed outcomes."
    )

# Event definitions
major = outcomes["final_size"] >= major_threshold
early_fadeout = outcomes["final_size"] < major_threshold

# Probability of extinction by each common time
extinction_probability = (all_I_paths == 0).mean(axis=0)

if (np.diff(extinction_probability) < -1e-12).any():
    raise ValueError("Extinction probability decreased.")

# Communication-ready probability table
table_rows = []

for event_name, definition, indicator in [
    (
        "Major outbreak",
        f"Final size >= {major_threshold}",
        major
    ),
    (
        "Early fade-out",
        f"Final size < {major_threshold}",
        early_fadeout
    )
]:
    count = int(indicator.sum())
    estimate = count / number_of_simulations
    mc_se = np.sqrt(
        estimate * (1 - estimate) / number_of_simulations
    )
    lower, upper = wilson_interval(
        count, number_of_simulations
    )

table_rows.append({
        "event": event_name,
        "definition": definition,
        "count": count,
        "M": number_of_simulations,
        "estimate": estimate,
        "MC SE": mc_se,
        "95% interval": f"[{lower:.3f}, {upper:.3f}]"
    })

probability_table = pd.DataFrame(table_rows)

outcome_table = pd.DataFrame({
    "outcome": [
        "Final epidemic size",
        "Peak infectious count",
        "Extinction time (days)"
    ],
    "median": [
        outcomes["final_size"].median(),
        outcomes["peak_I"].median(),
        outcomes["extinction_time"].median()
    ],
    "10th percentile": [
        outcomes["final_size"].quantile(0.10),
        outcomes["peak_I"].quantile(0.10),
        outcomes["extinction_time"].quantile(0.10)
    ],
    "90th percentile": [
        outcomes["final_size"].quantile(0.90),
        outcomes["peak_I"].quantile(0.90),
        outcomes["extinction_time"].quantile(0.90)
    ]
})

print("Probability table")
print(probability_table.round(4).to_string(index=False))
print()
print("Outcome-distribution table")
print(outcome_table.round(2).to_string(index=False))

# Figure 1: selected sample paths
plt.figure(figsize=(9, 5))
for simulation in range(15):
    plt.step(
        time_values,
        all_I_paths[simulation],
        where="post",
        alpha=0.45,
        linewidth=1
    )
plt.xlabel("Time (days)")
plt.ylabel("Infectious people")
plt.title("Fifteen individual SIR DTMC sample paths")
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: final-size histogram with threshold
bins = np.arange(0.5, N + 1.5, 4)

plt.figure(figsize=(9, 4.8))
plt.hist(
    outcomes["final_size"],
    bins=bins,
    color="#39815e",
    edgecolor="white"
)
plt.axvline(
    major_threshold,
    color="#b91c1c",
    linestyle="--",
    linewidth=2,
    label=f"Major threshold = {major_threshold}"
)
plt.xlabel("Completed final epidemic size")
plt.ylabel("Number of trajectories")
plt.title("Empirical distribution of final epidemic size")
plt.legend()
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

# Figure 3: extinction probability through time
plt.figure(figsize=(9, 4.8))
plt.plot(
    time_values,
    extinction_probability,
    color="#2563eb",
    linewidth=2.5
)
plt.xlabel("Time (days)")
plt.ylabel("Estimated probability")
plt.title(
    f"Probability of extinction by time "
    f"({number_of_simulations} simulations)"
)
plt.ylim(0, 1.03)
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the communication code

Code	Communication purpose
`alpha=0.45`	Makes overlapping paths partly transparent so several remain visible.
`where="post"`	Displays fixed-step DTMC states without suggesting continuous fractional changes.
`bins=...`	Defines the outcome intervals counted by the histogram.
`plt.axvline(...)`	Marks the stated major-outbreak threshold on the outcome distribution.
`(all_I_paths == 0).mean(axis=0)`	Estimates extinction probability separately at every common time.
`plt.ylim(0, 1.03)`	Keeps a probability axis anchored to its meaningful 0-to-1 range.
`.quantile(0.10)`	Calculates an empirical lower percentile for the outcome table.

Common display mistakes

Mistake	Why it misleads	Correction
One path labelled “the prediction”	Hides stochastic variability.	Call it one realisation and show distributional summaries.
A smooth line through DTMC states	Suggests continuous fractional movement.	Use a step plot.
Histogram with no vertical-axis definition	Counts and densities have different meanings.	Label count, frequency or density explicitly.
Probability axis truncated near the estimate	Can exaggerate small differences.	Normally show the meaningful 0-to-1 probability range.
Threshold used but not displayed	The event definition is hidden.	State it in text, table and graph.
Too many coloured paths	Creates clutter without adding meaning.	Use a restrained selection and transparency.
Incomplete paths treated as final outcomes	Understates possible later infections.	Extend the horizon or report censoring.
Excessive decimal places	Implies unsupported precision.	Round consistently with the Monte Carlo error.

A five-part biological interpretation

1. ScenarioState the model, population, initial state and parameters.

2. Individual pathsDescribe visible variation without selecting a favourite path.

3. DistributionExplain the range, concentration or multiple outcome groups.

4. ProbabilityQuote the event definition, estimate and Monte Carlo uncertainty.

5. LimitationState what the model or simulation does not include.

Example interpretation: “For a closed population of 100 beginning with two infectious people, the simulations produced both early fade-outs and large outbreaks. Final epidemic sizes were therefore highly variable rather than concentrated around one deterministic outcome. Using a pre-specified threshold of 20 infections, the major-outbreak probability was estimated from 600 independent trajectories and reported with its Monte Carlo interval. These results are conditional on fixed parameters, homogeneous mixing and the one-event DTMC approximation.”

Final checklist before publishing a DTMC result

Is the biological scenario stated?
Are parameter values and units supplied?
Is the initial state supplied?
Is the time step stated and valid?
Is the number of simulations reported?
Is the random-seed policy reproducible?
Are paths distinguished from distributions?
Are all graph axes labelled with units?
Is a step plot used for state trajectories?
Is every threshold defined and displayed?
Are incomplete trajectories identified?
Is Monte Carlo uncertainty reported?
Are probability axes interpreted between 0 and 1?
Are tables rounded consistently?
Does the conclusion remain conditional on model assumptions?

What this lesson has added

Choose the display from the scientific question.
Use sample paths for individual histories.
Use histograms for empirical outcome distributions.
Use probability curves for event probability through time.
Use tables for precise definitions, estimates and uncertainty.
Keep trajectories, distribution summaries and probabilities conceptually separate.
Interpret results through the biological scenario and stated limitations.

This completes the DTMC programming pathway. The next section introduces continuous-time Markov chains, where event times are random rather than restricted to a fixed time grid.