ferric 0.2.1

A Probablistic Programming Language with a declarative syntax for random variables.
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

// Copyright 2022 The Ferric AI Project Developers

use rand_distr::Distribution as Distribution2;
use rand_distr::Exp as Exp2;

use crate::distributions::Distribution;
use rand::Rng;

/// Exponential distribution over non-negative reals.
///
/// The PDF is
///
/// $$p(x \mid \lambda) = \lambda e^{-\lambda x}, \quad x \geq 0$$
///
/// where $\lambda > 0$ is the rate parameter (mean $= 1/\lambda$).
///
/// See [Exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution)
/// on Wikipedia for further details.
///
/// # Examples
///
/// ```
/// use ferric::distributions::{Distribution, Exponential};
/// use rand::thread_rng;
///
/// let dist = Exponential::new(2.0).unwrap();
/// let x: f64 = dist.sample(&mut thread_rng());
/// println!("sample = {:.4}", x);
/// ```
pub struct Exponential {
    rate: f64,
}

impl Exponential {
    /// Construct an Exponential distribution with rate parameter `rate` ($\lambda$).
    ///
    /// # Errors
    ///
    /// Returns `Err` if `rate` is not strictly positive.
    pub fn new(rate: f64) -> Result<Exponential, String> {
        if rate <= 0.0 {
            Err(format!(
                "Exponential: illegal rate `{}` should be greater than 0",
                rate
            ))
        } else {
            Ok(Exponential { rate })
        }
    }
}

impl<R: Rng + ?Sized> Distribution<R> for Exponential {
    type Domain = f64;

    fn sample(&self, rng: &mut R) -> f64 {
        Exp2::new(self.rate).unwrap().sample(rng)
    }

    /// Returns $\ln\lambda - \lambda x$ for $x \geq 0$, or $-\infty$ otherwise.
    fn log_prob(&self, x: &f64) -> f64 {
        if *x < 0.0 {
            f64::NEG_INFINITY
        } else {
            self.rate.ln() - self.rate * x
        }
    }

    fn is_discrete(&self) -> bool {
        false
    }
}

impl std::fmt::Display for Exponential {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        write!(f, "Exponential {{ rate = {} }}", self.rate)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use rand::rngs::ThreadRng;
    use rand::thread_rng;

    #[test]
    fn exponential_sample() {
        let mut rng = thread_rng();
        let rate = 2.0f64;
        let dist = Exponential::new(rate).unwrap();
        println!("dist = {}", dist);
        let mut total = 0f64;
        let trials = 10000;
        for _ in 0..trials {
            total += dist.sample(&mut rng);
        }
        let empirical_mean = total / (trials as f64);
        let expected_mean = 1.0 / rate;
        let expected_std = 1.0 / rate;
        let err = 5.0 * expected_std / (trials as f64).sqrt();
        assert!((empirical_mean - expected_mean).abs() < err);
    }

    #[test]
    fn exponential_log_prob() {
        // Exp(1): log_prob(1) = ln(1) - 1*1 = -1
        let dist = Exponential::new(1.0).unwrap();
        let lp = <Exponential as Distribution<ThreadRng>>::log_prob(&dist, &1.0);
        assert!((lp - (-1.0f64)).abs() < 1e-10);
        // outside support
        let lp_out = <Exponential as Distribution<ThreadRng>>::log_prob(&dist, &-0.1);
        assert_eq!(lp_out, f64::NEG_INFINITY);
        assert!(!<Exponential as Distribution<ThreadRng>>::is_discrete(
            &dist
        ));
    }

    #[test]
    #[should_panic]
    fn exponential_zero_rate() {
        Exponential::new(0.0).unwrap();
    }

    #[test]
    #[should_panic]
    fn exponential_negative_rate() {
        Exponential::new(-1.0).unwrap();
    }
}