Skip to main content

plotters_statistical/stats/
histogram.rs

1//! Histogram binning with several standard bin-count rules.
2
3use super::{quartiles, sorted_finite, StatsError};
4
5/// How to choose the number/width of histogram bins.
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub enum BinRule {
8    /// A fixed number of equal-width bins.
9    Count(usize),
10    /// A fixed bin width (bin count derived from the data range).
11    Width(f64),
12    /// Sturges' rule: `ceil(log2(n)) + 1`.
13    Sturges,
14    /// Freedman–Diaconis rule: width `= 2 * IQR * n^(-1/3)`.
15    FreedmanDiaconis,
16    /// Scott's rule: width `= 3.49 * std * n^(-1/3)`.
17    Scott,
18}
19
20/// A computed histogram: bin edges (length `bins + 1`), per-bin counts, and
21/// per-bin density (`count / (n * width)`, so the areas sum to 1).
22#[derive(Debug, Clone, PartialEq)]
23pub struct Histogram {
24    /// Bin edges, ascending; `edges.len() == counts.len() + 1`.
25    pub edges: Vec<f64>,
26    /// Count of samples in each bin.
27    pub counts: Vec<usize>,
28    /// Density of each bin (`count / (n * bin_width)`).
29    pub density: Vec<f64>,
30}
31
32impl Histogram {
33    /// The center of each bin.
34    pub fn centers(&self) -> Vec<f64> {
35        self.edges.windows(2).map(|w| (w[0] + w[1]) / 2.0).collect()
36    }
37
38    /// The uniform bin width.
39    pub fn bin_width(&self) -> f64 {
40        if self.edges.len() < 2 {
41            0.0
42        } else {
43            self.edges[1] - self.edges[0]
44        }
45    }
46}
47
48/// Build a [`Histogram`] from `data` using the chosen [`BinRule`]. Non-finite
49/// values are ignored.
50///
51/// # Errors
52/// * [`StatsError::EmptyInput`] if no finite values remain.
53/// * [`StatsError::ZeroVariance`] if all values are identical (no range to bin).
54/// * [`StatsError::InvalidBandwidth`] if an explicit [`BinRule::Width`] or
55///   [`BinRule::Count`] is non-positive.
56pub fn histogram(data: &[f64], rule: BinRule) -> Result<Histogram, StatsError> {
57    let sorted = sorted_finite(data);
58    if sorted.is_empty() {
59        return Err(StatsError::EmptyInput);
60    }
61    let n = sorted.len();
62    let min = sorted[0];
63    let max = sorted[n - 1];
64    if max <= min {
65        return Err(StatsError::ZeroVariance);
66    }
67    let range = max - min;
68    let nf = n as f64;
69
70    let bins = match rule {
71        BinRule::Count(c) => {
72            if c == 0 {
73                return Err(StatsError::InvalidBandwidth);
74            }
75            c
76        }
77        BinRule::Width(w) => {
78            if !(w.is_finite() && w > 0.0) {
79                return Err(StatsError::InvalidBandwidth);
80            }
81            (range / w).ceil().max(1.0) as usize
82        }
83        BinRule::Sturges => (nf.log2().ceil() as usize) + 1,
84        BinRule::FreedmanDiaconis => {
85            let q = quartiles(&sorted)?;
86            let w = 2.0 * q.iqr * nf.powf(-1.0 / 3.0);
87            if w > 0.0 {
88                (range / w).ceil().max(1.0) as usize
89            } else {
90                (nf.log2().ceil() as usize) + 1 // fall back to Sturges on zero IQR
91            }
92        }
93        BinRule::Scott => {
94            let mean = sorted.iter().sum::<f64>() / nf;
95            let std = (sorted.iter().map(|x| (x - mean).powi(2)).sum::<f64>()
96                / (nf - 1.0).max(1.0))
97            .sqrt();
98            let w = 3.49 * std * nf.powf(-1.0 / 3.0);
99            if w > 0.0 {
100                (range / w).ceil().max(1.0) as usize
101            } else {
102                (nf.log2().ceil() as usize) + 1
103            }
104        }
105    }
106    .max(1);
107
108    let width = range / bins as f64;
109    let edges: Vec<f64> = (0..=bins).map(|i| min + width * i as f64).collect();
110    let mut counts = vec![0usize; bins];
111    for &v in &sorted {
112        let mut idx = ((v - min) / width).floor() as usize;
113        if idx >= bins {
114            idx = bins - 1; // the maximum lands in the last bin
115        }
116        counts[idx] += 1;
117    }
118    let density = counts.iter().map(|&c| c as f64 / (nf * width)).collect();
119
120    Ok(Histogram {
121        edges,
122        counts,
123        density,
124    })
125}