Scaling: never let an observation's own fluctuation set its weight

A weighted mean is only unbiased while the weights are independent of the values
being averaged. The IUCr's own nomenclature report (Schwarzenbach et al., Acta
Cryst A45 (1989) 63-75) puts it directly: weights in averaging "should not be
based on the counting statistics of the individual observations whose estimated
variances are biased and result in larger weights for accidentally low
intensities". Two places in the rotation pipeline were doing exactly that, and
between them they drove whole resolution shells of merged intensity negative.

1. The profile fit computed its non-signal variance as

       var_bkg = max(0, 1/den - max(0, I) + bkg-estimate term)

   The point of a separate var_bkg is that it does NOT move with the
   reflection's own fluctuation, and 1/den - I is the quantity that does not:
   1/den is the fit variance taken at the fitted intensity and grows with it
   roughly one for one. Clamping the subtrahend at zero left a down-fluctuated
   reflection's own deflated variance standing as its background variance.
   Measured over 6.9 M partials of one weak rotation dataset, var_bkg/bkg came
   out at 3.7-5.4 for observations with I < 0 against 11.4-13.7 for I > 0 - the
   down-fluctuated half of every reflection carried a variance ~2.7x too small
   and was weighted up by the same factor, first in the 3D combine and then
   again in the merge. Removing the clamp makes var_bkg flat in I (~13 x bkg
   across the whole range).

2. The merge then weighted each combined full by 1/sigma_full^2, and sigma_full
   is by construction a function of the full's own answer: the combine's
   variance carries a corr*max(0, F) signal term, so every full with F <= 0 got
   the smallest variance the model allows while the strongest quartile got
   2.26x more. The merge now rebuilds that variance at the reflection's mean
   instead, from a linear model var(I) = var_bkg + var_per_I * I that the
   combine measures and stores on the full. This mirrors
   MergeOnTheFly::CorrectedSigma, whose comment already claimed to mirror the
   rotation combine.

Verified against an estimator that cannot see the fluctuation - summing the
partials and dividing by the summed partiality, the classical construction every
other program uses (Greenhough & Suddath, J. Appl. Cryst. 19 (1986) 400-409, via
Leslie, Acta Cryst D55 (1999) 1696-1702: profile fitting biases the individual
partials but not their sum). Reproducing the merge on dumped observations, the
shipped weighting sat ~1.9 sigma below that reference in the noise shells; the
two changes recover most of it, and every intensity-independent weighting
scheme agrees with the reference once (1) is in.

Four-crystal probe, XDS resolution limits, branch fingerprint identical on all
four (so none of these is a two-pass branch flip):

  weak cubic case   last shell <I/sig> -1.6 -> +0.2 (XDS +0.10), last shell
                    R_meas 478% -> 250% (XDS 246%), overall <I/sig> 6.1 -> 7.5
                    (XDS 7.18), R_meas 18.3% -> 18.1%, CC1/2_hi 38.2% -> 43.7%
  tetragonal case   outer shells <I/sig> -0.4/-0.8/-0.9/-1.0 -> +1.8/+1.2/
                    +0.9/+0.4, R_meas 184%/595%/7614%/nan -> 95%/119%/135%/232%
                    (the nan was the shell mean crossing zero), R_meas 33.3% ->
                    32.9%, CC1/2_hi 38.3% -> 56.5%
  trigonal case     R_meas 13.0% -> 12.5%, CC1/2_hi 14.4% -> 16.5%
  strong control    unchanged to every printed digit but ISa

Cost: ISa falls (17.2 -> 14.0 and 16.7 -> 14.9 on the two mid-strength cases,
28.3 -> 27.8 on the control). Strong reflections are untouched by (1) - their
partials are all positive, so var_bkg is bit-identical - but the joint a/b fit
redistributes: honest weak sigmas lower a, and b rises to keep the strong bins
fitted. The median reduced chi^2 improves (1.25 -> 1.14, 1.35 -> 1.28) so the
new split describes the scatter better, but ISa is the one headline metric that
moves the wrong way and it should be watched over the full battery.

The integrator change is shared, so the stills merge sees it too; there it feeds
GetExpectedVarianceMerge, which had been handed the same contaminated var_bkg.
That path is untested here.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
2026-08-11 17:40:15 +02:00
co-authored by Claude Opus 5
parent 1eb9d5a2fa
commit d8029524e7
7 changed files with 82 additions and 24 deletions
@@ -511,7 +511,13 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's
// own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing.
double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var);
double var_bkg = std::max(0.0, 1.0 / den - std::max(0.0, I)
// var_bkg is the NON-SIGNAL part of that variance, and 1/den is the fit's variance taken at the
// fitted intensity, so the signal part to remove is I itself - not max(0, I). Clamping it leaves
// a down-fluctuated reflection's own (deflated) variance standing as its background variance,
// which is 2-3x too small; the merge then weights exactly the down-fluctuated observations up.
// The whole point of a separate var_bkg is that it does not move with the observation's own
// fluctuation, and 1/den - I is what does not (1/den grows with I one for one).
double var_bkg = std::max(0.0, 1.0 / den - I
+ (wsum / den) * (wsum / den) * rh.bkg_var);
// The seed is a sum over the disk the box sum actually read, so when Exclude has taken pixels
// out of both, the fit's full-profile intensity has to be scaled down to that same disk