fix: accumulate central moments relative to the first observation#408
Open
agene0001 wants to merge 2 commits into
Open
fix: accumulate central moments relative to the first observation#408agene0001 wants to merge 2 commits into
agene0001 wants to merge 2 commits into
Conversation
`OnlineMoments` (added in statrs-dev#394 and wired into `Statistics` in cf11836 for this issue) cured the catastrophic case statrs-dev#376 opened with, but Welford is still poorly conditioned when the data carries a large offset: `mean += delta / n` cannot represent a small increment against a large running mean, so the low bits of every update are dropped. On the dataset from the issue - `1e12 + U(0, 1)`, n = 1e6 - measured against a Neumaier-compensated two-pass reference of 8.336923e-2: before 2.5e-4 relative error after 5e-15 Central moments are invariant under a shift, so accumulating moments of `x - first_observation` keeps every magnitude small. It costs one subtraction per observation - 1107 us vs 1103 us over 1e6 elements, i.e. free - and is better conditioned than plain Welford, which carries the same ~2.5e-4 here because it has the same mean-update problem. Also adds `OnlineMoments::merge`, the Chan-Golub-LeVeque pairwise update, which is the parallelisability the issue asks for. It has to reconcile two different offsets, so it is reviewed alongside them. Folding into several accumulators and merging is also slightly better conditioned than one long chain, since each chain accumulates over fewer updates. Refs statrs-dev#376
The tests collected their inputs into `Vec`, which does not exist under `--no-default-features`: the crate is `no_std` without `alloc`. Folding straight off the iterators gives the same accumulations with no allocation; `merge_reconciles_different_offsets` needs the data three times, so it takes closures returning fresh iterators rather than one collected copy. CI did not catch this because the test job runs only default features, and the feature-powerset job passes `--no-dev-deps`, so test code is never compiled without std. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Follow-up to #394 / cf11836, which moved
Statisticsonto Welford for this issue. That cured the catastrophic case #376 opened with, but Welford is still poorly conditioned when the data carries a large offset:mean += delta / ncannot represent a small increment against a large running mean, so the low bits of every update are dropped.On the dataset from the issue —
1e12 + U(0,1), n = 1e6 — measured against a Neumaier-compensated two-pass reference:varianceCentral moments are invariant under a shift, so
OnlineMomentsnow accumulates moments ofx − first_observation, keeping every magnitude small. Cost is one subtraction per observation: 1107 µs vs 1103 µs over 1e6 elements, i.e. free. (For comparison, I also prototyped a double-double running mean: 2x slower and less accurate, 1.9e-8.)Also adds
OnlineMoments::merge(Chan–Golub–LeVeque pairwise update, Pébay's form for M3), which is the parallelisability this issue asks for. It reconciles the two accumulators' different offsets, so it is included here rather than separately. Folding into several accumulators and merging is also slightly better conditioned than one long chain.Tests: exact-data offset sweep at offsets 2^0 … 2^40 (powers of two + a 2^-10 step so every sample is exactly representable and the reference variance is exact), merge-vs-single-chain equality at every split point, and merge with differing offsets.
Closes #376
🤖 Generated with Claude Code