statrs/distribution/
gumbel.rs

1use super::{Continuous, ContinuousCDF};
2use crate::consts::EULER_MASCHERONI;
3use crate::statistics::*;
4use std::f64::{self, consts::PI};
5
6/// Implements the [Gumbel](https://en.wikipedia.org/wiki/Gumbel_distribution)
7/// distribution, also known as the type-I generalized extreme value distribution.
8///
9/// # Examples
10///
11/// ```
12/// use statrs::distribution::{Gumbel, Continuous};
13/// use statrs::{consts::EULER_MASCHERONI, statistics::Distribution};
14///
15/// let n = Gumbel::new(0.0, 1.0).unwrap();
16/// assert_eq!(n.location(), 0.0);
17/// assert_eq!(n.skewness().unwrap(), 1.13955);
18/// assert_eq!(n.pdf(0.0), 0.36787944117144233);
19/// ```
20#[derive(Copy, Clone, PartialEq, Debug)]
21pub struct Gumbel {
22    location: f64,
23    scale: f64,
24}
25
26/// Represents the errors that can occur when creating a [`Gumbel`]
27#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
28#[non_exhaustive]
29pub enum GumbelError {
30    /// The location is invalid (NAN)
31    LocationInvalid,
32
33    /// The scale is NAN, zero or less than zero
34    ScaleInvalid,
35}
36
37impl std::fmt::Display for GumbelError {
38    #[cfg_attr(coverage_nightly, coverage(off))]
39    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
40        match self {
41            GumbelError::LocationInvalid => write!(f, "Location is NAN"),
42            GumbelError::ScaleInvalid => write!(f, "Scale is NAN, zero or less than zero"),
43        }
44    }
45}
46
47impl std::error::Error for GumbelError {}
48
49impl Gumbel {
50    /// Constructs a new Gumbel distribution with the given
51    /// location and scale.
52    ///
53    /// # Errors
54    ///
55    /// Returns an error if location or scale are `NaN` or `scale <= 0.0`
56    ///
57    /// # Examples
58    ///
59    /// ```
60    /// use statrs::distribution::Gumbel;
61    ///
62    /// let mut result = Gumbel::new(0.0, 1.0);
63    /// assert!(result.is_ok());
64    ///
65    /// result = Gumbel::new(0.0, -1.0);
66    /// assert!(result.is_err());
67    /// ```
68    pub fn new(location: f64, scale: f64) -> Result<Self, GumbelError> {
69        if location.is_nan() {
70            return Err(GumbelError::LocationInvalid);
71        }
72
73        if scale.is_nan() || scale <= 0.0 {
74            return Err(GumbelError::ScaleInvalid);
75        }
76
77        Ok(Self { location, scale })
78    }
79
80    /// Returns the location of the Gumbel distribution
81    ///
82    /// # Examples
83    ///
84    /// ```
85    /// use statrs::distribution::Gumbel;
86    ///
87    /// let n = Gumbel::new(0.0, 1.0).unwrap();
88    /// assert_eq!(n.location(), 0.0);
89    /// ```
90    pub fn location(&self) -> f64 {
91        self.location
92    }
93
94    /// Returns the scale of the Gumbel distribution
95    ///
96    /// # Examples
97    ///
98    /// ```
99    /// use statrs::distribution::Gumbel;
100    ///
101    /// let n = Gumbel::new(0.0, 1.0).unwrap();
102    /// assert_eq!(n.scale(), 1.0);
103    /// ```
104    pub fn scale(&self) -> f64 {
105        self.scale
106    }
107}
108
109impl std::fmt::Display for Gumbel {
110    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
111        write!(f, "Gumbel({:?}, {:?})", self.location, self.scale)
112    }
113}
114
115#[cfg(feature = "rand")]
116impl ::rand::distributions::Distribution<f64> for Gumbel {
117    fn sample<R: rand::Rng + ?Sized>(&self, r: &mut R) -> f64 {
118        self.location - self.scale * ((-(r.gen::<f64>())).ln()).ln()
119    }
120}
121
122impl ContinuousCDF<f64, f64> for Gumbel {
123    /// Calculates the cumulative distribution function for the
124    /// Gumbel distribution at `x`
125    ///
126    /// # Formula
127    ///
128    /// ```text
129    /// e^(-e^(-(x - μ) / β))
130    /// ```
131    ///
132    /// where `μ` is the location and `β` is the scale
133    fn cdf(&self, x: f64) -> f64 {
134        (-(-(x - self.location) / self.scale).exp()).exp()
135    }
136
137    /// Calculates the inverse cumulative distribution function for the
138    /// Gumbel distribution at `x`
139    ///
140    /// # Formula
141    ///
142    /// ```text
143    /// μ - β ln(-ln(p)) where 0 < p < 1
144    /// -INF             where p <= 0
145    /// INF              otherwise
146    /// ```
147    ///
148    /// where `μ` is the location and `β` is the scale
149    fn inverse_cdf(&self, p: f64) -> f64 {
150        if p <= 0.0 {
151            f64::NEG_INFINITY
152        } else if p >= 1.0 {
153            f64::INFINITY
154        } else {
155            self.location - self.scale * ((-(p.ln())).ln())
156        }
157    }
158
159    /// Calculates the survival function for the
160    /// Gumbel distribution at `x`
161    ///
162    /// # Formula
163    ///
164    /// ```text
165    /// 1 - e^(-e^(-(x - μ) / β))
166    /// ```
167    ///
168    /// where `μ` is the location and `β` is the scale
169    fn sf(&self, x: f64) -> f64 {
170        -(-(-(x - self.location) / self.scale).exp()).exp_m1()
171    }
172}
173
174impl Min<f64> for Gumbel {
175    /// Returns the minimum value in the domain of the Gumbel
176    /// distribution representable by a double precision float
177    ///
178    /// # Formula
179    ///
180    /// ```text
181    /// NEG_INF
182    /// ```
183    fn min(&self) -> f64 {
184        f64::NEG_INFINITY
185    }
186}
187
188impl Max<f64> for Gumbel {
189    /// Returns the maximum value in the domain of the Gumbel
190    /// distribution representable by a double precision float
191    ///
192    /// # Formula
193    ///
194    /// ```text
195    /// f64::INFINITY
196    /// ```
197    fn max(&self) -> f64 {
198        f64::INFINITY
199    }
200}
201
202impl Distribution<f64> for Gumbel {
203    /// Returns the entropy of the Gumbel distribution
204    ///
205    /// # Formula
206    ///
207    /// ```text
208    /// ln(β) + γ + 1
209    /// ```
210    ///
211    /// where `β` is the scale
212    /// and `γ` is the Euler-Mascheroni constant (approx 0.57721)
213    fn entropy(&self) -> Option<f64> {
214        Some(1.0 + EULER_MASCHERONI + (self.scale).ln())
215    }
216
217    /// Returns the mean of the Gumbel distribution
218    ///
219    /// # Formula
220    ///
221    /// ```text
222    /// μ + γβ
223    /// ```
224    ///
225    /// where `μ` is the location, `β` is the scale
226    /// and `γ` is the Euler-Mascheroni constant (approx 0.57721)
227    fn mean(&self) -> Option<f64> {
228        Some(self.location + (EULER_MASCHERONI * self.scale))
229    }
230
231    /// Returns the skewness of the Gumbel distribution
232    ///
233    /// # Formula
234    ///
235    /// ```text
236    /// 12 * sqrt(6) * ζ(3) / π^3 ≈ 1.13955
237    /// ```
238    /// ζ(3) is the Riemann zeta function evaluated at 3 (approx 1.20206)
239    /// and π is the constant PI (approx 3.14159)
240    ///
241    /// This approximately evaluates to 1.13955
242    fn skewness(&self) -> Option<f64> {
243        Some(1.13955)
244    }
245
246    /// Returns the variance of the Gumbel distribution
247    ///
248    /// # Formula
249    ///
250    /// ```text
251    /// (π^2 / 6) * β^2
252    /// ```
253    ///
254    /// where `β` is the scale and `π` is the constant PI (approx 3.14159)
255    fn variance(&self) -> Option<f64> {
256        Some(((PI * PI) / 6.0) * self.scale * self.scale)
257    }
258
259    /// Returns the standard deviation of the Gumbel distribution
260    ///
261    /// # Formula
262    ///
263    /// ```text
264    /// β * π / sqrt(6)
265    /// ```
266    ///
267    /// where `β` is the scale and `π` is the constant PI (approx 3.14159)
268    fn std_dev(&self) -> Option<f64> {
269        Some(self.scale * PI / 6.0_f64.sqrt())
270    }
271}
272
273impl Median<f64> for Gumbel {
274    /// Returns the median of the Gumbel distribution
275    ///
276    /// # Formula
277    ///
278    /// ```text
279    /// μ - β ln(ln(2))
280    /// ```
281    ///
282    /// where `μ` is the location and `β` is the scale parameter
283    fn median(&self) -> f64 {
284        self.location - self.scale * (((2.0_f64).ln()).ln())
285    }
286}
287
288impl Mode<f64> for Gumbel {
289    /// Returns the mode of the Gumbel distribution
290    ///
291    /// # Formula
292    ///
293    /// ```text
294    /// μ
295    /// ```
296    ///
297    /// where `μ` is the location
298    fn mode(&self) -> f64 {
299        self.location
300    }
301}
302
303impl Continuous<f64, f64> for Gumbel {
304    /// Calculates the probability density function for the Gumbel
305    /// distribution at `x`
306    ///
307    /// # Formula
308    ///
309    /// ```text
310    /// (1/β) * exp(-(x - μ)/β) * exp(-exp(-(x - μ)/β))
311    /// ```
312    ///
313    /// where `μ` is the location, `β` is the scale
314    fn pdf(&self, x: f64) -> f64 {
315        (1.0_f64 / self.scale)
316            * (-(x - self.location) / (self.scale)).exp()
317            * (-((-(x - self.location) / self.scale).exp())).exp()
318    }
319
320    /// Calculates the log probability density function for the Gumbel
321    /// distribution at `x`
322    ///
323    /// # Formula
324    ///
325    /// ```text
326    /// ln((1/β) * exp(-(x - μ)/β) * exp(-exp(-(x - μ)/β)))
327    /// ```
328    ///
329    /// where `μ` is the location, `β` is the scale
330    fn ln_pdf(&self, x: f64) -> f64 {
331        ((1.0_f64 / self.scale)
332            * (-(x - self.location) / (self.scale)).exp()
333            * (-((-(x - self.location) / self.scale).exp())).exp())
334        .ln()
335    }
336}
337
338#[rustfmt::skip]
339#[cfg(test)]
340mod tests {
341    use super::*;
342    use crate::testing_boiler;
343
344    testing_boiler!(location: f64, scale: f64; Gumbel; GumbelError);
345
346    #[test]
347    fn test_create() {
348        create_ok(0.0, 0.1);
349        create_ok(0.0, 1.0);
350        create_ok(0.0, 10.0);
351        create_ok(10.0, 11.0);
352        create_ok(-5.0, 100.0);
353        create_ok(0.0, f64::INFINITY);
354    }
355
356    #[test]
357    fn test_bad_create() {
358        let invalid = [
359            (f64::NAN, 1.0, GumbelError::LocationInvalid),
360            (1.0, f64::NAN, GumbelError::ScaleInvalid),
361            (f64::NAN, f64::NAN, GumbelError::LocationInvalid),
362            (1.0, 0.0, GumbelError::ScaleInvalid),
363            (0.0, f64::NEG_INFINITY, GumbelError::ScaleInvalid)
364        ];
365
366        for (location, scale, err) in invalid {
367            test_create_err(location, scale, err);
368        }
369    }
370
371    #[test]
372    fn test_min_max() {
373        let min = |x: Gumbel| x.min();
374        let max = |x:Gumbel| x.max();
375
376        test_exact(0.0, 1.0, f64::NEG_INFINITY, min);
377        test_exact(0.0, 1.0, f64::INFINITY, max);
378    }
379
380    #[test]
381    fn test_entropy() {
382        let entropy = |x: Gumbel| x.entropy().unwrap();
383        test_exact(0.0, 2.0, 2.270362845461478, entropy);
384        test_exact(0.1, 4.0, 2.9635100260214235, entropy);
385        test_exact(1.0, 10.0, 3.8798007578955787, entropy);
386        test_exact(10.0, 11.0, 3.9751109376999034, entropy); 
387    }
388
389    #[test]
390    fn test_mean() {
391        let mean = |x: Gumbel| x.mean().unwrap();
392        test_exact(0.0, 2.0, 1.1544313298030658, mean);
393        test_exact(0.1, 4.0, 2.4088626596061316, mean);
394        test_exact(1.0, 10.0, 6.772156649015328, mean);
395        test_exact(10.0, 11.0, 16.34937231391686, mean);
396        test_exact(10.0, f64::INFINITY, f64::INFINITY, mean);
397    }
398
399    #[test]
400    fn test_skewness() {
401        let skewness = |x: Gumbel| x.skewness().unwrap();
402        test_exact(0.0, 2.0, 1.13955, skewness);
403        test_exact(0.1, 4.0, 1.13955, skewness);
404        test_exact(1.0, 10.0, 1.13955, skewness);
405        test_exact(10.0, 11.0, 1.13955, skewness);
406        test_exact(10.0, f64::INFINITY, 1.13955, skewness);
407    }
408
409    #[test]
410    fn test_variance() {
411        let variance = |x: Gumbel| x.variance().unwrap();
412        test_exact(0.0, 2.0, 6.579736267392906, variance);
413        test_exact(0.1, 4.0, 26.318945069571624, variance);
414        test_exact(1.0, 10.0, 164.49340668482265, variance);
415        test_exact(10.0, 11.0, 199.03702208863538, variance);
416    }
417
418    #[test]
419    fn test_std_dev() {
420        let std_dev = |x: Gumbel| x.std_dev().unwrap();
421        test_exact(0.0, 2.0, 2.565099660323728, std_dev);
422        test_exact(0.1, 4.0, 5.130199320647456, std_dev);
423        test_exact(1.0, 10.0, 12.82549830161864, std_dev);
424        test_exact(10.0, 11.0, 14.108048131780505, std_dev);
425    }
426
427    #[test]
428    fn test_median() {
429        let median = |x: Gumbel| x.median();
430        test_exact(0.0, 2.0, 0.7330258411633287, median);
431        test_exact(0.1, 4.0, 1.5660516823266574, median);
432        test_exact(1.0, 10.0, 4.665129205816644, median);
433        test_exact(10.0, 11.0, 14.031642126398307, median);
434        test_exact(10.0, f64::INFINITY, f64::INFINITY, median);
435    }
436
437    #[test]
438    fn test_mode() {
439        let mode = |x: Gumbel| x.mode();
440        test_exact(0.0, 2.0, 0.0, mode);
441        test_exact(0.1, 4.0, 0.1, mode);
442        test_exact(1.0, 10.0, 1.0, mode);
443        test_exact(10.0, 11.0, 10.0, mode);
444        test_exact(10.0, f64::INFINITY, 10.0, mode);
445    }
446
447    #[test]
448    fn test_cdf() {
449        let cdf = |a: f64| move |x: Gumbel| x.cdf(a);
450        test_exact(0.0, 0.1, 0.0, cdf(-5.0));
451        test_exact(0.0, 0.1, 0.0, cdf(-1.0));
452        test_exact(0.0, 0.1, 0.36787944117144233, cdf(0.0));
453        test_exact(0.0, 0.1, 0.9999546011007987, cdf(1.0));
454        test_absolute(0.0, 0.1, 0.99999999999999999, 1e-12, cdf(5.0));
455        test_absolute(0.0, 1.0, 0.06598803584531253, 1e-12, cdf(-1.0));
456        test_exact(0.0, 1.0, 0.36787944117144233, cdf(0.0));
457        test_absolute(0.0, 10.0, 0.192295645547964928, 1e-12, cdf(-5.0));
458        test_absolute(0.0, 10.0, 0.3311542771529088, 1e-12, cdf(-1.0));
459        test_exact(0.0, 10.0, 0.36787944117144233, cdf(0.0));
460        test_absolute(0.0, 10.0, 0.4046076616641318, 1e-12, cdf(1.0));
461        test_absolute(0.0, 10.0, 0.545239211892605, 1e-12, cdf(5.0));
462        test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(-5.0));
463        test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(-1.0));
464        test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(0.0));
465        test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(1.0));
466        test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(5.0));
467        test_exact(f64::INFINITY, 1.0, 0.0, cdf(-5.0));
468        test_exact(f64::INFINITY, 1.0, 0.0, cdf(-1.0));
469        test_exact(f64::INFINITY, 1.0, 0.0, cdf(0.0));
470        test_exact(f64::INFINITY, 1.0, 0.0, cdf(1.0));
471        test_exact(f64::INFINITY, 1.0, 0.0, cdf(5.0));
472    }
473
474    #[test]
475    fn test_inverse_cdf() {
476        let inv_cdf = |a: f64| move |x: Gumbel| x.inverse_cdf(a);
477        test_exact(0.0, 0.1, f64::NEG_INFINITY, inv_cdf(-5.0));
478        test_exact(0.0, 0.1, f64::NEG_INFINITY, inv_cdf(-1.0));
479        test_exact(0.0, 0.1, f64::NEG_INFINITY, inv_cdf(0.0));
480        test_exact(0.0, 0.1, f64::INFINITY, inv_cdf(1.0));
481        test_exact(0.0, 0.1, f64::INFINITY, inv_cdf(5.0));
482        test_absolute(0.0, 1.0, -0.8340324452479557, 1e-12, inv_cdf(0.1));
483        test_absolute(0.0, 10.0, 3.6651292058166436, 1e-12, inv_cdf(0.5));
484        test_absolute(0.0, 10.0, 22.503673273124456, 1e-12, inv_cdf(0.9));
485        test_exact(2.0, f64::INFINITY, f64::NEG_INFINITY, inv_cdf(0.1));
486        test_exact(-2.0, f64::INFINITY, f64::INFINITY, inv_cdf(0.5));
487        test_exact(f64::INFINITY, 1.0, f64::INFINITY, inv_cdf(0.1));
488    }
489
490    #[test]
491    fn test_sf() {
492        let sf = |a: f64| move |x: Gumbel| x.sf(a);
493        test_exact(0.0, 0.1, 1.0, sf(-5.0));
494        test_exact(0.0, 0.1, 1.0, sf(-1.0));
495        test_absolute(0.0, 0.1, 0.632120558828557678, 1e-12, sf(0.0));
496        test_absolute(0.0, 0.1, 0.000045398899201269, 1e-12, sf(1.0));
497        test_absolute(0.0, 1.0, 0.934011964154687462, 1e-12, sf(-1.0));
498        test_absolute(0.0, 1.0, 0.632120558828557678, 1e-12, sf(0.0));
499        test_absolute(0.0, 1.0, 0.3077993724446536, 1e-12, sf(1.0));
500        test_absolute(0.0, 10.0, 0.66884572284709110, 1e-12, sf(-1.0));
501        test_absolute(0.0, 10.0, 0.632120558828557678, 1e-12, sf(0.0));
502        test_absolute(0.0, 10.0, 0.595392338335868174, 1e-12, sf(1.0));
503        test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(-5.0));
504        test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(-1.0));
505        test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(0.0));
506        test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(1.0));
507        test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(5.0));
508        test_exact(f64::INFINITY, 1.0, 1.0, sf(-5.0));
509        test_exact(f64::INFINITY, 1.0, 1.0, sf(-1.0));
510        test_exact(f64::INFINITY, 1.0, 1.0, sf(0.0));
511        test_exact(f64::INFINITY, 1.0, 1.0, sf(1.0));
512        test_exact(f64::INFINITY, 1.0, 1.0, sf(5.0));
513        test_absolute(0.0, 1.0, 4.248354255291589e-18, 1e-32, sf(40.0));
514        test_absolute(0.0, 1.0, 1.804851387845415e-35, 1e-50, sf(80.0));
515    }
516
517    #[test]
518    fn test_pdf() {
519        let pdf = |a: f64| move |x: Gumbel| x.pdf(a);
520        test_exact(0.0, 0.1, 0.0, pdf(-5.0));
521        test_exact(0.0, 0.1, 3.678794411714423215, pdf(0.0));
522        test_absolute(0.0, 0.1, 0.0004539786865564, 1e-12, pdf(1.0));
523        test_absolute(0.0, 1.0, 0.1793740787340171, 1e-12, pdf(-1.0));
524        test_exact(0.0, 1.0, 0.36787944117144233, pdf(0.0));
525        test_absolute(0.0, 1.0, 0.25464638004358249, 1e-12, pdf(1.0));
526        test_absolute(0.0, 10.0, 0.031704192107794217, 1e-12, pdf(-5.0));
527        test_absolute(0.0, 10.0, 0.0365982076505757, 1e-12, pdf(-1.0));
528        test_exact(0.0, 10.0, 0.036787944117144233, pdf(0.0));
529        test_absolute(0.0, 10.0, 0.03661041518977401428, 1e-12, pdf(1.0));
530        test_absolute(0.0, 10.0, 0.033070429889041, 1e-12, pdf(5.0));
531        test_exact(-2.0, f64::INFINITY, 0.0, pdf(-5.0));
532        test_exact(-2.0, f64::INFINITY, 0.0, pdf(-1.0));
533        test_exact(-2.0, f64::INFINITY, 0.0, pdf(0.0));
534        test_exact(-2.0, f64::INFINITY, 0.0, pdf(1.0));
535        test_exact(-2.0, f64::INFINITY, 0.0, pdf(5.0));
536    }
537
538    #[test]
539    fn test_ln_pdf() {
540        let ln_pdf = |a: f64| move |x: Gumbel| x.ln_pdf(a);
541        test_exact(0.0, 0.1, 0.0_f64.ln(), ln_pdf(-5.0));
542        test_exact(0.0, 0.1, 3.678794411714423215_f64.ln(), ln_pdf(0.0));
543        test_absolute(0.0, 0.1, 0.0004539786865564_f64.ln(), 1e-12, ln_pdf(1.0));
544        test_absolute(0.0, 1.0, 0.1793740787340171_f64.ln(), 1e-12, ln_pdf(-1.0));
545        test_exact(0.0, 1.0, 0.36787944117144233_f64.ln(), ln_pdf(0.0));
546        test_absolute(0.0, 1.0, 0.25464638004358249_f64.ln(), 1e-12, ln_pdf(1.0));
547        test_absolute(0.0, 10.0, 0.031704192107794217_f64.ln(), 1e-12, ln_pdf(-5.0));
548        test_absolute(0.0, 10.0, 0.0365982076505757_f64.ln(), 1e-12, ln_pdf(-1.0));
549        test_exact(0.0, 10.0, 0.036787944117144233_f64.ln(), ln_pdf(0.0));
550        test_absolute(0.0, 10.0, 0.03661041518977401428_f64.ln(), 1e-12, ln_pdf(1.0));
551        test_absolute(0.0, 10.0, 0.033070429889041_f64.ln(), 1e-12, ln_pdf(5.0));
552        test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(-5.0));
553        test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(-1.0));
554        test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(0.0));
555        test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(1.0));
556        test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(5.0));
557    }
558}
statrs/distribution/gumbel.rs

statrs/distribution/
gumbel.rs
