""" Absence of Evidence: how many surveys before a run of blanks counts as evidence that a species is gone. Science Journaling Club, Volume 1 Issue 1, Fall 2024, "Populations and Chance". ------------------------------------------------------------------------------ THE QUESTION ------------------------------------------------------------------------------ A survey visit that records nothing has two possible explanations: the species is not there, or the species is there and the surveyor missed it. If the chance of detecting the species on a single visit to an occupied site is p, how many consecutive blank visits are needed before "we did not find it" becomes a defensible statement about the species rather than about the surveyor? ------------------------------------------------------------------------------ THE MODEL ------------------------------------------------------------------------------ This is the detection half of a single-season occupancy model (MacKenzie et al. 2002). The site is TAKEN AS OCCUPIED throughout. We are computing the sampling distribution of a detection history conditional on occupancy, which is the quantity that controls how much a run of blanks is worth. Model H (homogeneous). Each of k visits detects the species independently with the same probability p. The probability that all k visits come up blank at an occupied site is P_false-absence(k) = (1 - p)^k (1) and the visits needed for a false-absence rate of at most alpha = 0.05 is k* = ceil( ln(alpha) / ln(1 - p) ) (2) Model V (heterogeneity redrawn every visit). Detection probability for visit j is p_j ~ Beta(a, b), drawn afresh each visit, independent across visits. Because the visits are independent and the miss probability is linear in p, P_false-absence(k) = ( E[1 - p] )^k = (1 - pbar)^k (3) which is exactly Model H at p = pbar. Variation that reshuffles between visits costs nothing. We simulate this to confirm it rather than assert it. Model S (heterogeneity that persists at the site). One p ~ Beta(a, b) is drawn for the site and every visit to that site uses it. Now the miss probability is a nonlinear function of p and the expectation does not pass through: P_false-absence(k) = E[(1 - p)^k] = B(a, b + k) / B(a, b) (4) By Jensen's inequality (4) >= (1 - pbar)^k, with equality only at zero variance. We parameterise the Beta by its mean pbar and a concentration kappa = a + b, so a = kappa * pbar, b = kappa * (1 - pbar), and SD(p) = sqrt( pbar (1 - pbar) / (kappa + 1) ). Small kappa means strong heterogeneity. Equation (4) extended to real k via log-gamma is used to invert for k*. ------------------------------------------------------------------------------ VALIDATION (the part that makes this science rather than output) ------------------------------------------------------------------------------ 1. Every simulated false-absence rate in Model H is printed beside the closed form (1 - p)^k with the signed difference and the difference in units of the Monte Carlo standard error. 2. Model S is checked three ways: simulation, the Beta-function closed form (4), and an independent numerical integration of E[(1-p)^k] over the Beta density by tanh-sinh (double exponential) quadrature, which uses no gamma function at all. Closed form and quadrature must agree to ~1e-13; simulation must agree with both inside Monte Carlo error. The first quadrature we wrote was wrong and is documented at the point of failure rather than quietly deleted. 3. A convergence run tracks the running estimate against the exact value as replicates accumulate, with a +/- 2 SE envelope. ------------------------------------------------------------------------------ ASSUMPTIONS, STATED PLAINLY ------------------------------------------------------------------------------ A1. The site is occupied. Everything here is conditional on that. The unconditional question ("is the species present?") needs a prior on occupancy; we add that as a Bayesian postscript rather than smuggling it in. A2. The population does not change during the survey (closure). A3. Visits are independent given p. No learning by the surveyor, no animal becoming trap-shy, no weather run that ruins a whole fortnight. A4. Detection probability does not depend on abundance, season, time of day, observer, or on how many visits have already happened. A5. No false positives. A detection is always real. A6. In Model S the heterogeneity is a Beta distribution, which is a convenient choice, not a measured one. ------------------------------------------------------------------------------ LIMITATIONS ------------------------------------------------------------------------------ L1. This is a computation, not a survey. The club has no field site and made no observations. Every number below comes from arithmetic and from pseudo-random numbers with a stated seed. L2. p is treated as known. In real work p is itself estimated, usually from the same small dataset, and the uncertainty in p propagates into k* in a way we do not model here. L3. Real surveys violate A3 badly. Correlated visits (a bad season) behave more like Model S than Model H, which is one reason Model S matters. L4. The k* values for strongly heterogeneous cases run to millions of visits. They are correct arithmetic for the model and should be read as "this model says the question cannot be settled by more effort", not as a survey plan. L5. Beta heterogeneity is one shape. A two-point mixture (some sites easy, some nearly undetectable) gives a different and generally worse answer. Runtime: about 30 seconds. Python 3.12, numpy 2.x. """ from __future__ import annotations import math import sys import time import numpy as np # ---------------------------------------------------------------------------- # Reproducibility # ---------------------------------------------------------------------------- MASTER_SEED = 20240921 ALPHA = 0.05 # target false-absence rate N_REP = 500_000 # replicates per configuration (>= 50000) N_CONV = 2_000_000 # replicates in the convergence study CHUNK = 25_000 _ss = np.random.SeedSequence(MASTER_SEED) def new_rng(label: str) -> np.random.Generator: """One independent stream per configuration, all descended from MASTER_SEED.""" return np.random.default_rng(_ss.spawn(1)[0]) def rule(ch: str = "-", n: int = 78) -> str: return ch * n def head(title: str) -> None: print() print(rule("=")) print(title) print(rule("=")) # ---------------------------------------------------------------------------- # Analytic pieces # ---------------------------------------------------------------------------- def p_miss_hom(p: float, k: int) -> float: return (1.0 - p) ** k def k_star_hom_int(p: float, alpha: float = ALPHA) -> int: return int(math.ceil(math.log(alpha) / math.log(1.0 - p))) def k_star_hom_cont(p: float, alpha: float = ALPHA) -> float: return math.log(alpha) / math.log(1.0 - p) def beta_ab(pbar: float, kappa: float) -> tuple[float, float]: return kappa * pbar, kappa * (1.0 - pbar) def beta_sd(pbar: float, kappa: float) -> float: return math.sqrt(pbar * (1.0 - pbar) / (kappa + 1.0)) def log_p_miss_site(a: float, b: float, k: float) -> float: """log E[(1-p)^k] for p ~ Beta(a,b): log B(a, b+k) - log B(a, b).""" return (math.lgamma(a + b) - math.lgamma(b) + math.lgamma(b + k) - math.lgamma(a + b + k)) def p_miss_site(a: float, b: float, k: float) -> float: return math.exp(log_p_miss_site(a, b, k)) # --- independent numerical integration of E[(1-p)^k] ------------------------- # # FIRST ATTEMPT, WHICH FAILED AND IS RECORDED HERE ON PURPOSE. # The original quadrature split [0,1] at 1/2 and substituted p = (1/2)u^(1/a) # on the left and 1-p = (1/2)v^(1/(b+k)) on the right, then used 240-node # Gauss-Legendre on each half. That removes the endpoint singularity when the # exponent is below 1, but when the exponent is ABOVE 1 the substitution puts # an infinite derivative at the other end of the interval, and polynomial # quadrature cannot see it. The symptom was a disagreement with the closed form # of up to 8.6e-01 (!), entirely in the kappa = 100 rows where b is large. # Replaced with tanh-sinh (double exponential) quadrature below, which handles # algebraic endpoint singularities of either sign without any case analysis. # _TS_H = 0.001953125 # step in the doubly-transformed variable _TS_T = 5.2 # truncation _TS_t = np.arange(-_TS_T, _TS_T + 1e-12, _TS_H) _TS_u = (math.pi / 2.0) * np.sinh(_TS_t) _TS_logcosh = np.logaddexp(_TS_t, -_TS_t) - math.log(2.0) _TS_logp = -np.logaddexp(0.0, -2.0 * _TS_u) # log p at each node _TS_log1mp = -np.logaddexp(0.0, 2.0 * _TS_u) # log (1-p) at each node def log_beta_quad(a: float, c: float) -> float: """ log of Integral_0^1 p^(a-1) (1-p)^(c-1) dp by tanh-sinh quadrature. Substituting p = 1 / (1 + exp(-pi sinh t)) gives dp = pi cosh(t) p (1-p) dt, so the integrand becomes pi cosh(t) p^a (1-p)^c on t in (-inf, inf) and decays doubly exponentially at both ends. The whole sum is done in logs and shifted by its maximum so nothing overflows or underflows. No gamma function is used anywhere in this routine. """ L = math.log(math.pi) + _TS_logcosh + a * _TS_logp + c * _TS_log1mp M = float(L.max()) return M + math.log(float(np.sum(np.exp(L - M))) * _TS_H) def p_miss_site_quad(a: float, b: float, k: float) -> float: """ E[(1-p)^k] for p ~ Beta(a,b) by numerical integration, computed with no reference at all to the Beta-function identity: E = Integral p^(a-1)(1-p)^(b+k-1) dp / Integral p^(a-1)(1-p)^(b-1) dp Both integrals come from the same tanh-sinh routine, so the normalising constant is integrated too rather than taken from lgamma. """ return math.exp(log_beta_quad(a, b + k) - log_beta_quad(a, b)) def k_star_site_cont(a: float, b: float, alpha: float = ALPHA) -> float: """Smallest real k with E[(1-p)^k] <= alpha. Returns inf if unreachable.""" lo, hi = 0.0, 1.0 target = math.log(alpha) if log_p_miss_site(a, b, lo) <= target: return 0.0 ok = False for _ in range(400): if log_p_miss_site(a, b, hi) <= target: ok = True break hi *= 2.0 if hi > 1e18: return float("inf") if not ok: return float("inf") for _ in range(200): mid = 0.5 * (lo + hi) if log_p_miss_site(a, b, mid) <= target: hi = mid else: lo = mid return hi def fmt_k(x: float) -> str: if not math.isfinite(x): return "unreachable" if x < 1e5: return f"{x:,.1f}" return f"{x:.3e}" # ---------------------------------------------------------------------------- # Simulation core # ---------------------------------------------------------------------------- def simulate(rng: np.random.Generator, n_rep: int, k_max: int, p=None, beta=None, per_visit: bool = False, checkpoints=None, track_k=None): """ Simulate n_rep independent occupied sites, each visited k_max times. p : constant detection probability (Model H) beta : (a, b) Beta parameters. per_visit=False -> Model S (one p per site, reused every visit). per_visit=True -> Model V (a fresh p for every visit). Returns survival[k] = simulated P(no detection in first k visits), k=0..k_max, plus optional running estimates of survival[track_k] at the given checkpoints. """ hist = np.zeros(k_max + 2, dtype=np.int64) # first-detection visit; k_max+1 = none running = [] done = 0 miss_running = 0 cp = sorted(checkpoints) if checkpoints else [] cp_i = 0 while done < n_rep: m = min(CHUNK, n_rep - done) if beta is None: det = rng.random((m, k_max)) < p elif per_visit: pv = rng.beta(beta[0], beta[1], size=(m, k_max)) det = rng.random((m, k_max)) < pv else: ps = rng.beta(beta[0], beta[1], size=m) det = rng.random((m, k_max)) < ps[:, None] any_det = det.any(axis=1) first = np.where(any_det, det.argmax(axis=1) + 1, k_max + 1) hist += np.bincount(first, minlength=k_max + 2)[: k_max + 2] if track_k is not None: miss_running += int(np.count_nonzero(first > track_k)) done += m while cp_i < len(cp) and cp[cp_i] <= done: if cp[cp_i] == done: running.append((done, miss_running / done)) cp_i += 1 # survival[k] = P(first detection > k) = P(k blanks) tail = np.zeros(k_max + 1, dtype=np.int64) acc = 0 for k in range(k_max, -1, -1): acc += int(hist[k + 1]) tail[k] = acc survival = tail / n_rep return survival, running def mc_se(q: float, n: int) -> float: return math.sqrt(max(q * (1.0 - q), 0.0) / n) # ---------------------------------------------------------------------------- # MAIN # ---------------------------------------------------------------------------- def main() -> None: t0 = time.time() print(rule("=")) print("ABSENCE OF EVIDENCE: HOW MANY SURVEYS BEFORE A SPECIES IS 'GONE'") print("Science Journaling Club - occupancy detection model") print(rule("=")) print(f"master seed : {MASTER_SEED}") print(f"numpy version : {np.__version__}") print(f"python version : {sys.version.split()[0]}") print(f"replicates/config : {N_REP:,}") print(f"convergence run : {N_CONV:,}") print(f"target false-absence rate alpha : {ALPHA}") print("All random draws descend from the master seed via numpy SeedSequence.") print("No field data. Nothing here was observed; everything was computed.") # ======================================================================= head("PART 1 MODEL H (constant p). VISITS NEEDED FOR 95% CONFIDENCE") p_grid = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80] k_cap = 80 print() print("Closed form: P(k blanks | occupied) = (1-p)^k ; k* = ceil(ln 0.05 / ln(1-p))") print() print(f"{'p':>6} {'k* exact':>9} {'k* int':>7} {'(1-p)^k* exact':>15} " f"{'simulated':>11} {'SE':>9} {'sim-exact':>11} {'z':>7}") print(rule("-")) part1 = [] sims = {} for p in p_grid: rng = new_rng(f"H-p{p}") surv, _ = simulate(rng, N_REP, k_cap, p=p) sims[p] = surv ks_int = k_star_hom_int(p) ks_cont = k_star_hom_cont(p) exact = p_miss_hom(p, ks_int) sim = float(surv[ks_int]) se = mc_se(exact, N_REP) diff = sim - exact z = diff / se if se > 0 else float("nan") part1.append((p, ks_cont, ks_int, exact, sim, se, diff, z)) print(f"{p:>6.2f} {ks_cont:>9.2f} {ks_int:>7d} {exact:>15.9f} " f"{sim:>11.6f} {se:>9.6f} {diff:>+11.6f} {z:>+7.2f}") print() print("Reading: after k* visits with no detection, an occupied site still") print("produces that record with probability <= 0.05. The 'z' column is the") print("gap between simulation and closed form in Monte Carlo standard errors.") # ======================================================================= head("PART 2 VALIDATION GRID, MODEL H. SIMULATION vs (1-p)^k") print() print("Every cell with an expected blank count >= 25 is included in the") print("summary statistics. Cells rarer than that are printed but flagged,") print("because a Monte Carlo estimate of a probability near 1e-5 from 500,000") print("draws carries no useful resolution.") print() print(f"{'p':>6} {'k':>4} {'exact (1-p)^k':>15} {'simulated':>12} {'SE':>10} " f"{'sim-exact':>12} {'z':>7} flag") print(rule("-")) k_probe = [1, 2, 3, 5, 8, 13, 21, 34, 55] z_pool = [] anchor_h = [] for p in p_grid: surv = sims[p] last_usable = None for k in k_probe: if k > k_cap: continue exact = p_miss_hom(p, k) sim = float(surv[k]) se = mc_se(exact, N_REP) diff = sim - exact z = diff / se if se > 0 else float("nan") usable = exact * N_REP >= 25 flag = "" if usable else "thin" if usable: z_pool.append(z) last_usable = z print(f"{p:>6.2f} {k:>4d} {exact:>15.9f} {sim:>12.6f} {se:>10.6f} " f"{diff:>+12.6f} {z:>+7.2f} {flag}") if last_usable is not None: anchor_h.append(last_usable) zs = np.array(z_pool) ah = np.array(anchor_h) print(rule("-")) print(f"usable cells : {len(zs)}") print(f"mean z : {zs.mean():+.4f} (expected 0)") print(f"sd of z : {zs.std(ddof=1):.4f} (expected 1)") print(f"max |z| : {np.abs(zs).max():.4f}") print(f"cells with |z| > 3 : {int((np.abs(zs) > 3).sum())} of {len(zs)}") print(f"cells with |z| > 2 : {int((np.abs(zs) > 2).sum())} of {len(zs)} " f"(expect about {0.0455 * len(zs):.1f})") print() print("CAUTION ON THE MEAN. The cells in one row of this table all come from") print("the same 500,000 simulated sites, so their z values are strongly") print("positively correlated: a run of sites that happens to miss a little too") print("often pushes every k in that row the same way. Averaging all 76 cells") print("therefore does NOT have standard error 1/sqrt(76). The independent test") print("uses one anchor cell per p value (the largest k with enough resolution).") print(f"independent anchors : {len(ah)}") print(f"mean anchor z : {ah.mean():+.4f} " f"(SE {1/math.sqrt(len(ah)):.4f}, so {abs(ah.mean())*math.sqrt(len(ah)):.2f} SE from 0)") print(f"sd of anchor z : {ah.std(ddof=1):.4f}") print() if np.abs(zs).max() < 4.5 and abs(ah.mean()) * math.sqrt(len(ah)) < 3.0: print("VERDICT: simulation reproduces (1-p)^k within Monte Carlo error.") else: print("VERDICT: DISCREPANCY. Do not report; debug the code.") # ======================================================================= head("PART 3 CONVERGENCE OF THE MONTE CARLO ESTIMATE") conv_p, conv_k = 0.20, 13 conv_exact = p_miss_hom(conv_p, conv_k) checkpoints = [25_000, 50_000, 75_000, 100_000, 150_000, 200_000, 300_000, 400_000, 500_000, 700_000, 1_000_000, 1_250_000, 1_500_000, 1_750_000, 2_000_000] rng = new_rng("conv") _, running = simulate(rng, N_CONV, conv_k + 2, p=conv_p, checkpoints=checkpoints, track_k=conv_k) print() print(f"configuration: p = {conv_p}, k = {conv_k}") print(f"exact (1-p)^k = {conv_exact:.12f}") print() print(f"{'replicates':>12} {'running est.':>14} {'error':>12} {'2 SE':>11} " f"{'|err|/SE':>9} inside 2SE") print(rule("-")) conv_rows = [] for n, est in running: err = est - conv_exact se = mc_se(conv_exact, n) inside = "yes" if abs(err) <= 2 * se else "NO" conv_rows.append((n, est, err, se)) print(f"{n:>12,} {est:>14.8f} {err:>+12.8f} {2*se:>11.8f} " f"{abs(err)/se:>9.2f} {inside}") print() print("The error shrinks like 1/sqrt(N): a hundredfold increase in replicates") print("buys a tenfold reduction. This is the whole reason the closed form is") print("worth having.") # ======================================================================= head("PART 4 MODEL V. DETECTION REDRAWN FROM Beta EVERY VISIT") print() print("Claim to test: if p_j ~ Beta(a,b) independently at each visit, the") print("false-absence probability is exactly (1 - pbar)^k, i.e. heterogeneity") print("that reshuffles between visits costs nothing at all.") print() print(f"{'pbar':>6} {'kappa':>7} {'SD(p)':>8} {'k':>4} {'(1-pbar)^k':>13} " f"{'simulated':>12} {'diff':>12} {'z':>7}") print(rule("-")) zv = [] anchor_v = [] for pbar, kappa in [(0.20, 2.0), (0.20, 20.0), (0.40, 2.0), (0.40, 20.0)]: a, b = beta_ab(pbar, kappa) rng = new_rng(f"V-{pbar}-{kappa}") surv, _ = simulate(rng, N_REP, 30, beta=(a, b), per_visit=True) for k in [1, 3, 6, 10, 14]: exact = p_miss_hom(pbar, k) sim = float(surv[k]) se = mc_se(exact, N_REP) z = (sim - exact) / se zv.append(z) if k == 14: anchor_v.append(z) print(f"{pbar:>6.2f} {kappa:>7.1f} {beta_sd(pbar,kappa):>8.4f} {k:>4d} " f"{exact:>13.8f} {sim:>12.6f} {sim-exact:>+12.6f} {z:>+7.2f}") zva = np.array(zv) av = np.array(anchor_v) print(rule("-")) print(f"max |z| across Model V cells : {np.abs(zva).max():.3f}") print(f"mean z over all cells : {zva.mean():+.4f} (cells correlated within a row)") print(f"mean anchor z ({len(av)} independent): {av.mean():+.4f} " f"(SE {1/math.sqrt(len(av)):.3f}, {abs(av.mean())*math.sqrt(len(av)):.2f} SE from 0)") print("Confirmed. Per-visit variation in p is invisible in the false-absence") print("rate. Only the mean matters. That is a result about WHICH kind of") print("heterogeneity you should worry about.") # ======================================================================= head("PART 5 MODEL S. ONE p PER SITE, HELD FOR EVERY VISIT") print() print("Three independent routes to the same number:") print(" (a) simulation, 500,000 sites;") print(" (b) closed form E[(1-p)^k] = B(a, b+k) / B(a, b) via lgamma;") print(" (c) tanh-sinh (double exponential) quadrature of both the numerator") print(" and the normalising integral, using no gamma function at all.") print() print("Route (c) started life as a Gauss-Legendre scheme with the endpoint") print("singularity substituted away. It disagreed with route (b) by up to") print("8.6e-01 in the kappa = 100 rows, which is not a finding, it is a bug:") print("the substitution cures a singularity at one end and manufactures an") print("infinite derivative at the other whenever the exponent exceeds 1. The") print("comment block above p_miss_site_quad records the failure.") print() print("Quadrature self-test: log of the Beta integral against lgamma.") print(f"{'a':>8} {'b':>8} {'quadrature':>16} {'lgamma':>16} {'diff':>11}") print(rule("-")) for a_t, b_t in [(0.2, 1.8), (0.5, 4.5), (1.0, 9.0), (4.0, 36.0), (10.0, 90.0), (60.0, 40.0), (2.0, 3.0)]: q = log_beta_quad(a_t, b_t) g = math.lgamma(a_t) + math.lgamma(b_t) - math.lgamma(a_t + b_t) print(f"{a_t:>8.2f} {b_t:>8.2f} {q:>16.12f} {g:>16.12f} {abs(q-g):>11.2e}") print() s_configs = [(0.10, 2.0), (0.10, 5.0), (0.10, 20.0), (0.10, 100.0), (0.20, 2.0), (0.20, 5.0), (0.20, 20.0), (0.20, 100.0), (0.40, 2.0), (0.40, 5.0), (0.40, 20.0), (0.40, 100.0), (0.60, 2.0), (0.60, 5.0), (0.60, 20.0), (0.60, 100.0)] k_cap_s = 150 k_probe_s = [1, 3, 6, 10, 20, 40, 80, 150] print(f"{'pbar':>5} {'kap':>6} {'k':>4} {'closed form':>14} {'quadrature':>14} " f"{'|cf-quad|':>11} {'simulated':>11} {'SE':>9} {'z':>7} flag") print(rule("-")) zs_pool = [] anchor_s = [] max_cf_quad = 0.0 s_surv = {} for pbar, kappa in s_configs: a, b = beta_ab(pbar, kappa) rng = new_rng(f"S-{pbar}-{kappa}") surv, _ = simulate(rng, N_REP, k_cap_s, beta=(a, b), per_visit=False) s_surv[(pbar, kappa)] = surv last_usable_s = None for k in k_probe_s: cf = p_miss_site(a, b, k) qd = p_miss_site_quad(a, b, k) d = abs(cf - qd) max_cf_quad = max(max_cf_quad, d) sim = float(surv[k]) se = mc_se(cf, N_REP) z = (sim - cf) / se if se > 0 else float("nan") usable = cf * N_REP >= 25 if usable: zs_pool.append(z) last_usable_s = z print(f"{pbar:>5.2f} {kappa:>6.1f} {k:>4d} {cf:>14.9f} {qd:>14.9f} " f"{d:>11.2e} {sim:>11.6f} {se:>9.6f} {z:>+7.2f} " f"{'' if usable else 'thin'}") if last_usable_s is not None: anchor_s.append(last_usable_s) zsa = np.array(zs_pool) asr = np.array(anchor_s) print(rule("-")) print(f"max |closed form - quadrature| over all cells : {max_cf_quad:.3e}") print(f"usable simulation cells : {len(zsa)}") print(f"mean z : {zsa.mean():+.4f}") print(f"sd of z : {zsa.std(ddof=1):.4f}") print(f"max |z| : {np.abs(zsa).max():.4f}") print(f"cells with |z| > 3 : {int((np.abs(zsa)>3).sum())}") print(f"cells with |z| > 2 : {int((np.abs(zsa)>2).sum())} " f"(expect about {0.0455*len(zsa):.1f})") print() print("Same caution as Part 2: the eight k values in one configuration share") print("the same 500,000 simulated sites and are correlated, and here the") print("correlation is stronger still, because a site's p is fixed for life, so") print("one unlucky draw of site-level p values tilts the whole row. The mean of") print("all 110 cells is therefore not a 110-fold-precise test. The independent") print("statistic is one anchor cell per configuration.") print(f"independent anchors : {len(asr)}") print(f"mean anchor z : {asr.mean():+.4f} " f"(SE {1/math.sqrt(len(asr)):.4f}, {abs(asr.mean())*math.sqrt(len(asr)):.2f} SE from 0)") print(f"sd of anchor z : {asr.std(ddof=1):.4f}") print() if (max_cf_quad < 1e-11 and np.abs(zsa).max() < 4.5 and abs(asr.mean()) * math.sqrt(len(asr)) < 3.0): print("VERDICT: all three routes agree. Closed form and quadrature agree to") print("machine precision; the simulation agrees with both inside MC error.") else: print("VERDICT: DISCREPANCY. Do not report; debug the code.") # ======================================================================= head("PART 6 HOW MUCH DOES PERSISTENT HETEROGENEITY INFLATE THE EFFORT?") print() print("k* is the visits needed for a false-absence rate of 0.05. The naive") print("answer uses the mean detection probability and equation (2). The") print("honest answer inverts equation (4). The ratio is the inflation factor.") print() print(f"{'pbar':>5} {'kappa':>7} {'SD(p)':>8} {'a':>7} {'b':>7} " f"{'k* naive':>10} {'k* Model S':>14} {'inflation':>13}") print(rule("-")) infl_rows = [] for pbar in [0.10, 0.20, 0.40, 0.60]: kn = k_star_hom_cont(pbar) for kappa in [1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 1000.0]: a, b = beta_ab(pbar, kappa) ks = k_star_site_cont(a, b) ratio = ks / kn if math.isfinite(ks) else float("inf") infl_rows.append((pbar, kappa, beta_sd(pbar, kappa), a, b, kn, ks, ratio)) rs = f"{ratio:>13.2f}" if ratio < 1e6 else f"{ratio:>13.3e}" print(f"{pbar:>5.2f} {kappa:>7.1f} {beta_sd(pbar,kappa):>8.4f} " f"{a:>7.3f} {b:>7.3f} {kn:>10.2f} {fmt_k(ks):>14} {rs}") print() print("Why the numbers explode. For large k, B(a,b+k)/B(a,b) behaves like") print(" E[(1-p)^k] ~ [Gamma(a+b)/Gamma(b)] * k^(-a),") print("a POWER LAW in k, not a geometric decay. Once the Beta puts real weight") print("near p = 0, a fraction of sites are effectively undetectable and no") print("amount of repeat visiting drives the blank-record probability to zero at") print("the rate the naive formula promises.") print() print("Asymptotic check against the exact inversion:") print(f"{'pbar':>5} {'kappa':>7} {'a':>7} {'k* exact':>14} {'k* asymptotic':>15} " f"{'rel. err':>10}") print(rule("-")) for pbar in [0.10, 0.20, 0.40]: for kappa in [1.0, 2.0, 5.0]: a, b = beta_ab(pbar, kappa) ks = k_star_site_cont(a, b) c = math.lgamma(a + b) - math.lgamma(b) k_asym = math.exp((c - math.log(ALPHA)) / a) rel = (k_asym - ks) / ks if math.isfinite(ks) and ks > 0 else float("nan") print(f"{pbar:>5.2f} {kappa:>7.1f} {a:>7.3f} {fmt_k(ks):>14} " f"{fmt_k(k_asym):>15} {rel:>+10.4f}") # ======================================================================= head("PART 7 WHAT 20 BLANK VISITS ACTUALLY BUY YOU") print() print("Survey budgets are finite. Fix k = 20 blank visits and ask what the") print("false-absence probability really is, naive answer beside Model S.") print() print(f"{'pbar':>5} {'kappa':>7} {'naive (1-pbar)^20':>18} {'Model S closed':>16} " f"{'simulated':>11} {'SE':>9} {'ratio S/naive':>14}") print(rule("-")) for pbar, kappa in s_configs: a, b = beta_ab(pbar, kappa) naive = p_miss_hom(pbar, 20) cf = p_miss_site(a, b, 20) sim = float(s_surv[(pbar, kappa)][20]) se = mc_se(cf, N_REP) print(f"{pbar:>5.2f} {kappa:>7.1f} {naive:>18.9f} {cf:>16.9f} " f"{sim:>11.6f} {se:>9.6f} {cf/naive:>14.1f}") # ======================================================================= head("PART 8 CONVERGENCE, MODEL S") pbar_c, kappa_c, kc = 0.20, 5.0, 20 a_c, b_c = beta_ab(pbar_c, kappa_c) exact_c = p_miss_site(a_c, b_c, kc) rng = new_rng("convS") _, running_s = simulate(rng, N_CONV, kc + 2, beta=(a_c, b_c), per_visit=False, checkpoints=checkpoints, track_k=kc) print() print(f"configuration: pbar = {pbar_c}, kappa = {kappa_c} " f"(a = {a_c}, b = {b_c}), k = {kc}") print(f"closed form = {exact_c:.12f}") print(f"quadrature = {p_miss_site_quad(a_c, b_c, kc):.12f}") print() print(f"{'replicates':>12} {'running est.':>14} {'error':>12} {'2 SE':>11} " f"{'|err|/SE':>9} inside 2SE") print(rule("-")) conv_rows_s = [] for n, est in running_s: err = est - exact_c se = mc_se(exact_c, n) conv_rows_s.append((n, est, err, se)) print(f"{n:>12,} {est:>14.8f} {err:>+12.8f} {2*se:>11.8f} " f"{abs(err)/se:>9.2f} {'yes' if abs(err)<=2*se else 'NO'}") # ======================================================================= head("PART 9 POSTSCRIPT: FROM 'A BLANK RECORD IS UNLIKELY' TO 'IT IS GONE'") print() print("Everything above is conditional on the site being occupied. The") print("statement a conservation manager wants is the other way round. With a") print("prior probability psi that the site is occupied, Bayes gives") print() print(" P(occupied | k blanks) = psi (1-p)^k / [ psi (1-p)^k + (1 - psi) ]") print() print("Visits needed for P(occupied | k blanks) <= 0.05, by prior:") print() print(f"{'p':>6} {'psi=0.9':>9} {'psi=0.7':>9} {'psi=0.5':>9} {'psi=0.3':>9} " f"{'psi=0.1':>9} {'conditional k*':>15}") print(rule("-")) for p in [0.05, 0.10, 0.20, 0.30, 0.50, 0.80]: cells = [] for psi in [0.9, 0.7, 0.5, 0.3, 0.1]: thresh = 0.05 * (1 - psi) / (0.95 * psi) if thresh >= 1.0: cells.append("0") else: cells.append(str(int(math.ceil(math.log(thresh) / math.log(1 - p))))) print(f"{p:>6.2f} " + " ".join(f"{c:>9}" for c in cells) + f" {k_star_hom_int(p):>15}") print() print("A low prior does most of the work. If you already believed the species") print("was probably gone, a few blanks finish the argument; if you believed it") print("was probably there, no realistic number of blanks will move you. That is") print("not a flaw in the arithmetic, it is what the arithmetic is for.") # ======================================================================= head("PART 10 CHECK AGAINST PUBLISHED SURVEY NUMBERS") print() print("Kery (2002, J. Wildl. Manage. 66:330-338) surveyed three European snake") print("species over 645 visits to 87 sites and published both the per-visit") print("detection probabilities for small populations and the number of visits") print("needed to infer absence with 95 percent confidence. Those are exactly the") print("two quantities equation (2) links, so his table is a free external test of") print("our arithmetic. We did not fit anything to his data; we put his p into our") print("formula and compare the answer with the number he printed.") print() kery = [("Vipera aspis (asp viper)", 0.23, 12), ("Coronella austriaca (smooth snake)", 0.09, 34), ("Natrix natrix (grass snake)", 0.11, 26)] print(f"{'species':>36} {'p published':>12} {'our k*':>8} {'Kery k*':>9} " f"{'diff':>6} {'p band consistent with Kery k*':>31}") print(rule("-")) for name, pk, kk in kery: ours = k_star_hom_int(pk) p_lo = 1.0 - ALPHA ** (1.0 / kk) # any p above this needs <= kk visits p_hi = 1.0 - ALPHA ** (1.0 / (kk - 1)) # any p at or below this needs >= kk band = f"{p_lo:.4f} to {p_hi:.4f}" print(f"{name:>36} {pk:>12.2f} {ours:>8d} {kk:>9d} {ours-kk:>+6d} " f"{band:>31}") print() print("Two of the three land exactly on the published value. The smooth snake is") print("two visits short: our formula at p = 0.09 gives 32, Kery prints 34.") print("Inverting his 34 says the detection probability behind it was somewhere in") print("0.0843 to 0.0868. The midpoint of that band, 0.0855, rounds to the 0.09") print("printed in his table. So the two-visit gap is consistent with rounding in") print("the published detection probability rather than a disagreement about the") print("model, but we cannot prove that from the printed table alone, and we would") print("rather show the gap than quietly drop the one row that did not match.") # ======================================================================= head("PART 11 FIGURE DATA") print() print("[FIG1] false-absence curves (1-p)^k, closed form, k = 0..60") print("p," + ",".join(str(k) for k in range(0, 61, 5))) for p in [0.05, 0.10, 0.20, 0.40, 0.80]: print(f"{p}," + ",".join(f"{p_miss_hom(p,k):.6f}" for k in range(0, 61, 5))) print() print("[FIG2] k* against p (continuous and integer)") print("p,k_star_cont,k_star_int") for p in [0.05, 0.075, 0.10, 0.125, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 0.70, 0.80]: print(f"{p},{k_star_hom_cont(p):.4f},{k_star_hom_int(p)}") print() print("[FIG3] convergence, Model H, p=0.20 k=13") print(f"exact,{conv_exact:.10f}") print("n,estimate,error,se") for n, est, err, se in conv_rows: print(f"{n},{est:.8f},{err:+.8f},{se:.8f}") print("[FIG3b] convergence, Model S, pbar=0.20 kappa=5 k=20") print(f"exact,{exact_c:.10f}") print("n,estimate,error,se") for n, est, err, se in conv_rows_s: print(f"{n},{est:.8f},{err:+.8f},{se:.8f}") print() print("[FIG4] Model S survival curves, pbar = 0.20, k = 0..60") print("kappa,k,closed_form,simulated") for kappa in [2.0, 5.0, 20.0, 100.0]: a, b = beta_ab(0.20, kappa) surv = s_surv[(0.20, kappa)] for k in range(0, 61, 5): print(f"{kappa},{k},{p_miss_site(a,b,k):.6f},{surv[k]:.6f}") for k in range(0, 61, 5): print(f"inf,{k},{p_miss_hom(0.20,k):.6f},{sims[0.20][k]:.6f}") print() print("[FIG5] inflation factor k*(Model S) / k*(naive) against kappa") print("pbar,kappa,sd,k_naive,k_site,inflation") for pbar, kappa, sd, a, b, kn, ks, ratio in infl_rows: print(f"{pbar},{kappa},{sd:.5f},{kn:.4f}," f"{ks if math.isfinite(ks) else 'inf'}," f"{ratio if math.isfinite(ratio) else 'inf'}") print() print("[FIG6] Bayesian posterior P(occupied | k blanks), psi = 0.5") print("p," + ",".join(str(k) for k in range(0, 41, 2))) for p in [0.05, 0.10, 0.20, 0.40, 0.80]: row = [] for k in range(0, 41, 2): s = (1 - p) ** k row.append(f"{0.5*s/(0.5*s+0.5):.6f}") print(f"{p}," + ",".join(row)) head("SUMMARY") print() for p, ks_cont, ks_int, exact, sim, se, diff, z in part1: print(f" p = {p:<5} -> {ks_int:>3d} blank visits for 95% confidence " f"(exact k* = {ks_cont:6.2f}); simulated miss rate {sim:.5f} " f"vs {exact:.5f} exact, z = {z:+.2f}") print() print(f"Model H validation : max |z| = {np.abs(zs).max():.3f} over {len(zs)} cells; " f"mean anchor z = {ah.mean():+.4f} over {len(ah)} independent anchors " f"({abs(ah.mean())*math.sqrt(len(ah)):.2f} SE from 0)") print(f"Model V validation : max |z| = {np.abs(zva).max():.3f}; mean anchor z = " f"{av.mean():+.4f} ({abs(av.mean())*math.sqrt(len(av)):.2f} SE from 0). " f"Per-visit heterogeneity costs nothing, as predicted.") print(f"Model S validation : closed form vs quadrature max diff " f"{max_cf_quad:.2e}; simulation max |z| = {np.abs(zsa).max():.3f}; " f"mean anchor z = {asr.mean():+.4f} " f"({abs(asr.mean())*math.sqrt(len(asr)):.2f} SE from 0)") print() ks20_5 = k_star_site_cont(*beta_ab(0.20, 5.0)) print("Headline: at p = 0.20 a clean run of " f"{k_star_hom_int(0.20)} blanks is 95% evidence under Model H. Let") print(f"detection vary between sites with kappa = 5 (SD " f"{beta_sd(0.20,5.0):.3f}) and the same standard") print(f"needs {fmt_k(ks20_5)} visits, an inflation of " f"{ks20_5/k_star_hom_cont(0.20):.1f}x.") print() print(f"total runtime: {time.time()-t0:.1f} s") print(rule("=")) if __name__ == "__main__": main()