plotters_statistical/stats/
kde.rs

1//! Gaussian kernel density estimation, used to shape violin plots.
2
3use super::{sorted_finite, StatsError};
4use std::f64::consts::PI;
5
6/// An evaluated kernel density estimate: paired `x` positions and density
7/// values, plus the bandwidth that produced them.
8#[derive(Debug, Clone, PartialEq)]
9pub struct KdeCurve {
10    /// Evaluation positions along the value axis, ascending.
11    pub xs: Vec<f64>,
12    /// Estimated density at each corresponding `xs` position. Integrates to ~1.
13    pub density: Vec<f64>,
14    /// The bandwidth (kernel standard deviation) actually used.
15    pub bandwidth: f64,
16}
17
18impl KdeCurve {
19    /// The largest density value on the curve — handy for scaling a violin's
20    /// half-width so its widest point maps to a fixed pixel/coordinate extent.
21    pub fn max_density(&self) -> f64 {
22        self.density.iter().copied().fold(0.0_f64, f64::max)
23    }
24}
25
26/// Silverman's rule-of-thumb bandwidth:
27/// `h = 0.9 * min(std, IQR/1.34) * n^(-1/5)`.
28///
29/// Uses the sample standard deviation (n − 1 denominator). Returns
30/// [`StatsError::InvalidBandwidth`] when the estimate is not finite and
31/// positive — e.g. fewer than two points, or an all-identical sample where both
32/// the spread and the IQR are zero.
33pub fn silverman_bandwidth(data: &[f64]) -> Result<f64, StatsError> {
34    let sorted = sorted_finite(data);
35    let n = sorted.len();
36    if n < 2 {
37        return Err(StatsError::InvalidBandwidth);
38    }
39    let nf = n as f64;
40    let mean = sorted.iter().sum::<f64>() / nf;
41    let var = sorted.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / (nf - 1.0);
42    let std = var.sqrt();
43
44    // IQR via the same type-7 percentile used for box plots.
45    let q1 = super::percentile_sorted(&sorted, 0.25);
46    let q3 = super::percentile_sorted(&sorted, 0.75);
47    let iqr = q3 - q1;
48
49    let spread = if iqr > 0.0 { std.min(iqr / 1.34) } else { std };
50    let h = 0.9 * spread * nf.powf(-1.0 / 5.0);
51    if h.is_finite() && h > 0.0 {
52        Ok(h)
53    } else {
54        Err(StatsError::InvalidBandwidth)
55    }
56}
57
58fn gaussian_kernel(u: f64) -> f64 {
59    (-0.5 * u * u).exp() / (2.0 * PI).sqrt()
60}
61
62/// Evaluate a Gaussian KDE of `data` on a uniform grid.
63///
64/// * `bandwidth` — `Some(h)` to force a bandwidth, or `None` to derive one via
65///   [`silverman_bandwidth`].
66/// * `n_points` — number of grid points; higher gives a smoother violin outline
67///   (values around 200 render cleanly). Clamped to at least 2.
68/// * `cut` — how far past the data range to extend, in bandwidths (a value of
69///   `3.0` captures the Gaussian tails without a visible cutoff).
70///
71/// Returns [`StatsError::EmptyInput`] for no finite data, or
72/// [`StatsError::InvalidBandwidth`] for a non-positive bandwidth (explicit or
73/// derived).
74pub fn kde_curve(
75    data: &[f64],
76    bandwidth: Option<f64>,
77    n_points: usize,
78    cut: f64,
79) -> Result<KdeCurve, StatsError> {
80    let sorted = sorted_finite(data);
81    if sorted.is_empty() {
82        return Err(StatsError::EmptyInput);
83    }
84    let h = match bandwidth {
85        Some(b) if b.is_finite() && b > 0.0 => b,
86        Some(_) => return Err(StatsError::InvalidBandwidth),
87        None => silverman_bandwidth(&sorted)?,
88    };
89
90    let n_points = n_points.max(2);
91    let lo = sorted[0] - cut * h;
92    let hi = sorted[sorted.len() - 1] + cut * h;
93    let step = (hi - lo) / (n_points as f64 - 1.0);
94
95    let nf = sorted.len() as f64;
96    let mut xs = Vec::with_capacity(n_points);
97    let mut density = Vec::with_capacity(n_points);
98    for i in 0..n_points {
99        let x = lo + step * i as f64;
100        let d = sorted
101            .iter()
102            .map(|&xi| gaussian_kernel((x - xi) / h))
103            .sum::<f64>()
104            / (nf * h);
105        xs.push(x);
106        density.push(d);
107    }
108
109    Ok(KdeCurve {
110        xs,
111        density,
112        bandwidth: h,
113    })
114}
plotters_statistical/stats/kde.rs

plotters_statistical/stats/
kde.rs
