Extinction and major-outbreak probabilities

Monte Carlo simulation becomes scientifically meaningful only after outcomes are defined precisely. We now define early fade-out, major outbreak and extinction by a chosen time, then estimate their probabilities and simulation uncertainty.

Why “extinction probability” is ambiguous

In a closed SIR model with no births or imported infections, every completed epidemic eventually reaches \(I=0\). Therefore:

\[ \Pr(\text{eventual extinction})=1. \]

This fact does not distinguish an introduction that disappears after infecting two people from an outbreak that infects most of the population and then ends.

Always specify whether “extinction” means eventual extinction, extinction before a threshold is reached, or extinction by a stated time.

Define the outcomes before running the simulations

Let \(Z\) be the completed epidemic size—the total number of people infected, including those infectious initially. Choose a major-outbreak threshold \(H\).

Major outbreak

\[Z\ge H.\]

The epidemic reaches or exceeds the pre-specified size threshold.

Early fade-out

\[Z<H.\]

Infection becomes extinct before reaching the major-outbreak threshold.

Extinct by time \(T\)

\[I(T)=0.\]

The absorbing extinction state has been reached no later than the stated time.

Which events form a partition?

For a completed closed SIR trajectory, exactly one of these is true:

\[ Z<H \qquad\text{or}\qquad Z\ge H. \]

Therefore:

\[ \Pr(\text{early fade-out}) + \Pr(\text{major outbreak}) =1. \]

“Extinct by day 20” is a different time-specific event. It may overlap with early fade-out, and it is not the complement of major outbreak.

Choosing the major-outbreak threshold

In the worked example, \(N=100\) and:

\[ H=20. \]

A major outbreak is therefore defined as at least 20 people infected, or 20% of the population.

This threshold is illustrative. In an applied study, it might be based on hospital demand, public-health workload, a percentage of the population, or a scientifically justified outbreak definition.

Choose the threshold before examining the simulated results.
State whether the threshold concerns final size, peak size, cumulative cases or another quantity.
Report sensitivity to reasonable alternative thresholds.

Changing \(H\) changes the event being measured, so it can change the estimated probability even when the epidemic model is unchanged.

Convert simulated outcomes into indicators

For simulation \(m\), define the major-outbreak indicator:

\[ Y_m= \begin{cases} 1,&Z_m\ge H,\\ 0,&Z_m<H. \end{cases} \]

The indicator turns a verbal event into data that can be counted. Its sample mean is the estimated probability:

\[ \widehat p = \frac{1}{M}\sum_{m=1}^{M}Y_m = \frac{\text{number of major outbreaks}}{M}. \]

The early-fade-out and time-specific extinction indicators are constructed in the same way.

Monte Carlo uncertainty

Repeating the entire experiment with another seed generally gives a slightly different estimate. This sampling variation is measured by the Monte Carlo standard error:

\[ \operatorname{SE}_{\mathrm{MC}}(\widehat p) = \sqrt{\frac{\widehat p(1-\widehat p)}{M}}. \]

Biological uncertainty

The model assigns multiple possible epidemic outcomes even when parameters are fixed.

Monte Carlo uncertainty

A finite number of simulations estimates the model probability imperfectly.

Running more simulations reduces Monte Carlo uncertainty but does not remove the stochastic variability inherent in the epidemic model.

Why four times as many runs halve the standard error

The standard error is proportional to \(1/\sqrt M\). If \(M\) is multiplied by 4:

\[ \frac{1}{\sqrt{4M}} = \frac{1}{2\sqrt M}. \]

Thus approximately four times as much simulation work is required to halve the Monte Carlo standard error.

A 95% interval for a simulated probability

The code reports a Wilson 95% interval. It is designed for a binomial proportion and behaves better than the simple \(\widehat p\pm1.96\operatorname{SE}\) interval when probabilities are near 0 or 1 or the simulation count is modest.

For \(x\) successes among \(M\) simulations and \(z=1.96\), the Wilson centre and half-width are:

\[ \text{centre} = \frac{\widehat p+z^2/(2M)}{1+z^2/M}, \] \[ \text{half-width} = \frac{z}{1+z^2/M} \sqrt{ \frac{\widehat p(1-\widehat p)}{M} + \frac{z^2}{4M^2} }. \]

This interval describes numerical uncertainty from a finite Monte Carlo sample under the chosen model. It does not include uncertainty in \(\beta\), \(\gamma\), initial conditions or model assumptions.

Worked scenario

Use the same SIR model as Lesson 12:

\[ N=100,\quad I_0=2,\quad \beta=0.3,\quad\gamma=0.1, \quad\Delta t=0.02\text{ day}. \]

Run \(M=1000\) trajectories up to 160 days. Define:

