Skip to main content

stats_claw/distributions/positive/
lognormal.rs

1//! Log-normal distribution numerics, for the
2//! [`LogNormalDistribution`].
3//!
4//! Equivalent to `scipy.stats.lognorm(s = std_log_value, scale =
5//! exp(mean_log_value))`: a variable whose logarithm is normal with mean
6//! `mean_log_value` and standard deviation `std_log_value`. The CDF and quantile
7//! delegate to a standard [`NormalDistribution`] on the log scale, reusing its
8//! erf-based machinery; sampling exponentiates a normal draw.
9//!
10//! # Examples
11//!
12//! ```
13//! use stats_claw::distributions::Quantile;
14//! use stats_claw::distributions::LogNormalDistribution;
15//!
16//! let d = LogNormalDistribution { mean_log_value: 0.0, std_log_value: 1.0, ..Default::default() };
17//! // Median is exp(mean_log_value) = exp(0) = 1.
18//! assert!((d.quantile(0.5) - 1.0).abs() < 1e-9, "median was {}", d.quantile(0.5));
19//! ```
20
21use super::super::{Cdf, Moments, Pdf, Quantile, Sample};
22use crate::distributions::LogNormalDistribution;
23use crate::distributions::NormalDistribution;
24use crate::rng::SplitMix64;
25use std::f64::consts::PI;
26
27impl LogNormalDistribution {
28    /// The unit normal on the log scale this distribution is built from.
29    fn log_normal(&self) -> NormalDistribution {
30        NormalDistribution {
31            mean: self.mean_log_value,
32            standard_deviation: self.std_log_value,
33            ..Default::default()
34        }
35    }
36}
37
38impl Pdf for LogNormalDistribution {
39    fn pdf(&self, x: f64) -> f64 {
40        if x <= 0.0 {
41            return 0.0;
42        }
43        let (mu, sigma) = (self.mean_log_value, self.std_log_value);
44        let z = (x.ln() - mu) / sigma;
45        (-0.5 * z * z).exp() / (x * sigma * (2.0 * PI).sqrt())
46    }
47}
48
49impl Cdf for LogNormalDistribution {
50    fn cdf(&self, x: f64) -> f64 {
51        if x <= 0.0 {
52            0.0
53        } else {
54            self.log_normal().cdf(x.ln())
55        }
56    }
57}
58
59impl Quantile for LogNormalDistribution {
60    fn quantile(&self, p: f64) -> f64 {
61        self.log_normal().quantile(p).exp()
62    }
63}
64
65impl Moments for LogNormalDistribution {
66    fn mean(&self) -> Option<f64> {
67        let (mu, sigma) = (self.mean_log_value, self.std_log_value);
68        Some((0.5 * sigma).mul_add(sigma, mu).exp())
69    }
70    fn variance(&self) -> Option<f64> {
71        let (mu, sigma) = (self.mean_log_value, self.std_log_value);
72        let s2 = sigma * sigma;
73        Some(s2.exp_m1() * 2.0f64.mul_add(mu, s2).exp())
74    }
75}
76
77impl Sample for LogNormalDistribution {
78    fn sample(&self, rng: &mut SplitMix64) -> f64 {
79        self.std_log_value
80            .mul_add(rng.standard_normal(), self.mean_log_value)
81            .exp()
82    }
83}
84
85#[cfg(test)]
86mod tests {
87    use super::*;
88
89    /// The median is `exp(mean_log_value)`.
90    #[test]
91    fn median_is_exp_mu() {
92        let d = LogNormalDistribution {
93            mean_log_value: 0.3,
94            std_log_value: 0.6,
95            ..Default::default()
96        };
97        assert!((d.quantile(0.5) - 0.3_f64.exp()).abs() < 1e-9);
98    }
99
100    /// The mean is `exp(μ + σ²/2)`.
101    #[test]
102    fn mean_matches_closed_form() {
103        let d = LogNormalDistribution {
104            mean_log_value: 0.0,
105            std_log_value: 1.0,
106            ..Default::default()
107        };
108        // exp(0.5) ≈ 1.648721
109        assert!(matches!(d.mean(), Some(m) if (m - 1.648_721_271).abs() < 1e-6));
110    }
111}