From 9403688f00f1dfa897617aef49cf449f01a9ac5f Mon Sep 17 00:00:00 2001 From: alxkm <19151554+alxkm@users.noreply.github.com> Date: Thu, 10 Sep 2026 17:57:34 +0200 Subject: [PATCH] feat: add WelfordAlgorithm, online mean and variance in one pass Welford's recurrence keeps a running mean and the sum of squared deviations from it, so it never forms the large nearly equal intermediate values that make the textbook variance formula lose its significant digits, and it needs O(1) time per sample and O(1) memory regardless of the stream length. Beyond the plain accumulation it supports removal, which runs the recurrence backwards and turns the accumulator into the statistics of a sliding window, and a static merge implementing Chan's parallel update so partial results from different shards combine exactly. Signed-off-by: alxkm <19151554+alxkm@users.noreply.github.com> --- .../streaming/WelfordAlgorithm.java | 245 ++++++++++++++++++ .../streaming/WelfordAlgorithmTest.java | 243 +++++++++++++++++ 2 files changed, 488 insertions(+) create mode 100644 src/main/java/com/thealgorithms/streaming/WelfordAlgorithm.java create mode 100644 src/test/java/com/thealgorithms/streaming/WelfordAlgorithmTest.java diff --git a/src/main/java/com/thealgorithms/streaming/WelfordAlgorithm.java b/src/main/java/com/thealgorithms/streaming/WelfordAlgorithm.java new file mode 100644 index 000000000000..62752df3d7da --- /dev/null +++ b/src/main/java/com/thealgorithms/streaming/WelfordAlgorithm.java @@ -0,0 +1,245 @@ +package com.thealgorithms.streaming; + +/** + * Online (single pass) mean and variance using Welford's algorithm. + * + *

The textbook formula {@code Var = (sum(x^2) - n * mean^2) / (n - 1)} is fast but numerically + * treacherous: {@code sum(x^2)} and {@code n * mean^2} may be huge and nearly equal, so their + * difference loses most of its significant digits and can even come out negative. Welford's + * recurrence never forms those large intermediate values. It keeps only the running mean and the sum + * of squared deviations from that running mean, {@code M2}: + * + * 

+ * n     <- n + 1
+ * delta <- x - mean
+ * mean  <- mean + delta / n
+ * M2    <- M2 + delta * (x - mean)   // note: the second factor uses the *updated* mean
+ *
+ * + *

Both {@link #add(double)} and {@link #remove(double)} run in O(1) time and the accumulator + * occupies O(1) memory no matter how many samples pass through it. + * + *

Sliding windows and map-reduce

+ * + * + * + *

Usage

+ * + * 
{@code
+ * WelfordAlgorithm stats = new WelfordAlgorithm();
+ * stats.add(2.0);
+ * stats.add(4.0);
+ * stats.add(4.0);
+ * stats.mean();                        // 3.3333...
+ * stats.populationStandardDeviation();  // 0.9428...
+ * }
+ * + *

This class is not thread-safe. + * + * @see ExponentialMovingAverage + * @see Algorithms for calculating variance + */ +public final class WelfordAlgorithm { + + private long count; + private double mean; + private double sumOfSquaredDeviations; + + /** + * Creates an empty accumulator. + */ + public WelfordAlgorithm() { + clear(); + } + + /** + * Incorporates one sample. + * + * @param value the sample to add + * @throws IllegalArgumentException if {@code value} is NaN or infinite + */ + public void add(double value) { + requireFinite(value); + count++; + double delta = value - mean; + mean += delta / count; + sumOfSquaredDeviations += delta * (value - mean); + } + + /** + * Incorporates every given sample, in order. + * + * @param values the samples to add + * @throws IllegalArgumentException if any value is NaN or infinite + * @throws NullPointerException if {@code values} is {@code null} + */ + public void addAll(double... values) { + for (double value : values) { + add(value); + } + } + + /** + * Removes a previously added sample, reversing {@link #add(double)}. This is what makes the + * accumulator usable for a sliding window. + * + * @param value the sample to remove; it must genuinely have been added before + * @throws IllegalStateException if the accumulator is empty + * @throws IllegalArgumentException if {@code value} is NaN or infinite + */ + public void remove(double value) { + requireFinite(value); + if (count == 0) { + throw new IllegalStateException("Cannot remove a sample from an empty accumulator"); + } + if (count == 1) { + clear(); + return; + } + double previousMean = mean; + mean = (count * mean - value) / (count - 1); + sumOfSquaredDeviations -= (value - previousMean) * (value - mean); + count--; + if (sumOfSquaredDeviations < 0.0) { + sumOfSquaredDeviations = 0.0; + } + } + + /** + * Combines two independently accumulated summaries using Chan's parallel variance update. + * + * @param left summary of the first batch of samples + * @param right summary of the second batch of samples + * @return a new summary describing the concatenation of both batches + * @throws NullPointerException if either argument is {@code null} + */ + public static WelfordAlgorithm merge(WelfordAlgorithm left, WelfordAlgorithm right) { + WelfordAlgorithm merged = new WelfordAlgorithm(); + merged.count = left.count + right.count; + if (merged.count == 0) { + return merged; + } + double delta = right.mean - left.mean; + merged.mean = left.mean + delta * right.count / merged.count; + merged.sumOfSquaredDeviations = left.sumOfSquaredDeviations + right.sumOfSquaredDeviations + delta * delta * left.count * right.count / merged.count; + return merged; + } + + /** + * Returns the number of samples seen so far. + * + * @return the sample count + */ + public long count() { + return count; + } + + /** + * Tells whether any sample has been added. + * + * @return {@code true} if no sample is currently accounted for + */ + public boolean isEmpty() { + return count == 0; + } + + /** + * Returns the arithmetic mean of the samples. + * + * @return the mean, or {@link Double#NaN} if no sample has been added + */ + public double mean() { + return count == 0 ? Double.NaN : mean; + } + + /** + * Returns the sum of the samples, reconstructed from the mean. + * + * @return {@code count * mean}, or {@code 0} if no sample has been added + */ + public double sum() { + return count == 0 ? 0.0 : mean * count; + } + + /** + * Returns the sum of squared deviations from the mean, {@code M2}. + * + * @return the sum of squared deviations, {@code 0} for an empty accumulator + */ + public double sumOfSquaredDeviations() { + return sumOfSquaredDeviations; + } + + /** + * Returns the unbiased sample variance, normalised by {@code count - 1}. + * + * @return the sample variance, or {@link Double#NaN} if fewer than two samples were added + */ + public double sampleVariance() { + return count < 2 ? Double.NaN : sumOfSquaredDeviations / (count - 1); + } + + /** + * Returns the population variance, normalised by {@code count}. + * + * @return the population variance, or {@link Double#NaN} if no sample has been added + */ + public double populationVariance() { + return count == 0 ? Double.NaN : sumOfSquaredDeviations / count; + } + + /** + * Returns the square root of {@link #sampleVariance()}. + * + * @return the sample standard deviation, or {@link Double#NaN} if fewer than two samples were added + */ + public double sampleStandardDeviation() { + return Math.sqrt(sampleVariance()); + } + + /** + * Returns the square root of {@link #populationVariance()}. + * + * @return the population standard deviation, or {@link Double#NaN} if no sample has been added + */ + public double populationStandardDeviation() { + return Math.sqrt(populationVariance()); + } + + /** + * Returns the standard error of the mean, {@code sampleStandardDeviation / sqrt(count)}. + * + * @return the standard error, or {@link Double#NaN} if fewer than two samples were added + */ + public double standardError() { + return sampleStandardDeviation() / Math.sqrt(count); + } + + /** + * Forgets every sample. + */ + public void clear() { + count = 0; + mean = 0.0; + sumOfSquaredDeviations = 0.0; + } + + @Override + public String toString() { + return "WelfordAlgorithm{count=" + count + ", mean=" + mean() + ", sampleStandardDeviation=" + sampleStandardDeviation() + '}'; + } + + private static void requireFinite(double value) { + if (!Double.isFinite(value)) { + throw new IllegalArgumentException("Samples must be finite, but was " + value); + } + } +} diff --git a/src/test/java/com/thealgorithms/streaming/WelfordAlgorithmTest.java b/src/test/java/com/thealgorithms/streaming/WelfordAlgorithmTest.java new file mode 100644 index 000000000000..5031e9de6e62 --- /dev/null +++ b/src/test/java/com/thealgorithms/streaming/WelfordAlgorithmTest.java @@ -0,0 +1,243 @@ +package com.thealgorithms.streaming; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; + +import java.util.Random; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.ValueSource; + +class WelfordAlgorithmTest { + + private static final double[] TEXTBOOK_SAMPLE = {2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0}; + + private static double twoPassPopulationVariance(double... values) { + double mean = 0.0; + for (double value : values) { + mean += value; + } + mean /= values.length; + + double sumOfSquares = 0.0; + for (double value : values) { + sumOfSquares += (value - mean) * (value - mean); + } + return sumOfSquares / values.length; + } + + @Test + void emptyAccumulatorReportsNothing() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + assertTrue(stats.isEmpty()); + assertEquals(0L, stats.count()); + assertEquals(0.0, stats.sum()); + assertEquals(0.0, stats.sumOfSquaredDeviations()); + assertTrue(Double.isNaN(stats.mean())); + assertTrue(Double.isNaN(stats.populationVariance())); + assertTrue(Double.isNaN(stats.sampleVariance())); + assertTrue(Double.isNaN(stats.sampleStandardDeviation())); + assertTrue(Double.isNaN(stats.standardError())); + } + + @Test + void singleSampleHasNoSpread() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.add(42.0); + assertFalse(stats.isEmpty()); + assertEquals(1L, stats.count()); + assertEquals(42.0, stats.mean()); + assertEquals(42.0, stats.sum()); + assertEquals(0.0, stats.populationVariance()); + assertTrue(Double.isNaN(stats.sampleVariance())); + } + + @Test + @DisplayName("reproduces the textbook values of a known sample") + void matchesKnownValues() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.addAll(TEXTBOOK_SAMPLE); + + assertEquals(8L, stats.count()); + assertEquals(5.0, stats.mean(), 1e-12); + assertEquals(40.0, stats.sum(), 1e-12); + assertEquals(32.0, stats.sumOfSquaredDeviations(), 1e-12); + assertEquals(4.0, stats.populationVariance(), 1e-12); + assertEquals(2.0, stats.populationStandardDeviation(), 1e-12); + assertEquals(32.0 / 7.0, stats.sampleVariance(), 1e-12); + assertEquals(Math.sqrt(32.0 / 7.0), stats.sampleStandardDeviation(), 1e-12); + assertEquals(Math.sqrt(32.0 / 7.0 / 8.0), stats.standardError(), 1e-12); + } + + @Test + @DisplayName("stays accurate where the naive sum of squares collapses") + void survivesLargeOffsets() { + double offset = 1e9; + double[] values = {offset + 4.0, offset + 7.0, offset + 13.0, offset + 16.0}; + + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.addAll(values); + + // The naive formula computes sum(x^2) - n * mean^2 on numbers of order 1e18 and loses every + // significant digit of the answer; Welford's recurrence never forms them. + double naive = 0.0; + for (double value : values) { + naive += value * value; + } + naive = (naive - values.length * stats.mean() * stats.mean()) / (values.length - 1); + + assertEquals(30.0, stats.sampleVariance(), 1e-6); + assertTrue(Math.abs(naive - 30.0) > 1.0, "the naive formula is expected to be far off, but returned " + naive); + } + + @ParameterizedTest + @ValueSource(ints = {2, 3, 10, 1_000}) + void agreesWithATwoPassComputation(int sampleCount) { + Random random = new Random(4242L + sampleCount); + double[] values = new double[sampleCount]; + for (int i = 0; i < sampleCount; i++) { + values[i] = random.nextGaussian() * 17.0 + 3.0; + } + + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.addAll(values); + + double expectedVariance = twoPassPopulationVariance(values); + assertEquals(expectedVariance, stats.populationVariance(), 1e-9 * Math.max(1.0, expectedVariance)); + } + + @Test + void removeUndoesAdd() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.addAll(TEXTBOOK_SAMPLE); + stats.add(1000.0); + stats.remove(1000.0); + + assertEquals(8L, stats.count()); + assertEquals(5.0, stats.mean(), 1e-9); + assertEquals(4.0, stats.populationVariance(), 1e-9); + } + + @Test + void removingTheLastSampleEmptiesTheAccumulator() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.add(3.0); + stats.remove(3.0); + assertTrue(stats.isEmpty()); + assertEquals(0.0, stats.sumOfSquaredDeviations()); + } + + @Test + void removeOnAnEmptyAccumulatorFails() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + assertThrows(IllegalStateException.class, () -> stats.remove(1.0)); + } + + @Test + @DisplayName("add plus remove turns the accumulator into a sliding window") + void supportsSlidingWindows() { + int windowSize = 20; + Random random = new Random(1234L); + double[] signal = new double[500]; + for (int i = 0; i < signal.length; i++) { + signal[i] = random.nextGaussian() * 5.0 + 100.0; + } + + WelfordAlgorithm window = new WelfordAlgorithm(); + for (int i = 0; i < signal.length; i++) { + window.add(signal[i]); + if (window.count() > windowSize) { + window.remove(signal[i - windowSize]); + } + if (i >= windowSize) { + double[] expectedWindow = new double[windowSize]; + System.arraycopy(signal, i - windowSize + 1, expectedWindow, 0, windowSize); + + WelfordAlgorithm reference = new WelfordAlgorithm(); + reference.addAll(expectedWindow); + assertEquals(windowSize, window.count()); + assertEquals(reference.mean(), window.mean(), 1e-9); + assertEquals(reference.populationVariance(), window.populationVariance(), 1e-8); + } + } + } + + @Test + void mergeCombinesPartialSummaries() { + Random random = new Random(20240517L); + double[] values = new double[1_000]; + for (int i = 0; i < values.length; i++) { + values[i] = random.nextGaussian(); + } + + WelfordAlgorithm whole = new WelfordAlgorithm(); + WelfordAlgorithm left = new WelfordAlgorithm(); + WelfordAlgorithm right = new WelfordAlgorithm(); + for (int i = 0; i < values.length; i++) { + whole.add(values[i]); + if (i < 337) { + left.add(values[i]); + } else { + right.add(values[i]); + } + } + + WelfordAlgorithm merged = WelfordAlgorithm.merge(left, right); + assertEquals(whole.count(), merged.count()); + assertEquals(whole.mean(), merged.mean(), 1e-12); + assertEquals(whole.sampleVariance(), merged.sampleVariance(), 1e-12); + } + + @Test + void mergeHandlesEmptyOperands() { + WelfordAlgorithm empty = new WelfordAlgorithm(); + WelfordAlgorithm filled = new WelfordAlgorithm(); + filled.addAll(TEXTBOOK_SAMPLE); + + assertTrue(WelfordAlgorithm.merge(empty, empty).isEmpty()); + assertEquals(filled.mean(), WelfordAlgorithm.merge(empty, filled).mean(), 1e-12); + assertEquals(filled.mean(), WelfordAlgorithm.merge(filled, empty).mean(), 1e-12); + assertEquals(filled.sampleVariance(), WelfordAlgorithm.merge(filled, empty).sampleVariance(), 1e-12); + } + + @Test + void clearForgetsEverything() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.addAll(TEXTBOOK_SAMPLE); + stats.clear(); + assertTrue(stats.isEmpty()); + assertTrue(Double.isNaN(stats.mean())); + } + + @ParameterizedTest + @ValueSource(doubles = {Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}) + void rejectsNonFiniteSamples(double value) { + WelfordAlgorithm stats = new WelfordAlgorithm(); + assertThrows(IllegalArgumentException.class, () -> stats.add(value)); + stats.add(1.0); + assertThrows(IllegalArgumentException.class, () -> stats.remove(value)); + } + + @Test + void toStringMentionsTheSummary() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.addAll(TEXTBOOK_SAMPLE); + assertTrue(stats.toString().contains("count=8"), stats.toString()); + } + + @Test + @DisplayName("removing a sample far from the mean cannot push the variance below zero") + void removalNeverYieldsANegativeVariance() { + WelfordAlgorithm stats = new WelfordAlgorithm(); + stats.add(1e16); + stats.add(1.0); + stats.remove(1e16); + + assertEquals(1L, stats.count()); + assertTrue(stats.sumOfSquaredDeviations() >= 0.0, "got " + stats.sumOfSquaredDeviations()); + assertTrue(stats.populationVariance() >= 0.0); + } +}