Simulating many independent trajectories

One stochastic trajectory shows one possible epidemic. Monte Carlo simulation repeats the model independently to reveal the range, frequency and empirical distribution of possible outcomes.

Why one trajectory is insufficient

Suppose two simulations use the same biological parameters and initial state. Different random numbers can produce different:

early infection or recovery sequences;
peak infectious counts and peak times;
outbreak durations; and
final epidemic sizes.

A Monte Carlo experiment approximates stochastic-model outcomes by running many independent simulated trajectories and summarising the resulting sample.

The model does not change between runs. Only the random realisation changes.

Two levels of repetition

Inner loop: one trajectory

Repeatedly apply infection, recovery or no change through time until extinction or the observation horizon.

Outer loop: many trajectories

Call the complete-trajectory function again for simulation 0, simulation 1, and so on.

\[ \underbrace{X_0^{(m)},X_1^{(m)},\ldots,X_K^{(m)}}_ {\text{trajectory }m}, \qquad m=1,2,\ldots,M. \]

The superscript \((m)\) identifies the simulation. It is not a mathematical power.

The Monte Carlo workflow

1. Fix model inputsUse the same parameters and initial state.

2. Create one generatorSet one seed before the outer loop.

3. Simulate independentlyLet the generator advance between runs.

4. SummariseBuild empirical tables and distributions.

Independence and random seeds

Create the random generator once:

rng = np.random.default_rng(seed)

for simulation in range(number_of_simulations):
    result = simulate_one_trajectory(..., rng)

The generator supplies a continuing sequence of pseudo-random numbers. Different runs receive different portions of that sequence and are treated as independent simulations.

Do not recreate np.random.default_rng(seed) with the same seed inside the outer loop. That would restart the same sequence for every simulation and could produce identical trajectories rather than independent replicates.

Using the seed once makes the entire Monte Carlo experiment reproducible. Changing the seed creates a different reproducible sample.

Two kinds of stored results

Path arrays

Store \(S\), \(I\) and \(R\) at every common time. These arrays support trajectory graphs, pointwise means and percentile bands.

One-row summaries

Store one peak, peak time, final size and completion indicator per simulation. These values support empirical outcome tables and histograms.

Keeping these two data structures separate prevents confusion between “one value per time” and “one value per simulation.”

Why trajectories must be aligned

Different outbreaks become extinct at different times. To calculate quantities across simulations at time \(t_n\), every path array must use the same time grid and length.

After extinction in a closed SIR model, the state is absorbing. Therefore it is biologically correct to extend the remaining array with the final state:

\[ (S_n,I_n,R_n)=(S_\mathrm{final},0,R_\mathrm{final}) \quad\text{for later stored times}. \]

This is called absorbing-state padding. It aligns paths without inventing new epidemic events.

Empirical distributions

If the completed final sizes are \(Z_1,Z_2,\ldots,Z_M\), their empirical distribution is the observed pattern across the simulations.

For a possible value \(z\), its empirical relative frequency is:

\[ \widehat p(z) = \frac{\text{number of simulations with }Z_m=z}{M}. \]

A histogram groups nearby values into bins and displays how frequently different ranges occurred. It approximates the outcome distribution; it is not one epidemic trajectory.

Mean and percentile bands across paths

Pointwise meanAverage the \(M\) infectious counts separately at every time.

10th percentileAt each time, about 10% of simulated counts lie at or below this value.

90th percentileAt each time, about 90% of simulated counts lie at or below this value.

The region between the 10th and 90th percentiles describes the central 80% of the simulated counts at each time.

This percentile band describes variability among trajectories. It is not a confidence interval for an unknown parameter, and it does not represent one path that simultaneously follows the lower or upper boundary.

Scenario

Simulate \(M=400\) independent SIR trajectories in a population of 100:

\[ (S_0,I_0,R_0)=(98,2,0), \quad \beta=0.3, \quad \gamma=0.1, \quad \Delta t=0.02\text{ day}. \]

The common observation horizon is 160 days. The program reports whether any trajectory is still infectious at that time, because its completed final size would not yet be known.

Interactive Python laboratory

Click Run code. The calculation may take longer than a single-path lesson because it simulates millions of individual DTMC steps.

Interactive Python400 independent SIR trajectories

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

