# Statistical verification: detection, pulse shape and a red upper limit

This chapter specifies the probability experiment before giving numerical checks. The measurements below are synthetic receiver-input photon counts. A “30 μs”, “80 μs” or “200 μs” label describes the expected received waveform. It is not an independently known cloud depth, discharge current or physical source class.

The notebook **03_statistics.ipynb** executes a new worked likelihood fit and rebuilds the reported decisions from all 81 saved trial summaries. It does not rerun the original 81,000-pulse generator by default. The original scripts are retained under reference/statistics for comparison; the verification module does not import their fitting or interval functions.

## Reference planes and units

The optical entrance is the reference plane after atmospheric propagation. Let Φ_ap be event photon fluence at that plane, in photons m⁻², and A_ent the effective entrance area in m². For a specified spectral, angular and temporal response η, the expected primary detected count is

\[
E_\gamma(\lambda)=hc/\lambda,\qquad
s=\eta A_{\rm ent}\Phi_{\rm ap}. \tag{S01}
\]

Here s is an expectation in photoelectrons or registered primary photon events, before PMT/APD charge multiplication. It is not a realized count and does not include background. For a real broadband/extended scene, replace the scalar η by the appropriate spectral-angular-time integral. Internal charge gain is not another photon-detection efficiency.

The reverse conversion is

\[
N_{\rm ap}=s/\eta,\qquad
\Phi_{\rm ap}=s/(\eta A_{\rm ent}),\qquad
F_{E,{\rm ap}}=E_\gamma\Phi_{\rm ap}. \tag{S02}
\]

These equations do not infer total emitted photons. A spacecraft-plane observed fluence has already undergone cloud and atmospheric transport; do not attenuate it again. The module imports h, c and the two photon energies from modules/common.py.

## State the on/off experiment

In this section b is the expected background **count in the on window**, not a count rate. Let s≥0 be the expected source count in that window and r=T_off/T_on>0 the known relative exposure. A stationary background and disjoint, uncontaminated windows give

\[
N\sim{\rm Pois}(s+b),\quad
M\sim{\rm Pois}(rb),\quad N\perp M. \tag{S03}
\]

The joint likelihood, with integer observations n and m, is

\[
{\cal L}(s,b)=
\frac{e^{-(s+b)}(s+b)^n}{n!}
\frac{e^{-rb}(rb)^m}{m!}. \tag{S04}
\]

This model excludes deadtime, saturation, afterpulses, correlated avalanches, background drift and source contamination of the off window. An effective exposure ratio can include differing throughput or live time only when that correction is justified. A finite prewindow is not an externally known background.

### Known-background detection and power

If b is known exactly and fixed before the test, choose the smallest integer k≥1 satisfying

\[
P_b(N\ge k)\le\alpha. \tag{S05}
\]

Declare a positive excess when N≥k. The count discreteness generally makes the actual size smaller than α. At b=0, k=1 for every 0<α<1. This zero-background control is mathematical, not a claim that a spacecraft detector has no false counts.

For a requested detection probability p_det, solve for the expected source count:

\[
P_{b+s_{p_{\rm det}}}(N\ge k)=p_{\rm det}. \tag{S06}
\]

The notebook uses p_det=0.95. It verifies the returned integer, the preceding integer and the achieved power. This is a specified-gate signal-count requirement, not a mission completeness curve.

### Unknown background: exact conditional on/off test

Under the null s=0, condition on the observed total t=N+M:

\[
N\mid(N+M=t)\sim{\rm Binom}\left(t,\frac{1}{1+r}\right),\quad
p_{\rm on/off}=P(X\ge n\mid t). \tag{S07}
\]

For every fixed t the rejection probability is at most α; averaging over t retains that bound for every constant unknown b. The condition is on **N+M**, not on M alone. A displayed minimum on-count at a fixed observed off-count is a convenient boundary representation, not a conditional false-alarm guarantee given M.

The review independently sums the finite binomial tail for selected examples and compares it with SciPy's survival function. It does not insert M/r into a known-background test.

### Search trials

If K is a predetermined number of tested hypotheses, or a deterministic upper bound, the union bound and linearity of expectation give

\[
P(\hbox{one or more false candidates})\le K\alpha,\qquad
E[N_{\rm false\ candidates}]\le K\alpha. \tag{S08}
\]

Independence between overlapping tests is unnecessary for this bound. A retrospectively observed stopping count is not automatically a valid K. The illustrative allocation α=0.1/(50,000×86,400×3) assumes that particular number of evaluated starts and window families. It is neither an empirically measured event rate nor a complete model of an adaptive spacecraft trigger, clustering or particle bursts.

