/// Percentile algorithm.
-#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+#[derive(Copy, Clone, Debug, Default, PartialEq, Eq)]
pub enum Algorithm {
/// Weighted average at lower observation.
WeightedAverage,
/// 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,
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 => {
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?}"
+ );
+ }
}
}
}
projects / pspp / commitdiff