diff --git a/meson.build b/meson.build index 8160418..babc82d 100644 --- a/meson.build +++ b/meson.build @@ -43,6 +43,7 @@ sources_list = [ 'mir/stat/internal/one_minus_exp', 'mir/stat/internal/shape_transform', 'mir/stat/internal/standardize', + 'mir/stat/internal/summation', 'mir/stat/descriptive/histogram/package', 'mir/stat/descriptive/histogram/api/rc', 'mir/stat/descriptive/histogram/api/gc', diff --git a/source/mir/stat/descriptive/univariate/dispersions.d b/source/mir/stat/descriptive/univariate/dispersions.d index ee3a2a3..4ba5d90 100644 --- a/source/mir/stat/descriptive/univariate/dispersions.d +++ b/source/mir/stat/descriptive/univariate/dispersions.d @@ -1087,6 +1087,19 @@ By default, if `F` is not floating point type or complex type, then the result will have a `double` type if `F` is implicitly convertible to a floating point type or a type for which `isComplex!F` is true. +Large integer observations that differ by small amounts can lose those +differences when converted to the floating-point accumulator type. For example, +the five `long` values `(1L << 60) + [0, 1, 2, 3, 4]` have sample variance +`2.5`, but the default calculation returns zero because their conversions to +`double` are identical. A more accurate summation algorithm cannot recover +differences already lost during conversion. + +When the data has a known common offset, subtracting it in integer arithmetic +before calling `variance` can preserve the differences. Ensure that the +subtractions fit the integer type. This translation preserves variance about +the mean, but must not be used with `VarianceAlgo.assumeZeroMean`, which +calculates squared distances from zero. + Params: F = controls type of output varianceAlgo = algorithm for calculating variance (default: VarianceAlgo.hybrid) @@ -1430,12 +1443,101 @@ unittest assert(x.sliced.variance!float.approxEqual(54.76562 / 11)); } +/* +Recover a representable standard deviation when squaring the observations or +residuals overflowed or underflowed. Scaling by a power of two keeps nearby +observations distinct; subtracting a fixed normalized observation avoids a +rounded mean creating a spread in constant inputs. The original variance is +computed first, so its contracts and ordinary results remain unchanged. +*/ +private T recoverStandardDeviation(T, VarianceAlgo algo, Range)( + Range values, bool isPopulation, T originalVariance) +{ + import std.math : abs, frexp, ldexp, isFinite, sqrt; + import std.range.primitives : empty, front, save; + + if (values.empty) + return sqrt(originalVariance); + const first = cast(T) values.front; + // Zero can mean either identical observations or underflowed squares. + // Equality is cheaper to check than computing the scale for every value. + if (originalVariance == 0 && isFinite(first)) + { + static if (algo == VarianceAlgo.assumeZeroMean) + const canBeConstant = first == 0; + else + enum canBeConstant = true; + if (canBeConstant) + { + bool equal = true; + foreach (x; values.save) + if (cast(T) x != first) + { + equal = false; + break; + } + if (equal) + return 0; + } + } + + T maximum = 0; + foreach (x; values.save) + { + const value = cast(T) x; + if (!isFinite(value)) + return sqrt(originalVariance); + if (abs(value) > maximum) + maximum = abs(value); + } + if (maximum == 0) + return sqrt(originalVariance); + int exponent; + frexp(maximum, exponent); + const scale = ldexp(T(0.5), exponent); + static if (algo == VarianceAlgo.assumeZeroMean) + enum anchor = T(0); + else + const anchor = first / scale; + size_t count; + Summator!(T, Summation.kbn) sum; + foreach (x; values.save) + { + sum.put(cast(T) x / scale - anchor); + ++count; + } + static if (algo == VarianceAlgo.assumeZeroMean) + enum center = T(0); + else + const center = sum.sum / count; + Summator!(T, Summation.kbn) squares; + foreach (x; values) + { + const residual = (cast(T) x / scale - anchor) - center; + squares.put(residual * residual); + } + const v = squares.sum / (count + isPopulation - 1); + return sqrt(v) * scale; +} + /++ Calculates the standard deviation of the input By default, if `F` is not floating point type, then the result will have a `double` type if `F` is implicitly convertible to a floating point type. +For replayable floating-point inputs, a nonfinite, negative, zero, or +subnormal variance triggers a scaled recalculation. This can preserve a +standard deviation even when its square is not representable. Recovery uses +compensated two-pass variance, except that `assumeZeroMean` retains its +zero-centered definition. Ordinary positive normal variances use the selected +algorithm unchanged. Single-pass inputs do not receive this recovery. + +Integer inputs retain the conversion limitation described for $(LREF variance): +large integers separated by small differences can produce a zero standard +deviation. The scaled floating-point recovery does not repair this loss of +integer precision. + Params: F = controls type of output varianceAlgo = algorithm for calculating variance (default: VarianceAlgo.hybrid) @@ -1460,7 +1562,25 @@ template standardDeviation( { import core.lifetime: move; alias G = typeof(return); - return r.move.variance!(G, varianceAlgo, ResolveSummationType!(summation, Range, G))(isPopulation).sqrt; + import std.range.primitives : isForwardRange, ElementType, save; + import std.traits : isFloatingPoint; + static if (isSlice!Range) + { + import mir.ndslice.topology : flattened; + auto flat = r.flattened; + static if (isForwardRange!(typeof(flat)) && isFloatingPoint!(ElementType!(typeof(flat)))) + auto saved = flat.save; + } + else static if (isForwardRange!Range && isFloatingPoint!(ElementType!Range)) + auto saved = r.save; + const v = r.move.variance!(G, varianceAlgo, ResolveSummationType!(summation, Range, G))(isPopulation); + static if (__traits(compiles, saved)) + { + import std.math : isFinite; + if (!isFinite(v) || v < G.min_normal) + return recoverStandardDeviation!(G, varianceAlgo)(saved.move, isPopulation, v); + } + return v.sqrt; } /++ @@ -1470,10 +1590,120 @@ template standardDeviation( @fmamath stdevType!F standardDeviation(scope const F[] ar...) { alias G = typeof(return); - return ar.variance!(G, varianceAlgo, ResolveSummationType!(summation, const(G)[], G)).sqrt; + const v = ar.variance!(G, varianceAlgo, ResolveSummationType!(summation, const(G)[], G)); + static if (isFloatingPoint!F) + { + import std.math : isFinite; + if (!isFinite(v) || v < G.min_normal) + return recoverStandardDeviation!(G, varianceAlgo)(ar, false, v); + } + return v.sqrt; } } +// Save the forward cursor independently before computing variance. +version(mir_stat_test) +unittest +{ + static struct Forward + { + const(double)[] data; + size_t* position; + size_t ownPosition; + bool empty() const { return (position ? *position : ownPosition) == data.length; } + double front() const { return data[position ? *position : ownPosition]; } + void popFront() { if (position) ++*position; else ++ownPosition; } + auto save() { return Forward(data, null, position ? *position : ownPosition); } + } + static struct Single + { + const(double)[] data; + size_t* position; + bool empty() const { return *position == data.length; } + double front() const { return data[*position]; } + void popFront() { ++*position; } + } + size_t position; + // Standard-deviation recovery also needs an independent saved cursor. + double[3] spread = [-1e200, 0, 1e200]; + position = 0; + assert(standardDeviation!(VarianceAlgo.online)(Forward(spread[], &position)) == 1e200); + assert(position == 3); + position = 0; + assert(standardDeviation!(VarianceAlgo.online)(Single(spread[], &position)) == double.infinity); + assert(position == 3); +} + +// Preserve assumeZeroMean's existing boundary results when recovery has no scale. +// Unlike the centered algorithms, this algorithm has no minimum-count contract. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.math : isNaN; + + double[] empty; + // Population variance divides zero by zero; recovery must not read front. + assert(empty.standardDeviation!(VarianceAlgo.assumeZeroMean)(true).isNaN); + // A singleton sample also divides zero by zero. All observations are zero, + // so scaling cannot recover a standard deviation from the NaN variance. + double[1] zero = [0.0]; + assert(zero[].standardDeviation!(VarianceAlgo.assumeZeroMean).isNaN); + assert(zero[].standardDeviation!(VarianceAlgo.assumeZeroMean)(true) == 0); +} + +// Recover standard deviations whose squared values are outside the type's range. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.meta : AliasSeq; + import std.math : sqrt, abs, isNaN; + import mir.ndslice.slice : sliced; + import mir.ndslice.dynamic : transposed; + + static foreach (T; AliasSeq!(float, double, real)) + {{ + foreach (magnitude; [T.max / 4, T.min_normal, T.min_normal * T.epsilon]) + { + T[3] values = [-magnitude, 0, magnitude]; + static foreach (algo; AliasSeq!(VarianceAlgo.naive, VarianceAlgo.online, + VarianceAlgo.twoPass, VarianceAlgo.hybrid, VarianceAlgo.assumeZeroMean)) + {{ + assert(values[].standardDeviation!(T, algo) == magnitude); + const population = values[].standardDeviation!(T, algo)(true); + const expected = magnitude * sqrt(T(2) / 3); + assert(abs(population - expected) <= 4 * T.epsilon * magnitude + + T.min_normal * T.epsilon); + }} + } + T[3] extremes = [-T.max, 0, T.max]; + assert(extremes[].standardDeviation!(T, VarianceAlgo.naive) == T.max); + assert(standardDeviation!T(-T.max, T(0), T.max) == T.max); + T[3] largeConstant; + largeConstant[] = T.max / 4; + const rms = largeConstant[].standardDeviation!(T, VarianceAlgo.assumeZeroMean)(true); + assert(abs(rms / (T.max / 4) - 1) <= 2 * T.epsilon); + T[5] identical; + identical[] = T.max; + assert(identical[].standardDeviation == 0); + T[3] zero; + zero[] = 0; + assert(zero[].standardDeviation == 0); + T[3] constant = [T(3), T(3), T(3)]; + assert(constant[].standardDeviation!(T, VarianceAlgo.assumeZeroMean)(true) == 3); + T[2] invalid = [T.infinity, T.infinity]; + assert(invalid[].standardDeviation.isNaN); + }} + double[3] tiny = [-2e-162, 0, 2e-162]; + assert(abs(tiny[].standardDeviation / 2e-162 - 1) <= 2 * double.epsilon); + const double[6] matrix = [-1e200, 0, 1e200, -1e200, 0, 1e200]; + assert(matrix[].sliced(2, 3).transposed[0].standardDeviation == 0); + const expected = 1e200 * sqrt(4.0 / 5); + assert(abs(matrix[].sliced(2, 3).transposed.standardDeviation / expected - 1) + <= 4 * double.epsilon); +} + /// ditto template standardDeviation( VarianceAlgo varianceAlgo = VarianceAlgo.hybrid, @@ -1605,11 +1835,9 @@ unittest auto y = x.standardDeviation; assert(y.approxEqual(sqrt(54.76562 / 11))); - // The naive algorithm is numerically unstable in this case - auto z0 = x.standardDeviation!"naive"; - assert(!z0.approxEqual(y)); - - // But the two-pass algorithm provides a consistent answer + // Prefer two-pass for large offsets. Naive variance can lose the spread; + // recovery catches zero or negative results, but not an inaccurate positive result. + // The two-pass algorithm provides a consistent answer auto z1 = x.standardDeviation!"twoPass"; assert(z1.approxEqual(y)); } diff --git a/source/mir/stat/descriptive/univariate/means.d b/source/mir/stat/descriptive/univariate/means.d index 9b59949..d214f50 100644 --- a/source/mir/stat/descriptive/univariate/means.d +++ b/source/mir/stat/descriptive/univariate/means.d @@ -102,8 +102,10 @@ struct MeanAccumulator(T, Summation summation) /// void put(F = T)(MeanAccumulator!(F, summation) m) { + import mir.stat.internal.summation : mergeSummator; + count += m.count; - summator.put(cast(T) m.summator); + mergeSummator(summator, m.summator); } } @@ -150,6 +152,137 @@ unittest assert(m0.mean == 29.25 / 12); } +// Merging compensated sums must retain contributions below the rounded sum. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.meta : AliasSeq; + + static foreach (mode; AliasSeq!(Summation.kbn, Summation.kb2)) + {{ + MeanAccumulator!(double, mode) source, merged; + source.put(0x1p54); + source.put(1.0); + merged.put(source); + merged.put(-0x1p54); + assert(merged.count == 3); + assert(merged.sum == 1); + assert(merged.mean == 1.0 / 3); + + // The merge must not consume or alter the source's corrections. + source.put(-0x1p54); + assert(source.sum == 1); + assert(source.count == 3); + }} +} + +// Retain scalar merging for summators without a direct merge operation. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + alias Accumulator = MeanAccumulator!(double, Summation.decimal); + Accumulator source, merged; + source.put(1.0); + source.put(3.0); + merged.put(2.0); + merged.put(source); + assert(merged.count == 3 && merged.sum == 6); + assert(merged.mean == 2); + assert(source.count == 2 && source.sum == 4); + merged.put(Accumulator.init); + assert(merged.count == 3 && merged.sum == 6); + merged.put(merged); + assert(merged.count == 6 && merged.sum == 12); + assert(merged.mean == 2); +} + +// Precise summation also retains overflow internally until later cancellation. +// Copying this summator by value currently requires @system in mir-algorithm. +version(mir_stat_test) +pure nothrow @nogc +unittest +{ + alias Accumulator = MeanAccumulator!(double, Summation.precise); + Accumulator source, merged; + source.put(double.max); + source.put(double.max); + source.put(1.0); + merged.put(source); + merged.put(-double.max); + merged.put(-double.max); + assert(merged.count == 5); + assert(merged.sum == 1); + assert(merged.mean == 0.2); + + source.put(-double.max); + source.put(-double.max); + assert(source.sum == 1); + + // Empty and self merges preserve the accumulated state and count. + merged.put(Accumulator.init); + assert(merged.count == 5 && merged.sum == 1); + merged.put(merged); + assert(merged.count == 10 && merged.sum == 2); + assert(merged.mean == 0.2); +} + +// Retry a nonfinite floating-point mean without overflowing the raw sum. +// The caller supplies an independently saved traversal of the same observations. +private T recoverOverflowingMean(T, Range)(scope Range values, size_t count, T original) +{ + import std.math : abs, isFinite; + + if (!count) + return original; + + // A power of two at least twice the count leaves headroom for the large + // terms' sum. Multiplication by its reciprocal is exact unless it enters + // the subnormal range; keep those small terms in a separate unscaled sum. + T scale = 2; + for (size_t n = count - 1; n; n >>= 1) + scale *= 2; + const T inverse = 1 / scale; + const T threshold = T.min_normal * scale; + Summator!(T, Summation.precise) high, low; + T lower = T.infinity, upper = -T.infinity; + foreach (value; values) + { + const T x = cast(T) value; + // Genuine NaN/Inf inputs (including overflow on conversion to T) are + // not finite-input sum overflow. Preserve the original result. + if (!isFinite(x)) + return original; + if (x < lower) lower = x; + if (x > upper) upper = x; + if (abs(x) < threshold) + low.put(x); + else + high.put(x * inverse); + } + + const T highSum = high.sum; + const T unscaled = highSum * scale; + T result; + if (isFinite(unscaled)) + { + // Cancellation can cross the two sums. Add the large subtotal to the + // still-precise small sum before rounding, so a tiny remainder is not + // lost by rounding the two opposing subtotals independently. + low.put(unscaled); + result = low.sum / cast(T) count; + } + else + { + // The unscaled subtotal itself overflows: divide before restoring its + // scale. At this magnitude the small sum cannot cancel the large one. + result = highSum * (scale / cast(T) count) + low.sum / cast(T) count; + } + // Rounding on rescaling must not leave the finite input interval. + return result < lower ? lower : result > upper ? upper : result; +} + /++ Computes the mean of the input. @@ -157,6 +290,18 @@ By default, if `F` is not floating point type or complex type, then the result will have a `double` type if `F` is implicitly convertible to a floating point type or a type for which `isComplex!F` is true. +For built-in floating-point results, a nonfinite result is retried with +scaled precise summation when the input can be traversed again. This supports +arrays, slices whose flattened view is a forward range, and other forward +ranges. For example, the mean of two maximum finite values is still that +maximum value, even though their sum overflows. A forward range's saved +traversal must yield the same observations as the first traversal. + +The retry preserves the original result if an observation is nonfinite after +conversion to the result type. Finite results use the selected summation +algorithm without a retry; this is not a general correction for cancellation +error. Single-pass inputs and $(LREF MeanAccumulator) do not retry. + Params: F = controls type of output summation = algorithm for calculating sums (default: Summation.appropriate) @@ -169,7 +314,7 @@ See_also: template mean(F, Summation summation = Summation.appropriate) { import core.lifetime: move; - import std.traits: isIterable; + import std.traits: isIterable, isFloatingPoint; /++ Params: @@ -180,8 +325,28 @@ template mean(F, Summation summation = Summation.appropriate) { alias G = typeof(return); MeanAccumulator!(G, ResolveSummationType!(summation, Range, G)) mean; + static if (isFloatingPoint!G) + { + import std.range.primitives : isForwardRange, save; + static if (isSlice!Range) + { + import mir.ndslice.topology : flattened; + auto flat = r.flattened; + static if (isForwardRange!(typeof(flat))) + auto saved = flat.save; + } + else static if (isForwardRange!Range) + auto saved = r.save; + } mean.put(r.move); - return mean.mean; + const result = mean.mean; + static if (__traits(compiles, saved)) + { + import std.math : isFinite; + if (!isFinite(result)) + return recoverOverflowingMean!G(saved.move, mean.count, result); + } + return result; } /++ @@ -193,7 +358,14 @@ template mean(F, Summation summation = Summation.appropriate) alias G = typeof(return); MeanAccumulator!(G, ResolveSummationType!(summation, const(G)[], G)) mean; mean.put(ar); - return mean.mean; + const result = mean.mean; + static if (isFloatingPoint!G) + { + import std.math : isFinite; + if (!isFinite(result)) + return recoverOverflowingMean!G(ar, mean.count, result); + } + return result; } } @@ -237,6 +409,93 @@ template mean(string summation) mixin("alias mean = .mean!(Summation." ~ summation ~ ");"); } +// Finite observations can have a representable mean but an overflowing sum. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + import std.meta : AliasSeq; + import std.math : isNaN; + import mir.ndslice.slice : sliced; + import mir.ndslice.dynamic : transposed; + + static foreach (T; AliasSeq!(float, double, real)) + static foreach (mode; AliasSeq!(Summation.naive, Summation.pairwise, + Summation.kbn, Summation.kb2, Summation.precise)) + {{ + T[4] values = [T.max, T.max, T.max, T.max]; + assert(mean!(T, mode)(values[]) == T.max); + assert(mean!(T, mode)(T.max, T.max) == T.max); + assert(mean!(T, mode)(values[].sliced) == T.max); + assert(mean!(T, mode)(values[].sliced(2, 2).transposed) == T.max); + assert(mean!(T, mode)(values[].sliced(2, 2).transposed[0]) == T.max); + values[] = -T.max; + assert(mean!(T, mode)(values[]) == -T.max); + + T[] empty; + assert(isNaN(mean!(T, mode)(empty))); + assert(isNaN(mean!(T, mode)(T.nan, T(1)))); + // Preserve the selected summator's handling of infinities. + MeanAccumulator!(T, mode) original; + original.put(T.infinity); + original.put(T(1)); + const actual = mean!(T, mode)(T.infinity, T(1)); + assert(isNaN(original.mean) ? isNaN(actual) : original.mean == actual); + }} + + static foreach (T; AliasSeq!(float, double, real)) + {{ + // Dividing every input first would erase this subnormal remainder. + const T tiny = T.min_normal * T.epsilon; + T[9] values = [T.max, T.max, -T.max, -T.max, + tiny, tiny, tiny, tiny, tiny]; + assert(mean!(T, Summation.naive)(values[]) == tiny); + + // Cancellation between the scaled and unscaled groups must happen + // before their subtotals are rounded. The exact mean is 5/9 * tiny, + // which rounds to tiny rather than zero. + const T boundary = T.min_normal * 32; + T[9] cross = [T.max, T.max, -T.max, -T.max, + boundary, -boundary / 2, -boundary / 2, 5 * tiny, 0]; + assert(mean!(T, Summation.naive)(cross[]) == tiny); + cross[] = -cross[]; + assert(mean!(T, Summation.naive)(cross[]) == -tiny); + }} + // Conversion overflow is not repaired as if it were sum overflow. + assert(mean!float(double.max, double.max) == float.infinity); +} + +// Save forward ranges before the first pass advances a shared cursor. +version(mir_stat_test) +unittest +{ + static struct Forward + { + const(double)[] data; + size_t* position; + size_t ownPosition; + bool empty() const { return (position ? *position : ownPosition) == data.length; } + double front() const { return data[position ? *position : ownPosition]; } + void popFront() { if (position) ++*position; else ++ownPosition; } + auto save() { return Forward(data, null, position ? *position : ownPosition); } + } + static struct Single + { + const(double)[] data; + size_t* position; + bool empty() const { return *position == data.length; } + double front() const { return data[*position]; } + void popFront() { ++*position; } + } + double[2] values = [double.max, double.max]; + size_t position; + assert(mean(Forward(values[], &position)) == double.max); + assert(position == 2); + position = 0; + assert(mean(Single(values[], &position)) == double.infinity); + assert(position == 2); +} + /// version(mir_stat_test) @safe pure nothrow diff --git a/source/mir/stat/internal/summation.d b/source/mir/stat/internal/summation.d new file mode 100644 index 0000000..756ed6b --- /dev/null +++ b/source/mir/stat/internal/summation.d @@ -0,0 +1,68 @@ +/++ +Internal summator merging for statistical accumulators. + +`mergeSummator` combines accumulated summation state. Matching summator types +are merged directly when supported, preserving partial sums and compensation +terms. Otherwise, the source sum is converted to the destination element type +before being added. Callers manage observation counts and any additional +statistical corrections. + +The helper has `package(mir.stat)` visibility and is intended for use inside +mir.stat; it is not a public application API. See +$(MREF mir,stat,descriptive,univariate) for the public statistical functions. + +License: $(HTTP www.apache.org/licenses/LICENSE-2.0, Apache-2.0) +Authors: John Michael Hall +Copyright: 2026 Mir Stat Authors. ++/ +module mir.stat.internal.summation; + +import mir.math.sum : Summator, Summation; + +/* +Preserve partial sums and corrections when matching summator types support a +direct merge. Otherwise retain scalar conversion to the destination element +type. In particular, decimal summators currently lack a direct merge, and +cross-type merges retain the existing conversion path and its constraints. + +Capability detection is compile-time only. Counts and any statistical merge +corrections remain the caller's responsibility. Source and target must be +distinct summator objects; callers supporting self merges must first copy the +source, as MeanAccumulator's by-value merge parameter does. The source is not +modified. +*/ +package(mir.stat) void mergeSummator(T, U, Summation mode)( + ref Summator!(T, mode) target, ref const Summator!(U, mode) source) +{ + static if (is(T == U) && __traits(compiles, target += source)) + target += source; + else + target.put(cast(T) source); +} + +// Both paths accept a const source without changing its accumulated state. +version(mir_stat_test) +@safe pure nothrow @nogc +unittest +{ + Summator!(double, Summation.kbn) source; + source.put(0x1p54); + source.put(1.0); + const saved = source; + Summator!(double, Summation.kbn) target; + mergeSummator(target, saved); + target.put(-0x1p54); + assert(target.sum == 1); + // The source's saved correction remains intact. + source.put(-0x1p54); + assert(source.sum == 1); + + Summator!(double, Summation.decimal) decimalSource; + decimalSource.put(2.0); + const savedDecimal = decimalSource; + Summator!(double, Summation.decimal) decimalTarget; + decimalTarget.put(1.0); + mergeSummator(decimalTarget, savedDecimal); + assert(decimalTarget.sum == 3); + assert(savedDecimal.sum == 2); +}