major outbreak: \(Z\ge20\);
early fade-out: \(Z<20\); and
extinct by day 20: extinction time \(\le20\) days.

Interactive Python laboratory

Click Run code to simulate the outcomes, estimate the three probabilities, calculate Monte Carlo uncertainty and examine sensitivity to the major-outbreak threshold.

Interactive PythonOutbreak probabilities and uncertainty

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

def simulate_sir_outcome(
    N, initial_I, beta, gamma, dt, max_days, rng
):
    S = N - initial_I
    I = initial_I
    R = 0

peak_I = I
    max_steps = int(max_days / dt)
    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

peak_I = max(peak_I, I)

if I == 0:
            extinction_time = (step + 1) * dt
            break

completed = extinction_time is not None
    final_size_so_far = N - S

return {
        "completed": completed,
        "final_size": final_size_so_far,
        "peak_I": peak_I,
        "extinction_time": extinction_time
    }

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

def probability_summary(name, indicator):
    successes = int(indicator.sum())
    total = len(indicator)
    estimate = successes / total

standard_error = np.sqrt(
        estimate * (1 - estimate) / total
    )

lower, upper = wilson_interval(successes, total)

return {
        "event": name,
        "count": successes,
        "simulations": total,
        "estimate": estimate,
        "MC standard error": standard_error,
        "95% lower": lower,
        "95% upper": upper
    }

# Model and Monte Carlo settings
N = 100
initial_I = 2
beta = 0.30
gamma = 0.10
dt = 0.02
max_days = 160

number_of_simulations = 1000
major_threshold = 20
time_cutoff = 20
seed = 31

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

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

rng = np.random.default_rng(seed)
rows = []

for simulation in range(number_of_simulations):
    result = simulate_sir_outcome(
        N, initial_I, beta, gamma, dt, max_days, rng
    )
    result["simulation"] = simulation
    rows.append(result)

outcomes = pd.DataFrame(rows)
incomplete_count = (~outcomes["completed"]).sum()

if incomplete_count > 0:
    raise ValueError(
        "Some paths are incomplete. Increase max_days "
        "before classifying final-size outcomes."
    )

# Define event indicators after using pre-chosen definitions
major_indicator = (
    outcomes["final_size"] >= major_threshold
)
fadeout_indicator = (
    outcomes["final_size"] < major_threshold
)
extinct_by_time_indicator = (
    outcomes["extinction_time"] <= time_cutoff
)

if not (major_indicator | fadeout_indicator).all():
    raise ValueError("An outcome was not classified.")

if (major_indicator & fadeout_indicator).any():
    raise ValueError("The two final-size events overlap.")

probability_table = pd.DataFrame([
    probability_summary(
        f"Major outbreak: Z >= {major_threshold}",
        major_indicator
    ),
    probability_summary(
        f"Early fade-out: Z < {major_threshold}",
        fadeout_indicator
    ),
    probability_summary(
        f"Extinct by day {time_cutoff}",
        extinct_by_time_indicator
    )
])

print("Definitions")
print(f"Major outbreak: final size >= {major_threshold}")
print(f"Early fade-out: final size < {major_threshold}")
print(f"Time-specific event: extinct by day {time_cutoff}")
print()
print("Monte Carlo probability estimates")
print(
    probability_table.round(5).to_string(index=False)
)
print()
print(
    "Fade-out + major-outbreak estimate: "
    f"{fadeout_indicator.mean() + major_indicator.mean():.6f}"
)

# Sensitivity to the chosen major-outbreak threshold
thresholds = np.array([10, 20, 30, 40, 50])
threshold_estimates = np.array([
    (outcomes["final_size"] >= threshold).mean()
    for threshold in thresholds
])

threshold_table = pd.DataFrame({
    "major threshold": thresholds,
    "estimated major-outbreak probability": threshold_estimates
})

print()
print("Threshold sensitivity")
print(threshold_table.round(5).to_string(index=False))

# Running estimate as simulations accumulate
running_count = np.cumsum(major_indicator.to_numpy())
running_numbers = np.arange(1, number_of_simulations + 1)
running_estimate = running_count / running_numbers

# Figure 1: estimates with Wilson intervals
estimates = probability_table["estimate"].to_numpy()
lower_errors = (
    estimates - probability_table["95% lower"].to_numpy()
)
upper_errors = (
    probability_table["95% upper"].to_numpy() - estimates
)

labels = ["Major outbreak", "Early fade-out", "Extinct by day 20"]
x = np.arange(len(labels))

plt.figure(figsize=(9, 4.9))
plt.bar(x, estimates, color=["#dc2626", "#2563eb", "#16a34a"])
plt.errorbar(
    x,
    estimates,
    yerr=[lower_errors, upper_errors],
    fmt="none",
    ecolor="black",
    capsize=6,
    linewidth=1.5
)
plt.xticks(x, labels)
plt.ylabel("Estimated probability")
plt.ylim(0, 1.05)
plt.title("Monte Carlo estimates with Wilson 95% intervals")
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

