03 — Poisson detection, pulse morphology and red upper limits

This is an executable, lightweight review of the statistical experiment. It uses synthetic receiver-input counts, not measured source-depth labels. It does not rerun the original 81,000-pulse generator.

Full derivations, equation IDs S01–S21, units, source functions and selection limits are in methods/statistics.md. The original scripts are retained for comparison; this notebook calls a new implementation.

Required: NumPy and SciPy in the portable runtime. Optional Node checks are labelled NOT_RUN when Node is unavailable. No downloads, installs or large simulations occur here.

from pathlib import Path
import sys, json, shutil, time
start = Path.cwd().resolve()
ROOT = next((p for p in [start, *start.parents] if (p / "modules/common.py").is_file()), None)
if ROOT is None:
    raise RuntimeError("Open this notebook from within the verification package.")
for p in [ROOT / "runtime/deps", ROOT / "modules"]:
    if p.is_dir():
        sys.path.insert(0, str(p))
from common import table, H, C, PHOTON_ENERGY_337, PHOTON_ENERGY_777
import statistics_review as st
import numpy as np
import scipy
print("Package:", ROOT)
print("NumPy:", np.__version__, "SciPy:", scipy.__version__)
print("Source:", Path(st.__file__).relative_to(ROOT))

Package: D:\Photometer_research_20261022\verification_20261022_r1
NumPy: 2.3.5 SciPy: 1.18.1
Source: modules\statistics_review.py

Execute the bounded checks

This call checks count tails and intervals, generates one explicit-bin pulse, independently fits its likelihood, and rebuilds decisions and intervals from the 81 saved synthetic trial summaries. It also checks source-free and background-free profile boundaries.

A software check does not establish physical source depth, hardware sensitivity or useful orbital event rate.

tick = time.perf_counter()
audit = st.run_checks(write=True)
elapsed_s = time.perf_counter() - tick
assert audit["passed"]
print("Checks completed in", round(elapsed_s, 3), "s")
print("Output: checks/statistics_results.json")
print("Saved scenarios:", audit["confusions"]["scenarios"])
print("Recomputed metrics:", audit["confusions"]["metrics_recomputed"])
print("Optional JavaScript status:", audit["javascript"]["status"])

Checks completed in 3.927 s
Output: checks/statistics_results.json
Saved scenarios: 81
Recomputed metrics: 648
Optional JavaScript status: EXECUTED_OPTIONAL_NODE

Photon reference planes — S01–S02

For an assumed efficiency $\eta$ and entrance area $A$, expected primary counts are $s=\eta A\Phi_{\rm ap}$. The aperture fluence already includes atmospheric delivery. Multiplication gain is not another photon efficiency.

The following is an assumed conversion, not a threshold or emitted-source yield.

table([audit["photon_planes"]], max_rows=1)

Detection with known background — S05–S06

The critical integer is the smallest $k$ such that $P_{\operatorname{Pois}(b)}(N\ge k)\le\alpha$. Required signal is an expected source count solving $P_{\operatorname{Pois}(b+s)}(N\ge k)=0.95$.

The two alphas are a single-gate tail and a hypothetical fixed search allocation. Neither validates an adaptive spacecraft trigger or a physical false-event rate.

table(audit["counts"]["threshold_power_checks"], max_rows=20)
print("Zero-background critical count:", st.critical_known(0))
print("Zero-background source mean for 95% detection:", st.required_signal_known(0)["source_mean_PE"])

Zero-background critical count: 1
Zero-background source mean for 95% detection: 2.995732273553991

Downstream threshold bindings

All 216 FOV and 192 auroral rows are checked against their declared known-background mean and false-alarm allocation. The FOV variable named alpha is a population slope, not the false-alarm alpha. The latter is recovered from the retained source constants without executing that source.

This verifies only the statistical columns. It does not validate a brightness population, sky scene, field geometry or useful mission rate.

bindings = audit["downstream_threshold_bindings"]
print(bindings["status"])
print("Rows:", bindings["rows_checked"], "Unique (background, alpha) computations:",
      bindings["unique_background_alpha_computations"])
