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); + } +}