Skip to content

mcstat: folded systematics and k-fold protections - #231

Draft
davidwalter2 wants to merge 6 commits into
WMass:mainfrom
davidwalter2:261110_mcStatFoldSyst
Draft

davidwalter2 wants to merge 6 commits into
WMass:mainfrom
davidwalter2:261110_mcStatFoldSyst

Conversation

@davidwalter2

Copy link
Copy Markdown
Collaborator

Follow-up to #141: folded (per-fold) systematics in the k-fold MC-stat de-biasing, plus optional protections against the runaways seen in the W-mass fit with split logk (WMass/WRemnants#731, #748).

TensorWriter (folded systematics)

  • add_systematic(..., fold_axis=...) with an up/down pair now accepts symmetrize="linear"/"quadratic" and splits every fold into SymAvg/SymDiff (the default for QCD scale, PDF, αs, quark-mass systematics in WRemnants). TensorWriter.fold_symmetrize_modes lists the supported modes so callers can fall back to a shared systematic. The sparse asymmetric case is refused before anything is registered.
  • In bins where a folded process has no fold structure (negative fold, or no MC-stat variance), every fold gets the full logk (1/k for the normal form). Previously the nominal's own fold noise entered the variation there.
  • Folds with fewer effective MC events (sumw²/sumw2) than mcstat_fold_min_kstat=3 take the full logk. In a W-mass test tensor such folds had logk up to 25, compared with about 0.05 elsewhere. The count is logged per process.
  • Empty folds (zero nominal) get logk 0, as for an empty full template. Before, they got log(1e-3), which gives 0·inf = NaN at large θ.

Fitter (all off by default)

  • --mcStatFoldBound κ (--mcStatFoldBoundForm taper|tanh): bounds the fold deviation x = E_f − E of the log_normal split-logk exponent by x/(1+x²/κ²). Value, gradient and curvature at the nominal point are unchanged. Far out the objective goes back to the ordinary likelihood. The jackknife combination is otherwise unbounded below.
  • --mcStatFoldMaxDeviation τ: (bin, process) templates whose fold deviation norm max_f |logk_f − logk|₂ exceeds τ take the shared logk in every fold. The k=2 noise-floor estimate over-subtracts on these templates.
  • --mcStatFoldMinKstat n: templates with kstat < n get the full-sample template in every fold (nominal, variance and logk), so they enter like the ordinary likelihood.

The W-mass fits use κ = 0.5, τ = 0.1 (and --mcStatFoldMinKstat 3 with the photon-induced samples). The design and the W-mass numbers are in docs/mcstat/RABBIT_MCSTAT_DESIGN.md. Tests are in tests/test_mcstat_kfold.py.

Draft: the history is rebased linearly onto main, with the same tree as 261019_mcStatFoldSyst (which WRemnants 261021_dev pins).

🤖 Generated with Claude Code

https://claude.ai/code/session_01GdFuNuyGWJYvaaMNBUGRmX

davidwalter2 and others added 6 commits November 10, 2026 13:46
…cture in flat bins

TensorWriter.add_systematic(..., fold_axis=...) with an up/down pair and
symmetrize="linear" or "quadratic" now splits every fold into SymAvg and
SymDiff, as for the full pair, instead of raising. These are the default
symmetrizations of the QCD scale, PDF, alphaS and quark-mass systematics in
WRemnants.

In the bins a folded process is written without fold structure (a negative
fold, or no MC-stat variance at all), every fold of a folded systematic now
gets the full logk (1/k of it for the normal form). Its per-fold logk was
computed against the flattened nominal (full / k), so the nominal's own fold
noise entered the systematic variation there.

The class attribute fold_symmetrize_modes lists the symmetrize modes of folded
pairs, so callers can fall back to a shared systematic. The sparse asymmetric
case is refused before anything is registered.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_01WmiWMzJgT1jHGqqLQ7dxro
…C events

In a fold with fewer effective MC events kstat = sumw^2/sumw2 than
TensorWriter(mcstat_fold_min_kstat=3) the per-fold variation of a folded
systematic is the ratio of two sums dominated by their noise, e.g. cancelling
negative weights: in a W mass test tensor (4 folds, 5 files per sample) such
folds had logk up to 25 against ~0.05 in the other folds, a variation by a
factor 1e6 at one sigma. Those folds now get the full logk (1/k of it for the
normal form), as the bins without fold structure; the nominal folds are
unchanged. Above 3 effective events the per-fold logk stayed within 0.7 of the
full one in that tensor. The number of such folds is logged per process.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_01WmiWMzJgT1jHGqqLQ7dxro
get_logk gives an empty fold (zero nominal) the sign-flip value log(1e-3),
while book_logk zeroes the logk of an empty full template. In a W mass test
tensor 40% of the per-fold entries were such folds. They do not change the
prediction (0 * exp(logk theta)), but overflow to 0 * inf = NaN for a large
trial value of the nuisance. Their per-fold logk is now zero, as for the full
template.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_01WmiWMzJgT1jHGqqLQ7dxro
…ates

The jackknife combination 2 L_full - 1/2 L_A - 1/2 L_B is a bias correction to
second order in the fold deviations x = E_f - E of a template's exponent. Beyond
that it is unbounded below, and templates of a few MC events with large weights
let the minimizer run away (W-mass fit: recoil_stat on Z->tautau and W->taunu;
with the isolation/mT sidebands the recoil calibration and isolation efficiency
stat variations on W->taunu templates of ~7 effective MC events).

- --mcStatFoldBound kappa (default 0.5): per (fold, bin, process) x is replaced
  by the taper x / (1 + x^2 / kappa^2). Value, gradient and curvature at the
  nominal point are unchanged, and far out a fold follows the full-sample
  variation, so the objective becomes the ordinary likelihood there.
  --mcStatFoldBoundForm tanh keeps kappa tanh(x / kappa) for comparison.
- --mcStatFoldMaxDeviation tau (default 0.1): the (bin, process) templates whose
  fold deviation norm max_f |logk_f - logk|_2 exceeds tau take the shared logk
  in every fold, since their k = 2 noise-floor estimate is not reliable.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_01WmiWMzJgT1jHGqqLQ7dxro
The (bin, process) templates with fewer effective MC events than the threshold,
kstat = (sum w)^2 / sum w^2, get the full-sample template in every fold
(nominal, MC-stat variance and logk), so they enter the de-biased objective like
the ordinary likelihood, without bias correction. With a few events a fold
estimate is noise (with negative weights a fold can cancel to zero while its
variations do not), which the jackknife cannot correct. Default 0 (off). The
number of templates taken out and their yield share are logged per process.

With --mcStatFoldMaxDeviation as well, the dense path holds two copies of the
folded logk for a moment.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_01WmiWMzJgT1jHGqqLQ7dxro
--mcStatFoldBound and --mcStatFoldMaxDeviation default to 0, like
--mcStatFoldMinKstat, so the k-fold de-biasing is the plain jackknife unless a
protection is requested. The W-mass fits use kappa = 0.5 and tau = 0.1. Without
the bound, the warning about the unbounded split-logk objective now names the
option that bounds it. The tests of the bound and of the shared-logk fallback
set them explicitly.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>
Claude-Session: https://claude.ai/code/session_01WmiWMzJgT1jHGqqLQ7dxro
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant