Skip to content

integrators: stabilize the MC error estimate (tail diagnostic, replicas, disclosed budgets) - #37

Closed
oshaughnessy-junior wants to merge 3 commits into
rift_O4dfrom
rift_O4d_mc_error_stabilization
Closed

integrators: stabilize the MC error estimate (tail diagnostic, replicas, disclosed budgets)#37
oshaughnessy-junior wants to merge 3 commits into
rift_O4dfrom
rift_O4d_mc_error_stabilization

Conversation

@oshaughnessy-junior

@oshaughnessy-junior oshaughnessy-junior commented Jul 29, 2026

Copy link
Copy Markdown
Owner

The reported sigma_lnL is the pooled weight sample-variance — algebraically 1/ESS_hat − 1/n restated, computed from the SAME weights as the integral. It is tail-blind: a run that misses the likelihood peak reports BOTH a low lnL AND a small sigma, so the error bar fails conditionally on the run being wrong. A controlled lognormal-weight study shows the median underestimate grows 1.0→2.4× (tau 1→5), 90% coverage collapses to 0.24, and the reported sigma is ANTI-correlated with the actual error; dslice exports corroborate (sigmaL understates redraw scatter 3–18× at high SNR).

Changes

  • statutils: pareto_khat_from_log (Zhang–Stephens GPD tail fit, PSIS xi convention), ess_from_log_weights, block_scatter_sigma (delete-one jackknife over per-chunk lnZ), bootstrap_lnZ_quantiles.
  • mcsampler / mcsamplerGPU: record per-chunk lnZ; at finalize emit pareto_khat, n_ESS, sigma_lnZ_block, and (relative error > 0.3) bootstrap lnZ quantiles via dict_return. Fixes the _rvs[weights] append-onto-joint_s_prior bug.
  • mcsamplerAdaptiveVolume: adds the missing Var(lnV) (stochastic live volume) and trunc_p^2 terms to the returned relative variance, components disclosed.
  • ILE batchmode: floors sigma at the between-chunk scatter; prints k-hat (>0.7 ⇒ sigma is a lower bound) and bootstrap CI; opt-in --mc-error-replicas N (+sigma/khat/ESS triggers) reruns cold replicas in-process and reports the LINEAR mean with sigma = max(propagated, replica scatter). Default off — production behavior unchanged unless enabled.
  • util_CleanILE: combines repeated evaluations by the ntot-weighted LINEAR mean with sigma = max(propagated, between-replica scatter). The previous 1/sigma^2 weighting overweighted exactly the replica that missed the peak (toy: −0.86 nat combine bias vs −0.26 linear).
  • mcsamplerGPU adaptive-proposal bias fix (ff0a04ba) — root-caused and closed the "−0.23 nat invisible bias" flagged in the first version of this PR. The 1D marginals of the drawn samples matched the claimed p_s exactly, but E[prior/p_s] = 0.37 on fresh draws: ~63% of the prior volume sat in histogram bins with exactly zero proposal probability. A zero bin is an absorbing state (it can never be re-drawn), so the sampled support shrinks irreversibly and the integral silently loses the mass outside it — invisible to every within-run error estimate. Four changes:
    • compute_hist: clamp the uniform-mixture floor to HIST_FLOOR_LEVEL_MIN = 1e-2 so no bin reaches zero (production ILE's --adapt-floor-level 0.1 already exceeds this; the clamp only binds for smaller/unset floors).
    • integrate_log: adaptation weights are now the stored tempered importance weights exp(tempering_exp*lnL + ln p − ln p_s), so the weighted histogram estimates the FIXED target L^beta * prior — the contract integrate() and the ILE driver comment always documented. The old lnL + max(maxlnL,200) weights ignored tempering_exp and 1/p_s; being near-flat they made each histogram replay the previous proposal's sampling noise, a multiplicative random walk that collapsed the proposal onto a comb of surviving bins (61/100 per dim) and drove the truncation.
    • integrate_log: the n_adapt freeze test double-multiplied by n (already scaled at parse time), so adaptation never froze at any setting; no-kwarg default scaled consistently.
    • pdf_from_hist: bin-index clamp off-by-one (right-edge sample indexed out of range).

Validation

  • 19/19 in the roboto-paper demos/error_estimate/run_validation.py (k-hat recovers GPD k=0.3/0.7 to ±0.04; all samplers emit diagnostics; CleanILE matches the new formula; non-adaptive AC unbiased).
  • Adaptive AC bias (2D Gaussian, n=2000, neff=1000, 16 runs/config): −0.321 ± 0.011 → −0.003 ± 0.004 for all n_adapt in {5..100} and tempering_exp in {0, 0.1, 1.0}; final proposal support 56% → 100%; E[prior/p_s] on fresh draws 0.373 → 1.0002; runs now reach neff≈1000 instead of exhausting nmax at neff≈350. T4b now reports delta = −0.016 against an honest sigma_block ≈ 0.009.
  • CI gate test_mcsamplerEnsemble_extended --as-test passes both variants unchanged.

🤖 Generated with Claude Code

oshaughnessy-junior and others added 2 commits July 29, 2026 08:31
…as, disclosed budgets)

