statrs/distribution/
log_normal.rs

1use crate::consts;
2use crate::distribution::{Continuous, ContinuousCDF};
3use crate::function::erf;
4use crate::statistics::*;
5use std::f64;
6
7/// Implements the
8/// [Log-normal](https://en.wikipedia.org/wiki/Log-normal_distribution)
9/// distribution
10///
11/// # Examples
12///
13/// ```
14/// use statrs::distribution::{LogNormal, Continuous};
15/// use statrs::statistics::Distribution;
16/// use statrs::prec;
17///
18/// let n = LogNormal::new(0.0, 1.0).unwrap();
19/// assert_eq!(n.mean().unwrap(), (0.5f64).exp());
20/// assert!(prec::almost_eq(n.pdf(1.0), 0.3989422804014326779399, 1e-16));
21/// ```
22#[derive(Copy, Clone, PartialEq, Debug)]
23pub struct LogNormal {
24    location: f64,
25    scale: f64,
26}
27
28/// Represents the errors that can occur when creating a [`LogNormal`].
29#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
30#[non_exhaustive]
31pub enum LogNormalError {
32    /// The location is NaN.
33    LocationInvalid,
34
35    /// The scale is NaN, zero or less than zero.
36    ScaleInvalid,
37}
38
39impl std::fmt::Display for LogNormalError {
40    #[cfg_attr(coverage_nightly, coverage(off))]
41    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
42        match self {
43            LogNormalError::LocationInvalid => write!(f, "Location is NaN"),
44            LogNormalError::ScaleInvalid => write!(f, "Scale is NaN, zero or less than zero"),
45        }
46    }
47}
48
49impl std::error::Error for LogNormalError {}
50
51impl LogNormal {
52    /// Constructs a new log-normal distribution with a location of `location`
53    /// and a scale of `scale`
54    ///
55    /// # Errors
56    ///
57    /// Returns an error if `location` or `scale` are `NaN`.
58    /// Returns an error if `scale <= 0.0`
59    ///
60    /// # Examples
61    ///
62    /// ```
63    /// use statrs::distribution::LogNormal;
64    ///
65    /// let mut result = LogNormal::new(0.0, 1.0);
66    /// assert!(result.is_ok());
67    ///
68    /// result = LogNormal::new(0.0, 0.0);
69    /// assert!(result.is_err());
70    /// ```
71    pub fn new(location: f64, scale: f64) -> Result<LogNormal, LogNormalError> {
72        if location.is_nan() {
73            return Err(LogNormalError::LocationInvalid);
74        }
75
76        if scale.is_nan() || scale <= 0.0 {
77            return Err(LogNormalError::ScaleInvalid);
78        }
79
80        Ok(LogNormal { location, scale })
81    }
82}
83
84impl std::fmt::Display for LogNormal {
85    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
86        write!(f, "LogNormal({}, {}^2)", self.location, self.scale)
87    }
88}
89
90#[cfg(feature = "rand")]
91#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
92impl ::rand::distributions::Distribution<f64> for LogNormal {
93    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
94        super::normal::sample_unchecked(rng, self.location, self.scale).exp()
95    }
96}
97
98impl ContinuousCDF<f64, f64> for LogNormal {
99    /// Calculates the cumulative distribution function for the log-normal
100    /// distribution
101    /// at `x`
102    ///
103    /// # Formula
104    ///
105    /// ```text
106    /// (1 / 2) + (1 / 2) * erf((ln(x) - μ) / sqrt(2) * σ)
107    /// ```
108    ///
109    /// where `μ` is the location, `σ` is the scale, and `erf` is the
110    /// error function
111    fn cdf(&self, x: f64) -> f64 {
112        if x <= 0.0 {
113            0.0
114        } else if x.is_infinite() {
115            1.0
116        } else {
117            0.5 * erf::erfc((self.location - x.ln()) / (self.scale * f64::consts::SQRT_2))
118        }
119    }
120
121    /// Calculates the survival function for the log-normal
122    /// distribution at `x`
123    ///
124    /// # Formula
125    ///
126    /// ```text
127    /// (1 / 2) + (1 / 2) * erf(-(ln(x) - μ) / sqrt(2) * σ)
128    /// ```
129    ///
130    /// where `μ` is the location, `σ` is the scale, and `erf` is the
131    /// error function
132    ///
133    /// note that this calculates the complement due to flipping
134    /// the sign of the argument error function with respect to the cdf.
135    ///
136    /// the normal cdf Φ (and internal error function) as the following property:
137    /// ```text
138    ///  Φ(-x) + Φ(x) = 1
139    ///  Φ(-x)        = 1 - Φ(x)
140    /// ```
141    fn sf(&self, x: f64) -> f64 {
142        if x <= 0.0 {
143            1.0
144        } else if x.is_infinite() {
145            0.0
146        } else {
147            0.5 * erf::erfc((x.ln() - self.location) / (self.scale * f64::consts::SQRT_2))
148        }
149    }
150
151    /// Calculates the inverse cumulative distribution function for the
152    /// log-normal distribution at `p`
153    ///
154    /// # Panics
155    ///
156    /// If `p < 0.0` or `p > 1.0`
157    ///
158    /// # Formula
159    ///
160    /// ```text
161    /// μ - σ * sqrt(2) * erfc_inv(2p)
162    /// ```
163    ///
164    /// where `μ` is the location, `σ` is the scale and `erfc_inv` is
165    /// the inverse of the complementary error function
166    fn inverse_cdf(&self, p: f64) -> f64 {
167        if p == 0.0 {
168            0.0
169        } else if p < 1.0 {
170            (self.location - (self.scale * f64::consts::SQRT_2 * erf::erfc_inv(2.0 * p))).exp()
171        } else if p == 1.0 {
172            f64::INFINITY
173        } else {
174            panic!("p must be within [0.0, 1.0]");
175        }
176    }
177}
178
179impl Min<f64> for LogNormal {
180    /// Returns the minimum value in the domain of the log-normal
181    /// distribution representable by a double precision float
182    ///
183    /// # Formula
184    ///
185    /// ```text
186    /// 0
187    /// ```
188    fn min(&self) -> f64 {
189        0.0
190    }
191}
192
193impl Max<f64> for LogNormal {
194    /// Returns the maximum value in the domain of the log-normal
195    /// distribution representable by a double precision float
196    ///
197    /// # Formula
198    ///
199    /// ```text
200    /// f64::INFINITY
201    /// ```
202    fn max(&self) -> f64 {
203        f64::INFINITY
204    }
205}
206
207impl Distribution<f64> for LogNormal {
208    /// Returns the mean of the log-normal distribution
209    ///
210    /// # Formula
211    ///
212    /// ```text
213    /// e^(μ + σ^2 / 2)
214    /// ```
215    ///
216    /// where `μ` is the location and `σ` is the scale
217    fn mean(&self) -> Option<f64> {
218        Some((self.location + self.scale * self.scale / 2.0).exp())
219    }
220
221    /// Returns the variance of the log-normal distribution
222    ///
223    /// # Formula
224    ///
225    /// ```text
226    /// (e^(σ^2) - 1) * e^(2μ + σ^2)
227    /// ```
228    ///
229    /// where `μ` is the location and `σ` is the scale
230    fn variance(&self) -> Option<f64> {
231        let sigma2 = self.scale * self.scale;
232        Some((sigma2.exp() - 1.0) * (self.location + self.location + sigma2).exp())
233    }
234
235    /// Returns the entropy of the log-normal distribution
236    ///
237    /// # Formula
238    ///
239    /// ```text
240    /// ln(σe^(μ + 1 / 2) * sqrt(2π))
241    /// ```
242    ///
243    /// where `μ` is the location and `σ` is the scale
244    fn entropy(&self) -> Option<f64> {
245        Some(0.5 + self.scale.ln() + self.location + consts::LN_SQRT_2PI)
246    }
247
248    /// Returns the skewness of the log-normal distribution
249    ///
250    /// # Formula
251    ///
252    /// ```text
253    /// (e^(σ^2) + 2) * sqrt(e^(σ^2) - 1)
254    /// ```
255    ///
256    /// where `μ` is the location and `σ` is the scale
257    fn skewness(&self) -> Option<f64> {
258        let expsigma2 = (self.scale * self.scale).exp();
259        Some((expsigma2 + 2.0) * (expsigma2 - 1.0).sqrt())
260    }
261}
262
263impl Median<f64> for LogNormal {
264    /// Returns the median of the log-normal distribution
265    ///
266    /// # Formula
267    ///
268    /// ```text
269    /// e^μ
270    /// ```
271    ///
272    /// where `μ` is the location
273    fn median(&self) -> f64 {
274        self.location.exp()
275    }
276}
277
278impl Mode<Option<f64>> for LogNormal {
279    /// Returns the mode of the log-normal distribution
280    ///
281    /// # Formula
282    ///
283    /// ```text
284    /// e^(μ - σ^2)
285    /// ```
286    ///
287    /// where `μ` is the location and `σ` is the scale
288    fn mode(&self) -> Option<f64> {
289        Some((self.location - self.scale * self.scale).exp())
290    }
291}
292
293impl Continuous<f64, f64> for LogNormal {
294    /// Calculates the probability density function for the log-normal
295    /// distribution at `x`
296    ///
297    /// # Formula
298    ///
299    /// ```text
300    /// (1 / xσ * sqrt(2π)) * e^(-((ln(x) - μ)^2) / 2σ^2)
301    /// ```
302    ///
303    /// where `μ` is the location and `σ` is the scale
304    fn pdf(&self, x: f64) -> f64 {
305        if x <= 0.0 || x.is_infinite() {
306            0.0
307        } else {
308            let d = (x.ln() - self.location) / self.scale;
309            (-0.5 * d * d).exp() / (x * consts::SQRT_2PI * self.scale)
310        }
311    }
312
313    /// Calculates the log probability density function for the log-normal
314    /// distribution at `x`
315    ///
316    /// # Formula
317    ///
318    /// ```text
319    /// ln((1 / xσ * sqrt(2π)) * e^(-((ln(x) - μ)^2) / 2σ^2))
320    /// ```
321    ///
322    /// where `μ` is the location and `σ` is the scale
323    fn ln_pdf(&self, x: f64) -> f64 {
324        if x <= 0.0 || x.is_infinite() {
325            f64::NEG_INFINITY
326        } else {
327            let d = (x.ln() - self.location) / self.scale;
328            (-0.5 * d * d) - consts::LN_SQRT_2PI - (x * self.scale).ln()
329        }
330    }
331}
332
333#[rustfmt::skip]
334#[cfg(test)]
335mod tests {
336    use super::*;
337    use crate::distribution::internal::*;
338    use crate::testing_boiler;
339
340    testing_boiler!(location: f64, scale: f64; LogNormal; LogNormalError);
341
342    #[test]
343    fn test_create() {
344        create_ok(10.0, 0.1);
345        create_ok(-5.0, 1.0);
346        create_ok(0.0, 10.0);
347        create_ok(10.0, 100.0);
348        create_ok(-5.0, f64::INFINITY);
349    }
350
351    #[test]
352    fn test_bad_create() {
353        test_create_err(f64::NAN, 1.0, LogNormalError::LocationInvalid);
354        test_create_err(1.0, f64::NAN, LogNormalError::ScaleInvalid);
355        create_err(0.0, 0.0);
356        create_err(f64::NAN, f64::NAN);
357        create_err(1.0, -1.0);
358    }
359
360    #[test]
361    fn test_mean() {
362        let mean = |x: LogNormal| x.mean().unwrap();
363        test_exact(-1.0, 0.1, 0.369723444544058982601, mean);
364        test_exact(-1.0, 1.5, 1.133148453066826316829, mean);
365        test_exact(-1.0, 2.5, 8.372897488127264663205, mean);
366        test_exact(-1.0, 5.5, 1362729.18425285481771, mean);
367        test_exact(-0.1, 0.1, 0.9093729344682314204933, mean);
368        test_exact(-0.1, 1.5, 2.787095460565850768514, mean);
369        test_exact(-0.1, 2.5, 20.59400471119602917533, mean);
370        test_absolute(-0.1, 5.5, 3351772.941252693807591, 1e-9, mean);
371        test_exact(0.1, 0.1, 1.110710610355705232259, mean);
372        test_exact(0.1, 1.5, 3.40416608279081898632, mean);
373        test_absolute(0.1, 2.5, 25.15357415581836182776, 1e-14, mean);
374        test_absolute(0.1, 5.5, 4093864.715172665106863, 1e-8, mean);
375        test_absolute(1.5, 0.1, 4.50415363028848413209, 1e-15, mean);
376        test_exact(1.5, 1.5, 13.80457418606709491926, mean);
377        test_exact(1.5, 2.5, 102.0027730826996844534, mean);
378        test_exact(1.5, 5.5, 16601440.05723477471392, mean);
379        test_absolute(2.5, 0.1, 12.24355896580102707724, 1e-14, mean);
380        test_absolute(2.5, 1.5, 37.52472315960099891407, 1e-11, mean);
381        test_exact(2.5, 2.5, 277.2722845231339804081, mean);
382        test_exact(2.5, 5.5, 45127392.83383337999291, mean);
383        test_absolute(5.5, 0.1, 245.9184556788219446833, 1e-13, mean);
384        test_exact(5.5, 1.5, 753.7042125545612656606, mean);
385        test_exact(5.5, 2.5, 5569.162708566004074422, mean);
386        test_exact(5.5, 5.5, 906407915.0111549133446, mean);
387    }
388
389    #[test]
390    fn test_variance() {
391        let variance = |x: LogNormal| x.variance().unwrap();
392        test_absolute(-1.0, 0.1, 0.001373811865368952608715, 1e-16, variance);
393        test_exact(-1.0, 1.5, 10.898468544015731954, variance);
394        test_exact(-1.0, 2.5, 36245.39726189994988081, variance);
395        test_absolute(-1.0, 5.5, 2.5481629178024539E+25, 1e10, variance);
396        test_absolute(-0.1, 0.1, 0.008311077467909703803238, 1e-16, variance);
397        test_exact(-0.1, 1.5, 65.93189259328902509552, variance);
398        test_absolute(-0.1, 2.5, 219271.8756420929704707, 1e-10, variance);
399        test_absolute(-0.1, 5.5, 1.541548733459471E+26, 1e12, variance);
400        test_absolute(0.1, 0.1, 0.01239867063063756838894, 1e-15, variance);
401        test_absolute(0.1, 1.5, 98.35882573290010981464, 1e-13, variance);
402        test_absolute(0.1, 2.5, 327115.1995809995715014, 1e-10, variance);
403        test_absolute(0.1, 5.5, 2.299720473192458E+26, 1e12, variance);
404        test_absolute(1.5, 0.1, 0.2038917589520099120699, 1e-14, variance);
405        test_absolute(1.5, 1.5, 1617.476145997433210727, 1e-12, variance);
406        test_absolute(1.5, 2.5, 5379293.910566451644527, 1e-9, variance);
407        test_absolute(1.5, 5.5, 3.7818090853910142E+27, 1e12, variance);
408        test_absolute(2.5, 0.1, 1.506567645006046841936, 1e-13, variance);
409        test_absolute(2.5, 1.5, 11951.62198145717670088, 1e-11, variance);
410        test_exact(2.5, 2.5, 39747904.47781154725843, variance);
411        test_absolute(2.5, 5.5, 2.7943999487399818E+28, 1e13, variance);
412        test_absolute(5.5, 0.1, 607.7927673399807484235, 1e-11, variance);
413        test_exact(5.5, 1.5, 4821628.436260521100027, variance);
414        test_exact(5.5, 2.5, 16035449147.34799637823, variance);
415        test_exact(5.5, 5.5, 1.127341399856331737823E+31, variance);
416    }
417
418    #[test]
419    fn test_entropy() {
420        let entropy = |x: LogNormal| x.entropy().unwrap();
421        test_exact(-1.0, 0.1, -1.8836465597893728867265104870209210873020761202386, entropy);
422        test_exact(-1.0, 1.5, 0.82440364131283712375834285186996677643338789710028, entropy);
423        test_exact(-1.0, 2.5, 1.335229265078827806963856948173628711311498693546, entropy);
424        test_exact(-1.0, 5.5, 2.1236866254430979764250411929125703716076041932149, entropy);
425        test_absolute(-0.1, 0.1, -0.9836465597893728922776256101467037894202344606927, 1e-15, entropy);
426        test_exact(-0.1, 1.5, 1.7244036413128371182072277287441840743152295566462, entropy);
427        test_exact(-0.1, 2.5, 2.2352292650788278014127418250478460091933403530919, entropy);
428        test_exact(-0.1, 5.5, 3.0236866254430979708739260697867876694894458527608, entropy);
429        test_absolute(0.1, 0.1, -0.7836465597893728811753953638951383851839177797845, 1e-15, entropy);
430        test_absolute(0.1, 1.5, 1.9244036413128371293094579749957494785515462375544, 1e-15, entropy);
431        test_exact(0.1, 2.5, 2.4352292650788278125149720712994114134296570340001, entropy);
432        test_exact(0.1, 5.5, 3.223686625443097981976156316038353073725762533669, entropy);
433        test_absolute(1.5, 0.1, 0.6163534402106271132734895129790789126979238797614, 1e-15, entropy);
434        test_exact(1.5, 1.5, 3.3244036413128371237583428518699667764333878971003, entropy);
435        test_exact(1.5, 2.5, 3.835229265078827806963856948173628711311498693546, entropy);
436        test_exact(1.5, 5.5, 4.6236866254430979764250411929125703716076041932149, entropy);
437        test_exact(2.5, 0.1, 1.6163534402106271132734895129790789126979238797614, entropy);
438        test_absolute(2.5, 1.5, 4.3244036413128371237583428518699667764333878971003, 1e-15, entropy);
439        test_exact(2.5, 2.5, 4.835229265078827806963856948173628711311498693546, entropy);
440        test_exact(2.5, 5.5, 5.6236866254430979764250411929125703716076041932149, entropy);
441        test_exact(5.5, 0.1, 4.6163534402106271132734895129790789126979238797614, entropy);
442        test_absolute(5.5, 1.5, 7.3244036413128371237583428518699667764333878971003, 1e-15, entropy);
443        test_exact(5.5, 2.5, 7.835229265078827806963856948173628711311498693546, entropy);
444        test_exact(5.5, 5.5, 8.6236866254430979764250411929125703716076041932149, entropy);
445    }
446
447    #[test]
448    fn test_skewness() {
449        let skewness = |x: LogNormal| x.skewness().unwrap();
450        test_absolute(-1.0, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness);
451        test_exact(-1.0, 1.5, 33.46804679732172529147579024311650645764144530123, skewness);
452        test_absolute(-1.0, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness);
453        test_absolute(-1.0, 5.5, 50829064464591483629.132631635472412625371367420496, 1e4, skewness);
454        test_absolute(-0.1, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness);
455        test_exact(-0.1, 1.5, 33.46804679732172529147579024311650645764144530123, skewness);
456        test_absolute(-0.1, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness);
457        test_absolute(-0.1, 5.5, 50829064464591483629.132631635472412625371367420496, 1e4, skewness);
458        test_absolute(0.1, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness);
459        test_exact(0.1, 1.5, 33.46804679732172529147579024311650645764144530123, skewness);
460        test_absolute(0.1, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness);
461        test_absolute(0.1, 5.5, 50829064464591483629.132631635472412625371367420496, 1e4, skewness);
462        test_absolute(1.5, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness);
463        test_exact(1.5, 1.5, 33.46804679732172529147579024311650645764144530123, skewness);
464        test_absolute(1.5, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness);
465        test_absolute(1.5, 5.5, 50829064464591483629.132631635472412625371367420496, 1e4, skewness);
466        test_absolute(2.5, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness);
467        test_exact(2.5, 1.5, 33.46804679732172529147579024311650645764144530123, skewness);
468        test_absolute(2.5, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness);
469        test_absolute(2.5, 5.5, 50829064464591483629.132631635472412625371367420496, 1e4, skewness);
470        test_absolute(5.5, 0.1, 0.30175909933883402945387113824982918009810212213629, 1e-14, skewness);
471        test_exact(5.5, 1.5, 33.46804679732172529147579024311650645764144530123, skewness);
472        test_absolute(5.5, 2.5, 11824.007933610287521341659465200553739278936344799, 1e-11, skewness);
473        test_absolute(5.5, 5.5, 50829064464591483629.132631635472412625371367420496, 1e4, skewness);
474    }
475
476    #[test]
477    fn test_mode() {
478        let mode = |x: LogNormal| x.mode().unwrap();
479        test_exact(-1.0, 0.1, 0.36421897957152331652213191863106773137983085909534, mode);
480        test_exact(-1.0, 1.5, 0.03877420783172200988689983526759614326014406193602, mode);
481        test_exact(-1.0, 2.5, 0.0007101743888425490635846003705775444086763023873619, mode);
482        test_exact(-1.0, 5.5, 0.000000000000026810038677818032221548731163905979029274677187036, mode);
483        test_exact(-0.1, 0.1, 0.89583413529652823774737070060865897390995185639633, mode);
484        test_exact(-0.1, 1.5, 0.095369162215549610417813418326627245539514227574881, mode);
485        test_exact(-0.1, 2.5, 0.0017467471362611196181003627521060283221112106850165, mode);
486        test_exact(-0.1, 5.5, 0.00000000000006594205454219929159167575814655534255162059017114, mode);
487        test_exact(0.1, 0.1, 1.0941742837052103542285651753780976842292770841345, mode);
488        test_exact(0.1, 1.5, 0.11648415777349696821514223131929465848700730137808, mode);
489        test_exact(0.1, 2.5, 0.0021334817700377079925027678518795817076296484352472, mode);
490        test_exact(0.1, 5.5, 0.000000000000080541807296590798973741710866097756565304960216803, mode);
491        test_exact(1.5, 0.1, 4.4370955190036645692996309927420381428715912422597, mode);
492        test_exact(1.5, 1.5, 0.47236655274101470713804655094326791297020357913648, mode);
493        test_exact(1.5, 2.5, 0.008651695203120634177071503957250390848166331197708, mode);
494        test_exact(1.5, 5.5, 0.00000000000032661313427874471360158184468030186601222739665225, mode);
495        test_exact(2.5, 0.1, 12.061276120444720299113038763305617245808510584994, mode);
496        test_exact(2.5, 1.5, 1.2840254166877414840734205680624364583362808652815, mode);
497        test_exact(2.5, 2.5, 0.023517745856009108236151185100432939470067655273072, mode);
498        test_exact(2.5, 5.5, 0.00000000000088782654784596584473099190326928541185172970391855, mode);
499        test_exact(5.5, 0.1, 242.2572068579541371904816252345031593584721473492, mode);
500        test_exact(5.5, 1.5, 25.790339917193062089080107669377221876655268848954, mode);
501        test_exact(5.5, 2.5, 0.47236655274101470713804655094326791297020357913648, mode);
502        test_exact(5.5, 5.5, 0.000000000017832472908146389493511850431527026413424899198327, mode);
503    }
504
505    #[test]
506    fn test_median() {
507        let median = |x: LogNormal| x.median();
508        test_exact(-1.0, 0.1, 0.36787944117144232159552377016146086744581113103177, median);
509        test_exact(-1.0, 1.5, 0.36787944117144232159552377016146086744581113103177, median);
510        test_exact(-1.0, 2.5, 0.36787944117144232159552377016146086744581113103177, median);
511        test_exact(-1.0, 5.5, 0.36787944117144232159552377016146086744581113103177, median);
512        test_exact(-0.1, 0.1, 0.90483741803595956814139238421693559530906465375738, median);
513        test_exact(-0.1, 1.5, 0.90483741803595956814139238421693559530906465375738, median);
514        test_exact(-0.1, 2.5, 0.90483741803595956814139238421693559530906465375738, median);
515        test_exact(-0.1, 5.5, 0.90483741803595956814139238421693559530906465375738, median);
516        test_exact(0.1, 0.1, 1.1051709180756476309466388234587796577416634163742, median);
517        test_exact(0.1, 1.5, 1.1051709180756476309466388234587796577416634163742, median);
518        test_exact(0.1, 2.5, 1.1051709180756476309466388234587796577416634163742, median);
519        test_exact(0.1, 5.5, 1.1051709180756476309466388234587796577416634163742, median);
520        test_exact(1.5, 0.1, 4.4816890703380648226020554601192758190057498683697, median);
521        test_exact(1.5, 1.5, 4.4816890703380648226020554601192758190057498683697, median);
522        test_exact(1.5, 2.5, 4.4816890703380648226020554601192758190057498683697, median);
523        test_exact(1.5, 5.5, 4.4816890703380648226020554601192758190057498683697, median);
524        test_exact(2.5, 0.1, 12.182493960703473438070175951167966183182767790063, median);
525        test_exact(2.5, 1.5, 12.182493960703473438070175951167966183182767790063, median);
526        test_exact(2.5, 2.5, 12.182493960703473438070175951167966183182767790063, median);
527        test_exact(2.5, 5.5, 12.182493960703473438070175951167966183182767790063, median);
528        test_exact(5.5, 0.1, 244.6919322642203879151889495118393501842287101075, median);
529        test_exact(5.5, 1.5, 244.6919322642203879151889495118393501842287101075, median);
530        test_exact(5.5, 2.5, 244.6919322642203879151889495118393501842287101075, median);
531        test_exact(5.5, 5.5, 244.6919322642203879151889495118393501842287101075, median);
532    }
533
534    #[test]
535    fn test_min_max() {
536        let min = |x: LogNormal| x.min();
537        let max = |x: LogNormal| x.max();
538        test_exact(0.0, 0.1, 0.0, min);
539        test_exact(-3.0, 10.0, 0.0, min);
540        test_exact(0.0, 0.1, f64::INFINITY, max);
541        test_exact(-3.0, 10.0, f64::INFINITY, max);
542    }
543
544    #[test]
545    fn test_pdf() {
546        let pdf = |arg: f64| move |x: LogNormal| x.pdf(arg);
547        test_absolute(-0.1, 0.1, 1.7968349035073582236359415565799753846986440127816e-104, 1e-118, pdf(0.1));
548        test_absolute(-0.1, 0.1, 0.00000018288923328441197822391757965928083462391836798722, 1e-21, pdf(0.5));
549        test_exact(-0.1, 0.1, 2.3363114904470413709866234247494393485647978367885, pdf(0.8));
550        test_absolute(-0.1, 1.5, 0.90492497850024368541682348133921492204585092983646, 1e-15, pdf(0.1));
551        test_absolute(-0.1, 1.5, 0.49191985207660942803818797602364034466489243416574, 1e-16, pdf(0.5));
552        test_exact(-0.1, 1.5, 0.33133347214343229148978298237579567194870525187207, pdf(0.8));
553        test_exact(-0.1, 2.5, 1.0824698632626565182080576574958317806389057196768, pdf(0.1));
554        test_absolute(-0.1, 2.5, 0.31029619474753883558901295436486123689563749784867, 1e-16, pdf(0.5));
555        test_absolute(-0.1, 2.5, 0.19922929916156673799861939824205622734205083805245, 1e-16, pdf(0.8));
556
557// Test removed because it was causing compiler issues (see issue 31407 for rust)
558// test_absolute(1.5, 0.1, 4.1070141770545881694056265342787422035256248474059e-313, 1e-322, pdf(0.1));
559//
560
561        test_absolute(1.5, 0.1, 2.8602688726477103843476657332784045661507239533567e-104, 1e-116, pdf(0.5));
562        test_exact(1.5, 0.1, 1.6670425710002183246335601541889400558525870482613e-64, pdf(0.8));
563        test_absolute(1.5, 1.5, 0.10698412103361841220076392503406214751353235895732, 1e-16, pdf(0.1));
564        test_absolute(1.5, 1.5, 0.18266125308224685664142384493330155315630876975024, 1e-16, pdf(0.5));
565        test_absolute(1.5, 1.5, 0.17185785323404088913982425377565512294017306418953, 1e-16, pdf(0.8));
566        test_absolute(1.5, 2.5, 0.50186885259059181992025035649158160252576845315332, 1e-15, pdf(0.1));
567        test_absolute(1.5, 2.5, 0.21721369314437986034957451699565540205404697589349, 1e-16, pdf(0.5));
568        test_exact(1.5, 2.5, 0.15729636000661278918949298391170443742675565300598, pdf(0.8));
569        test_exact(2.5, 0.1, 5.6836826548848916385760779034504046896805825555997e-500, pdf(0.1));
570        test_absolute(2.5, 0.1, 3.1225608678589488061206338085285607881363155340377e-221, 1e-233, pdf(0.5));
571        test_absolute(2.5, 0.1, 4.6994713794671660918554320071312374073172560048297e-161, 1e-173, pdf(0.8));
572        test_absolute(2.5, 1.5, 0.015806486291412916772431170442330946677601577502353, 1e-16, pdf(0.1));
573        test_absolute(2.5, 1.5, 0.055184331257528847223852028950484131834529030116388, 1e-16, pdf(0.5));
574        test_exact(2.5, 1.5, 0.063982134749859504449658286955049840393511776984362, pdf(0.8));
575        test_absolute(2.5, 2.5, 0.25212505662402617595900822552548977822542300480086, 1e-15, pdf(0.1));
576        test_absolute(2.5, 2.5, 0.14117186955911792460646517002386088579088567275401, 1e-16, pdf(0.5));
577        test_absolute(2.5, 2.5, 0.11021452580363707866161369621432656293405065561317, 1e-16, pdf(0.8));
578    }
579
580    #[test]
581    fn test_neg_pdf() {
582        let pdf = |arg: f64| move |x: LogNormal| x.pdf(arg);
583        test_exact(0.0, 1.0, 0.0, pdf(0.0));
584    }
585
586    #[test]
587    fn test_ln_pdf() {
588        let ln_pdf = |arg: f64| move |x: LogNormal| x.ln_pdf(arg);
589        test_exact(-0.1, 0.1, -238.88282294119596467794686179588610665317241097599, ln_pdf(0.1));
590        test_absolute(-0.1, 0.1, -15.514385149961296196003163062199569075052113039686, 1e-14, ln_pdf(0.5));
591        test_exact(-0.1, 0.1, 0.84857339958981283964373051826407417105725729082041, ln_pdf(0.8));
592        test_absolute(-0.1, 1.5, -0.099903235403144611051953094864849327288457482212211, 1e-15, ln_pdf(0.1));
593        test_absolute(-0.1, 1.5, -0.70943947804316122682964396008813828577195771418027, 1e-15, ln_pdf(0.5));
594        test_absolute(-0.1, 1.5, -1.1046299420497998262946038709903250420774183529995, 1e-15, ln_pdf(0.8));
595        test_absolute(-0.1, 2.5, 0.07924534056485078867266307735371665927517517183681, 1e-16, ln_pdf(0.1));
596        test_exact(-0.1, 2.5, -1.1702279707433794860424967893989374511050637417043, ln_pdf(0.5));
597        test_exact(-0.1, 2.5, -1.6132988605030400828957768752511536087538109996183, ln_pdf(0.8));
598        test_exact(1.5, 0.1, -719.29643782024317312262673764204041218720576249741, ln_pdf(0.1));
599        test_absolute(1.5, 0.1, -238.41793403955250272430898754048547661932857086122, 1e-13, ln_pdf(0.5));
600        test_exact(1.5, 0.1, -146.85439481068371057247137024006716189469284256628, ln_pdf(0.8));
601        test_absolute(1.5, 1.5, -2.2350748570877992856465076624973458117562108140674, 1e-15, ln_pdf(0.1));
602        test_absolute(1.5, 1.5, -1.7001219175524556705452882616787223585705662860012, 1e-15, ln_pdf(0.5));
603        test_absolute(1.5, 1.5, -1.7610875785399045023354101841009649273236721172008, 1e-15, ln_pdf(0.8));
604        test_absolute(1.5, 2.5, -0.68941644324162489418137656699398207513321602763104, 1e-15, ln_pdf(0.1));
605        test_exact(1.5, 2.5, -1.5268736489667254857801287379715477173125628275598, ln_pdf(0.5));
606        test_exact(1.5, 2.5, -1.8496236096394777662704671479709839674424623547308, ln_pdf(0.8));
607        test_absolute(2.5, 0.1, -1149.5549471196476523788026360929146688367845019398, 1e-12, ln_pdf(0.1));
608        test_absolute(2.5, 0.1, -507.73265209554698134113704985174959301922196605736, 1e-12, ln_pdf(0.5));
609        test_absolute(2.5, 0.1, -369.16874994210463740474549611573497379941224077335, 1e-13, ln_pdf(0.8));
610        test_absolute(2.5, 1.5, -4.1473348984184862316495477617980296904955324113457, 1e-15, ln_pdf(0.1));
611        test_absolute(2.5, 1.5, -2.8970762200235424747307247601045786110485663457169, 1e-15, ln_pdf(0.5));
612        test_exact(2.5, 1.5, -2.7491513791239977024488074547907467152956602019989, ln_pdf(0.8));
613        test_absolute(2.5, 2.5, -1.3778300581206721947424710027422282714793718026513, 1e-15, ln_pdf(0.1));
614        test_exact(2.5, 2.5, -1.9577771978563167352868858774048559682046428490575, ln_pdf(0.5));
615        test_exact(2.5, 2.5, -2.2053265778497513183112901654193054111123780652581, ln_pdf(0.8));
616    }
617
618    #[test]
619    fn test_neg_ln_pdf() {
620        let ln_pdf = |arg: f64| move |x: LogNormal| x.ln_pdf(arg);
621        test_exact(0.0, 1.0, f64::NEG_INFINITY, ln_pdf(0.0));
622    }
623
624    #[test]
625    fn test_cdf() {
626        cdf_tests(false);
627    }
628
629    #[test]
630    fn test_inverse_cdf() {
631        cdf_tests(true)
632    }
633
634    // we can reuse the (input, output) pairs from the CDF unit test
635    // and verify that passing an 'output' to .inverse_cdf gives 'input',
636    // except in cases where output would be 0.0 (the inverse_cdf is defined to
637    // always give 0.0 in this case).
638    fn cdf_tests(inverse: bool) {
639        let f = |arg: f64| move |x: LogNormal| if inverse {
640            x.inverse_cdf(arg)
641        } else {
642            x.cdf(arg)
643        };
644
645        // given some cdf_input and cdf_output, returns a tuple (input, output) where
646        // input is what we will provide to cdf/inverse_cdf, and output is expected return
647        // value
648        let arrange_input_output = |cdf_input: f64, cdf_output: f64| {
649            if inverse {
650                (cdf_output, cdf_input)
651            } else {
652                (cdf_input, cdf_output)
653            }
654        };
655
656        // calls test_almost after re-arranging the input/output arguments and calling f with input
657        let almost = |mean: f64, std_dev: f64, cdf_input: f64, cdf_output: f64, acc: f64| {
658            let (input, output) = arrange_input_output(cdf_input, cdf_output);
659            test_absolute(mean, std_dev, output, acc, f(input));
660        };
661
662        // calls test_case after re-arranging the input/output arguments and calling f with input
663        let case = |mean: f64, std_dev: f64, cdf_input: f64, cdf_output: f64| {
664            let (input, output) = arrange_input_output(cdf_input, cdf_output);
665            test_exact(mean, std_dev, output, f(input));
666        };
667
668        // we skip cases where the CDF outputs 0.0 when testing the inverse CDF because
669        // there are multiple inputs to the CDF which give an answer of 0.0, therefore testing whether
670        // inputting 0.0 to the inverse cdf will give the same answer is not a valid test
671        // the inverse cdf for log-normal is defined to give answer 0.0 for input 0.0
672        if inverse {
673            case(-0.1, 0.1, 0.0, 0.0);
674        }
675
676        if !inverse {
677            almost(-0.1, 0.1, 0.1, 0.0, 1e-107);
678        }
679        almost(-0.1, 0.1, 0.5, 0.0000000015011556178148777579869633555518882664666520593658, 1e-16);
680        almost(-0.1, 0.1, 0.8, 0.10908001076375810900224507908874442583171381706127, 1e-11);
681        almost(-0.1, 1.5, 0.1, 0.070999149762464508991968731574953594549291668468349, 1e-11);
682        case(-0.1, 1.5, 0.5, 0.34626224992888089297789445771047690175505847991946);
683        case(-0.1, 1.5, 0.8, 0.46728530589487698517090261668589508746353129242404);
684        almost(-0.1, 2.5, 0.1, 0.18914969879695093477606645992572208111152994999076, 1e-10);
685        case(-0.1, 2.5, 0.5, 0.40622798321378106125020505907901206714868922279347);
686        case(-0.1, 2.5, 0.8, 0.48035707589956665425068652807400957345208517749893);
687
688        // input to inverse would be 0.0
689        if !inverse {
690            almost(1.5, 0.1, 0.1, 0.0, 1e-315);
691            almost(1.5, 0.1, 0.5, 0.0, 1e-106);
692            almost(1.5, 0.1, 0.8, 0.0, 1e-66);
693        }
694
695        almost(1.5, 1.5, 0.1, 0.005621455876973168709588070988239748831823850202953, 1e-12);
696        almost(1.5, 1.5, 0.8, 0.12532699044614938400496547188720940854423187977236, 1e-11);
697        almost(1.5, 2.5, 0.1, 0.064125647996943514411570834861724406903677144126117, 1e-11);
698        almost(1.5, 2.5, 0.5, 0.19017302281590810871719754032332631806011441356498, 1e-10);
699        almost(1.5, 2.5, 0.8, 0.24533064397555500690927047163085419096928289095201, 1e-16);
700
701        // input to inverse would be 0.0
702        if !inverse {
703            case(2.5, 0.1, 0.1, 0.0);
704            almost(2.5, 0.1, 0.5, 0.0, 1e-223);
705            almost(2.5, 0.1, 0.8, 0.0, 1e-162);
706        }
707
708        almost(2.5, 1.5, 0.1, 0.00068304052220788502001572635016579586444611070077399, 1e-13);
709        almost(2.5, 1.5, 0.5, 0.016636862816580533038130583128179878924863968664206, 1e-12);
710        almost(2.5, 1.5, 0.8, 0.034729001282904174941366974418836262996834852343018, 1e-11);
711        almost(2.5, 2.5, 0.1, 0.027363708266690978870139978537188410215717307180775, 1e-11);
712        almost(2.5, 2.5, 0.5, 0.10075543423327634536450625420610429181921642201567, 1e-11);
713        almost(2.5, 2.5, 0.8, 0.13802019192453118732001307556787218421918336849121, 1e-11);
714    }
715
716    #[test]
717    fn test_sf() {
718        let sf = |arg: f64| move |x: LogNormal| x.sf(arg);
719
720        // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(-0.1, 0.1), 0.1]
721        test_absolute(-0.1, 0.1, 1.0, 1e-107, sf(0.1));
722
723        // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(-0.1, 0.1), 0.8]
724        test_absolute(-0.1, 0.1, 0.890919989231123, 1e-14, sf(0.8));
725
726        // Wolfram Alpha:: SurvivalFunction[LogNormalDistribution[1.5, 1], 0.8]
727        test_absolute(1.5, 1.0, 0.957568715612642, 1e-14, sf(0.8));
728
729        // Wolfram Alpha:: SurvivalFunction[ LogNormalDistribution(2.5, 1.5), 0.1]
730        test_absolute(2.5, 1.5, 0.9993169594777358, 1e-14, sf(0.1));
731    }
732
733    #[test]
734    fn test_neg_cdf() {
735        let cdf = |arg: f64| move |x: LogNormal| x.cdf(arg);
736        test_exact(0.0, 1.0, 0.0, cdf(0.0));
737    }
738
739
740    #[test]
741    fn test_neg_sf() {
742        let sf = |arg: f64| move |x: LogNormal| x.sf(arg);
743        test_exact(0.0, 1.0, 1.0, sf(0.0));
744    }
745
746    #[test]
747    fn test_continuous() {
748        test::check_continuous_distribution(&create_ok(0.0, 0.25), 0.0, 10.0);
749        test::check_continuous_distribution(&create_ok(0.0, 0.5), 0.0, 10.0);
750    }
751}
statrs/distribution/log_normal.rs

statrs/distribution/
log_normal.rs
