"""Exact ideal-Poisson count inference; no deadtime correction or model fitting. Run with bundled Python 3.12, or any Python with numpy and scipy. exposure_ratio always means T_off/T_on; E[M_off]=exposure_ratio*b_on. """ from __future__ import annotations import csv import hashlib import json import math import numbers import pathlib import platform import sys from datetime import datetime, timezone sys.dont_write_bytecode = True _PRIVATE_SCIPY = pathlib.Path("D:/Photometer_research_20261022/research/cloud_transport/private_readers") if _PRIVATE_SCIPY.is_dir(): sys.path.insert(0, str(_PRIVATE_SCIPY)) import numpy as np import scipy from scipy.stats import beta, binom, chi2, poisson SCHEMA = "IDEAL_POISSON_COUNT_INFERENCE_V1" START_UTC = "2026-10-22T08:03:13Z" SEED = 20261022 MC_FAMILY_ALPHA = .01 MC_COMPARISONS = 96 # 64 coverage + 32 null-detection probabilities EXPOSURE_RATIO_DEFINITION = "T_off/T_on; E[M_off]=exposure_ratio*b_on" def _count(value, name="count"): if isinstance(value, (bool, np.bool_)) or not isinstance(value, numbers.Integral) or value < 0: raise ValueError(f"{name} must be a nonnegative integer") return int(value) def _positive(value, name): if isinstance(value, (bool, np.bool_)): raise ValueError(f"{name} must be finite and positive") value = float(value) if not math.isfinite(value) or value <= 0: raise ValueError(f"{name} must be finite and positive") return value def _alpha(value): value = _positive(value, "alpha") if value >= 1: raise ValueError("alpha must lie strictly between zero and one") return value def _background(value): if isinstance(value, (bool, np.bool_)): raise ValueError("b_on must be finite and nonnegative") value = float(value) if not math.isfinite(value) or value < 0: raise ValueError("b_on must be finite and nonnegative") return value def alpha_for_search(alpha_family, number_of_tests): """Bonferroni per-test alpha; arbitrary dependence is allowed. number_of_tests is the actual number of hypotheses evaluated, including all channels/windows/templates whose union defines the false alarm event. """ a = _alpha(alpha_family) m = _count(number_of_tests, "number_of_tests") if m == 0: raise ValueError("number_of_tests must be positive") return a / m def count_threshold_known_background(b_on, alpha_per_test): """Smallest k>=1 with P(Poisson(b_on)>=k)<=alpha_per_test.""" b = _background(b_on) a = _alpha(alpha_per_test) loga = math.log(a) hi = max(1, int(math.ceil(b + 8 * math.sqrt(b + 1) + 32))) while poisson.logsf(hi - 1, b) > loga: hi *= 2 lo = 0 # k=0 has tail probability one, and cannot satisfy alpha<1. while hi - lo > 1: mid = (lo + hi) // 2 if poisson.logsf(mid - 1, b) <= loga: hi = mid else: lo = mid return hi def exact_onoff_pvalue(n_on, m_off, exposure_ratio): """One-sided conditional-binomial p value for no source. N~Pois(s+b), M~Pois(r*b), independent. Under s=0: N | (N+M=t) ~ Binomial(t,1/(1+r)). No estimated b is plugged in. """ n = _count(n_on, "n_on") m = _count(m_off, "m_off") r = _positive(exposure_ratio, "exposure_ratio") return float(binom.sf(n - 1, n + m, 1 / (1 + r))) def count_threshold_onoff(m_off, exposure_ratio, alpha_per_test): """Minimum n for the exact conditional test at the observed off count. This is an equivalent decision-rule convenience, not a conditional false-alarm guarantee with M fixed. Unconditional validity follows by conditioning on N+M, not by conditioning on M. """ m = _count(m_off, "m_off") r = _positive(exposure_ratio, "exposure_ratio") a = _alpha(alpha_per_test) p = 1 / (1 + r) loga = math.log(a) hi = max(1, int(math.ceil(m / r + 8 * math.sqrt(m / r + 1) + 32))) while binom.logsf(hi - 1, hi + m, p) > loga: hi *= 2 lo = 0 while hi - lo > 1: mid = (lo + hi) // 2 if binom.logsf(mid - 1, mid + m, p) <= loga: hi = mid else: lo = mid return hi def poisson_total_mean_upper(n, alpha=0.05): """One-sided Garwood upper bound: P_mu(X<=n)=alpha at mu=U.""" n = _count(n, "n") a = _alpha(alpha) return float(chi2.isf(a, 2 * (n + 1)) / 2) def poisson_total_mean_lower(m, alpha=0.05): """One-sided lower bound, zero for m=0.""" m = _count(m, "m") a = _alpha(alpha) return 0.0 if m == 0 else float(chi2.ppf(a, 2 * m) / 2) def red_source_upper_known_background(n_on, b_on, alpha=0.05): """At-least-(1-alpha) frequentist upper limit on expected source counts. b_on is a known fixed mean in this same on window. No source prior, plug-in background, conditioning on nondetection, or unified interval. """ return max(0.0, poisson_total_mean_upper(n_on, alpha) - _background(b_on)) def red_source_upper_finite_baseline(n_on, m_off, exposure_ratio, alpha=0.05, split=0.5): """Conservative source upper bound by Bonferroni subtraction. T_off/T_on=r; alpha_total=alpha*split, alpha_off=alpha*(1-split). split must be fixed before seeing these counts. The union-bound proof does not require independence of confidence-bound failure events. The Poisson on/off observation model itself assumes independent windows. """ n = _count(n_on, "n_on") m = _count(m_off, "m_off") r = _positive(exposure_ratio, "exposure_ratio") a = _alpha(alpha) f = _alpha(split) return max(0.0, poisson_total_mean_upper(n, a * f) - poisson_total_mean_lower(m, a * (1 - f)) / r) def _upper_vec(n, a): return chi2.isf(a, 2 * (np.asarray(n) + 1)) / 2 def _lower_vec(m, a): m = np.asarray(m) # Avoid chi-square df=0 evaluations; lower bound there is exactly zero. out = np.zeros(m.shape, dtype=float) use = m > 0 out[use] = chi2.ppf(a, 2 * m[use]) / 2 return out def _clopper_pearson(successes, trials, confidence=0.95): a = 1 - confidence lo = 0.0 if successes == 0 else float(beta.ppf(a / 2, successes, trials - successes + 1)) hi = 1.0 if successes == trials else float(beta.isf(a / 2, successes + 1, trials - successes)) return [lo, hi] def _known_exact_coverage(s, b, alpha): if s == 0: return 1.0 target = s + b hi = max(1, int(math.ceil(target + 8 * math.sqrt(target + 1) + 32))) if poisson_total_mean_upper(0, alpha) >= target: return 1.0 lo = 0 while hi - lo > 1: mid = (lo + hi) // 2 if poisson_total_mean_upper(mid, alpha) >= target: hi = mid else: lo = mid return float(poisson.sf(hi - 1, target)) def _finite_exact_coverage_bounds(s, b, r, alpha, split=0.5, omitted_mass_tolerance=1e-12): """Analytic finite sum with explicit omitted-Poisson-tail enclosure. No Monte Carlo is used here. The enclosure is numerical, not interval arithmetic certified against every possible floating-point error. """ nmax = int(poisson.isf(omitted_mass_tolerance / 4, s + b)) mmax = int(poisson.isf(omitted_mass_tolerance / 4, r * b)) ns = np.arange(nmax + 1) ms = np.arange(mmax + 1) un = _upper_vec(ns, alpha * split) lm = _lower_vec(ms, alpha * (1 - split)) / r limits = np.maximum(0, un[:, None] - lm[None, :]) # Tiny tolerance only prevents loss of exact-boundary points to roundoff. covered = limits >= s - 1e-12 mass = poisson.pmf(ns, s + b)[:, None] * poisson.pmf(ms, r * b)[None, :] inside = float(np.sum(mass[covered])) missing_n = float(poisson.sf(nmax, s + b)) missing_m = float(poisson.sf(mmax, r * b)) missing = missing_n + missing_m - missing_n * missing_m return {"coverage_lower": min(1.0, inside), "coverage_upper": min(1.0, inside + missing + 2e-13), "omitted_joint_mass": missing, "n_max": nmax, "m_max": mmax} def _conditional_exact_false_alarm_bounds(b, r, alpha, omitted_mass_tolerance=1e-12): tmax = int(poisson.isf(omitted_mass_tolerance / 4, (1 + r) * b)) ts = np.arange(tmax + 1) p = 1 / (1 + r) critical = binom.isf(alpha, ts, p) + 1 conditional_false = binom.sf(critical - 1, ts, p) total = float(np.sum(poisson.pmf(ts, (1 + r) * b) * conditional_false)) omitted = float(poisson.sf(tmax, (1 + r) * b)) return {"probability_lower": total, "probability_upper": min(1.0, total + omitted + 2e-13), "omitted_total_count_mass": omitted, "largest_conditional_false_alarm": float(np.max(conditional_false))} def _tests(): tests = [] def check(name, condition, detail=None): if not condition: raise AssertionError(name) tests.append({"name": name, "passed": True, "detail": detail}) check("zero_count_known_zero_background_closed_form", abs(red_source_upper_known_background(0, 0) + math.log(.05)) < 2e-14) check("zero_count_known_background_subtraction_and_floor", abs(red_source_upper_known_background(0, 1) - (-math.log(.05) - 1)) < 2e-14 and red_source_upper_known_background(0, 4) == 0) check("zero_on_and_off_finite_baseline_closed_form", abs(red_source_upper_finite_baseline(0, 0, 20) + math.log(.025)) < 2e-14) check("zero_off_lower_bound", poisson_total_mean_lower(0, .025) == 0) check("zero_background_detection_is_one_count", count_threshold_known_background(0, 1e-12) == 1) for b in [0, .01, .2, 3, 30, 1000]: for a in [.05, 1e-3, 1e-9, 1e-12]: k = count_threshold_known_background(b, a) check(f"known_threshold_minimal_b{b}_a{a}", poisson.sf(k - 1, b) <= a * (1 + 1e-11) and (k == 1 or poisson.sf(k - 2, b) > a)) for m in [0, 1, 10, 100]: for r in [2, 20]: for a in [.05, 1e-6]: k = count_threshold_onoff(m, r, a) check(f"onoff_threshold_minimal_m{m}_r{r}_a{a}", exact_onoff_pvalue(k, m, r) <= a * (1 + 1e-11) and (k == 1 or exact_onoff_pvalue(k - 1, m, r) > a)) check("onoff_zero_off_closed_form", abs(exact_onoff_pvalue(4, 0, 20) - (1 / 21)**4) < 1e-17) check("onoff_no_counts_p_is_one", exact_onoff_pvalue(0, 0, 20) == 1) check("Bonferroni_union_bound_arbitrary_dependence", alpha_for_search(.01, 1000000) == 1e-8) # Total-mean upper inversion and lower inversion. for n in [0, 1, 10, 100]: u = poisson_total_mean_upper(n, .05) check(f"upper_inversion_n{n}", abs(poisson.cdf(n, u) - .05) < 2e-13) if n: lo = poisson_total_mean_lower(n, .025) check(f"lower_inversion_m{n}", abs(poisson.sf(n - 1, lo) - .025) < 2e-13) invalid = [ lambda: count_threshold_known_background(-1, .05), lambda: count_threshold_known_background(1, 0), lambda: exact_onoff_pvalue(1.2, 0, 2), lambda: exact_onoff_pvalue(1, -1, 2), lambda: exact_onoff_pvalue(1, 0, 0), lambda: red_source_upper_finite_baseline(1, 0, 2, split=1), lambda: poisson_total_mean_upper(True), lambda: alpha_for_search(.01, 0)] for i, fn in enumerate(invalid): try: fn() except ValueError: check(f"invalid_input_rejected_{i}", True) else: check(f"invalid_input_rejected_{i}", False) return tests def main(output_directory=None, trials=200000): out = pathlib.Path(output_directory) if output_directory else pathlib.Path(__file__).resolve().parent out.mkdir(parents=True, exist_ok=True) tests = _tests() rng = np.random.default_rng(SEED) alpha = .05 coverage_rows = [] for s in [0.0, .3, 3.0, 30.0]: for b in [0.0, .2, 3.0, 30.0]: for r in [2.0, 20.0]: n = rng.poisson(s + b, trials) m = rng.poisson(r * b, trials) un, inverse_n = np.unique(n, return_inverse=True) um, inverse_m = np.unique(m, return_inverse=True) known = np.maximum(0, _upper_vec(un, alpha) - b)[inverse_n] finite = np.maximum(0, _upper_vec(un, alpha / 2)[inverse_n] - (_lower_vec(um, alpha / 2) / r)[inverse_m]) exact_k = _known_exact_coverage(s, b, alpha) exact_f = _finite_exact_coverage_bounds(s, b, r, alpha) for name, values, truth_lo, truth_hi in [ ("known_background", known, exact_k, exact_k), ("finite_baseline_bonferroni", finite, exact_f["coverage_lower"], exact_f["coverage_upper"])]: successes = int(np.count_nonzero(values >= s - 1e-12)) ci95 = _clopper_pearson(successes, trials) ci999 = _clopper_pearson(successes, trials, .999) consistency = ci999[0] <= truth_hi + 1e-12 and ci999[1] >= truth_lo - 1e-12 ci_family = _clopper_pearson(successes, trials, 1 - MC_FAMILY_ALPHA / MC_COMPARISONS) family_consistency = ci_family[0] <= truth_hi + 1e-12 and ci_family[1] >= truth_lo - 1e-12 if truth_lo < 1 - alpha - 2e-12: raise AssertionError("Exact numerical coverage below nominal confidence") coverage_rows.append({ "method": name, "source_mean_on": s, "background_mean_on": b, "exposure_ratio_T_off_over_T_on": r, "alpha": alpha, "trials": trials, "covered": successes, "empirical_coverage": successes / trials, "MC_Clopper_Pearson_95_lower": ci95[0], "MC_Clopper_Pearson_95_upper": ci95[1], "MC_Clopper_Pearson_99p9_lower": ci999[0], "MC_Clopper_Pearson_99p9_upper": ci999[1], "analytic_coverage_lower": truth_lo, "analytic_coverage_upper": truth_hi, "MC_99p9_interval_overlaps_analytic_enclosure": consistency, "MC_family99_lower": ci_family[0], "MC_family99_upper": ci_family[1], "MC_family99_interval_overlaps_analytic_enclosure": family_consistency}) detections = [] for b in [0.0, .2, 3.0, 30.0]: for r in [2.0, 20.0]: for a in [.05, .001]: n = rng.poisson(b, trials) m = rng.poisson(r * b, trials) k = count_threshold_known_background(b, a) exact_k = float(poisson.sf(k - 1, b)) hit_k = int(np.count_nonzero(n >= k)) hit_o = int(np.count_nonzero(binom.sf(n - 1, n + m, 1 / (1 + r)) <= a)) exact_o = _conditional_exact_false_alarm_bounds(b, r, a) ci_k = _clopper_pearson(hit_k, trials, .999) ci_o = _clopper_pearson(hit_o, trials, .999) family_k = _clopper_pearson(hit_k, trials, 1 - MC_FAMILY_ALPHA / MC_COMPARISONS) family_o = _clopper_pearson(hit_o, trials, 1 - MC_FAMILY_ALPHA / MC_COMPARISONS) family_consistency_k = family_k[0] <= exact_k + 1e-12 <= family_k[1] + 1e-12 family_consistency_o = family_o[0] <= exact_o["probability_upper"] + 1e-12 and family_o[1] >= exact_o["probability_lower"] - 1e-12 if exact_o["probability_upper"] > a + 2e-12: raise AssertionError("Conditional null size exceeds alpha") detections.append({ "background_mean_on": b, "exposure_ratio_T_off_over_T_on": r, "alpha_per_test": a, "trials": trials, "known_b_threshold_count": k, "known_b_exact_false_alarm": exact_k, "known_b_empirical_false_alarm": hit_k / trials, "known_b_MC_99p9_interval": ci_k, "onoff_empirical_false_alarm": hit_o / trials, "onoff_MC_99p9_interval": ci_o, "onoff_analytic_false_alarm_bounds": exact_o, "known_b_MC_family99_interval": family_k, "onoff_MC_family99_interval": family_o, "known_b_MC_family99_consistent": family_consistency_k, "onoff_MC_family99_consistent": family_consistency_o}) examples = [] for n, m, r, b in [(0, 0, 20, 0), (0, 20, 20, 1), (1, 20, 20, 1), (5, 100, 20, 5), (10, 200, 20, 10), (50, 1000, 20, 50), (5, 10, 2, 5), (10, 20, 2, 10)]: examples.append({ "n_on": n, "m_off": m, "exposure_ratio_T_off_over_T_on": r, "known_b_on_for_comparison": b, "plug_in_m_off_over_ratio": m / r, "red_source_upper95_known_b_counts": red_source_upper_known_background(n, b), "red_source_upper95_finite_baseline_counts": red_source_upper_finite_baseline(n, m, r), "exact_onoff_pvalue": exact_onoff_pvalue(n, m, r), "known_b_detection_threshold_alpha0p001": count_threshold_known_background(b, .001), "onoff_detection_threshold_alpha0p001_at_m_off": count_threshold_onoff(m, r, .001)}) search = [] for count in [1, 1000, 1000000]: a = alpha_for_search(.01, count) for b in [0, .2, 3, 30]: k = count_threshold_known_background(b, a) search.append({"alpha_family": .01, "actual_number_of_tests_assumed": count, "alpha_per_test": a, "background_mean_on": b, "threshold_count": k, "actual_per_test_tail": float(poisson.sf(k - 1, b)), "Bonferroni_family_bound": min(1.0, count * float(poisson.sf(k - 1, b)))}) result = { "schema": SCHEMA, "status": "EXECUTED_IDEAL_POISSON_STATISTICAL_CHECKS_NOT_INSTRUMENT_VALIDATION", "work_start_utc": START_UTC, "execution_finished_utc": datetime.now(timezone.utc).isoformat(), "seed": SEED, "generator": "numpy.default_rng/PCG64", "trials_per_configuration": trials, "MC_audit_family": {"alpha_family": MC_FAMILY_ALPHA, "comparisons": MC_COMPARISONS, "per_interval_alpha": MC_FAMILY_ALPHA / MC_COMPARISONS, "construction": "Bonferroni exact Clopper-Pearson intervals; no independence between diagnostics assumed", "seed_not_changed": True, "initial_diagnostic_observation": "The first run stopped at known-b s=30,b=3,r=20 with191686/200000 covered, exact.9567073249868225; its individual99.9%interval excludes exact. This observation is retained; no seed reroll. Simultaneous intervals are implementation diagnostics, not the coverage proof."}, "versions": {"python": platform.python_version(), "numpy": np.__version__, "scipy": scipy.__version__}, "exposure_ratio_definition": EXPOSURE_RATIO_DEFINITION, "assumptions": [ "Registered counts in nonoverlapping on/off windows are independent Poisson variables, before any deadtime/censoring issue.", "Source expectation s>=0 is confined to on window; background stationary with common rate and known exposure ratio.", "Known-background upper bound requires externally known fixed b_on, not a finite-baseline plug-in estimate.", "Bonferroni alpha split is fixed before observing counts; total nominal alpha=.05, default equal split.", "Limits cover expected detected source counts in the selected window, not emitted photons, current, peak rate, or source class.", "Per-test alpha is not a search false-alarm probability. A family bound needs the actual number of evaluated tests; operational overlapping-window trial count is unspecified.", "Rate times assumed17ns deadtime is only a diagnostic flag in the parent analysis. Flagging does not restore ideal Poisson counts; no Bernoulli-bin detector model is imposed here.", "Applying limits only after data-dependent selection or nondetection can change conditional coverage; unconditional bounds do not automatically survive that conditioning.", "Unknown backgrounds, drift, overdispersion, overlapping on/off windows, optical calibration and detector artefacts require separate models." ], "closed_form_controls": { "zero_count_known_zero_background_upper95": -math.log(.05), "zero_on_zero_off_finite_baseline_equal_split_upper95": -math.log(.025), "zero_background_detection_threshold_for_any_alpha_between_0_and_1": 1, "onoff_zero_off_pvalue": "(1/(1+exposure_ratio))**n_on"}, "deterministic_tests": tests, "coverage_MC_rows": coverage_rows, "null_detection_MC_rows": detections, "count_examples": examples, "search_trial_examples": search, "coverage_summary": { "coverage_MC_rows": len(coverage_rows), "minimum_known_background_analytic_coverage": min(z["analytic_coverage_lower"] for z in coverage_rows if z["method"]=="known_background"), "minimum_finite_baseline_analytic_coverage": min(z["analytic_coverage_lower"] for z in coverage_rows if z["method"]=="finite_baseline_bonferroni"), "all_MC_99p9_intervals_overlap_analytic_enclosures": all(z["MC_99p9_interval_overlaps_analytic_enclosure"] for z in coverage_rows), "MC_family99_all_overlap_analytic": all(z["MC_family99_interval_overlaps_analytic_enclosure"] for z in coverage_rows) and all(z["known_b_MC_family99_consistent"] and z["onoff_MC_family99_consistent"] for z in detections), "individual_99p9_coverage_disagreements": [z for z in coverage_rows if not z["MC_99p9_interval_overlaps_analytic_enclosure"]], "MC_intervals_are_diagnostics_not_proof": True, "proof_is_in_README": True}, "no_new_atmospheric_or_detector_simulation": True, "input_source_script_sha256": hashlib.sha256(pathlib.Path(__file__).read_bytes()).hexdigest() } (out/"count_inference_results.json").write_text(json.dumps(result, indent=2, allow_nan=False)+"\n", encoding="utf-8") for name, data in [("coverage_checks.csv",coverage_rows),("count_examples.csv",examples),("search_trial_examples.csv",search)]: with (out/name).open("w",newline="",encoding="utf-8") as f: w=csv.DictWriter(f,fieldnames=list(data[0])) w.writeheader();w.writerows(data) print(json.dumps({"status":result["status"],"tests":len(tests),"coverage_summary":result["coverage_summary"], "output_directory":str(out),"zero_controls":result["closed_form_controls"]},indent=2)) return result if __name__ == "__main__": main()