#!/usr/bin/env python3 """ Branching process extinction probabilities and Otter-Dwass cluster size distributions under negative binomial transmission overdispersion. MIT License. """ import json import math import pathlib import random def pgf_negbin(s, r0, k): """Probability generating function of NegBin(mean=r0, dispersion=k).""" return (1.0 + (r0 / k) * (1.0 - s)) ** (-k) def solve_extinction_prob(r0, k, tol=1e-14, max_iter=200): """ Solve the fixed-point equation q = G(q) for q in [0, 1) using Newton-Raphson. For r0 <= 1, the unique extinction probability is 1.0. """ if r0 <= 1.0: return 1.0 # Newton-Raphson starting at 0.1 q = 0.1 for _ in range(max_iter): # f(q) = q - (1 + (r0/k)*(1 - q))^(-k) # f'(q) = 1 - r0 * (1 + (r0/k)*(1 - q))^(-(k + 1)) term = 1.0 + (r0 / k) * (1.0 - q) g = term ** (-k) dg = r0 * (term ** (-(k + 1.0))) f = q - g df = 1.0 - dg step = f / df q_next = q - step if q_next < 0.0: q_next = q / 2.0 elif q_next >= 1.0: q_next = (q + 1.0) / 2.0 if abs(q_next - q) < tol: return q_next q = q_next return q def solve_extinction_poisson(r0, tol=1e-14, max_iter=200): """Solve q = exp(-r0 * (1 - q)) for Poisson(r0) transmission.""" if r0 <= 1.0: return 1.0 q = 0.1 for _ in range(max_iter): g = math.exp(-r0 * (1.0 - q)) dg = r0 * g f = q - g df = 1.0 - dg step = f / df q_next = q - step if q_next < 0.0: q_next = q / 2.0 elif q_next >= 1.0: q_next = (q + 1.0) / 2.0 if abs(q_next - q) < tol: return q_next q = q_next return q def critical_introductions(q1, alpha=0.05): """ Minimum number of independent index cases j* required to achieve at least 1 - alpha probability of a sustained epidemic: (q1)^j <= alpha. """ if q1 >= 1.0: return math.inf return math.ceil(math.log(alpha) / math.log(q1)) def otter_dwass_prob(s, r0, k): """ Exact probability that a branching process initiated by 1 case dies out with total final size S = s (Otter-Dwass / Lagrange inversion formula): P(S = s) = (1/s) * binom(s - 1 + s*k - 1, s - 1) * (k/(r0+k))^(s*k) * (r0/(r0+k))^(s-1) """ log_p = ( math.lgamma(s * k + s - 1.0) - math.lgamma(s) - math.lgamma(s * k) + (s * k) * math.log(k / (r0 + k)) + (s - 1.0) * math.log(r0 / (r0 + k)) ) return math.exp(log_p) / s def sample_poisson(lam, rng): """Exact Knuth Poisson sampler with zero and underflow guards.""" if lam <= 0.0: return 0 L = math.exp(-lam) if L >= 1.0: return 0 count = 0 p = 1.0 while p > L: count += 1 p *= rng.random() return count - 1 def sample_negbin(r0, k, rng): """Sample an integer from NegBin(r0, k) using Poisson(Gamma(k, r0/k)).""" lam = rng.gammavariate(k, r0 / k) return sample_poisson(lam, rng) def run_monte_carlo(r0, k, n_sims=50000, max_cases=2000, seed=42): """Simulate n_sims branching processes to estimate extinction probability.""" rng = random.Random(seed) extinct_count = 0 solitary_count = 0 for _ in range(n_sims): active = 1 total = 1 while active > 0 and total < max_cases: new_active = 0 for _ in range(active): new_active += sample_negbin(r0, k, rng) active = new_active total += active if active == 0: extinct_count += 1 if total == 1: solitary_count += 1 emp_q1 = extinct_count / n_sims se_q1 = math.sqrt(emp_q1 * (1.0 - emp_q1) / n_sims) emp_p_s1 = solitary_count / n_sims return emp_q1, se_q1, emp_p_s1 def main(): root = pathlib.Path(__file__).resolve().parent.parent results_dir = root / "results" results_dir.mkdir(parents=True, exist_ok=True) # Primary parameter regime: R0 = 2.50 r0_base = 2.50 k_sars = 0.10 k_geometric = 1.00 # 1. Solve extinction probabilities q1_sars = solve_extinction_prob(r0_base, k_sars) q1_geometric = solve_extinction_prob(r0_base, k_geometric) q1_poisson = solve_extinction_poisson(r0_base) # Verification: for k = 1.0, exact theoretical root is 1 / r0 exact_geometric = 1.0 / r0_base assert abs(q1_geometric - exact_geometric) < 1e-12, "Geometric root mismatch" # 2. Critical introduction thresholds for 95% epidemic certainty (alpha = 0.05) j_sars_95 = critical_introductions(q1_sars, alpha=0.05) j_geom_95 = critical_introductions(q1_geometric, alpha=0.05) j_poiss_95 = critical_introductions(q1_poisson, alpha=0.05) # 3. Additional reference regimes from literature # MERS-like (R0 = 1.20, k = 0.16) q1_mers = solve_extinction_prob(1.20, 0.16) j_mers_95 = critical_introductions(q1_mers, alpha=0.05) # Smallpox / Pneumonic Plague-like (R0 = 1.50, k = 0.50) q1_smallpox = solve_extinction_prob(1.50, 0.50) j_smallpox_95 = critical_introductions(q1_smallpox, alpha=0.05) # 4. Otter-Dwass distribution of finite outbreak sizes for SARS-like regime p_s1_sars = otter_dwass_prob(1, r0_base, k_sars) # Theoretical solitary probability is (k / (r0 + k))^k theory_p_s1 = (k_sars / (r0_base + k_sars)) ** k_sars assert abs(p_s1_sars - theory_p_s1) < 1e-12, "Otter-Dwass solitary prob mismatch" p_s1_fraction_of_extinctions = p_s1_sars / q1_sars # Cumulative sum over s = 1..500 otter_sum_500 = sum(otter_dwass_prob(s, r0_base, k_sars) for s in range(1, 501)) residual_500 = abs(q1_sars - otter_sum_500) assert residual_500 < 1e-5, f"Residual too large: {residual_500}" # 5. Monte Carlo validation n_sims = 50000 emp_q1, se_q1, emp_p_s1 = run_monte_carlo(r0_base, k_sars, n_sims=n_sims, max_cases=2000, seed=42) # Check within 3 standard errors z_score = abs(emp_q1 - q1_sars) / se_q1 assert z_score < 3.0, f"Monte Carlo deviation exceeds threshold: z={z_score}" out = { "r0_base": round(r0_base, 2), "k_sars": round(k_sars, 2), "k_geometric": round(k_geometric, 2), "confidence_pct": 95, "max_cluster_size": 500, # Claim C1 outputs "sars_q1": round(q1_sars, 6), "geometric_q1": round(q1_geometric, 6), "poisson_q1": round(q1_poisson, 6), "sars_j_star_95": int(j_sars_95), "geometric_j_star_95": int(j_geom_95), "poisson_j_star_95": int(j_poiss_95), # Comparative literature regimes "mers_q1": round(q1_mers, 6), "mers_j_star_95": int(j_mers_95), "smallpox_q1": round(q1_smallpox, 6), "smallpox_j_star_95": int(j_smallpox_95), # Claim C2 outputs (Otter-Dwass & Monte Carlo) "sars_p_s1": round(p_s1_sars, 6), "sars_p_s1_extinction_pct": round(p_s1_fraction_of_extinctions * 100.0, 2), "sars_sum_500": round(otter_sum_500, 6), "sars_sum_500_residual": round(residual_500, 8), "mc_sims": int(n_sims), "mc_empirical_q1": round(emp_q1, 6), "mc_se_q1": round(se_q1, 6), "mc_z_score": round(z_score, 3), # Table 1 summary rows "table_regimes": [ { "pathogen": "SARS-like (high superspreading)", "r0": 2.50, "k": 0.10, "extinction_prob": round(q1_sars, 6), "epidemic_prob": round(1.0 - q1_sars, 6), "j_star_95": int(j_sars_95) }, { "pathogen": "Geometric baseline (moderate variance)", "r0": 2.50, "k": 1.00, "extinction_prob": round(q1_geometric, 6), "epidemic_prob": round(1.0 - q1_geometric, 6), "j_star_95": int(j_geom_95) }, { "pathogen": "Poisson limit (homogeneous)", "r0": 2.50, "k": "inf", "extinction_prob": round(q1_poisson, 6), "epidemic_prob": round(1.0 - q1_poisson, 6), "j_star_95": int(j_poiss_95) }, { "pathogen": "MERS-like (near-critical overdispersed)", "r0": 1.20, "k": 0.16, "extinction_prob": round(q1_mers, 6), "epidemic_prob": round(1.0 - q1_mers, 6), "j_star_95": int(j_mers_95) }, { "pathogen": "Pneumonic plague / Smallpox-like", "r0": 1.50, "k": 0.50, "extinction_prob": round(q1_smallpox, 6), "epidemic_prob": round(1.0 - q1_smallpox, 6), "j_star_95": int(j_smallpox_95) } ] } (results_dir / "R1.json").write_text(json.dumps(out, indent=2) + "\n") print(f"Generated R1.json successfully. SARS q1={q1_sars:.6f}, j*(95%)={j_sars_95}, residual={residual_500:.2e}") if __name__ == "__main__": main()