Skip to main content

solvr/stats/cpu/
descriptive.rs

1//! CPU implementation of descriptive statistics algorithms.
2
3use crate::stats::TensorDescriptiveStats;
4use crate::stats::impl_generic::{
5    describe_impl, iqr_impl, kurtosis_impl, percentile_impl, sem_impl, skewness_impl, zscore_impl,
6};
7use crate::stats::traits::DescriptiveStatisticsAlgorithms;
8use numr::error::Result;
9use numr::runtime::cpu::{CpuClient, CpuRuntime};
10use numr::tensor::Tensor;
11
12impl DescriptiveStatisticsAlgorithms<CpuRuntime> for CpuClient {
13    fn describe(&self, x: &Tensor<CpuRuntime>) -> Result<TensorDescriptiveStats<CpuRuntime>> {
14        describe_impl(self, x)
15    }
16
17    fn percentile(&self, x: &Tensor<CpuRuntime>, p: f64) -> Result<Tensor<CpuRuntime>> {
18        percentile_impl(self, x, p)
19    }
20
21    fn iqr(&self, x: &Tensor<CpuRuntime>) -> Result<Tensor<CpuRuntime>> {
22        iqr_impl(self, x)
23    }
24
25    fn skewness(&self, x: &Tensor<CpuRuntime>) -> Result<Tensor<CpuRuntime>> {
26        skewness_impl(self, x)
27    }
28
29    fn kurtosis(&self, x: &Tensor<CpuRuntime>) -> Result<Tensor<CpuRuntime>> {
30        kurtosis_impl(self, x)
31    }
32
33    fn zscore(&self, x: &Tensor<CpuRuntime>) -> Result<Tensor<CpuRuntime>> {
34        zscore_impl(self, x)
35    }
36
37    fn sem(&self, x: &Tensor<CpuRuntime>) -> Result<Tensor<CpuRuntime>> {
38        sem_impl(self, x)
39    }
40}
41
42#[cfg(test)]
43mod tests {
44    use super::*;
45    use crate::stats::distribution::{ContinuousDistribution, DiscreteDistribution};
46    use crate::stats::helpers::extract_scalar;
47    use numr::runtime::cpu::CpuDevice;
48
49    fn setup() -> (CpuClient, CpuDevice) {
50        let device = CpuDevice::new();
51        let client = CpuClient::new(device.clone());
52        (client, device)
53    }
54
55    #[test]
56    fn test_describe() {
57        let (client, device) = setup();
58        let data = Tensor::<CpuRuntime>::from_slice(&[1.0f64, 2.0, 3.0, 4.0, 5.0], &[5], &device);
59
60        let stats = client.describe(&data).unwrap();
61
62        assert_eq!(stats.nobs, 5);
63
64        let min_val = extract_scalar(&stats.min).unwrap();
65        let max_val = extract_scalar(&stats.max).unwrap();
66        let mean_val = extract_scalar(&stats.mean).unwrap();
67
68        assert!((min_val - 1.0).abs() < 1e-10);
69        assert!((max_val - 5.0).abs() < 1e-10);
70        assert!((mean_val - 3.0).abs() < 1e-10);
71    }
72
73    #[test]
74    fn test_percentile() {
75        let (client, device) = setup();
76        let data = Tensor::<CpuRuntime>::from_slice(&[1.0f64, 2.0, 3.0, 4.0, 5.0], &[5], &device);
77
78        let p50 = DescriptiveStatisticsAlgorithms::percentile(&client, &data, 50.0).unwrap();
79        let p50_val = extract_scalar(&p50).unwrap();
80        assert!((p50_val - 3.0).abs() < 1e-10);
81    }
82
83    #[test]
84    fn test_zscore() {
85        let (client, device) = setup();
86        let data = Tensor::<CpuRuntime>::from_slice(&[1.0f64, 2.0, 3.0, 4.0, 5.0], &[5], &device);
87
88        let z = client.zscore(&data).unwrap();
89        let z_data: Vec<f64> = z.to_vec();
90
91        // Mean of z-scores should be ~0
92        let z_mean: f64 = z_data.iter().sum::<f64>() / z_data.len() as f64;
93        assert!(z_mean.abs() < 1e-10);
94    }
95
96    #[test]
97    fn test_skewness() {
98        let (client, device) = setup();
99        // Symmetric data should have ~0 skewness
100        let data = Tensor::<CpuRuntime>::from_slice(&[1.0f64, 2.0, 3.0, 4.0, 5.0], &[5], &device);
101
102        let skew = client.skewness(&data).unwrap();
103        let skew_val = extract_scalar(&skew).unwrap();
104        assert!(skew_val.abs() < 0.1);
105    }
106
107    #[test]
108    fn test_kurtosis() {
109        let (client, device) = setup();
110        let data = Tensor::<CpuRuntime>::from_slice(
111            &[1.0f64, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0],
112            &[10],
113            &device,
114        );
115
116        let kurt = client.kurtosis(&data).unwrap();
117        let kurt_val = extract_scalar(&kurt).unwrap();
118        // Uniform-ish data has negative excess kurtosis
119        assert!(kurt_val < 0.0);
120    }
121
122    #[test]
123    fn test_iqr() {
124        let (client, device) = setup();
125        let data = Tensor::<CpuRuntime>::from_slice(&[1.0f64, 2.0, 3.0, 4.0, 5.0], &[5], &device);
126
127        let iqr = client.iqr(&data).unwrap();
128        let iqr_val = extract_scalar(&iqr).unwrap();
129        // Q3 - Q1 for [1,2,3,4,5] is 4-2 = 2
130        assert!((iqr_val - 2.0).abs() < 1e-10);
131    }
132
133    #[test]
134    fn test_sem() {
135        let (client, device) = setup();
136        let data = Tensor::<CpuRuntime>::from_slice(&[1.0f64, 2.0, 3.0, 4.0, 5.0], &[5], &device);
137
138        let sem = client.sem(&data).unwrap();
139        let sem_val = extract_scalar(&sem).unwrap();
140        // std / sqrt(n) for this data
141        assert!(sem_val > 0.0);
142    }
143
144    // ============================================================================
145    // Distribution Batch Tensor Method Tests
146    // ============================================================================
147
148    #[test]
149    fn test_normal_pdf_tensor() {
150        let (client, device) = setup();
151        let n = crate::stats::Normal::standard();
152
153        let x = Tensor::<CpuRuntime>::from_slice(&[-2.0, -1.0, 0.0, 1.0, 2.0], &[5], &device);
154        let pdf = n.pdf_tensor(&x, &client).unwrap();
155        let pdf_data: Vec<f64> = pdf.to_vec();
156
157        // Compare with scalar results
158        assert!((pdf_data[0] - n.pdf(-2.0)).abs() < 1e-10);
159        assert!((pdf_data[1] - n.pdf(-1.0)).abs() < 1e-10);
160        assert!((pdf_data[2] - n.pdf(0.0)).abs() < 1e-10);
161        assert!((pdf_data[3] - n.pdf(1.0)).abs() < 1e-10);
162        assert!((pdf_data[4] - n.pdf(2.0)).abs() < 1e-10);
163
164        // PDF at 0 for standard normal
165        assert!((pdf_data[2] - 0.3989422804014327).abs() < 1e-10);
166    }
167
168    #[test]
169    fn test_normal_cdf_tensor() {
170        let (client, device) = setup();
171        let n = crate::stats::Normal::standard();
172
173        let x = Tensor::<CpuRuntime>::from_slice(&[-2.0, 0.0, 2.0], &[3], &device);
174        let cdf = n.cdf_tensor(&x, &client).unwrap();
175        let cdf_data: Vec<f64> = cdf.to_vec();
176
177        // CDF at 0 should be 0.5
178        assert!((cdf_data[1] - 0.5).abs() < 1e-10);
179
180        // Symmetry: CDF(-2) + CDF(2) should be ~1
181        assert!((cdf_data[0] + cdf_data[2] - 1.0).abs() < 1e-10);
182
183        // Compare with scalar
184        assert!((cdf_data[0] - n.cdf(-2.0)).abs() < 1e-10);
185        assert!((cdf_data[2] - n.cdf(2.0)).abs() < 1e-10);
186    }
187
188    #[test]
189    fn test_normal_ppf_tensor() {
190        let (client, device) = setup();
191        let n = crate::stats::Normal::standard();
192
193        let p = Tensor::<CpuRuntime>::from_slice(&[0.1, 0.5, 0.9], &[3], &device);
194        let ppf = n.ppf_tensor(&p, &client).unwrap();
195        let ppf_data: Vec<f64> = ppf.to_vec();
196
197        // PPF(0.5) = 0 for standard normal
198        assert!(ppf_data[1].abs() < 1e-10);
199
200        // Symmetry: PPF(0.1) = -PPF(0.9)
201        assert!((ppf_data[0] + ppf_data[2]).abs() < 1e-10);
202
203        // Compare with scalar
204        assert!((ppf_data[0] - n.ppf(0.1).unwrap()).abs() < 1e-10);
205        assert!((ppf_data[2] - n.ppf(0.9).unwrap()).abs() < 1e-10);
206    }
207
208    #[test]
209    fn test_normal_sf_tensor() {
210        let (client, device) = setup();
211        let n = crate::stats::Normal::standard();
212
213        let x = Tensor::<CpuRuntime>::from_slice(&[0.0, 1.96], &[2], &device);
214        let sf = n.sf_tensor(&x, &client).unwrap();
215        let sf_data: Vec<f64> = sf.to_vec();
216
217        // SF(0) = 0.5
218        assert!((sf_data[0] - 0.5).abs() < 1e-10);
219
220        // SF(1.96) ~ 0.025
221        assert!((sf_data[1] - 0.025).abs() < 0.001);
222    }
223
224    #[test]
225    fn test_uniform_tensor() {
226        let (client, device) = setup();
227        let u = crate::stats::Uniform::new(0.0, 1.0).unwrap();
228
229        let x = Tensor::<CpuRuntime>::from_slice(&[0.0, 0.25, 0.5, 0.75, 1.0], &[5], &device);
230
231        // PDF should be constant = 1 over [0, 1]
232        let pdf = u.pdf_tensor(&x, &client).unwrap();
233        let pdf_data: Vec<f64> = pdf.to_vec();
234        for &p in &pdf_data {
235            assert!((p - 1.0).abs() < 1e-10);
236        }
237
238        // CDF should equal x for x in [0, 1]
239        let cdf = u.cdf_tensor(&x, &client).unwrap();
240        let cdf_data: Vec<f64> = cdf.to_vec();
241        let x_data: Vec<f64> = x.to_vec();
242        for (c, x) in cdf_data.iter().zip(x_data.iter()) {
243            assert!((c - x).abs() < 1e-10);
244        }
245    }
246
247    #[test]
248    fn test_binomial_pmf_tensor() {
249        let (client, device) = setup();
250        let b = crate::stats::Binomial::new(10, 0.5).unwrap();
251
252        let k = Tensor::<CpuRuntime>::from_slice(&[0.0, 5.0, 10.0], &[3], &device);
253        let pmf = b.pmf_tensor(&k, &client).unwrap();
254        let pmf_data: Vec<f64> = pmf.to_vec();
255
256        // Compare with scalar
257        assert!((pmf_data[0] - b.pmf(0)).abs() < 1e-10);
258        assert!((pmf_data[1] - b.pmf(5)).abs() < 1e-10);
259        assert!((pmf_data[2] - b.pmf(10)).abs() < 1e-10);
260
261        // PMF(5) should be maximum for fair coin (most likely outcome)
262        assert!(pmf_data[1] > pmf_data[0]);
263        assert!(pmf_data[1] > pmf_data[2]);
264    }
265
266    #[test]
267    fn test_binomial_cdf_tensor() {
268        let (client, device) = setup();
269        let b = crate::stats::Binomial::new(10, 0.5).unwrap();
270
271        let k = Tensor::<CpuRuntime>::from_slice(&[0.0, 5.0, 10.0], &[3], &device);
272        let cdf = b.cdf_tensor(&k, &client).unwrap();
273        let cdf_data: Vec<f64> = cdf.to_vec();
274
275        // CDF is monotonic
276        assert!(cdf_data[0] < cdf_data[1]);
277        assert!(cdf_data[1] < cdf_data[2]);
278
279        // CDF(10) = 1
280        assert!((cdf_data[2] - 1.0).abs() < 1e-10);
281
282        // Compare with scalar
283        assert!((cdf_data[0] - b.cdf(0)).abs() < 1e-10);
284        assert!((cdf_data[1] - b.cdf(5)).abs() < 1e-6);
285    }
286
287    #[test]
288    fn test_binomial_ppf_tensor() {
289        let (client, device) = setup();
290        let b = crate::stats::Binomial::new(10, 0.5).unwrap();
291
292        let p = Tensor::<CpuRuntime>::from_slice(&[0.0, 0.5, 1.0], &[3], &device);
293        let ppf = b.ppf_tensor(&p, &client).unwrap();
294        let ppf_data: Vec<f64> = ppf.to_vec();
295
296        // PPF(0) = 0
297        assert!((ppf_data[0] - 0.0).abs() < 1e-10);
298
299        // PPF(1) = n = 10
300        assert!((ppf_data[2] - 10.0).abs() < 1e-10);
301
302        // Compare with scalar
303        assert!((ppf_data[1] - b.ppf(0.5).unwrap() as f64).abs() < 1e-10);
304    }
305
306    #[test]
307    fn test_2d_tensor_batch() {
308        let (client, device) = setup();
309        let n = crate::stats::Normal::standard();
310
311        // Test with 2D tensor
312        let x =
313            Tensor::<CpuRuntime>::from_slice(&[-1.0, 0.0, 1.0, -2.0, 0.0, 2.0], &[2, 3], &device);
314
315        let pdf = n.pdf_tensor(&x, &client).unwrap();
316        let pdf_data: Vec<f64> = pdf.to_vec();
317
318        // Shape should be preserved
319        assert_eq!(pdf.shape(), &[2, 3]);
320
321        // Values should match scalar computation
322        assert!((pdf_data[0] - n.pdf(-1.0)).abs() < 1e-10);
323        assert!((pdf_data[1] - n.pdf(0.0)).abs() < 1e-10);
324        assert!((pdf_data[2] - n.pdf(1.0)).abs() < 1e-10);
325    }
326}