From fa2e8e32a06234dce20fc7904daf6ef61241ab82 Mon Sep 17 00:00:00 2001 From: Ben Pfaff Date: Fri, 16 Oct 2026 09:46:56 -0700 Subject: [PATCH] Work on percentile calculation --- rust/pspp/src/lib.rs | 1 + rust/pspp/src/percentiles.rs | 201 +++++++++++++++++++++++++++++++++++ 2 files changed, 202 insertions(+) create mode 100644 rust/pspp/src/percentiles.rs diff --git a/rust/pspp/src/lib.rs b/rust/pspp/src/lib.rs index acc1e8d4d0..4e1d83cdd3 100644 --- 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 index 0000000000..64fa44270c --- /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 { + 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, Option); +} + +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, Option) { + 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 { + fn total_weight(&self) -> f64 { + self.as_slice().total_weight() + } + + fn observations_around(&self, cum_weight: f64) -> (Option, Option) { + 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)); + } + } +} -- 2.47.3