statrs/distribution/
laplace.rs

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

statrs/distribution/
laplace.rs
