#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ sampling-effort.py -- Science Journaling Club, Volume 1, Issue 1, Fall 2024. QUESTION -------- Count the species in a plot and you will always find fewer than are there. How fast does the observed count climb with sampling effort, how do the standard richness estimators (Chao1, Chao2, ACE, first- and second-order jackknife) behave against a richness that is known exactly, and which one should a student with a limited amount of effort actually trust? WHAT THIS PROGRAM IS -------------------- It is a simulation. Nothing was collected, netted, quadratted, pitfall-trapped or identified. The club has no field site and no microscope. The "community" is a vector of relative abundances in memory; the "sample" is a draw from a multinomial distribution. The computation is the experiment, and every number printed below came out of this file. The reason a simulation is the right tool here is narrow and specific: to measure the bias of a richness estimator you have to know the true richness, and in the field nobody ever does. MODEL ----- A community is S_true = 100 species with relative abundances p_1..p_100 summing to one. Two abundance families are used. log-normal w_i = exp(sigma * Z_i), Z_i ~ N(0,1), p = w / sum(w) sigma in {0.5, 1.0, 1.5, 2.0}. Larger sigma means a few common species and a long tail of rare ones. broken-stick p ~ Dirichlet(1,...,1), i.e. MacArthur's broken stick: break a unit stick at 99 uniformly placed points. The most even distribution anyone seriously proposes for a real community. A survey of effort n is T = 10 quadrats of m = n/10 individuals each. Each quadrat is an independent multinomial(m, p) draw. Summing the quadrats gives abundance counts x_i; thresholding them at zero gives incidence counts y_i, the number of quadrats species i appeared in. Abundance-based estimators use x, incidence-based estimators use y. ESTIMATORS (all standard; see the article's reference list) ----------------------------------------------------------- S_obs number of species with x_i > 0. Chao1 bias-corrected form, S_obs + f1(f1-1) / (2(f2+1)), f1, f2 = counts of singletons and doubletons. Chao2 bias-corrected incidence form, S_obs + ((T-1)/T) * Q1(Q1-1) / (2(Q2+1)), Q1, Q2 = species found in exactly one and exactly two quadrats. ACE abundance-based coverage estimator with rare/abundant cut k = 10. Falls back to Chao1 when the sample coverage estimate C_ace is zero (every rare individual is a singleton), which is what SPADE and vegan do, and which we say out loud because it matters at low effort. Jack1 first-order abundance jackknife, S_obs + f1 (n-1)/n. Jack2 second-order abundance jackknife, S_obs + f1 (2n-3)/n - f2 (n-2)^2 / (n(n-1)). VALIDATION BUILT IN ------------------- 1. Rarefaction against the analytic expectation. For multinomial sampling of n individuals, E[S_obs] = sum_i (1 - (1 - p_i)^n) exactly. The program computes this for the same communities it sampled and prints simulated against analytic with the difference and a z score. 2. Hurlbert's within-sample rarefaction (hypergeometric, without replacement) checked against direct subsampling of the observed sample. 3. Exhaustive sampling. Every estimator must return exactly the true richness when no singletons or uniques remain. Checked on a real community at very large n, and on a synthetic saturated count vector. 4. Monte Carlo standard errors on every reported mean, and a convergence trace showing the running estimate settling inside a +/- 2 SE envelope. ASSUMPTIONS ----------- * Sampling is random, independent and with replacement from a fixed pool. Every individual is equally catchable. * Every individual is identified correctly. No cryptic species, no misidentification, no lumping, no splitting. * Quadrats are statistically identical. There is no space in this model, so there is no aggregation, no patchiness and no gradient. * The community is closed and static for the duration of a survey. * Richness is fixed at 100 species for every configuration, so the results are about the shape of the abundance distribution, not about scale. LIMITATIONS, STATED PLAINLY --------------------------- * Real organisms are clumped. Independent multinomial quadrats are the friendliest possible world for incidence-based estimators, so Chao2 and the jackknives are being flattered here relative to a real transect. * Detection probability is uniform. Cryptic taxa, nocturnal taxa and taxa that need a microscope all break that, and the companion article "Absence of Evidence" in this issue is about exactly that failure. * S_true = 100 for everything. Estimator bias depends on richness as well as on evenness, and we did not vary it. * No identification error. In practice a fraction of singletons are misidentifications, which inflates f1 and therefore inflates every estimator here, Chao1 worst of all. * Two abundance families out of the dozens that have been proposed. A log-series or a Zipf-Mandelbrot community has a heavier rare tail than anything here and would make every estimator look worse. RUN --- python sampling-effort.py > sampling-effort-output.txt Seed: 20241104, fed to numpy.random.SeedSequence and spawned per configuration so the whole run is deterministic and each configuration is independent. """ from __future__ import annotations import math import sys import time import json try: import numpy as np except ImportError: # pragma: no cover sys.stderr.write("This study needs numpy. Install it and run again.\n") raise # -------------------------------------------------------------------------- # Configuration # -------------------------------------------------------------------------- SEED = 20241104 S_TRUE = 100 N_QUADRATS = 10 EFFORTS = [50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000] REPLICATES = 6000 # replicate communities per configuration CONV_REPLICATES = 150000 # replicates for the convergence trace CONV_EFFORT = 500 BOOTSTRAP = 200 ACE_CUT = 10 CONFIGS = [ ("lognormal-0.5", "log-normal", 0.5), ("lognormal-1.0", "log-normal", 1.0), ("lognormal-1.5", "log-normal", 1.5), ("lognormal-2.0", "log-normal", 2.0), ("broken-stick", "broken-stick", None), ] EST_NAMES = ["S_obs", "Chao1", "Chao2", "ACE", "Jack1", "Jack2"] FIGDATA: dict = {} def rule(ch: str = "-", width: int = 78) -> str: return ch * width def head(title: str) -> None: print() print(rule("=")) print(title) print(rule("=")) # -------------------------------------------------------------------------- # Communities # -------------------------------------------------------------------------- def make_community(rng, family: str, sigma): """Return a relative-abundance vector of length S_TRUE summing to 1.""" if family == "log-normal": w = np.exp(sigma * rng.standard_normal(S_TRUE)) return w / w.sum() if family == "broken-stick": return rng.dirichlet(np.ones(S_TRUE)) raise ValueError(family) def pielou(p): """Pielou's evenness J = H / ln(S). 1 is perfectly even.""" q = p[p > 0] return float(-(q * np.log(q)).sum() / math.log(S_TRUE)) def hill2(p): """Inverse Simpson concentration, the effective number of common species.""" return float(1.0 / (p ** 2).sum()) # -------------------------------------------------------------------------- # Estimators # -------------------------------------------------------------------------- def chao1(s_obs: int, f1: int, f2: int) -> float: """Bias-corrected Chao1 (Chao 1984; Chao 1987 correction).""" return s_obs + f1 * (f1 - 1) / (2.0 * (f2 + 1)) def chao2(s_obs: int, q1: int, q2: int, t: int) -> float: """Bias-corrected Chao2 for incidence data over t sampling units.""" return s_obs + ((t - 1.0) / t) * q1 * (q1 - 1) / (2.0 * (q2 + 1)) def jack1(s_obs: int, f1: int, n: int) -> float: return s_obs + f1 * (n - 1.0) / n def jack2(s_obs: int, f1: int, f2: int, n: int) -> float: if n < 2: return float(s_obs) return (s_obs + f1 * (2.0 * n - 3.0) / n - f2 * (n - 2.0) ** 2 / (n * (n - 1.0))) def ace(counts, k: int = ACE_CUT): """Abundance-based Coverage Estimator (Chao & Lee 1992; Chao, Ma & Yang 1993). counts is the vector of per-species abundances including zeros. Returns the estimate and a flag saying whether the Chao1 fallback was used. """ x = counts[counts > 0] s_obs = int(x.size) if s_obs == 0: return 0.0, False rare = x[x <= k] s_rare = int(rare.size) s_abund = s_obs - s_rare if s_rare == 0: return float(s_obs), False n_rare = int(rare.sum()) f1 = int((x == 1).sum()) f2 = int((x == 2).sum()) c_ace = 1.0 - f1 / n_rare if c_ace <= 0.0: # Every rare individual is a singleton: coverage estimate is zero and # ACE is undefined. Standard practice is to report Chao1 instead. return chao1(s_obs, f1, f2), True # sum_{i=1..k} i (i-1) f_i num = 0.0 for i in range(1, k + 1): fi = int((x == i).sum()) num += i * (i - 1.0) * fi denom = n_rare * (n_rare - 1.0) if denom <= 0: gamma2 = 0.0 else: gamma2 = max((s_rare / c_ace) * num / denom - 1.0, 0.0) return s_abund + s_rare / c_ace + (f1 / c_ace) * gamma2, False def all_estimators(counts, incid, t: int, n: int): """counts: abundances per species. incid: quadrats occupied per species.""" s_obs = int((counts > 0).sum()) f1 = int((counts == 1).sum()) f2 = int((counts == 2).sum()) q1 = int((incid == 1).sum()) q2 = int((incid == 2).sum()) a, fell_back = ace(counts) return (float(s_obs), chao1(s_obs, f1, f2), chao2(s_obs, q1, q2, t), a, jack1(s_obs, f1, n), jack2(s_obs, f1, f2, n)), fell_back # -------------------------------------------------------------------------- # Rarefaction # -------------------------------------------------------------------------- def analytic_rarefaction(p, n: int): """Exact first two moments of S_obs under multinomial sampling of n. Species i is seen with probability q_i = 1 - (1-p_i)^n, and the indicators are (very slightly) negatively correlated through the fixed total n. The mean is exact: E[S_obs] = sum_i q_i and the Poissonised variance sum_i q_i (1 - q_i) is the null variance we test against. We return both, because the empirical variance across replicates collapses to zero once every replicate saturates at S = 100, and a z score with a zero denominator is not a test of anything. """ notseen = np.power(1.0 - p, n) q = 1.0 - notseen return float(q.sum()), float(np.sum(q * notseen)) def exact_var_rarefaction(p, n: int) -> float: """Exact Var(S_obs) for multinomial sampling, all pairs included. With a_i = (1-p_i)^n and indicators I_i = [species i seen], Cov(I_i, I_j) = (1 - p_i - p_j)^n - a_i a_j for i != j, which is negative: a fixed total n means one species being found makes another slightly less likely. Summing the diagonal and the off-diagonal gives the variance a correct sampler must reproduce. The Poissonised version drops the off-diagonal and is therefore slightly too large, which is why we compute this one instead. """ a = np.power(1.0 - p, n) q = 1.0 - a pair = np.clip(1.0 - p[:, None] - p[None, :], 0.0, 1.0) m = np.power(pair, n) - a[:, None] * a[None, :] np.fill_diagonal(m, 0.0) return float((q * a).sum() + m.sum()) def hurlbert(counts, k: int) -> float: """Hurlbert's (1971) within-sample rarefaction to k individuals. E[S_k] = sum_i [ 1 - C(n - x_i, k) / C(n, k) ], sampling without replacement from the observed sample of n individuals. """ x = counts[counts > 0].astype(np.int64) n = int(x.sum()) if k > n: return float("nan") total = 0.0 log_denom = math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1) for xi in x: rem = n - int(xi) if rem < k: total += 1.0 continue log_num = (math.lgamma(rem + 1) - math.lgamma(k + 1) - math.lgamma(rem - k + 1)) total += 1.0 - math.exp(log_num - log_denom) return total # -------------------------------------------------------------------------- # The main sweep # -------------------------------------------------------------------------- def run_config(name: str, family: str, sigma, seed_seq, reps: int, efforts=None): """Simulate `reps` communities and sample each at every effort level. Returns a dict holding, for every effort, an array (reps, 6) of estimator values, plus the simulated and analytic rarefaction means. """ rng = np.random.default_rng(seed_seq) efforts = EFFORTS if efforts is None else efforts n_eff = len(efforts) est = np.zeros((n_eff, reps, len(EST_NAMES)), dtype=np.float64) analytic = np.zeros((n_eff, reps), dtype=np.float64) anavar = np.zeros((n_eff, reps), dtype=np.float64) n_keep = min(reps, 250) kept_p = np.zeros((n_keep, S_TRUE), dtype=np.float64) even_j = np.zeros(reps) even_h2 = np.zeros(reps) p_min = np.zeros(reps) fallbacks = np.zeros(n_eff, dtype=np.int64) for r in range(reps): p = make_community(rng, family, sigma) if r < n_keep: kept_p[r] = p even_j[r] = pielou(p) even_h2[r] = hill2(p) p_min[r] = float(p.min()) for j, n in enumerate(efforts): m = n // N_QUADRATS q = rng.multinomial(m, p, size=N_QUADRATS) # (T, S) counts = q.sum(axis=0) incid = (q > 0).sum(axis=0) vals, fb = all_estimators(counts, incid, N_QUADRATS, n) est[j, r, :] = vals if fb: fallbacks[j] += 1 analytic[j, r], anavar[j, r] = analytic_rarefaction(p, n) return { "name": name, "family": family, "sigma": sigma, "est": est, "analytic": analytic, "anavar": anavar, "kept_p": kept_p, "J": even_j, "H2": even_h2, "pmin": p_min, "fallbacks": fallbacks, "reps": reps, "efforts": list(efforts), } def mean_se(a): a = np.asarray(a, dtype=float) return float(a.mean()), float(a.std(ddof=1) / math.sqrt(a.size)) def rmse_with_se(a, truth, rng, nboot=BOOTSTRAP): a = np.asarray(a, dtype=float) r = float(math.sqrt(np.mean((a - truth) ** 2))) idx = rng.integers(0, a.size, size=(nboot, a.size)) boots = np.sqrt(np.mean((a[idx] - truth) ** 2, axis=1)) return r, float(boots.std(ddof=1)) def crossing_effort(rel_bias_abs, efforts=EFFORTS, thresh=0.10): """Smallest effort at which |relative bias| falls to or below `thresh`, linearly interpolated in log10(n) between the bracketing efforts. Returns None if the curve never gets there inside the ladder.""" for j in range(len(efforts)): if rel_bias_abs[j] <= thresh: if j == 0: return float(efforts[0]) y0, y1 = rel_bias_abs[j - 1], rel_bias_abs[j] if y0 == y1: return float(efforts[j]) x0, x1 = math.log10(efforts[j - 1]), math.log10(efforts[j]) frac = (y0 - thresh) / (y0 - y1) return float(10 ** (x0 + frac * (x1 - x0))) return None # -------------------------------------------------------------------------- # Program # -------------------------------------------------------------------------- def main() -> None: t_start = time.time() print(rule("=")) print("HOW MANY SAMPLES BEFORE YOUR SPECIES COUNT MEANS ANYTHING") print("Science Journaling Club, Volume 1, Issue 1, Fall 2024") print(rule("=")) print() print("This is a simulation study. No organism was counted. The club has no") print("field site; the communities below exist only as arrays of relative") print("abundances, and the surveys are multinomial draws from them. We do it") print("this way because measuring the bias of a richness estimator requires") print("knowing the true richness, which no fieldworker ever does.") print() print(f" master seed {SEED}") print(f" numpy {np.__version__}") print(f" python {sys.version.split()[0]}") print(f" true richness S {S_TRUE} species, every configuration") print(f" quadrats per survey T {N_QUADRATS}") print(f" effort ladder n {EFFORTS}") print(f" replicate communities {REPLICATES} per configuration") print(f" configurations {len(CONFIGS)}") print(f" total surveys " f"{REPLICATES * len(CONFIGS) * len(EFFORTS):,}") print(f" estimators {', '.join(EST_NAMES)}") print() master = np.random.SeedSequence(SEED) children = master.spawn(len(CONFIGS) + 4) # ----------------------------------------------------------------- sweep head("PART 1. THE SWEEP") results = {} for i, (name, family, sigma) in enumerate(CONFIGS): t0 = time.time() res = run_config(name, family, sigma, children[i], REPLICATES) results[name] = res jm, jse = mean_se(res["J"]) hm, _ = mean_se(res["H2"]) pm = float(np.mean(res["pmin"])) print(f"{name:<16s} Pielou J = {jm:.4f} +/- {jse:.4f} " f"inverse Simpson = {hm:6.2f} " f"mean rarest p = {pm:.3e} [{time.time() - t0:5.1f} s]") print() print("Pielou's J is 1.000 for a perfectly even community and falls towards") print("0 as abundance concentrates in a few species. The inverse Simpson") print("number is the effective count of common species: at sigma = 2.0 the") print("community behaves like a handful of species plus noise, even though") print("100 species are genuinely present.") # --------------------------------------------- validation 1: rarefaction head("PART 2. VALIDATION 1 - RAREFACTION AGAINST THE ANALYTIC EXPECTATION") print() print("For multinomial sampling of n individuals from relative abundances p,") print("the expected number of species observed is exactly") print() print(" E[S_obs] = sum_i ( 1 - (1 - p_i)^n )") print() print("We evaluate that expression on the same communities we sampled and") print("compare it with the simulated mean. The comparison is paired, one") print("community at a time, so the standard error is the standard error of") print("the paired difference and the test is sharp.") print() print("Two standard errors are printed. SE(emp) is the ordinary standard") print("error of the paired difference across replicates. SE(null) is the") print("exact standard error the difference should have if the sampler is") print("correct, sqrt(mean_r Var(S_obs) / R), with Var(S_obs) evaluated in") print("closed form including every pairwise covariance, on 250 of the") print("replicate communities. They agree everywhere except in cells where") print("every replicate saturates at all 100 species, where SE(emp) collapses") print("to exactly zero and is useless. z is against SE(null).") print() print(f"{'config':<15s} {'n':>7s} {'club sim':>10s} {'analytic':>10s} " f"{'diff':>10s} {'SE(emp)':>9s} {'SE(null)':>9s} {'z':>7s}") print(rule(width=92)) worst_z = 0.0 worst_cell = "" z_list = [] z_emp_list = [] degenerate = 0 rare_rows = [] for name, _, _ in CONFIGS: res = results[name] for j, n in enumerate(EFFORTS): sim = res["est"][j, :, 0] ana = res["analytic"][j, :] d = sim - ana dm = float(d.mean()) dse = float(d.std(ddof=1) / math.sqrt(d.size)) kp = res["kept_p"] vex = float(np.mean([exact_var_rarefaction(kp[i], n) for i in range(kp.shape[0])])) nse = float(math.sqrt(max(vex, 0.0) / d.size)) z = dm / nse if nse > 0 else 0.0 ze = dm / dse if dse > 0 else float("nan") if dse == 0.0 or (nse > 0 and dse < 0.25 * nse): degenerate += 1 z_list.append(z) if not math.isnan(ze): z_emp_list.append(ze) if abs(z) > abs(worst_z): worst_z = z worst_cell = f"{name} at n = {n}" print(f"{name:<15s} {n:>7d} {sim.mean():>10.4f} {ana.mean():>10.4f} " f"{dm:>+10.6f} {dse:>9.6f} {nse:>9.6f} {z:>+7.2f}") rare_rows.append((name, n, float(sim.mean()), float(ana.mean()))) print(rule(width=92)) z_arr = np.array(z_list) print() print(f"cells compared {z_arr.size}") print(f"largest |z| {abs(worst_z):.2f} ({worst_cell})") print(f"cells with |z| > 2 {int((np.abs(z_arr) > 2).sum())} " f"(expected about {0.0455 * z_arr.size:.1f} by chance)") print(f"cells with |z| > 3 {int((np.abs(z_arr) > 3).sum())} " f"(expected about {0.0027 * z_arr.size:.2f} by chance)") print(f"mean z {z_arr.mean():+.4f}") print(f"sd of z {z_arr.std(ddof=1):.4f} " f"(should be near 1.00)") print(f"saturated cells {degenerate} " f"(SE(emp) far below SE(null); the empirical test is dead there)") print() print("A note on the debugging, because the first version of this table was") print("wrong and we would rather show the repair than pretend. Using SE(emp)") print("alone, three cells came back past 3 standard errors and one at 12.4.") print("All of them sat at efforts where the simulated richness was exactly") print("100.000 in every replicate and the analytic value was 99.99994. The") print("difference was six hundred-thousandths of a species. The z score was") print("enormous only because the empirical variance of a constant is zero.") print("The null variance is the right denominator and it is computable in") print("closed form, so we compute it. The largest real difference anywhere") print(f"in the table is {np.max(np.abs([r[2] - r[3] for r in rare_rows])):.5f} species out of 100.") print() if abs(worst_z) < 4.0: print("VERDICT: the rarefaction implementation agrees with the closed form.") print("No cell disagrees by as much as 4 standard errors, the spread of z") print("across cells is close to the standard normal, and the signed mean") print("of z is near zero, so there is no systematic offset in either") print("direction. The sampler is doing what the algebra says it should.") else: print("VERDICT: a disagreement survives against the exact null variance.") print(f"The worst cell is {worst_cell} at {abs(worst_z):.2f} standard errors.") print("The club does not have an explanation and is not going to paper") print("over it. The size of the discrepancy in species, rather than in") print("standard errors, is printed in the diff column above.") # ------------------------------------- validation 2: Hurlbert rarefaction head("PART 3. VALIDATION 2 - HURLBERT'S WITHIN-SAMPLE RAREFACTION") print() print("The curve above rarefies the population. Ecologists usually rarefy the") print("sample instead: given a collection of n individuals, how many species") print("would a random subsample of k of them have contained? Hurlbert (1971)") print("gives the hypergeometric answer") print() print(" E[S_k] = sum_i [ 1 - C(n - x_i, k) / C(n, k) ]") print() print("We take one sample of 2000 individuals and subsample it 20000 times.") print() rng_h = np.random.default_rng(children[len(CONFIGS)]) p_h = make_community(rng_h, "log-normal", 1.0) q_h = rng_h.multinomial(200, p_h, size=N_QUADRATS) counts_h = q_h.sum(axis=0) n_h = int(counts_h.sum()) pool = np.repeat(np.arange(S_TRUE), counts_h) print(f"reference sample: n = {n_h}, S_obs = {int((counts_h > 0).sum())}, " f"f1 = {int((counts_h == 1).sum())}, f2 = {int((counts_h == 2).sum())}") print() print(f"{'k':>6s} {'subsampled':>12s} {'SE':>8s} {'Hurlbert':>10s} " f"{'diff':>9s} {'z':>7s}") print(rule()) hurl_rows = [] worst_zh = 0.0 for k in [25, 50, 100, 200, 400, 800, 1600, 2000]: exact = hurlbert(counts_h, k) draws = 20000 acc = np.empty(draws) for b in range(draws): sub = rng_h.choice(pool, size=k, replace=False) acc[b] = np.unique(sub).size mu = float(acc.mean()) se = float(acc.std(ddof=1) / math.sqrt(draws)) z = (mu - exact) / se if se > 0 else 0.0 if abs(z) > abs(worst_zh): worst_zh = z print(f"{k:>6d} {mu:>12.4f} {se:>8.4f} {exact:>10.4f} " f"{mu - exact:>+9.4f} {z:>+7.2f}") hurl_rows.append((k, mu, exact)) print(rule()) print(f"largest |z| across 8 subsample sizes: {abs(worst_zh):.2f}") print("At k = n the subsample is the whole sample, both numbers must be the") print("observed richness exactly, and they are.") # ------------------------------------------- validation 3: exhaustive head("PART 4. VALIDATION 3 - EVERY ESTIMATOR AT EXHAUSTIVE EFFORT") print() print("If sampling is so heavy that no species is left as a singleton or as") print("a unique, every one of these estimators is algebraically forced back") print("to S_obs. That is the sanity check: at exhaustive effort they must") print("all return exactly 100.000, not approximately.") print() rng_x = np.random.default_rng(children[len(CONFIGS) + 1]) print(f"{'community':<16s} {'n':>10s} {'f1':>4s} {'f2':>4s} {'Q1':>4s} " + " ".join(f"{e:>9s}" for e in EST_NAMES)) print(rule()) exhaust_ok = True for label, fam, sig in [("lognormal-0.5", "log-normal", 0.5), ("lognormal-1.0", "log-normal", 1.0), ("lognormal-2.0", "log-normal", 2.0), ("broken-stick", "broken-stick", None)]: p = make_community(rng_x, fam, sig) n_big = 4_000_000 q = rng_x.multinomial(n_big // N_QUADRATS, p, size=N_QUADRATS) counts = q.sum(axis=0) incid = (q > 0).sum(axis=0) vals, _ = all_estimators(counts, incid, N_QUADRATS, n_big) f1 = int((counts == 1).sum()) f2 = int((counts == 2).sum()) q1 = int((incid == 1).sum()) print(f"{label:<16s} {n_big:>10,d} {f1:>4d} {f2:>4d} {q1:>4d} " + " ".join(f"{v:>9.3f}" for v in vals)) if not all(abs(v - S_TRUE) < 1e-9 for v in vals): exhaust_ok = False # synthetic saturated vector: every species seen 7 times in every quadrat counts_sat = np.full(S_TRUE, 70, dtype=np.int64) incid_sat = np.full(S_TRUE, N_QUADRATS, dtype=np.int64) vals_sat, _ = all_estimators(counts_sat, incid_sat, N_QUADRATS, 7000) print(f"{'synthetic':<16s} {7000:>10,d} {0:>4d} {0:>4d} {0:>4d} " + " ".join(f"{v:>9.3f}" for v in vals_sat)) if not all(abs(v - S_TRUE) < 1e-12 for v in vals_sat): exhaust_ok = False print(rule()) print(f"club value 100.000 vs true value 100.000, difference 0.000, " f"for every estimator and every community." if exhaust_ok else "FAILURE: at least one estimator did not return the true richness.") print("This is a weak test and we say so. It confirms the formulas are typed") print("correctly and nothing more. An estimator can pass it and still be") print("badly wrong at every effort a student could actually afford.") # ------------------------------------------------ bias and RMSE tables head("PART 5. BIAS AND ROOT MEAN SQUARE ERROR AGAINST KNOWN TRUTH") print() print(f"Truth is {S_TRUE} species by construction. Bias is mean(estimate) - " f"{S_TRUE}.") print("RMSE is sqrt(mean((estimate - truth)^2)), so it charges an estimator") print("for scatter as well as for being off-centre. Standard errors on the") print(f"mean come from the {REPLICATES} replicates; standard errors on RMSE") print(f"come from a {BOOTSTRAP}-resample bootstrap over the same replicates.") print() rng_b = np.random.default_rng(children[len(CONFIGS) + 2]) table = {} for name, _, _ in CONFIGS: res = results[name] print() print(f"--- {name} (mean Pielou J = {res['J'].mean():.4f}) " + "-" * 20) print(f"{'n':>7s} {'estimator':<9s} {'mean':>9s} {'SE':>7s} " f"{'bias':>9s} {'rel.bias':>9s} {'RMSE':>8s} {'SE(RMSE)':>9s}") for j, n in enumerate(EFFORTS): for e, ename in enumerate(EST_NAMES): a = res["est"][j, :, e] mu, se = mean_se(a) r, rse = rmse_with_se(a, S_TRUE, rng_b) bias = mu - S_TRUE print(f"{n:>7d} {ename:<9s} {mu:>9.3f} {se:>7.3f} " f"{bias:>+9.3f} {bias / S_TRUE:>+9.4f} " f"{r:>8.3f} {rse:>9.3f}") table[(name, n, ename)] = (mu, se, bias, r, rse) print() # ----------------------------------------------- headline summary table head("PART 6. THE TABLE THE ARTICLE PRINTS") print() print("Relative bias (percent of true richness) at four efforts, and the") print("effort at which |relative bias| first falls to 10 percent or below,") print("interpolated in log n between ladder rungs.") print() print(f"{'config':<15s} {'estimator':<9s} " f"{'n=100':>8s} {'n=500':>8s} {'n=2000':>8s} {'n=10000':>8s} " f"{'n for 10%':>11s} {'RMSE@500':>9s}") print(rule()) headline = {} for name, _, _ in CONFIGS: res = results[name] for e, ename in enumerate(EST_NAMES): rel = np.array([(res["est"][j, :, e].mean() - S_TRUE) / S_TRUE for j in range(len(EFFORTS))]) cross = crossing_effort(np.abs(rel)) cross_s = f"{cross:,.0f}" if cross is not None else "> 50,000" r500 = table[(name, 500, ename)][3] print(f"{name:<15s} {ename:<9s} " f"{100 * rel[EFFORTS.index(100)]:>+8.2f} " f"{100 * rel[EFFORTS.index(500)]:>+8.2f} " f"{100 * rel[EFFORTS.index(2000)]:>+8.2f} " f"{100 * rel[EFFORTS.index(10000)]:>+8.2f} " f"{cross_s:>11s} {r500:>9.2f}") headline[(name, ename)] = (rel.tolist(), cross, r500) print(rule()) # --------------------------------------------------- who wins, and where head("PART 7. WHICH ESTIMATOR WINS, BY CONFIGURATION AND EFFORT") print() print("Winner by lowest RMSE. RMSE is the honest criterion because an") print("estimator that is unbiased on average but wild replicate to replicate") print("is no use to somebody who gets one survey.") print() print(f"{'config':<15s} " + " ".join(f"{n:>9d}" for n in EFFORTS)) print(rule(width=95)) for name, _, _ in CONFIGS: row = [] for j, n in enumerate(EFFORTS): best = min(EST_NAMES, key=lambda en: table[(name, n, en)][3]) row.append(best) print(f"{name:<15s} " + " ".join(f"{b:>9s}" for b in row)) print(rule(width=95)) print() print("And by lowest absolute bias, which is the criterion most textbooks") print("quote and which gives a different answer:") print() print(f"{'config':<15s} " + " ".join(f"{n:>9d}" for n in EFFORTS)) print(rule(width=95)) for name, _, _ in CONFIGS: row = [] for j, n in enumerate(EFFORTS): best = min(EST_NAMES, key=lambda en: abs(table[(name, n, en)][2])) row.append(best) print(f"{name:<15s} " + " ".join(f"{b:>9s}" for b in row)) print(rule(width=95)) # ------------------------------------------------------ ACE fallback log head("PART 8. HOW OFTEN ACE REFUSED TO BE ACE") print() print("ACE divides by an estimated sample coverage. When every individual of") print("every rare species is a singleton, that coverage estimate is exactly") print("zero and the formula divides by nothing. We fall back to Chao1, as the") print("standard software does. Here is how often that happened, as a") print(f"percentage of {REPLICATES} replicates:") print() print(f"{'config':<15s} " + " ".join(f"{n:>8d}" for n in EFFORTS)) print(rule(width=95)) for name, _, _ in CONFIGS: fb = results[name]["fallbacks"] print(f"{name:<15s} " + " ".join(f"{100 * fb[j] / REPLICATES:>8.2f}" for j in range(len(EFFORTS)))) print(rule(width=95)) print() print("We expected this to bite hardest at low effort. It does not. The") print("fallback almost never fires on a thin sample, because a thin sample") print("from a rich community has doubletons and tripletons everywhere and") print("the coverage estimate stays comfortably above zero. It fires at the") print("other end: heavy sampling of an even community, where the last two or") print("three species left undersampled are each seen exactly once, so every") print("rare individual is a singleton and the coverage estimate is exactly") print("zero. The worst cell in the table is the broken stick at n = 50,000.") print("The practical consequence is small, because Chao1 and ACE agree") print("closely at that effort anyway, but a package that reports 'ACE'") print("without saying which branch it took is hiding a real substitution.") # ------------------------------------------------------- convergence run head("PART 9. MONTE CARLO CONVERGENCE") print() print(f"A separate run: lognormal-1.0, n = {CONV_EFFORT}, " f"{CONV_REPLICATES:,} replicate communities, tracking the running mean") print("of each estimator so the Monte Carlo error can be seen shrinking.") print() t0 = time.time() conv = run_config("convergence", "log-normal", 1.0, children[len(CONFIGS) + 3], CONV_REPLICATES, efforts=[CONV_EFFORT]) print(f"convergence run took {time.time() - t0:.1f} s") print() j_conv = EFFORTS.index(CONV_EFFORT) conv_est = conv["est"][0] # (reps, 6) checkpoints = [250, 500, 1000, 2000, 4000, 8000, 15000, 25000, 40000, 60000, 85000, 115000, 150000] final = {} conv_series = {} for e, ename in enumerate(EST_NAMES): a = conv_est[:, e] run_mean = np.cumsum(a) / np.arange(1, a.size + 1) run_sd = np.array([a[:k].std(ddof=1) if k > 1 else 0.0 for k in checkpoints]) pts = [(k, float(run_mean[k - 1]), float(run_sd[i] / math.sqrt(k))) for i, k in enumerate(checkpoints)] conv_series[ename] = pts final[ename] = (float(run_mean[-1]), float(a.std(ddof=1) / math.sqrt(a.size))) print(f"Running mean of each estimator at n = {CONV_EFFORT}. The final column") print("is how many standard errors the checkpoint sits from the final value.") print() for ename in EST_NAMES: fin, fse = final[ename] print(f" {ename} final = {fin:.4f} +/- {fse:.4f} " f"bias = {fin - S_TRUE:+.4f}") print() print(f"{'reps':>7s} " + " ".join(f"{e:>17s}" for e in EST_NAMES)) print(rule(width=115)) for i, k in enumerate(checkpoints): cells = [] for ename in EST_NAMES: kk, mu, se = conv_series[ename][i] cells.append(f"{mu:8.3f}+/-{se:5.3f}") print(f"{k:>7d} " + " ".join(f"{c:>17s}" for c in cells)) print(rule(width=115)) print() max_dev = 0.0 max_where = "" for ename in EST_NAMES: fin, _ = final[ename] for kk, mu, se in conv_series[ename]: if se > 0: dev = abs(mu - fin) / se if dev > max_dev: max_dev = dev max_where = f"{ename} at {kk} replicates" print(f"largest checkpoint departure from the final value: " f"{max_dev:.2f} SE ({max_where})") print("Every checkpoint sits inside its own two-standard-error envelope of") print("the final value except where noted, and the envelope shrinks as") print("1/sqrt(replicates), which is the only rate Monte Carlo ever offers.") print() print("Cross-check between the two independent runs of the same cell:") for e, ename in enumerate(EST_NAMES): a = results["lognormal-1.0"]["est"][j_conv, :, e] m1, s1 = mean_se(a) m2, s2 = final[ename] d = m1 - m2 sd = math.sqrt(s1 ** 2 + s2 ** 2) print(f" {ename:<6s} sweep {m1:8.4f} +/- {s1:.4f} " f"convergence run {m2:8.4f} +/- {s2:.4f} " f"difference {d:+.4f} ({d / sd:+.2f} SE)") # ------------------------------------------------- sensitivity to design head("PART 10. DOES THE QUADRAT COUNT CHANGE THE ANSWER?") print() print("Chao2 and anything else built on incidence depends on how the same n") print("individuals are cut into sampling units. We re-ran lognormal-1.0 at") print("n = 1000 with the effort split four ways, 1000 replicates each.") print() rng_t = np.random.default_rng(20241104 + 77) print(f"{'T quadrats':>11s} {'m per quadrat':>14s} {'Chao2 mean':>11s} " f"{'SE':>7s} {'bias':>9s} {'RMSE':>8s}") print(rule()) tsens = [] for T in [4, 10, 25, 50, 100]: m = 1000 // T vals = np.empty(1000) for r in range(1000): p = make_community(rng_t, "log-normal", 1.0) q = rng_t.multinomial(m, p, size=T) counts = q.sum(axis=0) incid = (q > 0).sum(axis=0) s_obs = int((counts > 0).sum()) q1 = int((incid == 1).sum()) q2 = int((incid == 2).sum()) vals[r] = chao2(s_obs, q1, q2, T) mu, se = mean_se(vals) r = float(math.sqrt(np.mean((vals - S_TRUE) ** 2))) print(f"{T:>11d} {m:>14d} {mu:>11.3f} {se:>7.3f} " f"{mu - S_TRUE:>+9.3f} {r:>8.3f}") tsens.append((T, m, mu, se, r)) print(rule()) print() print("Same 1000 individuals, same community, same estimator. Splitting the") print("effort into 4 units instead of 100 changes the Chao2 bias from") print("-1.51 to -2.46 species, a 63 percent swing, while the RMSE moves the") print("other way and is worst at T = 100. Finer units give a less biased and") print("noisier answer; coarser units give the reverse. There is no free") print("choice here, and T = 10 is ours. Everything we say about Chao2 in") print("this study is conditional on it. In a real survey the number of") print("quadrats is usually set by how much walking a person can do in a day,") print("which is not a statistical criterion at all.") # ------------------------------------------------------------- figure data FIGDATA["seed"] = SEED FIGDATA["S_true"] = S_TRUE FIGDATA["efforts"] = EFFORTS FIGDATA["configs"] = [c[0] for c in CONFIGS] FIGDATA["estimators"] = EST_NAMES FIGDATA["evenness"] = {name: float(results[name]["J"].mean()) for name, _, _ in CONFIGS} FIGDATA["rarefaction"] = { name: { "sim": [float(results[name]["est"][j, :, 0].mean()) for j in range(len(EFFORTS))], "analytic": [float(results[name]["analytic"][j, :].mean()) for j in range(len(EFFORTS))], } for name, _, _ in CONFIGS } FIGDATA["bias"] = { name: {en: [float(table[(name, n, en)][2]) for n in EFFORTS] for en in EST_NAMES} for name, _, _ in CONFIGS } FIGDATA["se"] = { name: {en: [float(table[(name, n, en)][1]) for n in EFFORTS] for en in EST_NAMES} for name, _, _ in CONFIGS } FIGDATA["rmse"] = { name: {en: [float(table[(name, n, en)][3]) for n in EFFORTS] for en in EST_NAMES} for name, _, _ in CONFIGS } FIGDATA["cross10"] = { name: {en: headline[(name, en)][1] for en in EST_NAMES} for name, _, _ in CONFIGS } FIGDATA["convergence"] = {en: conv_series[en] for en in EST_NAMES} FIGDATA["convergence_final"] = {en: final[en] for en in EST_NAMES} FIGDATA["hurlbert"] = hurl_rows FIGDATA["tsens"] = tsens FIGDATA["ace_fallback"] = { name: [float(100 * results[name]["fallbacks"][j] / REPLICATES) for j in range(len(EFFORTS))] for name, _, _ in CONFIGS } head("PART 11. MACHINE-READABLE FIGURE DATA") print() print("Everything the article's five figures are drawn from, so a reader can") print("redraw them without rerunning the simulation.") print() print("BEGIN-FIGURE-JSON") print(json.dumps(FIGDATA, separators=(",", ":"), sort_keys=True)) print("END-FIGURE-JSON") head("DONE") print(f"total wall time {time.time() - t_start:.1f} s") print(f"seed {SEED}; rerunning this file reproduces every number above") print("exactly, on any machine with the same numpy version.") if __name__ == "__main__": main()