The reported sigma_lnL is the pooled weight sample-variance -- algebraically
1/ESS_hat - 1/n restated, computed from the SAME weights as the integral.  It
is tail-blind: a run that misses the likelihood peak reports BOTH a low lnL
AND a small sigma, so the error bar fails conditionally on the run being wrong
(toy model: median underestimate 1.0-2.4x for tau=1-5, 90% coverage collapses
to 0.24, reported sigma ANTI-correlated with actual error; corroborated by
dslice exports where sigmaL understates redraw scatter 3-18x at high SNR).

Defenses, all disclosed in the job log:

* statutils: new helpers pareto_khat_from_log (Zhang-Stephens GPD tail fit,
  PSIS xi convention), ess_from_log_weights, block_scatter_sigma (delete-one
  jackknife over per-chunk lnZ), bootstrap_lnZ_quantiles.
* mcsampler + mcsamplerGPU: record per-chunk lnZ during the run; at finalize
  emit pareto_khat, n_ESS, sigma_lnZ_block, and (when relative error > 0.3)
  bootstrap lnZ quantiles via dict_return.  Also fix a bug where the
  _rvs['weights'] record appended new weights onto joint_s_prior.
* mcsamplerAdaptiveVolume: the naive term is blind to the stochastic live
  volume (Z ~ V*mean(w)) and to the truncated probability; accumulate
  Var(lnV) per cycle and add var_lnV + trunc_p^2 to the returned relative
  variance, with components in dict_return.
* ILE (integrate_likelihood_extrinsic_batchmode): floor the reported sigma at
  the between-chunk scatter; print k-hat (k>0.7: sigma is a lower bound) and
  bootstrap CI; new opt-in --mc-error-replicas N (+ --mc-error-sigma-trigger
  0.4, --mc-error-khat-trigger 0.7, --mc-error-ess-trigger 30): when
  triggered, rerun cold replicas in-process and report the LINEAR mean of the
  replica integrals with sigma = max(propagated, replica scatter), neff = sum;
  posterior export keeps the best-converged replica.  Default off.
* util_CleanILE: combine repeated evaluations by the ntot-weighted LINEAR
  mean with sigma = max(propagated, between-replica scatter).  The previous
  1/sigma^2 weighting overweighted exactly the replica that missed the peak
  (low lnL AND low sigma): toy model shows -0.86 nat combine bias vs -0.26
  for the linear mean.

