Skip to main content

numra_stats/distributions/
lognormal.rs

1//! Log-normal distribution.
2//!
3//! Author: Moussa Leblouba
4//! Date: 9 February 2026
5//! Modified: 2 May 2026
6
7use numra_core::Scalar;
8use numra_special::{erf, erfinv};
9use rand::RngCore;
10
11use super::normal::Normal;
12use super::ContinuousDistribution;
13
14/// Log-normal distribution: ln(X) ~ N(mu, sigma^2).
15#[derive(Clone, Debug)]
16pub struct LogNormal<S: Scalar> {
17    pub mu: S,
18    pub sigma: S,
19}
20
21impl<S: Scalar> LogNormal<S> {
22    pub fn new(mu: S, sigma: S) -> Self {
23        Self { mu, sigma }
24    }
25}
26
27impl<S: Scalar> ContinuousDistribution<S> for LogNormal<S> {
28    fn pdf(&self, x: S) -> S {
29        if x <= S::ZERO {
30            return S::ZERO;
31        }
32        let two = S::TWO;
33        let pi2 = S::from_f64(core::f64::consts::TAU);
34        let z = (x.ln() - self.mu) / self.sigma;
35        (S::ZERO - z * z / two).exp() / (x * self.sigma * pi2.sqrt())
36    }
37
38    fn cdf(&self, x: S) -> S {
39        if x <= S::ZERO {
40            return S::ZERO;
41        }
42        let sqrt2 = S::from_f64(core::f64::consts::SQRT_2);
43        let half = S::HALF;
44        let z = (x.ln() - self.mu) / (self.sigma * sqrt2);
45        half * (S::ONE + erf(z))
46    }
47
48    fn quantile(&self, p: S) -> S {
49        let sqrt2 = S::from_f64(core::f64::consts::SQRT_2);
50        let two = S::TWO;
51        (self.mu + self.sigma * sqrt2 * erfinv(two * p - S::ONE)).exp()
52    }
53
54    fn mean(&self) -> S {
55        let half = S::HALF;
56        (self.mu + half * self.sigma * self.sigma).exp()
57    }
58
59    fn variance(&self) -> S {
60        let s2 = self.sigma * self.sigma;
61        ((S::TWO * self.mu + s2).exp()) * (s2.exp() - S::ONE)
62    }
63
64    fn sample(&self, rng: &mut dyn RngCore) -> S {
65        let normal = Normal::new(self.mu, self.sigma);
66        normal.sample(rng).exp()
67    }
68}
69
70#[cfg(test)]
71mod tests {
72    use super::*;
73
74    #[test]
75    fn test_lognormal_pdf_positive_only() {
76        let ln = LogNormal::new(0.0_f64, 1.0);
77        assert!(ln.pdf(-1.0).abs() < 1e-14);
78        assert!(ln.pdf(0.0).abs() < 1e-14);
79        assert!(ln.pdf(1.0) > 0.0);
80    }
81
82    #[test]
83    fn test_lognormal_cdf() {
84        let ln = LogNormal::new(0.0_f64, 1.0);
85        // CDF at x=1 should be 0.5 (since ln(1) = 0 = mean)
86        assert!((ln.cdf(1.0) - 0.5).abs() < 1e-10);
87    }
88
89    #[test]
90    fn test_lognormal_quantile_roundtrip() {
91        let ln = LogNormal::new(1.0_f64, 0.5);
92        for &p in &[0.1, 0.5, 0.9] {
93            let x = ln.quantile(p);
94            let p2 = ln.cdf(x);
95            assert!((p - p2).abs() < 1e-8, "p={}, p2={}", p, p2);
96        }
97    }
98
99    #[test]
100    fn test_lognormal_mean() {
101        let ln = LogNormal::new(0.0_f64, 1.0);
102        // E[X] = exp(mu + sigma^2/2) = exp(0.5)
103        assert!((ln.mean() - 0.5_f64.exp()).abs() < 1e-12);
104    }
105}