""" extinction-by-bad-luck.py Science Journaling Club, Volume 1 Issue 1 (Fall 2024), theme "Populations and Chance". THE QUESTION ------------ A population whose average offspring number is above replacement is supposed to be safe. It is not. Every individual reproduces by drawing from a distribution, and a run of bad draws can take the whole population to zero before the average has any chance to assert itself. We ask two things: (1) How large must a population be before random reproductive variation stops killing it, where "stops killing it" is defined as extinction probability below 5 percent within a fixed horizon? (2) How does that threshold size move as the VARIANCE in offspring number grows, with the mean held exactly fixed? THE MODELS (both are computations, not observations) ---------------------------------------------------- Model A. Galton-Watson branching process in discrete generations. Every individual alive in generation g independently produces a random number of offspring drawn from a fixed distribution with mean m and then dies. We hold m = 1.15 across every run and change only the offspring distribution, so the five variance levels differ in nothing except spread: Poisson(m) variance = m = 1.15 Negative binomial r = 4 variance = m + m^2/4 = 1.48 Geometric (NB r = 1) variance = m + m^2 = 2.47 Negative binomial r = 1/4 variance = m + 4 m^2 = 6.44 Lottery: K = 30 with prob m/K, else 0 variance = m(K - m) = 33.18 A sum of n independent draws from each of these families is itself a draw from a single named distribution (Poisson(n m); negative binomial with shape n r; K times a binomial(n, m/K)), so a whole generation is produced with one vectorised call. That is an exact identity, not an approximation. Model B. Linear birth-death process in continuous time, overlapping generations. Each individual gives birth at per-capita rate lambda and dies at per-capita rate mu, both constant. The Malthusian growth rate r = lambda - mu is held fixed at 0.15 while the turnover lambda + mu is raised, which raises demographic variance (the variance of the change in population size accumulates at rate n(lambda + mu)) with the mean growth untouched. Simulated exactly by its embedded jump chain with exponential waiting times. VALIDATION ---------- For a branching process with offspring probability generating function f(s), the probability of eventual extinction starting from one individual is the smallest non-negative root of f(s) = s, and the probability of extinction within G generations starting from one individual is exactly the G-fold composition f_G(0). Starting from N0 independent individuals both quantities are raised to the power N0. So every simulated cell has an exact analytic partner, computed here two independent ways: * the finite-horizon value by iterating the generating function G times from zero; * the eventual-extinction root by bisection on f(s) - s, cross-checked against the closed form 1/m for geometric offspring and against the Lambert-W form q = -W(-m e^{-m})/m for Poisson offspring, with W found by its own Newton solve. For the birth-death process the probability of extinction by time t from n0 individuals has the closed form [ mu (1 - e^{-rt}) / (lambda - mu e^{-rt}) ]^{n0}, and the simulated value is printed beside it. Every comparison prints measured, analytic, difference, and the difference in units of the Monte Carlo standard error. ASSUMPTIONS, STATED SO THEY CAN BE ATTACKED ------------------------------------------- * Individuals reproduce independently of one another. No mate finding, no Allee effect, no crowding, no competition for a limited resource. * The offspring distribution never changes. No good years and bad years, which is environmental stochasticity; everything here is demographic stochasticity alone. * No spatial structure, no age structure, no sex, no immigration, no rescue effect. * The population is unbounded above. There is no carrying capacity, so survivors grow for ever at 15 percent a generation, which no real population does. * Generations do not overlap in Model A. Model B exists partly to show that the conclusion survives that assumption being dropped. LIMITATIONS OF THE COMPUTATION ITSELF ------------------------------------- * A replicate whose population reaches CAP individuals is recorded as a survivor and stopped. The bias this introduces is at most q^CAP per replicate, where q is the eventual extinction probability from one individual; the script prints that bound for every parameter set. It is below 1e-9 for every branching-process run and below 1e-6 for every birth-death run, in both cases far under the Monte Carlo error. * "Extinction probability" throughout means extinction within the stated horizon, not eventual extinction. For the highest-variance level those two differ noticeably and the script prints both. Run: python extinction-by-bad-luck.py Seed: 21091847 (numpy PCG64 via SeedSequence, one spawned child stream per cell) """ import math import sys import time import numpy as np SEED = 21091847 M_MEAN = 1.15 # mean offspring per individual per generation, Model A HORIZON_G = 100 # generations, Model A CAP_A = 20000 # escape threshold, Model A R_REPS = 20000 # replicates per cell R_CONV = 200000 # replicates for the convergence study Z95 = 1.959963984540054 START_SIZES = [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 24, 26, 28, 30, 33, 36, 40, 45, 50, 56, 63, 70, 80, 90, 100, 120, 150, 200, 250, 300, 400, 500] # ---------------------------------------------------------------------------------- # offspring families # ---------------------------------------------------------------------------------- class Offspring: """One offspring distribution with mean m, plus the total-offspring sampler.""" def __init__(self, key, label, kind, param): self.key = key self.label = label self.kind = kind self.param = param self.m = M_MEAN if kind == "poisson": self.var = self.m elif kind == "nb": r = param self.var = self.m + self.m * self.m / r elif kind == "lottery": K = param self.alpha = self.m / K self.var = self.m * (K - self.m) else: raise ValueError(kind) def pgf(self, s): """Probability generating function E[s^X] of a single individual's offspring.""" if self.kind == "poisson": return math.exp(self.m * (s - 1.0)) if self.kind == "nb": r = self.param return (r / (r + self.m - self.m * s)) ** r K = self.param return 1.0 - self.alpha + self.alpha * (s ** K) def total_offspring(self, rng, pop): """Sum of `pop[i]` independent offspring draws, for every replicate i at once.""" if self.kind == "poisson": return rng.poisson(self.m * pop) if self.kind == "nb": r = self.param p = r / (r + self.m) # numpy: mean = n(1-p)/p = pop*r*m/r = pop*m return rng.negative_binomial(r * pop, p) K = self.param return K * rng.binomial(pop, self.alpha) FAMILIES = [ Offspring("pois", "Poisson", "poisson", None), Offspring("nb4", "Neg. binomial r=4", "nb", 4.0), Offspring("geom", "Geometric (NB r=1)", "nb", 1.0), Offspring("nb025", "Neg. binomial r=0.25", "nb", 0.25), Offspring("lot30", "Lottery K=30", "lottery", 30), ] # ---------------------------------------------------------------------------------- # analytic machinery # ---------------------------------------------------------------------------------- def finite_horizon_single(fam, gens): """P(lineage from ONE individual is extinct by generation `gens`) = f_G(0), exact.""" u = 0.0 for _ in range(gens): u = fam.pgf(u) return u def extinction_root(fam): """Smallest non-negative root of f(s) = s, by bisection. Supercritical case.""" lo, hi = 0.0, 1.0 - 1e-14 # f(s) - s is positive at 0 (f(0) > 0) and negative just below 1 when m > 1 for _ in range(200): mid = 0.5 * (lo + hi) if fam.pgf(mid) - mid > 0.0: lo = mid else: hi = mid return 0.5 * (lo + hi) def lambert_w0(z): """Principal branch of Lambert W by Newton iteration, for z in (-1/e, 0).""" w = -0.5 if z < 0 else 0.5 for _ in range(200): ew = math.exp(w) f = w * ew - z dw = f / (ew * (w + 1.0) - (w + 2.0) * f / (2.0 * w + 2.0)) w -= dw if abs(dw) < 1e-16: break return w def wilson(k, n, z=Z95): if n == 0: return (0.0, 1.0) p = k / n d = 1.0 + z * z / n c = p + z * z / (2.0 * n) s = z * math.sqrt(p * (1.0 - p) / n + z * z / (4.0 * n * n)) return ((c - s) / d, (c + s) / d) def binom_se(p, n): return math.sqrt(max(p * (1.0 - p), 0.0) / n) def poisson_tail_p(k, expected): """Two-sided exact Poisson tail probability, for cells too rare for a z-score.""" if expected <= 0.0: return 1.0 if k == 0 else 0.0 term = math.exp(-expected) cdf = term below = term for i in range(1, max(k, 1) + 1): term *= expected / i cdf += term if i < k: below = cdf if k >= expected: upper = 1.0 - (cdf - term) return min(1.0, 2.0 * upper) return min(1.0, 2.0 * cdf) # ---------------------------------------------------------------------------------- # Model A: branching process # ---------------------------------------------------------------------------------- def run_branching(rng, fam, n0, reps, gens=HORIZON_G, cap=CAP_A): """Return array ext_gen[reps]: generation of extinction, or 0 if never extinct.""" pop = np.full(reps, n0, dtype=np.int64) idx = np.arange(reps) ext_gen = np.zeros(reps, dtype=np.int32) for g in range(1, gens + 1): pop = fam.total_offspring(rng, pop) dead = pop == 0 if dead.any(): ext_gen[idx[dead]] = g keep = ~dead & (pop < cap) if not keep.all(): idx = idx[keep] pop = pop[keep] if idx.size == 0: break return ext_gen # ---------------------------------------------------------------------------------- # Model B: continuous-time linear birth-death, exact jump-chain simulation # ---------------------------------------------------------------------------------- def run_birth_death(rng, lam, mu, n0, reps, horizon, cap): """Return (extinct bool array, extinction times array with nan for survivors).""" n = np.full(reps, n0, dtype=np.int64) t = np.zeros(reps, dtype=np.float64) idx = np.arange(reps) ext = np.zeros(reps, dtype=bool) ext_t = np.full(reps, np.nan, dtype=np.float64) pbirth = lam / (lam + mu) rate = lam + mu steps = 0 while idx.size: steps += 1 t = t + rng.exponential(1.0 / (n * rate)) over = t > horizon if over.any(): keep = ~over idx, n, t = idx[keep], n[keep], t[keep] if idx.size == 0: break birth = rng.random(idx.size) < pbirth n = n + np.where(birth, 1, -1) dead = n == 0 if dead.any(): ext[idx[dead]] = True ext_t[idx[dead]] = t[dead] keep = ~dead & (n < cap) if not keep.all(): idx, n, t = idx[keep], n[keep], t[keep] return ext, ext_t, steps def bd_analytic(lam, mu, t, n0): r = lam - mu e = math.exp(-r * t) single = mu * (1.0 - e) / (lam - mu * e) return single ** n0 # ---------------------------------------------------------------------------------- # weighted log-linear threshold estimate # ---------------------------------------------------------------------------------- def threshold_from_sim(sizes, phat, reps, target=0.05): """Weighted least squares of ln(p) on N over the usable window; return N*, SE.""" xs, ys, ws = [], [], [] for N, p in zip(sizes, phat): if 0.004 < p < 0.45: se_ln = binom_se(p, reps) / p xs.append(float(N)) ys.append(math.log(p)) ws.append(1.0 / (se_ln * se_ln)) if len(xs) < 3: return None, None, 0 x = np.array(xs) y = np.array(ys) w = np.array(ws) X = np.column_stack([np.ones_like(x), x]) W = np.diag(w) XtW = X.T @ W cov = np.linalg.inv(XtW @ X) beta = cov @ (XtW @ y) a, b = beta[0], beta[1] lt = math.log(target) nstar = (lt - a) / b # delta method on N* = (lt - a)/b da = -1.0 / b db = -(lt - a) / (b * b) var = (da * da * cov[0, 0] + 2 * da * db * cov[0, 1] + db * db * cov[1, 1]) return nstar, math.sqrt(max(var, 0.0)), len(xs) # ---------------------------------------------------------------------------------- # reporting helpers # ---------------------------------------------------------------------------------- def rule(ch="-", n=96): print(ch * n) def header(title): print() rule("=") print(title) rule("=") def main(): t_start = time.time() ss = np.random.SeedSequence(SEED) stream = iter(ss.spawn(4000)) print("EXTINCTION BY BAD LUCK") print("Science Journaling Club, Volume 1 Issue 1, Fall 2024") print("A population that grows on average can still die out.") print() print("master seed : %d" % SEED) print("numpy version : %s" % np.__version__) print("python : %s" % sys.version.split()[0]) print("replicates per cell : %d" % R_REPS) print("replicates, convergence run: %d" % R_CONV) print("mean offspring m : %.4f (held fixed across all variance levels)" % M_MEAN) print("horizon, Model A : %d generations" % HORIZON_G) print("escape cap, Model A : %d individuals" % CAP_A) # ---------------------------------------------------------------- section 1 header("SECTION 1. ANALYTIC EXTINCTION PROBABILITIES, AND THE CHECKS ON THEM") print("For each offspring law: q is the smallest non-negative root of f(s) = s,") print("found by bisection; q_G = f_G(0) is the exact probability that ONE lineage is") print("extinct by generation %d. Closed forms are solved independently where they exist." % HORIZON_G) print() print("%-22s %8s %8s %12s %12s %12s %10s" % ("offspring law", "mean", "var", "q (bisect)", "q_G=f_G(0)", "resid f(q)-q", "cap bias")) rule() ANA = {} for fam in FAMILIES: q = extinction_root(fam) qG = finite_horizon_single(fam, HORIZON_G) resid = fam.pgf(q) - q capbias = q ** CAP_A ANA[fam.key] = {"q": q, "qG": qG} print("%-22s %8.4f %8.4f %12.8f %12.8f %12.2e %10.1e" % (fam.label, fam.m, fam.var, q, qG, resid, capbias)) rule() print() print("How fast does extinction-within-G-generations approach eventual extinction?") print("f_G(0) for one starting individual, by horizon length:") print() print("%-22s %11s %11s %11s %11s %11s %14s" % ("offspring law", "G=5", "G=10", "G=20", "G=50", "G=100", "q - f_100(0)")) rule() HORIZ = {} for fam in FAMILIES: vals = [finite_horizon_single(fam, G) for G in (5, 10, 20, 50, 100)] HORIZ[fam.key] = vals print("%-22s %11.7f %11.7f %11.7f %11.7f %11.7f %14.2e" % (fam.label, vals[0], vals[1], vals[2], vals[3], vals[4], ANA[fam.key]["q"] - vals[4])) rule() print() print("Independent closed-form checks on the same roots:") print() print("%-26s %14s %14s %12s" % ("check", "club value", "closed form", "difference")) rule() geom = [f for f in FAMILIES if f.key == "geom"][0] q_geom_cf = 1.0 / M_MEAN print("%-26s %14.10f %14.10f %12.2e" % ("geometric: q = 1/m", ANA["geom"]["q"], q_geom_cf, ANA["geom"]["q"] - q_geom_cf)) pois = [f for f in FAMILIES if f.key == "pois"][0] z = -M_MEAN * math.exp(-M_MEAN) q_pois_cf = -lambert_w0(z) / M_MEAN print("%-26s %14.10f %14.10f %12.2e" % ("Poisson: q = -W(-m e^-m)/m", ANA["pois"]["q"], q_pois_cf, ANA["pois"]["q"] - q_pois_cf)) nb4 = [f for f in FAMILIES if f.key == "nb4"][0] # negative binomial shape r: f(s)=s has a root found here by bisection; verify by # substituting back at full double precision (printed above as resid). print("%-26s %14.10f %14.10f %12.2e" % ("geometric q_G vs 1/m limit", ANA["geom"]["qG"], q_geom_cf, ANA["geom"]["qG"] - q_geom_cf)) rule() print() print("The last line is negative by construction: q_G is extinction within %d" % HORIZON_G) print("generations and must sit below the eventual value q. At G = %d the gap has" % HORIZON_G) print("closed to less than a part in a million for every law we use, so the horizon is") print("not what sets our answer. At G = 20 it would be, which is why the table above") print("is printed rather than assumed.") # ---------------------------------------------------------------- section 2 header("SECTION 2. MODEL A, BRANCHING PROCESS: SIMULATED VS ANALYTIC, EVERY CELL") print("P_sim is the fraction of %d replicates extinct within %d generations." % (R_REPS, HORIZON_G)) print("P_ana = (f_G(0))^N0 is exact. z = (P_sim - P_ana) / SE, SE from the binomial.") print("CI is the 95 percent Wilson interval on P_sim.") RESULTS = {} all_z = [] rare_a = [] for fam in FAMILIES: print() print("--- %s (offspring variance %.4f, q = %.6f, q_G = %.6f) ---" % (fam.label, fam.var, ANA[fam.key]["q"], ANA[fam.key]["qG"])) print("%6s %10s %20s %10s %11s %7s %12s %10s" % ("N0", "P_sim", "95% CI", "P_ana", "diff", "z", "mean T_ext", "n_ext")) rule() rows = [] for n0 in START_SIZES: rng = np.random.default_rng(next(stream)) ext_gen = run_branching(rng, fam, n0, R_REPS) k = int((ext_gen > 0).sum()) p = k / R_REPS lo, hi = wilson(k, R_REPS) pa = ANA[fam.key]["qG"] ** n0 se = binom_se(pa, R_REPS) zz = (p - pa) / se if se > 0 else 0.0 te = ext_gen[ext_gen > 0] mt = float(te.mean()) if te.size else float("nan") mt_se = float(te.std(ddof=1) / math.sqrt(te.size)) if te.size > 1 else float("nan") rows.append({"n0": n0, "k": k, "p": p, "lo": lo, "hi": hi, "pa": pa, "z": zz, "mt": mt, "mt_se": mt_se, "next": int(te.size)}) if pa * R_REPS >= 10.0: all_z.append(zz) else: rare_a.append((fam.label, n0, k, pa * R_REPS, poisson_tail_p(k, pa * R_REPS))) print("%6d %10.5f [%8.5f,%8.5f] %10.5f %11.5f %7.2f %12s %10d" % (n0, p, lo, hi, pa, p - pa, zz, ("%.2f" % mt) if te.size else " -", te.size)) rule() RESULTS[fam.key] = rows # ---------------------------------------------------------------- section 3 header("SECTION 3. DID THE SIMULATOR PASS?") az = np.array(all_z) print("Every simulated cell has an exact analytic partner, so the z-scores above are") print("a sample from the standard normal if and only if the simulator is correct.") print("Cells where fewer than 10 extinctions are expected are held back from this") print("pool, because the normal approximation to a binomial is worthless there. They") print("are checked separately with an exact Poisson tail probability.") print() print("cells compared : %d" % az.size) print("mean z : %+.4f (expected 0, SE of the mean %.4f)" % (az.mean(), 1.0 / math.sqrt(az.size))) print("sd of z : %.4f (expected 1)" % az.std(ddof=1)) print("max |z| : %.3f" % np.abs(az).max()) print("cells with |z| > 1.96 : %d of %d (expected about %.1f)" % (int((np.abs(az) > 1.96).sum()), az.size, 0.05 * az.size)) print("cells with |z| > 3 : %d of %d (expected about %.1f)" % (int((np.abs(az) > 3.0).sum()), az.size, 0.0027 * az.size)) ks = float(np.abs(az).max()) worst = min(rare_a, key=lambda r: r[4]) if rare_a else None print("rare cells held back : %d (fewer than 10 extinctions expected)" % len(rare_a)) if worst: print("smallest Poisson tail p : %.4f (%s, N0 = %d: %d seen, %.2f expected)" % (worst[4], worst[0], worst[1], worst[2], worst[3])) print("rare cells with p < 0.001 : %d of %d" % (sum(1 for r in rare_a if r[4] < 0.001), len(rare_a))) print() print("Verdict: %s" % ("PASS, the measured curve is the analytic curve within Monte Carlo error." if abs(az.mean()) < 4.0 / math.sqrt(az.size) + 0.05 and ks < 5.0 else "FAIL, investigate before trusting anything below.")) # ---------------------------------------------------------------- section 4 header("SECTION 4. THE THRESHOLD: HOW BIG IS BIG ENOUGH FOR 5 PERCENT?") print("N* is the starting size at which extinction probability within %d generations" % HORIZON_G) print("falls to 0.05. Analytic: N* = ln(0.05) / ln(q_G). Simulated: weighted") print("least-squares fit of ln(P_sim) on N0 over the cells where P_sim is between") print("0.004 and 0.45, extrapolated to 0.05, with a delta-method standard error.") print("N_grid is the smallest size on our grid whose entire 95 percent CI sits below 0.05.") print() print("%-22s %9s %11s %11s %9s %10s %8s" % ("offspring law", "variance", "N* analytic", "N* sim", "SE", "diff", "N_grid")) rule() THRESH = {} for fam in FAMILIES: rows = RESULTS[fam.key] qG = ANA[fam.key]["qG"] n_ana = math.log(0.05) / math.log(qG) sizes = [r["n0"] for r in rows] ps = [r["p"] for r in rows] n_sim, n_se, npts = threshold_from_sim(sizes, ps, R_REPS) n_grid = None for r in rows: if r["hi"] < 0.05: n_grid = r["n0"] break THRESH[fam.key] = {"ana": n_ana, "sim": n_sim, "se": n_se, "grid": n_grid, "var": fam.var, "q": ANA[fam.key]["q"], "qG": qG, "ana_ultimate": math.log(0.05) / math.log(ANA[fam.key]["q"])} print("%-22s %9.4f %11.2f %11.2f %9.2f %10.2f %8s" % (fam.label, fam.var, n_ana, n_sim, n_se, n_sim - n_ana, str(n_grid) if n_grid else ">500")) rule() print() print("Same threshold, but for EVENTUAL extinction rather than extinction within the") print("horizon. This is where the choice of horizon shows up in the answer.") print() print("%-22s %9s %14s %16s %10s" % ("offspring law", "variance", "N* (100 gen)", "N* (eventual)", "ratio")) rule() for fam in FAMILIES: t = THRESH[fam.key] print("%-22s %9.4f %14.2f %16.2f %10.3f" % (fam.label, t["var"], t["ana"], t["ana_ultimate"], t["ana_ultimate"] / t["ana"])) rule() print() print("Threshold against variance, as a scaling: N* rises close to linearly in the") print("offspring variance once the variance is well above the mean.") print() print("%-22s %9s %13s %14s" % ("offspring law", "variance", "N* (100 gen)", "N*/variance")) rule() for fam in FAMILIES: t = THRESH[fam.key] print("%-22s %9.4f %13.2f %14.3f" % (fam.label, t["var"], t["ana"], t["ana"] / t["var"])) rule() xs = np.log(np.array([THRESH[f.key]["var"] for f in FAMILIES])) ys = np.log(np.array([THRESH[f.key]["ana"] for f in FAMILIES])) A = np.column_stack([np.ones_like(xs), xs]) coef, *_ = np.linalg.lstsq(A, ys, rcond=None) pred = A @ coef ssr = float(((ys - pred) ** 2).sum()) sst = float(((ys - ys.mean()) ** 2).sum()) print() print("log-log fit over the five levels: ln N* = %.4f + %.4f ln(variance), R^2 = %.5f" % (coef[0], coef[1], 1.0 - ssr / sst)) print("A slope of 1 would mean the threshold is exactly proportional to variance.") # ---------------------------------------------------------------- section 5 header("SECTION 5. TIME TO EXTINCTION, AMONG THE POPULATIONS THAT DIED") print("Conditioned on dying inside the horizon. Generations, mean and its standard") print("error, with quartiles and the 95th percentile of the same conditional") print("distribution. Blank rows are cells where fewer than 30 replicates died.") print() print("%-22s %6s %10s %8s %8s %8s %8s %8s %9s" % ("offspring law", "N0", "n_extinct", "mean", "SE", "q25", "median", "q75", "p95")) rule() TIMES = {} for fam in FAMILIES: TIMES[fam.key] = [] for n0 in [2, 5, 10, 20, 50, 100, 200, 500]: rng = np.random.default_rng(next(stream)) ext_gen = run_branching(rng, fam, n0, R_REPS) te = ext_gen[ext_gen > 0].astype(float) if te.size < 30: print("%-22s %6d %10d %8s %8s %8s %8s %8s %9s" % (fam.label, n0, te.size, "-", "-", "-", "-", "-", "-")) TIMES[fam.key].append({"n0": n0, "n": int(te.size), "mean": None}) continue mean = float(te.mean()) se = float(te.std(ddof=1) / math.sqrt(te.size)) q25, med, q75, p95 = [float(v) for v in np.percentile(te, [25, 50, 75, 95])] TIMES[fam.key].append({"n0": n0, "n": int(te.size), "mean": mean, "se": se, "q25": q25, "med": med, "q75": q75, "p95": p95}) print("%-22s %6d %10d %8.3f %8.3f %8.1f %8.1f %8.1f %9.1f" % (fam.label, n0, te.size, mean, se, q25, med, q75, p95)) rule() # ---------------------------------------------------------------- section 6 header("SECTION 6. CONVERGENCE OF THE MONTE CARLO ESTIMATE") print("One cell run to %d replicates, estimate recorded as trials accumulate." % R_CONV) print("halfwidth is the 95 percent normal-approximation halfwidth at that trial count.") CONV = {} conv_cells = [("geom", 20), ("pois", 10), ("lot30", 120)] checkpoints = [100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000, 100000, 200000] for key, n0 in conv_cells: fam = [f for f in FAMILIES if f.key == key][0] rng = np.random.default_rng(next(stream)) ext_gen = run_branching(rng, fam, n0, R_CONV) flags = (ext_gen > 0).astype(np.float64) csum = np.cumsum(flags) pa = ANA[key]["qG"] ** n0 print() print("--- %s, N0 = %d, analytic P = %.6f ---" % (fam.label, n0, pa)) print("%10s %12s %12s %12s %10s" % ("trials", "estimate", "halfwidth", "error", "err/SE")) rule() rowlist = [] for c in checkpoints: est = csum[c - 1] / c hw = Z95 * binom_se(est if 0 < est < 1 else pa, c) err = est - pa se = binom_se(pa, c) rowlist.append({"n": c, "est": est, "hw": hw, "err": err, "z": err / se}) print("%10d %12.6f %12.6f %+12.6f %10.2f" % (c, est, hw, err, err / se)) rule() CONV[key] = {"n0": n0, "pa": pa, "rows": rowlist} # ---------------------------------------------------------------- section 7 header("SECTION 7. MODEL B, CONTINUOUS-TIME BIRTH-DEATH WITH RISING TURNOVER") lam_mu = [(0.65, 0.50), (1.15, 1.00), (2.15, 2.00), (3.15, 3.00)] T_HORIZON = 100.0 bd_sizes = [2, 3, 5, 8, 12, 20, 30, 45, 60, 80] print("Every level has the same Malthusian growth rate r = lambda - mu = 0.15 per unit") print("time. Only the turnover lambda + mu changes, which is the demographic variance") print("rate per individual. Horizon t = %.0f, %d replicates per cell." % (T_HORIZON, R_REPS)) print("The escape cap is set per level so that the extinction probability still available") print("to a capped replicate is below 1e-6, and never below twice the largest start size.") print("P_ana = [mu(1 - e^{-rt}) / (lambda - mu e^{-rt})]^{N0}, exact.") BD = {} all_zb = [] rare_b = [] for lam, mu in lam_mu: q_ult = mu / lam CAP_B = max(2 * max(bd_sizes), int(math.ceil(math.log(1e-6) / math.log(q_ult)))) print() print("--- lambda = %.2f, mu = %.2f, turnover %.2f, eventual q = %.6f, cap %d, cap bias %.1e ---" % (lam, mu, lam + mu, q_ult, CAP_B, q_ult ** CAP_B)) print("%6s %10s %20s %10s %11s %7s %12s %10s" % ("N0", "P_sim", "95% CI", "P_ana", "diff", "z", "mean t_ext", "n_ext")) rule() rows = [] for n0 in bd_sizes: rng = np.random.default_rng(next(stream)) ext, ext_t, steps = run_birth_death(rng, lam, mu, n0, R_REPS, T_HORIZON, CAP_B) k = int(ext.sum()) p = k / R_REPS lo, hi = wilson(k, R_REPS) pa = bd_analytic(lam, mu, T_HORIZON, n0) se = binom_se(pa, R_REPS) zz = (p - pa) / se if se > 0 else 0.0 tt = ext_t[ext] mt = float(tt.mean()) if tt.size else float("nan") rows.append({"n0": n0, "p": p, "lo": lo, "hi": hi, "pa": pa, "z": zz, "mt": mt, "next": int(tt.size), "cap": CAP_B}) if pa * R_REPS >= 10.0: all_zb.append(zz) else: rare_b.append(("lam %.2f" % lam, n0, k, pa * R_REPS, poisson_tail_p(k, pa * R_REPS))) print("%6d %10.5f [%8.5f,%8.5f] %10.5f %11.5f %7.2f %12s %10d" % (n0, p, lo, hi, pa, p - pa, zz, ("%.3f" % mt) if tt.size else " -", tt.size)) rule() n_ana_b = math.log(0.05) / math.log(bd_analytic(lam, mu, T_HORIZON, 1)) sizes = [r["n0"] for r in rows] ps = [r["p"] for r in rows] n_sim_b, n_se_b, _ = threshold_from_sim(sizes, ps, R_REPS) BD[(lam, mu)] = {"rows": rows, "n_ana": n_ana_b, "n_sim": n_sim_b, "n_se": n_se_b, "q": q_ult, "turnover": lam + mu, "cap": CAP_B} print("threshold N* for 5 percent: analytic %.2f, simulated %.2f (SE %.2f)" % (n_ana_b, n_sim_b if n_sim_b else float("nan"), n_se_b if n_se_b else float("nan"))) azb = np.array(all_zb) worst_b = min(rare_b, key=lambda r: r[4]) if rare_b else None print() print("Model B validation, %d cells with at least 10 expected extinctions:" % azb.size) print("mean z %+.4f, sd %.4f, max |z| %.3f, |z| above 1.96 in %d cells." % (azb.mean(), azb.std(ddof=1), np.abs(azb).max(), int((np.abs(azb) > 1.96).sum()))) if worst_b: print("%d rarer cells checked by exact Poisson tail. Smallest p = %.4f (%s, N0 = %d:" % (len(rare_b), worst_b[4], worst_b[0], worst_b[1])) print("%d seen against %.2f expected). Rare cells with p < 0.001: %d." % (worst_b[2], worst_b[3], sum(1 for r in rare_b if r[4] < 0.001))) print() print("Model B thresholds against turnover:") print() print("%10s %10s %10s %12s %12s %10s" % ("lambda", "mu", "turnover", "eventual q", "N* analytic", "N* sim")) rule() for lam, mu in lam_mu: b = BD[(lam, mu)] print("%10.2f %10.2f %10.2f %12.6f %12.2f %10.2f" % (lam, mu, b["turnover"], b["q"], b["n_ana"], b["n_sim"])) rule() # ---------------------------------------------------------------- section 8 header("SECTION 8. SENSITIVITY: WHAT IF THE MEAN IS DIFFERENT?") print("Everything above fixes m = %.2f. Here is the analytic threshold N* for" % M_MEAN) print("extinction within %d generations, across a range of mean growth rates, for" % HORIZON_G) print("three of the offspring laws. Analytic only, no simulation needed.") print() means = [1.02, 1.05, 1.10, 1.15, 1.25, 1.50, 2.00] print("%10s %16s %16s %16s" % ("mean m", "Poisson N*", "Geometric N*", "Lottery K=30 N*")) rule() SENS = [] for mm in means: row = [mm] for kind, param in [("poisson", None), ("nb", 1.0), ("lottery", 30)]: f = Offspring("tmp", "tmp", kind, param) f.m = mm if kind == "lottery": f.alpha = mm / param qg = finite_horizon_single(f, HORIZON_G) row.append(math.log(0.05) / math.log(qg) if qg < 1.0 else float("inf")) SENS.append(row) print("%10.2f %16.2f %16.2f %16.2f" % (row[0], row[1], row[2], row[3])) rule() print() print("The threshold explodes as the mean approaches replacement. That is the honest") print("caveat on any single number we quote: N* depends on how far above 1 the mean is") print("at least as strongly as it depends on the variance.") # ---------------------------------------------------------------- section 8b header("SECTION 8B. AGAINST THE CLASSICAL APPROXIMATION (HALDANE 1927)") print("Haldane's branching-process argument for a rare advantageous type gives the") print("survival probability of a single lineage as roughly 2s / V, with s = m - 1 the") print("excess growth and V the offspring variance. It is derived for small s. Our") print("exact roots let us say how wrong it is at s = %.2f." % (M_MEAN - 1.0)) print() print("%-22s %9s %13s %13s %9s %11s %11s" % ("offspring law", "variance", "1 - q exact", "2s/V approx", "ratio", "N* exact", "N* approx")) rule() HALD = [] for fam in FAMILIES: q = ANA[fam.key]["q"] surv = 1.0 - q hald = 2.0 * (M_MEAN - 1.0) / fam.var n_exact = math.log(0.05) / math.log(q) n_appr = -math.log(0.05) / hald HALD.append((fam.key, fam.var, surv, hald, n_exact, n_appr)) print("%-22s %9.4f %13.6f %13.6f %9.4f %11.2f %11.2f" % (fam.label, fam.var, surv, hald, hald / surv, n_exact, n_appr)) rule() print() print("The approximation sits within %.1f to %.1f percent of the exact survival" % (min(abs(h[3] / h[2] - 1.0) for h in HALD) * 100, max(abs(h[3] / h[2] - 1.0) for h in HALD) * 100)) print("probability across a variance range of nearly thirty, which is why the linear") print("scaling of N* with variance that we measure was predictable from theory a") print("century old. What the simulation adds is the size of the error in that theory,") print("and the shape of the whole extinction curve rather than its rate alone.") # ---------------------------------------------------------------- section 9 header("SECTION 9. NUMBERS THE ARTICLE QUOTES") g = THRESH["geom"] p2 = RESULTS["geom"][0] print("Mean offspring fixed at : %.2f" % M_MEAN) print("Poisson threshold N* (5%%, 100 generations) : %.1f" % THRESH["pois"]["ana"]) print("Geometric threshold N* : %.1f" % THRESH["geom"]["ana"]) print("Lottery K=30 threshold N* : %.1f" % THRESH["lot30"]["ana"]) print("Ratio, highest to lowest variance threshold : %.2f" % (THRESH["lot30"]["ana"] / THRESH["pois"]["ana"])) print("Variance ratio over the same span : %.2f" % (THRESH["lot30"]["var"] / THRESH["pois"]["var"])) print("Extinction probability, geometric, N0 = 2 : %.4f" % p2["p"]) print("Extinction probability, geometric, N0 = 50 : %.4f" % [r for r in RESULTS["geom"] if r["n0"] == 50][0]["p"]) print("Extinction probability, Poisson, N0 = 20 : %.4f" % [r for r in RESULTS["pois"] if r["n0"] == 20][0]["p"]) print("Extinction probability, lottery, N0 = 100 : %.4f" % [r for r in RESULTS["lot30"] if r["n0"] == 100][0]["p"]) print("Largest |z| anywhere in Model A : %.3f over %d cells" % (np.abs(az).max(), az.size)) print("Largest |z| anywhere in Model B : %.3f over %d cells" % (np.abs(azb).max(), azb.size)) print("Runtime so far : %.1f s" % (time.time() - t_start)) # machine-readable block for the figures header("SECTION 10. MACHINE-READABLE BLOCK (used to draw the figures)") print("#BEGIN CURVES key,variance,N0,P_sim,lo,hi,P_ana") for fam in FAMILIES: for r in RESULTS[fam.key]: print("%s,%.4f,%d,%.6f,%.6f,%.6f,%.6f" % (fam.key, fam.var, r["n0"], r["p"], r["lo"], r["hi"], r["pa"])) print("#END CURVES") print("#BEGIN THRESH key,variance,N_ana,N_sim,SE,N_ultimate") for fam in FAMILIES: t = THRESH[fam.key] print("%s,%.4f,%.4f,%.4f,%.4f,%.4f" % (fam.key, t["var"], t["ana"], t["sim"], t["se"], t["ana_ultimate"])) print("#END THRESH") print("#BEGIN CONV key,n0,trials,estimate,halfwidth,analytic") for key, c in CONV.items(): for r in c["rows"]: print("%s,%d,%d,%.6f,%.6f,%.6f" % (key, c["n0"], r["n"], r["est"], r["hw"], c["pa"])) print("#END CONV") print("#BEGIN TIMES key,N0,n_extinct,mean,se,q25,median,q75,p95") for fam in FAMILIES: for r in TIMES[fam.key]: if r.get("mean") is None: continue print("%s,%d,%d,%.4f,%.4f,%.1f,%.1f,%.1f,%.1f" % (fam.key, r["n0"], r["n"], r["mean"], r["se"], r["q25"], r["med"], r["q75"], r["p95"])) print("#END TIMES") print("#BEGIN BD lambda,mu,turnover,N0,P_sim,lo,hi,P_ana,mean_t_ext") for (lam, mu), b in BD.items(): for r in b["rows"]: print("%.2f,%.2f,%.2f,%d,%.6f,%.6f,%.6f,%.6f,%.4f" % (lam, mu, b["turnover"], r["n0"], r["p"], r["lo"], r["hi"], r["pa"], r["mt"] if r["mt"] == r["mt"] else float("nan"))) print("#END BD") print("#BEGIN BDTHRESH lambda,mu,turnover,q,N_ana,N_sim,SE,cap") for (lam, mu), b in BD.items(): print("%.2f,%.2f,%.2f,%.6f,%.4f,%.4f,%.4f,%d" % (lam, mu, b["turnover"], b["q"], b["n_ana"], b["n_sim"], b["n_se"], b["cap"])) print("#END BDTHRESH") print("#BEGIN HORIZON key,G5,G10,G20,G50,G100,q") for fam in FAMILIES: v = HORIZ[fam.key] print("%s,%.8f,%.8f,%.8f,%.8f,%.8f,%.8f" % (fam.key, v[0], v[1], v[2], v[3], v[4], ANA[fam.key]["q"])) print("#END HORIZON") print("#BEGIN HALDANE key,variance,surv_exact,surv_haldane,N_exact,N_haldane") for h in HALD: print("%s,%.4f,%.6f,%.6f,%.4f,%.4f" % h) print("#END HALDANE") print("#BEGIN SENS mean,Poisson,Geometric,Lottery30") for row in SENS: print("%.2f,%.4f,%.4f,%.4f" % tuple(row)) print("#END SENS") print() rule("=") print("TOTAL RUNTIME: %.1f seconds" % (time.time() - t_start)) rule("=") if __name__ == "__main__": main()