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