rv/dist/
gev.rs

1#[cfg(feature = "serde1")]
2use serde::{Deserialize, Serialize};
3
4use crate::consts;
5use crate::impl_display;
6use crate::misc::gammafn;
7use crate::traits::*;
8use rand::Rng;
9use std::f32;
10use std::f64;
11use std::f64::consts::{LN_2, PI};
12use std::fmt;
13
14/// [Generalized Extreme Value Distribution](https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution)
15/// Gev(μ, σ, ξ) where the parameters are
16/// μ is location
17/// σ is the scale
18/// ξ is the shape
19///
20/// ```math
21/// f(x|μ, σ, ξ) = \frac{1}{σ} t(x)^{ξ + 1} e^{-t(x)}
22///
23/// t(x) = ⎰ (1 + ξ ((x - μ) / σ))^(-1/ξ) if ξ ≠ 0
24///        ⎱ e^{(μ - x) / σ}              if ξ = 0
25/// ```
26#[derive(Debug, Clone, PartialEq)]
27#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
28#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
29pub struct Gev {
30    loc: f64,
31    scale: f64,
32    shape: f64,
33}
34
35pub struct GevParameters {
36    pub loc: f64,
37    pub scale: f64,
38    pub shape: f64,
39}
40
41impl Parameterized for Gev {
42    type Parameters = GevParameters;
43
44    fn emit_params(&self) -> Self::Parameters {
45        Self::Parameters {
46            loc: self.loc(),
47            scale: self.scale(),
48            shape: self.shape(),
49        }
50    }
51
52    fn from_params(params: Self::Parameters) -> Self {
53        Self::new_unchecked(params.loc, params.scale, params.shape)
54    }
55}
56
57#[derive(Debug, Clone, PartialEq)]
58#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
59#[cfg_attr(feature = "serde1", serde(rename_all = "snake_case"))]
60pub enum GevError {
61    /// The location parameter is infinite or NaN
62    LocNotFinite { loc: f64 },
63    /// The shape parameter is infinite or NaN
64    ShapeNotFinite { shape: f64 },
65    /// The scale parameter is infinite or NaN
66    ScaleNotFinite { scale: f64 },
67    /// The scale parameter is less than or equal to zero
68    ScaleTooLow { scale: f64 },
69}
70
71impl Gev {
72    /// Create a new `Gev` distribution with location, scale, and shape.
73    pub fn new(loc: f64, scale: f64, shape: f64) -> Result<Self, GevError> {
74        if scale <= 0.0 {
75            Err(GevError::ScaleTooLow { scale })
76        } else if !scale.is_finite() {
77            Err(GevError::ScaleNotFinite { scale })
78        } else if !shape.is_finite() {
79            Err(GevError::ShapeNotFinite { shape })
80        } else if !loc.is_finite() {
81            Err(GevError::LocNotFinite { loc })
82        } else {
83            Ok(Gev { loc, scale, shape })
84        }
85    }
86
87    /// Creates a new Gev without checking whether the parameters are valid.
88    #[inline]
89    pub fn new_unchecked(loc: f64, scale: f64, shape: f64) -> Self {
90        Gev { loc, scale, shape }
91    }
92
93    /// Get the location parameter
94    ///
95    /// # Example
96    ///
97    /// ```rust
98    /// # use rv::dist::Gev;
99    /// let gev = Gev::new(1.2, 2.3, 3.4).unwrap();
100    ///
101    /// assert_eq!(gev.loc(), 1.2);
102    /// ```
103    #[inline]
104    pub fn loc(&self) -> f64 {
105        self.loc
106    }
107
108    /// Set the loc parameter without input validation
109    ///
110    /// # Example
111    ///
112    /// ```rust
113    /// # use rv::dist::Gev;
114    /// let mut gev = Gev::new(1.2, 2.3, 3.4).unwrap();
115    /// assert_eq!(gev.loc(), 1.2);
116    ///
117    /// gev.set_loc(2.8).unwrap();
118    /// assert_eq!(gev.loc(), 2.8);
119    /// ```
120    ///
121    /// Will error for invalid values
122    ///
123    /// ```rust
124    /// # use rv::dist::Gev;
125    /// # let mut gev = Gev::new(1.2, 2.3, 3.4).unwrap();
126    /// assert!(gev.set_loc(2.8).is_ok());
127    /// assert!(gev.set_loc(f64::INFINITY).is_err());
128    /// assert!(gev.set_loc(f64::NEG_INFINITY).is_err());
129    /// assert!(gev.set_loc(f64::NAN).is_err());
130    /// ```
131    #[inline]
132    pub fn set_loc(&mut self, loc: f64) -> Result<(), GevError> {
133        if !loc.is_finite() {
134            Err(GevError::LocNotFinite { loc })
135        } else {
136            self.set_loc_unchecked(loc);
137            Ok(())
138        }
139    }
140
141    /// Set the loc parameter without input validation
142    #[inline]
143    pub fn set_loc_unchecked(&mut self, loc: f64) {
144        self.loc = loc
145    }
146
147    /// Get the shape parameter
148    ///
149    /// # Example
150    ///
151    /// ```rust
152    /// # use rv::dist::Gev;
153    /// let gev = Gev::new(1.2, 2.3, 3.4).unwrap();
154    ///
155    /// assert_eq!(gev.shape(), 3.4);
156    /// ```
157    #[inline]
158    pub fn shape(&self) -> f64 {
159        self.shape
160    }
161
162    /// Set the shape parameter without input validation
163    ///
164    /// # Example
165    ///
166    /// ```rust
167    /// # use rv::dist::Gev;
168    /// let mut gev = Gev::new(1.2, 2.3, 3.4).unwrap();
169    /// assert_eq!(gev.shape(), 3.4);
170    ///
171    /// gev.set_shape(2.8).unwrap();
172    /// assert_eq!(gev.shape(), 2.8);
173    /// ```
174    ///
175    /// Will error for invalid values
176    ///
177    /// ```rust
178    /// # use rv::dist::Gev;
179    /// # let mut gev = Gev::new(1.2, 2.3, 3.4).unwrap();
180    /// assert!(gev.set_shape(2.8).is_ok());
181    /// assert!(gev.set_shape(f64::INFINITY).is_err());
182    /// assert!(gev.set_shape(f64::NEG_INFINITY).is_err());
183    /// assert!(gev.set_shape(f64::NAN).is_err());
184    /// ```
185    #[inline]
186    pub fn set_shape(&mut self, shape: f64) -> Result<(), GevError> {
187        if !shape.is_finite() {
188            Err(GevError::ShapeNotFinite { shape })
189        } else {
190            self.set_shape_unchecked(shape);
191            Ok(())
192        }
193    }
194
195    /// Set the shape parameter without input validation
196    #[inline]
197    pub fn set_shape_unchecked(&mut self, shape: f64) {
198        self.shape = shape
199    }
200
201    /// Get the scale parameter
202    ///
203    /// # Example
204    ///
205    /// ```rust
206    /// # use rv::dist::Gev;
207    /// let gev = Gev::new(1.2, 2.3, 3.4).unwrap();
208    ///
209    /// assert_eq!(gev.scale(), 2.3);
210    /// ```
211    #[inline]
212    pub fn scale(&self) -> f64 {
213        self.scale
214    }
215
216    /// Set the scale parameter without input validation
217    ///
218    /// # Example
219    ///
220    /// ```rust
221    /// # use rv::dist::Gev;
222    /// let mut gev = Gev::new(1.2, 2.3, 3.4).unwrap();
223    /// assert_eq!(gev.scale(), 2.3);
224    ///
225    /// gev.set_scale(2.8).unwrap();
226    /// assert_eq!(gev.scale(), 2.8);
227    /// ```
228    ///
229    /// Will error for invalid values
230    ///
231    /// ```rust
232    /// # use rv::dist::Gev;
233    /// # let mut gev = Gev::new(1.2, 2.3, 3.4).unwrap();
234    /// assert!(gev.set_scale(2.8).is_ok());
235    /// assert!(gev.set_scale(0.0).is_err());
236    /// assert!(gev.set_scale(-1.0).is_err());
237    /// assert!(gev.set_scale(f64::INFINITY).is_err());
238    /// assert!(gev.set_scale(f64::NEG_INFINITY).is_err());
239    /// assert!(gev.set_scale(f64::NAN).is_err());
240    /// ```
241    #[inline]
242    pub fn set_scale(&mut self, scale: f64) -> Result<(), GevError> {
243        if !scale.is_finite() {
244            Err(GevError::ScaleNotFinite { scale })
245        } else if scale <= 0.0 {
246            Err(GevError::ScaleTooLow { scale })
247        } else {
248            self.set_scale_unchecked(scale);
249            Ok(())
250        }
251    }
252
253    /// Set the scale parameter without input validation
254    #[inline]
255    pub fn set_scale_unchecked(&mut self, scale: f64) {
256        self.scale = scale
257    }
258}
259
260fn t(loc: f64, shape: f64, scale: f64, x: f64) -> f64 {
261    if shape == 0.0 {
262        ((loc - x) / scale).exp()
263    } else {
264        (1.0 + shape * (x - loc) / scale).powf(-1.0 / shape)
265    }
266}
267
268impl From<&Gev> for String {
269    fn from(gev: &Gev) -> String {
270        format!("GEV(μ: {}, α: {}, ξ: {})", gev.loc, gev.scale, gev.shape)
271    }
272}
273
274impl_display!(Gev);
275
276macro_rules! impl_traits {
277    ($kind: ty) => {
278        impl HasDensity<$kind> for Gev {
279            fn ln_f(&self, x: &$kind) -> f64 {
280                // TODO: could cache ln(scale)
281                let tv = t(self.loc, self.shape, self.scale, f64::from(*x));
282                (self.shape + 1.0).mul_add(tv.ln(), -self.scale.ln()) - tv
283            }
284        }
285
286        impl Sampleable<$kind> for Gev {
287            fn draw<R: Rng>(&self, rng: &mut R) -> $kind {
288                let uni = rand_distr::Open01;
289                let u: f64 = rng.sample(uni);
290                let lnu = -u.ln();
291                if self.shape == 0.0 {
292                    self.scale.mul_add(-lnu.ln(), self.loc) as $kind
293                } else {
294                    (self.loc
295                        + self.scale * (lnu.powf(-self.shape) - 1.0)
296                            / self.shape) as $kind
297                }
298            }
299        }
300
301        impl ContinuousDistr<$kind> for Gev {}
302
303        #[allow(clippy::cmp_owned)]
304        impl Support<$kind> for Gev {
305            fn supports(&self, x: &$kind) -> bool {
306                if self.shape > 0.0 {
307                    x.is_finite()
308                        && f64::from(*x) >= self.loc - self.scale / self.shape
309                } else if self.shape == 0.0 {
310                    x.is_finite()
311                } else {
312                    x.is_finite()
313                        && f64::from(*x) <= self.loc - self.scale / self.shape
314                }
315            }
316        }
317
318        impl Cdf<$kind> for Gev {
319            fn cdf(&self, x: &$kind) -> f64 {
320                (-t(self.loc, self.shape, self.scale, f64::from(*x))).exp()
321            }
322        }
323
324        impl Mean<$kind> for Gev {
325            fn mean(&self) -> Option<$kind> {
326                if self.shape == 0.0 {
327                    Some(self.scale.mul_add(consts::EULER_MASCERONI, self.loc)
328                        as $kind)
329                } else if self.shape >= 1.0 {
330                    Some(f64::INFINITY as $kind)
331                } else {
332                    let g1 = gammafn(1.0 - self.shape);
333                    Some(
334                        (self.loc + self.scale * (g1 - 1.0) / self.shape)
335                            as $kind,
336                    )
337                }
338            }
339        }
340
341        impl Mode<$kind> for Gev {
342            fn mode(&self) -> Option<$kind> {
343                if self.shape == 0.0 {
344                    Some(self.loc as $kind)
345                } else {
346                    Some(
347                        (self.loc
348                            + self.scale.mul_add(
349                                (1.0 + self.shape).powf(-self.shape),
350                                -1.0,
351                            ) / self.shape) as $kind,
352                    )
353                }
354            }
355        }
356
357        impl Median<$kind> for Gev {
358            fn median(&self) -> Option<$kind> {
359                if self.shape == 0.0 {
360                    Some(self.scale.mul_add(-consts::LN_LN_2, self.loc) as $kind)
361                } else {
362                    Some(
363                        (self.loc
364                            + self.scale * (LN_2.powf(-self.shape) - 1.0)
365                                / self.shape) as $kind,
366                    )
367                }
368            }
369        }
370    };
371}
372
373impl Variance<f64> for Gev {
374    fn variance(&self) -> Option<f64> {
375        if self.shape == 0.0 {
376            Some(self.scale * self.scale * PI * PI / 6.0)
377        } else if self.shape >= 0.5 {
378            Some(f64::INFINITY)
379        } else {
380            let g1 = gammafn(1.0 - self.shape);
381            let g2 = gammafn(2.0_f64.mul_add(-self.shape, 1.0));
382            Some(
383                self.scale * self.scale * g1.mul_add(-g1, g2)
384                    / (self.shape * self.shape),
385            )
386        }
387    }
388}
389
390impl Entropy for Gev {
391    fn entropy(&self) -> f64 {
392        consts::EULER_MASCERONI.mul_add(self.shape, self.scale.ln())
393            + consts::EULER_MASCERONI
394            + 1.0
395    }
396}
397
398impl_traits!(f32);
399impl_traits!(f64);
400
401impl std::error::Error for GevError {}
402
403impl fmt::Display for GevError {
404    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
405        match self {
406            Self::LocNotFinite { loc } => write!(f, "non-finite loc: {}", loc),
407            Self::ShapeNotFinite { shape } => {
408                write!(f, "non-finite shape: {}", shape)
409            }
410            Self::ScaleNotFinite { scale } => {
411                write!(f, "non-finite scale: {}", scale)
412            }
413            Self::ScaleTooLow { scale } => {
414                write!(f, "scale ({}) must be greater than zero", scale)
415            }
416        }
417    }
418}
419
420#[cfg(test)]
421mod tests {
422    use super::*;
423    use crate::misc::ks_test;
424    use crate::misc::linspace;
425    use crate::test_basic_impls;
426
427    const TOL: f64 = 1E-12;
428    const KS_PVAL: f64 = 0.2;
429    const N_TRIES: usize = 5;
430
431    test_basic_impls!(f64, Gev, Gev::new(0.0, 1.0, 2.0).unwrap());
432
433    #[test]
434    fn new() {
435        let gev = Gev::new(0.0, 1.0, 0.0).unwrap();
436        assert::close(gev.loc, 0.0, TOL);
437        assert::close(gev.scale, 1.0, TOL);
438        assert::close(gev.shape, 0.0, TOL);
439    }
440
441    #[test]
442    fn ln_pdf_0() {
443        let gev = Gev::new(0.0, 1.0, 0.0).unwrap();
444        let xs: Vec<f64> = linspace(-10.0, 10.0, 50);
445
446        let this_ln_pdf: Vec<f64> = xs.iter().map(|x| gev.ln_f(x)).collect();
447        let known_ln_pdf: Vec<f64> = vec![
448            -22_016.465_794_806_718,
449            -14_635.151_518_216_397,
450            -9_727.671_522_925_775,
451            -6_464.970_517_138_57,
452            -4_295.834_243_945_598,
453            -2_853.776_704_129_995_3,
454            -1_895.132_234_225_268_8,
455            -1_257.894_766_661_472_9,
456            -834.351_275_499_323_7,
457            -552.886_566_746_239_7,
458            -365.885_823_476_682_3,
459            -241.691_367_138_702_4,
460            -159.254_946_872_366_7,
461            -104.582_205_399_212_58,
462            -68.368_709_921_451_17,
463            -44.428_219_230_148_72,
464            -28.647_685_154_952_825,
465            -18.292_464_043_881_665,
466            -11.544_372_497_764_506,
467            -7.194_554_338_713_512_5,
468            -4.439_276_973_259_749,
469            -2.744_162_455_026_668,
470            -1.753_918_747_963_923_8,
471            -1.232_322_722_474_304_5,
472            -1.022_316_630_860_929_3,
473            -1.019_477_438_200_524,
474            -1.154_377_367_879_148_5,
475            -1.380_855_951_863_127_6,
476            -1.668_222_465_013_204_5,
477            -1.996_071_555_094_084_2,
478            -2.350_836_309_041_406,
479            -2.723_496_489_028_663,
480            -3.108_054_806_648_338_4,
481            -3.500_523_842_839_567,
482            -3.898_252_481_016_556,
483            -4.299_478_072_447_337,
484            -4.703_028_684_305_955,
485            -5.108_125_133_239_749,
486            -5.514_249_363_363_929,
487            -5.921_056_934_696_743,
488            -6.328_318_839_317_411,
489            -6.735_882_816_656_426,
490            -7.143_647_633_180_262_5,
491            -7.551_545_981_717_12,
492            -7.959_533_111_726_17,
493            -8.367_579_269_901_015,
494            -8.775_664_674_151_324,
495            -9.183_776_171_952_376,
496            -9.591_905_018_580_851,
497            -10.000_045_399_929_762,
498        ];
499
500        assert::close(known_ln_pdf, this_ln_pdf, TOL);
501    }
502
503    #[test]
504    fn ln_pdf_one_half() {
505        let gev = Gev::new(0.0, 1.0, 0.5).unwrap();
506        let xs: Vec<f64> = linspace(-1.9, 10.0, 50);
507
508        let this_ln_pdf: Vec<f64> = xs.iter().map(|x| gev.ln_f(x)).collect();
509        let known_ln_pdf: Vec<f64> = vec![
510            -391.012_803_179_337_25,
511            -28.737_012_000_993_683,
512            -7.975_515_285_646_1,
513            -3.182_798_910_065_802_3,
514            -1.611_981_517_225_461_7,
515            -1.056_128_444_415_616_3,
516            -0.898_809_165_661_975_9,
517            -0.918_486_354_261_089,
518            -1.022_088_700_315_005_2,
519            -1.166_219_177_873_567_8,
520            -1.329_140_366_485_869_2,
521            -1.499_425_187_893_265,
522            -1.670_888_816_300_113_7,
523            -1.840_148_707_985_855_8,
524            -2.005_377_976_051_166,
525            -2.165_637_389_658_979,
526            -2.320_503_404_006_605_5,
527            -2.469_854_528_307_160_5,
528            -2.613_745_589_117_496_8,
529            -2.752_332_330_080_753_4,
530            -2.885_825_601_856_779_6,
531            -3.014_463_329_165_149,
532            -3.138_493_349_819_397,
533            -3.258_162_997_160_523,
534            -3.373_712_908_918_762,
535            -3.485_373_502_391_386_3,
536            -3.593_363_135_370_378,
537            -3.697_887_329_478_015_2,
538            -3.799_138_656_162_781_6,
539            -3.897_297_027_438_662,
540            -3.992_530_224_449_923,
541            -4.084_994_555_888_228,
542            -4.174_835_576_763_844_5,
543            -4.262_188_823_288_923,
544            -4.347_180_536_268_123,
545            -4.429_928_356_362_477,
546            -4.510_541_981_812_075,
547            -4.589_123_783_926_562,
548            -4.665_769_378_708_396,
549            -4.740_568_154_914_112,
550            -4.813_603_760_052_456,
551            -4.884_954_546_513_983,
552            -4.954_693_980_392_284_5,
553            -5.022_891_015_706_101,
554            -5.089_610_436_741_23,
555            -5.154_913_171_153_501,
556            -5.218_856_576_344_224_5,
557            -5.281_494_701_461_209,
558            -5.342_878_527_206_992,
559            -5.403_056_185_461_942,
560        ];
561
562        assert::close(known_ln_pdf, this_ln_pdf, TOL);
563    }
564
565    #[test]
566    fn ln_pdf_minus_one_half() {
567        let gev = Gev::new(0.0, 1.0, -0.5).unwrap();
568        let xs: Vec<f64> = linspace(-10.0, 1.9, 50);
569
570        let this_ln_pdf: Vec<f64> = xs.iter().map(|x| gev.ln_f(x)).collect();
571        let known_ln_pdf: Vec<f64> = vec![
572            -34.208_240_530_771_945,
573            -32.786_288_262_748_58,
574            -31.394_252_557_653_69,
575            -30.032_151_611_967_51,
576            -28.700_004_811_360_04,
577            -27.397_832_836_625_59,
578            -26.125_657_781_682_3,
579            -24.883_503_285_324_355,
580            -23.671_394_678_696_16,
581            -22.489_359_150_795_124,
582            -21.337_425_934_714_08,
583            -20.215_626_517_822_667,
584            -19.123_994_879_677_34,
585            -18.062_567_762_169_61,
586            -17.031_384_977_300_48,
587            -16.030_489_759_050_92,
588            -15.059_929_167_153_456,
589            -14.119_754_552_231_905,
590            -13.210_022_093_853_19,
591            -12.330_793_425_651_896,
592            -11.482_136_365_003_441,
593            -10.664_125_768_956_731,
594            -9.876_844_543_585_701,
595            -9.120_384_840_989_5,
596            -8.394_849_487_426_416,
597            -7.700_353_698_296_974,
598            -7.037_027_152_018_743,
599            -6.405_016_516_870_358,
600            -5.804_488_554_972_564,
601            -5.235_633_969_193_058,
602            -4.698_672_217_128_37,
603            -4.193_857_599_416_090_5,
604            -3.721_487_049_917_897,
605            -3.281_910_232_624_673,
606            -2.875_542_816_807_392_7,
607            -2.502_884_212_064_577_3,
608            -2.164_541_691_614_048_5,
609            -1.861_263_880_969_972_4,
610            -1.593_988_345_178_63,
611            -1.363_911_057_382_413_6,
612            -1.172_591_056_355_069_7,
613            -1.022_114_118_880_009_1,
614            -0.915_360_515_657_826_6,
615            -0.856_468_009_767_915_7,
616            -0.851_690_580_254_141_6,
617            -0.911_144_104_991_361_2,
618            -1.052_832_065_124_108_8,
619            -1.313_835_662_027_446,
620            -1.792_976_347_363_397_7,
621            -2.998_232_273_553_990_4,
622        ];
623
624        assert::close(known_ln_pdf, this_ln_pdf, TOL);
625    }
626
627    #[test]
628    fn cdf() {
629        let gev_a = Gev::new(0.0, 1.0, 0.0).unwrap();
630        let gev_b = Gev::new(0.0, 1.0, 0.5).unwrap();
631        let gev_c = Gev::new(0.0, 1.0, -0.5).unwrap();
632
633        assert::close(gev_a.cdf(&0.0), 0.367_879_441_171_442_33, TOL);
634        assert::close(gev_a.cdf(&2.0), 0.873_423_018_493_116_7, TOL);
635        assert::close(gev_a.cdf(&-2.0), 0.000_617_978_989_331_093_4, TOL);
636
637        assert::close(gev_b.cdf(&0.0), 0.367_879_441_171_442_33, TOL);
638        assert::close(gev_b.cdf(&2.0), 0.778_800_783_071_404_9, TOL);
639        assert::close(gev_b.cdf(&-2.0), 0.0, TOL);
640
641        assert::close(gev_c.cdf(&0.0), 0.367_879_441_171_442_33, TOL);
642        assert::close(gev_c.cdf(&2.0), 1.0, TOL);
643        assert::close(gev_c.cdf(&-2.0), 0.018_315_638_888_734_18, TOL);
644    }
645
646    #[test]
647    fn entropy() {
648        let gev_a = Gev::new(0.0, 1.0, 0.0).unwrap();
649        let gev_b = Gev::new(0.0, 1.0, 0.5).unwrap();
650        let gev_c = Gev::new(0.0, 1.0, -0.5).unwrap();
651
652        assert::close(gev_a.entropy(), 1.577_215_664_901_532_8, TOL);
653        assert::close(gev_b.entropy(), 1.865_823_497_352_299_4, TOL);
654        assert::close(gev_c.entropy(), 1.288_607_832_450_766_4, TOL);
655    }
656
657    #[test]
658    fn draw_0() {
659        let mut rng = rand::thread_rng();
660        let gev = Gev::new(0.0, 1.0, 0.0).unwrap();
661        let cdf = |x: f64| gev.cdf(&x);
662
663        let passes = (0..N_TRIES).fold(0, |acc, _| {
664            let xs: Vec<f64> = gev.sample(1000, &mut rng);
665            let (_, p) = ks_test(&xs, cdf);
666            if p > KS_PVAL {
667                acc + 1
668            } else {
669                acc
670            }
671        });
672
673        assert!(passes > 0);
674    }
675
676    #[test]
677    fn draw_one_half() {
678        let mut rng = rand::thread_rng();
679        let gev = Gev::new(0.0, 1.0, 0.5).unwrap();
680        let cdf = |x: f64| gev.cdf(&x);
681
682        let passes = (0..N_TRIES).fold(0, |acc, _| {
683            let xs: Vec<f64> = gev.sample(1000, &mut rng);
684            let (_, p) = ks_test(&xs, cdf);
685            if p > KS_PVAL {
686                acc + 1
687            } else {
688                acc
689            }
690        });
691
692        assert!(passes > 0);
693    }
694
695    #[test]
696    fn draw_negative_one_half() {
697        let mut rng = rand::thread_rng();
698        let gev = Gev::new(0.0, 1.0, -0.5).unwrap();
699        let cdf = |x: f64| gev.cdf(&x);
700
701        let passes = (0..N_TRIES).fold(0, |acc, _| {
702            let xs: Vec<f64> = gev.sample(1000, &mut rng);
703            let (_, p) = ks_test(&xs, cdf);
704            if p > KS_PVAL {
705                acc + 1
706            } else {
707                acc
708            }
709        });
710
711        assert!(passes > 0);
712    }
713
714    #[test]
715    fn mode() {
716        let m1: f64 = Gev::new(0.0, 1.0, 0.0).unwrap().mode().unwrap();
717        let m2: f64 = Gev::new(0.0, 1.0, 0.5).unwrap().mode().unwrap();
718        let m3: f64 = Gev::new(0.0, 1.0, -0.5).unwrap().mode().unwrap();
719
720        assert::close(m1, 0.0, TOL);
721        assert::close(m2, -0.367_006_838_144_547_93, TOL);
722        assert::close(m3, 0.585_786_437_626_905, TOL);
723    }
724
725    #[test]
726    fn median() {
727        let m1: f64 = Gev::new(0.0, 1.0, 0.0).unwrap().median().unwrap();
728        let m2: f64 = Gev::new(0.0, 1.0, 0.5).unwrap().median().unwrap();
729        let m3: f64 = Gev::new(0.0, 1.0, -0.5).unwrap().median().unwrap();
730
731        assert::close(m1, 0.366_512_920_581_664_35, TOL);
732        assert::close(m2, 0.402_244_817_572_899_6, TOL);
733        assert::close(m3, 0.334_890_777_684_604_5, TOL);
734    }
735
736    #[test]
737    fn variance() {
738        let m1: f64 = Gev::new(0.0, 1.0, 0.0).unwrap().variance().unwrap();
739        let m2: f64 = Gev::new(0.0, 1.0, 0.5).unwrap().variance().unwrap();
740        let m3: f64 = Gev::new(0.0, 1.0, -0.5).unwrap().variance().unwrap();
741
742        assert::close(m1, 1.644_934_066_848_226_4, TOL);
743        assert::close(m2, f64::INFINITY, TOL);
744        assert::close(m3, 0.858_407_346_410_206_8, TOL);
745    }
746}
rv/dist/gev.rs

rv/dist/
gev.rs
