rv/dist/
laplace.rs

1//! Laplace (double exponential) distribution
2#[cfg(feature = "serde1")]
3use serde::{Deserialize, Serialize};
4
5use crate::impl_display;
6use crate::traits::{
7    Cdf, ContinuousDistr, Entropy, HasDensity, Kurtosis, Mean, Median, Mode,
8    Parameterized, Sampleable, Scalable, Shiftable, Skewness, Support,
9    Variance,
10};
11use rand::Rng;
12use std::f64::consts::{E, FRAC_1_SQRT_2, LN_2};
13use std::fmt;
14
15/// [Laplace](https://en.wikipedia.org/wiki/Laplace_distribution), or double
16/// exponential, distribution over x in (-∞, ∞).
17///
18/// # Example
19///
20/// ```
21/// use rv::prelude::*;
22///
23/// let laplace = Laplace::new(0.0, 1.0).expect("Invalid params");
24///
25/// // 100 draws from Laplace
26/// let mut rng = rand::rng();
27/// let xs: Vec<f64> = laplace.sample(100, &mut rng);
28/// assert_eq!(xs.len(), 100);
29/// ```
30#[derive(Debug, Clone, PartialEq)]
31#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
32#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
33pub struct Laplace {
34    /// Location in (-∞, ∞)
35    mu: f64,
36    /// Scale in (0, ∞)
37    b: f64,
38}
39
40impl Shiftable for Laplace {
41    type Output = Laplace;
42    type Error = LaplaceError;
43
44    fn shifted(self, shift: f64) -> Result<Self::Output, Self::Error>
45    where
46        Self: Sized,
47    {
48        Laplace::new(self.mu() + shift, self.b())
49    }
50
51    fn shifted_unchecked(self, shift: f64) -> Self::Output
52    where
53        Self: Sized,
54    {
55        Laplace::new_unchecked(self.mu() + shift, self.b())
56    }
57}
58
59impl Scalable for Laplace {
60    type Output = Laplace;
61    type Error = LaplaceError;
62
63    fn scaled(self, scale: f64) -> Result<Self::Output, Self::Error>
64    where
65        Self: Sized,
66    {
67        Laplace::new(self.mu() * scale, self.b() * scale)
68    }
69
70    fn scaled_unchecked(self, scale: f64) -> Self::Output
71    where
72        Self: Sized,
73    {
74        Laplace::new_unchecked(self.mu() * scale, self.b() * scale)
75    }
76}
77
78pub struct LaplaceParameters {
79    pub mu: f64,
80    pub b: f64,
81}
82
83impl Parameterized for Laplace {
84    type Parameters = LaplaceParameters;
85
86    fn emit_params(&self) -> Self::Parameters {
87        Self::Parameters {
88            mu: self.mu(),
89            b: self.b(),
90        }
91    }
92
93    fn from_params(params: Self::Parameters) -> Self {
94        Self::new_unchecked(params.mu, params.b)
95    }
96}
97
98#[derive(Debug, Clone, PartialEq)]
99#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
100#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
101pub enum LaplaceError {
102    /// The mu parameter is infinite or NaN
103    MuNotFinite { mu: f64 },
104    /// The b parameter less than or equal to zero
105    BTooLow { b: f64 },
106    /// The b parameter is infinite or NaN
107    BNotFinite { b: f64 },
108}
109
110impl Laplace {
111    /// Create a new Laplace distribution.
112    pub fn new(mu: f64, b: f64) -> Result<Self, LaplaceError> {
113        if !mu.is_finite() {
114            Err(LaplaceError::MuNotFinite { mu })
115        } else if !b.is_finite() {
116            Err(LaplaceError::BNotFinite { b })
117        } else if b <= 0.0 {
118            Err(LaplaceError::BTooLow { b })
119        } else {
120            Ok(Laplace { mu, b })
121        }
122    }
123
124    /// Creates a new Laplace without checking whether the parameters are
125    /// valid.
126    #[inline]
127    #[must_use]
128    pub fn new_unchecked(mu: f64, b: f64) -> Self {
129        Laplace { mu, b }
130    }
131
132    /// Get the mu parameter
133    ///
134    /// # Examples
135    ///
136    /// ```rust
137    /// # use rv::dist::Laplace;
138    /// let laplace = Laplace::new(-1.0, 2.0).unwrap();
139    /// assert_eq!(laplace.mu(), -1.0);
140    /// ```
141    #[inline]
142    #[must_use]
143    pub fn mu(&self) -> f64 {
144        self.mu
145    }
146
147    /// Set the value of the mu parameter
148    ///
149    /// # Example
150    /// ```rust
151    /// # use rv::dist::Laplace;
152    /// let mut laplace = Laplace::new(-1.0, 2.0).unwrap();
153    /// assert_eq!(laplace.mu(), -1.0);
154    ///
155    /// laplace.set_mu(2.3).unwrap();
156    /// assert_eq!(laplace.mu(), 2.3);
157    /// ```
158    ///
159    /// Will error for invalid values
160    ///
161    /// ```rust
162    /// # use rv::dist::Laplace;
163    /// # let mut laplace = Laplace::new(-1.0, 2.0).unwrap();
164    /// assert!(laplace.set_mu(0.0).is_ok());
165    /// assert!(laplace.set_mu(f64::INFINITY).is_err());
166    /// assert!(laplace.set_mu(f64::NEG_INFINITY).is_err());
167    /// assert!(laplace.set_mu(f64::NAN).is_err());
168    /// ```
169    #[inline]
170    pub fn set_mu(&mut self, mu: f64) -> Result<(), LaplaceError> {
171        if mu.is_finite() {
172            self.set_mu_unchecked(mu);
173            Ok(())
174        } else {
175            Err(LaplaceError::MuNotFinite { mu })
176        }
177    }
178
179    /// Set the value of the mu parameter without input validation
180    #[inline]
181    pub fn set_mu_unchecked(&mut self, mu: f64) {
182        self.mu = mu;
183    }
184
185    /// Get the b parameter
186    ///
187    /// # Examples
188    ///
189    /// ```rust
190    /// # use rv::dist::Laplace;
191    /// let laplace = Laplace::new(-1.0, 2.0).unwrap();
192    /// assert_eq!(laplace.b(), 2.0);
193    /// ```
194    #[inline]
195    #[must_use]
196    pub fn b(&self) -> f64 {
197        self.b
198    }
199
200    /// Set the value of the b parameter
201    ///
202    /// # Example
203    /// ```rust
204    /// # use rv::dist::Laplace;
205    /// let mut laplace = Laplace::new(-1.0, 2.0).unwrap();
206    /// assert_eq!(laplace.b(), 2.0);
207    ///
208    /// laplace.set_b(2.3).unwrap();
209    /// assert_eq!(laplace.b(), 2.3);
210    /// ```
211    ///
212    /// Will error for invalid values
213    ///
214    /// ```rust
215    /// # use rv::dist::Laplace;
216    /// # let mut laplace = Laplace::new(-1.0, 2.0).unwrap();
217    /// assert!(laplace.set_b(2.3).is_ok());
218    /// assert!(laplace.set_b(0.0).is_err());
219    /// assert!(laplace.set_b(f64::INFINITY).is_err());
220    /// assert!(laplace.set_b(f64::NEG_INFINITY).is_err());
221    /// assert!(laplace.set_b(f64::NAN).is_err());
222    /// ```
223    #[inline]
224    pub fn set_b(&mut self, b: f64) -> Result<(), LaplaceError> {
225        if b <= 0.0 {
226            Err(LaplaceError::BTooLow { b })
227        } else if !b.is_finite() {
228            Err(LaplaceError::BNotFinite { b })
229        } else {
230            self.set_b_unchecked(b);
231            Ok(())
232        }
233    }
234
235    /// Set the value of the b parameter without input validation
236    #[inline]
237    pub fn set_b_unchecked(&mut self, b: f64) {
238        self.b = b;
239    }
240}
241
242/// Laplace with mean 0 and variance 1
243impl Default for Laplace {
244    fn default() -> Self {
245        Laplace::new_unchecked(0.0, FRAC_1_SQRT_2)
246    }
247}
248
249impl From<&Laplace> for String {
250    fn from(laplace: &Laplace) -> String {
251        format!("Laplace(μ: {}, b: {})", laplace.mu, laplace.b)
252    }
253}
254
255impl_display!(Laplace);
256
257#[inline]
258fn laplace_partial_draw(u: f64) -> f64 {
259    let r = u - 0.5;
260    r.signum() * 2.0_f64.mul_add(-r.abs(), 1.0).ln()
261}
262
263macro_rules! impl_traits {
264    ($kind:ty) => {
265        impl HasDensity<$kind> for Laplace {
266            fn ln_f(&self, x: &$kind) -> f64 {
267                // TODO: could cache ln(b)
268                -(f64::from(*x) - self.mu).abs() / self.b - self.b.ln() - LN_2
269            }
270        }
271
272        impl Sampleable<$kind> for Laplace {
273            fn draw<R: Rng>(&self, rng: &mut R) -> $kind {
274                let u = rng.sample(rand_distr::OpenClosed01);
275                self.b.mul_add(-laplace_partial_draw(u), self.mu) as $kind
276            }
277        }
278
279        impl Support<$kind> for Laplace {
280            fn supports(&self, x: &$kind) -> bool {
281                x.is_finite()
282            }
283        }
284
285        impl ContinuousDistr<$kind> for Laplace {}
286
287        impl Cdf<$kind> for Laplace {
288            fn cdf(&self, x: &$kind) -> f64 {
289                let xf: f64 = f64::from(*x);
290                if xf < self.mu {
291                    0.5 * ((xf - self.mu) / self.b).exp()
292                } else {
293                    0.5_f64.mul_add(-(-(xf - self.mu) / self.b).exp(), 1.0)
294                }
295            }
296        }
297
298        impl Mean<$kind> for Laplace {
299            fn mean(&self) -> Option<$kind> {
300                Some(self.mu as $kind)
301            }
302        }
303
304        impl Median<$kind> for Laplace {
305            fn median(&self) -> Option<$kind> {
306                Some(self.mu as $kind)
307            }
308        }
309
310        impl Mode<$kind> for Laplace {
311            fn mode(&self) -> Option<$kind> {
312                Some(self.mu as $kind)
313            }
314        }
315
316        impl Variance<$kind> for Laplace {
317            fn variance(&self) -> Option<$kind> {
318                Some((2.0 * self.b * self.b) as $kind)
319            }
320        }
321    };
322}
323
324impl Skewness for Laplace {
325    fn skewness(&self) -> Option<f64> {
326        Some(0.0)
327    }
328}
329
330impl Kurtosis for Laplace {
331    fn kurtosis(&self) -> Option<f64> {
332        Some(3.0)
333    }
334}
335
336impl Entropy for Laplace {
337    fn entropy(&self) -> f64 {
338        (2.0 * self.b * E).ln()
339    }
340}
341
342impl_traits!(f64);
343impl_traits!(f32);
344
345impl std::error::Error for LaplaceError {}
346
347#[cfg_attr(coverage_nightly, coverage(off))]
348impl fmt::Display for LaplaceError {
349    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
350        match self {
351            Self::MuNotFinite { mu } => write!(f, "non-finite mu: {mu}"),
352            Self::BTooLow { b } => {
353                write!(f, "b ({b}) must be greater than zero")
354            }
355            Self::BNotFinite { b } => write!(f, "non-finite b: {b}"),
356        }
357    }
358}
359
360#[cfg(test)]
361mod tests {
362    use super::*;
363    use crate::misc::ks_test;
364    use crate::test_basic_impls;
365    use std::f64;
366
367    const TOL: f64 = 1E-12;
368    const KS_PVAL: f64 = 0.2;
369    const N_TRIES: usize = 5;
370
371    test_basic_impls!(f64, Laplace);
372
373    #[test]
374    fn new() {
375        let laplace = Laplace::new(0.0, 1.0).unwrap();
376        assert::close(laplace.mu, 0.0, TOL);
377        assert::close(laplace.b, 1.0, TOL);
378    }
379
380    #[test]
381    fn new_should_reject_non_finite_mu() {
382        assert!(Laplace::new(f64::NEG_INFINITY, 1.0).is_err());
383        assert!(Laplace::new(f64::INFINITY, 1.0).is_err());
384        assert!(Laplace::new(f64::NAN, 1.0).is_err());
385    }
386
387    #[test]
388    fn new_should_reject_negative_b() {
389        assert!(Laplace::new(0.0, 0.0).is_err());
390        assert!(Laplace::new(0.0, -1e-12).is_err());
391        assert!(Laplace::new(0.0, -1e12).is_err());
392    }
393
394    #[test]
395    fn new_should_reject_non_finite_b() {
396        assert!(Laplace::new(0.0, f64::NAN).is_err());
397        assert!(Laplace::new(0.0, f64::INFINITY).is_err());
398    }
399
400    #[test]
401    fn mean() {
402        let m: f64 = Laplace::new(1.2, 3.4).unwrap().mean().unwrap();
403        assert::close(m, 1.2, TOL);
404    }
405
406    #[test]
407    fn median() {
408        let m: f64 = Laplace::new(1.2, 3.4).unwrap().median().unwrap();
409        assert::close(m, 1.2, TOL);
410    }
411
412    #[test]
413    fn mode() {
414        let m: f64 = Laplace::new(1.2, 3.4).unwrap().mode().unwrap();
415        assert::close(m, 1.2, TOL);
416    }
417
418    #[test]
419    fn variance() {
420        let v: f64 = Laplace::new(1.2, 3.4).unwrap().variance().unwrap();
421        assert::close(v, 23.119_999_999_999_997, TOL);
422    }
423
424    #[test]
425    fn entropy() {
426        let h: f64 = Laplace::new(1.2, 3.4).unwrap().entropy();
427        assert::close(h, 2.916_922_612_182_061, TOL);
428    }
429
430    #[test]
431    fn skewness() {
432        let s: f64 = Laplace::new(1.2, 3.4).unwrap().skewness().unwrap();
433        assert::close(s, 0.0, TOL);
434    }
435
436    #[test]
437    fn kurtosis() {
438        let k: f64 = Laplace::new(1.2, 3.4).unwrap().kurtosis().unwrap();
439        assert::close(k, 3.0, TOL);
440    }
441
442    #[test]
443    fn cdf_at_mu() {
444        let laplace = Laplace::new(1.2, 3.4).unwrap();
445        let cdf = laplace.cdf(&1.2_f64);
446        assert::close(cdf, 0.5, TOL);
447    }
448
449    #[test]
450    fn cdf_below_mu() {
451        let laplace = Laplace::new(1.2, 3.4).unwrap();
452        let cdf = laplace.cdf(&0.0_f64);
453        assert::close(cdf, 0.351_309_261_331_497_76, TOL);
454    }
455
456    #[test]
457    fn cdf_above_mu() {
458        let laplace = Laplace::new(1.2, 3.4).unwrap();
459        let cdf = laplace.cdf(&3.0_f64);
460        assert::close(cdf, 0.705_524_345_124_723_3, TOL);
461    }
462
463    #[test]
464    fn ln_pdf() {
465        let laplace = Laplace::new(1.2, 3.4).unwrap();
466        assert::close(laplace.ln_pdf(&1.2), -1.916_922_612_182_061, TOL);
467        assert::close(laplace.ln_pdf(&0.2), -2.211_040_259_240_884_4, TOL);
468    }
469
470    #[test]
471    fn draw_test() {
472        // Since we've had to implement the laplace draw ourselves, we have to
473        // make sure the thing works, so we use the Kolmogorov-Smirnov test.
474        let mut rng = rand::rng();
475        let laplace = Laplace::new(1.2, 3.4).unwrap();
476        let cdf = |x: f64| laplace.cdf(&x);
477
478        // test is flaky, try a few times
479        let passes = (0..N_TRIES).fold(0, |acc, _| {
480            let xs: Vec<f64> = laplace.sample(1000, &mut rng);
481            let (_, p) = ks_test(&xs, cdf);
482            if p > KS_PVAL { acc + 1 } else { acc }
483        });
484        assert!(passes > 0);
485    }
486    use crate::test_shiftable_cdf;
487    use crate::test_shiftable_density;
488    use crate::test_shiftable_entropy;
489    use crate::test_shiftable_method;
490
491    test_shiftable_method!(Laplace::new(2.0, 1.0).unwrap(), mean);
492    test_shiftable_method!(Laplace::new(2.0, 1.0).unwrap(), median);
493    test_shiftable_method!(Laplace::new(2.0, 1.0).unwrap(), mode);
494    test_shiftable_method!(Laplace::new(2.0, 1.0).unwrap(), variance);
495    test_shiftable_method!(Laplace::new(2.0, 1.0).unwrap(), skewness);
496    test_shiftable_method!(Laplace::new(2.0, 1.0).unwrap(), kurtosis);
497    test_shiftable_density!(Laplace::new(2.0, 1.0).unwrap());
498    test_shiftable_entropy!(Laplace::new(2.0, 1.0).unwrap());
499    test_shiftable_cdf!(Laplace::new(2.0, 1.0).unwrap());
500
501    use crate::test_scalable_cdf;
502    use crate::test_scalable_density;
503    use crate::test_scalable_entropy;
504    use crate::test_scalable_method;
505
506    test_scalable_method!(Laplace::new(2.0, 1.0).unwrap(), mean);
507    test_scalable_method!(Laplace::new(2.0, 1.0).unwrap(), median);
508    test_scalable_method!(Laplace::new(2.0, 1.0).unwrap(), mode);
509    test_scalable_method!(Laplace::new(2.0, 1.0).unwrap(), variance);
510    test_scalable_method!(Laplace::new(2.0, 1.0).unwrap(), skewness);
511    test_scalable_method!(Laplace::new(2.0, 1.0).unwrap(), kurtosis);
512    test_scalable_density!(Laplace::new(2.0, 1.0).unwrap());
513    test_scalable_entropy!(Laplace::new(2.0, 1.0).unwrap());
514    test_scalable_cdf!(Laplace::new(2.0, 1.0).unwrap());
515
516    #[test]
517    fn emit_and_from_params_are_identity() {
518        let dist_a = Laplace::new(3.0, 4.0).unwrap();
519        let dist_b = Laplace::from_params(dist_a.emit_params());
520        assert_eq!(dist_a, dist_b);
521    }
522}
rv/dist/laplace.rs

rv/dist/
laplace.rs