# Figure 2: sensitivity to threshold
plt.figure(figsize=(8.5, 4.8))
plt.plot(
    thresholds,
    threshold_estimates,
    marker="o",
    linewidth=2.2,
    color="#7c3aed"
)
plt.xlabel("Major-outbreak final-size threshold H")
plt.ylabel("Estimated Pr(Z >= H)")
plt.title("Major-outbreak probability depends on its definition")
plt.ylim(0, 1.05)
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 3: Monte Carlo convergence
plt.figure(figsize=(9, 4.8))
plt.plot(
    running_numbers,
    running_estimate,
    color="#f59e0b",
    linewidth=1.5
)
plt.axhline(
    major_indicator.mean(),
    color="black",
    linestyle="--",
    label="Estimate using all simulations"
)
plt.xlabel("Number of simulations used")
plt.ylabel("Running major-outbreak estimate")
plt.title("Monte Carlo estimate stabilises as simulations accumulate")
plt.ylim(0, 1.05)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the new probability code

Code	Meaning
`outcomes["final_size"] >= major_threshold`	Creates one Boolean major-outbreak indicator per simulation.
`indicator.sum()`	Counts `True` values because Python treats `True` as 1.
`indicator.mean()`	Calculates the proportion of `True` values—the Monte Carlo probability estimate.
`major_indicator \| fadeout_indicator`	Uses elementwise logical “or” to check that every completed outcome belongs to at least one category.
`major_indicator & fadeout_indicator`	Uses elementwise logical “and” to check that the categories do not overlap.
`np.cumsum(...)`	Calculates cumulative successes after 1, 2, 3, and later simulations.
`yerr=[lower_errors, upper_errors]`	Adds possibly unequal lower and upper Wilson interval lengths to the bar chart.

Read the probability table correctly

Column	Interpretation
`count`	Number of simulations satisfying the event definition.
`estimate`	Count divided by the total number of independent simulations.
`MC standard error`	Estimated numerical sampling variability of the Monte Carlo proportion.
`95% lower/upper`	Wilson interval endpoints for the model probability.

Report the event definition, number of simulations, estimate and uncertainty together. A bare probability without its definition or simulation count is incomplete.

Why incomplete paths must not be silently classified

If the horizon is reached while \(I>0\) and the current final size is below \(H\), the path might cross the threshold later. Calling it an early fade-out would be wrong.

The code therefore stops with an error if any trajectory is incomplete. Other defensible approaches include extending the horizon or reporting an explicit unclassified/censored category.

Interpret threshold sensitivity

As \(H\) increases, the event \(\{Z\ge H\}\) becomes more demanding. Therefore its probability cannot increase.

If conclusions change strongly across reasonable thresholds, report that sensitivity rather than presenting one threshold as uniquely correct.

Interpret the running estimate

The first few simulations can make the running probability jump substantially. Each additional trajectory has less influence as the accumulated sample grows.

Stabilisation is evidence about Monte Carlo precision, not proof that the epidemic model is biologically correct. A precisely simulated poor model remains a poor model.

Try these experiments

Change seed = 31. Compare the new estimate with the original Wilson interval.
Use only 100 simulations. Observe the wider intervals and more irregular running estimate.
Increase the simulations to 2,000 if your browser can complete the calculation comfortably.
Change major_threshold to 10, 30 or 40. State the changed event definition before interpreting the probability.
Change time_cutoff from 20 to 10 and then 40 days. The extinction-by-time probability should be non-decreasing with the cutoff.
Shorten max_days enough to produce incomplete paths and read the protective error message.

Biological interpretation

The major-outbreak estimate describes how frequently the model crosses a specified public-health-relevant size under fixed parameters and initial conditions. The early-fade-out estimate describes failure to reach that size because stochastic recovery and transmission events unfold differently across introductions.

These are conditional model probabilities. They change if parameters, population size, starting infections, time step, intervention assumptions or outcome definitions change.

What this lesson has added

Event definitions are fixed before probability estimation.
Early fade-out and major outbreak partition completed final-size outcomes.
Extinction by a stated time is a separate event.
Boolean indicators convert definitions into countable simulation data.
Monte Carlo standard error measures finite-simulation uncertainty.
Wilson intervals quantify uncertainty in estimated event probabilities.
Threshold sensitivity reveals dependence on the operational outbreak definition.
Incomplete trajectories must be extended or explicitly treated as unclassified.

Lesson 14 completes the DTMC pathway by choosing appropriate graphs, tables and written interpretations for trajectories, distributions and probability estimates.