Skip to main content

stats_claw/distributions/positive/
lognormal.rs

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