peroxide/statistics/stat.rs
1//! Basic statistics
2//!
3//! ## `Statistics` trait
4//!
5//! * To make generic code, there is `Statistics` trait
6//! * `mean`: just mean
7//! * `var` : variance
8//! * `sd` : standard deviation (R-like notation)
9//! * `cov` : covariance
10//! * `cor` : correlation coefficient
11//! ```rust
12//! pub trait Statistics {
13//! type Array;
14//! type Value;
15//!
16//! fn mean(&self) -> Self::Value;
17//! fn var(&self) -> Self::Value;
18//! fn sd(&self) -> Self::Value;
19//! fn cov(&self) -> Self::Array;
20//! fn cor(&self) -> Self::Array;
21//! }
22//! ```
23//!
24//! ### For `Vec<f64>`
25//!
26//! * Caution: For `Vec<f64>`, `cov` & `cor` are unimplemented (those for `Matrix`)
27//!
28//! ```rust
29//! #[macro_use]
30//! extern crate peroxide;
31//! use peroxide::fuga::*;
32//!
33//! fn main() {
34//! let a = c!(1,2,3,4,5);
35//! a.mean().print(); // 3
36//! a.var().print(); // 2.5
37//! a.sd().print(); // 1.5811388300841898
38//! }
39//! ```
40//!
41//! * But there are other functions to calculate `cov` & `cor`
42//!
43//! ```rust
44//! #[macro_use]
45//! extern crate peroxide;
46//! use peroxide::fuga::*;
47//!
48//! fn main() {
49//! let v1 = c!(1,2,3);
50//! let v2 = c!(3,2,1);
51//!
52//! cov(&v1, &v2).print(); // -0.9999999999999998
53//! cor(&v1, &v2).print(); // -0.9999999999999993
54//! }
55//! ```
56//!
57//! ### For `Matrix`
58//!
59//! * For `Matrix`, `mean, var, sd` means column operations
60//! * `cov` means covariance matrix & `cor` means also correlation coefficient matrix
61//!
62//! ```rust
63//! #[macro_use]
64//! extern crate peroxide;
65//! use peroxide::fuga::*;
66//!
67//! fn main() {
68//! let m = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
69//!
70//! m.mean().print(); // [2, 2]
71//! m.var().print(); // [1.0000, 1.0000]
72//! m.sd().print(); // [1.0000, 1.0000]
73//!
74//! m.cov().print();
75//! // c[0] c[1]
76//! // r[0] 1.0000 -1.0000
77//! // r[1] -1.0000 1.0000
78//!
79//! m.cor().print();
80//! // c[0] c[1]
81//! // r[0] 1 -1.0000
82//! // r[1] -1.0000 1
83//! }
84//! ```
85//!
86//! ### For `DataFrame`
87//!
88//! * Similar to Matrix but, `Value` is `DataFrame`
89//! * `cov` means covariance matrix.
90//!
91//! ```rust
92//! #[macro_use]
93//! extern crate peroxide;
94//! use peroxide::fuga::*;
95//!
96//! fn main() {
97//! #[cfg(feature = "dataframe")]
98//! {
99//! let mut m = DataFrame::with_header(vec!["x", "y"]);
100//! m["x"] = c!(1,2,3);
101//! m["y"] = c!(3,2,1);
102//!
103//! m.cov().print();
104//! // c[0] c[1]
105//! // r[0] 1.0000 -1.0000
106//! // r[1] -1.0000 1.0000
107//! }
108//! }
109//! ```
110//!
111//! ## Confusion Matrix
112//!
113//! * `ConfusionMatrix` is a struct to calculate confusion matrix
114//! * The reference is [here](https://en.wikipedia.org/wiki/Confusion_matrix)
115//!
116//! ### Example
117//!
118//! ```rust
119//! use peroxide::fuga::*;
120//!
121//! fn main() {
122//! let y = vec![1usize, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0];
123//! let y_hat = vec![0usize, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0];
124//! let true_val = 1usize;
125//!
126//! let cm = ConfusionMatrix::new(&y, &y_hat, true_val);
127//! cm.print();
128//! // c[0] c[1]
129//! // r[0] 6 2
130//! // r[1] 1 3
131//!
132//! // to matrix
133//! let cm_mat = cm.to_matrix();
134//!
135//! // Calculate accuracy
136//! cm.ACC().print(); // 0.75
137//!
138//! // Calculate TPR (Sensitivity or Recall)
139//! cm.TPR().print(); // 0.6666....
140//!
141//! // Calculate some metrics
142//! let metrics = cm.calc_metrics(&[ACC, TPR, TNR, F1]);
143//!
144//! // Print some metrics
145//! cm.summary(&[ACC, TPR, TNR, F1]);
146//! }
147//! ```
148
149use std::fmt;
150
151use self::QType::*;
152//use crate::structure::dataframe::*;
153use crate::structure::matrix::*;
154use crate::traits::matrix::{LinearAlgebra, MatrixTrait};
155use order_stat::kth_by;
156
157/// Statistics Trait
158///
159/// It contains `mean`, `var`, `sd`, `cov`
160pub trait Statistics {
161 type Array;
162 type Value;
163
164 fn mean(&self) -> Self::Value;
165 fn var(&self) -> Self::Value;
166 fn sd(&self) -> Self::Value;
167 fn cov(&self) -> Self::Array;
168 fn cor(&self) -> Self::Array;
169}
170
171impl Statistics for Vec<f64> {
172 type Array = Vec<f64>;
173 type Value = f64;
174
175 /// Mean
176 ///
177 /// Uses welfords online algorithm for numerically stable computation.
178 ///
179 /// # Examples
180 /// ```
181 /// #[macro_use]
182 /// extern crate peroxide;
183 /// use peroxide::fuga::*;
184 ///
185 /// fn main() {
186 /// let a = c!(1,2,3,4,5);
187 /// assert_eq!(a.mean(), 3.0);
188 /// }
189 /// ```
190 fn mean(&self) -> f64 {
191 let mut xn = 0f64;
192 let mut n = 0f64;
193
194 for x in self.iter() {
195 n += 1f64;
196 xn += (x - xn) / n;
197 }
198 xn
199 }
200
201 /// Variance
202 ///
203 /// Uses welfords online algorithm for numerically stable computation.
204 ///
205 /// # Examples
206 /// ```
207 /// #[macro_use]
208 /// extern crate peroxide;
209 /// use peroxide::fuga::*;
210 ///
211 /// fn main() {
212 /// let a = c!(1,2,3,4,5);
213 /// assert_eq!(a.var(), 2.5);
214 /// }
215 /// ```
216 fn var(&self) -> f64 {
217 let mut xn = 0f64;
218 let mut n = 0f64;
219 let mut m2n: f64 = 0f64;
220
221 for x in self.iter() {
222 n += 1f64;
223 let diff_1 = x - xn;
224 xn += diff_1 / n;
225 m2n += diff_1 * (x - xn);
226 }
227 assert_ne!(n, 1f64);
228 m2n / (n - 1f64)
229 }
230
231 /// Standard Deviation
232 ///
233 /// # Examples
234 /// ```
235 /// #[macro_use]
236 /// extern crate peroxide;
237 /// use peroxide::fuga::*;
238 ///
239 /// fn main() {
240 /// let a = c!(1,2,3);
241 /// assert!(nearly_eq(a.sd(), 1f64)); // Floating Number Error
242 /// }
243 /// ```
244 fn sd(&self) -> f64 {
245 self.var().sqrt()
246 }
247
248 fn cov(&self) -> Vec<f64> {
249 unimplemented!()
250 }
251 fn cor(&self) -> Vec<f64> {
252 unimplemented!()
253 }
254}
255
256impl Statistics for Vec<f32> {
257 type Array = Vec<f32>;
258 type Value = f32;
259
260 /// Mean
261 ///
262 /// Uses welfords online algorithm for numerically stable computation.
263 ///
264 /// # Examples
265 /// ```
266 /// #[macro_use]
267 /// extern crate peroxide;
268 /// use peroxide::fuga::*;
269 ///
270 /// fn main() {
271 /// let a = c!(1,2,3,4,5);
272 /// assert_eq!(a.mean(), 3.0);
273 /// }
274 /// ```
275 fn mean(&self) -> f32 {
276 let mut xn = 0f32;
277 let mut n = 0f32;
278
279 for x in self.iter() {
280 n += 1f32;
281 xn += (x - xn) / n;
282 }
283 xn
284 }
285
286 /// Variance
287 ///
288 /// Uses welfords online algorithm for numerically stable computation.
289 ///
290 /// # Examples
291 /// ```
292 /// #[macro_use]
293 /// extern crate peroxide;
294 /// use peroxide::fuga::*;
295 ///
296 /// fn main() {
297 /// let a = c!(1,2,3,4,5);
298 /// assert_eq!(a.var(), 2.5);
299 /// }
300 /// ```
301 fn var(&self) -> f32 {
302 let mut xn = 0f32;
303 let mut n = 0f32;
304 let mut m2n: f32 = 0f32;
305
306 for x in self.iter() {
307 n += 1f32;
308 let diff_1 = x - xn;
309 xn += diff_1 / n;
310 m2n += diff_1 * (x - xn);
311 }
312 assert_ne!(n, 1f32);
313 m2n / (n - 1f32)
314 }
315
316 /// Standard Deviation
317 ///
318 /// # Examples
319 /// ```
320 /// #[macro_use]
321 /// extern crate peroxide;
322 /// use peroxide::fuga::*;
323 ///
324 /// fn main() {
325 /// let a = c!(1,2,3);
326 /// assert!(nearly_eq(a.sd(), 1f64)); // Floating Number Error
327 /// }
328 /// ```
329 fn sd(&self) -> f32 {
330 self.var().sqrt()
331 }
332
333 fn cov(&self) -> Vec<f32> {
334 unimplemented!()
335 }
336 fn cor(&self) -> Vec<f32> {
337 unimplemented!()
338 }
339}
340
341impl Statistics for Matrix {
342 type Array = Matrix;
343 type Value = Vec<f64>;
344
345 /// Column Mean
346 ///
347 /// # Examples
348 /// ```
349 /// #[macro_use]
350 /// extern crate peroxide;
351 /// use peroxide::fuga::*;
352 ///
353 /// fn main() {
354 /// let m = matrix(c!(1,3,3,1), 2, 2, Col);
355 /// assert_eq!(m.mean(), c!(2,2));
356 /// }
357 /// ```
358 fn mean(&self) -> Vec<f64> {
359 let mut container: Vec<f64> = Vec::new();
360 let c = self.col;
361
362 for i in 0..c {
363 container.push(self.col(i).mean());
364 }
365 container
366 }
367
368 /// Column variance
369 ///
370 /// # Examples
371 /// ```
372 /// #[macro_use]
373 /// extern crate peroxide;
374 /// use peroxide::fuga::*;
375 ///
376 /// fn main() {
377 /// let m = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
378 /// assert!(nearly_eq(m.var()[0], 1));
379 /// }
380 /// ```
381 fn var(&self) -> Vec<f64> {
382 let mut container: Vec<f64> = Vec::new();
383 let c = self.col;
384
385 for i in 0..c {
386 container.push(self.col(i).var());
387 }
388 container
389 }
390
391 /// Column Standard Deviation
392 ///
393 /// # Examples
394 /// ```
395 /// #[macro_use]
396 /// extern crate peroxide;
397 /// use peroxide::fuga::*;
398 ///
399 /// fn main() {
400 /// let m = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
401 /// assert!(nearly_eq(m.sd()[0], 1));
402 /// }
403 /// ```
404 fn sd(&self) -> Vec<f64> {
405 let mut container: Vec<f64> = Vec::new();
406 let c = self.col;
407
408 for i in 0..c {
409 container.push(self.col(i).sd());
410 }
411 container
412 }
413
414 /// Covariance Matrix (Column based)
415 ///
416 /// # Examples
417 /// ```
418 /// #[macro_use]
419 /// extern crate peroxide;
420 /// use peroxide::fuga::*;
421 ///
422 /// fn main() {
423 /// let m = matrix(c!(1,2,3,3,2,1), 3, 2, Col);
424 /// println!("{}", m.cov());
425 ///
426 /// // c[0] c[1]
427 /// // r[0] 1.0000 -1.0000
428 /// // r[1] -1.0000 1.0000
429 /// }
430 /// ```
431 fn cov(&self) -> Self {
432 let c = self.col;
433
434 let mut m: Self = matrix(vec![0f64; c * c], c, c, self.shape);
435
436 for i in 0..c {
437 let m1 = self.col(i);
438 for j in 0..c {
439 let m2 = self.col(j);
440 m[(i, j)] = cov(&m1, &m2);
441 }
442 }
443 m
444 }
445
446 fn cor(&self) -> Self {
447 let c = self.col;
448
449 let mut m: Self = matrix(vec![0f64; c * c], c, c, self.shape);
450
451 for i in 0..c {
452 let m1 = self.col(i);
453 for j in 0..c {
454 let m2 = self.col(j);
455 m[(i, j)] = cor(&m1, &m2);
456 }
457 }
458 m
459 }
460}
461
462//impl Statistics for DataFrame {
463// type Array = Matrix;
464// type Value = Self;
465//
466// fn mean(&self) -> Self::Value {
467// let mut df = DataFrame::with_header(self.headers().map(|x| x.as_str()).collect());
468// for k in self.headers() {
469// df[k] = vec![self[k].mean()];
470// }
471// df
472// }
473//
474// fn var(&self) -> Self::Value {
475// let mut df = DataFrame::with_header(self.headers().map(|x| x.as_str()).collect());
476// for k in self.headers() {
477// df[k] = vec![self[k].var()];
478// }
479// df
480// }
481//
482// fn sd(&self) -> Self::Value {
483// let mut df = DataFrame::with_header(self.headers().map(|x| x.as_str()).collect());
484// for k in self.headers() {
485// df[k] = vec![self[k].sd()];
486// }
487// df
488// }
489//
490// fn cov(&self) -> Self::Array {
491// let m: Matrix = self.into();
492// m.cov()
493// }
494//
495// fn cor(&self) -> Self::Array {
496// self.to_matrix().cor()
497// }
498//}
499
500/// Covariance (to Value)
501///
502/// # Examples
503/// ```
504/// #[macro_use]
505/// extern crate peroxide;
506/// use peroxide::fuga::*;
507///
508/// fn main() {
509/// let v1 = c!(1,2,3);
510/// let v2 = c!(3,2,1);
511/// assert!(nearly_eq(cov(&v1, &v2), -1f64));
512/// }
513/// ```
514pub fn cov(v1: &Vec<f64>, v2: &Vec<f64>) -> f64 {
515 let mut ss = 0f64;
516 let mut sx = 0f64;
517 let mut sy = 0f64;
518 let mut l = 0f64;
519
520 for (x, y) in v1.into_iter().zip(v2) {
521 ss += x * y;
522 sx += *x;
523 sy += *y;
524 l += 1f64;
525 }
526 assert_ne!(l, 1f64);
527 (ss / l - (sx * sy) / (l.powf(2f64))) * l / (l - 1f64)
528}
529
530/// Pearson's correlation coefficient
531///
532/// # Examples
533/// ```
534/// #[macro_use]
535/// extern crate peroxide;
536/// use peroxide::fuga::*;
537///
538/// fn main() {
539/// let a = c!(1,2,3);
540/// let b = c!(3,2,1);
541/// assert!(nearly_eq(cor(&a, &b),-1));
542/// }
543/// ```
544pub fn cor(v1: &Vec<f64>, v2: &Vec<f64>) -> f64 {
545 cov(v1, v2) / (v1.sd() * v2.sd())
546}
547
548/// R like linear regression
549///
550/// # Examples
551/// ```
552/// #[macro_use]
553/// extern crate peroxide;
554/// use peroxide::fuga::*;
555///
556/// fn main() {
557/// let a: Matrix = c!(1,2,3,4,5).into();
558/// let noise: Matrix = Normal(0,1).sample(5).into();
559/// let b: Matrix = &a + &noise;
560/// lm(&a, &b).print();
561///
562/// // c[0]
563/// // r[0] 0.7219
564/// // r[1] 0.8058
565/// }
566/// ```
567pub fn lm(input: &Matrix, target: &Matrix) -> Matrix {
568 let mut ones = vec![1f64; input.row * input.col];
569 ones.extend(&input.data);
570 let x = matrix(ones, input.row, input.col + 1, input.shape);
571 &x.pseudo_inv() * target
572}
573
574// =============================================================================
575// Ordered Statistics (Use `order-stat`)
576// =============================================================================
577/// Trait for Ordered Statistics
578///
579/// * `median`
580/// * `quantile`
581pub trait OrderedStat {
582 type Array;
583 type Value;
584
585 fn median(&self) -> Self::Value;
586 fn quantile(&self, q: f64, qtype: QType) -> Self::Value;
587 fn quantiles(&self, q: Vec<f64>, qtype: QType) -> Self::Array;
588}
589
590/// R Quantile Type enums
591#[derive(Debug, Copy, Clone)]
592pub enum QType {
593 Type1,
594 Type2,
595 Type3,
596 Type4,
597 Type5,
598 Type6,
599 Type7,
600 Type8,
601 Type9,
602}
603
604impl OrderedStat for Vec<f64> {
605 type Array = Self;
606 type Value = f64;
607
608 fn median(&self) -> Self::Value {
609 self.quantile(0.5, Type2)
610 }
611
612 fn quantile(&self, q: f64, qtype: QType) -> Self::Value {
613 let mut m = self.clone();
614 quantile_mut(&mut m, q, qtype)
615 }
616
617 fn quantiles(&self, q: Vec<f64>, qtype: QType) -> Self::Array {
618 let mut v = vec![0f64; q.len()];
619 let mut m = self.clone();
620 for i in 0..q.len() {
621 v[i] = quantile_mut(&mut m, q[i], qtype);
622 }
623 v
624 }
625}
626
627fn quantile_mut(v: &mut [f64], q: f64, t: QType) -> f64 {
628 let l = v.len();
629 let p = 1f64 / (l as f64);
630 let k = (q / p) as usize;
631 match t {
632 Type1 => {
633 let k = if q == 0f64 {
634 0
635 } else if q == 1f64 {
636 l - 1
637 } else if q - (k as f64) * p > 0f64 {
638 k
639 } else {
640 k - 1
641 };
642 *kth_by(v, k, |x, y| x.partial_cmp(y).unwrap())
643 }
644 Type2 => {
645 if q == 0f64 {
646 let k = 0;
647 *kth_by(v, k, |x, y| x.partial_cmp(y).unwrap())
648 } else if q == 1f64 {
649 let k = l - 1;
650 *kth_by(v, k, |x, y| x.partial_cmp(y).unwrap())
651 } else if q - (k as f64) * p > 0f64 {
652 *kth_by(v, k, |x, y| x.partial_cmp(y).unwrap())
653 } else {
654 let prev = *kth_by(v, k - 1, |x, y| x.partial_cmp(y).unwrap());
655 let next = *kth_by(v, k, |x, y| x.partial_cmp(y).unwrap());
656 (prev + next) / 2f64
657 }
658 }
659 _ => unimplemented!(),
660 }
661}
662
663pub fn quantile(v: &Vec<f64>, qtype: QType) -> Vec<f64> {
664 let q_vec = vec![0.0, 0.25, 0.5, 0.75, 1.0];
665 v.quantiles(q_vec, qtype)
666}
667
668// =============================================================================
669// Confusion Matrix
670// =============================================================================
671/// Confusion Matrix
672///
673/// * `TP` : True Positive
674/// * `TN` : True Negative
675/// * `FP` : False Positive
676/// * `FN` : False Negative
677///
678/// # Examples
679/// ```
680/// use peroxide::fuga::*;
681///
682/// fn main() {
683/// let y = vec![1usize, 1, 1, 0, 0, 0];
684/// let y_hat = vec![1usize, 0, 1, 0, 0, 1];
685/// let true_val = 1usize;
686///
687/// // Create Confusion Matrix
688/// let cm = ConfusionMatrix::new(&y, &y_hat, true_val);
689///
690/// // Print
691/// cm.print();
692/// // c[0] c[1]
693/// // r[0] 2.0000 1.0000
694/// // r[1] 1.0000 2.0000
695///
696/// // To Matrix
697/// let cm_mat = cm.to_matrix();
698///
699/// // Calculate Accuracy
700/// let acc = cm.ACC();
701///
702/// // Calculate for some metrics
703/// let metrics = cm.calc_metrics(&[ACC, TPR, FPR, F1]);
704///
705/// // Print summary for some metrics
706/// cm.summary(&[ACC, TPR, FPR, F1]);
707/// }
708/// ```
709#[allow(non_snake_case)]
710#[derive(Debug, Clone, PartialEq)]
711pub struct ConfusionMatrix {
712 pub TP: usize,
713 pub TN: usize,
714 pub FP: usize,
715 pub FN: usize,
716}
717
718impl ConfusionMatrix {
719 /// Create Confusion Matrix
720 ///
721 /// # Examples
722 /// ```
723 /// use peroxide::fuga::*;
724 ///
725 /// fn main() {
726 /// let y = vec![1usize, 1, 1, 0, 0, 0];
727 /// let y_hat = vec![1usize, 0, 1, 0, 0, 1];
728 /// let true_val = 1usize;
729 ///
730 /// let true_val = 1usize;
731 /// let cm = ConfusionMatrix::new(&y, &y_hat, true_val);
732 /// cm.print();
733 /// // c[0] c[1]
734 /// // r[0] 2.0000 1.0000
735 /// // r[1] 1.0000 2.0000
736 /// }
737 /// ```
738 #[allow(non_snake_case)]
739 pub fn new<T: PartialEq + Clone + Copy>(y: &Vec<T>, y_hat: &Vec<T>, true_val: T) -> Self {
740 let mut TP = 0;
741 let mut TN = 0;
742 let mut FP = 0;
743 let mut FN = 0;
744
745 for (&y, &y_hat) in y.iter().zip(y_hat.iter()) {
746 if y == true_val && y_hat == true_val {
747 TP += 1;
748 } else if y != true_val && y_hat != true_val {
749 TN += 1;
750 } else if y != true_val && y_hat == true_val {
751 FP += 1;
752 } else if y == true_val && y_hat != true_val {
753 FN += 1;
754 }
755 }
756
757 Self { TP, TN, FP, FN }
758 }
759
760 /// Condition Positive
761 #[allow(non_snake_case)]
762 pub fn P(&self) -> usize {
763 self.TP + self.FN
764 }
765
766 /// Condition Negative
767 #[allow(non_snake_case)]
768 pub fn N(&self) -> usize {
769 self.TN + self.FP
770 }
771
772 /// True Positive Rate (Sensitivity, Recall, Hit-rate)
773 #[allow(non_snake_case)]
774 pub fn TPR(&self) -> f64 {
775 self.TP as f64 / (self.TP + self.FN) as f64
776 }
777
778 /// True Negative Rate (Specificity, Selectivity)
779 #[allow(non_snake_case)]
780 pub fn TNR(&self) -> f64 {
781 self.TN as f64 / (self.TN + self.FP) as f64
782 }
783
784 /// Positive Predictive Value (Precision)
785 #[allow(non_snake_case)]
786 pub fn PPV(&self) -> f64 {
787 self.TP as f64 / (self.TP + self.FP) as f64
788 }
789
790 /// Negative Predictive Value
791 #[allow(non_snake_case)]
792 pub fn NPV(&self) -> f64 {
793 self.TN as f64 / (self.TN + self.FN) as f64
794 }
795
796 /// False Negative Rate (Miss-rate)
797 #[allow(non_snake_case)]
798 pub fn FNR(&self) -> f64 {
799 self.FN as f64 / (self.FN + self.TP) as f64
800 }
801
802 /// False Positive Rate (Fall-out)
803 #[allow(non_snake_case)]
804 pub fn FPR(&self) -> f64 {
805 self.FP as f64 / (self.FP + self.TN) as f64
806 }
807
808 /// False Discovery Rate
809 #[allow(non_snake_case)]
810 pub fn FDR(&self) -> f64 {
811 self.FP as f64 / (self.FP + self.TP) as f64
812 }
813
814 /// False Omission Rate
815 #[allow(non_snake_case)]
816 pub fn FOR(&self) -> f64 {
817 self.FN as f64 / (self.FN + self.TN) as f64
818 }
819
820 /// Positive Likelihood Ratio
821 #[allow(non_snake_case)]
822 pub fn LR_plus(&self) -> f64 {
823 self.TPR() / self.FPR()
824 }
825
826 /// Negative Likelihood Ratio
827 #[allow(non_snake_case)]
828 pub fn LR_minus(&self) -> f64 {
829 self.FNR() / self.TNR()
830 }
831
832 /// Prevalence Threshold
833 #[allow(non_snake_case)]
834 pub fn PT(&self) -> f64 {
835 self.FPR().sqrt() / (self.TPR().sqrt() + self.FPR().sqrt())
836 }
837
838 /// Threat Score (Critical Success Index)
839 #[allow(non_snake_case)]
840 pub fn TS(&self) -> f64 {
841 self.TP as f64 / (self.TP + self.FP + self.FN) as f64
842 }
843
844 /// Prevalence
845 #[allow(non_snake_case)]
846 pub fn prevalence(&self) -> f64 {
847 self.P() as f64 / (self.P() + self.N()) as f64
848 }
849
850 /// Accuracy
851 #[allow(non_snake_case)]
852 pub fn ACC(&self) -> f64 {
853 (self.TP + self.TN) as f64 / (self.P() + self.N()) as f64
854 }
855
856 /// Balanced Accuracy
857 #[allow(non_snake_case)]
858 pub fn BA(&self) -> f64 {
859 (self.TPR() + self.TNR()) / 2f64
860 }
861
862 /// F1 Score
863 #[allow(non_snake_case)]
864 pub fn F1(&self) -> f64 {
865 2f64 * self.TP as f64 / (2f64 * self.TP as f64 + self.FP as f64 + self.FN as f64)
866 }
867
868 /// Matthews Correlation Coefficient (Phi Coefficient)
869 #[allow(non_snake_case)]
870 pub fn MCC(&self) -> f64 {
871 let a = self.TP as f64;
872 let b = self.FP as f64;
873 let c = self.FN as f64;
874 let d = self.TN as f64;
875 (a * d - b * c) / ((a + b) * (a + c) * (d + b) * (d + c)).sqrt()
876 }
877
878 /// Fowlkes-Mallows Index
879 #[allow(non_snake_case)]
880 pub fn FM(&self) -> f64 {
881 (self.PPV() * self.TPR()).sqrt()
882 }
883
884 /// Bookmaker Informedness (Informedness)
885 #[allow(non_snake_case)]
886 pub fn BM(&self) -> f64 {
887 self.TPR() + self.TNR() - 1f64
888 }
889
890 /// Markedness (deltaP)
891 #[allow(non_snake_case)]
892 pub fn MK(&self) -> f64 {
893 self.PPV() + self.NPV() - 1f64
894 }
895
896 /// Diagnostic Odds Ratio
897 #[allow(non_snake_case)]
898 pub fn DOR(&self) -> f64 {
899 self.LR_plus() / self.LR_minus()
900 }
901
902 /// To Matrix
903 pub fn to_matrix(&self) -> Matrix {
904 let mut m = matrix(vec![0f64; 4], 2, 2, Row);
905 m[(0, 0)] = self.TP as f64;
906 m[(0, 1)] = self.FP as f64;
907 m[(1, 0)] = self.FN as f64;
908 m[(1, 1)] = self.TN as f64;
909 m
910 }
911
912 /// Calculate a specific metric
913 pub fn calc_metric(&self, metric: Metric) -> f64 {
914 match metric {
915 Metric::TPR => self.TPR(),
916 Metric::TNR => self.TNR(),
917 Metric::PPV => self.PPV(),
918 Metric::NPV => self.NPV(),
919 Metric::FNR => self.FNR(),
920 Metric::FPR => self.FPR(),
921 Metric::FDR => self.FDR(),
922 Metric::FOR => self.FOR(),
923 Metric::LR_plus => self.LR_plus(),
924 Metric::LR_minus => self.LR_minus(),
925 Metric::PT => self.PT(),
926 Metric::TS => self.TS(),
927 Metric::prevalence => self.prevalence(),
928 Metric::ACC => self.ACC(),
929 Metric::BA => self.BA(),
930 Metric::F1 => self.F1(),
931 Metric::MCC => self.MCC(),
932 Metric::FM => self.FM(),
933 Metric::BM => self.BM(),
934 Metric::MK => self.MK(),
935 Metric::DOR => self.DOR(),
936 }
937 }
938
939 /// Calculate for some metrics
940 pub fn calc_metrics(&self, metrics: &[Metric]) -> Vec<f64> {
941 metrics.iter().map(|m| self.calc_metric(*m)).collect()
942 }
943
944 /// Summarize some metrics
945 pub fn summary(&self, metrics: &[Metric]) {
946 let width = metrics.iter().fold(0, |acc, m| {
947 if m.to_string().len() > acc {
948 m.to_string().len()
949 } else {
950 acc
951 }
952 });
953 println!("============================================================");
954 println!("Summary of metrics");
955 println!("============================================================");
956 for m in metrics {
957 println!("{:width$}:\t{:.4}", m.to_string(), self.calc_metric(*m));
958 }
959 println!("============================================================");
960 }
961}
962
963#[allow(non_camel_case_types)]
964#[derive(Debug, Clone, Copy)]
965pub enum Metric {
966 TPR,
967 TNR,
968 PPV,
969 NPV,
970 FNR,
971 FPR,
972 FDR,
973 FOR,
974 LR_plus,
975 LR_minus,
976 PT,
977 TS,
978 prevalence,
979 ACC,
980 BA,
981 F1,
982 MCC,
983 FM,
984 BM,
985 MK,
986 DOR,
987}
988
989impl fmt::Display for Metric {
990 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
991 match self {
992 Metric::TPR => write!(f, "TPR"),
993 Metric::TNR => write!(f, "TNR"),
994 Metric::PPV => write!(f, "PPV"),
995 Metric::NPV => write!(f, "NPV"),
996 Metric::FNR => write!(f, "FNR"),
997 Metric::FPR => write!(f, "FPR"),
998 Metric::FDR => write!(f, "FDR"),
999 Metric::FOR => write!(f, "FOR"),
1000 Metric::LR_plus => write!(f, "LR+"),
1001 Metric::LR_minus => write!(f, "LR-"),
1002 Metric::PT => write!(f, "PT"),
1003 Metric::TS => write!(f, "TS"),
1004 Metric::prevalence => write!(f, "prevalence"),
1005 Metric::ACC => write!(f, "ACC"),
1006 Metric::BA => write!(f, "BA"),
1007 Metric::F1 => write!(f, "F1"),
1008 Metric::MCC => write!(f, "MCC"),
1009 Metric::FM => write!(f, "FM"),
1010 Metric::BM => write!(f, "BM"),
1011 Metric::MK => write!(f, "MK"),
1012 Metric::DOR => write!(f, "DOR"),
1013 }
1014 }
1015}