print("Largest S95 difference (PE):", bindings["maximum_S95_absolute_difference_PE"])
print("Largest power difference:", bindings["maximum_power_absolute_difference"])
show = [r for r in bindings["rows"] if r["study"] == "fov_population"][:2]
show += [r for r in bindings["rows"] if r["study"] == "auroral_counts"][:2]
table(show, max_rows=4)

ALL408_STATISTICAL_COLUMN_BINDINGS_CHECKED
Rows: 408 Unique (background, alpha) computations: 228
Largest S95 difference (PE): 1.2287055142223835e-06
Largest power difference: 1.2996648202090455e-10

Finite baseline and red limits — S03–S12

For $N\sim\operatorname{Pois}(s+b)$, $M\sim\operatorname{Pois}(rb)$, $r=T_{\rm off}/T_{\rm on}$, condition on N+M under the null:

N\mid N+M \sim \operatorname{Binom}(N+M,1/(1+r)).

Known-background upper limits subtract a known fixed $b$. Finite-baseline limits use a conservative union-bound construction with a prechosen alpha split. Plugging $M/r$ into a known-background limit is a different procedure.

Selecting only nondetections using the same counts can alter conditional coverage.

table(audit["counts"]["onoff_checks"], max_rows=10)
table(audit["counts"]["upper_inversions"], max_rows=10)
print(json.dumps(audit["counts"]["zero_count_controls"], indent=2))
table(audit["counts"]["coverage_finite_enumerations"], max_rows=10)

{
  "known_zero_background95": 2.9957322735539913,
  "finite_baseline_zero_zero95": 3.6888794541139363
}

Worked likelihood with every bin exposed — S13–S16

The unit-peak ramp has an 80 µs 10–90% rise and ends 3 ms after onset. Early bins, one late plateau and an independent 10 ms prewindow form this example.

For each bin, $\lambda_j=Aq_j+\beta d_j$. Both $d_j$ and $q_j$ are in seconds; $A,\beta$ are primary count rates; $n_j,\lambda_j$ are counts.

With $Q=\sum q_j$, $D=\sum d_j$, $h=q/Q$, $w=d/D$, the total-scale MLE is $\widehat\eta_{\rm tot}=\sum n_j$. A concave profile gives $\widehat A=\widehat\theta N/Q$, $\widehat\beta=(1-\widehat\theta)N/D$.

A separate two-rate optimizer checks the original likelihood. It does not call the old R4 fitter. The table includes all bins in early/late/prewindow order, not chronological order.

example = audit["worked_example"]
print("Seed:", example["seed"], "Observed total:", example["observed_total_PE"])
table([
    {"method": "concave mixture profile", **example["profile"]},
    {"method": "direct two-rate SciPy optimization", **example["direct_SciPy_two_parameter_fit"]},
], max_rows=2)
print("Absolute log-likelihood difference:", example["absolute_loglikelihood_difference"])
print("Fitted total minus observed counts:", example["fitted_total_minus_observed_PE"])
table(example["bins"], max_rows=40)

Seed: 202610221110 Observed total: 743

Absolute log-likelihood difference: 2.1316282072803006e-14
Fitted total minus observed counts: 0.0

The next table searches three rise templates and an onset grid with unknown amplitude/background. Its winning template is a waveform morphology, not cloud depth. An individual noisy outcome may differ from its generating label.

table(example["top_grid_candidates"], max_rows=9)

Rebuild denominators and Monte Carlo intervals — S17–S19

A fixed profile-score gap $\ln19$ is a rejection heuristic, not calibrated 95% confidence. Detection/all trials, acceptance/all trials and wrong/accepted trials have different denominators.

The summaries contain fitted scores and gate/prewindow counts, not every original bin count. Reconstructing them verifies decisions, tables and intervals; it does not independently refit 81,000 pulses.

The table shows the 50 kS/s low-background subset. Full results remain in the machine-readable output. Error is undefined when no trial is accepted.

