stats_claw/distributions/positive/
lognormal.rs1use 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 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 #[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 #[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 assert!(matches!(d.mean(), Some(m) if (m - 1.648_721_271).abs() < 1e-6));
108 }
109}