Validation: 19/19 in the roboto-paper demos/error_estimate/run_validation.py
(k-hat recovers GPD k=0.3/0.7 to +-0.04; all samplers emit the diagnostics;
CleanILE matches the new formula) and .travis test_mcsamplerEnsemble_extended
--as-test passes both variants unchanged.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
On a mild 2D Gaussian test (T4b, floor_level unset), the adapted sampler was
biased by -0.32 nats while every within-run error estimate read ~0.02: the
1D marginals of the drawn samples matched the claimed p_s exactly, but
E[prior/p_s] = 0.37, i.e. ~63% of the prior volume sat in histogram bins with
exactly zero proposal probability.  A zero bin is an absorbing state (it can
never be re-drawn), so the sampled support shrinks irreversibly and the
integral silently loses the mass outside it.

Four fixes:
- compute_hist: clamp the uniform-mixture floor to HIST_FLOOR_LEVEL_MIN=1e-2
  so no bin can reach zero probability (production ILE already passes 0.1;
  the clamp only binds for smaller/unset floors).
- integrate_log: build adaptation weights from the stored tempered importance
  weights exp(tempering_exp*lnL + ln p - ln p_s) so the weighted histogram
  estimates the fixed target L^beta * prior -- the documented contract
  (see integrate() and the ILE driver comment).  The old lnL+max(maxlnL,200)
  weights ignored tempering_exp and 1/p_s; being near-flat, each histogram
  replayed the previous proposal's sampling noise, a multiplicative random
  walk that collapsed the proposal onto a comb of surviving bins (61/100 per
  dim in the test) and drove the truncation above.
- integrate_log: n_adapt freeze test double-multiplied by n (n_adapt was
  already scaled at parse time), so adaptation never froze regardless of the
  requested chunk count; also scale the no-kwarg default consistently.
- pdf_from_hist: clamp bin index to n_bins-1 (right-edge sample previously
  indexed out of range).

Validation (2D Gaussian, n=2000, neff=1000, 16 runs/config): bias
-0.321 +/- 0.011 -> -0.003 +/- 0.004 for all n_adapt in {5..100} and
tempering_exp in {0, 0.1, 1.0}; final proposal support 56% -> 100%;
runs now reach neff~1000 instead of exhausting nmax at neff~350.
E[prior/p_s] over fresh draws from the adapted proposal: 0.373 -> 1.0002.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@oshaughnessy-junior

Copy link
Copy Markdown
Owner Author

Auto-review: not clear to merge yet. This is a core integrator change (MC-error estimator paths in the standard, AV, and GPU samplers), so the mandatory CPU base-vs-candidate shape-recovery test is required. I verified a clean merge, changed-Python compilation, and test_mcsamplerEnsemble_extended.py --as-test locally.

However, this PR’s base predates the hard shape-gate suite and still has the eager optional-NF entry-point import issue, which prevents a valid strict portfolio gate without torch. GitHub also shows no CI checks for this head. Please rebase/depend on the branch with the hard gate and optional-plugin isolation, then provide a green CI run and standard base/candidate shape result before merge.

@oshaughnessy-junior
oshaughnessy-junior changed the base branch from rift_O4d_junior_calmarg_in_loop to rift_O4d July 29, 2026 16:19
Conflict: mcsamplerAdaptiveVolume.py integrate_log preamble -- kept both sides
(this branch's var_lnV accumulator init for the disclosed variance budget, and
rift_O4d's AV warm-start state re-application block).

Validation on the merged tree: 19/19 in the roboto-paper
demos/error_estimate/run_validation.py; adaptive-AC 2D Gaussian bias
+0.001 +/- 0.006 (8 runs).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@oshaughnessy-junior

Copy link
Copy Markdown
Owner Author

Superseded by #38: same single commit rebased onto rift_O4d (which now carries the hard shape-gate suite and optional-plugin isolation), retargeted to rift_O4d per review, with the required base-vs-candidate shape-recovery gate attached there (0 blocking regressions; base 8 strict failures vs PR 4 -- the delta is honest error bars widening sigma-scaled tolerances, biases byte-identical; details in #38). Closing in favor of #38.

🤖 Generated with Claude Code

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