probability/distribution/
lognormal.rs

1use alloc::{vec, vec::Vec};
2#[allow(unused_imports)]
3use special::Primitive;
4
5use distribution::{self, Gaussian};
6use source::Source;
7
8/// A lognormal distribution.
9#[derive(Clone, Copy, Debug)]
10pub struct Lognormal {
11    mu: f64,
12    sigma: f64,
13    gaussian: Gaussian,
14}
15
16impl Lognormal {
17    /// Create a lognormal distribution with location `mu` and scale `sigma`.
18    ///
19    /// It should hold that `sigma > 0`.
20    #[inline]
21    pub fn new(mu: f64, sigma: f64) -> Self {
22        should!(sigma > 0.0);
23        Lognormal {
24            mu,
25            sigma,
26            gaussian: Gaussian::new(mu, sigma),
27        }
28    }
29
30    /// Return the location parameter.
31    #[inline(always)]
32    pub fn mu(&self) -> f64 {
33        self.mu
34    }
35
36    /// Return the scale parameter.
37    #[inline(always)]
38    pub fn sigma(&self) -> f64 {
39        self.sigma
40    }
41}
42
43impl Default for Lognormal {
44    #[inline]
45    fn default() -> Self {
46        Lognormal::new(0.0, 1.0)
47    }
48}
49
50impl distribution::Continuous for Lognormal {
51    fn density(&self, x: f64) -> f64 {
52        use core::f64::consts::PI;
53        if x <= 0.0 {
54            0.0
55        } else {
56            let &Lognormal { mu, sigma, .. } = self;
57            (-(x.ln() - mu).powi(2) / (2.0 * sigma * sigma)).exp() / (x * sigma * (2.0 * PI).sqrt())
58        }
59    }
60}
61
62impl distribution::Distribution for Lognormal {
63    type Value = f64;
64
65    fn distribution(&self, x: f64) -> f64 {
66        if x <= 0.0 {
67            0.0
68        } else {
69            self.gaussian.distribution(x.ln())
70        }
71    }
72}
73
74impl distribution::Entropy for Lognormal {
75    #[inline]
76    fn entropy(&self) -> f64 {
77        use core::f64::consts::PI;
78        (self.sigma * (self.mu + 0.5).exp() * (2.0 * PI).sqrt()).ln()
79    }
80}
81
82impl distribution::Inverse for Lognormal {
83    fn inverse(&self, p: f64) -> f64 {
84        self.gaussian.inverse(p).exp()
85    }
86}
87
88impl distribution::Kurtosis for Lognormal {
89    #[inline]
90    fn kurtosis(&self) -> f64 {
91        let s2 = self.sigma * self.sigma;
92        (4.0 * s2).exp() + 2.0 * (3.0 * s2).exp() + 3.0 * (2.0 * s2).exp() - 6.0
93    }
94}
95
96impl distribution::Mean for Lognormal {
97    #[inline]
98    fn mean(&self) -> f64 {
99        (self.mu + self.sigma * self.sigma / 2.0).exp()
100    }
101}
102
103impl distribution::Median for Lognormal {
104    #[inline]
105    fn median(&self) -> f64 {
106        self.mu.exp()
107    }
108}
109
110impl distribution::Modes for Lognormal {
111    #[inline]
112    fn modes(&self) -> Vec<f64> {
113        vec![(self.mu - self.sigma * self.sigma).exp()]
114    }
115}
116
117impl distribution::Sample for Lognormal {
118    #[inline]
119    fn sample<S>(&self, source: &mut S) -> f64
120    where
121        S: Source,
122    {
123        self.gaussian.sample(source).exp()
124    }
125}
126
127impl distribution::Skewness for Lognormal {
128    #[inline]
129    fn skewness(&self) -> f64 {
130        let es2 = (self.sigma * self.sigma).exp();
131        (es2 - 1.0).sqrt() * (2.0 + es2)
132    }
133}
134
135impl distribution::Variance for Lognormal {
136    #[inline]
137    fn variance(&self) -> f64 {
138        let s2 = self.sigma * self.sigma;
139        (s2.exp() - 1.0) * (2.0 * self.mu + s2).exp()
140    }
141}
142
143#[cfg(test)]
144mod tests {
145    use alloc::{vec, vec::Vec};
146    use assert;
147    use prelude::*;
148
149    macro_rules! new(
150        ($mu:expr, $sigma:expr) => (Lognormal::new($mu, $sigma));
151    );
152
153    #[test]
154    fn density() {
155        let d = new!(1.0, 2.0);
156        let x = vec![0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0];
157        let p = vec![
158            0.0000000000000000e+00,
159            2.7879404629273086e-01,
160            1.7603266338214976e-01,
161            1.2723305581441105e-01,
162            9.8568580344013113e-02,
163            7.9718599555316239e-02,
164            6.6409606924506773e-02,
165            5.6538422820400766e-02,
166            4.8946227003151078e-02,
167            4.2941143217487855e-02,
168            3.8084403129689012e-02,
169        ];
170
171        assert::close(
172            &x.iter().map(|&x| d.density(x)).collect::<Vec<_>>(),
173            &p,
174            1e-15,
175        );
176    }
177
178    #[test]
179    fn distribution() {
180        let d = new!(1.0, 2.0);
181        let x = vec![0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0];
182        let p = vec![
183            0.0000000000000000e+00,
184            1.9861641975736130e-01,
185            3.0853753872598694e-01,
186            3.8313116661630492e-01,
187            4.3903100974768944e-01,
188            4.8330729072740009e-01,
189            5.1966233849751675e-01,
190            5.5028502097208276e-01,
191            5.7657814823924480e-01,
192            5.9949442394950303e-01,
193            6.1970989457732906e-01,
194        ];
195
196        assert::close(
197            &x.iter().map(|&x| d.distribution(x)).collect::<Vec<_>>(),
198            &p,
199            1e-15,
200        );
201    }
202
203    #[test]
204    fn entropy() {
205        use core::f64::consts::PI;
206        assert_eq!(new!(-0.5, 1.0 / (2.0 * PI).sqrt()).entropy(), 0.0);
207    }
208
209    #[test]
210    fn inverse() {
211        use core::f64::INFINITY;
212        let d = new!(1.0, 2.0);
213        let p = vec![
214            0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65,
215            0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00,
216        ];
217        let x = vec![
218            0.0000000000000000e+00,
219            1.0129611155505908e-01,
220            2.0948500212405705e-01,
221            3.4202659595680435e-01,
222            5.0497696371871126e-01,
223            7.0540759070071157e-01,
224            9.5237060839269883e-01,
225            1.2577935903399797e+00,
226            1.6377212497125082e+00,
227            2.1142017250556107e+00,
228            2.7182818284590451e+00,
229            3.4949626666945868e+00,
230            4.5117910634839467e+00,
231            5.8746173900706555e+00,
232            7.7585931714136498e+00,
233            1.0474874663016855e+01,
234            1.4632461735515125e+01,
235            2.1603747153814467e+01,
236            3.5272482631261830e+01,
237            7.2945110977081981e+01,
238            INFINITY,
239        ];
240
241        assert::close(
242            &p.iter().map(|&p| d.inverse(p)).collect::<Vec<_>>(),
243            &x,
244            1e-12,
245        );
246    }
247
248    #[test]
249    fn kurtosis() {
250        assert::close(new!(0.0, 1.0).kurtosis(), 1.1093639217631153e+02, 1e-15);
251    }
252
253    #[test]
254    fn mean() {
255        assert_eq!(new!(-2.0, 2.0).mean(), 1.0);
256    }
257
258    #[test]
259    fn median() {
260        assert_eq!(new!(0.0, 1.0).median(), 1.0);
261    }
262
263    #[test]
264    fn modes() {
265        assert_eq!(new!(1.0, 1.0).modes(), vec![1.0]);
266    }
267
268    #[test]
269    fn skewness() {
270        assert!(4.0 - new!(0.0, 2f64.ln().sqrt()).skewness() < 1e-10);
271    }
272
273    #[test]
274    fn variance() {
275        assert!(2.0 - new!(0.0, 2f64.ln().sqrt()).variance() < 1e-10);
276    }
277
278    #[test]
279    fn deviation() {
280        assert!(2f64.sqrt() - new!(0.0, 2f64.ln().sqrt()).variance() < 1e-10);
281    }
282}
probability/distribution/lognormal.rs

probability/distribution/
lognormal.rs
