"""Analysis 3: reproduce the core of Snyder-Beattie, Sandberg, Drexler & Bonsall (2021, Astrobiology 21:265) and test its sensitivity. Model (their Section 3): n sequential transitions with exponential waiting times of mean beta_i; likelihood = prod (1/beta_i) exp(-t_i/beta_i) / P(all n done within L); prior: log-uniform on each beta_i between 1e-14 Ga and 1e10 Ga. Their Table 2 reports, for nine combinations of transitions, the posterior weight on parameter values for which P(final transition within L) exceeds "10^-2, 10^-12, 10^-24". """ from __future__ import annotations import os import numpy as np import pandas as pd import config as cfg from timing_lib import ess, log_hypoexp_cdf, log_hypoexp_cdf2, loglike_steps LOG_LO, LOG_HI = -14.0, 10.0 # log10 beta (Ga), as stated in the paper # Table 2 of the paper: (transition dates in Gya incl. the 0th, window L, published weights) TABLE2 = [ ([4.3, 4.2, 0.0], 5.8, (0.068, 0.40, 0.83)), ([4.2, 3.8, 0.0], 5.2, (0.026, 0.25, 0.68)), ([4.2, 3.8, 1.8], 5.2, (0.041, 0.27, 0.70)), ([4.2, 1.8, 0.0], 5.2, (0.012, 0.18, 0.61)), ([4.3, 4.2, 1.8, 0.0], 5.8, (0.0064, 0.11, 0.43)), ([4.2, 3.8, 1.8, 0.0], 5.2, (0.0022, 0.059, 0.30)), ([4.2, 3.8, 1.8, 1.2], 5.2, (0.0056, 0.080, 0.34)), ([4.2, 3.8, 1.8, 1.2, 0.0], 5.2, (0.0003, 0.014, 0.11)), ([4.3, 4.2, 1.8, 1.7, 0.0], 5.8, (0.0027, 0.056, 0.26)), ] PRIOR_PUBLISHED = {2: (0.42, 0.66, 0.94), 3: (0.27, 0.51, 0.84), 4: (0.18, 0.39, 0.71)} TABLE3_PUBLISHED = { # log10 of the upper bound (Ga) needed 2: {"10%|P>1%": 4.8, "10%|P>10%": 3.1, "50%|P>1%": 1.8, "50%|P>10%": 1.2}, 3: {"10%|P>1%": 2.4, "10%|P>10%": 1.6, "50%|P>1%": 1.2, "50%|P>10%": 0.8}, 4: {"10%|P>1%": 1.8, "10%|P>10%": 1.2, "50%|P>1%": 1.0, "50%|P>10%": 0.6}, } def logF(lam, L): if lam.shape[1] == 2: return log_hypoexp_cdf2(lam[:, 0], lam[:, 1], L) return log_hypoexp_cdf(lam, L) def posterior_samples(waits, L, rng, M, lo=LOG_LO, hi=LOG_HI): """Importance samples of log10(beta) from the prior restricted to where the likelihood is not numerically zero (beta_i > wait_i/400, where exp(-t/beta) < e^-400), with weights.""" waits = np.asarray(waits, float) n = len(waits) lob = np.maximum(lo, np.log10(waits / 400.0)) u = lob + (hi - lob) * rng.random((M, n)) lam = 10.0 ** (-u) lf = logF(lam, L) ll = np.sum(np.log(lam) - lam * waits, axis=1) - lf w = np.exp(ll - ll.max()) return u, lf, w def reproduce_table2(rng, M=200_000): rows = [] thr10 = np.log(10.0) * np.array([-2, -12, -24]) thre = np.array([-2.0, -12.0, -24.0]) for times, L, pub in TABLE2: waits = -np.diff(times) n = len(waits) u, lf, w = posterior_samples(waits, L, rng, M) post10 = [(w * (lf > t)).sum() / w.sum() for t in thr10] poste = [(w * (lf > t)).sum() / w.sum() for t in thre] # reading C: the prior range (not the thresholds) is in natural-log units, i.e. the prior # actually spans 10^(-14 ln10) .. 10^(10 ln10) Ga, thresholds as labelled s = np.log(10.0) u2, lf2, w2 = posterior_samples(waits, L, rng, M, lo=LOG_LO * s, hi=LOG_HI * s) postC = [(w2 * (lf2 > t)).sum() / w2.sum() for t in thr10] # prior weights up = LOG_LO + (LOG_HI - LOG_LO) * rng.random((150_000, n)) lfp = logF(10.0 ** (-up), L) rows.append(dict(transitions_Gya=" ".join(str(t) for t in times), L=L, n=n, pub_1=pub[0], pub_2=pub[1], pub_3=pub[2], ours_10pow_1=post10[0], ours_10pow_2=post10[1], ours_10pow_3=post10[2], ours_epow_1=poste[0], ours_epow_2=poste[1], ours_epow_3=poste[2], ours_widerprior_1=postC[0], ours_widerprior_2=postC[1], ours_widerprior_3=postC[2], prior_10pow_1=(lfp > thr10[0]).mean(), prior_10pow_2=(lfp > thr10[1]).mean(), prior_10pow_3=(lfp > thr10[2]).mean(), prior_epow_1=(lfp > thre[0]).mean(), prior_epow_2=(lfp > thre[1]).mean(), prior_epow_3=(lfp > thre[2]).mean(), ESS=ess(w))) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a3_table2_reproduction.csv"), index=False) def rms(prefix): a = np.log10(df[[f"{prefix}_1", f"{prefix}_2", f"{prefix}_3"]].values.clip(1e-6)) b = np.log10(df[["pub_1", "pub_2", "pub_3"]].values) return float(np.sqrt(np.mean((a - b) ** 2))) summ = {"rms_log10_dev_thresholds_as_labelled_10^-k": rms("ours_10pow"), "rms_log10_dev_thresholds_e^-k": rms("ours_epow"), "rms_log10_dev_prior_range_in_e-folds": rms("ours_widerprior"), "max_abs_ratio_dev_e^-k": float(np.max(np.abs( df[["ours_epow_1", "ours_epow_2", "ours_epow_3"]].values / df[["pub_1", "pub_2", "pub_3"]].values - 1))), "prior_weights_published": PRIOR_PUBLISHED, "prior_weights_ours_10^-k": {int(n): [float(g[f"prior_10pow_{i}"].mean()) for i in (1, 2, 3)] for n, g in df.groupby("n")}, "prior_weights_ours_e^-k": {int(n): [float(g[f"prior_epow_{i}"].mean()) for i in (1, 2, 3)] for n, g in df.groupby("n")}, # smallest completion probability the stated prior can produce at all (n = 2, L = 5.2) "min_possible_P_n2_stated_prior": float(np.exp(log_hypoexp_cdf2(np.array([1e-10]), np.array([1e-10]), 5.2))[0]), } return df, summ def table3_check(rng, M=200_000): """Upper prior bound needed for 10% / 50% of posterior weight to have P(life) > 1% / 10%.""" out = {} rows_for_n = {2: TABLE2[1], 3: TABLE2[5], 4: TABLE2[7]} bgrid = np.arange(0.0, 10.0001, 0.05) for n, (times, L, _) in rows_for_n.items(): waits = -np.diff(times) u, lf, w = posterior_samples(waits, L, rng, M) umax = u.max(axis=1) res = {} for plab, thr in [("P>1%", np.log(0.01)), ("P>10%", np.log(0.10)), ("P>e^-2", -2.0), ("P>e^-1", -1.0)]: frac = np.array([(w * (lf > thr) * (umax <= b)).sum() / max((w * (umax <= b)).sum(), 1e-300) for b in bgrid]) for tgt in (0.10, 0.50): ok = np.where(frac >= tgt)[0] res[f"{int(tgt * 100)}%|{plab}"] = float(bgrid[ok[-1]]) if len(ok) else None out[n] = res return out def _updated_row(rng, model, ages, fut, lo, hi, first_bound, M): """One posterior under the Snyder-Beattie et al. prior (log10 beta uniform on [lo, hi], Ga) with 2024-26 dates marginalised. Returns a table row.""" n = len(ages) t0 = cfg.draw_t0(rng, M) T = t0 + cfg.draw_future(rng, M, fut) tl = t0 - np.minimum(cfg.split_normal(rng, *cfg.LUCA_PAPER, M), t0 - 0.005) s_later = (t0[:, None] - np.array(ages[1:])[None, :]) # proposal = prior restricted to where the likelihood is not numerically zero lob = np.full(n, float(lo)) waits_c = np.diff(np.concatenate([[0.0], cfg.T0_ZIRCON - np.array(ages)])) lob[1:] = np.maximum(lo, np.log10(waits_c[1:] / 400.0)) if not first_bound: lob[0] = max(lo, np.log10(waits_c[0] / 400.0)) u = lob + (hi - lob) * rng.random((M, n)) lam = 10.0 ** (-u) ll = loglike_steps(lam, tl, s_later, T, first_is_bound=first_bound) # importance correction for the restricted proposal is a constant -> cancels w = np.exp(ll - np.nanmax(ll)) w = np.nan_to_num(w) lf = logF(lam, T) if n > 2 else log_hypoexp_cdf2(lam[:, 0], lam[:, 1], T) row = dict(model=model, n=n, future=fut, log10_beta_min=lo, log10_beta_max=hi, abiogenesis="bounded by LUCA" if first_bound else "exactly at LUCA", post_P_gt_1em2=(w * (lf > np.log(1e-2))).sum() / w.sum(), post_P_gt_1em6=(w * (lf > np.log(1e-6))).sum() / w.sum(), post_P_gt_1em12=(w * (lf > np.log(1e-12))).sum() / w.sum(), post_median_log10_P=float(np.interp(0.5, np.cumsum(w[np.argsort(lf)]) / w.sum(), np.sort(lf))) / np.log(10), ESS=ess(w)) for i in range(n): row[f"post_P_beta{i + 1}_gt_T"] = (w * (u[:, i] > np.log10(T))).sum() / w.sum() row["post_P_any_beta_gt_T"] = (w * (u > np.log10(T)[:, None]).any(axis=1)).sum() / w.sum() # prior probability of the same event, for comparison (beta_1 > T under the log-uniform prior) row["prior_P_beta1_gt_T"] = float(np.mean(np.clip((hi - np.log10(T)) / (hi - lo), 0, 1))) return row MODELS = [("2 steps: life, humans", [4.20, 0.0]), ("3 steps: life, eukaryotes, humans", [4.20, 1.735, 0.0]), ("4 steps: life, eukaryotes, animals, humans", [4.20, 1.735, 0.60, 0.0]), ("5 steps: life, O2 photosynthesis, eukaryotes, animals, humans", [4.20, 2.99, 1.735, 0.60, 0.0])] def updated_analysis(rng, M=150_000): """Same machinery with 2024-26 dates: abiogenesis only bounded (before LUCA), eukaryotes = LECA, animals, humans; habitability from t0; several windows and prior upper bounds.""" rows = [] for model, ages in MODELS: for fut in ("short", "long", "mix"): for hi in (2.0, 4.0, 10.0, 20.0): if fut != "mix" and hi != 10.0: # short/long futures only at the published bound continue for first_bound in (True, False): rows.append(_updated_row(rng, model, ages, fut, LOG_LO, hi, first_bound, M)) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a3_updated_dates.csv"), index=False) return df def lower_bound_sensitivity(rng, M=150_000): """The statement 'abiogenesis probably does not have an expected time longer than the window' depends on the LOWER prior bound: 13 of the 24 decades of the published prior are faster than the data can resolve, and all of that volume counts as 'fast' (referee report). Scan it.""" rows = [] for model, ages in MODELS: for lo in (-14.0, -8.0, -4.0, -2.0, -1.0): rows.append(_updated_row(rng, model, ages, "mix", lo, LOG_HI, True, M)) df = pd.DataFrame(rows) df.to_csv(os.path.join(cfg.DATA, "a3_lower_bound_sensitivity.csv"), index=False) return df def table3_deviation(ours): """Recomputed minus published Table 3 bounds (dex) under the two readings of the thresholds.""" pairs = {"10%|P>1%": "10%|P>e^-2", "10%|P>10%": "10%|P>e^-1", "50%|P>1%": "50%|P>e^-2", "50%|P>10%": "50%|P>e^-1"} de, d10 = [], [] for n, pub in TABLE3_PUBLISHED.items(): for k, v in pub.items(): if ours[n].get(pairs[k]) is not None: de.append(ours[n][pairs[k]] - v) if ours[n].get(k) is not None: d10.append(ours[n][k] - v) de, d10 = np.array(de), np.array(d10) return {"reading_e^-k": {"n_cells": int(len(de)), "mean_dex": float(de.mean()), "min_dex": float(de.min()), "max_dex": float(de.max())}, "reading_10^-k_as_printed": {"n_cells": int(len(d10)), "mean_dex": float(d10.mean()), "min_dex": float(d10.min()), "max_dex": float(d10.max())}, "grid_resolution_dex": 0.05} def run(rng): df, summ = reproduce_table2(rng) out = {"table2": df.to_dict(orient="records"), "table2_summary": summ, "table3_published_log10_bound": TABLE3_PUBLISHED, "table3_ours_log10_bound": table3_check(rng)} out["table3_deviation_dex"] = table3_deviation(out["table3_ours_log10_bound"]) # what the paper's verbal statements become when the stated prior is taken literally out["verbal_statements"] = { f"{n}_steps": {"post_weight_below_1pct": [float(1 - g.ours_10pow_1.max()), float(1 - g.ours_10pow_1.min())], "post_weight_below_1e-12": [float(1 - g.ours_10pow_2.max()), float(1 - g.ours_10pow_2.min())], "post_weight_below_1e-24": [float(1 - g.ours_10pow_3.max()), float(1 - g.ours_10pow_3.min())]} for n, g in df.groupby("n")} upd = updated_analysis(rng) out["updated"] = upd.to_dict(orient="records") out["lower_bound_sensitivity"] = lower_bound_sensitivity(rng).to_dict(orient="records") return out if __name__ == "__main__": import json r = run(np.random.default_rng(cfg.SEED)) pd.set_option("display.width", 250) print(pd.DataFrame(r["table2"]).round(4).to_string()) print(json.dumps(r["table2_summary"], indent=1, default=float)) print(json.dumps(r["table3_published_log10_bound"])) print(json.dumps(r["table3_ours_log10_bound"])) u = pd.DataFrame(r["updated"]) print(u.round(4).to_string())