From 94a5b8a13d0f27543038829d5f249917768242f8 Mon Sep 17 00:00:00 2001 From: Soumya Date: Mon, 9 Sep 2024 23:44:06 +0200 Subject: [PATCH 01/10] Add support for Gumbel Distribution - Added basic struct and implementation of Gumbel Distribution and GumbelError --- src/distribution/gumbel.rs | 56 ++++++++++++++++++++++++++++++++++++++ src/distribution/mod.rs | 1 + 2 files changed, 57 insertions(+) create mode 100644 src/distribution/gumbel.rs diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs new file mode 100644 index 00000000..d21386a3 --- /dev/null +++ b/src/distribution/gumbel.rs @@ -0,0 +1,56 @@ +/// https://en.wikipedia.org/wiki/Gumbel_distribution +#[derive(Copy, Clone, PartialEq, Debug)] +pub struct Gumbel { + location: f64, + scale: f64, +} + +#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)] +#[non_exhaustive] +pub enum GumbelError { + /// The location is invalid (NAN) + LocationInvalid, + + /// The scale is NAN, zero or less than zero + ScaleInvalid, +} + +impl std::fmt::Display for GumbelError { + #[cfg_attr(coverage_nightly, coverage(off))] + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + GumbelError::LocationInvalid => write!(f, "Location is NAN"), + GumbelError::ScaleInvalid => write!(f, "Scale is NAN, zero or less than zero"), + } + } +} + +impl std::error::Error for GumbelError {} + +impl Gumbel { + pub fn new(location: f64, scale: f64) -> Result { + if location.is_nan() { + return Err(GumbelError::LocationInvalid); + } + + if scale.is_nan() && scale <= 0.0 { + return Err(GumbelError::ScaleInvalid); + } + + Ok(Self { location, scale }) + } + + pub fn location(&self) -> f64 { + self.location + } + + pub fn scale(&self) -> f64 { + self.scale + } +} + +impl std::fmt::Display for Gumbel { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "Gumbel({:?}, {:?})", self.location(), self.scale()) + } +} diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index e8c9574c..df8f5e68 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -58,6 +58,7 @@ mod exponential; mod fisher_snedecor; mod gamma; mod geometric; +mod gumbel; mod hypergeometric; #[macro_use] mod internal; From d865190eca1f6971d2da4db03a8bd94f62f9bb45 Mon Sep 17 00:00:00 2001 From: Soumya Date: Tue, 10 Sep 2024 20:21:30 +0200 Subject: [PATCH 02/10] Implemented Continuous CDF trait for Gumbel - Added cdf, inverse_cdf and sf for Gumbel - Added Min and Max traits for Gumbel - Added documentation for the same --- src/distribution/gumbel.rs | 87 +++++++++++++++++++++++++++++++++++++- 1 file changed, 86 insertions(+), 1 deletion(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index d21386a3..ca5e168d 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -1,3 +1,9 @@ +use core::f64; + +use crate::statistics::{Max, Min}; + +use super::ContinuousCDF; + /// https://en.wikipedia.org/wiki/Gumbel_distribution #[derive(Copy, Clone, PartialEq, Debug)] pub struct Gumbel { @@ -51,6 +57,85 @@ impl Gumbel { impl std::fmt::Display for Gumbel { fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { - write!(f, "Gumbel({:?}, {:?})", self.location(), self.scale()) + write!(f, "Gumbel({:?}, {:?})", self.location, self.scale) + } +} + +impl ::rand::distributions::Distribution for Gumbel { + fn sample(&self, r: &mut R) -> f64 { + // Check: Quantile formula: mu - beta*ln(-ln(p)) + self.location - self.scale * (((-(r.gen::())).ln()).ln()) + } +} + +impl ContinuousCDF for Gumbel { + /// Calculates the cumulative distribution function for the + /// gumbel distribution at `x` + /// + /// # Formula + /// + /// ```text + /// e^(-e^(-(x - μ) / β)) + /// ``` + /// + /// where `μ` is the location and `β` is the scale + fn cdf(&self, x: f64) -> f64 { + (-(-(x - self.location) / self.scale).exp()).exp() + } + + /// Calculates the inverse cumulative distribution function for the + /// gumbel distribution at `x` + /// + /// # Formula + /// + /// ```text + /// μ - β ln(-ln(p)) + /// ``` + /// + /// where `μ` is the location and `β` is the scale + fn inverse_cdf(&self, p: f64) -> f64 { + self.location - self.scale * ((-(p.ln())).ln()) + } + + /// Calculates the survival function for the + /// gumbel distribution at `x` + /// + /// # Formula + /// + /// ```text + /// 1 - e^(-e^(-(x - μ) / β)) + /// ``` + /// + /// where `μ` is the location and `β` is the scale + fn sf(&self, x: f64) -> f64 { + 1.0 - (-(-(x - self.location) / self.scale).exp()).exp() + } +} + +impl Min for Gumbel { + /// Returns the minimum value in the domain of the gumbel + /// distribution representable by a double precision float + /// + /// # Formula + /// + /// ```text + /// NEG_INF + /// ``` + fn min(&self) -> f64 { + f64::NEG_INFINITY + } +} + +impl Max for Gumbel { + /// Returns the maximum value in the domain of the gumbel + /// distribution representable by a double precision float + /// + /// # Formula + /// + /// ```text + /// f64::INFINITY + /// ``` + fn max(&self) -> f64 { + f64::INFINITY } } From 9e62e34dc06f313e35f3bd9fdbd6ae45bf9954c2 Mon Sep 17 00:00:00 2001 From: Soumya Date: Tue, 10 Sep 2024 20:53:39 +0200 Subject: [PATCH 03/10] Implement Distribution trait for Gumbel - Added mean, skewness, variance and entropy for Gumbel distribution - Added documentation for the same --- src/distribution/gumbel.rs | 66 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 64 insertions(+), 2 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index ca5e168d..d49eb597 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -1,6 +1,10 @@ use core::f64; +use std::f64::consts::PI; -use crate::statistics::{Max, Min}; +use crate::{ + consts::EULER_MASCHERONI, + statistics::{Distribution, Max, Min}, +}; use super::ContinuousCDF; @@ -64,7 +68,7 @@ impl std::fmt::Display for Gumbel { impl ::rand::distributions::Distribution for Gumbel { fn sample(&self, r: &mut R) -> f64 { // Check: Quantile formula: mu - beta*ln(-ln(p)) - self.location - self.scale * (((-(r.gen::())).ln()).ln()) + self.location - self.scale * ((-(r.gen::())).ln()).ln() } } @@ -139,3 +143,61 @@ impl Max for Gumbel { f64::INFINITY } } + +impl Distribution for Gumbel { + /// Returns the entropy of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// ln(β) + γ + 1 + /// ``` + /// + /// where `β` is the scale + /// and `γ` is the Euler-Mascheroni constant (approx 0.57721) + fn entropy(&self) -> Option { + Some(1.0 + EULER_MASCHERONI + (self.scale).ln()) + } + + /// Returns the mean of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// μ + γβ + /// ``` + /// + /// where `μ` is the location, `β` is the scale + /// and `γ` is the Euler-Mascheroni constant (approx 0.57721) + fn mean(&self) -> Option { + Some(self.location + (EULER_MASCHERONI * self.scale)) + } + + /// Returns the skewness of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// 12 * sqrt(6) * ζ(3) / π^3 ≈ 1.13955 + /// ``` + /// ζ(3) is the Riemann zeta function evaluated at 3 (approx 1.20206) + /// π is the constant PI (approx 3.14159) + /// This approximately evaluates to 1.13955\ + /// Ref: https://www.statsref.com/HTML/gumbel_extreme_value_distribut.html + fn skewness(&self) -> Option { + Some(1.13955) + } + + /// Returns the variance of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// (π^2 / 6) * β^2 + /// ``` + /// + /// where `β` is the scale and `π` is the constant PI (approx 3.14159) + fn variance(&self) -> Option { + Some(((PI * PI) / 6.0) * self.scale * self.scale) + } +} From fbc0def57d4fbd107e774395275e62da3080f364 Mon Sep 17 00:00:00 2001 From: Soumya Date: Tue, 10 Sep 2024 21:04:28 +0200 Subject: [PATCH 04/10] Implemented median and mode for Gumbel - Added median and mode traits - Added std deviation for the Distribution trait - Minor cleanups --- src/distribution/gumbel.rs | 48 +++++++++++++++++++++++++++++++++++--- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index d49eb597..af193200 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -3,7 +3,7 @@ use std::f64::consts::PI; use crate::{ consts::EULER_MASCHERONI, - statistics::{Distribution, Max, Min}, + statistics::{Distribution, Max, Median, Min, Mode}, }; use super::ContinuousCDF; @@ -182,8 +182,7 @@ impl Distribution for Gumbel { /// ``` /// ζ(3) is the Riemann zeta function evaluated at 3 (approx 1.20206) /// π is the constant PI (approx 3.14159) - /// This approximately evaluates to 1.13955\ - /// Ref: https://www.statsref.com/HTML/gumbel_extreme_value_distribut.html + /// This approximately evaluates to 1.13955 fn skewness(&self) -> Option { Some(1.13955) } @@ -200,4 +199,47 @@ impl Distribution for Gumbel { fn variance(&self) -> Option { Some(((PI * PI) / 6.0) * self.scale * self.scale) } + + /// Returns the standard deviation of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// β * π / sqrt(6) + /// ``` + /// + /// where `β` is the scale and `π` is the constant PI (approx 3.14159) + fn std_dev(&self) -> Option { + Some(self.scale * PI / 6.0_f64.sqrt()) + } +} + +impl Median for Gumbel { + /// Returns the median of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// μ - β ln(ln(2)) + /// ``` + /// + /// where `μ` is the location and `β` is the scale parameter + fn median(&self) -> f64 { + self.location - self.scale * (((2.0_f64).ln()).ln()) + } +} + +impl Mode for Gumbel { + /// Returns the mode of the gumbel distribution + /// + /// # Formula + /// + /// ```text + /// μ + /// ``` + /// + /// where `μ` is the location + fn mode(&self) -> f64 { + self.location + } } From 7ad585ad382ca86b85d487f662095500073e1141 Mon Sep 17 00:00:00 2001 From: Soumya Date: Tue, 10 Sep 2024 21:39:09 +0200 Subject: [PATCH 05/10] Implemented Continuous trait for Gumbel - Added pdf and ln_pdf for Gumbel --- src/distribution/gumbel.rs | 39 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index af193200..7f73353e 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -3,10 +3,10 @@ use std::f64::consts::PI; use crate::{ consts::EULER_MASCHERONI, - statistics::{Distribution, Max, Median, Min, Mode}, + statistics::{Distribution, Max, MeanN, Median, Min, Mode}, }; -use super::ContinuousCDF; +use super::{Continuous, ContinuousCDF}; /// https://en.wikipedia.org/wiki/Gumbel_distribution #[derive(Copy, Clone, PartialEq, Debug)] @@ -243,3 +243,38 @@ impl Mode for Gumbel { self.location } } + +impl Continuous for Gumbel { + /// Calculates the probability density function for the gumbel + /// distribution at `x` + /// + /// # Formula + /// + /// ```text + /// (1/β) * exp(-(x - μ)/β) * exp(-exp(-(x - μ)/β)) + /// ``` + /// + /// where `μ` is the location, `β` is the scale + fn pdf(&self, x: f64) -> f64 { + (1.0_f64 / self.scale) + * (-(x - self.location) / (self.scale)).exp() + * (-((-(x - self.location) / self.scale).exp())).exp() + } + + /// Calculates the log probability density function for the gumbel + /// distribution at `x` + /// + /// # Formula + /// + /// ```text + /// ln((1/β) * exp(-(x - μ)/β) * exp(-exp(-(x - μ)/β))) + /// ``` + /// + /// where `μ` is the location, `β` is the scale + fn ln_pdf(&self, x: f64) -> f64 { + ((1.0_f64 / self.scale) + * (-(x - self.location) / (self.scale)).exp() + * (-((-(x - self.location) / self.scale).exp())).exp()) + .ln() + } +} From 443c9feda3c364df0ea215a76092225a55340b82 Mon Sep 17 00:00:00 2001 From: Soumya Date: Tue, 10 Sep 2024 23:38:01 +0200 Subject: [PATCH 06/10] Commence testing for Gumbel - Added tests for mean, median, mode, min-max, std dev, variance, entropy - Added valid and invalid creation tests --- src/distribution/gumbel.rs | 115 ++++++++++++++++++++++++++++++++++++- 1 file changed, 113 insertions(+), 2 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index 7f73353e..4ff43362 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -3,7 +3,7 @@ use std::f64::consts::PI; use crate::{ consts::EULER_MASCHERONI, - statistics::{Distribution, Max, MeanN, Median, Min, Mode}, + statistics::{Distribution, Max, Median, Min, Mode}, }; use super::{Continuous, ContinuousCDF}; @@ -43,7 +43,7 @@ impl Gumbel { return Err(GumbelError::LocationInvalid); } - if scale.is_nan() && scale <= 0.0 { + if scale.is_nan() || scale <= 0.0 { return Err(GumbelError::ScaleInvalid); } @@ -278,3 +278,114 @@ impl Continuous for Gumbel { .ln() } } + +#[rustfmt::skip] +#[cfg(test)] +mod tests { + use super::*; + use crate::distribution::internal::*; + use crate::testing_boiler; + + testing_boiler!(location: f64, scale: f64; Gumbel; GumbelError); + + #[test] + fn test_create() { + create_ok(0.0, 0.1); + create_ok(0.0, 1.0); + create_ok(0.0, 10.0); + create_ok(10.0, 11.0); + create_ok(-5.0, 100.0); + create_ok(0.0, f64::INFINITY); + } + + #[test] + fn test_bad_create() { + let invalid = [ + (f64::NAN, 1.0, GumbelError::LocationInvalid), + (1.0, f64::NAN, GumbelError::ScaleInvalid), + (f64::NAN, f64::NAN, GumbelError::LocationInvalid), + (1.0, 0.0, GumbelError::ScaleInvalid), + (0.0, f64::NEG_INFINITY, GumbelError::ScaleInvalid) + ]; + + for (location, scale, err) in invalid { + test_create_err(location, scale, err); + } + } + + #[test] + fn test_min_max() { + let min = |x: Gumbel| x.min(); + let max = |x:Gumbel| x.max(); + + test_exact(0.0, 1.0, f64::NEG_INFINITY, min); + test_exact(0.0, 1.0, f64::INFINITY, max); + } + + #[test] + fn test_entropy() { + let entropy = |x: Gumbel| x.entropy().unwrap(); + test_exact(0.0, 2.0, 2.270362845461478, entropy); + test_exact(0.1, 4.0, 2.9635100260214235, entropy); + test_exact(1.0, 10.0, 3.8798007578955787, entropy); + test_exact(10.0, 11.0, 3.9751109376999034, entropy); + } + + #[test] + fn test_mean() { + let mean = |x: Gumbel| x.mean().unwrap(); + test_exact(0.0, 2.0, 1.1544313298030658, mean); + test_exact(0.1, 4.0, 2.4088626596061316, mean); + test_exact(1.0, 10.0, 6.772156649015328, mean); + test_exact(10.0, 11.0, 16.34937231391686, mean); + test_exact(10.0, f64::INFINITY, f64::INFINITY, mean); + } + + #[test] + fn test_skewness() { + let skewness = |x: Gumbel| x.skewness().unwrap(); + test_exact(0.0, 2.0, 1.13955, skewness); + test_exact(0.1, 4.0, 1.13955, skewness); + test_exact(1.0, 10.0, 1.13955, skewness); + test_exact(10.0, 11.0, 1.13955, skewness); + test_exact(10.0, f64::INFINITY, 1.13955, skewness); + } + + #[test] + fn test_variance() { + let variance = |x: Gumbel| x.variance().unwrap(); + test_exact(0.0, 2.0, 6.579736267392906, variance); + test_exact(0.1, 4.0, 26.318945069571624, variance); + test_exact(1.0, 10.0, 164.49340668482265, variance); + test_exact(10.0, 11.0, 199.03702208863538, variance); + } + + #[test] + fn test_std_dev() { + let std_dev = |x: Gumbel| x.std_dev().unwrap(); + test_exact(0.0, 2.0, 2.565099660323728, std_dev); + test_exact(0.1, 4.0, 5.130199320647456, std_dev); + test_exact(1.0, 10.0, 12.82549830161864, std_dev); + test_exact(10.0, 11.0, 14.108048131780505, std_dev); + } + + #[test] + fn test_median() { + let median = |x: Gumbel| x.median(); + test_exact(0.0, 2.0, 0.7330258411633287, median); + test_exact(0.1, 4.0, 1.5660516823266574, median); + test_exact(1.0, 10.0, 4.665129205816644, median); + test_exact(10.0, 11.0, 14.031642126398307, median); + test_exact(10.0, f64::INFINITY, f64::INFINITY, median); + } + + #[test] + fn test_mode() { + let mode = |x: Gumbel| x.mode(); + test_exact(0.0, 2.0, 0.0, mode); + test_exact(0.1, 4.0, 0.1, mode); + test_exact(1.0, 10.0, 1.0, mode); + test_exact(10.0, 11.0, 10.0, mode); + test_exact(10.0, f64::INFINITY, 10.0, mode); + } +} From a0b766848b0c17107de34e5481657d04ab8b91d8 Mon Sep 17 00:00:00 2001 From: Soumya Date: Wed, 11 Sep 2024 11:23:48 +0200 Subject: [PATCH 07/10] Added tests for Gumbel distrbution - Added test for CDF and inverse_CDF --- src/distribution/gumbel.rs | 55 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 53 insertions(+), 2 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index 4ff43362..d923bff7 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -93,12 +93,20 @@ impl ContinuousCDF for Gumbel { /// # Formula /// /// ```text - /// μ - β ln(-ln(p)) + /// μ - β ln(-ln(p)) where 0 < p < 1 + /// -INF where p <= 0 + /// INF otherwise /// ``` /// /// where `μ` is the location and `β` is the scale fn inverse_cdf(&self, p: f64) -> f64 { - self.location - self.scale * ((-(p.ln())).ln()) + if p <= 0.0 { + f64::NEG_INFINITY + } else if p >= 1.0 { + f64::INFINITY + } else { + self.location - self.scale * ((-(p.ln())).ln()) + } } /// Calculates the survival function for the @@ -388,4 +396,47 @@ mod tests { test_exact(10.0, 11.0, 10.0, mode); test_exact(10.0, f64::INFINITY, 10.0, mode); } + + #[test] + fn test_cdf() { + let cdf = |a: f64| move |x: Gumbel| x.cdf(a); + test_exact(0.0, 0.1, 0.0, cdf(-5.0)); + test_exact(0.0, 0.1, 0.0, cdf(-1.0)); + test_exact(0.0, 0.1, 0.36787944117144233, cdf(0.0)); + test_exact(0.0, 0.1, 0.9999546011007987, cdf(1.0)); + test_absolute(0.0, 0.1, 0.99999999999999999, 1e-12, cdf(5.0)); + test_absolute(0.0, 1.0, 0.06598803584531253, 1e-12, cdf(-1.0)); + test_exact(0.0, 1.0, 0.36787944117144233, cdf(0.0)); + test_absolute(0.0, 10.0, 0.192295645547964928, 1e-12, cdf(-5.0)); + test_absolute(0.0, 10.0, 0.3311542771529088, 1e-12, cdf(-1.0)); + test_exact(0.0, 10.0, 0.36787944117144233, cdf(0.0)); + test_absolute(0.0, 10.0, 0.4046076616641318, 1e-12, cdf(1.0)); + test_absolute(0.0, 10.0, 0.545239211892605, 1e-12, cdf(5.0)); + test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(-5.0)); + test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(-1.0)); + test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(0.0)); + test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(1.0)); + test_exact(-2.0, f64::INFINITY, 0.36787944117144233, cdf(5.0)); + test_exact(f64::INFINITY, 1.0, 0.0, cdf(-5.0)); + test_exact(f64::INFINITY, 1.0, 0.0, cdf(-1.0)); + test_exact(f64::INFINITY, 1.0, 0.0, cdf(0.0)); + test_exact(f64::INFINITY, 1.0, 0.0, cdf(1.0)); + test_exact(f64::INFINITY, 1.0, 0.0, cdf(5.0)); + } + + #[test] + fn test_inverse_cdf() { + let inv_cdf = |a: f64| move |x: Gumbel| x.inverse_cdf(a); + test_exact(0.0, 0.1, f64::NEG_INFINITY, inv_cdf(-5.0)); + test_exact(0.0, 0.1, f64::NEG_INFINITY, inv_cdf(-1.0)); + test_exact(0.0, 0.1, f64::NEG_INFINITY, inv_cdf(0.0)); + test_exact(0.0, 0.1, f64::INFINITY, inv_cdf(1.0)); + test_exact(0.0, 0.1, f64::INFINITY, inv_cdf(5.0)); + test_absolute(0.0, 1.0, -0.8340324452479557, 1e-12, inv_cdf(0.1)); + test_absolute(0.0, 10.0, 3.6651292058166436, 1e-12, inv_cdf(0.5)); + test_absolute(0.0, 10.0, 22.503673273124456, 1e-12, inv_cdf(0.9)); + test_exact(2.0, f64::INFINITY, f64::NEG_INFINITY, inv_cdf(0.1)); + test_exact(-2.0, f64::INFINITY, f64::INFINITY, inv_cdf(0.5)); + test_exact(f64::INFINITY, 1.0, f64::INFINITY, inv_cdf(0.1)); + } } From 90527672b2b662f402653198f85d97e5fb0218d2 Mon Sep 17 00:00:00 2001 From: Soumya Date: Wed, 11 Sep 2024 12:32:47 +0200 Subject: [PATCH 08/10] Added tests for Gumbel - Added tests for pdf, inverse_pdf and sf --- src/distribution/gumbel.rs | 73 +++++++++++++++++++++++++++++++++++--- 1 file changed, 69 insertions(+), 4 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index d923bff7..d2b640bd 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -8,7 +8,6 @@ use crate::{ use super::{Continuous, ContinuousCDF}; -/// https://en.wikipedia.org/wiki/Gumbel_distribution #[derive(Copy, Clone, PartialEq, Debug)] pub struct Gumbel { location: f64, @@ -67,7 +66,6 @@ impl std::fmt::Display for Gumbel { impl ::rand::distributions::Distribution for Gumbel { fn sample(&self, r: &mut R) -> f64 { - // Check: Quantile formula: mu - beta*ln(-ln(p)) self.location - self.scale * ((-(r.gen::())).ln()).ln() } } @@ -189,7 +187,8 @@ impl Distribution for Gumbel { /// 12 * sqrt(6) * ζ(3) / π^3 ≈ 1.13955 /// ``` /// ζ(3) is the Riemann zeta function evaluated at 3 (approx 1.20206) - /// π is the constant PI (approx 3.14159) + /// and π is the constant PI (approx 3.14159) + /// /// This approximately evaluates to 1.13955 fn skewness(&self) -> Option { Some(1.13955) @@ -291,7 +290,6 @@ impl Continuous for Gumbel { #[cfg(test)] mod tests { use super::*; - use crate::distribution::internal::*; use crate::testing_boiler; testing_boiler!(location: f64, scale: f64; Gumbel; GumbelError); @@ -439,4 +437,71 @@ mod tests { test_exact(-2.0, f64::INFINITY, f64::INFINITY, inv_cdf(0.5)); test_exact(f64::INFINITY, 1.0, f64::INFINITY, inv_cdf(0.1)); } + + #[test] + fn test_sf() { + let sf = |a: f64| move |x: Gumbel| x.sf(a); + test_exact(0.0, 0.1, 1.0, sf(-5.0)); + test_exact(0.0, 0.1, 1.0, sf(-1.0)); + test_absolute(0.0, 0.1, 0.632120558828557678, 1e-12, sf(0.0)); + test_absolute(0.0, 0.1, 0.000045398899201269, 1e-12, sf(1.0)); + test_absolute(0.0, 1.0, 0.934011964154687462, 1e-12, sf(-1.0)); + test_absolute(0.0, 1.0, 0.632120558828557678, 1e-12, sf(0.0)); + test_absolute(0.0, 1.0, 0.3077993724446536, 1e-12, sf(1.0)); + test_absolute(0.0, 10.0, 0.66884572284709110, 1e-12, sf(-1.0)); + test_absolute(0.0, 10.0, 0.632120558828557678, 1e-12, sf(0.0)); + test_absolute(0.0, 10.0, 0.595392338335868174, 1e-12, sf(1.0)); + test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(-5.0)); + test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(-1.0)); + test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(0.0)); + test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(1.0)); + test_exact(-2.0, f64::INFINITY, 0.6321205588285576784, sf(5.0)); + test_exact(f64::INFINITY, 1.0, 1.0, sf(-5.0)); + test_exact(f64::INFINITY, 1.0, 1.0, sf(-1.0)); + test_exact(f64::INFINITY, 1.0, 1.0, sf(0.0)); + test_exact(f64::INFINITY, 1.0, 1.0, sf(1.0)); + test_exact(f64::INFINITY, 1.0, 1.0, sf(5.0)); + } + + #[test] + fn test_pdf() { + let pdf = |a: f64| move |x: Gumbel| x.pdf(a); + test_exact(0.0, 0.1, 0.0, pdf(-5.0)); + test_exact(0.0, 0.1, 3.678794411714423215, pdf(0.0)); + test_absolute(0.0, 0.1, 0.0004539786865564, 1e-12, pdf(1.0)); + test_absolute(0.0, 1.0, 0.1793740787340171, 1e-12, pdf(-1.0)); + test_exact(0.0, 1.0, 0.36787944117144233, pdf(0.0)); + test_absolute(0.0, 1.0, 0.25464638004358249, 1e-12, pdf(1.0)); + test_absolute(0.0, 10.0, 0.031704192107794217, 1e-12, pdf(-5.0)); + test_absolute(0.0, 10.0, 0.0365982076505757, 1e-12, pdf(-1.0)); + test_exact(0.0, 10.0, 0.036787944117144233, pdf(0.0)); + test_absolute(0.0, 10.0, 0.03661041518977401428, 1e-12, pdf(1.0)); + test_absolute(0.0, 10.0, 0.033070429889041, 1e-12, pdf(5.0)); + test_exact(-2.0, f64::INFINITY, 0.0, pdf(-5.0)); + test_exact(-2.0, f64::INFINITY, 0.0, pdf(-1.0)); + test_exact(-2.0, f64::INFINITY, 0.0, pdf(0.0)); + test_exact(-2.0, f64::INFINITY, 0.0, pdf(1.0)); + test_exact(-2.0, f64::INFINITY, 0.0, pdf(5.0)); + } + + #[test] + fn test_ln_pdf() { + let ln_pdf = |a: f64| move |x: Gumbel| x.ln_pdf(a); + test_exact(0.0, 0.1, 0.0_f64.ln(), ln_pdf(-5.0)); + test_exact(0.0, 0.1, 3.678794411714423215_f64.ln(), ln_pdf(0.0)); + test_absolute(0.0, 0.1, 0.0004539786865564_f64.ln(), 1e-12, ln_pdf(1.0)); + test_absolute(0.0, 1.0, 0.1793740787340171_f64.ln(), 1e-12, ln_pdf(-1.0)); + test_exact(0.0, 1.0, 0.36787944117144233_f64.ln(), ln_pdf(0.0)); + test_absolute(0.0, 1.0, 0.25464638004358249_f64.ln(), 1e-12, ln_pdf(1.0)); + test_absolute(0.0, 10.0, 0.031704192107794217_f64.ln(), 1e-12, ln_pdf(-5.0)); + test_absolute(0.0, 10.0, 0.0365982076505757_f64.ln(), 1e-12, ln_pdf(-1.0)); + test_exact(0.0, 10.0, 0.036787944117144233_f64.ln(), ln_pdf(0.0)); + test_absolute(0.0, 10.0, 0.03661041518977401428_f64.ln(), 1e-12, ln_pdf(1.0)); + test_absolute(0.0, 10.0, 0.033070429889041_f64.ln(), 1e-12, ln_pdf(5.0)); + test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(-5.0)); + test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(-1.0)); + test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(0.0)); + test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(1.0)); + test_exact(-2.0, f64::INFINITY, 0.0_f64.ln(), ln_pdf(5.0)); + } } From dbbec67b72ef21494c28d48f166e01499fca645f Mon Sep 17 00:00:00 2001 From: Soumya Date: Wed, 11 Sep 2024 16:24:18 +0200 Subject: [PATCH 09/10] feat: Add support for Gumbel distribution - Improved documentation for Gumbel - Added missing examples - Added Gumbel to mod.rs --- src/distribution/gumbel.rs | 94 ++++++++++++++++++++++++++++---------- src/distribution/mod.rs | 1 + 2 files changed, 72 insertions(+), 23 deletions(-) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index d2b640bd..273ea874 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -1,19 +1,29 @@ -use core::f64; -use std::f64::consts::PI; - -use crate::{ - consts::EULER_MASCHERONI, - statistics::{Distribution, Max, Median, Min, Mode}, -}; - use super::{Continuous, ContinuousCDF}; - +use crate::consts::EULER_MASCHERONI; +use crate::statistics::*; +use std::f64::{self, consts::PI}; + +/// Implements the [Gumbel](https://en.wikipedia.org/wiki/Gumbel_distribution) +/// distribution, also known as the type-I generalized extreme value distribution. +/// +/// # Examples +/// +/// ``` +/// use statrs::distribution::{Gumbel, Continuous}; +/// use statrs::{consts::EULER_MASCHERONI, statistics::Distribution}; +/// +/// let n = Gumbel::new(0.0, 1.0).unwrap(); +/// assert_eq!(n.location(), 0.0); +/// assert_eq!(n.skewness().unwrap(), 1.13955); +/// assert_eq!(n.pdf(0.0), 0.36787944117144233); +/// ``` #[derive(Copy, Clone, PartialEq, Debug)] pub struct Gumbel { location: f64, scale: f64, } +/// Represents the errors that can occur when creating a [`Gumbel`] #[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)] #[non_exhaustive] pub enum GumbelError { @@ -37,6 +47,24 @@ impl std::fmt::Display for GumbelError { impl std::error::Error for GumbelError {} impl Gumbel { + /// Constructs a new Gumbel distribution with the given + /// location and scale. + /// + /// # Errors + /// + /// Returns an error if location or scale are `NaN` or `scale <= 0.0` + /// + /// # Examples + /// + /// ``` + /// use statrs::distribution::Gumbel; + /// + /// let mut result = Gumbel::new(0.0, 1.0); + /// assert!(result.is_ok()); + /// + /// result = Gumbel::new(0.0, -1.0); + /// assert!(result.is_err()); + /// ``` pub fn new(location: f64, scale: f64) -> Result { if location.is_nan() { return Err(GumbelError::LocationInvalid); @@ -49,10 +77,30 @@ impl Gumbel { Ok(Self { location, scale }) } + /// Returns the location of the Gumbel distribution + /// + /// # Examples + /// + /// ``` + /// use statrs::distribution::Gumbel; + /// + /// let n = Gumbel::new(0.0, 1.0).unwrap(); + /// assert_eq!(n.location(), 0.0); + /// ``` pub fn location(&self) -> f64 { self.location } + /// Returns the scale of the Gumbel distribution + /// + /// # Examples + /// + /// ``` + /// use statrs::distribution::Gumbel; + /// + /// let n = Gumbel::new(0.0, 1.0).unwrap(); + /// assert_eq!(n.scale(), 1.0); + /// ``` pub fn scale(&self) -> f64 { self.scale } @@ -72,7 +120,7 @@ impl ::rand::distributions::Distribution for Gumbel { impl ContinuousCDF for Gumbel { /// Calculates the cumulative distribution function for the - /// gumbel distribution at `x` + /// Gumbel distribution at `x` /// /// # Formula /// @@ -86,7 +134,7 @@ impl ContinuousCDF for Gumbel { } /// Calculates the inverse cumulative distribution function for the - /// gumbel distribution at `x` + /// Gumbel distribution at `x` /// /// # Formula /// @@ -108,7 +156,7 @@ impl ContinuousCDF for Gumbel { } /// Calculates the survival function for the - /// gumbel distribution at `x` + /// Gumbel distribution at `x` /// /// # Formula /// @@ -123,7 +171,7 @@ impl ContinuousCDF for Gumbel { } impl Min for Gumbel { - /// Returns the minimum value in the domain of the gumbel + /// Returns the minimum value in the domain of the Gumbel /// distribution representable by a double precision float /// /// # Formula @@ -137,7 +185,7 @@ impl Min for Gumbel { } impl Max for Gumbel { - /// Returns the maximum value in the domain of the gumbel + /// Returns the maximum value in the domain of the Gumbel /// distribution representable by a double precision float /// /// # Formula @@ -151,7 +199,7 @@ impl Max for Gumbel { } impl Distribution for Gumbel { - /// Returns the entropy of the gumbel distribution + /// Returns the entropy of the Gumbel distribution /// /// # Formula /// @@ -165,7 +213,7 @@ impl Distribution for Gumbel { Some(1.0 + EULER_MASCHERONI + (self.scale).ln()) } - /// Returns the mean of the gumbel distribution + /// Returns the mean of the Gumbel distribution /// /// # Formula /// @@ -179,7 +227,7 @@ impl Distribution for Gumbel { Some(self.location + (EULER_MASCHERONI * self.scale)) } - /// Returns the skewness of the gumbel distribution + /// Returns the skewness of the Gumbel distribution /// /// # Formula /// @@ -194,7 +242,7 @@ impl Distribution for Gumbel { Some(1.13955) } - /// Returns the variance of the gumbel distribution + /// Returns the variance of the Gumbel distribution /// /// # Formula /// @@ -207,7 +255,7 @@ impl Distribution for Gumbel { Some(((PI * PI) / 6.0) * self.scale * self.scale) } - /// Returns the standard deviation of the gumbel distribution + /// Returns the standard deviation of the Gumbel distribution /// /// # Formula /// @@ -222,7 +270,7 @@ impl Distribution for Gumbel { } impl Median for Gumbel { - /// Returns the median of the gumbel distribution + /// Returns the median of the Gumbel distribution /// /// # Formula /// @@ -237,7 +285,7 @@ impl Median for Gumbel { } impl Mode for Gumbel { - /// Returns the mode of the gumbel distribution + /// Returns the mode of the Gumbel distribution /// /// # Formula /// @@ -252,7 +300,7 @@ impl Mode for Gumbel { } impl Continuous for Gumbel { - /// Calculates the probability density function for the gumbel + /// Calculates the probability density function for the Gumbel /// distribution at `x` /// /// # Formula @@ -268,7 +316,7 @@ impl Continuous for Gumbel { * (-((-(x - self.location) / self.scale).exp())).exp() } - /// Calculates the log probability density function for the gumbel + /// Calculates the log probability density function for the Gumbel /// distribution at `x` /// /// # Formula diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index df8f5e68..b0b50f67 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -22,6 +22,7 @@ pub use self::exponential::{Exp, ExpError}; pub use self::fisher_snedecor::{FisherSnedecor, FisherSnedecorError}; pub use self::gamma::{Gamma, GammaError}; pub use self::geometric::{Geometric, GeometricError}; +pub use self::gumbel::{Gumbel, GumbelError}; pub use self::hypergeometric::{Hypergeometric, HypergeometricError}; pub use self::inverse_gamma::{InverseGamma, InverseGammaError}; pub use self::laplace::{Laplace, LaplaceError}; From 3dd22f4a0b4247e245dafe10658b8e5fd04c39a4 Mon Sep 17 00:00:00 2001 From: Soumya Date: Fri, 13 Sep 2024 00:38:38 +0200 Subject: [PATCH 10/10] feat: Add support for Gumbel distribution - Fix added: rand feature added to distribution to fix CI error --- src/distribution/gumbel.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/distribution/gumbel.rs b/src/distribution/gumbel.rs index 273ea874..a71f0dcb 100644 --- a/src/distribution/gumbel.rs +++ b/src/distribution/gumbel.rs @@ -112,6 +112,7 @@ impl std::fmt::Display for Gumbel { } } +#[cfg(feature = "rand")] impl ::rand::distributions::Distribution for Gumbel { fn sample(&self, r: &mut R) -> f64 { self.location - self.scale * ((-(r.gen::())).ln()).ln()