"""Analysis 2: Carter's hard-steps model. (a) posterior for the number of hard steps n from the time remaining in the window, (b) the same using also the early start (first step) - jointly, not separately, (c) effect of fixed Earth-system delays (Watson 2025), (d) goodness of fit of candidate transition sets to the uniform-order-statistics prediction, with a power calculation that shows how little such a test can exclude. Hard-steps facts used (all rates << 1/T, conditioned on success before T): the n step times are distributed as the order statistics of n U(0,T) variates, i.e. the joint density is n!/T^n on the ordered simplex; the last step has density n t^(n-1)/T^n. """ from __future__ import annotations import os import numpy as np import pandas as pd from scipy import stats import config as cfg NMAX = 20 def like_last(n, tobs, T, delay=0.0): x = (tobs - delay) / (T - delay) return n * x ** (n - 1) / (T - delay) def like_first_last(n, tLp, tobs, T): """abiogenesis is the first of n hard steps and happened before tLp; last step at tobs.""" n = np.asarray(n, float) out = n * (tobs ** (n - 1) - (tobs - tLp) ** (n - 1)) / T ** n return np.where(n >= 2, out, 0.0) def summarise_posterior(ns, post): post = post / post.sum() cdf = np.cumsum(post) q = lambda p: int(ns[np.searchsorted(cdf, p)]) pos = post[:10][post[:10] > 0] # n = 1 has zero likelihood in the first+last case return dict(mode=int(ns[np.argmax(post)]), mean=float((ns * post).sum()), median=q(0.5), lo90=q(0.05), hi90=q(0.95), max_over_min_like_n_le_10=float(pos.max() / pos.min()), post=[float(p) for p in post]) def n_posteriors(rng): ns = np.arange(1, NMAX + 1) t0 = cfg.T0_ZIRCON tobs = t0 rows = [] out = {} M = 200_000 cases = {} for lab, F in cfg.FUTURE.items(): cases[f"future {lab}"] = (np.full(1, tobs), np.full(1, tobs + F), np.full(1, tobs - 4.20)) t0d = cfg.draw_t0(rng, M) for which in ("short", "long", "mix"): Fd = cfg.draw_future(rng, M, which) luca = np.minimum(cfg.split_normal(rng, *cfg.LUCA_PAPER, M), t0d - 0.005) cases[f"marginalised: t0~U(4.404,4.50), future '{which}'"] = (t0d, t0d + Fd, t0d - luca) for lab, (to, T, tl) in cases.items(): # likelihoods averaged over date draws (marginal likelihood for each n) L_last = np.array([like_last(n, to, T).mean() for n in ns]) L_last_delay = np.array([like_last(n, to, T, 1.5).mean() for n in ns]) L_fl = np.array([like_first_last(n, tl, to, T).mean() for n in ns]) x = float(np.mean(to / T)) entry = {"x_mean": x, "n_mle_last_step": float(-1 / np.log(x)), "n_moments_last": float(np.mean(T) / np.mean(T - to) - 1), "n_moments_first": float(np.mean(T) / np.mean(tl) - 1)} for plab, prior in [("flat prior n=1..20", np.ones(NMAX)), ("1/n prior n=1..20", 1.0 / ns)]: entry[f"last step only | {plab}"] = summarise_posterior(ns, L_last * prior) entry[f"last step, 1.5 Gyr delay | {plab}"] = summarise_posterior(ns, L_last_delay * prior) entry[f"first (abiogenesis hard) + last | {plab}"] = summarise_posterior(ns, L_fl * prior) out[lab] = entry for n in ns: rows.append(dict(case=lab, n=n, like_last=L_last[n - 1], like_last_delay1p5=L_last_delay[n - 1], like_first_last=L_fl[n - 1])) pd.DataFrame(rows).to_csv(os.path.join(cfg.DATA, "a2_n_likelihoods.csv"), index=False) return out # --------------------------------------------------------------------------------- # goodness of fit of candidate transition sets # --------------------------------------------------------------------------------- CANDIDATE_SETS = { # label: (start Ga, [ages Ga of the steps, oldest first], future Gyr or None -> scan) "2 steps: life (LUCA 4.2), humans": (cfg.T0_ZIRCON, [4.20, 0.0]), "3 steps (Watson 2025): O2 photosynthesis 2.99, eukaryotes 1.735, humans": (cfg.T0_ZIRCON, [2.99, 1.735, 0.0]), "4 steps (Snyder-Beattie+ 2021 dates): life 3.8, eukaryotes 1.8, sex 1.2, intelligence": (4.2, [3.8, 1.8, 1.2, 0.0]), "4 steps (2024-26 dates): life 4.2, eukaryotes 1.735, animals 0.60, humans": (cfg.T0_ZIRCON, [4.20, 1.735, 0.60, 0.0]), "5 steps (Mills+ 2025 candidates): life 4.2, O2 photosynthesis 2.99, eukaryotes 1.735, animals 0.60, humans": (cfg.T0_ZIRCON, [4.20, 2.99, 1.735, 0.60, 0.0]), } def gof_stats(u): """u: (M, n) sorted points in (0,1). Returns dict of statistics, each (M,).""" M, n = u.shape sp = np.diff(np.concatenate([np.zeros((M, 1)), u, np.ones((M, 1))], axis=1), axis=1) # n+1 spacings k = np.arange(1, n + 1) dplus = np.max(k / n - u, axis=1) dminus = np.max(u - (k - 1) / n, axis=1) return dict(first=u[:, 0], last_gap=1 - u[:, -1], greenwood=np.sum(sp ** 2, axis=1), ks=np.maximum(dplus, dminus), max_spacing=sp.max(axis=1)) def two_sided_p(sim, obs): lo = np.mean(sim <= obs) hi = np.mean(sim >= obs) return float(min(1.0, 2 * min(lo, hi))) def goodness_of_fit(rng): M = 200_000 rows = [] for lab, (t0, ages) in CANDIDATE_SETS.items(): n = len(ages) sim = gof_stats(np.sort(rng.random((M, n)), axis=1)) for flab, F in cfg.FUTURE.items(): T = t0 + F u = (t0 - np.array(ages)) / T obs = gof_stats(u[None, :]) row = dict(candidate_set=lab, n=n, future=flab, T=T, u=" ".join(f"{x:.3f}" for x in u), expected_u=" ".join(f"{(k) / (n + 1):.3f}" for k in range(1, n + 1))) row["p_first_step_this_early"] = float(1 - (1 - u[0]) ** n) # exact, one-sided row["p_last_gap_this_large"] = float(u[-1] ** n) # exact, one-sided for s in ("greenwood", "ks", "max_spacing"): row[f"p_two_sided_{s}"] = two_sided_p(sim[s], obs[s][0]) # one-sided: are the steps MORE evenly spaced than uniform order statistics would be? row["p_too_even_greenwood"] = float(np.mean(sim["greenwood"] <= obs["greenwood"][0])) rows.append(row) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a2_goodness_of_fit.csv"), index=False) # ---- power: how often would these tests reject the hard-steps null (alpha = 0.05) if the # truth were something else? n = 5 transitions, T = 5.3 Gyr. n, T = 5, 5.304 null = gof_stats(np.sort(rng.random((M, n)), axis=1)) crit = {s: (np.quantile(null[s], 0.025), np.quantile(null[s], 0.975)) for s in ("greenwood", "ks")} alts = {} K = 50_000 # (i) environmental pacing with Earth-like gates: transitions at Earth's observed positions +- 10% jitter base = (cfg.T0_ZIRCON - np.array([4.20, 2.99, 1.735, 0.60, 0.0])) / T jit = np.clip(base[None, :] * np.exp(rng.normal(0, 0.10, (K, n))), 0, 0.999) alts["pacing, Earth-like gate times, 10% scatter"] = np.sort(jit, axis=1) # (ii) pacing with evenly spaced gates even = np.clip((np.arange(1, n + 1) / (n + 1))[None, :] + rng.normal(0, 0.03, (K, n)), 0.001, 0.999) alts["pacing, evenly spaced gates, 3% scatter"] = np.sort(even, axis=1) # (iii) all steps easy: exponential waits of mean 0.1 Gyr each (finish early) easy = np.cumsum(rng.exponential(0.1, (K, n)), axis=1) / T alts["all steps fast (mean wait 0.1 Gyr)"] = np.clip(easy, 0, 0.999) # (iv) one long deterministic delay then quick steps (gate opens at 80% of the window) late = 0.80 + np.cumsum(rng.exponential(0.02, (K, n)), axis=1) alts["single late gate at 0.8T then fast steps"] = np.sort(np.clip(late, 0, 0.999), axis=1) power = {} for lab, u in alts.items(): st = gof_stats(u) rej_g = (st["greenwood"] < crit["greenwood"][0]) | (st["greenwood"] > crit["greenwood"][1]) rej_k = (st["ks"] < crit["ks"][0]) | (st["ks"] > crit["ks"][1]) power[lab] = dict(power_greenwood=float(rej_g.mean()), power_ks=float(rej_k.mean())) return df, power def run(rng): out = {"n_posteriors": n_posteriors(rng)} df, power = goodness_of_fit(rng) out["gof"] = df.to_dict(orient="records") out["gof_power_n5_alpha0.05"] = power out["published_n"] = {"Carter 1983": "1-2", "Hanson 1998": "up to ~5-7", "Watson 2008": "4 (<10)", "Carter 2008": "5-6", "Hanson+ 2021": "3.9 from remaining time, 12.5 from abiogenesis delay", "Graham+ 2024": "~2.4 with 1.6-1.86 Gyr remaining", "Watson 2025": "2-3 with 1.5 Gyr delay", "Mills+ 2025": "0 (no hard steps)"} return out if __name__ == "__main__": import json r = run(np.random.default_rng(cfg.SEED)) for k, v in r["n_posteriors"].items(): print(k, {a: (b if not isinstance(b, dict) else {x: y for x, y in b.items() if x != "post"}) for a, b in v.items()}) print(pd.DataFrame(r["gof"]).drop(columns=["candidate_set"]).round(3).to_string()) print(json.dumps(r["gof_power_n5_alpha0.05"], indent=1))