#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ When Most People Infect Nobody: What Overdispersion Does to an Outbreak Science Journaling Club, Volume 1 Issue 2, Winter 2025 THE QUESTION ------------ The average number of people an infectious case goes on to infect, R, is the number everybody quotes. It hides how unevenly that transmission is spread across cases. In a real outbreak most infected people infect nobody at all, and a small minority produce most of the onward chains, sometimes dozens of them at a single event. Two epidemics can share the same R and behave nothing alike. This script asks, holding R fixed, what changing only that unevenness does to whether an outbreak takes off, how big it gets if it goes extinct on its own, how concentrated its transmission is, and which of two control strategies -- lowering R everywhere, or specifically capping the largest transmission events -- gets more extinction for the same amount of transmission removed. THE MODEL --------- A Galton-Watson branching process. One index case starts generation zero. Every case in a generation independently produces a number of secondary cases drawn from a negative binomial distribution with mean R and dispersion parameter k: Var(offspring) = R + R^2 / k Small k means heavy overdispersion: most draws are 0, a few are large. As k -> infinity the negative binomial converges to Poisson(R), the classical homogeneous-mixing case with no superspreading at all. R is held fixed at 2.5 throughout the main sweep; k is swept from 0.01 (extreme overdispersion) to 10,000 (indistinguishable from Poisson) on a log grid of 17 values, chosen to include several dispersion estimates reported for real pathogens (k around 0.1-0.16 for SARS and SARS-CoV-2, k of tens or more for measles and influenza). The process has no depletion of susceptibles, no network structure, no spatial component and no time axis; a generation is a transmission generation, not a calendar interval. It is deliberately the simplest object that has an R and a variance around that R. Its companion article in this issue, "The Same Disease on Four Different Networks," studies the same underlying idea, contact heterogeneity, from the opposite end: fixed offspring distribution, varying contact network. This script fixes the contact structure implicitly (well-mixed, unlimited susceptible pool) and varies only the offspring distribution's shape. Because generation sizes only ever need to be summed, not enumerated individual by individual, the whole main sweep uses one algebraic shortcut: the sum of n independent NegBin(k, p) draws is itself NegBin(n*k, p) for fixed p. So an entire generation of a branching process can be advanced with a single call to a vectorised negative-binomial sampler, one per generation, regardless of how many individuals are alive in it. A run is stopped and classed EXTINCT the generation a lineage produces zero new cases, or ESCAPED once its cumulative case count reaches a large-outbreak cutoff (5,000 for the main sweep). No branching process with R > 1 that clears a cutoff that large is coming back down; a supercritical Galton-Watson process that has not gone extinct grows to infinity with probability 1 (Athreya & Ney 1972), so "escaped" and "did not go extinct" are the same event for practical purposes here. The one part of the study that needs individual-level detail is the control comparison, because capping a distribution's largest values is not a negative-binomial operation and breaks the generation-sum shortcut. That part samples every individual's offspring count separately (still fully vectorised across all currently-alive individuals in all replicates at once), at a smaller cutoff (1,000) chosen to keep the extra cost bounded. VALIDATION (the part that makes this science and not output) -------------------------------------------------------------- For a Galton-Watson process with offspring probability generating function G(s) = E[s^offspring], the extinction probability starting from one individual is the smallest non-negative root of G(s) = s. For the negative binomial offspring distribution used here, G(s) = (1 + (R/k) * (1 - s)) ** (-k) which has no closed algebraic root but is trivial to solve numerically: G is convex and increasing on [0, 1], so the fixed-point iteration q -> G(q), started at q = 0, increases monotonically to the smallest fixed point. Every row of Part 1 below prints this analytic root beside the simulated extinction fraction from 40,000 independent replicates, together with the Monte Carlo standard error and the resulting z-score, so agreement (or its absence) is stated in the same units as the noise. Part 1 also computes the k -> infinity limit by solving the Poisson analogue, G(s) = exp(-R(1-s)), the classical branching-process result, and confirms the largest-k row of the negative-binomial sweep has converged to it. ASSUMPTIONS AND LIMITATIONS, NAMED ----------------------------------- - No depletion of susceptibles. R is treated as constant regardless of how many people have already been infected, which only describes the early, exponential phase of a real epidemic, before immunity or behaviour change bite. Nothing here models a full epidemic curve. - No network structure, no repeated contacts, no clustering. Every secondary case is an independent, freshly chosen individual. The negative binomial's overdispersion is used here purely as a stand-in for "some cases transmit far more than others," and is agnostic about *why* -- biology (viral load, duration of shedding), behaviour (contact count), or setting (indoor gatherings, choirs, wards). Real superspreading mixes all three, and this model cannot and does not distinguish them. - Discrete, synchronous transmission generations, not calendar time. A "generation" here is a step in the chain of transmission, not a fixed number of days, so nothing in this script makes a claim about outbreak duration or doubling time. - The control-strategy comparison models "limiting large gatherings" as a deterministic cap on any one case's secondary-infection count, with a fractional randomised rounding so the resulting mean can be matched exactly to the mean-reduction strategy (documented at the point it is used, Part 3). Real gathering limits are patchy, evadable, and do not touch every large event; this is a best-case idealisation of that control, not a forecast of what any specific policy would achieve. - "Effort" for the two control strategies is defined identically as the fraction of total raw transmission (the uncontrolled mean R) removed. This is one honest way to make the comparison fair; it says nothing about which strategy is cheaper, more enforceable, or more politically feasible in the real world, none of which this script has any data on. - Every random number in this file comes from numpy's PCG64 generator, seeded once from the integer below and split by numpy's SeedSequence into independent streams per sweep cell, so cells cannot leak correlation into each other and the whole run is exactly reproducible. Master seed: 20251217. Python 3.12, numpy (no other numerical dependency). Expected runtime: well under two minutes on a single laptop core. """ import sys import time import numpy as np SEED = 20251217 R0 = 2.5 K_GRID = np.array([0.01, 0.02, 0.05, 0.1, 0.16, 0.3, 0.5, 1.0, 2.0, 5.0, 10.0, 30.0, 100.0, 300.0, 1000.0, 3000.0, 10000.0]) M_MAIN = 40000 # replicates per dispersion cell, main sweep CAP_MAIN = 5000 # cumulative-size cutoff classed as "escaped" MAXGEN_MAIN = 500 # safety cap on generations SHARE_N = 3000000 # sample size for the top-10%-of-transmission estimate BOOT_B = 2000 # bootstrap resamples for the variance standard error K0_CTRL = 0.1 # baseline dispersion for the control comparison EFFORT_GRID = np.round(np.arange(0.0, 0.651, 0.05), 3) M_CTRL = 20000 # replicates per (strategy, effort) cell CAP_CTRL = 1000 MAXGEN_CTRL = 500 CONV_K = 0.1 # which main-sweep row doubles as the convergence figure def rule(ch="-", n=78): print(ch * n) def section(title): print() rule("=") print(title) rule("=") # -------------------------------------------------------------------------- # Analytic extinction probability: smallest non-negative root of G(s) = s. # -------------------------------------------------------------------------- def pgf_nb(s, R, k): return (1.0 + (R / k) * (1.0 - s)) ** (-k) def pgf_poisson(s, R): return np.exp(-R * (1.0 - s)) def extinction_root_nb(R, k, tol=1e-15, maxit=200000): q = 0.0 for i in range(maxit): qn = pgf_nb(q, R, k) if abs(qn - q) < tol: return qn, i + 1 q = qn return q, maxit def extinction_root_poisson(R, tol=1e-15, maxit=200000): q = 0.0 for i in range(maxit): qn = pgf_poisson(q, R) if abs(qn - q) < tol: return qn, i + 1 q = qn return q, maxit # -------------------------------------------------------------------------- # Part 1 simulator: generation-vectorised branching process, exact NB. # Uses the closure property that a sum of n iid NegBin(k, p) draws is itself # NegBin(n*k, p) for fixed p, so one call per generation advances every # still-alive replicate regardless of how many individuals it holds. # -------------------------------------------------------------------------- def simulate_branching_fast(rng, R, k, M, cap, maxgen): p = k / (k + R) cs = np.ones(M, dtype=np.float64) # current-generation size total = np.ones(M, dtype=np.float64) # cumulative outbreak size active = np.ones(M, dtype=bool) gens_used = np.zeros(M, dtype=np.int64) gen = 0 while active.any() and gen < maxgen: idx = np.nonzero(active)[0] n_param = cs[idx] * k draws = rng.negative_binomial(n_param, p).astype(np.float64) cs[idx] = draws total[idx] += draws gens_used[idx] += 1 newly_extinct = idx[draws == 0.0] newly_escaped = idx[(draws > 0.0) & (total[idx] >= cap)] active[newly_extinct] = False active[newly_escaped] = False gen += 1 unresolved = int(active.sum()) if unresolved: # Never observed in testing; if it ever fires, count as escaped # rather than silently mis-classing a lineage that is still growing. total[active] = np.maximum(total[active], cap) extinct_mask = (~active) & (cs == 0.0) escaped_mask = ~extinct_mask return extinct_mask, escaped_mask, total, unresolved, gen # -------------------------------------------------------------------------- # Part 3 simulator: per-individual, because capping breaks the NB closure. # Every currently-alive individual across every still-active replicate draws # its own raw offspring count; values above the (possibly fractional) cap C # are replaced by floor(C) or floor(C)+1, chosen at random with probability # equal to the fractional part, so E[capped] hits C exactly on average. # -------------------------------------------------------------------------- def simulate_branching_capped(rng, R, k, cap_C, M, size_cap, maxgen): p = k / (k + R) C0 = int(np.floor(cap_C)) frac = cap_C - C0 sizes = np.ones(M, dtype=np.int64) # current-generation size per replicate total = np.ones(M, dtype=np.int64) # cumulative outbreak size per replicate active = np.ones(M, dtype=bool) gen = 0 while active.any() and gen < maxgen: idx = np.nonzero(active)[0] counts = sizes[idx] n_indiv = int(counts.sum()) rep_of = np.repeat(idx, counts) # which replicate each individual belongs to raw = rng.negative_binomial(k, p, size=n_indiv) over = raw > C0 n_over = int(over.sum()) if n_over: bump = rng.random(n_over) < frac capped_vals = np.where(bump, C0 + 1, C0) raw = raw.copy() raw[over] = capped_vals next_sizes = np.bincount(rep_of, weights=raw, minlength=M).astype(np.int64) sizes[idx] = next_sizes[idx] total[idx] += next_sizes[idx] newly_extinct = idx[next_sizes[idx] == 0] newly_escaped = idx[(next_sizes[idx] > 0) & (total[idx] >= size_cap)] active[newly_extinct] = False active[newly_escaped] = False gen += 1 unresolved = int(active.sum()) if unresolved: total[active] = np.maximum(total[active], size_cap) extinct_mask = (~active) & (sizes == 0) escaped_mask = ~extinct_mask return extinct_mask, escaped_mask, total, unresolved, gen def se_prop(p_hat, n): return np.sqrt(max(p_hat * (1.0 - p_hat), 0.0) / n) def bootstrap_var_se(x, B, rng): n = len(x) if n < 2: return float("nan") idx = rng.integers(0, n, size=(B, n)) resampled = x[idx] vars_ = resampled.var(axis=1, ddof=1) return float(vars_.std(ddof=1)) def mean_cap_target(k, R, C): """Exact E[min(X, C)] for X ~ NegBin(k, p), p = k/(k+R), for real C >= 0, via E[min(X,C)] = sum_{j=0}^{floor(C)-1} P(X > j) + frac(C) * P(X > floor(C)). P(X > j) computed from the NB survival function using the same PGF-free recursive pmf update used nowhere else in this script -- built once here because scipy is not a dependency of this file.""" p = k / (k + R) C0 = int(np.floor(C)) frac = C - C0 # pmf(0), then recursive ratio pmf(j+1)/pmf(j) = (j+k)/(j+1) * (1-p) pmf0 = p ** k pmf = pmf0 cdf = pmf0 total = 0.0 j = 0 while j < C0: total += 1.0 - cdf # P(X > j) pmf = pmf * (j + k) / (j + 1) * (1.0 - p) cdf += pmf j += 1 total += frac * (1.0 - cdf) # P(X > C0) return total def find_cap_for_effort(k, R, target_mean, lo=0.0, hi=200000.0, tol=1e-9, maxit=200): # mean_cap_target is monotonically increasing in C, from 0 at C=0 to R as # C -> infinity, so ordinary bisection applies. if target_mean >= R - 1e-12: return hi a, b = lo, hi fa = mean_cap_target(k, R, a) - target_mean for _ in range(maxit): mid = 0.5 * (a + b) fm = mean_cap_target(k, R, mid) - target_mean if abs(fm) < tol or (b - a) < 1e-9: return mid if (fa < 0) == (fm < 0): a, fa = mid, fm else: b = mid return 0.5 * (a + b) def nb_pmf_row(k, R, xmax): p = k / (k + R) pmf = np.empty(xmax + 1) pmf[0] = p ** k for j in range(xmax): pmf[j + 1] = pmf[j] * (j + k) / (j + 1) * (1.0 - p) return pmf def main(): T0 = time.time() master = np.random.SeedSequence(SEED) section("WHEN MOST PEOPLE INFECT NOBODY: OVERDISPERSION IN A BRANCHING-PROCESS OUTBREAK") print("Science Journaling Club, Volume 1 Issue 2, Winter 2025") print("Master seed : %d" % SEED) print("Mean reproduction number R : %.2f" % R0) print("Dispersion grid (k), %d values : %s" % (len(K_GRID), ", ".join("%g" % k for k in K_GRID))) print("Replicates per dispersion cell : %d" % M_MAIN) print("Large-outbreak cutoff (main) : %d cumulative cases" % CAP_MAIN) print("numpy : %s" % np.__version__) print("python : %s" % sys.version.split()[0]) # ------------------------------------------------------------------ # PART 1: main dispersion sweep # ------------------------------------------------------------------ section("PART 1. THE DISPERSION SWEEP: EXTINCTION, FINAL SIZE, CONCENTRATION") rows = [] conv_extinct_trace = None for i, k in enumerate(K_GRID): cell_seed = np.random.SeedSequence(entropy=SEED, spawn_key=(1, i)) rng = np.random.default_rng(cell_seed) q_analytic, iters = extinction_root_nb(R0, k) extinct_mask, escaped_mask, total, unresolved, gens_run = simulate_branching_fast( rng, R0, k, M_MAIN, CAP_MAIN, MAXGEN_MAIN) n_extinct = int(extinct_mask.sum()) p_sim = n_extinct / M_MAIN se = se_prop(p_sim, M_MAIN) diff = p_sim - q_analytic z = diff / se if se > 0 else float("nan") final_sizes = total[extinct_mask] mean_final = float(final_sizes.mean()) se_mean_final = float(final_sizes.std(ddof=1) / np.sqrt(len(final_sizes))) var_final = float(final_sizes.var(ddof=1)) se_var_final = bootstrap_var_se(final_sizes, BOOT_B, rng) # Top-10%-of-cases transmission share, from a large direct sample of # the marginal offspring distribution (every case independently # draws from the same NegBin(k, R0), so this needs no branching # simulation of its own). p_share = k / (k + R0) sample = rng.negative_binomial(k, p_share, size=SHARE_N) sample.sort() total_transmission = sample.sum(dtype=np.float64) top10_n = int(round(0.10 * SHARE_N)) top20_n = int(round(0.20 * SHARE_N)) top10_share = sample[-top10_n:].sum(dtype=np.float64) / total_transmission top20_share = sample[-top20_n:].sum(dtype=np.float64) / total_transmission zero_frac_sim = float((sample == 0).mean()) zero_frac_exact = p_share ** k rows.append(dict(k=k, q_analytic=q_analytic, p_sim=p_sim, se=se, diff=diff, z=z, n_extinct=n_extinct, unresolved=unresolved, mean_final=mean_final, se_mean_final=se_mean_final, var_final=var_final, se_var_final=se_var_final, top10_share=top10_share, top20_share=top20_share, zero_frac_sim=zero_frac_sim, zero_frac_exact=zero_frac_exact, gens_run=gens_run)) if abs(k - CONV_K) < 1e-9: conv_extinct_trace = extinct_mask.astype(np.float64) print() print("%-9s %10s %12s %10s %10s %8s | %10s %6s" % ("k", "q_analytic", "p_sim", "SE", "diff", "z", "n_extinct", "unresolv")) rule() max_abs_z = 0.0 for r in rows: print("%-9g %10.6f %12.6f %10.6f %10.6f %8.3f | %10d %6d" % (r["k"], r["q_analytic"], r["p_sim"], r["se"], r["diff"], r["z"], r["n_extinct"], r["unresolved"])) max_abs_z = max(max_abs_z, abs(r["z"])) print() print("Largest |z| across all %d dispersion cells: %.3f" % (len(rows), max_abs_z)) n_beyond_3 = sum(1 for r in rows if abs(r["z"]) > 3.0) print("Cells with |z| > 3: %d of %d (expect roughly %.2f by chance at this count)" % (n_beyond_3, len(rows), len(rows) * 0.0027)) # Poisson limit check q_pois, it_pois = extinction_root_poisson(R0) q_largek = rows[-1]["q_analytic"] print() print("Poisson (k -> infinity) analytic extinction root : %.8f (%d iterations)" % (q_pois, it_pois)) print("Largest-k NB row (k = %g) analytic root : %.8f" % (K_GRID[-1], q_largek)) print("Absolute difference : %.2e" % abs(q_pois - q_largek)) print("Relative difference : %.2e" % (abs(q_pois - q_largek) / q_pois)) section("PART 1b. FINAL SIZE AND TRANSMISSION CONCENTRATION") print("Final outbreak size is measured only among replicates that went extinct on") print("their own (a well-defined finite random variable); replicates that escaped") print("the %d-case cutoff are excluded from this table, not treated as size %d." % (CAP_MAIN, CAP_MAIN)) print() print("%-9s %8s %14s %14s %10s %10s %10s %10s" % ("k", "n_ext", "mean_final", "var_final", "top10%sh", "top20%sh", "P(X=0)sim", "P(X=0)exact")) rule() for r in rows: print("%-9g %8d %8.2f+/-%4.2f %10.1f+/-%6.1f %9.4f%% %9.4f%% %9.5f %10.5f" % (r["k"], r["n_extinct"], r["mean_final"], r["se_mean_final"], r["var_final"], r["se_var_final"], 100 * r["top10_share"], 100 * r["top20_share"], r["zero_frac_sim"], r["zero_frac_exact"])) # ------------------------------------------------------------------ # PART 2: Monte Carlo convergence, one representative dispersion value # ------------------------------------------------------------------ section("PART 2. MONTE CARLO CONVERGENCE (k = %g, reusing its Part 1 replicates)" % CONV_K) row_conv = next(r for r in rows if abs(r["k"] - CONV_K) < 1e-9) q_conv = row_conv["q_analytic"] checkpoints = np.unique(np.round(np.logspace(np.log10(100), np.log10(M_MAIN), 24)).astype(int)) print("Analytic root at k = %g: %.6f" % (CONV_K, q_conv)) print() print("%10s %12s %10s %10s" % ("n_used", "running_p", "run_SE", "z")) rule() conv_rows = [] for n in checkpoints: p_run = float(conv_extinct_trace[:n].mean()) se_run = se_prop(p_run, n) z_run = (p_run - q_conv) / se_run if se_run > 0 else float("nan") conv_rows.append((int(n), p_run, se_run, z_run)) print("%10d %12.6f %10.6f %10.3f" % (n, p_run, se_run, z_run)) # ------------------------------------------------------------------ # PART 3: matched-effort control comparison # ------------------------------------------------------------------ section("PART 3. CONTROL COMPARISON AT MATCHED EFFORT (k0 = %g)" % K0_CTRL) print("Baseline: R0 = %.2f, dispersion k0 = %.2f (order of magnitude reported for" % (R0, K0_CTRL)) print("SARS-CoV and SARS-CoV-2 by Lloyd-Smith et al. 2005 and Endo et al. 2020).") print("Effort e is the fraction of the UNCONTROLLED mean transmission R0 removed,") print("defined identically for both strategies so the comparison is at equal cost") print("in that one specific currency (see docstring for what this does not claim).") print() print("Strategy A, MEAN reduction : thin every case's transmission uniformly by") print(" retention (1-e). This is a binomial thinning of a NegBin(k0,.) variable,") print(" which is itself NegBin(k0,.) with the same k0 and mean R0*(1-e) -- thinning") print(" leaves the SHAPE of the heterogeneity untouched and only scales it down.") print("Strategy B, VARIANCE reduction : cap any one case's secondary infections at") print(" a threshold C(e), solved numerically so E[min(X,C)] = R0*(1-e) exactly,") print(" same target mean as Strategy A at the same e, by construction.") print() q0_analytic, _ = extinction_root_nb(R0, K0_CTRL) print("Uncontrolled (e=0) analytic extinction probability: %.6f" % q0_analytic) print() ctrl_rows = [] for j, e in enumerate(EFFORT_GRID): R_eff = R0 * (1.0 - e) seedA = np.random.SeedSequence(entropy=SEED, spawn_key=(2, j, 0)) seedB = np.random.SeedSequence(entropy=SEED, spawn_key=(2, j, 1)) rngA = np.random.default_rng(seedA) rngB = np.random.default_rng(seedB) qA_analytic, _ = extinction_root_nb(R_eff, K0_CTRL) if R_eff > 0 else (1.0, 0) extA, escA, totA, unresA, gensA = simulate_branching_fast( rngA, R_eff, K0_CTRL, M_CTRL, CAP_CTRL, MAXGEN_CTRL) if R_eff > 0 else ( np.ones(M_CTRL, dtype=bool), np.zeros(M_CTRL, dtype=bool), np.ones(M_CTRL), 0, 0) pA = float(extA.mean()) seA = se_prop(pA, M_CTRL) cap_C = find_cap_for_effort(K0_CTRL, R0, R_eff) extB, escB, totB, unresB, gensB = simulate_branching_capped( rngB, R0, K0_CTRL, cap_C, M_CTRL, CAP_CTRL, MAXGEN_CTRL) pB = float(extB.mean()) seB = se_prop(pB, M_CTRL) achieved_mean_B = mean_cap_target(K0_CTRL, R0, cap_C) diff = pA - pB se_diff = np.sqrt(seA ** 2 + seB ** 2) z_diff = diff / se_diff if se_diff > 0 else float("nan") ctrl_rows.append(dict(e=e, R_eff=R_eff, qA_analytic=qA_analytic, pA=pA, seA=seA, cap_C=cap_C, achieved_mean_B=achieved_mean_B, pB=pB, seB=seB, diff=diff, se_diff=se_diff, z_diff=z_diff, unresA=unresA, unresB=unresB)) print("%6s %8s %10s | %10s %8s %8s | %8s %10s %8s | %8s %6s" % ("effort", "R_eff", "qA_theory", "pA_sim", "SE_A", "cap_C", "meanB", "pB_sim", "SE_B", "A-B", "z")) rule() for r in ctrl_rows: print("%6.2f %8.4f %10.6f | %10.6f %8.6f %8.2f | %8.4f %10.6f %8.6f | %8.4f %6.2f" % (r["e"], r["R_eff"], r["qA_analytic"], r["pA"], r["seA"], r["cap_C"], r["achieved_mean_B"], r["pB"], r["seB"], r["diff"], r["z_diff"])) print() n_favors_A = sum(1 for r in ctrl_rows[1:] if r["diff"] > 0 and abs(r["z_diff"]) > 2) n_favors_B = sum(1 for r in ctrl_rows[1:] if r["diff"] < 0 and abs(r["z_diff"]) > 2) print("Cells (excluding e=0) where A beats B by >2 SE : %d" % n_favors_A) print("Cells (excluding e=0) where B beats A by >2 SE : %d" % n_favors_B) e_guarantee = 1.0 - 1.0 / R0 print("Effort at which Strategy A alone drives R_eff to 1 (guaranteed eventual") print("extinction regardless of dispersion): e = 1 - 1/R0 = %.4f" % e_guarantee) # ------------------------------------------------------------------ # PART 4: distribution snapshots for the figure (two matched-mean cases) # ------------------------------------------------------------------ section("PART 4. OFFSPRING DISTRIBUTION SNAPSHOTS (same mean R = %.2f, two k)" % R0) for k_show in (0.1, 10000.0): pmf = nb_pmf_row(k_show, R0, 10) print("k = %-8g P(X=0..10): " % k_show + ", ".join("%.4f" % v for v in pmf) + " [P(X>10) = %.4f]" % (1.0 - pmf.sum())) section("HEADLINE NUMBERS (for the article; all read directly from the tables above)") row_headline = next(r for r in rows if abs(r["k"] - 0.1) < 1e-9) row_poisson_end = rows[-1] row_mid_k1 = next(r for r in rows if abs(r["k"] - 1.0) < 1e-9) row_low = rows[0] print("R0 fixed at : %.2f" % R0) print("k = 0.1 (SARS/SARS-CoV-2-like) simulated extinction : %.4f +/- %.4f" % (row_headline["p_sim"], row_headline["se"])) print("k = 0.1 analytic root : %.6f" % row_headline["q_analytic"]) print("k = 10000 (near-Poisson) simulated extinction : %.4f +/- %.4f" % (row_poisson_end["p_sim"], row_poisson_end["se"])) print("k = 10000 analytic root vs true Poisson root, abs diff : %.2e" % abs(q_pois - q_largek)) print("k = 0.01 (extreme) simulated extinction : %.4f +/- %.4f" % (row_low["p_sim"], row_low["se"])) print("k = 0.1 top-10%%-of-cases transmission share : %.2f%%" % (100 * row_headline["top10_share"])) print("k = 0.1 top-20%%-of-cases transmission share : %.2f%%" % (100 * row_headline["top20_share"])) print("k = 0.1 P(a case infects exactly nobody) : %.2f%%" % (100 * row_headline["zero_frac_sim"])) print("k = 10000 P(a case infects exactly nobody) : %.2f%%" % (100 * row_poisson_end["zero_frac_sim"])) print("k = 1 top-10%% share : %.2f%%" % (100 * row_mid_k1["top10_share"])) ctrl_mid = ctrl_rows[6] print("Control comparison at effort e = %.2f: A - B = %.4f (z = %.2f)" % (ctrl_mid["e"], ctrl_mid["diff"], ctrl_mid["z_diff"])) section("DONE") print("Total replicates, main sweep : %d" % (M_MAIN * len(K_GRID))) print("Total replicates, control sweep : %d" % (M_CTRL * len(EFFORT_GRID) * 2)) print("Total wall-clock time : %.1f s" % (time.time() - T0)) print("Seed: %d. Every number above is reproducible by rerunning this file." % SEED) rule("=") if __name__ == "__main__": main()