Work on percentile calculation

diff --git a/rust/pspp/src/lib.rs b/rust/pspp/src/lib.rs
index acc1e8d4d03041978bb2967aa3c91b65df7cf417..4e1d83cdd307925d8de787700edfbd8db6d96a9d 100644 (file)
--- a/rust/pspp/src/lib.rs
+++ b/rust/pspp/src/lib.rs
@@ -118,6 +118,7 @@ pub mod macros;
 pub mod message;
 pub mod output;
 pub mod pc;
+pub mod percentiles;
 pub mod por;
 pub mod prompt;
 pub mod settings;

diff --git a/rust/pspp/src/percentiles.rs b/rust/pspp/src/percentiles.rs
new file mode 100644 (file)
index 0000000..64fa442
--- /dev/null
+++ b/rust/pspp/src/percentiles.rs
@@ -0,0 +1,201 @@
+/// Percentile algorithm.
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+pub enum Algorithm {
+    /// Weighted average at lower observation.
+    WeightedAverage,
+    /// Closest observation.
+    Round,
+    /// Empirical distribution function.
+    Empirical,
+    /// Weighted average at higher observation.
+    HighAverage,
+    /// Empirical distribution function with averaging.
+    AverageEmpirical,
+}
+
+impl Algorithm {
+    fn weight_adjustment(self) -> f64 {
+        match self {
+            Algorithm::WeightedAverage
+            | Algorithm::Round
+            | Algorithm::Empirical
+            | Algorithm::AverageEmpirical => 0.0,
+            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;
+
+        let (a, b) = match data.observations_around(tc) {
+            (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;
+
+        Some(match self {
+            Algorithm::WeightedAverage | Algorithm::HighAverage => {
+                if gs >= 1.0 {
+                    b.y
+                } else if b.c >= 1.0 {
+                    (1.0 - gs) * a.y + gs * b.y
+                } else {
+                    (1.0 - g) * a.y + g * b.y
+                }
+            }
+            Algorithm::Round => {
+                let threshold = if b.c >= 1.0 { gs } else { g };
+                if threshold < 0.5 { a.y } else { b.y }
+            }
+            Algorithm::Empirical => {
+                if gs == 0.0 {
+                    a.y
+                } else {
+                    b.y
+                }
+            }
+            Algorithm::AverageEmpirical => {
+                if gs == 0.0 {
+                    a.y.midpoint(b.y)
+                } else {
+                    b.y
+                }
+            }
+        })
+    }
+}
+
+#[derive(Copy, Clone, Debug)]
+pub struct Observation {
+    /// Observed value.
+    y: f64,
+
+    /// Weight.
+    ///
+    /// Must be positive.
+    c: f64,
+
+    /// Sum of weights up to and including this 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>);
+}
+
+struct Observation2 {
+    /// Observed value.
+    y: f64,
+
+    /// Sum of weights up to and including this observation.
+    cc: 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)));
+            }
+        };
+
+        if k + 1 >= self.len() {
+            (Some(get(self, k)), None)
+        } else {
+            (Some(get(self, k)), Some(get(self, k + 1)))
+        }
+    }
+}
+
+impl Observations for Vec<Observation2> {
+    fn total_weight(&self) -> f64 {
+        self.as_slice().total_weight()
+    }
+
+    fn observations_around(&self, cum_weight: f64) -> (Option<Observation>, Option<Observation>) {
+        self.as_slice().observations_around(cum_weight)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::percentiles::{Algorithm::HighAverage, Observation2};
+
+    #[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));
+        }
+        //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));
+        }
+        //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));
+        }
+        //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 },
+        ];
+        for p in (0..=20).map(|p| (p * 5) as f64 / 100.0) {
+            println!("{p:.2}: {:.2?}", HighAverage.percentile(&data, p));
+        }
+    }
+}