fix(tdigest): interpolate quantiles toward the nearer centroid and normalize the left tail

CHANGELOG.md

@@ -47,6 +47,7 @@ All significant changes to this project will be documented in this file.
 * HLL, Theta, and Tuple deserializers now return `InvalidData` for malformed payload sizes and entry counts instead of risking oversized allocations or decoding failures.
 * Malformed CPC images now return `InvalidData` instead of panicking.
 * Seeded deserializers now return `InvalidData` rather than panicking when the caller supplies a seed whose hash is the reserved zero value.
+* Fix T-Digest interpolation and tail calculations that could produce non-monotonic or out-of-range quantiles and invalid rank, CDF, or PMF values.
 
 ## v0.4.0 (2026-08-18)
 

datasketches/src/tdigest/sketch.rs

@@ -1278,8 +1278,9 @@ impl TDigestView<'_> {
                 return Some(if value == self.min {
                     0.5 / centroids_weight
                 } else {
-                    1. + (((value - self.min) / (first_mean - self.min))
-                        * ((self.centroids[0].weight() / 2.) - 1.))
+                    (1. + (((value - self.min) / (first_mean - self.min))
+                        * ((self.centroids[0].weight() / 2.) - 1.)))
+                        / centroids_weight
                 });
             }
             return Some(0.); // should never happen
@@ -1376,9 +1377,12 @@ impl TDigestView<'_> {
         }
         let last_weight = self.centroids[num_centroids - 1].weight();
         if last_weight > 1. && (centroids_weight - weight <= last_weight / 2.) {
+            if last_weight == 2. {
+                return Some(self.max);
+            }
             return Some(
                 self.max
-                    + (((centroids_weight - weight - 1.) / ((last_weight / 2.) - 1.))
+                    - (((centroids_weight - weight - 1.) / ((last_weight / 2.) - 1.))
                         * (self.max - self.centroids[num_centroids - 1].mean)),
             );
         }
@@ -1403,13 +1407,15 @@ impl TDigestView<'_> {
                     }
                     right_weight = 0.5;
                 }
-                let w1 = weight - weight_so_far - left_weight;
-                let w2 = weight_so_far + dw - weight - right_weight;
+                // Each centroid is weighted by the distance from the target to the *other*
+                // centroid, so the estimate approaches the nearer one.
+                let distance_from_left = weight - weight_so_far - left_weight;
+                let distance_to_right = weight_so_far + dw - weight - right_weight;
                 return Some(weighted_average(
                     self.centroids[i].mean,
-                    w1,
+                    distance_to_right,
                     self.centroids[i + 1].mean,
-                    w2,
+                    distance_from_left,
                 ));
             }
             weight_so_far += dw;

tests-integration/tests/tdigest_test/main.rs

@@ -15,4 +15,5 @@
 // specific language governing permissions and limitations
 // under the License.
 
+mod property;
 mod sketch;

tests-integration/tests/tdigest_test/property.rs

@@ -0,0 +1,89 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements.  See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership.  The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License.  You may obtain a copy of the License at
+//
+//   http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied.  See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! Property-based t-digest tests.
+
+use datasketches::tdigest::TDigestMut;
+use quickcheck::Gen;
+use quickcheck::QuickCheck;
+use quickcheck::TestResult;
+
+const RANK_STEPS: usize = 500;
+
+fn digest_of(values: &[u32]) -> TDigestMut {
+    let mut tdigest = TDigestMut::new(100).unwrap();
+    for value in values {
+        tdigest.update(f64::from(*value) / 4096.0);
+    }
+    tdigest
+}
+
+fn assert_quantiles_are_monotonic(tdigest: &mut TDigestMut) {
+    let min = tdigest.min_value().unwrap();
+    let max = tdigest.max_value().unwrap();
+    let mut previous = min;
+
+    for step in 0..=RANK_STEPS {
+        let rank = step as f64 / RANK_STEPS as f64;
+        let quantile = tdigest.quantile(rank).unwrap();
+        assert!(
+            (previous..=max).contains(&quantile),
+            "quantile {quantile} at rank {rank} is outside [{previous}, {max}]"
+        );
+        previous = quantile;
+    }
+}
+
+#[test]
+fn prop_quantile_is_non_decreasing_and_within_the_observed_range() {
+    fn property(values: Vec<u32>) -> TestResult {
+        if !(500..1500).contains(&values.len()) {
+            return TestResult::discard();
+        }
+
+        assert_quantiles_are_monotonic(&mut digest_of(&values));
+
+        TestResult::passed()
+    }
+
+    QuickCheck::new()
+        .tests(128)
+        .min_tests_passed(128)
+        .rng(Gen::new(1200))
+        .quickcheck(property as fn(Vec<u32>) -> TestResult);
+}
+
+#[test]
+fn prop_merged_quantile_is_non_decreasing() {
+    fn property(left: Vec<u32>, right: Vec<u32>) -> TestResult {
+        if left.len() < 300 || right.len() < 300 {
+            return TestResult::discard();
+        }
+
+        let mut tdigest = digest_of(&left);
+        tdigest.merge(&digest_of(&right));
+        assert_quantiles_are_monotonic(&mut tdigest);
+
+        TestResult::passed()
+    }
+
+    QuickCheck::new()
+        .tests(64)
+        .min_tests_passed(64)
+        .rng(Gen::new(900))
+        .quickcheck(property as fn(Vec<u32>, Vec<u32>) -> TestResult);
+}

