From: Ben Pfaff <blp@cs.stanford.edu>
Date: Tue, 20 Oct 2026 01:55:44 +0000 (-0700)
Subject: work on percentiles
X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=13ca2ef180ef104bdc5236a03ff9d0f2b63e7ecc;p=pspp

work on percentiles
---

diff --git a/rust/pspp/src/percentiles.rs b/rust/pspp/src/percentiles.rs
index 64fa44270c..7a8d3ef63a 100644
--- a/rust/pspp/src/percentiles.rs
+++ b/rust/pspp/src/percentiles.rs
@@ -1,5 +1,5 @@
 /// Percentile algorithm.
-#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+#[derive(Copy, Clone, Debug, Default, PartialEq, Eq)]
 pub enum Algorithm {
     /// Weighted average at lower observation.
     WeightedAverage,
@@ -8,6 +8,10 @@ pub enum Algorithm {
     /// Empirical distribution function.
     Empirical,
     /// Weighted average at higher observation.
+    ///
+    /// The `FREQUENCIES` command uses this algorithm, and this is the default
+    /// algorithm for `EXAMINE`.
+    #[default]
     HighAverage,
     /// Empirical distribution function with averaging.
     AverageEmpirical,
@@ -23,24 +27,41 @@ impl Algorithm {
             Algorithm::HighAverage => 1.0,
         }
     }
-    pub fn percentile(self, data: &impl Observations, p: f64) -> Option<f64> {
-        let w = data.total_weight();
-
-        // Percentile must be between 0 and 1.
-        let p = p.clamp(0.0, 1.0);
-
-        // Target cumulative weight.
-        let tc = (w + self.weight_adjustment()) * p;
+    pub fn target_cumulative_weight(self, total_weight: f64, p: f64) -> f64 {
+        (total_weight + self.weight_adjustment()) * p.clamp(0.0, 1.0)
+    }
 
-        let (a, b) = match data.observations_around(tc) {
+    /// Calculates a percentile for this algorithm at  target cumulative weight `tcw`,
+    /// given that `(a, b)` brackets the observation with that cumulative weight:
+    ///
+    /// - `a` is the observation with the least cumulative weight whose
+    ///   cumulative weight is at least `tcw`.
+    ///
+    /// - `b` is the observation with cumulative weight just above `a`.
+    ///
+    /// Corner cases:
+    ///
+    /// - `a` and `b` are both `None` if there are no observations.
+    ///
+    /// - If `a` is the observation with the greatest cumulative weight, then
+    ///   `b` is `None`.
+    ///
+    /// - If `tcw` is less than the least cumulative weight among the
+    ///   observations, then `b` is the observation with the least cumulative
+    ///   weight and `a` is `None`.
+    pub fn percentile(
+        self,
+        tcw: f64,
+        (a, b): (Option<Observation>, Option<Observation>),
+    ) -> Option<f64> {
+        let (a, b) = match (a, b) {
             (None, None) => return None,
             (None, Some(x)) | (Some(x), None) => return Some(x.y),
             (Some(a), Some(b)) => (a, b),
         };
-        //dbg!([a, b]);
 
-        let g = (tc - a.cc) / b.c;
-        let gs = tc - a.cc;
+        let g = (tcw - a.cc) / b.c;
+        let gs = tcw - a.cc;
 
         Some(match self {
             Algorithm::WeightedAverage | Algorithm::HighAverage => {
@@ -88,114 +109,212 @@ pub struct Observation {
     cc: f64,
 }
 
-/// Observations suitable for calculating percentiles.
-///
-/// Conceptually, this trait represents an array of [Observation] with strictly
-/// increasing values of `y`.
-pub trait Observations {
-    /// Returns the sum of weights of all the observations.
-    fn total_weight(&self) -> f64;
-
-    /// Finds and returns the first observation whose cumulative weight is at
-    /// least `cum_weight`, and the observation after that (if it is not the
-    /// last observation).
-    ///
-    /// If `cum_weight` is greater than the total weight, returns an empty
-    /// vector.
-    fn observations_around(&self, cum_weight: f64) -> (Option<Observation>, Option<Observation>);
+pub trait Percentile {
+    /// Calculates percentile `0 <= p <= 1` using `algorithm` on the data
+    /// represented by this object.
+    fn percentile(&self, algorithm: Algorithm, p: f64) -> Option<f64>;
 }
 
-struct Observation2 {
+trait CumulativeFrequency {
     /// Observed value.
-    y: f64,
+    fn y(&self) -> f64;
 
     /// Sum of weights up to and including this observation.
-    cc: f64,
+    fn cc(&self) -> f64;
 }
 
-impl Observations for [Observation2] {
-    fn total_weight(&self) -> f64 {
-        self.last().map_or(0.0, |last| last.cc)
-    }
-    fn observations_around(&self, cum_weight: f64) -> (Option<Observation>, Option<Observation>) {
-        fn get(this: &[Observation2], k: usize) -> Observation {
-            let Observation2 { y, cc } = this[k];
-            let c = if k == 0 { cc } else { cc - this[k - 1].cc };
-            Observation { y, c, cc }
-        }
-
-        let Some(k) = self
-            .partition_point(|obs| obs.cc <= cum_weight)
-            .checked_sub(1)
-        else {
-            if self.is_empty() {
-                return (None, None);
-            } else {
-                return (None, Some(get(self, 0)));
+impl<T> Percentile for [T]
+where
+    T: CumulativeFrequency,
+{
+    /// Calculates percentile `0 <= p <= 1` using `algorithm` on the data
+    /// represented by this object.
+    ///
+    /// The frequencies in the slice must be in strictly increasing order of
+    /// y-value and of cumulative weight.
+    fn percentile(&self, algorithm: Algorithm, p: f64) -> Option<f64> {
+        fn observations_around<T>(
+            this: &[T],
+            cum_weight: f64,
+        ) -> (Option<Observation>, Option<Observation>)
+        where
+            T: CumulativeFrequency,
+        {
+            fn get<T>(this: &[T], k: usize) -> Observation
+            where
+                T: CumulativeFrequency,
+            {
+                let y = this[k].y();
+                let cc = this[k].cc();
+                let c = if k == 0 { cc } else { cc - this[k - 1].cc() };
+                Observation { y, c, cc }
             }
-        };
 
-        if k + 1 >= self.len() {
-            (Some(get(self, k)), None)
-        } else {
-            (Some(get(self, k)), Some(get(self, k + 1)))
+            match this
+                .partition_point(|obs| obs.cc() <= cum_weight)
+                .checked_sub(1)
+            {
+                Some(k) => {
+                    if k + 1 >= this.len() {
+                        (Some(get(this, k)), None)
+                    } else {
+                        (Some(get(this, k)), Some(get(this, k + 1)))
+                    }
+                }
+                None => {
+                    if this.is_empty() {
+                        (None, None)
+                    } else {
+                        (None, Some(get(this, 0)))
+                    }
+                }
+            }
         }
-    }
-}
 
-impl Observations for Vec<Observation2> {
-    fn total_weight(&self) -> f64 {
-        self.as_slice().total_weight()
+        let total_weight = self.last().map_or(0.0, |last| last.cc());
+        let tcw = algorithm.target_cumulative_weight(total_weight, p);
+        algorithm.percentile(tcw, observations_around(self, tcw))
     }
+}
 
-    fn observations_around(&self, cum_weight: f64) -> (Option<Observation>, Option<Observation>) {
-        self.as_slice().observations_around(cum_weight)
+impl<T> Percentile for Vec<T>
+where
+    T: CumulativeFrequency,
+{
+    /// Calculates percentile `0 <= p <= 1` using `algorithm` on the data
+    /// represented by this object.
+    ///
+    /// The frequencies in the slice must be in strictly increasing order of
+    /// y-value and of cumulative weight.
+    fn percentile(&self, algorithm: Algorithm, p: f64) -> Option<f64> {
+        self.as_slice().percentile(algorithm, p)
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use crate::percentiles::{Algorithm::HighAverage, Observation2};
+    use crate::percentiles::{Algorithm, CumulativeFrequency, Percentile as _};
+
+    struct Item {
+        y: f64,
+        cc: f64,
+    }
+    impl CumulativeFrequency for Item {
+        fn y(&self) -> f64 {
+            self.y
+        }
+
+        fn cc(&self) -> f64 {
+            self.cc
+        }
+    }
 
     #[test]
     fn zero_observations() {
-        let data = vec![];
-        for p in (0..=20).map(|p| (p * 5) as f64 / 100.0) {
-            println!("{p:.2}: {:.2?}", HighAverage.percentile(&data, p));
+        let data = Vec::<Item>::new();
+        for p in (0..=20).map(|p| p as f64 / 20.0) {
+            assert_eq!(data.percentile(Algorithm::HighAverage, p), None);
         }
-        //assert_eq!(HighAverage.percentile(&data, 0.05), Some(5.0));
     }
 
     #[test]
     fn one_observation() {
-        let data = vec![Observation2 { y: 5.0, cc: 1.0 }];
-        for p in (0..=20).map(|p| (p * 5) as f64 / 100.0) {
-            println!("{p:.2}: {:.2?}", HighAverage.percentile(&data, p));
+        let data = vec![Item { y: 5.0, cc: 1.0 }];
+        for p in (0..=20).map(|p| p as f64 / 20.0) {
+            assert_eq!(data.percentile(Algorithm::HighAverage, p), Some(5.0));
         }
-        //assert_eq!(HighAverage.percentile(&data, 0.05), Some(5.0));
     }
 
     #[test]
     fn two_observations() {
-        let data = vec![
-            Observation2 { y: 2.0, cc: 1.0 },
-            Observation2 { y: 5.0, cc: 2.0 },
-        ];
-        for p in (0..=20).map(|p| (p * 5) as f64 / 100.0) {
-            println!("{p:.2}: {:.2?}", HighAverage.percentile(&data, p));
+        let data = vec![Item { y: 2.0, cc: 1.0 }, Item { y: 5.0, cc: 2.0 }];
+
+        for (p, expected) in [
+            (0.00, "2.00"),
+            (0.05, "2.00"),
+            (0.10, "2.00"),
+            (0.15, "2.00"),
+            (0.20, "2.00"),
+            (0.25, "2.00"),
+            (0.30, "2.00"),
+            (0.35, "2.15"),
+            (0.40, "2.60"),
+            (0.45, "3.05"),
+            (0.50, "3.50"),
+            (0.55, "3.95"),
+            (0.60, "4.40"),
+            (0.65, "4.85"),
+            (0.70, "5.00"),
+            (0.75, "5.00"),
+            (0.80, "5.00"),
+            (0.85, "5.00"),
+            (0.90, "5.00"),
+            (0.95, "5.00"),
+            (1.00, "5.00"),
+        ] {
+            let percentile = data.percentile(Algorithm::HighAverage, p).unwrap();
+            assert_eq!(format!("{:.2}", percentile), expected);
         }
-        //assert_eq!(HighAverage.percentile(&data, 0.05), Some(5.0));
     }
 
     #[test]
     fn three_observations() {
         let data = vec![
-            Observation2 { y: 2.0, cc: 1.0 },
-            Observation2 { y: 5.0, cc: 2.0 },
-            Observation2 { y: 8.0, cc: 3.0 },
+            Item { y: 2.0, cc: 1.0 },
+            Item { y: 5.0, cc: 2.0 },
+            Item { y: 8.0, cc: 3.0 },
         ];
-        for p in (0..=20).map(|p| (p * 5) as f64 / 100.0) {
-            println!("{p:.2}: {:.2?}", HighAverage.percentile(&data, p));
+
+        let algorithms = [
+            Algorithm::HighAverage,
+            Algorithm::WeightedAverage,
+            Algorithm::Round,
+            Algorithm::Empirical,
+            Algorithm::AverageEmpirical,
+        ];
+        /*
+               for p in (0..=20).map(|p| p as f64 / 20.0) {
+                   println!(
+                       "({p:.2}, [{}]),",
+                       algorithms
+                           .into_iter()
+                           .map(|algorithm| format!("\"{:.2}\"", data.percentile(algorithm, p).unwrap()))
+                           .join(", ")
+                   );
+               }
+        */
+
+        for (p, expected) in [
+            (0.00, ["2.00", "2.00", "2.00", "2.00", "2.00"]),
+            (0.05, ["2.00", "2.00", "2.00", "2.00", "2.00"]),
+            (0.10, ["2.00", "2.00", "2.00", "2.00", "2.00"]),
+            (0.15, ["2.00", "2.00", "2.00", "2.00", "2.00"]),
+            (0.20, ["2.00", "2.00", "2.00", "2.00", "2.00"]),
+            (0.25, ["2.00", "2.00", "2.00", "2.00", "2.00"]),
+            (0.30, ["2.60", "2.00", "2.00", "2.00", "2.00"]),
+            (0.35, ["3.20", "2.15", "2.00", "5.00", "5.00"]),
+            (0.40, ["3.80", "2.60", "2.00", "5.00", "5.00"]),
+            (0.45, ["4.40", "3.05", "2.00", "5.00", "5.00"]),
+            (0.50, ["5.00", "3.50", "5.00", "5.00", "5.00"]),
+            (0.55, ["5.60", "3.95", "5.00", "5.00", "5.00"]),
+            (0.60, ["6.20", "4.40", "5.00", "5.00", "5.00"]),
+            (0.65, ["6.80", "4.85", "5.00", "5.00", "5.00"]),
+            (0.70, ["7.40", "5.30", "5.00", "8.00", "8.00"]),
+            (0.75, ["8.00", "5.75", "5.00", "8.00", "8.00"]),
+            (0.80, ["8.00", "6.20", "5.00", "8.00", "8.00"]),
+            (0.85, ["8.00", "6.65", "8.00", "8.00", "8.00"]),
+            (0.90, ["8.00", "7.10", "8.00", "8.00", "8.00"]),
+            (0.95, ["8.00", "7.55", "8.00", "8.00", "8.00"]),
+            (1.00, ["8.00", "8.00", "8.00", "8.00", "8.00"]),
+        ] {
+            for (algorithm, expected) in algorithms.into_iter().zip(expected) {
+                let percentile = data.percentile(algorithm, p).unwrap();
+                assert_eq!(
+                    format!("{:.2}", percentile),
+                    expected,
+                    "{algorithm:?} percentile {p:.2?}"
+                );
+            }
         }
     }
 }