feat: add GeneralizedEsd, Rosner's test for up to r outliers (#7612)

Grubbs' test asks whether the single most extreme point is an outlier and breaks as soon as there are two of them, because each one inflates the standard deviation the other is measured against. The generalized extreme Studentized deviate test removes the most extreme point, recomputes the statistics on what is left and repeats r times, and the answer is the largest step whose statistic exceeds its critical value rather than the first, which is what defeats that masking.

The critical values need a Student t quantile, which the repository did not have, so the class computes it from the regularized incomplete beta function through a Lentz continued fraction and a Lanczos log gamma, inverted by bisection. No table is needed. The tests check the critical values and the number of outliers against Rosner's published 54 point example.

Signed-off-by: alxkm <19151554+alxkm@users.noreply.github.com>
Co-authored-by: alxkm <19151554+alxkm@users.noreply.github.com>

Author: alxkm

src/main/java/com/thealgorithms/streaming/GeneralizedEsd.java

@@ -0,0 +1,334 @@
+package com.thealgorithms.streaming;
+
+import java.util.Arrays;
+
+/**
+ * The <b>generalized extreme Studentized deviate test</b> of Rosner: finds up to {@code r} outliers
+ * in a sample, without being told how many are there.
+ *
+ * <p>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:
+ *
+ * <pre>
+ * 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
+ * </pre>
+ *
+ * <p>The number of outliers is the <i>largest</i> {@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.
+ *
+ * <p>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.
+ *
+ * <h2>Usage</h2>
+ *
+ * <pre>{@code
+ * GeneralizedEsd test = new GeneralizedEsd(0.05);
+ * int[] outliers = test.findOutliers(sample, 10);
+ * }</pre>
+ *
+ * <p>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 <a href="https://en.wikipedia.org/wiki/Grubbs%27s_test">Grubbs's test and its generalization</a>
+ */
+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);
+    }
+}