Analytical functions

Closed-form solutions for offspring distribution properties — no simulation needed.

Extinction and epidemic probability

For a branching process with Negative Binomial offspring, the extinction probability is computed via fixed-point iteration on the probability generating function:

using EpiBranch
using Distributions
using DataFrames

# For Geometric offspring (k=1), extinction prob = 1/R
q = extinction_probability(3.0, 1.0)
println("Geometric(R=3): P(ext) = $(round(q, digits=4)) (analytical: $(round(1/3, digits=4)))")

Geometric(R=3): P(ext) = 0.3333 (analytical: 0.3333)

The epidemic probability is the complement of the extinction probability:

p = epidemic_probability(2.5, 0.16)
println("NegBin(R=2.5, k=0.16): P(epidemic) = $(round(p, digits=3))")

NegBin(R=2.5, k=0.16): P(epidemic) = 0.205

Effect of dispersion

Higher overdispersion (lower k) increases extinction probability — even with the same R:

for k in [0.01, 0.1, 0.5, 1.0, 10.0]
    q = extinction_probability(2.5, k)
    println("k = $(lpad(k, 5)): P(ext) = $(round(q, digits=3))")
end

┌ Warning: Assignment to `q` in soft scope is ambiguous because a global variable by the same name exists: `q` will be treated as a new local. Disambiguate by using `local q` to suppress this warning or `global q` to assign to the existing global variable.
└ @ analytical.md:35
k =  0.01: P(ext) = 0.984
k =   0.1: P(ext) = 0.861
k =   0.5: P(ext) = 0.558
k =   1.0: P(ext) = 0.4
k =  10.0: P(ext) = 0.144

Dispatch on distributions and models

Distribution objects and BranchingProcess models are also accepted:

println("Poisson(2.0):   P(ext) = $(round(extinction_probability(Poisson(2.0)), digits=4))")
println("NegBin(2, 0.5): P(ext) = $(round(extinction_probability(NegBin(2.0, 0.5)), digits=4))")

model = BranchingProcess(NegBin(2.5, 0.16), LogNormal(1.6, 0.5))
println("Model:          P(ext) = $(round(extinction_probability(model), digits=4))")
println("Model:          top 20% = $(round(proportion_transmission(model; prop_cases=0.2) * 100, digits=1))%")
println("Model:          chain size dist = $(chain_size_distribution(model))")

Poisson(2.0):   P(ext) = 0.2032
NegBin(2, 0.5): P(ext) = 0.6404
Model:          P(ext) = 0.7945
Model:          top 20% = 88.3%
Model:          chain size dist = GammaBorel{Float64}(k=0.16, R=2.5)

All analytical functions accept R/k values, Distribution objects, or BranchingProcess models.

Superspreading: proportion of transmission

What fraction of transmission is caused by the most infectious fraction of cases? This is the "80/20 rule" metric:

# SARS-CoV-2: R = 2.5, k = 0.16
prop = proportion_transmission(2.5, 0.16; prop_cases = 0.2)
println("Top 20% of cases cause $(round(prop * 100, digits=1))% of transmission")

Top 20% of cases cause 88.3% of transmission

The result depends almost entirely on the dispersion parameter k:

for k in [0.01, 0.1, 0.16, 0.5, 1.0, 10.0, 1000.0]
    prop = proportion_transmission(2.5, k; prop_cases = 0.2)
    println("k = $(lpad(k, 6)): top 20% → $(round(prop * 100, digits=1))%")
end

┌ Warning: Assignment to `prop` in soft scope is ambiguous because a global variable by the same name exists: `prop` will be treated as a new local. Disambiguate by using `local prop` to suppress this warning or `global prop` to assign to the existing global variable.
└ @ analytical.md:72
k =   0.01: top 20% → 100.0%
k =    0.1: top 20% → 95.1%
k =   0.16: top 20% → 88.3%
k =    0.5: top 20% → 65.0%
k =    1.0: top 20% → 52.2%
k =   10.0: top 20% → 29.5%
k = 1000.0: top 20% → 20.9%

As k → ∞ (no overdispersion), the top 20% cause ≈ 20% of transmission.

Chain size distributions

Analytical chain size distributions for specific offspring families:

# Borel distribution (chain sizes from Poisson offspring)
d = Borel(0.8)
println("Borel(0.8):")
println("  Mean chain size: $(round(mean(d), digits=2))")
for n in 1:5
    println("  P(size=$n) = $(round(pdf(d, n), digits=4))")
end

Borel(0.8):
  Mean chain size: 5.0
  P(size=1) = 0.4493
  P(size=2) = 0.1615
  P(size=3) = 0.0871
  P(size=4) = 0.0557
  P(size=5) = 0.0391

chain_size_distribution dispatches on the offspring type:

d_pois = chain_size_distribution(Poisson(0.8))
d_nb = chain_size_distribution(NegBin(0.8, 0.5))
println("Poisson(0.8) chain sizes: $(typeof(d_pois))")
println("NegBin(0.8, 0.5) chain sizes: $(typeof(d_nb))")

Poisson(0.8) chain sizes: Borel{Float64}
NegBin(0.8, 0.5) chain sizes: GammaBorel{Float64}

Validation against simulation

using StableRNGs

R, k = 1.5, 0.5
q_exact = extinction_probability(R, k)

model = BranchingProcess(NegBin(R, k), Exponential(5.0))
results = simulate_batch(model, 1000;
    sim_opts = SimOpts(max_cases = 10_000, max_generations = 200),
    rng = StableRNG(42),
)
q_sim = containment_probability(results)

println("R=$R, k=$k:")
println("  Analytical: $(round(q_exact, digits=4))")
println("  Simulated:  $(round(q_sim, digits=4))")

R=1.5, k=0.5:
  Analytical: 0.7676
  Simulated:  0.767