"""Analysis 5: posterior predictive for f_bio under four labelled prior families. f_bio := probability that every step up to technological intelligence completes within the habitable window T of an Earth-like planet (same window as Earth's), given the model parameters. The posterior is over the parameters, given Earth's timeline and the fact that we can only be looking at a timeline on which every step succeeded (selection-conditioned likelihood). ANTHROPIC CONVENTION (stated explicitly after the referee report). The main ("post_") columns divide the likelihood by P(all steps succeed within T | parameters). That is the convention of Carter, Spiegel & Turner, Kipping and Snyder-Beattie et al.: hypotheses are NOT rewarded for producing more observers. The "sia_" columns drop the division: posterior proportional to prior x p(timeline AND success | parameters). That is what the self-indication assumption gives (prior weight proportional to the number of observers, which cancels the division exactly when the rates are universal), and equivalently what one gets by treating Earth as a random draw from Earth-like planets (Whitmire 2025; Baak & Snoek 2026). The two sets of quantiles answer different questions and must not be mixed. Families A optimist / pacing, Earth typical: the time to observers is set by planetary pacing, log10(time) ~ N(log10 t_obs, s^2) across Earth-like planets, s ~ U(0.05, 0.4) dex. (Earth is ASSUMED to be the median planet: nothing is learned from the timeline.) A' pacing, centre unknown: log10 of the median time ~ U(-1, 2) [0.1 .. 100 Gyr]; Earth's t_obs updates it with the truncated-normal likelihood phi(.)/Phi(.). B few hard steps, moderate priors: n = 2 or 3 sequential Poisson steps, each rate log-uniform on [1e-3, 1e3] per Gyr. C several hard steps, wide priors: n = 4 or 5 steps, each rate log-uniform on [1e-10, 1e14] per Gyr (the Snyder-Beattie et al. range). D Sandberg-Drexler-Ord style: abiogenesis rate log-normal with sigma = 50 dex around one event per window; the remaining steps collapse into q = f_i * f_c with f_i ~ LU[1e-3, 1], f_c ~ LU[1e-2, 1], q = P(intelligence within T | life). K reference: Kipping (2020) two-step model with Jeffreys priors in f for both rates. """ from __future__ import annotations import os import numpy as np import pandas as pd from scipy import stats import config as cfg from timing_lib import ess, log_hypoexp_cdf, log_hypoexp_cdf2, loglike_steps, weighted_quantile QS = [0.05, 0.16, 0.25, 0.50, 0.75, 0.84, 0.95] LN10 = np.log(10.0) STEP_AGES = {2: [], 3: ["LECA"], 4: ["LECA", "METAZOA"], 5: ["OXYPHOTO", "LECA", "METAZOA"]} def draw_dates(rng, M, n, fut="mix"): t0 = cfg.draw_t0(rng, M) T = t0 + cfg.draw_future(rng, M, fut) luca = np.minimum(cfg.split_normal(rng, *cfg.LUCA_PAPER, M), t0 - 0.005) cols = [t0 - cfg.draw_uniform(rng, getattr(cfg, k), M) for k in STEP_AGES[n]] + [t0] return t0, T, t0 - luca, np.stack(cols, axis=1) def summarise(label, logf_prior, logf_post, w, extra=None): """logf_* are natural logs of f_bio; returns one table row (log10 quantiles).""" row = {"family": label} lp = logf_prior / LN10 lq = logf_post / LN10 for q, v in zip(QS, np.quantile(lp, QS)): row[f"prior_q{int(q * 100):02d}"] = float(v) for q, v in zip(QS, weighted_quantile(lq, w, QS)): row[f"post_q{int(q * 100):02d}"] = float(v) for thr in (-1, -2, -6, -12, -22): row[f"prior_P_log10f_gt_{thr}"] = float((lp > thr).mean()) row[f"post_P_log10f_gt_{thr}"] = float((w * (lq > thr)).sum() / w.sum()) # information gained about log10 f_bio (bits): Kullback-Leibler divergence of the posterior # from the prior of this ONE marginal, estimated on prior-quantile bins. A coarse estimate # (binning can only lower the true value; finer bins add sampling noise), so it is given for # 20, 40 and 80 bins. "KL_post_vs_prior_bits" is the 40-bin value used in the first pass. def kl(nb): edges = np.unique(np.quantile(lp, np.linspace(0, 1, nb + 1))) if len(edges) <= 3: return 0.0 hp = np.histogram(lp, edges)[0].astype(float) hq = np.histogram(np.clip(lq, edges[0], edges[-1]), edges, weights=w)[0] hp, hq = hp / hp.sum(), hq / hq.sum() m = (hq > 0) & (hp > 0) return float(np.sum(hq[m] * np.log2(hq[m] / hp[m]))) row["KL_post_vs_prior_bits"] = kl(40) row["KL_bits_20bins"] = kl(20) row["KL_bits_80bins"] = kl(80) row["ESS"] = float(ess(w)) if extra: row.update(extra) return row def family_A(rng, M=200_000): t0, T, _, _ = draw_dates(rng, M, 2) s = rng.uniform(0.05, 0.4, M) d = np.log10(T / t0) / s # Earth's margin in sigma units logf_A = stats.norm.logcdf(d) # centre = Earth rowA = summarise("A pacing, Earth typical (assumed; no update)", logf_A, logf_A, np.ones(M)) c = rng.uniform(-1.0, 2.0, M) # log10 median time to observers (Gyr) z = (np.log10(T) - c) / s logf = stats.norm.logcdf(z) loglike = stats.norm.logpdf((np.log10(t0) - c) / s) - np.log(s) - logf w = np.exp(loglike - loglike.max()) rowAp = summarise("A' pacing, centre unknown (0.1-100 Gyr)", logf, logf, w) return [rowA, rowAp], {"A": (logf_A, logf_A, np.ones(M)), "A'": (logf, logf, w)} def steps_family(rng, n, lo, hi, label, M=200_000, Mprior=100_000): """n sequential Poisson steps with rates log-uniform on [10^lo, 10^hi] per Gyr.""" t0, T, tl, s_later = draw_dates(rng, M, n) # proposal = prior restricted to rates for which exp(-rate * gap) is not numerically zero gaps = np.diff(np.concatenate([np.zeros((M, 1)), s_later], axis=1), axis=1) gaps[:, 0] = s_later[:, 0] - tl his = np.full(n, float(hi)) his[1:] = np.minimum(hi, np.log10(400.0 / gaps.min(axis=0))) u = lo + (his - lo)[None, :] * rng.random((M, n)) lam = 10.0 ** u ll, logF = loglike_steps(lam, tl, s_later, T, first_is_bound=True, return_logF=True) w = np.nan_to_num(np.exp(ll - np.nanmax(ll))) # prior sample Tp = cfg.draw_t0(rng, Mprior) + cfg.draw_future(rng, Mprior, "mix") lamp = 10.0 ** (lo + (hi - lo) * rng.random((Mprior, n))) logFp = log_hypoexp_cdf2(lamp[:, 0], lamp[:, 1], Tp) if n == 2 else log_hypoexp_cdf(lamp, Tp) extra = {"post_P_flife_gt_0.5": float((w * (lam[:, 0] * T > np.log(2))).sum() / w.sum()), "prior_P_flife_gt_0.5": float((lamp[:, 0] * Tp > np.log(2)).mean()), "post_P_any_later_step_slower_than_window": float((w * (lam[:, 1:] * T[:, None] < 1).any(axis=1)).sum() / w.sum())} return summarise(label, logFp, logF, w, extra), (logFp, logF, w) def family_D(rng, M=300_000): def prior(m): T = cfg.draw_t0(rng, m) + cfg.draw_future(rng, m, "mix") x = np.clip(rng.normal(0.0, 50.0, m), -290, 290) # log10(lambda_L T) lamL = 10.0 ** x / T q = 10.0 ** rng.uniform(-3, 0, m) * 10.0 ** rng.uniform(-2, 0, m) lamI = -np.log1p(-np.clip(q, 0, 1 - 1e-12)) / T return T, lamL, lamI Tp, lLp, lIp = prior(M) logFp = log_hypoexp_cdf2(lLp, lIp, Tp) t0, T, tl, s_later = draw_dates(rng, M, 2) _, lL, lI = prior(M) lL, lI = lL * Tp[:M] / T, lI * Tp[:M] / T # same dimensionless draws, this sample's T ll, logF = loglike_steps(np.stack([lL, lI], axis=1), tl, s_later, T, True, return_logF=True) w = np.nan_to_num(np.exp(ll - np.nanmax(ll))) extra = {"post_P_flife_gt_0.5": float((w * (lL * T > np.log(2))).sum() / w.sum()), "prior_P_flife_gt_0.5": float((lLp * Tp > np.log(2)).mean())} return summarise("D Sandberg-Drexler-Ord style (sigma = 50 dex abiogenesis)", logFp, logF, w, extra), (logFp, logF, w) def family_K(rng, M=300_000): """Kipping 2020: Jeffreys prior in f = 1 - exp(-lambda T) for each of two rates.""" def prior(m): T = cfg.draw_t0(rng, m) + cfg.draw_future(rng, m, "mix") f = np.clip(rng.beta(0.5, 0.5, (m, 2)), 1e-300, 1 - 1e-15) return T, -np.log1p(-f) / T[:, None] Tp, lamp = prior(M) logFp = log_hypoexp_cdf2(lamp[:, 0], lamp[:, 1], Tp) t0, T, tl, s_later = draw_dates(rng, M, 2) f = np.clip(rng.beta(0.5, 0.5, (M, 2)), 1e-300, 1 - 1e-15) lam = -np.log1p(-f) / T[:, None] ll, logF = loglike_steps(lam, tl, s_later, T, True, return_logF=True) w = np.nan_to_num(np.exp(ll - np.nanmax(ll))) return summarise("K reference: Kipping 2020 two-step, Jeffreys priors in f", logFp, logF, w), (logFp, logF, w) # --------------------------------------------------------------------------------- # Sensitivity to the anthropic convention: posterior WITHOUT the division by P(success) # --------------------------------------------------------------------------------- def _sia_row(logf, logw): w = np.nan_to_num(np.exp(logw - np.nanmax(logw))) q = weighted_quantile(logf / LN10, w, QS) row = {f"sia_q{int(a * 100):02d}": float(v) for a, v in zip(QS, q)} for thr in (-1, -2, -6): row[f"sia_P_log10f_lt_{thr}"] = float((w * (logf / LN10 < thr)).sum() / w.sum()) row["sia_ESS"] = float(ess(w)) return row, w def steps_family_sia(rng, n, lo, hi, M=200_000): """n Poisson steps, rates log-uniform on [10^lo, 10^hi] per Gyr, posterior proportional to prior x p(timeline | rates) with NO division by P(success). With a log-uniform prior the posterior density of each later rate is ~ exp(-rate x gap), so the proposal draws later rates from exponentials (rate parameter = smallest gap over the date draws) and the abiogenesis rate log-uniformly; importance weights correct the rest exactly.""" t0, T, tl, s_later = draw_dates(rng, M, n) gaps = np.diff(np.concatenate([np.zeros((M, 1)), s_later], axis=1), axis=1) gaps[:, 0] = s_later[:, 0] - tl # smallest possible gap after abiogenesis gmin = gaps.min(axis=0) lam = np.empty((M, n)) lo1 = max(lo, -7.0) # below this the weight is < 1e-7 of the peak lam[:, 0] = 10.0 ** (lo1 + (hi - lo1) * rng.random(M)) lam[:, 1:] = rng.exponential(1.0, (M, n - 1)) / gmin[None, :] ll, logF = loglike_steps(lam, tl, s_later, T, first_is_bound=True, return_logF=True) num = ll + logF # log p(timeline | rates), unconditioned later = lam[:, 1:] logw = num + np.sum(-np.log(later) + later * gmin[None, :], axis=1) inside = ((later >= 10.0 ** lo) & (later <= 10.0 ** hi)).all(axis=1) logw = np.where(inside, logw, -np.inf) row, w = _sia_row(logF, logw) row["sia_P_flife_gt_0.5"] = float((w * (lam[:, 0] * T > np.log(2))).sum() / w.sum()) return row, (logF, w) def sia_grid_check(lo=-3.0, hi=3.0): """Deterministic cross-check of steps_family_sia for two steps at central dates (t0 = 4.404, LUCA 4.20, window +0.9 and +1.73 Gyr): 2-D grid in log rates, no Monte Carlo.""" out = {} g1 = np.linspace(lo, hi, 1201) g2 = np.linspace(lo, hi, 1201) L1, L2 = np.meshgrid(10.0 ** g1, 10.0 ** g2, indexing="ij") lam = np.stack([L1.ravel(), L2.ravel()], axis=1) for F in (0.9, 1.73): t0 = cfg.T0_ZIRCON ll, logF = loglike_steps(lam, t0 - 4.20, np.full((1, 1), t0), t0 + F, True, return_logF=True) for lab, lw in (("conditioned_on_success", ll), ("unconditioned_sia", ll + logF)): w = np.exp(lw - lw.max()) out[f"future {F} | {lab}"] = [float(v) for v in weighted_quantile(logF / LN10, w, [0.05, 0.5, 0.95])] return out def sia_all(rng, samples): """Unconditioned ('self-indication') posterior quantiles of log10 f_bio for every family.""" rows = {} keep = {} M = 200_000 # A: Earth assumed typical -> nothing is learned under either convention lpA = samples["A"][0] rows["A "] = {f"sia_q{int(a * 100):02d}": float(v) for a, v in zip(QS, np.quantile(lpA / LN10, QS))} for thr in (-1, -2, -6): rows["A "][f"sia_P_log10f_lt_{thr}"] = float((lpA / LN10 < thr).mean()) rows["A "]["sia_ESS"] = float(len(lpA)) keep["A"] = (lpA, np.ones(len(lpA))) # A': pacing with unknown centre; likelihood phi((log t_obs - c)/s)/s without the 1/Phi factor t0, T, _, _ = draw_dates(rng, M, 2) sg = rng.uniform(0.05, 0.4, M) c = rng.uniform(-1.0, 2.0, M) logf = stats.norm.logcdf((np.log10(T) - c) / sg) rows["A'"], w = _sia_row(logf, stats.norm.logpdf((np.log10(t0) - c) / sg) - np.log(sg)) keep["A'"] = (logf, w) for n, lo, hi, key in [(2, -3, 3, "B2"), (3, -3, 3, "B3"), (4, -10, 14, "C4"), (5, -10, 14, "C5")]: rows[key], keep[key] = steps_family_sia(rng, n, lo, hi) # D: sigma = 50 dex abiogenesis prior, q = f_i f_c; importance-sample the fast half of the # log-normal (weights for lambda_L T < 1e-6 are < 1e-6 of the peak) and q log-uniform on [1e-5, 1] MD = 400_000 t0, T, tl, s_later = draw_dates(rng, MD, 2) x = rng.uniform(-6.0, 290.0, MD) # log10(lambda_L T) y = rng.uniform(-5.0, 0.0, MD) # log10 q lamL = 10.0 ** x / T lamI = -np.log1p(-np.clip(10.0 ** y, 0, 1 - 1e-12)) / T # prior density of y = log10 f_i + log10 f_c with U(-3,0) + U(-2,0): trapezoid on [-5, 0] py = np.where(y < -3, (y + 5) / 6, np.where(y < -2, 2 / 6, -y / 6)) ll, logF = loglike_steps(np.stack([lamL, lamI], axis=1), tl, s_later, T, True, return_logF=True) rows["D "], w = _sia_row(logF, ll + logF + stats.norm.logpdf(x, 0.0, 50.0) + np.log(np.maximum(py, 1e-300))) keep["D"] = (logF, w) # K: Kipping two-step, Jeffreys priors in f t0, T, tl, s_later = draw_dates(rng, M, 2) f = np.clip(rng.beta(0.5, 0.5, (M, 2)), 1e-300, 1 - 1e-15) lam = -np.log1p(-f) / T[:, None] ll, logF = loglike_steps(lam, tl, s_later, T, True, return_logF=True) rows["K "], w = _sia_row(logF, ll + logF) keep["K"] = (logF, w) # insensitivity to the prior bounds (compare a5_fbio_prior_bound_sensitivity.csv) sens = [] for n, lo, hi in [(2, -2, 2), (2, -6, 3), (3, -6, 3), (4, -3, 3), (4, -20, 14), (5, -3, 3), (5, -20, 14)]: r, _ = steps_family_sia(rng, n, lo, hi, M=120_000) sens.append(dict(family=f"{n} steps, rates LU[1e{lo},1e{hi}]/Gyr", **r)) return rows, keep, sens def profile_likelihood(): """Profile likelihood of f_bio in the two-step model: sup over (lambda_L, lambda_I) with F(T; lambda_L, lambda_I) = f of the selection-conditioned likelihood, relative to its global supremum. Central dates (t0 = 4.404, LUCA 4.2).""" targets = np.array([1e-12, 1e-6, 1e-3, 1e-2, 0.1, 0.3, 0.5, 0.7, 0.9, 0.99, 0.999]) gL = np.linspace(-9, 6, 601) gI = np.linspace(-14, 4, 2401) rows = [] curves = {} for F in (0.9, 1.73): t0 = cfg.T0_ZIRCON T, tl = t0 + F, t0 - 4.20 LL, LI = np.meshgrid(10.0 ** gL, 10.0 ** gI, indexing="ij") lam = np.stack([LL.ravel(), LI.ravel()], axis=1) ll, logF = loglike_steps(lam, tl, np.full((1, 1), t0), T, True, return_logF=True) ll = ll.reshape(LL.shape) logF = logF.reshape(LL.shape) sup = ll.max() # attained for lambda_I -> 0 and lambda_L ~ 30-40 / Gyr, # 0.4% above the lambda_L -> inf value 1/T fgrid = np.concatenate([10.0 ** np.linspace(-12, -1, 111), np.linspace(0.1, 0.999, 300)[1:]]) prof = np.full(len(fgrid), -np.inf) for i in range(len(gL)): # F is increasing in lambda_I along each row ok = np.log(fgrid) < logF[i].max() v = np.interp(np.log(fgrid[ok]), logF[i], ll[i]) prof[ok] = np.maximum(prof[ok], v) rel = np.exp(prof - sup) curves[F] = (fgrid, rel) for t in targets: rows.append(dict(future=F, f_bio=t, path="sup over both rates at fixed f_bio", relative_profile_likelihood=float(np.interp(np.log(t), np.log(fgrid), rel)))) # the fully hard corner (both rates -> 0, f_bio -> 0 along the diagonal): 2 t_L' / T^2 rows.append(dict(future=F, f_bio=0.0, path="both rates -> 0 together (fully hard corner)", relative_profile_likelihood=float(2 * tl / T ** 2 / np.exp(sup)))) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a5_profile_likelihood.csv"), index=False) return df, curves def run(rng): rows, samples = family_A(rng) for n, lo, hi, lab in [(2, -3, 3, "B2 2 hard steps, rates LU[1e-3,1e3]/Gyr"), (3, -3, 3, "B3 3 hard steps, rates LU[1e-3,1e3]/Gyr"), (4, -10, 14, "C4 4 hard steps, rates LU[1e-10,1e14]/Gyr"), (5, -10, 14, "C5 5 hard steps, rates LU[1e-10,1e14]/Gyr")]: r, s = steps_family(rng, n, lo, hi, lab) rows.append(r) samples[lab[:2]] = s r, s = family_D(rng) rows.append(r) samples["D"] = s r, s = family_K(rng) rows.append(r) samples["K"] = s # sensitivity of B and C to their own prior bounds (posterior median and 5-95%) sens = [] for n, lo, hi in [(2, -2, 2), (2, -6, 3), (3, -6, 3), (4, -3, 3), (4, -20, 14), (5, -3, 3)]: r, _ = steps_family(rng, n, lo, hi, f"{n} steps, rates LU[1e{lo},1e{hi}]/Gyr", M=120_000, Mprior=60_000) sens.append(r) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a5_fbio_quantiles.csv"), index=False) pd.DataFrame(sens).to_csv(os.path.join(cfg.DATA, "a5_fbio_prior_bound_sensitivity.csv"), index=False) prof, curves = profile_likelihood() # ---- sensitivity to the anthropic convention (uses the generator AFTER everything above, so # the first-pass numbers are unchanged) ---- sia_rows, sia_samples, sia_sens = sia_all(rng, samples) for r in rows: r.update(sia_rows[r["family"][:2]]) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a5_fbio_quantiles.csv"), index=False) pd.DataFrame(sia_sens).to_csv(os.path.join(cfg.DATA, "a5_fbio_sia_prior_bound_sensitivity.csv"), index=False) # data dictionary for the integrator (the two conventions must not be mixed) qtxt = "quantile of log10 f_bio" dd = [("family", "prior family label (A, A', B2, B3, C4, C5, D, K)"), ("prior_qXX", f"XX% {qtxt} under the prior alone"), ("post_qXX", f"XX% {qtxt}, posterior with the likelihood CONDITIONED ON SUCCESS: " "p(timeline | parameters) / P(all steps finish within the window | parameters). " "Carter / Kipping / Snyder-Beattie convention; no reward for producing more observers. " "Use only with a no-update (self-sampling) treatment of our own existence."), ("sia_qXX", f"XX% {qtxt}, posterior WITHOUT the division (self-indication, or Earth as a random draw " "from Earth-like planets). Use only inside a full self-indication calculation that also " "weights every other factor and the Fermi observation. Never combine with post_ columns."), ("prior_P_log10f_gt_T / post_P_log10f_gt_T", "probability that log10 f_bio exceeds T (prior / conditioned posterior)"), ("sia_P_log10f_lt_T", "probability that log10 f_bio is below T, unconditioned posterior"), ("KL_post_vs_prior_bits, KL_bits_20bins, KL_bits_80bins", "binned Kullback-Leibler divergence (bits) of the conditioned posterior from the prior of log10 f_bio, " "40 / 20 / 80 prior-quantile bins; a coarse estimate for this one marginal"), ("ESS, sia_ESS", "importance-sampling effective sample size of the conditioned / unconditioned posterior"), ("prior_P_flife_gt_0.5, post_P_flife_gt_0.5, sia_P_flife_gt_0.5", "probability that abiogenesis succeeds within the window more often than not"), ("post_P_any_later_step_slower_than_window", "conditioned posterior probability that at least one later step has rate < 1/T")] pd.DataFrame(dd, columns=["column", "meaning"]).to_csv( os.path.join(cfg.DATA, "a5_fbio_quantiles_columns.csv"), index=False) for k, v in sia_samples.items(): samples[k] = samples[k] + (v,) out = {"families": rows, "prior_bound_sensitivity": sens, "sia_prior_bound_sensitivity": sia_sens, "sia_grid_check_B2_central_dates": sia_grid_check(), "profile_likelihood": prof.to_dict(orient="records")} return out, samples, curves if __name__ == "__main__": out, samples, curves = run(np.random.default_rng(cfg.SEED)) pd.set_option("display.width", 250) df = pd.DataFrame(out["families"]).set_index("family") print(df.round(2).T.to_string()) print(pd.DataFrame(out["prior_bound_sensitivity"]).set_index("family").round(2).T.to_string()) print(pd.DataFrame(out["profile_likelihood"]).to_string())