rebuilt = audit["confusions"]
print("Scenarios:", rebuilt["scenarios"], "Synthetic pulses:", rebuilt["signal_trials"])
print("Metrics:", rebuilt["metrics_recomputed"])
print("Largest interval difference from retained report:", rebuilt["max_interval_absolute_difference"])
rows = []
for x in rebuilt["rows"]:
    if x["sample_khz"] == 50 and x["background_cps"] == 50:
        detected, a3, w3, a2, w2, r3, r2, allerr = x["rebuilt_metrics"]
        rows.append({
            "plateau_equivalent_PE": x["N_plateau_equivalent_PE"],
            "true_receiver_rise_us": x["truth_rise_us"],
            "detected_n": detected["k"],
            "accepted_n": a2["k"],
            "wrong_n": w2["k"],
            "wrong_denominator": w2["n"],
            "wrong_fraction": w2["fraction"],
            "individual_exact95_interval": w2["Clopper_Pearson"],
            "confusion_fast_slow_reject_undetected": x["confusion_fastslow"],
        })
table(rows, max_rows=12)

Scenarios: 81 Synthetic pulses: 81000
Metrics: 648
Largest interval difference from retained report: 2.220446049250313e-16

Independent numerical precision checks

Selected tails are summed from Poisson masses; selected on/off tails are finite binomial sums. No Gaussian approximation replaces the count probabilities.

The high-mean pins come from a retained 60-digit Decimal/Stirling/geometric-remainder calculation. They expose a small SciPy incomplete-gamma precision limitation near a mean of $10^6$. The tested critical integers agree. This is bounded floating-point QA, not formal directed interval arithmetic.

table(audit["counts"]["pinned_high_mean_controls"], max_rows=4)
table(audit["counts"]["tail_checks"], max_rows=12)

Known-parameter oracle — S20–S21

For known independent Poisson means, ${\rm BC}=\exp[-\tfrac12\sum_j(\sqrt{\mu_j}-\sqrt{\nu_j})^2]$. Under equal priors and zero-one loss:

(1-\sqrt{1-{\rm BC}^2})/2\le P_e^\star\le{\rm BC}/2.

Three new small two-bin joint-PMF controls test the formula and bounds. The original 320-row grid is not independently rebuilt here. Unknown nuisance fits and accepted-only practical errors have different assumptions and loss/selection rules.

table(audit["oracle_affinity"]["rows"], max_rows=3)
print("Original 320-row grid:", audit["oracle_affinity"]["retained320_grid"])

Original 320-row grid: NOT_INDEPENDENTLY_REBUILT_IN_THIS_LIGHT_REVIEW

Source-function audit and optional browser execution

Retaining original code is provenance; running it twice is not independent theory validation. The likelihood is derived in the methods chapter, and this new profile is checked against direct two-parameter optimization.

The browser copy is linked by hash/function lines. Optional actual Node checks exercise thresholds, red limits and the corrected tiny-background boundary. NOT_RUN means JavaScript was not executed in this environment.

for row in audit["source_map"]:
    print(json.dumps(row, indent=2))
print(json.dumps({k:v for k,v in audit["javascript"].items() if k != "checks"}, indent=2))
table(audit["javascript"].get("checks", []), max_rows=12)
equations = []
for eq, funcs in audit["equation_function_map"].items():
    for f in funcs:
        equations.append({"equation": eq, **f})
table(equations, max_rows=40)

