Skip to main content

ggplot_rs/stat/
ecdf.rs

1use crate::aes::Aesthetic;
2use crate::data::{DataFrame, Value};
3use crate::scale::ScaleSet;
4
5use super::Stat;
6
7/// Empirical cumulative distribution function.
8/// Sorts x values and assigns y = rank / n.
9pub struct StatEcdf;
10
11impl Default for StatEcdf {
12    fn default() -> Self {
13        StatEcdf
14    }
15}
16
17impl Stat for StatEcdf {
18    fn compute_group(&self, data: &DataFrame, _scales: &ScaleSet) -> DataFrame {
19        let x_col = match data.column("x") {
20            Some(c) => c,
21            None => return DataFrame::new(),
22        };
23
24        let mut values: Vec<f64> = x_col
25            .iter()
26            .filter_map(|v| v.as_f64())
27            .filter(|v| v.is_finite())
28            .collect();
29        if values.is_empty() {
30            return DataFrame::new();
31        }
32
33        values.sort_by(|a, b| a.total_cmp(b));
34        let n = values.len() as f64;
35
36        let mut x_vals = Vec::with_capacity(values.len() + 2);
37        let mut y_vals = Vec::with_capacity(values.len() + 2);
38
39        // ggplot2 pads the step to ±Inf (y = 0 before the first point, y = 1
40        // after the last) so it spans the panel. Scales ignore non-finite values
41        // when training, and geom_step clamps the ±Inf segments to the panel edge.
42        x_vals.push(Value::Float(f64::NEG_INFINITY));
43        y_vals.push(Value::Float(0.0));
44        for (i, &x) in values.iter().enumerate() {
45            x_vals.push(Value::Float(x));
46            y_vals.push(Value::Float((i + 1) as f64 / n));
47        }
48        x_vals.push(Value::Float(f64::INFINITY));
49        y_vals.push(Value::Float(1.0));
50
51        let mut result = DataFrame::new();
52        result.add_column("x".to_string(), x_vals);
53        result.add_column("y".to_string(), y_vals);
54
55        // Carry over grouping columns
56        let nrows = values.len() + 2;
57        for col_name in &["color", "fill", "group"] {
58            if let Some(col) = data.column(col_name) {
59                if let Some(first) = col.first() {
60                    result.add_column(col_name.to_string(), vec![first.clone(); nrows]);
61                }
62            }
63        }
64
65        result
66    }
67
68    fn required_aes(&self) -> Vec<Aesthetic> {
69        vec![Aesthetic::X]
70    }
71
72    fn computed_y(&self) -> bool {
73        true
74    }
75
76    fn name(&self) -> &str {
77        "ecdf"
78    }
79}
80
81/// ECDF with a simultaneous Dvoretzky–Kiefer–Wolfowitz confidence band:
82/// `F̂(x) ± ε`, `ε = √(ln(2 / (1 − level)) / (2n))`, clamped to `[0, 1]`.
83/// Output `x` (±Inf padded, as [`StatEcdf`]), `y`, `ymin`, `ymax` — draw it
84/// with `geom_stepribbon` (see `GGPlot::stat_ecdf_band`).
85#[derive(Clone, Debug)]
86pub struct StatEcdfBand {
87    /// Confidence level (default 0.95).
88    pub level: f64,
89}
90
91impl Default for StatEcdfBand {
92    fn default() -> Self {
93        StatEcdfBand { level: 0.95 }
94    }
95}
96
97impl StatEcdfBand {
98    pub fn new(level: f64) -> Self {
99        StatEcdfBand { level }
100    }
101
102    /// The DKW half-width `ε` for `n` observations at this level.
103    pub fn epsilon(&self, n: usize) -> f64 {
104        ((2.0 / (1.0 - self.level)).ln() / (2.0 * n as f64)).sqrt()
105    }
106}
107
108impl Stat for StatEcdfBand {
109    fn compute_group(&self, data: &DataFrame, scales: &ScaleSet) -> DataFrame {
110        if !(self.level > 0.0 && self.level < 1.0) {
111            return DataFrame::new();
112        }
113        // Only finite x count towards n (±Inf / NaN input rows are dropped).
114        let finite: Vec<Value> = data
115            .column("x")
116            .map(|c| {
117                c.iter()
118                    .filter(|v| v.as_f64().is_some_and(f64::is_finite))
119                    .cloned()
120                    .collect()
121            })
122            .unwrap_or_default();
123        let n = finite.len();
124        if n == 0 {
125            return DataFrame::new();
126        }
127        let mut input = DataFrame::new();
128        input.add_column("x".into(), finite);
129        for col in ["color", "fill", "group"] {
130            if let Some(first) = data.column(col).and_then(|c| c.first()) {
131                input.add_column(col.into(), vec![first.clone(); n]);
132            }
133        }
134        let mut out = StatEcdf.compute_group(&input, scales);
135        let eps = self.epsilon(n);
136        let y: Vec<f64> = out
137            .column("y")
138            .map(|c| c.iter().filter_map(|v| v.as_f64()).collect())
139            .unwrap_or_default();
140        out.add_column(
141            "ymin".into(),
142            y.iter()
143                .map(|&v| Value::Float((v - eps).max(0.0)))
144                .collect(),
145        );
146        out.add_column(
147            "ymax".into(),
148            y.iter()
149                .map(|&v| Value::Float((v + eps).min(1.0)))
150                .collect(),
151        );
152        out
153    }
154
155    fn required_aes(&self) -> Vec<Aesthetic> {
156        vec![Aesthetic::X]
157    }
158
159    fn name(&self) -> &str {
160        "ecdf_band"
161    }
162}