diff --git a/src/distribution/chi.rs b/src/distribution/chi.rs index 208047db..c087b42d 100644 --- a/src/distribution/chi.rs +++ b/src/distribution/chi.rs @@ -119,12 +119,17 @@ impl ContinuousCDF for Chi { /// where `k` is the degrees of freedom and `P` is /// the regularized lower incomplete Gamma function fn cdf(&self, x: f64) -> f64 { - if x == f64::INFINITY { - 1.0 - } else if x <= 0.0 { + let k = self.freedom() as f64; + let u = x * x / 2.0; + if x <= 0.0 { 0.0 + } else if u == f64::INFINITY { + 1.0 + } else if u < f64::MIN_POSITIVE { + // x * x / 2 has lost precision; use the leading term (x^2 / 2)^(k/2) / Gamma(k/2 + 1) + (k * x.ln() - k / 2.0 * f64_consts::LN_2 - gamma::ln_gamma(k / 2.0 + 1.0)).exp() } else { - gamma::gamma_lr(self.freedom() as f64 / 2.0, x * x / 2.0) + gamma::gamma_lr(k / 2.0, u) } } @@ -140,12 +145,13 @@ impl ContinuousCDF for Chi { /// where `k` is the degrees of freedom and `P` is /// the regularized upper incomplete Gamma function fn sf(&self, x: f64) -> f64 { - if x == f64::INFINITY { - 0.0 - } else if x <= 0.0 { + let u = x * x / 2.0; + if x <= 0.0 { 1.0 + } else if u == f64::INFINITY { + 0.0 } else { - gamma::gamma_ur(self.freedom() as f64 / 2.0, x * x / 2.0) + gamma::gamma_ur(self.freedom() as f64 / 2.0, u) } } @@ -503,6 +509,14 @@ mod tests { test_exact(2, 1.0, cdf(f64::INFINITY)); test_exact(2, 0.0, cdf(0.0)); test_exact(2, 1.0, cdf(f64::INFINITY)); + // x * x / 2 underflows or overflows; for k = 1, cdf(x) = erf(x / sqrt(2)) ~ x * sqrt(2 / pi) + for x in [1e-160, 1e-200] { + let expected = x * (2.0 / f64_consts::PI).sqrt(); + crate::prec::assert_relative_eq!(create_ok(1).cdf(x), expected, epsilon = 0.0, max_relative = 1e-12); + } + test_exact(3, 0.0, cdf(1e-200)); + test_exact(1, 1.0, cdf(1e160)); + test_exact(3, 1.0, cdf(f64::MAX)); } #[test] @@ -520,18 +534,23 @@ mod tests { test_exact(2, 0.0, sf(f64::INFINITY)); test_exact(2, 1.0, sf(0.0)); test_exact(2, 0.0, sf(f64::INFINITY)); + test_exact(1, 1.0, sf(1e-200)); + test_exact(1, 0.0, sf(1e160)); + test_exact(3, 0.0, sf(f64::MAX)); } #[test] fn test_neg_cdf() { let cdf = |arg: f64| move |x: Chi| x.cdf(arg); test_exact(1, 0.0, cdf(-1.0)); + test_exact(1, 0.0, cdf(f64::NEG_INFINITY)); } #[test] fn test_neg_sf() { let sf = |arg: f64| move |x: Chi| x.sf(arg); test_exact(1, 1.0, sf(-1.0)); + test_exact(1, 1.0, sf(f64::NEG_INFINITY)); } #[test] diff --git a/src/distribution/gamma.rs b/src/distribution/gamma.rs index 6df26e38..6b92cc20 100644 --- a/src/distribution/gamma.rs +++ b/src/distribution/gamma.rs @@ -152,7 +152,7 @@ impl ContinuousCDF for Gamma { 1.0 } else if self.rate.is_infinite() { 0.0 - } else if x.is_infinite() { + } else if x.is_infinite() || (x * self.rate).is_infinite() { 1.0 } else { gamma::gamma_lr(self.shape, x * self.rate) @@ -177,7 +177,7 @@ impl ContinuousCDF for Gamma { 0.0 } else if self.rate.is_infinite() { 1.0 - } else if x.is_infinite() { + } else if x.is_infinite() || (x * self.rate).is_infinite() { 0.0 } else { gamma::gamma_ur(self.shape, x * self.rate) @@ -745,6 +745,14 @@ mod tests { test_relative(1.0, 0.1, 0.0, |x| x.cdf(0.0)); } + #[test] + fn test_cdf_sf_when_x_times_rate_overflows_or_underflows() { + test_exact(2.5, 1.5, 1.0, |g| g.cdf(f64::MAX)); + test_exact(2.5, 1.5, 0.0, |g| g.sf(f64::MAX)); + test_exact(1.0, 0.5, 0.0, |g| g.cdf(5e-324)); + test_exact(1.0, 0.5, 1.0, |g| g.sf(5e-324)); + } + #[test] fn test_infinite_rate_special_point_is_exact() { let dist = create_ok(10.0, f64::INFINITY); diff --git a/src/function/gamma.rs b/src/function/gamma.rs index 51d5fed1..eaf793ad 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -14,7 +14,7 @@ pub enum GammaFuncError { /// `a` is infinite, zero or less than zero. AInvalid, - /// `x` is infinite, zero or less than zero. + /// `x` is infinite or less than zero. XInvalid, } @@ -22,7 +22,7 @@ impl core::fmt::Display for GammaFuncError { fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result { match self { GammaFuncError::AInvalid => write!(f, "a is infinite, zero or less than zero"), - GammaFuncError::XInvalid => write!(f, "x is infinite, zero or less than zero"), + GammaFuncError::XInvalid => write!(f, "x is infinite or less than zero"), } } } @@ -113,7 +113,7 @@ pub fn gamma(x: f64) -> f64 { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_ui(a: f64, x: f64) -> f64 { checked_gamma_ui(a, x).unwrap() } @@ -125,7 +125,7 @@ pub fn gamma_ui(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_ui(a: f64, x: f64) -> Result { checked_gamma_ur(a, x).map(|x| x * gamma(a)) } @@ -138,7 +138,7 @@ pub fn checked_gamma_ui(a: f64, x: f64) -> Result { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_li(a: f64, x: f64) -> f64 { checked_gamma_li(a, x).unwrap() } @@ -151,9 +151,10 @@ pub fn gamma_li(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_li(a: f64, x: f64) -> Result { - checked_gamma_lr(a, x).map(|x| x * gamma(a)) + // gamma(a) can overflow, and 0 * inf is NaN + checked_gamma_lr(a, x).map(|p| if x == 0.0 { 0.0 } else { p * gamma(a) }) } /// Computes the upper incomplete regularized gamma function @@ -167,7 +168,7 @@ pub fn checked_gamma_li(a: f64, x: f64) -> Result { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_ur(a: f64, x: f64) -> f64 { checked_gamma_ur(a, x).unwrap() } @@ -183,7 +184,7 @@ pub fn gamma_ur(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_ur(a: f64, x: f64) -> Result { if a.is_nan() || x.is_nan() { return Ok(f64::NAN); @@ -191,7 +192,7 @@ pub fn checked_gamma_ur(a: f64, x: f64) -> Result { if a <= 0.0 || a == f64::INFINITY { return Err(GammaFuncError::AInvalid); } - if x <= 0.0 || x == f64::INFINITY { + if x < 0.0 || x == f64::INFINITY { return Err(GammaFuncError::XInvalid); } @@ -261,7 +262,7 @@ pub fn checked_gamma_ur(a: f64, x: f64) -> Result { /// /// # Panics /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn gamma_lr(a: f64, x: f64) -> f64 { checked_gamma_lr(a, x).unwrap() } @@ -277,7 +278,7 @@ pub fn gamma_lr(a: f64, x: f64) -> f64 { /// /// # Errors /// -/// if `a` or `x` are not in `(0, +inf)` +/// if `a` is not in `(0, +inf)` or `x` is not in `[0, +inf)` pub fn checked_gamma_lr(a: f64, x: f64) -> Result { if a.is_nan() || x.is_nan() { return Ok(f64::NAN); @@ -285,7 +286,7 @@ pub fn checked_gamma_lr(a: f64, x: f64) -> Result { if a <= 0.0 || a == f64::INFINITY { return Err(GammaFuncError::AInvalid); } - if x <= 0.0 || x == f64::INFINITY { + if x < 0.0 || x == f64::INFINITY { return Err(GammaFuncError::XInvalid); } @@ -915,6 +916,16 @@ mod tests { assert!(super::checked_gamma_lr(1.0, f64::INFINITY).is_err()); } + #[test] + fn test_incomplete_gamma_at_x_zero() { + assert_eq!(gamma_lr(2.5, 0.0), 0.0); + assert_eq!(gamma_ur(2.5, 0.0), 1.0); + assert_eq!(super::checked_gamma_lr(0.5, 0.0), Ok(0.0)); + assert_eq!(super::checked_gamma_ur(0.5, 0.0), Ok(1.0)); + assert_eq!(gamma_li(172.0, 0.0), 0.0); + assert_eq!(super::checked_gamma_li(172.0, 0.0), Ok(0.0)); + } + #[test] fn test_gamma_li() { assert!(gamma_li(f64::NAN, f64::NAN).is_nan());