## A red nondetection is an upper confidence limit

For a Poisson total mean μ, the one-sided Garwood upper limit is

\[
U(n;\alpha)=\tfrac12\chi^2_{2(n+1),\,1-\alpha},
\quad P_{\mu=U}(X\le n)=\alpha. \tag{S09}
\]

Its repeated-sampling coverage is at least 1−α; discreteness can make it conservative. With externally known b,

\[
s_{\rm up}=\max[0,U(n;\alpha)-b]. \tag{S10}
\]

For s>0, this covers s whenever U(N;α) covers s+b. At s=0 the floor always covers. A downward fluctuation relative to a known b can produce an upper bound of zero; it is not a posterior assertion that a source is impossible.

For a finite independent off window, choose α_t+α_b=α **before examining the data**. Define the lower bound

\[
L(m;\alpha_b)=
\begin{cases}
0,&m=0,\\
\tfrac12\chi^2_{2m,\,\alpha_b},&m>0 .
\end{cases} \tag{S11}
\]

Then

\[
s_{\rm up,off}=
\max\left[0,U(n;\alpha_t)-L(m;\alpha_b)/r\right]
\quad\hbox{covers with probability at least }1-\alpha_t-\alpha_b. \tag{S12}
\]

On the event that both component bounds cover, U≥s+b and L/r≤b, so their difference is at least s. The probability either component fails is bounded by α_t+α_b. Independence of the two bound-failure events is unnecessary for this union-bound proof; the observation model in S03 still specifies independent windows.

At n=0 and known b=0, the 95% upper mean is −ln(0.05)=2.995732 counts. At n=m=0, the equal-split conservative finite-baseline construction gives −ln(0.025)=3.688879 counts. These are different stated procedures, not conflicting physical measurements. The latter is a conservative bound, not an “exact joint 95%” interval or a Bayesian probability.

Unconditional coverage need not remain nominal after selecting only nondetections, accepted shapes or extreme ratios using the same data. An independent trigger can avoid some selection effects, but its independence must be demonstrated. A red upper limit is on the expected registered source count in the specified window. It requires calibrated response and live time to become an aperture-photon bound, and a separate physical lower red/UV bound to exclude a thermal alternative. No universal positive thermal 777/337 bound is supplied by these statistics.

## Unknown amplitude, onset and background

For the waveform calculation, use β for a background **rate** in counts s⁻¹. For hardware bin j of duration d_j in seconds, let q_j=∫_bin f(t−t₀;R)dt, also in seconds, where f is a unit-peak pulse template. A is the peak source count rate in counts s⁻¹. The bin expectation is

\[
\lambda_j=Aq_j+\beta d_j,\qquad
n_j\sim{\rm Pois}(\lambda_j)\quad\hbox{independently}. \tag{S13}
\]

The background prewindow is another bin with q_pre=0. The unit-peak template rises linearly over D_ramp=R_10–90/0.8, stays on a plateau, and ends 3 ms after onset. Exact primitive differences give q_j; a midpoint sample is not substituted for the integral.

For each candidate onset and rise, define Q=Σq_j, D_tot=Σd_j, h_j=q_j/Q and w_j=d_j/D_tot. With η_tot=AQ+βD_tot and θ=AQ/η_tot,

\[
\ell=N\log\eta_{\rm tot}-\eta_{\rm tot}
+\sum_j n_j\log[(1-\theta)w_j+\theta h_j]
-\sum_j\log(n_j!),\quad N=\sum_j n_j. \tag{S14}
\]

Both h and w sum to one. Therefore, for N>0, the maximum in total scale is η_hat=N. The remaining one-dimensional objective is concave:

\[
g'(\theta)=\sum_j\frac{n_j(h_j-w_j)}
{(1-\theta)w_j+\theta h_j},\qquad
g''(\theta)=-\sum_j\frac{n_j(h_j-w_j)^2}
{[(1-\theta)w_j+\theta h_j]^2}\le0 . \tag{S15}
\]

Check both endpoints and solve the interior derivative root where appropriate. At N=0 the mixture is unidentifiable and A=β=0 suffices. The fitted rates are

\[
\widehat A=\widehat\theta N/Q,\qquad
\widehat\beta=(1-\widehat\theta)N/D_{\rm tot}. \tag{S16}
\]

The new worked example saves every n_j, d_j, q_j and fitted λ_j. It compares this profile solution with a separate two-parameter SciPy optimization of the original A q_j+β d_j likelihood. The original R4 grouped/vectorized implementation is not imported. Agreement checks implementation of the algebra; it cannot validate the physical pulse family.

