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

use crate::stats::LogProb;

pub mod evidence {
    /// Scale of evidence as defined by
    /// [Kass and Raftery 1995](http://www.andrew.cmu.edu/user/kk3n/simplicity/KassRaftery1995.pdf).
    #[derive(
        Display,
        Debug,
        Clone,
        Copy,
        PartialEq,
        Eq,
        PartialOrd,
        Ord,
        Serialize,
        Deserialize,
        EnumString,
        EnumIter,
        IntoStaticStr,
        EnumVariantNames,
    )]
    pub enum KassRaftery {
        #[strum(serialize = "none")]
        None,
        #[strum(serialize = "barely")]
        Barely,
        #[strum(serialize = "positive")]
        Positive,
        #[strum(serialize = "strong")]
        Strong,
        #[strum(serialize = "very-strong")]
        VeryStrong,
    }
}

custom_derive! {
    /// A newtype for Bayes factors.
    #[derive(NewtypeFrom, NewtypeDeref, NewtypeAdd(*), NewtypeSub(*), Default, Copy, Clone, PartialEq, PartialOrd, Debug)]
    pub struct BayesFactor(pub f64);
}

impl BayesFactor {
    /// Calculate Bayes factor from given probabilities.
    pub fn new(a: LogProb, b: LogProb) -> Self {
        BayesFactor((a - b).exp())
    }

    /// Calculate strength of evidence as defined by
    /// [Kass and Raftery 1995](http://www.andrew.cmu.edu/user/kk3n/simplicity/KassRaftery1995.pdf).
    pub fn evidence_kass_raftery(&self) -> evidence::KassRaftery {
        let k = **self;
        if k <= 1.0 {
            evidence::KassRaftery::None
        } else if k <= 3.0 {
            evidence::KassRaftery::Barely
        } else if k <= 20.0 {
            evidence::KassRaftery::Positive
        } else if k <= 150.0 {
            evidence::KassRaftery::Strong
        } else {
            evidence::KassRaftery::VeryStrong
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_bayes_factor() {
        let bf = BayesFactor::new(LogProb(0.5_f64.ln()), LogProb(0.1_f64.ln()));
        assert_relative_eq!(*bf, 5.0, epsilon = 1e-9);
        assert_eq!(bf.evidence_kass_raftery(), evidence::KassRaftery::Positive);
    }
}