# Fixed-length aligned arrays
    S_path = np.empty(max_steps + 1, dtype=int)
    I_path = np.empty(max_steps + 1, dtype=int)
    R_path = np.empty(max_steps + 1, dtype=int)

S_path[0] = S
    I_path[0] = I
    R_path[0] = R

peak_I = I
    peak_step = 0
    extinction_step = 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

S_path[step + 1] = S
        I_path[step + 1] = I
        R_path[step + 1] = R

if I > peak_I:
            peak_I = I
            peak_step = step + 1

if I == 0:
            extinction_step = step + 1

# Pad the absorbing final state
            S_path[step + 2:] = S
            I_path[step + 2:] = 0
            R_path[step + 2:] = R
            break

completed = extinction_step is not None

if completed:
        end_time = extinction_step * dt
    else:
        end_time = max_days

return {
        "S_path": S_path,
        "I_path": I_path,
        "R_path": R_path,
        "peak_I": peak_I,
        "peak_time": peak_step * dt,
        "final_size": N - S,
        "completed": completed,
        "end_time": end_time
    }

# Monte Carlo settings
N = 100
initial_I = 2
beta = 0.30
gamma = 0.10
dt = 0.02
max_days = 160
number_of_simulations = 400
seed = 21

if number_of_simulations <= 0:
    raise ValueError("Number of simulations must be positive.")

if initial_I <= 0 or initial_I > N:
    raise ValueError("initial_I must be between 1 and N.")

# Sufficient whole-state-space probability check
upper_probability_bound = (
    beta * N / 4 + gamma * N
) * dt

if upper_probability_bound > 1:
    raise ValueError("dt is not guaranteed valid.")

max_steps = int(max_days / dt)
time_values = np.arange(max_steps + 1) * dt

# One generator is created before the outer loop
rng = np.random.default_rng(seed)

all_I_paths = np.empty(
    (number_of_simulations, max_steps + 1),
    dtype=int
)
summary_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"]

summary_rows.append({
        "simulation": simulation,
        "peak_I": result["peak_I"],
        "peak_time": result["peak_time"],
        "final_size": result["final_size"],
        "completed": result["completed"],
        "end_time": result["end_time"]
    })

summaries = pd.DataFrame(summary_rows)

# Monte Carlo summaries across aligned paths
mean_I = all_I_paths.mean(axis=0)
lower_I = np.percentile(all_I_paths, 10, axis=0)
upper_I = np.percentile(all_I_paths, 90, axis=0)

completed_summaries = summaries[summaries["completed"]]
number_incomplete = (~summaries["completed"]).sum()

descriptive_table = completed_summaries[
    ["peak_I", "peak_time", "final_size", "end_time"]
].describe().loc[["mean", "std", "min", "50%", "max"]]

print("Monte Carlo experiment")
print(f"Seed used once: {seed}")
print(f"Independent trajectories: {number_of_simulations}")
print(f"Common observation horizon: {max_days} days")
print(f"Incomplete trajectories at horizon: {number_incomplete}")
print()
print("First 10 trajectory summaries")
print(summaries.head(10).to_string(index=False))
print()
print("Descriptive outcomes for completed trajectories")
print(descriptive_table.round(3).to_string())

# Figure 1: selected paths, mean and percentile band
plt.figure(figsize=(10, 5.5))

for simulation in range(15):
    plt.step(
        time_values,
        all_I_paths[simulation],
        where="post",
        color="#94a3b8",
        alpha=0.35,
        linewidth=0.8
    )

plt.fill_between(
    time_values,
    lower_I,
    upper_I,
    color="#bfdbfe",
    alpha=0.7,
    label="10th to 90th percentile"
)
plt.plot(
    time_values,
    mean_I,
    color="#1d4ed8",
    linewidth=2.5,
    label="Monte Carlo mean"
)
plt.xlabel("Time (days)")
plt.ylabel("Infectious people")
plt.title("Many stochastic trajectories and their pointwise summary")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: empirical outcome distributions
fig, axes = plt.subplots(1, 2, figsize=(11, 4.7))

axes[0].hist(
    completed_summaries["final_size"],
    bins=np.arange(0.5, N + 1.5, 4),
    color="#16a34a",
    edgecolor="white"
)
axes[0].set_xlabel("Completed final epidemic size")
axes[0].set_ylabel("Number of trajectories")
axes[0].set_title("Empirical final-size distribution")
axes[0].grid(axis="y", alpha=0.25)

