Extinction probability and outbreak probability

Stochastic models are valuable because identical initial conditions can lead to early extinction or a large epidemic. Monte Carlo simulation estimates how often each clearly defined event occurs.

Define the event before estimating its probability

Start an SIR epidemic with one infectious person in a population of 200. We define a major outbreak as reaching at least 20 people ever infected. If the epidemic becomes extinct before that threshold is reached, we call it early extinction.

Major outbreak

\[A=\{N-S(t)\ge20\text{ for some }t\}.\]

Early extinction

\[A^c=\{I(t)=0\text{ before }N-S(t)\ge20\}.\]

Under this stopping experiment, these events are complementary, so \(\Pr(A)+\Pr(A^c)=1\).

The threshold is part of the question. “Outbreak probability” is incomplete unless the model, initial state, time window and outbreak criterion are stated.

Monte Carlo estimator

For run \(r\), define \(Y_r=1\) if a major outbreak occurs and \(Y_r=0\) otherwise. For \(M\) independent simulations:

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

This is an estimate, not the exact probability. Repeating the whole Monte Carlo experiment with a different seed produces a slightly different value.

Quantify simulation uncertainty

An estimated Monte Carlo standard error is

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

A simple approximate 95% Monte Carlo interval is \(\widehat p\pm1.96\operatorname{SE}(\widehat p)\). Increasing \(M\) reduces simulation error at the slow rate \(1/\sqrt M\).

Interactive Python laboratory

Run 3,000 independent threshold experiments. Each simulation stops as soon as it reaches the major-outbreak threshold or becomes extinct, because later events cannot change its classification.

Interactive PythonExtinction and outbreak probabilities

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

def reaches_major_outbreak(
    rng, N, I0, beta, gamma, outbreak_threshold
):
    S = N - I0
    I = I0
    ever_infected = I0

while I > 0 and ever_infected < outbreak_threshold:
        infection_rate = beta * S * I / N
        recovery_rate = gamma * I
        total_rate = infection_rate + recovery_rate

# The event type is enough for this threshold question.
        # Waiting times do not affect which boundary is reached first.
        if rng.random() < infection_rate / total_rate:
            S -= 1
            I += 1
            ever_infected += 1
        else:
            I -= 1

return ever_infected >= outbreak_threshold

# Model and event definition
N = 200
I0 = 1
beta = 0.25
gamma = 0.10
outbreak_threshold = 20

number_of_runs = 3000
seed = 314
rng = np.random.default_rng(seed)

major = np.zeros(number_of_runs, dtype=bool)

for run in range(number_of_runs):
    major[run] = reaches_major_outbreak(
        rng, N, I0, beta, gamma, outbreak_threshold
    )

major_count = int(major.sum())
extinction_count = number_of_runs - major_count

p_major = major.mean()
p_extinction = 1 - p_major
standard_error = np.sqrt(
    p_major * (1 - p_major) / number_of_runs
)
lower = max(0, p_major - 1.96 * standard_error)
upper = min(1, p_major + 1.96 * standard_error)

results = pd.DataFrame({
    "event": ["Major outbreak", "Early extinction"],
    "count": [major_count, extinction_count],
    "estimated probability": [p_major, p_extinction]
})

print("Event definition")
print(f"Initial infectious count: {I0}")
print(f"Major-outbreak threshold: {outbreak_threshold} ever infected")
print(f"Independent runs: {number_of_runs}")
print()
print(results.round(4).to_string(index=False))
print()
print(f"Monte Carlo SE for major-outbreak estimate: {standard_error:.4f}")
print(f"Approximate 95% Monte Carlo interval: [{lower:.4f}, {upper:.4f}]")
print(f"Probabilities sum to: {p_major + p_extinction:.4f}")

# Running estimate shows Monte Carlo convergence
run_numbers = np.arange(1, number_of_runs + 1)
running_probability = np.cumsum(major) / run_numbers

plt.figure(figsize=(9.5, 4.9))
plt.plot(run_numbers, running_probability,
         color="#174c34", linewidth=1.5)
plt.axhline(p_major, color="#b42318", linestyle="--",
            label=f"Final estimate = {p_major:.3f}")
plt.xlabel("Number of simulated trajectories")
plt.ylabel("Estimated major-outbreak probability")
plt.title("Monte Carlo estimate stabilises as runs accumulate")
plt.ylim(0, 1)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Probability display with Monte Carlo uncertainty
plt.figure(figsize=(7.8, 4.8))
probabilities = [p_major, p_extinction]
colours = ["#b42318", "#39815e"]

plt.bar(results["event"], probabilities, color=colours)
plt.errorbar(
    0, p_major, yerr=1.96 * standard_error,
    fmt="none", ecolor="black", capsize=6
)
plt.ylabel("Estimated probability")
plt.title("Major outbreak versus early extinction")
plt.ylim(0, 1)
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Why waiting times are omitted here

This specific question asks which boundary is reached first, not when it is reached. In a CTMC, the next event type depends on rate proportions, while the exponential waiting time only places that event on the clock. Omitting waiting times makes this classification experiment faster without changing the boundary reached.

Do not generalise this omission. Waiting times are essential for duration, peak time, probability by a specified date, intervention timing and time-dependent rates. Those questions require the full Gillespie clock.

Early-phase analytical reference

When the population is large and susceptible depletion is negligible, an SIR epidemic starting with \(I_0\) infectious people is approximately a birth–death branching process. If \(\beta>\gamma\), its eventual fade-out probability is approximately

\[q^{I_0}\approx\left(\frac{\gamma}{\beta}\right)^{I_0}.\]

This is not exactly the same as extinction before a finite threshold in a finite SIR population. Use it as early-phase intuition, not as a replacement for the stated Monte Carlo event.

Interpretation and sensitivity to definitions

Increasing the threshold usually reduces the estimated probability of crossing it.
Starting with more infectious people usually reduces early fade-out probability.
A probability estimate applies only to the parameter values and initial state used.
Monte Carlo uncertainty is different from uncertainty in biological parameters.

What this lesson adds

You can now define complementary epidemic events, estimate their probabilities from independent CTMC simulations, report counts as well as proportions, calculate Monte Carlo standard error, and distinguish a finite-threshold estimate from an early branching-process approximation. The next lesson measures continuous outcomes such as size, duration and peak.