numra_stats/distributions/
lognormal.rs1use numra_core::Scalar;
8use numra_special::{erf, erfinv};
9use rand::RngCore;
10
11use super::normal::Normal;
12use super::ContinuousDistribution;
13
14#[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 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 assert!((ln.mean() - 0.5_f64.exp()).abs() < 1e-12);
104 }
105}