axes[1].hist(
    completed_summaries["peak_I"],
    bins=np.arange(0.5, N + 1.5, 3),
    color="#dc2626",
    edgecolor="white"
)
axes[1].set_xlabel("Peak infectious count")
axes[1].set_ylabel("Number of trajectories")
axes[1].set_title("Empirical peak-size distribution")
axes[1].grid(axis="y", alpha=0.25)

fig.tight_layout()

Output

Run the code to see the result.

Understand the new Monte Carlo code

Code	Meaning
`np.empty((M, K + 1), dtype=int)`	Reserves a rectangular array with one trajectory per row and one time point per column.
`all_I_paths[simulation]`	Selects the row assigned to one simulation.
`summary_rows.append({...})`	Adds one dictionary of outcome values for each completed call.
`extinction_step is not None`	Checks whether an extinction step was recorded.
`S_path[step + 2:] = S`	Fills every later position with the absorbing final susceptible count.
`all_I_paths.mean(axis=0)`	Averages down the simulation rows separately for every time column.
`np.percentile(..., axis=0)`	Calculates a pointwise percentile across simulations at each time.
`summaries[summaries["completed"]]`	Filters the table so completed-final-size summaries exclude unfinished trajectories.

Why use `axis=0`?

The aligned array has shape:

\[ (\text{number of simulations}) \times (\text{number of time points}). \]

Rows distinguish simulations and columns distinguish time. Therefore axis=0 collapses the simulation direction and leaves one result for every time column.

Using axis=1 would instead average each trajectory across its time points, answering a different question.

Interpret the figures correctly

The grey lines are 15 individual trajectories selected only to keep the figure readable.
The blue line is the pointwise mean across all 400 trajectories.
The blue band is the pointwise central 80% empirical range.
The histograms use one summary outcome from each completed trajectory.

A histogram bar height is a count of simulations in an outcome interval. It is not the number of infectious people at a particular time.

Computational cost

If \(M\) simulations each allow \(K\) time steps, the maximum number of one-step calculations grows approximately like:

\[ M\times K. \]

Doubling both the number of simulations and the number of time steps can make the work roughly four times larger. Early stopping reduces work when extinction occurs before the horizon.

Try these experiments

Change seed = 21. The empirical summaries change, but the biological model remains the same.
Reduce number_of_simulations to 25. Observe the less stable histograms and percentile band.
Increase it to 1,000 if your browser can run the larger experiment comfortably.
Set initial_I = 1. Look for stronger separation between small and large final sizes.
Shorten max_days. Examine the incomplete-trajectory count and explain why those rows are excluded from completed final-size summaries.
Move the generator creation inside the outer loop with the same seed. Observe why repeated identical random sequences are a mistake, then restore the correct code.

Biological interpretation

The collection of paths shows that identical starting conditions and parameters do not imply identical epidemic histories. Some introductions disappear early, while others produce larger peaks and final sizes.

The empirical distributions quantify this outcome variation. They do not yet define which outcomes count as extinction or a major outbreak. Those event definitions and their estimated probabilities are introduced carefully in Lesson 13.

What this lesson has—and has not—done

You can now:

distinguish the inner trajectory loop from the outer Monte Carlo loop;
generate independent, reproducible trajectories correctly;
align paths using absorbing-state padding;
store path arrays and one-row outcome summaries;
calculate pointwise means and empirical percentile bands; and
construct empirical peak and final-size distributions.

Not yet: this page has not chosen a major-outbreak threshold or estimated event probabilities and Monte Carlo standard errors. Lesson 13 defines those outcomes before estimating them.

Simulating many independent trajectories

Why one trajectory is insufficient

Two levels of repetition

Inner loop: one trajectory

Outer loop: many trajectories

The Monte Carlo workflow

Independence and random seeds

Two kinds of stored results

Path arrays

One-row summaries

Why trajectories must be aligned

Empirical distributions

Mean and percentile bands across paths

Scenario

Interactive Python laboratory

Understand the new Monte Carlo code

Why use axis=0?

Interpret the figures correctly

Computational cost

Try these experiments

Biological interpretation

What this lesson has—and has not—done

Why use `axis=0`?