The practical study searches a declared onset grid using early rise bins and a late plateau, with a 10 ms source-free prewindow. It excludes the final pulse cutoff from the shape fit. A coarse onset bracket is provided externally; blind search and its false-candidate rate are not established.

### Detection, acceptance and shape error have different denominators

The rule first detects a fixed candidate gate. It then selects a shape and rejects insufficiently separated profile scores:

\[
{\rm accept}={\rm detected}\ \land\
[\ell_{\rm winning\ profile}-\ell_{\rm competing\ profile}\ge\ln19]. \tag{S17}
\]

This margin was fixed in advance. It is not calibrated 95% confidence, a 95% correctness guarantee or a posterior probability. Maximizing over unknown parameters and unequal template families does not supply a normalized likelihood under a specified population prior.

The saved trials support three separate quantities:

- detected / all trials;
- accepted / all trials;
- wrong / accepted trials.

Always read error among accepted cases together with acceptance and rejection. An error fraction based on zero accepted cases is undefined, not zero. A sparse accepted subset has a wide uncertainty interval even if it contains no observed mistakes.

The review rebuilds all decisions, confusion columns and denominators from saved synthetic scores/count summaries. It recomputes 648 metric intervals across 81 scenarios. Those summaries do not retain every original bin count; rebuilding the confusion table is not an independent re-fit of all 81,000 pulses. The new explicit-bin worked example supplies a separate likelihood check.

### Monte Carlo intervals

For k errors out of a fixed n relevant trials, the Wilson score interval uses z=Φ⁻¹(1−α/2) and p_hat=k/n:

\[
{\rm Wilson}=
\frac{\widehat p+z^2/(2n)\ \pm\
z\sqrt{\widehat p(1-\widehat p)/n+z^2/(4n^2)}}{1+z^2/n}. \tag{S18}
\]

The Clopper–Pearson endpoints are

\[
p_L=B^{-1}_{k,n-k+1}(\alpha/2),\qquad
p_U=B^{-1}_{k+1,n-k}(1-\alpha/2), \tag{S19}
\]

with p_L=0 at k=0 and p_U=1 at k=n. This is an individual binomial interval; no simultaneous family coverage is claimed for displaying many such intervals. The intervals describe repeated synthetic trials under this specified generator, not uncertainty in real cloud populations or correctness of the source model. The review recomputes both kinds and compares them with the published values.

## Numerical accuracy and the browser cross-reference

The review compares SciPy tails with a positive Poisson-PMF recurrence, retaining a decreasing geometric bound for truncated terms. Small-count on/off tails are also summed directly. No normal/Gaussian approximation replaces the count probabilities.

“Exact Poisson” identifies the discrete statistical model and tail construction; ordinary floating-point code is not formal interval arithmetic. Near a mean of 10⁶ and a roughly 5σ upper tail, the retained independent high-precision audit found SciPy's incomplete-gamma survival value differed by about 4.6 parts per million from direct summation. The tested integer thresholds remain the same. The notebook reports this reference-precision qualification rather than hiding it.

The copied reference/statistics/browser_math.js provides the live web-engine functions sf, critical, threshold and upper. The source map gives their line numbers and hash. If Node is available, an optional check actually executes selected browser thresholds and compares them with Python. If Node is unavailable, the result is **NOT_RUN**, not a JavaScript validation claim. Node is not required for the Python notebook.

The previous tiny-background cancellation defect was corrected by computing P(N≥1)=−expm1(−b), and zero-count mass/CDF as exp(−b). The reference copied into this package is the corrected version. The browser has an explicitly bounded live domain; use the offline code beyond its stated limits.

## Known-hypothesis oracle: what the ideal bounds mean

For independent Poisson bins with two fully specified mean vectors μ and ν, the Bhattacharyya affinity is

\[
{\rm BC}(P_\mu,P_\nu)
=\prod_j\sum_{n=0}^{\infty}\sqrt{p_{\mu_j}(n)p_{\nu_j}(n)}
=\exp\left[-\frac12\sum_j(\sqrt{\mu_j}-\sqrt{\nu_j})^2\right]. \tag{S20}
\]

For one bin, the sum is the exponential series for \(\exp(\sqrt{\mu\nu})\), multiplied by \(\exp[-(\mu+\nu)/2]\). Independence multiplies affinities. No Gaussian approximation is involved.

Under equal prior probabilities and zero-one classification loss, the optimal binary error satisfies

