malevich 0.7.0

Terminal plotting: a small grammar of marks, honest axes, millions of points
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

//! Kernel density estimation: a smooth distribution from a sample.

/// The Gaussian KDE of `values`, evaluated at `points` positions across the data
/// extent (padded by three bandwidths). Returns `(positions, densities)`, or `None`
/// without finite values.
///
/// Bandwidth follows Silverman's rule of thumb —
/// `0.9 * min(σ, IQR / 1.34) * n^(-1/5)` — and evaluation runs on linearly binned
/// counts with a truncated Gaussian kernel: O(n + points × kernel), no FFT.
pub fn kde(values: &[f64], points: usize) -> Option<(Vec<f64>, Vec<f64>)> {
    let finite: Vec<f64> = values.iter().copied().filter(|v| v.is_finite()).collect();
    if finite.is_empty() || points < 2 {
        return None;
    }
    let n = finite.len() as f64;
    let mean = finite.iter().sum::<f64>() / n;
    let variance = finite.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / n;
    let sigma = variance.sqrt();

    let mut sorted = finite.clone();
    sorted.sort_by(f64::total_cmp);
    let quantile = |p: f64| -> f64 {
        let position = (sorted.len() - 1) as f64 * p;
        let index = position.floor() as usize;
        let fraction = position - index as f64;
        if index + 1 < sorted.len() {
            sorted[index] + fraction * (sorted[index + 1] - sorted[index])
        } else {
            sorted[index]
        }
    };
    let iqr = quantile(0.75) - quantile(0.25);
    let spread = if iqr > 0.0 {
        sigma.min(iqr / 1.34)
    } else {
        sigma
    };
    let bandwidth = if spread > 0.0 {
        0.9 * spread * n.powf(-0.2)
    } else {
        // A degenerate sample still deserves a bump rather than a spike.
        1.0
    };

    let (low, high) = (sorted[0], sorted[sorted.len() - 1]);
    let start = low - 3.0 * bandwidth;
    let end = high + 3.0 * bandwidth;
    let step = (end - start) / (points - 1) as f64;

    // Linear binning onto the evaluation grid.
    let mut binned = vec![0.0f64; points];
    for &value in &finite {
        let position = (value - start) / step;
        let index = position.floor() as usize;
        let fraction = position - position.floor();
        if index + 1 < points {
            binned[index] += 1.0 - fraction;
            binned[index + 1] += fraction;
        } else {
            binned[points - 1] += 1.0;
        }
    }

    // Truncated Gaussian kernel over the binned counts.
    let radius = ((3.0 * bandwidth / step).ceil() as usize).max(1);
    let kernel: Vec<f64> = (0..=radius)
        .map(|k| {
            let distance = k as f64 * step / bandwidth;
            (-0.5 * distance * distance).exp()
        })
        .collect();
    let normalization = 1.0 / (n * bandwidth * (2.0 * std::f64::consts::PI).sqrt());

    let densities: Vec<f64> = (0..points)
        .map(|i| {
            let mut sum = binned[i] * kernel[0];
            for k in 1..=radius {
                if i >= k {
                    sum += binned[i - k] * kernel[k];
                }
                if i + k < points {
                    sum += binned[i + k] * kernel[k];
                }
            }
            sum * normalization
        })
        .collect();
    let positions = (0..points).map(|i| start + i as f64 * step).collect();
    Some((positions, densities))
}

#[cfg(test)]
#[path = "tests/kde_tests.rs"]
mod tests;