rv/dist/
lognormal.rs

1//! Log Normal Distribution over x in (0, ∞)
2#[cfg(feature = "serde1")]
3use serde::{Deserialize, Serialize};
4
5use crate::consts::*;
6use crate::impl_display;
7use crate::traits::*;
8use rand::Rng;
9use special::Error as _;
10use std::f64::consts::SQRT_2;
11use std::fmt;
12
13/// [LogNormal Distribution](https://en.wikipedia.org/wiki/Log-normal_distribution)
14/// If x ~ Normal(μ, σ), then e^x ~ LogNormal(μ, σ).
15#[derive(Debug, Clone, PartialEq)]
16#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
17pub struct LogNormal {
18    /// log scale mean
19    mu: f64,
20    /// log scale standard deviation
21    sigma: f64,
22}
23
24pub struct LogNormalParameters {
25    pub mu: f64,
26    pub sigma: f64,
27}
28
29impl Parameterized for LogNormal {
30    type Parameters = LogNormalParameters;
31
32    fn emit_params(&self) -> Self::Parameters {
33        Self::Parameters {
34            mu: self.mu(),
35            sigma: self.sigma(),
36        }
37    }
38
39    fn from_params(params: Self::Parameters) -> Self {
40        Self::new_unchecked(params.mu, params.sigma)
41    }
42}
43
44#[derive(Debug, Clone, PartialEq)]
45#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
46#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
47pub enum LogNormalError {
48    /// The mu parameter is infinite or NaN
49    MuNotFinite { mu: f64 },
50    /// The sigma parameter is less than or equal to zero
51    SigmaTooLow { sigma: f64 },
52    /// The sigma parameter is infinite or NaN
53    SigmaNotFinite { sigma: f64 },
54}
55
56impl LogNormal {
57    /// Create a new LogNormal distribution
58    ///
59    /// # Arguments
60    /// - mu: log scale mean
61    /// - sigma: log scale standard deviation
62    #[inline]
63    pub fn new(mu: f64, sigma: f64) -> Result<Self, LogNormalError> {
64        if !mu.is_finite() {
65            Err(LogNormalError::MuNotFinite { mu })
66        } else if sigma <= 0.0 {
67            Err(LogNormalError::SigmaTooLow { sigma })
68        } else if !sigma.is_finite() {
69            Err(LogNormalError::SigmaNotFinite { sigma })
70        } else {
71            Ok(LogNormal { mu, sigma })
72        }
73    }
74
75    /// Creates a new LogNormal without checking whether the parameters are
76    /// valid.
77    #[inline]
78    pub fn new_unchecked(mu: f64, sigma: f64) -> Self {
79        LogNormal { mu, sigma }
80    }
81
82    /// LogNormal(0, 1)
83    ///
84    /// # Example
85    ///
86    /// ```rust
87    /// # use rv::dist::LogNormal;
88    /// let lognormal = LogNormal::standard();
89    /// assert_eq!(lognormal, LogNormal::new(0.0, 1.0).unwrap());
90    /// ```
91    #[inline]
92    pub fn standard() -> Self {
93        LogNormal {
94            mu: 0.0,
95            sigma: 1.0,
96        }
97    }
98
99    /// Get the mu parameter
100    ///
101    /// # Example
102    ///
103    /// ```rust
104    /// # use rv::dist::LogNormal;
105    /// let lognormal = LogNormal::new(-1.0, 2.0).unwrap();
106    /// assert_eq!(lognormal.mu(), -1.0);
107    /// ```
108    #[inline]
109    pub fn mu(&self) -> f64 {
110        self.mu
111    }
112
113    /// Set the value of mu
114    ///
115    /// # Example
116    ///
117    /// ```rust
118    /// # use rv::dist::LogNormal;
119    /// let mut lognormal = LogNormal::new(2.0, 1.5).unwrap();
120    /// assert_eq!(lognormal.mu(), 2.0);
121    ///
122    /// lognormal.set_mu(1.3).unwrap();
123    /// assert_eq!(lognormal.mu(), 1.3);
124    /// ```
125    ///
126    /// Will error for invalid values
127    ///
128    /// ```rust
129    /// # use rv::dist::LogNormal;
130    /// # let mut lognormal = LogNormal::new(2.0, 1.5).unwrap();
131    /// assert!(lognormal.set_mu(1.3).is_ok());
132    /// assert!(lognormal.set_mu(f64::NEG_INFINITY).is_err());
133    /// assert!(lognormal.set_mu(f64::INFINITY).is_err());
134    /// assert!(lognormal.set_mu(f64::NAN).is_err());
135    /// ```
136    #[inline]
137    pub fn set_mu(&mut self, mu: f64) -> Result<(), LogNormalError> {
138        if !mu.is_finite() {
139            Err(LogNormalError::MuNotFinite { mu })
140        } else {
141            self.set_mu_unchecked(mu);
142            Ok(())
143        }
144    }
145
146    /// Set the value of mu without input validation
147    #[inline]
148    pub fn set_mu_unchecked(&mut self, mu: f64) {
149        self.mu = mu;
150    }
151
152    /// Get the sigma parameter
153    ///
154    /// # Example
155    ///
156    /// ```rust
157    /// # use rv::dist::LogNormal;
158    /// let lognormal = LogNormal::new(-1.0, 2.0).unwrap();
159    /// assert_eq!(lognormal.sigma(), 2.0);
160    /// ```
161    #[inline]
162    pub fn sigma(&self) -> f64 {
163        self.sigma
164    }
165
166    /// Set the value of sigma
167    ///
168    /// # Example
169    ///
170    /// ```rust
171    /// # use rv::dist::LogNormal;
172    /// let mut lognormal = LogNormal::standard();
173    /// assert_eq!(lognormal.sigma(), 1.0);
174    ///
175    /// lognormal.set_sigma(2.3).unwrap();
176    /// assert_eq!(lognormal.sigma(), 2.3);
177    /// ```
178    ///
179    /// Will error for invalid values
180    ///
181    /// ```rust
182    /// # use rv::dist::LogNormal;
183    /// # let mut lognormal = LogNormal::standard();
184    /// assert!(lognormal.set_sigma(2.3).is_ok());
185    /// assert!(lognormal.set_sigma(0.0).is_err());
186    /// assert!(lognormal.set_sigma(-1.0).is_err());
187    /// assert!(lognormal.set_sigma(f64::INFINITY).is_err());
188    /// assert!(lognormal.set_sigma(f64::NEG_INFINITY).is_err());
189    /// assert!(lognormal.set_sigma(f64::NAN).is_err());
190    /// ```
191    #[inline]
192    pub fn set_sigma(&mut self, sigma: f64) -> Result<(), LogNormalError> {
193        if sigma <= 0.0 {
194            Err(LogNormalError::SigmaTooLow { sigma })
195        } else if !sigma.is_finite() {
196            Err(LogNormalError::SigmaNotFinite { sigma })
197        } else {
198            self.set_sigma_unchecked(sigma);
199            Ok(())
200        }
201    }
202
203    /// Set the value of sigma
204    #[inline]
205    pub fn set_sigma_unchecked(&mut self, sigma: f64) {
206        self.sigma = sigma;
207    }
208}
209
210impl Default for LogNormal {
211    fn default() -> Self {
212        LogNormal::standard()
213    }
214}
215
216impl From<&LogNormal> for String {
217    fn from(lognorm: &LogNormal) -> String {
218        format!("LogNormal(μ: {}, σ: {})", lognorm.mu, lognorm.sigma)
219    }
220}
221
222impl_display!(LogNormal);
223
224macro_rules! impl_traits {
225    ($kind: ty) => {
226        impl HasDensity<$kind> for LogNormal {
227            fn ln_f(&self, x: &$kind) -> f64 {
228                // TODO: cache ln(sigma)
229                let xk = f64::from(*x);
230                let xk_ln = xk.ln();
231                let d = (xk_ln - self.mu) / self.sigma;
232                (0.5 * d).mul_add(-d, -xk_ln - self.sigma.ln() - HALF_LN_2PI)
233            }
234        }
235
236        impl Sampleable<$kind> for LogNormal {
237            fn draw<R: Rng>(&self, rng: &mut R) -> $kind {
238                let g =
239                    rand_distr::LogNormal::new(self.mu, self.sigma).unwrap();
240                rng.sample(g) as $kind
241            }
242
243            fn sample<R: Rng>(&self, n: usize, rng: &mut R) -> Vec<$kind> {
244                let g =
245                    rand_distr::LogNormal::new(self.mu, self.sigma).unwrap();
246                (0..n).map(|_| rng.sample(g) as $kind).collect()
247            }
248        }
249
250        impl ContinuousDistr<$kind> for LogNormal {}
251
252        impl Support<$kind> for LogNormal {
253            fn supports(&self, x: &$kind) -> bool {
254                *x > 0.0 && x.is_finite()
255            }
256        }
257
258        impl Cdf<$kind> for LogNormal {
259            fn cdf(&self, x: &$kind) -> f64 {
260                let xk = f64::from(*x);
261                0.5_f64.mul_add(
262                    ((xk.ln() - self.mu) / (SQRT_2 * self.sigma)).error(),
263                    0.5,
264                )
265            }
266        }
267
268        impl InverseCdf<$kind> for LogNormal {
269            fn invcdf(&self, p: f64) -> $kind {
270                SQRT_2
271                    .mul_add(
272                        self.sigma * 2.0_f64.mul_add(p, -1.0).inv_error(),
273                        self.mu,
274                    )
275                    .exp() as $kind
276            }
277        }
278
279        impl Mean<$kind> for LogNormal {
280            fn mean(&self) -> Option<$kind> {
281                Some((self.mu + self.sigma * self.sigma / 2.0).exp() as $kind)
282            }
283        }
284
285        impl Median<$kind> for LogNormal {
286            fn median(&self) -> Option<$kind> {
287                Some(self.mu.exp() as $kind)
288            }
289        }
290
291        impl Mode<$kind> for LogNormal {
292            fn mode(&self) -> Option<$kind> {
293                Some(self.sigma.mul_add(-self.sigma, self.mu) as $kind)
294            }
295        }
296    };
297}
298
299impl Variance<f64> for LogNormal {
300    fn variance(&self) -> Option<f64> {
301        Some(
302            (self.sigma * self.sigma).exp_m1()
303                * self.sigma.mul_add(self.sigma, 2.0 * self.mu).exp(),
304        )
305    }
306}
307
308impl Entropy for LogNormal {
309    fn entropy(&self) -> f64 {
310        (self.mu + 0.5) + self.sigma.ln() + HALF_LN_2PI
311    }
312}
313
314impl Skewness for LogNormal {
315    fn skewness(&self) -> Option<f64> {
316        let e_sigma_2 = (self.sigma * self.sigma).exp();
317        Some((e_sigma_2 + 2.0) * (e_sigma_2 - 1.0).sqrt())
318    }
319}
320
321impl Kurtosis for LogNormal {
322    fn kurtosis(&self) -> Option<f64> {
323        let s2 = self.sigma * self.sigma;
324        Some(
325            3.0_f64.mul_add(
326                (2.0 * s2).exp(),
327                (3.0 * s2).exp().mul_add(2.0, (4.0 * s2).exp()),
328            ) - 6.0,
329        )
330    }
331}
332
333impl_traits!(f32);
334impl_traits!(f64);
335
336impl std::error::Error for LogNormalError {}
337
338impl fmt::Display for LogNormalError {
339    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
340        match self {
341            Self::MuNotFinite { mu } => write!(f, "non-finite mu: {}", mu),
342            Self::SigmaTooLow { sigma } => {
343                write!(f, "sigma ({}) must be greater than zero", sigma)
344            }
345            Self::SigmaNotFinite { sigma } => {
346                write!(f, "non-finite sigma: {}", sigma)
347            }
348        }
349    }
350}
351
352#[cfg(test)]
353mod tests {
354    use super::*;
355    use crate::test_basic_impls;
356    use proptest::prelude::*;
357    use std::f64;
358
359    const TOL: f64 = 1E-12;
360
361    test_basic_impls!(f64, LogNormal);
362
363    #[test]
364    fn new() {
365        let lognorm = LogNormal::new(1.2, 3.0).unwrap();
366        assert::close(lognorm.mu, 1.2, TOL);
367        assert::close(lognorm.sigma, 3.0, TOL);
368    }
369
370    #[test]
371    fn mean() {
372        let mu = 3.0;
373        let sigma = 2.0;
374        let mean: f64 = LogNormal::new(mu, sigma).unwrap().mean().unwrap();
375        assert::close(mean, 5.0_f64.exp(), TOL);
376    }
377
378    #[test]
379    fn median_should_be_exp_mu() {
380        let mu = 3.4;
381        let median: f64 = LogNormal::new(mu, 0.5).unwrap().median().unwrap();
382        assert::close(median, mu.exp(), TOL);
383    }
384
385    #[test]
386    fn mode() {
387        let mode: f64 = LogNormal::new(4.0, 2.0).unwrap().mode().unwrap();
388        assert::close(mode, 0.0, TOL);
389    }
390
391    #[test]
392    fn variance() {
393        let lognorm_1 = LogNormal::new(3.4, 1.0).unwrap();
394        let lognorm_2 = LogNormal::new(1.0, 3.0).unwrap();
395        assert::close(
396            lognorm_1.variance().unwrap(),
397            1.0_f64.exp_m1() * 7.8_f64.exp(),
398            TOL,
399        );
400        assert::close(
401            lognorm_2.variance().unwrap(),
402            9.0_f64.exp_m1() * 11.0_f64.exp(),
403            TOL,
404        );
405    }
406
407    #[test]
408    fn draws_should_be_finite() {
409        let mut rng = rand::thread_rng();
410        let lognorm = LogNormal::standard();
411        for _ in 0..100 {
412            let x: f64 = lognorm.draw(&mut rng);
413            assert!(x.is_finite())
414        }
415    }
416
417    #[test]
418    fn sample_length() {
419        let mut rng = rand::thread_rng();
420        let lognorm = LogNormal::standard();
421        let xs: Vec<f64> = lognorm.sample(10, &mut rng);
422        assert_eq!(xs.len(), 10);
423    }
424
425    #[test]
426    fn standard_ln_pdf_at_one() {
427        let lognorm = LogNormal::standard();
428        assert::close(lognorm.ln_pdf(&1.0_f64), -0.918_938_533_204_672_7, TOL);
429    }
430
431    #[test]
432    fn standard_ln_pdf_at_e() {
433        let lognorm = LogNormal::standard();
434        assert::close(
435            lognorm.ln_pdf(&f64::consts::E),
436            -2.418_938_533_204_672_7,
437            TOL,
438        );
439    }
440
441    #[test]
442    fn should_contain_positive_finite_values() {
443        let lognorm = LogNormal::standard();
444        assert!(lognorm.supports(&1E-8_f32));
445        assert!(lognorm.supports(&10E8_f64));
446    }
447
448    #[test]
449    fn should_not_contain_negative_or_zero() {
450        let lognorm = LogNormal::standard();
451        assert!(!lognorm.supports(&-1.0_f64));
452        assert!(!lognorm.supports(&0.0_f64));
453    }
454
455    #[test]
456    fn should_not_contain_nan() {
457        let lognorm = LogNormal::standard();
458        assert!(!lognorm.supports(&f64::NAN));
459    }
460
461    #[test]
462    fn should_not_contain_positive_or_negative_infinity() {
463        let lognorm = LogNormal::standard();
464        assert!(!lognorm.supports(&f64::INFINITY));
465        assert!(!lognorm.supports(&f64::NEG_INFINITY));
466    }
467
468    #[test]
469    fn skewness() {
470        let lognorm = LogNormal::new(-1.2, 3.4).unwrap();
471        assert::close(lognorm.skewness().unwrap(), 33_936_928.306_623_81, TOL);
472    }
473
474    #[test]
475    fn kurtosis() {
476        let lognorm = LogNormal::new(-1.2, 1.0).unwrap();
477        assert::close(lognorm.kurtosis().unwrap(), 110.936_392_176_311_53, TOL);
478    }
479
480    #[test]
481    fn cdf_standard_at_one_should_be_one_half() {
482        let lognorm1 = LogNormal::new(0.0, 1.0).unwrap();
483        assert::close(lognorm1.cdf(&1.0_f64), 0.5, TOL);
484    }
485
486    #[test]
487    fn cdf_standard_value_at_two() {
488        let lognorm = LogNormal::standard();
489        assert::close(lognorm.cdf(&2.0_f64), 0.755_891_404_214_417_3, TOL);
490    }
491
492    proptest! {
493        #[test]
494        fn quantile_agree_with_cdf(p in 0.0..1.0) {
495            prop_assume!(p > 0.0);
496            prop_assume!(p < 1.0);
497            let lognorm = LogNormal::standard();
498            let y: f64 = lognorm.quantile(p);
499            let p2 = lognorm.cdf(&y);
500            assert::close(p, p2, TOL);
501        }
502    }
503
504    #[test]
505    fn entropy() {
506        let lognorm = LogNormal::new(1.2, 3.4).unwrap();
507        assert::close(lognorm.entropy(), 3.842_713_964_826_788_5, TOL);
508    }
509
510    #[test]
511    fn entropy_standard() {
512        let lognorm = LogNormal::standard();
513        assert::close(lognorm.entropy(), 1.418_938_533_204_672_7, TOL);
514    }
515}
rv/dist/lognormal.rs

rv/dist/
lognormal.rs
