Review index · Download raw source
reference/statistics/count_inference/count_inference.py
1"""Exact ideal-Poisson count inference; no deadtime correction or model fitting.23Run with bundled Python 3.12, or any Python with numpy and scipy.4exposure_ratio always means T_off/T_on; E[M_off]=exposure_ratio*b_on.5"""6from __future__ import annotations7import csv8import hashlib9import json10import math11import numbers12import pathlib13import platform14import sys15from datetime import datetime, timezone1617sys.dont_write_bytecode = True18_PRIVATE_SCIPY = pathlib.Path("D:/Photometer_research_20261022/research/cloud_transport/private_readers")19if _PRIVATE_SCIPY.is_dir():20 sys.path.insert(0, str(_PRIVATE_SCIPY))2122import numpy as np23import scipy24from scipy.stats import beta, binom, chi2, poisson2526SCHEMA = "IDEAL_POISSON_COUNT_INFERENCE_V1"27START_UTC = "2026-10-22T08:03:13Z"28SEED = 2026102229MC_FAMILY_ALPHA = .0130MC_COMPARISONS = 96 # 64 coverage + 32 null-detection probabilities31EXPOSURE_RATIO_DEFINITION = "T_off/T_on; E[M_off]=exposure_ratio*b_on"323334def _count(value, name="count"):35 if isinstance(value, (bool, np.bool_)) or not isinstance(value, numbers.Integral) or value < 0:36 raise ValueError(f"{name} must be a nonnegative integer")37 return int(value)383940def _positive(value, name):41 if isinstance(value, (bool, np.bool_)):42 raise ValueError(f"{name} must be finite and positive")43 value = float(value)44 if not math.isfinite(value) or value <= 0:45 raise ValueError(f"{name} must be finite and positive")46 return value474849def _alpha(value):50 value = _positive(value, "alpha")51 if value >= 1:52 raise ValueError("alpha must lie strictly between zero and one")53 return value545556def _background(value):57 if isinstance(value, (bool, np.bool_)):58 raise ValueError("b_on must be finite and nonnegative")59 value = float(value)60 if not math.isfinite(value) or value < 0:61 raise ValueError("b_on must be finite and nonnegative")62 return value636465def alpha_for_search(alpha_family, number_of_tests):66 """Bonferroni per-test alpha; arbitrary dependence is allowed.6768 number_of_tests is the actual number of hypotheses evaluated, including69 all channels/windows/templates whose union defines the false alarm event.70 """71 a = _alpha(alpha_family)72 m = _count(number_of_tests, "number_of_tests")73 if m == 0:74 raise ValueError("number_of_tests must be positive")75 return a / m767778def count_threshold_known_background(b_on, alpha_per_test):79 """Smallest k>=1 with P(Poisson(b_on)>=k)<=alpha_per_test."""80 b = _background(b_on)81 a = _alpha(alpha_per_test)82 loga = math.log(a)83 hi = max(1, int(math.ceil(b + 8 * math.sqrt(b + 1) + 32)))84 while poisson.logsf(hi - 1, b) > loga:85 hi *= 286 lo = 0 # k=0 has tail probability one, and cannot satisfy alpha<1.87 while hi - lo > 1:88 mid = (lo + hi) // 289 if poisson.logsf(mid - 1, b) <= loga:90 hi = mid91 else:92 lo = mid93 return hi949596def exact_onoff_pvalue(n_on, m_off, exposure_ratio):97 """One-sided conditional-binomial p value for no source.9899 N~Pois(s+b), M~Pois(r*b), independent. Under s=0:100 N | (N+M=t) ~ Binomial(t,1/(1+r)). No estimated b is plugged in.101 """102 n = _count(n_on, "n_on")103 m = _count(m_off, "m_off")104 r = _positive(exposure_ratio, "exposure_ratio")105 return float(binom.sf(n - 1, n + m, 1 / (1 + r)))106107108def count_threshold_onoff(m_off, exposure_ratio, alpha_per_test):109 """Minimum n for the exact conditional test at the observed off count.110111 This is an equivalent decision-rule convenience, not a conditional112 false-alarm guarantee with M fixed. Unconditional validity follows by113 conditioning on N+M, not by conditioning on M.114 """115 m = _count(m_off, "m_off")116 r = _positive(exposure_ratio, "exposure_ratio")117 a = _alpha(alpha_per_test)118 p = 1 / (1 + r)119 loga = math.log(a)120 hi = max(1, int(math.ceil(m / r + 8 * math.sqrt(m / r + 1) + 32)))121 while binom.logsf(hi - 1, hi + m, p) > loga:122 hi *= 2123 lo = 0124 while hi - lo > 1:125 mid = (lo + hi) // 2126 if binom.logsf(mid - 1, mid + m, p) <= loga:127 hi = mid128 else:129 lo = mid130 return hi131132133def poisson_total_mean_upper(n, alpha=0.05):134 """One-sided Garwood upper bound: P_mu(X<=n)=alpha at mu=U."""135 n = _count(n, "n")136 a = _alpha(alpha)137 return float(chi2.isf(a, 2 * (n + 1)) / 2)138139140def poisson_total_mean_lower(m, alpha=0.05):141 """One-sided lower bound, zero for m=0."""142 m = _count(m, "m")143 a = _alpha(alpha)144 return 0.0 if m == 0 else float(chi2.ppf(a, 2 * m) / 2)145146147def red_source_upper_known_background(n_on, b_on, alpha=0.05):148 """At-least-(1-alpha) frequentist upper limit on expected source counts.149150 b_on is a known fixed mean in this same on window. No source prior,151 plug-in background, conditioning on nondetection, or unified interval.152 """153 return max(0.0, poisson_total_mean_upper(n_on, alpha) - _background(b_on))154155156def red_source_upper_finite_baseline(n_on, m_off, exposure_ratio, alpha=0.05, split=0.5):157 """Conservative source upper bound by Bonferroni subtraction.158159 T_off/T_on=r; alpha_total=alpha*split, alpha_off=alpha*(1-split).160 split must be fixed before seeing these counts. The union-bound proof161 does not require independence of confidence-bound failure events.162 The Poisson on/off observation model itself assumes independent windows.163 """164 n = _count(n_on, "n_on")165 m = _count(m_off, "m_off")166 r = _positive(exposure_ratio, "exposure_ratio")167 a = _alpha(alpha)168 f = _alpha(split)169 return max(0.0, poisson_total_mean_upper(n, a * f) -170 poisson_total_mean_lower(m, a * (1 - f)) / r)171172173def _upper_vec(n, a):174 return chi2.isf(a, 2 * (np.asarray(n) + 1)) / 2175176177def _lower_vec(m, a):178 m = np.asarray(m)179 # Avoid chi-square df=0 evaluations; lower bound there is exactly zero.180 out = np.zeros(m.shape, dtype=float)181 use = m > 0182 out[use] = chi2.ppf(a, 2 * m[use]) / 2183 return out184185186def _clopper_pearson(successes, trials, confidence=0.95):187 a = 1 - confidence188 lo = 0.0 if successes == 0 else float(beta.ppf(a / 2, successes, trials - successes + 1))189 hi = 1.0 if successes == trials else float(beta.isf(a / 2, successes + 1, trials - successes))190 return [lo, hi]191192193def _known_exact_coverage(s, b, alpha):194 if s == 0:195 return 1.0196 target = s + b197 hi = max(1, int(math.ceil(target + 8 * math.sqrt(target + 1) + 32)))198 if poisson_total_mean_upper(0, alpha) >= target:199 return 1.0200 lo = 0201 while hi - lo > 1:202 mid = (lo + hi) // 2203 if poisson_total_mean_upper(mid, alpha) >= target:204 hi = mid205 else:206 lo = mid207 return float(poisson.sf(hi - 1, target))208209210def _finite_exact_coverage_bounds(s, b, r, alpha, split=0.5, omitted_mass_tolerance=1e-12):211 """Analytic finite sum with explicit omitted-Poisson-tail enclosure.212213 No Monte Carlo is used here. The enclosure is numerical, not interval214 arithmetic certified against every possible floating-point error.215 """216 nmax = int(poisson.isf(omitted_mass_tolerance / 4, s + b))217 mmax = int(poisson.isf(omitted_mass_tolerance / 4, r * b))218 ns = np.arange(nmax + 1)219 ms = np.arange(mmax + 1)220 un = _upper_vec(ns, alpha * split)221 lm = _lower_vec(ms, alpha * (1 - split)) / r222 limits = np.maximum(0, un[:, None] - lm[None, :])223 # Tiny tolerance only prevents loss of exact-boundary points to roundoff.224 covered = limits >= s - 1e-12225 mass = poisson.pmf(ns, s + b)[:, None] * poisson.pmf(ms, r * b)[None, :]226 inside = float(np.sum(mass[covered]))227 missing_n = float(poisson.sf(nmax, s + b))228 missing_m = float(poisson.sf(mmax, r * b))229 missing = missing_n + missing_m - missing_n * missing_m230 return {"coverage_lower": min(1.0, inside),231 "coverage_upper": min(1.0, inside + missing + 2e-13),232 "omitted_joint_mass": missing, "n_max": nmax, "m_max": mmax}233234235def _conditional_exact_false_alarm_bounds(b, r, alpha, omitted_mass_tolerance=1e-12):236 tmax = int(poisson.isf(omitted_mass_tolerance / 4, (1 + r) * b))237 ts = np.arange(tmax + 1)238 p = 1 / (1 + r)239 critical = binom.isf(alpha, ts, p) + 1240 conditional_false = binom.sf(critical - 1, ts, p)241 total = float(np.sum(poisson.pmf(ts, (1 + r) * b) * conditional_false))242 omitted = float(poisson.sf(tmax, (1 + r) * b))243 return {"probability_lower": total,244 "probability_upper": min(1.0, total + omitted + 2e-13),245 "omitted_total_count_mass": omitted,246 "largest_conditional_false_alarm": float(np.max(conditional_false))}247248249def _tests():250 tests = []251 def check(name, condition, detail=None):252 if not condition:253 raise AssertionError(name)254 tests.append({"name": name, "passed": True, "detail": detail})255 check("zero_count_known_zero_background_closed_form",256 abs(red_source_upper_known_background(0, 0) + math.log(.05)) < 2e-14)257 check("zero_count_known_background_subtraction_and_floor",258 abs(red_source_upper_known_background(0, 1) - (-math.log(.05) - 1)) < 2e-14259 and red_source_upper_known_background(0, 4) == 0)260 check("zero_on_and_off_finite_baseline_closed_form",261 abs(red_source_upper_finite_baseline(0, 0, 20) + math.log(.025)) < 2e-14)262 check("zero_off_lower_bound", poisson_total_mean_lower(0, .025) == 0)263 check("zero_background_detection_is_one_count",264 count_threshold_known_background(0, 1e-12) == 1)265 for b in [0, .01, .2, 3, 30, 1000]:266 for a in [.05, 1e-3, 1e-9, 1e-12]:267 k = count_threshold_known_background(b, a)268 check(f"known_threshold_minimal_b{b}_a{a}",269 poisson.sf(k - 1, b) <= a * (1 + 1e-11)270 and (k == 1 or poisson.sf(k - 2, b) > a))271 for m in [0, 1, 10, 100]:272 for r in [2, 20]:273 for a in [.05, 1e-6]:274 k = count_threshold_onoff(m, r, a)275 check(f"onoff_threshold_minimal_m{m}_r{r}_a{a}",276 exact_onoff_pvalue(k, m, r) <= a * (1 + 1e-11)277 and (k == 1 or exact_onoff_pvalue(k - 1, m, r) > a))278 check("onoff_zero_off_closed_form",279 abs(exact_onoff_pvalue(4, 0, 20) - (1 / 21)**4) < 1e-17)280 check("onoff_no_counts_p_is_one", exact_onoff_pvalue(0, 0, 20) == 1)281 check("Bonferroni_union_bound_arbitrary_dependence", alpha_for_search(.01, 1000000) == 1e-8)282 # Total-mean upper inversion and lower inversion.283 for n in [0, 1, 10, 100]:284 u = poisson_total_mean_upper(n, .05)285 check(f"upper_inversion_n{n}", abs(poisson.cdf(n, u) - .05) < 2e-13)286 if n:287 lo = poisson_total_mean_lower(n, .025)288 check(f"lower_inversion_m{n}", abs(poisson.sf(n - 1, lo) - .025) < 2e-13)289 invalid = [290 lambda: count_threshold_known_background(-1, .05),291 lambda: count_threshold_known_background(1, 0),292 lambda: exact_onoff_pvalue(1.2, 0, 2),293 lambda: exact_onoff_pvalue(1, -1, 2),294 lambda: exact_onoff_pvalue(1, 0, 0),295 lambda: red_source_upper_finite_baseline(1, 0, 2, split=1),296 lambda: poisson_total_mean_upper(True),297 lambda: alpha_for_search(.01, 0)]298 for i, fn in enumerate(invalid):299 try:300 fn()301 except ValueError:302 check(f"invalid_input_rejected_{i}", True)303 else:304 check(f"invalid_input_rejected_{i}", False)305 return tests306307308def main(output_directory=None, trials=200000):309 out = pathlib.Path(output_directory) if output_directory else pathlib.Path(__file__).resolve().parent310 out.mkdir(parents=True, exist_ok=True)311 tests = _tests()312 rng = np.random.default_rng(SEED)313 alpha = .05314 coverage_rows = []315 for s in [0.0, .3, 3.0, 30.0]:316 for b in [0.0, .2, 3.0, 30.0]:317 for r in [2.0, 20.0]:318 n = rng.poisson(s + b, trials)319 m = rng.poisson(r * b, trials)320 un, inverse_n = np.unique(n, return_inverse=True)321 um, inverse_m = np.unique(m, return_inverse=True)322 known = np.maximum(0, _upper_vec(un, alpha) - b)[inverse_n]323 finite = np.maximum(0, _upper_vec(un, alpha / 2)[inverse_n] -324 (_lower_vec(um, alpha / 2) / r)[inverse_m])325 exact_k = _known_exact_coverage(s, b, alpha)326 exact_f = _finite_exact_coverage_bounds(s, b, r, alpha)327 for name, values, truth_lo, truth_hi in [328 ("known_background", known, exact_k, exact_k),329 ("finite_baseline_bonferroni", finite, exact_f["coverage_lower"], exact_f["coverage_upper"])]:330 successes = int(np.count_nonzero(values >= s - 1e-12))331 ci95 = _clopper_pearson(successes, trials)332 ci999 = _clopper_pearson(successes, trials, .999)333 consistency = ci999[0] <= truth_hi + 1e-12 and ci999[1] >= truth_lo - 1e-12334 ci_family = _clopper_pearson(successes, trials, 1 - MC_FAMILY_ALPHA / MC_COMPARISONS)335 family_consistency = ci_family[0] <= truth_hi + 1e-12 and ci_family[1] >= truth_lo - 1e-12336 if truth_lo < 1 - alpha - 2e-12:337 raise AssertionError("Exact numerical coverage below nominal confidence")338 coverage_rows.append({339 "method": name, "source_mean_on": s, "background_mean_on": b,340 "exposure_ratio_T_off_over_T_on": r, "alpha": alpha,341 "trials": trials, "covered": successes, "empirical_coverage": successes / trials,342 "MC_Clopper_Pearson_95_lower": ci95[0], "MC_Clopper_Pearson_95_upper": ci95[1],343 "MC_Clopper_Pearson_99p9_lower": ci999[0], "MC_Clopper_Pearson_99p9_upper": ci999[1],344 "analytic_coverage_lower": truth_lo, "analytic_coverage_upper": truth_hi,345 "MC_99p9_interval_overlaps_analytic_enclosure": consistency,346 "MC_family99_lower": ci_family[0], "MC_family99_upper": ci_family[1],347 "MC_family99_interval_overlaps_analytic_enclosure": family_consistency})348 detections = []349 for b in [0.0, .2, 3.0, 30.0]:350 for r in [2.0, 20.0]:351 for a in [.05, .001]:352 n = rng.poisson(b, trials)353 m = rng.poisson(r * b, trials)354 k = count_threshold_known_background(b, a)355 exact_k = float(poisson.sf(k - 1, b))356 hit_k = int(np.count_nonzero(n >= k))357 hit_o = int(np.count_nonzero(binom.sf(n - 1, n + m, 1 / (1 + r)) <= a))358 exact_o = _conditional_exact_false_alarm_bounds(b, r, a)359 ci_k = _clopper_pearson(hit_k, trials, .999)360 ci_o = _clopper_pearson(hit_o, trials, .999)361 family_k = _clopper_pearson(hit_k, trials, 1 - MC_FAMILY_ALPHA / MC_COMPARISONS)362 family_o = _clopper_pearson(hit_o, trials, 1 - MC_FAMILY_ALPHA / MC_COMPARISONS)363 family_consistency_k = family_k[0] <= exact_k + 1e-12 <= family_k[1] + 1e-12364 family_consistency_o = family_o[0] <= exact_o["probability_upper"] + 1e-12 and family_o[1] >= exact_o["probability_lower"] - 1e-12365 if exact_o["probability_upper"] > a + 2e-12:366 raise AssertionError("Conditional null size exceeds alpha")367 detections.append({368 "background_mean_on": b, "exposure_ratio_T_off_over_T_on": r,369 "alpha_per_test": a, "trials": trials,370 "known_b_threshold_count": k, "known_b_exact_false_alarm": exact_k,371 "known_b_empirical_false_alarm": hit_k / trials,372 "known_b_MC_99p9_interval": ci_k,373 "onoff_empirical_false_alarm": hit_o / trials,374 "onoff_MC_99p9_interval": ci_o, "onoff_analytic_false_alarm_bounds": exact_o,375 "known_b_MC_family99_interval": family_k, "onoff_MC_family99_interval": family_o,376 "known_b_MC_family99_consistent": family_consistency_k,377 "onoff_MC_family99_consistent": family_consistency_o})378 examples = []379 for n, m, r, b in [(0, 0, 20, 0), (0, 20, 20, 1), (1, 20, 20, 1),380 (5, 100, 20, 5), (10, 200, 20, 10), (50, 1000, 20, 50),381 (5, 10, 2, 5), (10, 20, 2, 10)]:382 examples.append({383 "n_on": n, "m_off": m, "exposure_ratio_T_off_over_T_on": r,384 "known_b_on_for_comparison": b, "plug_in_m_off_over_ratio": m / r,385 "red_source_upper95_known_b_counts": red_source_upper_known_background(n, b),386 "red_source_upper95_finite_baseline_counts": red_source_upper_finite_baseline(n, m, r),387 "exact_onoff_pvalue": exact_onoff_pvalue(n, m, r),388 "known_b_detection_threshold_alpha0p001": count_threshold_known_background(b, .001),389 "onoff_detection_threshold_alpha0p001_at_m_off": count_threshold_onoff(m, r, .001)})390 search = []391 for count in [1, 1000, 1000000]:392 a = alpha_for_search(.01, count)393 for b in [0, .2, 3, 30]:394 k = count_threshold_known_background(b, a)395 search.append({"alpha_family": .01, "actual_number_of_tests_assumed": count,396 "alpha_per_test": a, "background_mean_on": b,397 "threshold_count": k,398 "actual_per_test_tail": float(poisson.sf(k - 1, b)),399 "Bonferroni_family_bound": min(1.0, count * float(poisson.sf(k - 1, b)))})400 result = {401 "schema": SCHEMA, "status": "EXECUTED_IDEAL_POISSON_STATISTICAL_CHECKS_NOT_INSTRUMENT_VALIDATION",402 "work_start_utc": START_UTC, "execution_finished_utc": datetime.now(timezone.utc).isoformat(),403 "seed": SEED, "generator": "numpy.default_rng/PCG64", "trials_per_configuration": trials,404 "MC_audit_family": {"alpha_family": MC_FAMILY_ALPHA, "comparisons": MC_COMPARISONS,405 "per_interval_alpha": MC_FAMILY_ALPHA / MC_COMPARISONS,406 "construction": "Bonferroni exact Clopper-Pearson intervals; no independence between diagnostics assumed",407 "seed_not_changed": True,408 "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."},409 "versions": {"python": platform.python_version(), "numpy": np.__version__, "scipy": scipy.__version__},410 "exposure_ratio_definition": EXPOSURE_RATIO_DEFINITION,411 "assumptions": [412 "Registered counts in nonoverlapping on/off windows are independent Poisson variables, before any deadtime/censoring issue.",413 "Source expectation s>=0 is confined to on window; background stationary with common rate and known exposure ratio.",414 "Known-background upper bound requires externally known fixed b_on, not a finite-baseline plug-in estimate.",415 "Bonferroni alpha split is fixed before observing counts; total nominal alpha=.05, default equal split.",416 "Limits cover expected detected source counts in the selected window, not emitted photons, current, peak rate, or source class.",417 "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.",418 "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.",419 "Applying limits only after data-dependent selection or nondetection can change conditional coverage; unconditional bounds do not automatically survive that conditioning.",420 "Unknown backgrounds, drift, overdispersion, overlapping on/off windows, optical calibration and detector artefacts require separate models."421 ],422 "closed_form_controls": {423 "zero_count_known_zero_background_upper95": -math.log(.05),424 "zero_on_zero_off_finite_baseline_equal_split_upper95": -math.log(.025),425 "zero_background_detection_threshold_for_any_alpha_between_0_and_1": 1,426 "onoff_zero_off_pvalue": "(1/(1+exposure_ratio))**n_on"},427 "deterministic_tests": tests, "coverage_MC_rows": coverage_rows,428 "null_detection_MC_rows": detections, "count_examples": examples,429 "search_trial_examples": search,430 "coverage_summary": {431 "coverage_MC_rows": len(coverage_rows),432 "minimum_known_background_analytic_coverage": min(z["analytic_coverage_lower"] for z in coverage_rows if z["method"]=="known_background"),433 "minimum_finite_baseline_analytic_coverage": min(z["analytic_coverage_lower"] for z in coverage_rows if z["method"]=="finite_baseline_bonferroni"),434 "all_MC_99p9_intervals_overlap_analytic_enclosures": all(z["MC_99p9_interval_overlaps_analytic_enclosure"] for z in coverage_rows),435 "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),436 "individual_99p9_coverage_disagreements": [z for z in coverage_rows if not z["MC_99p9_interval_overlaps_analytic_enclosure"]],437 "MC_intervals_are_diagnostics_not_proof": True,438 "proof_is_in_README": True},439 "no_new_atmospheric_or_detector_simulation": True,440 "input_source_script_sha256": hashlib.sha256(pathlib.Path(__file__).read_bytes()).hexdigest()441 }442 (out/"count_inference_results.json").write_text(json.dumps(result, indent=2, allow_nan=False)+"\n", encoding="utf-8")443 for name, data in [("coverage_checks.csv",coverage_rows),("count_examples.csv",examples),("search_trial_examples.csv",search)]:444 with (out/name).open("w",newline="",encoding="utf-8") as f:445 w=csv.DictWriter(f,fieldnames=list(data[0]))446 w.writeheader();w.writerows(data)447 print(json.dumps({"status":result["status"],"tests":len(tests),"coverage_summary":result["coverage_summary"],448 "output_directory":str(out),"zero_controls":result["closed_form_controls"]},indent=2))449 return result450451452if __name__ == "__main__":453 main()454