Bin-by-bin (BBB) statistical treatment: numerical fixes + modularization

[davidwalter2]

This PR groups three lines of work on the bin-by-bin statistical treatment, all preserving the existing CLI and giving identical results to the pre-change baseline on tests/make_tensor.py (6/6 binByBinStatType × binByBinStatMode combinations) to all reported digits.

Numerical fixes

1. --binByBinStatType gamma --binByBinStatMode full postfit covariance was wrong. The Newton iteration over nbeta uses tf.Variable.assign_sub, which is not a differentiable op. After the loop, self.nbeta was treated as a leaf constant by the outer tape, so the implicit dependence dβ_profile/dx was dropped from the postfit Hessian — covariance and POI σ came out too small. Fix: take one more Newton step at the converged value with autodiff enabled, restoring the implicit-function-theorem gradient. On tests/make_tensor.py the sig POI σ goes from 0.00835 (buggy) to 0.01437, matching all five other type×mode combinations.

2. New --minBBKstat option. Extends betamask with kstat = sumw² / sumw2 < threshold so degenerate per-process bins (e.g. mixed-sign-weight cancellations) get their β nuisance fixed at β0 instead of being profiled. Useful in full mode when individual processes have very low effective stats per bin and otherwise produce a non-PD postfit Hessian. Default 0.0 preserves existing behaviour.

3. lite mode no longer scales zero-variance processes by β. For processes with sumw2 == 0 (e.g. data-driven backgrounds whose yield is determined by a paramModel), β is passed through unchanged. The merged kstat in lite mode is computed from finite-variance contributions only (data-driven processes were previously inflating sumw_total without contributing to sumw2_total, over-tightening β). Each lite-mode profile formula was re-derived to handle nexp(β) = β·n_active + n_zero_var; when no zero-variance process is present each formula reduces exactly to the original closed form, so the regression is exact. On a W-mass fit with one data-driven process (Fake), POI σ in lite matches full to 5 sig figs (0.09753 vs 0.09752), where previously lite was inflated by ~1%.

Modularization

The ~700–900 BBB-related lines previously interleaved through fitter.py move into a new rabbit/bbstat/ package:

  rabbit/bbstat/  
    __init__.py                               
    bbstat.py     # BinByBinStat: Variables, masks, kstat, profile + apply,
                  # constraint NLL, toy randomization, Newton helper                        
    formulas.py   # solve_quad_eq             

Fitter constructs a single self.bbstat = BinByBinStat(...) and delegates all β-related operations to it. _compute_yields_with_beta shrinks from ~614 lines to a 10-line wrapper. The legacy attribute / method shims on Fitter (binByBinStat, beta, kstat, betamask, set_beta0, …) have been removed; the two external callers (workspace.py, bin/rabbit_fit.py) now read fitter.bbstat. directly.

Net diff in fitter.py: 3118 → 2337 lines. New module: 942 lines.

The selective-lite naming uses physics-neutral terms: proc_zero_var_mask / n_active / n_zero_var (rather than "data-driven" / "MC", which were specific to one use case).

Tests + CI

• New tests/test_bbstat.py (15 tests, runs as a script and under pytest): solve_quad_eq, default β0 per type, lbeta = 0 at β = β0, normal-additive / normal-multiplicative constraint values, betamask construction, --minBBKstat extension, proc_zero_var_mask detection, full-mode skip, invalid type validation.
• unit-tests matrix in .github/workflows/main.yml gains the new test.
• All fit-running jobs (fitting, sparse-fits, alternative-fits, regularization, bsm, likelihoodscans, abcd-fits) now declare needs: [setenv, unit-tests, …] so a regression in the bbstat module short-circuits the pipeline before the long fits.

[bendavid]

Does 3 only handle the zero variance case, or does it also "distribute" the beta impact across processes with different (but non-zero) relative statistical uncertainty?

[davidwalter2]

Does 3 only handle the zero variance case, or does it also "distribute" the beta impact across processes with different (but non-zero) relative statistical uncertainty?

It distributes the beta impact across all processes with finite statistical uncertainty relative to their contribution, not considering their statistical uncertainty. This is the same way as it was done before just leaving out the zero variance contribution.

Intuitively I would expect that distributing the beta impact weighted with the relative statistical uncertainty of each process is more precise but I think it needs further study to quantify that. The BB full treats this exact.