statrs/distribution/
laplace.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::statistics::{Distribution, Max, Median, Min, Mode};
3use std::f64;
4
5/// Implements the [Laplace](https://en.wikipedia.org/wiki/Laplace_distribution)
6/// distribution.
7///
8/// # Examples
9///
10/// ```
11/// use statrs::distribution::{Laplace, Continuous};
12/// use statrs::statistics::Mode;
13///
14/// let n = Laplace::new(0.0, 1.0).unwrap();
15/// assert_eq!(n.mode().unwrap(), 0.0);
16/// assert_eq!(n.pdf(1.0), 0.18393972058572117);
17/// ```
18#[derive(Copy, Clone, PartialEq, Debug)]
19pub struct Laplace {
20    location: f64,
21    scale: f64,
22}
23
24/// Represents the errors that can occur when creating a [`Laplace`].
25#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
26#[non_exhaustive]
27pub enum LaplaceError {
28    /// The location is NaN.
29    LocationInvalid,
30
31    /// The scale is NaN, zero or less than zero.
32    ScaleInvalid,
33}
34
35impl std::fmt::Display for LaplaceError {
36    #[cfg_attr(coverage_nightly, coverage(off))]
37    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
38        match self {
39            LaplaceError::LocationInvalid => write!(f, "Location is NaN"),
40            LaplaceError::ScaleInvalid => write!(f, "Scale is NaN, zero or less than zero"),
41        }
42    }
43}
44
45impl std::error::Error for LaplaceError {}
46
47impl Laplace {
48    /// Constructs a new laplace distribution with the given
49    /// location and scale.
50    ///
51    /// # Errors
52    ///
53    /// Returns an error if location or scale are `NaN` or `scale <= 0.0`
54    ///
55    /// # Examples
56    ///
57    /// ```
58    /// use statrs::distribution::Laplace;
59    ///
60    /// let mut result = Laplace::new(0.0, 1.0);
61    /// assert!(result.is_ok());
62    ///
63    /// result = Laplace::new(0.0, -1.0);
64    /// assert!(result.is_err());
65    /// ```
66    pub fn new(location: f64, scale: f64) -> Result<Laplace, LaplaceError> {
67        if location.is_nan() {
68            return Err(LaplaceError::LocationInvalid);
69        }
70
71        if scale.is_nan() || scale <= 0.0 {
72            return Err(LaplaceError::ScaleInvalid);
73        }
74
75        Ok(Laplace { location, scale })
76    }
77
78    /// Returns the location of the laplace distribution
79    ///
80    /// # Examples
81    ///
82    /// ```
83    /// use statrs::distribution::Laplace;
84    ///
85    /// let n = Laplace::new(0.0, 1.0).unwrap();
86    /// assert_eq!(n.location(), 0.0);
87    /// ```
88    pub fn location(&self) -> f64 {
89        self.location
90    }
91
92    /// Returns the scale of the laplace distribution
93    ///
94    /// # Examples
95    ///
96    /// ```
97    /// use statrs::distribution::Laplace;
98    ///
99    /// let n = Laplace::new(0.0, 1.0).unwrap();
100    /// assert_eq!(n.scale(), 1.0);
101    /// ```
102    pub fn scale(&self) -> f64 {
103        self.scale
104    }
105}
106
107impl std::fmt::Display for Laplace {
108    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
109        write!(f, "Laplace({}, {})", self.location, self.scale)
110    }
111}
112
113#[cfg(feature = "rand")]
114#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
115impl ::rand::distributions::Distribution<f64> for Laplace {
116    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
117        let x: f64 = rng.gen_range(-0.5..0.5);
118        self.location - self.scale * x.signum() * (1. - 2. * x.abs()).ln()
119    }
120}
121
122impl ContinuousCDF<f64, f64> for Laplace {
123    /// Calculates the cumulative distribution function for the
124    /// laplace distribution at `x`
125    ///
126    /// # Formula
127    ///
128    /// ```text
129    /// (1 / 2) * (1 + signum(x - μ)) - signum(x - μ) * exp(-|x - μ| / b)
130    /// ```
131    ///
132    /// where `μ` is the location, `b` is the scale
133    fn cdf(&self, x: f64) -> f64 {
134        let y = (-(x - self.location).abs() / self.scale).exp() / 2.;
135        if x >= self.location {
136            1. - y
137        } else {
138            y
139        }
140    }
141
142    /// Calculates the survival function for the
143    /// laplace distribution at `x`
144    ///
145    /// # Formula
146    ///
147    /// ```text
148    /// 1 - [(1 / 2) * (1 + signum(x - μ)) - signum(x - μ) * exp(-|x - μ| / b)]
149    /// ```
150    ///
151    /// where `μ` is the location, `b` is the scale
152    fn sf(&self, x: f64) -> f64 {
153        let y = (-(x - self.location).abs() / self.scale).exp() / 2.;
154        if x >= self.location {
155            y
156        } else {
157            1. - y
158        }
159    }
160
161    /// Calculates the inverse cumulative distribution function for the
162    /// laplace distribution at `p`
163    ///
164    /// # Formula
165    ///
166    /// if p <= 1/2
167    /// ```text
168    /// μ + b * ln(2p)
169    /// ```
170    /// if p >= 1/2
171    /// ```text
172    /// μ - b * ln(2 - 2p)
173    /// ```
174    ///
175    /// where `μ` is the location, `b` is the scale
176    fn inverse_cdf(&self, p: f64) -> f64 {
177        if p <= 0. || 1. <= p {
178            panic!("p must be in [0, 1]");
179        };
180        if p <= 0.5 {
181            self.location + self.scale * (2. * p).ln()
182        } else {
183            self.location - self.scale * (2. - 2. * p).ln()
184        }
185    }
186}
187
188impl Min<f64> for Laplace {
189    /// Returns the minimum value in the domain of the laplace
190    /// distribution representable by a double precision float
191    ///
192    /// # Formula
193    ///
194    /// ```text
195    /// NEG_INF
196    /// ```
197    fn min(&self) -> f64 {
198        f64::NEG_INFINITY
199    }
200}
201
202impl Max<f64> for Laplace {
203    /// Returns the maximum value in the domain of the laplace
204    /// distribution representable by a double precision float
205    ///
206    /// # Formula
207    ///
208    /// ```text
209    /// f64::INFINITY
210    /// ```
211    fn max(&self) -> f64 {
212        f64::INFINITY
213    }
214}
215
216impl Distribution<f64> for Laplace {
217    /// Returns the mode of the laplace distribution
218    ///
219    /// # Formula
220    ///
221    /// ```text
222    /// μ
223    /// ```
224    ///
225    /// where `μ` is the location
226    fn mean(&self) -> Option<f64> {
227        Some(self.location)
228    }
229
230    /// Returns the variance of the laplace distribution
231    ///
232    /// # Formula
233    ///
234    /// ```text
235    /// 2*b^2
236    /// ```
237    ///
238    /// where `b` is the scale
239    fn variance(&self) -> Option<f64> {
240        Some(2. * self.scale * self.scale)
241    }
242
243    /// Returns the entropy of the laplace distribution
244    ///
245    /// # Formula
246    ///
247    /// ```text
248    /// ln(2be)
249    /// ```
250    ///
251    /// where `b` is the scale
252    fn entropy(&self) -> Option<f64> {
253        Some((2. * self.scale).ln() + 1.)
254    }
255
256    /// Returns the skewness of the laplace distribution
257    ///
258    /// # Formula
259    ///
260    /// ```text
261    /// 0
262    /// ```
263    fn skewness(&self) -> Option<f64> {
264        Some(0.)
265    }
266}
267
268impl Median<f64> for Laplace {
269    /// Returns the median of the laplace distribution
270    ///
271    /// # Formula
272    ///
273    /// ```text
274    /// μ
275    /// ```
276    ///
277    /// where `μ` is the location
278    fn median(&self) -> f64 {
279        self.location
280    }
281}
282
283impl Mode<Option<f64>> for Laplace {
284    /// Returns the mode of the laplace distribution
285    ///
286    /// # Formula
287    ///
288    /// ```text
289    /// μ
290    /// ```
291    ///
292    /// where `μ` is the location
293    fn mode(&self) -> Option<f64> {
294        Some(self.location)
295    }
296}
297
298impl Continuous<f64, f64> for Laplace {
299    /// Calculates the probability density function for the laplace
300    /// distribution at `x`
301    ///
302    /// # Formula
303    ///
304    /// ```text
305    /// (1 / 2b) * exp(-|x - μ| / b)
306    /// ```
307    /// where `μ` is the location and `b` is the scale
308    fn pdf(&self, x: f64) -> f64 {
309        (-(x - self.location).abs() / self.scale).exp() / (2. * self.scale)
310    }
311
312    /// Calculates the log probability density function for the laplace
313    /// distribution at `x`
314    ///
315    /// # Formula
316    ///
317    /// ```text
318    /// ln((1 / 2b) * exp(-|x - μ| / b))
319    /// ```
320    ///
321    /// where `μ` is the location and `b` is the scale
322    fn ln_pdf(&self, x: f64) -> f64 {
323        ((-(x - self.location).abs() / self.scale).exp() / (2. * self.scale)).ln()
324    }
325}
326
327#[cfg(test)]
328mod tests {
329    use super::*;
330
331    use crate::testing_boiler;
332
333    testing_boiler!(location: f64, scale: f64; Laplace; LaplaceError);
334
335    // A wrapper for the `assert_relative_eq!` macro from the approx crate.
336    //
337    // `rtol` is the accepable relative error.  This function is for testing
338    // relative tolerance *only*.  It should not be used with `expected = 0`.
339    //
340    fn test_rel_close<F>(location: f64, scale: f64, expected: f64, rtol: f64, get_fn: F)
341    where
342        F: Fn(Laplace) -> f64,
343    {
344        let x = create_and_get(location, scale, get_fn);
345        assert_relative_eq!(expected, x, epsilon = 0.0, max_relative = rtol);
346    }
347
348    #[test]
349    fn test_create() {
350        create_ok(1.0, 2.0);
351        create_ok(f64::NEG_INFINITY, 0.1);
352        create_ok(-5.0 - 1.0, 1.0);
353        create_ok(0.0, 5.0);
354        create_ok(1.0, 7.0);
355        create_ok(5.0, 10.0);
356        create_ok(f64::INFINITY, f64::INFINITY);
357    }
358
359    #[test]
360    fn test_bad_create() {
361        test_create_err(2.0, -1.0, LaplaceError::ScaleInvalid);
362        test_create_err(f64::NAN, 1.0, LaplaceError::LocationInvalid);
363        create_err(f64::NAN, -1.0);
364    }
365
366    #[test]
367    fn test_mean() {
368        let mean = |x: Laplace| x.mean().unwrap();
369        test_exact(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, mean);
370        test_exact(-5.0 - 1.0, 1.0, -6.0, mean);
371        test_exact(0.0, 5.0, 0.0, mean);
372        test_exact(1.0, 10.0, 1.0, mean);
373        test_exact(f64::INFINITY, f64::INFINITY, f64::INFINITY, mean);
374    }
375
376    #[test]
377    fn test_variance() {
378        let variance = |x: Laplace| x.variance().unwrap();
379        test_absolute(f64::NEG_INFINITY, 0.1, 0.02, 1E-12, variance);
380        test_absolute(-5.0 - 1.0, 1.0, 2.0, 1E-12, variance);
381        test_absolute(0.0, 5.0, 50.0, 1E-12, variance);
382        test_absolute(1.0, 7.0, 98.0, 1E-12, variance);
383        test_absolute(5.0, 10.0, 200.0, 1E-12, variance);
384        test_absolute(f64::INFINITY, f64::INFINITY, f64::INFINITY, 1E-12, variance);
385    }
386    #[test]
387    fn test_entropy() {
388        let entropy = |x: Laplace| x.entropy().unwrap();
389        test_absolute(
390            f64::NEG_INFINITY,
391            0.1,
392            (2.0 * f64::consts::E * 0.1).ln(),
393            1E-12,
394            entropy,
395        );
396        test_absolute(-6.0, 1.0, (2.0 * f64::consts::E).ln(), 1E-12, entropy);
397        test_absolute(1.0, 7.0, (2.0 * f64::consts::E * 7.0).ln(), 1E-12, entropy);
398        test_absolute(5., 10., (2. * f64::consts::E * 10.).ln(), 1E-12, entropy);
399        test_absolute(f64::INFINITY, f64::INFINITY, f64::INFINITY, 1E-12, entropy);
400    }
401
402    #[test]
403    fn test_skewness() {
404        let skewness = |x: Laplace| x.skewness().unwrap();
405        test_exact(f64::NEG_INFINITY, 0.1, 0.0, skewness);
406        test_exact(-6.0, 1.0, 0.0, skewness);
407        test_exact(1.0, 7.0, 0.0, skewness);
408        test_exact(5.0, 10.0, 0.0, skewness);
409        test_exact(f64::INFINITY, f64::INFINITY, 0.0, skewness);
410    }
411
412    #[test]
413    fn test_mode() {
414        let mode = |x: Laplace| x.mode().unwrap();
415        test_exact(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, mode);
416        test_exact(-6.0, 1.0, -6.0, mode);
417        test_exact(1.0, 7.0, 1.0, mode);
418        test_exact(5.0, 10.0, 5.0, mode);
419        test_exact(f64::INFINITY, f64::INFINITY, f64::INFINITY, mode);
420    }
421
422    #[test]
423    fn test_median() {
424        let median = |x: Laplace| x.median();
425        test_exact(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, median);
426        test_exact(-6.0, 1.0, -6.0, median);
427        test_exact(1.0, 7.0, 1.0, median);
428        test_exact(5.0, 10.0, 5.0, median);
429        test_exact(f64::INFINITY, f64::INFINITY, f64::INFINITY, median);
430    }
431
432    #[test]
433    fn test_min() {
434        test_exact(0.0, 1.0, f64::NEG_INFINITY, |l| l.min());
435    }
436
437    #[test]
438    fn test_max() {
439        test_exact(0.0, 1.0, f64::INFINITY, |l| l.max());
440    }
441
442    #[test]
443    fn test_density() {
444        let pdf = |arg: f64| move |x: Laplace| x.pdf(arg);
445        test_absolute(0.0, 0.1, 1.529511602509129e-06, 1E-12, pdf(1.5));
446        test_absolute(1.0, 0.1, 7.614989872356341e-08, 1E-12, pdf(2.8));
447        test_absolute(-1.0, 0.1, 3.8905661205668983e-19, 1E-12, pdf(-5.4));
448        test_absolute(5.0, 0.1, 5.056107463052243e-43, 1E-12, pdf(-4.9));
449        test_absolute(-5.0, 0.1, 1.9877248679543235e-30, 1E-12, pdf(2.0));
450        test_absolute(f64::INFINITY, 0.1, 0.0, 1E-12, pdf(5.5));
451        test_absolute(f64::NEG_INFINITY, 0.1, 0.0, 1E-12, pdf(-0.0));
452        test_absolute(0.0, 1.0, 0.0, 1E-12, pdf(f64::INFINITY));
453        test_absolute(1.0, 1.0, 0.00915781944436709, 1E-12, pdf(5.0));
454        test_absolute(-1.0, 1.0, 0.5, 1E-12, pdf(-1.0));
455        test_absolute(5.0, 1.0, 0.0012393760883331792, 1E-12, pdf(-1.0));
456        test_absolute(-5.0, 1.0, 0.0002765421850739168, 1E-12, pdf(2.5));
457        test_absolute(f64::INFINITY, 0.1, 0.0, 1E-12, pdf(2.0));
458        test_absolute(f64::NEG_INFINITY, 0.1, 0.0, 1E-12, pdf(15.0));
459        test_absolute(0.0, f64::INFINITY, 0.0, 1E-12, pdf(89.3));
460        test_absolute(1.0, f64::INFINITY, 0.0, 1E-12, pdf(-0.1));
461        test_absolute(-1.0, f64::INFINITY, 0.0, 1E-12, pdf(0.1));
462        test_absolute(5.0, f64::INFINITY, 0.0, 1E-12, pdf(-6.1));
463        test_absolute(-5.0, f64::INFINITY, 0.0, 1E-12, pdf(-10.0));
464        test_is_nan(f64::INFINITY, f64::INFINITY, pdf(2.0));
465        test_is_nan(f64::NEG_INFINITY, f64::INFINITY, pdf(-5.1));
466    }
467
468    #[test]
469    fn test_ln_density() {
470        let ln_pdf = |arg: f64| move |x: Laplace| x.ln_pdf(arg);
471        test_absolute(0.0, 0.1, -13.3905620875659, 1E-12, ln_pdf(1.5));
472        test_absolute(1.0, 0.1, -16.390562087565897, 1E-12, ln_pdf(2.8));
473        test_absolute(-1.0, 0.1, -42.39056208756591, 1E-12, ln_pdf(-5.4));
474        test_absolute(5.0, 0.1, -97.3905620875659, 1E-12, ln_pdf(-4.9));
475        test_absolute(-5.0, 0.1, -68.3905620875659, 1E-12, ln_pdf(2.0));
476        test_exact(f64::INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(5.5));
477        test_exact(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(-0.0));
478        test_exact(0.0, 1.0, f64::NEG_INFINITY, ln_pdf(f64::INFINITY));
479        test_absolute(1.0, 1.0, -4.693147180559945, 1E-12, ln_pdf(5.0));
480        test_absolute(-1.0, 1.0, -f64::consts::LN_2, 1E-12, ln_pdf(-1.0));
481        test_absolute(5.0, 1.0, -6.693147180559945, 1E-12, ln_pdf(-1.0));
482        test_absolute(-5.0, 1.0, -8.193147180559945, 1E-12, ln_pdf(2.5));
483        test_exact(f64::INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(2.0));
484        test_exact(f64::NEG_INFINITY, 0.1, f64::NEG_INFINITY, ln_pdf(15.0));
485        test_exact(0.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(89.3));
486        test_exact(1.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-0.1));
487        test_exact(-1.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(0.1));
488        test_exact(5.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-6.1));
489        test_exact(-5.0, f64::INFINITY, f64::NEG_INFINITY, ln_pdf(-10.0));
490        test_is_nan(f64::INFINITY, f64::INFINITY, ln_pdf(2.0));
491        test_is_nan(f64::NEG_INFINITY, f64::INFINITY, ln_pdf(-5.1));
492    }
493
494    #[test]
495    fn test_cdf() {
496        let cdf = |arg: f64| move |x: Laplace| x.cdf(arg);
497        let loc = 0.0f64;
498        let scale = 1.0f64;
499        let reltol = 1e-15f64;
500
501        // Expected value from Wolfram Alpha: CDF[LaplaceDistribution[0, 1], 1/2].
502        let expected = 0.69673467014368328819810023250440977f64;
503        test_rel_close(loc, scale, expected, reltol, cdf(0.5));
504
505        // Wolfram Alpha: CDF[LaplaceDistribution[0, 1], -1/2]
506        let expected = 0.30326532985631671180189976749559023f64;
507        test_rel_close(loc, scale, expected, reltol, cdf(-0.5));
508
509        // Wolfram Alpha: CDF[LaplaceDistribution[0, 1], -100]
510        let expected = 1.8600379880104179814798479019315592e-44f64;
511        test_rel_close(loc, scale, expected, reltol, cdf(-100.0));
512    }
513
514    #[test]
515    fn test_sf() {
516        let sf = |arg: f64| move |x: Laplace| x.sf(arg);
517        let loc = 0.0f64;
518        let scale = 1.0f64;
519        let reltol = 1e-15f64;
520
521        // Expected value from Wolfram Alpha: SurvivalFunction[LaplaceDistribution[0, 1], 1/2].
522        let expected = 0.30326532985631671180189976749559022f64;
523        test_rel_close(loc, scale, expected, reltol, sf(0.5));
524
525        // Wolfram Alpha: SurvivalFunction[LaplaceDistribution[0, 1], -1/2]
526        let expected = 0.69673467014368328819810023250440977f64;
527        test_rel_close(loc, scale, expected, reltol, sf(-0.5));
528
529        // Wolfram Alpha: SurvivalFunction[LaplaceDistribution[0, 1], 100]
530        let expected = 1.86003798801041798147984790193155916e-44;
531        test_rel_close(loc, scale, expected, reltol, sf(100.0));
532    }
533
534    #[test]
535    fn test_inverse_cdf() {
536        let inverse_cdf = |arg: f64| move |x: Laplace| x.inverse_cdf(arg);
537        let loc = 0.0f64;
538        let scale = 1.0f64;
539        let reltol = 1e-15f64;
540
541        // Wolfram Alpha: Inverse CDF[LaplaceDistribution[0, 1], 1/10000000000]
542        let expected = -22.3327037493805115307626824253854655f64;
543        test_rel_close(loc, scale, expected, reltol, inverse_cdf(1e-10));
544
545        // Wolfram Alpha: Inverse CDF[LaplaceDistribution[0, 1], 1/1000].
546        let expected = -6.2146080984221917426367422425949161f64;
547        test_rel_close(loc, scale, expected, reltol, inverse_cdf(0.001));
548
549        // Wolfram Alpha: Inverse CDF[LaplaceDistribution[0, 1], 95/100]
550        let expected = 2.3025850929940456840179914546843642f64;
551        test_rel_close(loc, scale, expected, reltol, inverse_cdf(0.95));
552    }
553
554    #[cfg(feature = "rand")]
555    #[test]
556    fn test_sample() {
557        use ::rand::distributions::Distribution;
558        use ::rand::thread_rng;
559
560        let l = create_ok(0.1, 0.5);
561        l.sample(&mut thread_rng());
562    }
563
564    #[cfg(feature = "rand")]
565    #[test]
566    fn test_sample_distribution() {
567        use ::rand::distributions::Distribution;
568        use ::rand::rngs::StdRng;
569        use ::rand::SeedableRng;
570
571        // sanity check sampling
572        let location = 0.0;
573        let scale = 1.0;
574        let n = create_ok(location, scale);
575        let trials = 10_000;
576        let tolerance = 250;
577
578        for seed in 0..10 {
579            let mut r: StdRng = SeedableRng::seed_from_u64(seed);
580
581            let result = (0..trials).map(|_| n.sample(&mut r)).fold(0, |sum, val| {
582                if val > 0.0 {
583                    sum + 1
584                } else if val < 0.0 {
585                    sum - 1
586                } else {
587                    0
588                }
589            });
590            assert!(
591                result > -tolerance && result < tolerance,
592                "Balance is {result} for seed {seed}"
593            );
594        }
595    }
596}
statrs/distribution/laplace.rs

statrs/distribution/
laplace.rs
