Computing an exact quantile requires storing everything seen so far, which is impossible for an + * unbounded stream. P² instead keeps five markers that track the minimum, the quantile + * {@code p/2}, the quantile {@code p} itself, the quantile {@code (1 + p)/2} and the maximum. Each + * marker has a height (its current value) and a position (how many samples are known to be below + * it). Every new sample shifts the positions of the markers it falls below, and each marker is then + * nudged back towards the position it is supposed to occupy. The nudge follows a piecewise + * parabolic prediction fitted through the three neighbouring markers, falling back to a linear one + * whenever the parabola would break the ordering of the heights. + * + *
Accuracy is typically a fraction of a percent in rank for stationary streams and improves as + * more samples arrive. The markers adapt continuously, so a drifting distribution is tracked rather + * than averaged away, but no accuracy bound is guaranteed for adversarial inputs. + * + *
| Operation | Complexity |
|---|---|
| {@link #add(double)} | O(1) |
| {@link #quantile()} | O(1) |
| memory | O(1) - five markers, independent of the stream length |
{@code
+ * P2QuantileEstimator p99 = new P2QuantileEstimator(0.99);
+ * for (double latency : latencies) {
+ * p99.add(latency);
+ * }
+ * double tail = p99.quantile();
+ * }
+ *
+ * The first five samples are stored verbatim, so the estimate is exact until the sixth one + * arrives. This class is not thread-safe. + * + * @see R. Jain, I. Chlamtac, The P-square algorithm (1985) + */ +public final class P2QuantileEstimator { + + private static final int MARKERS = 5; + + private final double probability; + + /** Marker heights, kept in non-decreasing order; also the raw buffer for the first five samples. */ + private final double[] heights = new double[MARKERS]; + + /** One-based rank of each marker, i.e. how many samples are known to sit at or below it. */ + private final int[] positions = new int[MARKERS]; + + /** Position each marker should ideally occupy. */ + private final double[] desiredPositions = new double[MARKERS]; + + /** How much {@link #desiredPositions} grows per sample. */ + private final double[] positionIncrements = new double[MARKERS]; + + private long count; + + /** + * Creates an estimator for a single quantile. + * + * @param probability the quantile to track, strictly between 0 and 1, e.g. {@code 0.5} for the median + * @throws IllegalArgumentException if {@code probability} is outside {@code (0, 1)} + */ + public P2QuantileEstimator(double probability) { + if (!(probability > 0.0) || !(probability < 1.0)) { + throw new IllegalArgumentException("The quantile probability must lie strictly between 0 and 1, but was " + probability); + } + this.probability = probability; + } + + /** + * Incorporates one sample. + * + * @param value the sample to add + * @throws IllegalArgumentException if {@code value} is NaN or infinite + */ + public void add(double value) { + if (!Double.isFinite(value)) { + throw new IllegalArgumentException("Samples must be finite, but was " + value); + } + if (count < MARKERS) { + collectInitialSample(value); + return; + } + + int cell = locate(value); + for (int i = cell + 1; i < MARKERS; i++) { + positions[i]++; + } + for (int i = 0; i < MARKERS; i++) { + desiredPositions[i] += positionIncrements[i]; + } + for (int i = 1; i < MARKERS - 1; i++) { + adjustMarker(i); + } + count++; + } + + /** + * 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); + } + } + + /** + * Returns the current estimate of the tracked quantile. + * + * @return the estimated quantile, exact while fewer than five samples have been seen + * @throws IllegalStateException if no sample has been added yet + */ + public double quantile() { + requireNonEmpty(); + if (count < MARKERS) { + return exactQuantileOfInitialSamples(); + } + return heights[2]; + } + + /** + * Returns the smallest sample seen so far, which P² tracks exactly. + * + * @return the running minimum + * @throws IllegalStateException if no sample has been added yet + */ + public double min() { + requireNonEmpty(); + return count < MARKERS ? Arrays.stream(heights, 0, (int) count).min().orElse(Double.NaN) : heights[0]; + } + + /** + * Returns the largest sample seen so far, which P² tracks exactly. + * + * @return the running maximum + * @throws IllegalStateException if no sample has been added yet + */ + public double max() { + requireNonEmpty(); + return count < MARKERS ? Arrays.stream(heights, 0, (int) count).max().orElse(Double.NaN) : heights[MARKERS - 1]; + } + + /** + * Returns the quantile this estimator was configured for. + * + * @return the probability given at construction time + */ + public double probability() { + return probability; + } + + /** + * 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 the estimator holds no samples + */ + public boolean isEmpty() { + return count == 0; + } + + /** + * Forgets every sample. + */ + public void reset() { + count = 0; + Arrays.fill(heights, 0.0); + Arrays.fill(positions, 0); + Arrays.fill(desiredPositions, 0.0); + Arrays.fill(positionIncrements, 0.0); + } + + @Override + public String toString() { + return "P2QuantileEstimator{p=" + probability + ", count=" + count + ", quantile=" + (count == 0 ? Double.NaN : quantile()) + '}'; + } + + /** + * Stores one of the first five samples and, on the fifth, lays out the markers. + */ + private void collectInitialSample(double value) { + heights[(int) count] = value; + count++; + if (count < MARKERS) { + return; + } + Arrays.sort(heights); + // Marker positions are one-based, exactly as in the paper. Only differences of positions are + // ever used, so the origin is arithmetically irrelevant - but the desired positions are grown + // by repeated addition, and starting them one lower would round differently and flip the + // strict comparisons below on the samples where a desired position lands on a whole number. + for (int i = 0; i < MARKERS; i++) { + positions[i] = i + 1; + } + desiredPositions[0] = 1.0; + desiredPositions[1] = 1.0 + 2.0 * probability; + desiredPositions[2] = 1.0 + 4.0 * probability; + desiredPositions[3] = 3.0 + 2.0 * probability; + desiredPositions[4] = 5.0; + positionIncrements[0] = 0.0; + positionIncrements[1] = probability / 2.0; + positionIncrements[2] = probability; + positionIncrements[3] = (1.0 + probability) / 2.0; + positionIncrements[4] = 1.0; + } + + /** + * Finds the cell the sample falls into, stretching the outer markers when it falls outside the + * range seen so far. + * + * @param value the incoming sample + * @return the index of the marker immediately below the sample + */ + private int locate(double value) { + if (value < heights[0]) { + heights[0] = value; + return 0; + } + for (int i = 1; i < MARKERS - 1; i++) { + if (value < heights[i]) { + return i - 1; + } + } + if (value > heights[MARKERS - 1]) { + heights[MARKERS - 1] = value; + } + return MARKERS - 2; + } + + /** + * Moves one inner marker by a single position if it has drifted too far from where it should be. + */ + private void adjustMarker(int i) { + double drift = desiredPositions[i] - positions[i]; + boolean shiftRight = drift >= 1.0 && positions[i + 1] - positions[i] > 1; + boolean shiftLeft = drift <= -1.0 && positions[i - 1] - positions[i] < -1; + if (!shiftRight && !shiftLeft) { + return; + } + + int direction = shiftRight ? 1 : -1; + double candidate = parabolicPrediction(i, direction); + heights[i] = heights[i - 1] < candidate && candidate < heights[i + 1] ? candidate : linearPrediction(i, direction); + positions[i] += direction; + } + + /** + * Piecewise parabolic prediction: fits a parabola through markers {@code i - 1}, {@code i} and + * {@code i + 1} and evaluates it one position away from the current one. + */ + private double parabolicPrediction(int i, int direction) { + double left = positions[i] - positions[i - 1]; + double right = positions[i + 1] - positions[i]; + // Grouped exactly as in the paper. The result then decides, by a strict comparison against the + // neighbouring heights, whether the parabola is used at all, so regrouping the arithmetic would + // flip that decision on knife-edge inputs and make this estimator drift away from the reference. + return heights[i] + direction / (left + right) * ((left + direction) * (heights[i + 1] - heights[i]) / right + (right - direction) * (heights[i] - heights[i - 1]) / left); + } + + /** + * Linear fallback used whenever the parabolic prediction would violate the ordering of the + * marker heights. + */ + private double linearPrediction(int i, int direction) { + int neighbour = i + direction; + return heights[i] + direction * (heights[neighbour] - heights[i]) / (positions[neighbour] - positions[i]); + } + + /** + * Exact quantile of the fewer than five samples buffered so far, by linear interpolation between + * the two order statistics surrounding the requested rank. + */ + private double exactQuantileOfInitialSamples() { + double[] sorted = Arrays.copyOf(heights, (int) count); + Arrays.sort(sorted); + if (sorted.length == 1) { + return sorted[0]; + } + double rank = probability * (sorted.length - 1); + int lower = (int) Math.floor(rank); + int upper = Math.min(lower + 1, sorted.length - 1); + return sorted[lower] + (rank - lower) * (sorted[upper] - sorted[lower]); + } + + private void requireNonEmpty() { + if (count == 0) { + throw new IllegalStateException("The estimator has not seen any sample yet"); + } + } +} diff --git a/src/test/java/com/thealgorithms/streaming/P2QuantileEstimatorTest.java b/src/test/java/com/thealgorithms/streaming/P2QuantileEstimatorTest.java new file mode 100644 index 000000000000..777e00d62d66 --- /dev/null +++ b/src/test/java/com/thealgorithms/streaming/P2QuantileEstimatorTest.java @@ -0,0 +1,206 @@ +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.Arrays; +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 P2QuantileEstimatorTest { + + /** + * Exact quantile of a sample, by linear interpolation between the surrounding order statistics. + */ + private static double exactQuantile(double[] sortedValues, double probability) { + double rank = probability * (sortedValues.length - 1); + int lower = (int) Math.floor(rank); + int upper = Math.min(lower + 1, sortedValues.length - 1); + return sortedValues[lower] + (rank - lower) * (sortedValues[upper] - sortedValues[lower]); + } + + private static double[] shuffledUniformSample(int count, long seed) { + Random random = new Random(seed); + double[] values = new double[count]; + for (int i = 0; i < count; i++) { + values[i] = random.nextDouble(); + } + return values; + } + + @ParameterizedTest + @ValueSource(doubles = {0.0, 1.0, -0.1, 1.5, Double.NaN}) + void rejectsInvalidProbabilities(double probability) { + assertThrows(IllegalArgumentException.class, () -> new P2QuantileEstimator(probability)); + } + + @Test + void queriesBeforeTheFirstSampleFail() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + assertTrue(estimator.isEmpty()); + assertEquals(0L, estimator.count()); + assertEquals(0.5, estimator.probability()); + assertThrows(IllegalStateException.class, estimator::quantile); + assertThrows(IllegalStateException.class, estimator::min); + assertThrows(IllegalStateException.class, estimator::max); + } + + @ParameterizedTest + @ValueSource(doubles = {Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}) + void rejectsNonFiniteSamples(double value) { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + assertThrows(IllegalArgumentException.class, () -> estimator.add(value)); + } + + @Test + @DisplayName("the first five samples are answered exactly") + void isExactDuringWarmUp() { + P2QuantileEstimator median = new P2QuantileEstimator(0.5); + median.add(3.0); + assertEquals(3.0, median.quantile()); + median.add(1.0); + assertEquals(2.0, median.quantile(), 1e-12); + median.add(2.0); + assertEquals(2.0, median.quantile(), 1e-12); + median.add(10.0); + assertEquals(2.5, median.quantile(), 1e-12); + assertEquals(1.0, median.min()); + assertEquals(10.0, median.max()); + assertEquals(4L, median.count()); + } + + @Test + void tracksTheExtremesExactly() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + double[] values = shuffledUniformSample(5_000, 7L); + estimator.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + assertEquals(sorted[0], estimator.min(), 0.0); + assertEquals(sorted[sorted.length - 1], estimator.max(), 0.0); + assertEquals(5_000L, estimator.count()); + } + + @Test + void aConstantStreamHasAConstantQuantile() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.9); + for (int i = 0; i < 1_000; i++) { + estimator.add(7.0); + } + assertEquals(7.0, estimator.quantile(), 1e-12); + } + + @ParameterizedTest + @ValueSource(doubles = {0.05, 0.25, 0.5, 0.75, 0.9, 0.99}) + @DisplayName("estimates a uniform stream within a percent of the exact quantile") + void approximatesUniformQuantiles(double probability) { + double[] values = shuffledUniformSample(100_000, 20240517L); + P2QuantileEstimator estimator = new P2QuantileEstimator(probability); + estimator.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + assertEquals(exactQuantile(sorted, probability), estimator.quantile(), 0.01); + } + + @Test + @DisplayName("estimates the tail of a skewed stream") + void approximatesTheTailOfAnExponentialStream() { + Random random = new Random(31337L); + double[] values = new double[100_000]; + for (int i = 0; i < values.length; i++) { + values[i] = -Math.log(1.0 - random.nextDouble()) * 10.0; + } + + P2QuantileEstimator estimator = new P2QuantileEstimator(0.95); + estimator.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + double expected = exactQuantile(sorted, 0.95); + assertEquals(expected, estimator.quantile(), 0.05 * expected); + } + + @Test + @DisplayName("stays sane even when the stream arrives already sorted, the worst case for P-square") + void handlesMonotonicInput() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + for (int i = 1; i <= 10_000; i++) { + estimator.add(i); + } + assertEquals(5_000.0, estimator.quantile(), 500.0); + assertEquals(1.0, estimator.min()); + assertEquals(10_000.0, estimator.max()); + } + + @Test + void estimatesOfDifferentQuantilesStayOrdered() { + double[] values = shuffledUniformSample(50_000, 99L); + P2QuantileEstimator low = new P2QuantileEstimator(0.25); + P2QuantileEstimator middle = new P2QuantileEstimator(0.5); + P2QuantileEstimator high = new P2QuantileEstimator(0.75); + low.addAll(values); + middle.addAll(values); + high.addAll(values); + + assertTrue(low.quantile() < middle.quantile()); + assertTrue(middle.quantile() < high.quantile()); + } + + @Test + void resetForgetsEverything() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + estimator.addAll(shuffledUniformSample(1_000, 5L)); + estimator.reset(); + + assertTrue(estimator.isEmpty()); + assertThrows(IllegalStateException.class, estimator::quantile); + + estimator.addAll(1.0, 2.0, 3.0); + assertEquals(2.0, estimator.quantile(), 1e-12); + } + + @Test + void toStringMentionsTheEstimate() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + assertTrue(estimator.toString().contains("count=0"), estimator.toString()); + estimator.addAll(1.0, 2.0, 3.0, 4.0, 5.0, 6.0); + assertFalse(estimator.isEmpty()); + assertTrue(estimator.toString().contains("quantile="), estimator.toString()); + } + + @Test + @DisplayName("violent jumps force the parabolic prediction to fall back to a linear one") + void survivesWildlyJumpingInput() { + P2QuantileEstimator estimator = new P2QuantileEstimator(0.5); + double value = 1.0; + for (int i = 0; i < 5_000; i++) { + value = i % 2 == 0 ? value * 3.0 + 1.0 : 1.0 / (i + 1); + estimator.add(value); + assertTrue(estimator.quantile() >= estimator.min(), "the estimate fell below the minimum"); + assertTrue(estimator.quantile() <= estimator.max(), "the estimate rose above the maximum"); + } + } + + @ParameterizedTest + @ValueSource(doubles = {0.01, 0.5, 0.99}) + @DisplayName("lands exactly on the textbook quantile of a ramp, which pins the marker convention") + void isExactOnARamp(double probability) { + // The markers are numbered from one, as in the paper. That origin is arithmetically irrelevant + // - only differences of positions are ever used - but the desired positions are grown by + // repeated addition, so a different origin rounds differently and flips the strict comparisons + // that pick between the parabolic and the linear prediction. On a ramp, whose quantiles are + // whole numbers, that shows up as an answer off by one. + P2QuantileEstimator estimator = new P2QuantileEstimator(probability); + for (int i = 0; i <= 1_700; i++) { + estimator.add(i); + } + assertEquals(probability * 1_700, estimator.quantile(), 1e-9); + } +}