tests-integration/tests/tdigest_test/sketch.rs

@@ -332,3 +332,80 @@ fn test_estimate_repeat_values() {
     }
     assert_eq!(tdigest.quantile(0.9), Some(1.0));
 }
+
+/// Builds a digest whose centroids carry the given weights.
+///
+/// Compression never merges the extreme centroids, so digests built through `update` and `merge`
+/// always keep unit-weight tails. Heavier tails arrive only through deserialization, including the
+/// reference implementation format, and they select the tail interpolation branches.
+fn deserialize_with_centroids(k: u16, min: f64, max: f64, centroids: &[(f64, u64)]) -> TDigestMut {
+    const PREAMBLE_LONGS: u8 = 2;
+    const SERIAL_VERSION: u8 = 1;
+    const FAMILY_TDIGEST: u8 = 20;
+
+    let mut bytes = vec![PREAMBLE_LONGS, SERIAL_VERSION, FAMILY_TDIGEST];
+    bytes.extend_from_slice(&k.to_le_bytes());
+    bytes.push(0); // flags
+    bytes.extend_from_slice(&0u16.to_le_bytes()); // unused
+    bytes.extend_from_slice(&(centroids.len() as u32).to_le_bytes());
+    bytes.extend_from_slice(&0u32.to_le_bytes()); // buffered values
+    bytes.extend_from_slice(&min.to_le_bytes());
+    bytes.extend_from_slice(&max.to_le_bytes());
+    for (mean, weight) in centroids {
+        bytes.extend_from_slice(&mean.to_le_bytes());
+        bytes.extend_from_slice(&weight.to_le_bytes());
+    }
+    TDigestMut::deserialize(&bytes).unwrap()
+}
+
+#[test]
+fn test_quantile_moves_toward_the_nearer_bracketing_centroid() {
+    let mut tdigest =
+        deserialize_with_centroids(100, -1.0, 21.0, &[(0.0, 4), (10.0, 4), (20.0, 4)]);
+
+    assert_eq!(tdigest.total_weight(), 12);
+    // Ranks 2/12 and 6/12 sit exactly on the two centroids bracketing the first interval.
+    assert_that!(tdigest.quantile(2.0 / 12.0).unwrap(), near(0.0, 1e-12));
+    assert_that!(tdigest.quantile(3.0 / 12.0).unwrap(), near(2.5, 1e-12));
+    assert_that!(tdigest.quantile(4.0 / 12.0).unwrap(), near(5.0, 1e-12));
+    assert_that!(tdigest.quantile(5.0 / 12.0).unwrap(), near(7.5, 1e-12));
+    assert_that!(tdigest.quantile(6.0 / 12.0).unwrap(), near(10.0, 1e-12));
+}
+
+#[test]
+fn test_quantile_right_tail_stays_within_max() {
+    let mut tdigest =
+        deserialize_with_centroids(100, 0.0, 100.0, &[(10.0, 10), (50.0, 10), (90.0, 10)]);
+
+    assert_eq!(tdigest.max_value(), Some(100.0));
+    assert_that!(tdigest.quantile(0.9).unwrap(), near(95.0, 1e-12));
+    assert_that!(tdigest.quantile(29.0 / 30.0).unwrap(), near(100.0, 1e-12));
+    // Mirrors the left tail, which interpolates from min up to the first centroid mean.
+    assert_that!(tdigest.quantile(1.0 / 30.0).unwrap(), near(0.0, 1e-12));
+    assert_that!(tdigest.quantile(5.0 / 30.0).unwrap(), near(10.0, 1e-12));
+}
+
+#[test]
+fn test_quantile_handles_two_sample_last_centroid() {
+    let mut tdigest =
+        deserialize_with_centroids(100, 0.0, 100.0, &[(0.0, 1), (50.0, 1), (90.0, 2)]);
+
+    assert_eq!(tdigest.quantile(0.75), Some(100.0));
+}
+
+#[test]
+fn test_rank_left_tail_is_a_fraction_of_the_total_weight() {
+    let mut tdigest =
+        deserialize_with_centroids(100, 0.0, 100.0, &[(10.0, 10), (50.0, 10), (90.0, 10)]);
+
+    assert_that!(tdigest.rank(5.0).unwrap(), near(0.1, 1e-12));
+    assert_that!(tdigest.rank(10.0).unwrap(), near(5.0 / 30.0, 1e-12));
+    // The right tail is the mirror image and pins the scale the left tail must match.
+    assert_that!(tdigest.rank(95.0).unwrap(), near(0.9, 1e-12));
+    assert_that!(tdigest.rank(90.0).unwrap(), near(25.0 / 30.0, 1e-12));
+
+    let pmf = tdigest.pmf(&[5.0, 95.0]).unwrap();
+    assert_that!(pmf[0], near(0.1, 1e-12));
+    assert_that!(pmf[1], near(0.8, 1e-12));
+    assert_that!(pmf[2], near(0.1, 1e-12));
+}