{
  "file": "reference\\statistics\\oracle_reference\\calculate_ideal_poisson_rise.py",
  "present": true,
  "sha256": "542d4cc4495b0e7d379b3d6071bc668d904533567638ccca35853ac56abc3269",
  "functions": {
    "experiment": 18,
    "affinity_bounds": 28
  }
}
{
  "file": "reference\\statistics\\practical_morphology.py",
  "present": true,
  "sha256": "08df8f61294885ba658d75ae23bc3a7bb6511ffd01f9608a0de5a143f33f477c",
  "functions": {
    "primitive": 25,
    "fit_one": 57,
    "draw_case": 138,
    "analyze_case": 154
  }
}
{
  "file": "reference\\statistics\\count_inference\\count_inference.py",
  "present": true,
  "sha256": "4e4e25ae4db608500fb51c8406b599a7ee8efb92a3adee0f49e8a6d2a45eba9c",
  "functions": {
    "count_threshold_known_background": 78,
    "exact_onoff_pvalue": 96,
    "poisson_total_mean_upper": 133,
    "poisson_total_mean_lower": 140,
    "red_source_upper_known_background": 147,
    "red_source_upper_finite_baseline": 156
  }
}
{
  "file": "reference\\statistics\\browser_math.js",
  "sha256": "69186918327b076b053769ad3597c3f472a93919e5cd73057303f8034704bd37",
  "function_lines": {
    "sf": 5,
    "critical": 7,
    "threshold": 8,
    "upper": 9
  },
  "status": "EXECUTED_OPTIONAL_NODE",
  "node_version": "v24.19.0",
  "red_upper_checks": [
    {
      "n_on": 0,
      "b_on_PE": 0,
      "JS_upper95": 2.995732273553991,
      "Python_upper95": 2.9957322735539913,
      "abs_difference": 4.440892098500626e-16
    },
    {
      "n_on": 0,
      "b_on_PE": 2,
      "JS_upper95": 0.9957322735539909,
      "Python_upper95": 0.9957322735539913,
      "abs_difference": 4.440892098500626e-16
    },
    {
      "n_on": 1,
      "b_on_PE": 0,
      "JS_upper95": 4.743864518390579,
      "Python_upper95": 4.743864518390579,
      "abs_difference": 0.0
    },
    {
      "n_on": 3,
      "b_on_PE": 2,
      "JS_upper95": 5.753656527932727,
      "Python_upper95": 5.753656527932727,
      "abs_difference": 0.0
    },
    {
      "n_on": 10,
      "b_on_PE": 5,
      "JS_upper95": 11.962219235721907,
      "Python_upper95": 11.962219235721896,
      "abs_difference": 1.0658141036401503e-14
    },
    {
      "n_on": 100,
      "b_on_PE": 80,
      "JS_upper95": 38.07927278209726,
      "Python_upper95": 38.079272782097064,
      "abs_difference": 1.9895196601282805e-13
    }
  ],
  "tiny_background_boundary_checks": [
    {
      "b_PE": 2.8665161009724685e-07,
      "alpha": 2.866515718791933e-07,
      "JS": {
        "sf1": 2.86651569012678e-07,
        "critical": 1
      },
      "expected_sf1": 2.86651569012678e-07,
      "expected_critical": 1
    },
    {
      "b_PE": 2.866516158302791e-07,
      "alpha": 2.866515718791933e-07,
      "JS": {
        "sf1": 2.8665157474570857e-07,
        "critical": 2
      },
      "expected_sf1": 2.8665157474570857e-07,
      "expected_critical": 2
    },
    {
      "b_PE": 7.716049305585323e-12,
      "alpha": 7.716049382716049e-12,
      "JS": {
        "sf1": 7.716049305555554e-12,
        "critical": 1
      },
      "expected_sf1": 7.716049305555554e-12,
      "expected_critical": 1
    },
    {
      "b_PE": 7.716049459906311e-12,
      "alpha": 7.716049382716049e-12,
      "JS": {
        "sf1": 7.716049459876542e-12,
        "critical": 2
      },
      "expected_sf1": 7.716049459876542e-12,
      "expected_critical": 2
    }
  ]
}

Interpretation limits

• kS/s specifies integration spacing, not analog noise bandwidth.
• Analog post-counting filtering generally correlates samples; it is absent from the independent-Poisson shape model.
• Deadtime flags do not restore missing photons or Poisson statistics.
• Red limits require calibration/live time and a separately justified thermal-alternative criterion.
• Pulse labels describe expected received rate after transport, not intrinsic histories or known depths.
• Physical populations, exposure, model mismatch and flight hardware remain unvalidated. Useful mission rate is not computable from these inputs.

Expected controls: 2.995732 PE known-zero-background upper limit; 3.688879 PE specified finite-baseline zero/zero bound; 81 scenarios, 81,000 saved pulses and 648 rebuilt metrics. Worked log-likelihood difference tolerance is $2\times10^{-7}$; interval differences must be below $2\times10^{-12}$.