Grubbs' test asks whether the single most extreme point is an outlier, and it breaks as soon as + * there are two of them: each one inflates the standard deviation the other is measured against, so + * a pair of outliers can hide one another completely. This is the masking problem, and the + * generalized test is the standard answer to it. It removes the most extreme point, recomputes the + * statistics on what is left, and repeats {@code r} times: + * + *
+ * R_i = max |x - mean| / sd over the points that are still in + * lambda_i = (n - i) * t / sqrt( (n - i - 1 + t^2)(n - i + 1) ) + * where t = the 1 - alpha / (2(n - i + 1)) quantile of a t distribution with n - i - 1 degrees of freedom + *+ * + *
The number of outliers is the largest {@code i} whose statistic exceeds its critical + * value, not the first one. That is what defeats masking: in a sample with three outliers the first + * two tests can easily fall short while the third one succeeds, and the test then reports all three. + * + *
This is a batch test and it needs the whole sample, which is exactly the trade it makes against + * the streaming detectors in this package: it has a stated significance level and it decides how many + * points are outliers, where {@link HampelFilter} answers one sample at a time against a threshold the + * caller has to choose. Rosner's own recommendation is to use it on at least 25 points, and {@code r} + * is an upper bound that may be set generously, because overstating it costs accuracy only in extreme + * cases. + * + *
{@code
+ * GeneralizedEsd test = new GeneralizedEsd(0.05);
+ * int[] outliers = test.findOutliers(sample, 10);
+ * }
+ *
+ * The test costs O(r * n) time and allocates O(n). The t quantile behind the critical value is + * computed from the regularized incomplete beta function, so no table is needed. + * + * @see HampelFilter + * @see Grubbs's test and its generalization + */ +public final class GeneralizedEsd { + + /** Significance level used when none is given. */ + public static final double DEFAULT_SIGNIFICANCE_LEVEL = 0.05; + + private static final int MINIMUM_SAMPLE_SIZE = 3; + private static final double[] LANCZOS = { + 0.99999999999980993, + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7, + }; + + private final double significanceLevel; + + /** + * Creates a test at the customary significance level of {@code 0.05}. + */ + public GeneralizedEsd() { + this(DEFAULT_SIGNIFICANCE_LEVEL); + } + + /** + * Creates a test. + * + * @param significanceLevel the probability of declaring an outlier in clean data, in {@code (0, 1)} + * @throws IllegalArgumentException if {@code significanceLevel} is not inside {@code (0, 1)} + */ + public GeneralizedEsd(double significanceLevel) { + if (!(significanceLevel > 0.0) || !(significanceLevel < 1.0)) { + throw new IllegalArgumentException("The significance level must lie in (0, 1), but was " + significanceLevel); + } + this.significanceLevel = significanceLevel; + } + + /** + * Finds the outliers of a sample. + * + * @param sample the data to inspect, left untouched + * @param maxOutliers upper bound on the number of outliers, at least one and at most + * {@code sample.length - 2} + * @return the indices of the outliers in ascending order, empty if the sample looks clean + * @throws IllegalArgumentException if the sample holds fewer than three points or a non-finite + * value, or if {@code maxOutliers} is out of range + * @throws NullPointerException if {@code sample} is {@code null} + */ + public int[] findOutliers(double[] sample, int maxOutliers) { + if (sample.length < MINIMUM_SAMPLE_SIZE) { + throw new IllegalArgumentException("The sample must hold at least " + MINIMUM_SAMPLE_SIZE + " points, but held " + sample.length); + } + if (maxOutliers < 1 || maxOutliers > sample.length - 2) { + throw new IllegalArgumentException("The number of outliers to look for must lie in [1, " + (sample.length - 2) + "], but was " + maxOutliers); + } + for (double value : sample) { + if (!Double.isFinite(value)) { + throw new IllegalArgumentException("Samples must be finite, but was " + value); + } + } + + int size = sample.length; + double[] values = Arrays.copyOf(sample, size); + int[] indices = new int[size]; + for (int i = 0; i < size; i++) { + indices[i] = i; + } + + int[] removed = new int[maxOutliers]; + int outliers = 0; + int remaining = size; + + for (int step = 1; step <= maxOutliers; step++) { + double mean = mean(values, remaining); + double deviation = standardDeviation(values, remaining, mean); + if (deviation == 0.0) { + break; + } + + int worst = indexOfLargestDeviation(values, remaining, mean); + double statistic = Math.abs(values[worst] - mean) / deviation; + removed[step - 1] = indices[worst]; + if (statistic > criticalValue(size, step)) { + outliers = step; + } + + remaining--; + values[worst] = values[remaining]; + indices[worst] = indices[remaining]; + } + + int[] result = Arrays.copyOf(removed, outliers); + Arrays.sort(result); + return result; + } + + /** + * Returns the critical value a statistic has to exceed at one step of the test. + * + * @param sampleSize the size of the whole sample + * @param step which removal this is, counting from one + * @return the critical value, that is {@code lambda_i} + * @throws IllegalArgumentException if {@code sampleSize} is smaller than three, or if {@code step} + * is outside {@code [1, sampleSize - 2]} + */ + public double criticalValue(int sampleSize, int step) { + if (sampleSize < MINIMUM_SAMPLE_SIZE) { + throw new IllegalArgumentException("The sample size must be at least " + MINIMUM_SAMPLE_SIZE + ", but was " + sampleSize); + } + if (step < 1 || step > sampleSize - 2) { + throw new IllegalArgumentException("The step must lie in [1, " + (sampleSize - 2) + "], but was " + step); + } + int size = sampleSize - step + 1; + int degreesOfFreedom = size - 2; + double probability = 1.0 - significanceLevel / (2.0 * size); + double t = studentTQuantile(probability, degreesOfFreedom); + return (size - 1) * t / Math.sqrt((degreesOfFreedom + t * t) * size); + } + + /** + * Returns the significance level of the test. + * + * @return the level given at construction time + */ + public double significanceLevel() { + return significanceLevel; + } + + @Override + public String toString() { + return "GeneralizedEsd{significanceLevel=" + significanceLevel + "}"; + } + + private static double mean(double[] values, int size) { + double sum = 0.0; + for (int i = 0; i < size; i++) { + sum += values[i]; + } + return sum / size; + } + + private static double standardDeviation(double[] values, int size, double mean) { + double sum = 0.0; + for (int i = 0; i < size; i++) { + double deviation = values[i] - mean; + sum += deviation * deviation; + } + return Math.sqrt(sum / (size - 1)); + } + + private static int indexOfLargestDeviation(double[] values, int size, double mean) { + int worst = 0; + double largest = -1.0; + for (int i = 0; i < size; i++) { + double deviation = Math.abs(values[i] - mean); + if (deviation > largest) { + largest = deviation; + worst = i; + } + } + return worst; + } + + /** + * Returns the {@code probability} quantile of a t distribution, found by bisecting its cumulative + * distribution function. + * + * @param probability the probability to invert, in {@code (0, 1)} + * @param degreesOfFreedom the degrees of freedom, at least one + * @return the quantile + */ + private static double studentTQuantile(double probability, int degreesOfFreedom) { + if (probability < 0.5) { + return -studentTQuantile(1.0 - probability, degreesOfFreedom); + } + double low = 0.0; + double high = 1.0; + while (studentTCumulative(high, degreesOfFreedom) < probability && high < 1e12) { + high *= 2.0; + } + for (int i = 0; i < 200; i++) { + double middle = 0.5 * (low + high); + if (studentTCumulative(middle, degreesOfFreedom) < probability) { + low = middle; + } else { + high = middle; + } + } + return 0.5 * (low + high); + } + + /** + * Returns the probability that a t distributed variable stays below {@code t}. + * + * @param t the point to evaluate at + * @param degreesOfFreedom the degrees of freedom, at least one + * @return the cumulative probability + */ + private static double studentTCumulative(double t, int degreesOfFreedom) { + double x = degreesOfFreedom / (degreesOfFreedom + t * t); + double tail = 0.5 * regularizedIncompleteBeta(0.5 * degreesOfFreedom, 0.5, x); + return t >= 0.0 ? 1.0 - tail : tail; + } + + /** + * Returns the regularized incomplete beta function, evaluated with the continued fraction of + * Lentz. + * + * @param a first shape parameter, strictly positive + * @param b second shape parameter, strictly positive + * @param x the point to evaluate at, in {@code [0, 1]} + * @return {@code I_x(a, b)} + */ + private static double regularizedIncompleteBeta(double a, double b, double x) { + if (x <= 0.0) { + return 0.0; + } + if (x >= 1.0) { + return 1.0; + } + double logBeta = logGamma(a) + logGamma(b) - logGamma(a + b); + if (x < (a + 1.0) / (a + b + 2.0)) { + return Math.exp(a * Math.log(x) + b * Math.log1p(-x) - logBeta) * betaContinuedFraction(a, b, x) / a; + } + return 1.0 - Math.exp(b * Math.log1p(-x) + a * Math.log(x) - logBeta) * betaContinuedFraction(b, a, 1.0 - x) / b; + } + + private static double betaContinuedFraction(double a, double b, double x) { + double tiny = 1e-300; + double c = 1.0; + double d = 1.0 - (a + b) * x / (a + 1.0); + if (Math.abs(d) < tiny) { + d = tiny; + } + d = 1.0 / d; + double fraction = d; + + for (int m = 1; m <= 300; m++) { + int even = 2 * m; + double numerator = m * (b - m) * x / ((a + even - 1.0) * (a + even)); + d = 1.0 + numerator * d; + if (Math.abs(d) < tiny) { + d = tiny; + } + c = 1.0 + numerator / c; + if (Math.abs(c) < tiny) { + c = tiny; + } + d = 1.0 / d; + fraction *= d * c; + + numerator = -(a + m) * (a + b + m) * x / ((a + even) * (a + even + 1.0)); + d = 1.0 + numerator * d; + if (Math.abs(d) < tiny) { + d = tiny; + } + c = 1.0 + numerator / c; + if (Math.abs(c) < tiny) { + c = tiny; + } + d = 1.0 / d; + double step = d * c; + fraction *= step; + + if (Math.abs(step - 1.0) < 1e-15) { + break; + } + } + return fraction; + } + + /** + * Returns the logarithm of the gamma function, using the Lanczos approximation. + * + * @param x the point to evaluate at, strictly positive + * @return {@code log(gamma(x))} + */ + private static double logGamma(double x) { + double shifted = x - 1.0; + double series = LANCZOS[0]; + for (int i = 1; i < LANCZOS.length; i++) { + series += LANCZOS[i] / (shifted + i); + } + double t = shifted + 7.5; + return 0.5 * Math.log(2.0 * Math.PI) + (shifted + 0.5) * Math.log(t) - t + Math.log(series); + } +} diff --git a/src/test/java/com/thealgorithms/streaming/GeneralizedEsdTest.java b/src/test/java/com/thealgorithms/streaming/GeneralizedEsdTest.java new file mode 100644 index 000000000000..cd5e1db492d9 --- /dev/null +++ b/src/test/java/com/thealgorithms/streaming/GeneralizedEsdTest.java @@ -0,0 +1,212 @@ +package com.thealgorithms.streaming; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; + +import java.util.Arrays; +import java.util.Random; +import org.junit.jupiter.api.Assertions; +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 GeneralizedEsdTest { + + /** + * The sample Rosner published with the test, also used as the worked example in the NIST + * engineering statistics handbook. + */ + private static final double[] ROSNER_SAMPLE = {-0.25, 0.68, 0.94, 1.15, 1.20, 1.26, 1.26, 1.34, 1.38, 1.43, 1.49, 1.49, 1.55, 1.56, 1.58, 1.65, 1.69, 1.70, 1.76, 1.77, 1.81, 1.91, 1.94, 1.96, 1.99, 2.06, 2.09, 2.10, 2.14, 2.15, 2.23, 2.24, 2.26, 2.35, 2.37, 2.40, 2.47, 2.54, 2.62, 2.64, 2.90, + 2.92, 2.92, 2.93, 3.21, 3.26, 3.30, 3.59, 3.68, 4.30, 4.64, 5.34, 5.42, 6.01}; + + private static double[] gaussianSample(int size, long seed) { + Random random = new Random(seed); + double[] sample = new double[size]; + for (int i = 0; i < size; i++) { + sample[i] = 100.0 + 5.0 * random.nextGaussian(); + } + return sample; + } + + @ParameterizedTest + @ValueSource(doubles = {0.0, 1.0, -0.1, 1.5, Double.NaN}) + void rejectsInvalidSignificanceLevels(double level) { + assertThrows(IllegalArgumentException.class, () -> new GeneralizedEsd(level)); + } + + @Test + void rejectsSamplesThatAreTooSmall() { + GeneralizedEsd test = new GeneralizedEsd(); + + assertThrows(IllegalArgumentException.class, () -> test.findOutliers(new double[] {1.0, 2.0}, 1)); + } + + @Test + void rejectsAnOutlierBoundOutOfRange() { + GeneralizedEsd test = new GeneralizedEsd(); + double[] sample = gaussianSample(20, 1L); + + assertThrows(IllegalArgumentException.class, () -> test.findOutliers(sample, 0)); + assertThrows(IllegalArgumentException.class, () -> test.findOutliers(sample, 19)); + } + + @ParameterizedTest + @ValueSource(doubles = {Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}) + void rejectsNonFiniteSamples(double value) { + GeneralizedEsd test = new GeneralizedEsd(); + double[] sample = gaussianSample(20, 2L); + sample[7] = value; + + assertThrows(IllegalArgumentException.class, () -> test.findOutliers(sample, 3)); + } + + @Test + void exposesItsConfiguration() { + assertEquals(0.05, new GeneralizedEsd().significanceLevel()); + assertEquals(0.01, new GeneralizedEsd(0.01).significanceLevel()); + } + + @Test + @DisplayName("the critical values match the table published with the test") + void reproducesThePublishedCriticalValues() { + GeneralizedEsd test = new GeneralizedEsd(0.05); + + assertEquals(3.159, test.criticalValue(54, 1), 5e-4); + assertEquals(3.151, test.criticalValue(54, 2), 5e-4); + assertEquals(3.144, test.criticalValue(54, 3), 5e-4); + assertEquals(3.085, test.criticalValue(54, 10), 5e-4); + } + + @Test + @DisplayName("Rosner's own sample of 54 points holds exactly three outliers") + void reproducesTheWorkedExample() { + GeneralizedEsd test = new GeneralizedEsd(0.05); + + int[] outliers = test.findOutliers(ROSNER_SAMPLE, 10); + + assertArrayEqualsSorted(new int[] {51, 52, 53}, outliers); + assertEquals(5.34, ROSNER_SAMPLE[outliers[0]]); + assertEquals(5.42, ROSNER_SAMPLE[outliers[1]]); + assertEquals(6.01, ROSNER_SAMPLE[outliers[2]]); + } + + @Test + @DisplayName("the answer is the last test that succeeds, not the first") + void findsOutliersThatMaskEachOther() { + double[] sample = gaussianSample(40, 0L); + sample[10] = 160.0; + sample[20] = 161.0; + sample[30] = 162.0; + + int[] outliers = new GeneralizedEsd(0.05).findOutliers(sample, 5); + + assertArrayEqualsSorted(new int[] {10, 20, 30}, outliers); + } + + @Test + void findsASingleOutlier() { + double[] sample = gaussianSample(30, 4L); + sample[13] = 200.0; + + int[] outliers = new GeneralizedEsd().findOutliers(sample, 5); + + assertArrayEqualsSorted(new int[] {13}, outliers); + } + + @Test + @DisplayName("clean data holds no outliers") + void findsNothingInCleanData() { + for (long seed = 0; seed < 12; seed++) { + double[] sample = gaussianSample(60, seed); + + int[] outliers = new GeneralizedEsd(0.001).findOutliers(sample, 6); + + assertEquals(0, outliers.length, "seed " + seed + " reported " + Arrays.toString(outliers)); + } + } + + @Test + @DisplayName("the significance level is a rate, not a promise: clean data does report the odd outlier") + void falseAlarmsFollowTheSignificanceLevel() { + int reporting = 0; + for (long seed = 0; seed < 12; seed++) { + if (new GeneralizedEsd(0.05).findOutliers(gaussianSample(60, seed), 6).length > 0) { + reporting++; + } + } + + assertTrue(reporting > 0, "at five percent some clean sample is expected to report something"); + assertTrue(reporting <= 4, "but it must stay rare, " + reporting + " of 12 samples reported"); + } + + @Test + @DisplayName("a sample with no spread at all has no outliers") + void handlesAConstantSample() { + double[] sample = new double[20]; + Arrays.fill(sample, 3.0); + + assertEquals(0, new GeneralizedEsd().findOutliers(sample, 5).length); + } + + @Test + void leavesTheSampleUntouched() { + double[] sample = gaussianSample(30, 5L); + sample[3] = 500.0; + double[] copy = Arrays.copyOf(sample, sample.length); + + new GeneralizedEsd().findOutliers(sample, 4); + + Assertions.assertArrayEquals(copy, sample); + } + + @Test + @DisplayName("a stricter level reports at most as many outliers as a looser one") + void aStricterLevelIsMoreConservative() { + double[] sample = gaussianSample(50, 6L); + sample[5] = 125.0; + sample[6] = 126.0; + + int loose = new GeneralizedEsd(0.10).findOutliers(sample, 8).length; + int strict = new GeneralizedEsd(0.001).findOutliers(sample, 8).length; + + assertTrue(strict <= loose, "strict found " + strict + ", loose found " + loose); + } + + @Test + @DisplayName("the bound on the number of outliers is respected") + void neverReportsMoreThanTheBound() { + double[] sample = gaussianSample(40, 7L); + for (int i = 0; i < 6; i++) { + sample[i] = 300.0 + i; + } + + int[] outliers = new GeneralizedEsd().findOutliers(sample, 2); + + assertTrue(outliers.length <= 2); + } + + @Test + void findsOutliersOnBothSides() { + double[] sample = gaussianSample(40, 8L); + sample[4] = 300.0; + sample[24] = -150.0; + + int[] outliers = new GeneralizedEsd().findOutliers(sample, 4); + + assertArrayEqualsSorted(new int[] {4, 24}, outliers); + } + + @Test + void toStringMentionsTheLevel() { + String text = new GeneralizedEsd(0.05).toString(); + + assertTrue(text.contains("GeneralizedEsd")); + assertTrue(text.contains("0.05")); + } + + private static void assertArrayEqualsSorted(int[] expected, int[] actual) { + Assertions.assertArrayEquals(expected, actual, "expected " + Arrays.toString(expected) + " but got " + Arrays.toString(actual)); + } +}