\[
\frac{1-\sqrt{1-{\rm BC}^2}}{2}
\le P_e^\star=\frac{1-{\rm TV}(P,Q)}{2}
\le \frac{{\rm BC}}{2}. \tag{S21}
\]

The upper bound follows from \(\min(P,Q)\le\sqrt{PQ}\). For the lower bound, write \(|P-Q|=|\sqrt P-\sqrt Q|(\sqrt P+\sqrt Q)\) and apply Cauchy–Schwarz; then TV≤sqrt(1−BC²).

The review independently enumerates three small two-bin joint PMFs and checks their affinity and exact finite-grid Bayes errors against these expressions. The finite grids contain probability within 1e-13 of one. **The original 320-row oracle grid is retained but not independently rebuilt by this light notebook.**

These bounds concern known simple hypotheses, equal priors and a specified loss. Substituting best-fit unknown amplitudes/onsets/backgrounds does not retain them. They also do not rank a practical accepted-only error fraction without accounting for rejected events and its different decision problem.

## Statistical bindings in the FOV and auroral tables

The downstream audit recomputes critical integers and source means for 95% detection for all 216 FOV rows and 192 auroral rows. It caches repeated background/false-alarm pairs, while comparing every stored row.

Auroral rows state their known background and false-alarm alpha directly. The FOV table's variable called alpha is the **hypothetical population slope**, not the false-alarm allocation. Its retained source uses count_requirements.ALPHA_OPERATIONAL=0.1/(50,000×86,400×3) and B=(sky count rate+50 dark counts/s)×window. These constants are read through the Python AST; the retained source is not executed.

This checks only the statistical-column binding to the declared background, allocation and power. It does not independently validate FOV geometry, sky radiance, luminosity law, orbital coverage, real counter response or mission event rate. Source hashes and the maximum numerical differences are recorded in the check output.



## What the verification does not establish

There is no analog post-counting low-pass in the Poisson shape generator. A deterministic mean RC illustration elsewhere does not make filtered samples independent Poisson variables. Deadtime occupancy flags do not restore missing arrivals or Poisson likelihoods.

A 10 kS/s stream has 100 μs bins; 50 kS/s gives 20 μs; 100 kS/s gives 10 μs. These are sampling/integration rates, not analog noise bandwidths. The input shapes refer to an expected received count rate after transport, not intrinsic nitrogen excitation histories.

No absolute source luminosity distribution, physical-depth labels, end-of-life hardware response, thermal alternative ratio, polar background climatology, operational false-event rate or mission useful-event rate is determined by these checks.

## Source and function crosswalk

Equation IDs S01–S02 map to photon_planes and common.py constants; S03–S12 to critical_known, required_signal_known, onoff_pvalue, source_upper_known and source_upper_off; S13–S16 to primitive, q_integrals, poisson_loglik, nuisance_profile and worked_example; S17–S19 to confusion_checks and mc_intervals; S20–S21 to affinity_bound_values and oracle_affinity_controls.

checks/statistics_results.json records exact source hashes and function lines for retained originals. In the original scripts:

| Original file | Function | Original line |
|---|---|---:|
| practical_morphology.py | fit_one | See generated source map |
| practical_morphology.py | draw_case | 138 |
| practical_morphology.py | analyze_case | 154 |
| count_inference/count_inference.py | count_threshold_known_background | 78 |
| count_inference/count_inference.py | exact_onoff_pvalue | 96 |
| count_inference/count_inference.py | poisson_total_mean_upper/lower | 133 / 140 |
| count_inference/count_inference.py | red_source_upper_known_background | 147 |
| count_inference/count_inference.py | red_source_upper_finite_baseline | 156 |

Line numbers are for the retained versions; the recorded hashes and generated map take precedence if formatting changes. Retaining old code is provenance, not independent theory validation.

Primary implementation documentation: [SciPy Poisson](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.poisson.html), [SciPy binomial](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom.html), [SciPy chi-square](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html), and [NumPy Poisson generator](https://numpy.org/doc/stable/reference/random/generated/numpy.random.Generator.poisson.html). The mathematical likelihood and coverage arguments are derived above, not inferred from software documentation.

## Execute

Run notebook 03 from top to bottom for the light review. To generate the machine-readable checks without a notebook:

    python modules/statistics_review.py --write

Optional browser execution, with a known Node executable:

    python modules/statistics_review.py --write --node /path/to/node

The default does not rerun the original large pulse simulation. It uses packaged data/statistics_results.json and the 81 own-generated data/trials/case_XXX.npz summaries. Output is confined to checks/statistics_results.json, checks/statistics_summary.json and checks/statistics_worked_bins.csv.

