Merge pull request #40 from mstoeckl/norm-invg-edge-case

Avoid panic on large alpha for NormalInverseGaussian

Author: dhardy

CHANGELOG.md

@@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixes
 - Fix `Geometric::new` for small `p > 0` where `1 - p` rounds to 1 (#36)
 - Fix panic in `FisherF::new` on almost zero parameters (#39)
+- Fix panic in `NormalInverseGaussian::new` with very large `alpha`; this is a Value-breaking change (#40)
 
 ## [0.5.2]
 

src/normal_inverse_gaussian.rs

@@ -1,4 +1,4 @@
-use crate::{Distribution, InverseGaussian, StandardNormal, StandardUniform};
+use crate::{Distribution, InverseGaussian, InverseGaussianError, StandardNormal, StandardUniform};
 use core::fmt;
 use num_traits::Float;
 use rand::Rng;
@@ -8,6 +8,8 @@ use rand::Rng;
 pub enum Error {
     /// `alpha <= 0` or `nan`.
     AlphaNegativeOrNull,
+    /// `alpha` is `inf` (or, if subnormals are disabled, too close to the maximum finite float value)
+    AlphaInfinite,
     /// `|beta| >= alpha` or `nan`.
     AbsoluteBetaNotLessThanAlpha,
 }
@@ -18,6 +20,9 @@ impl fmt::Display for Error {
             Error::AlphaNegativeOrNull => {
                 "alpha <= 0 or is NaN in normal inverse Gaussian distribution"
             }
+            Error::AlphaInfinite => {
+                "alpha is +infinity (or too close to the maximum finite value, if subnormal numbers are not supported) in normal inverse Gaussian distribution"
+            }
             Error::AbsoluteBetaNotLessThanAlpha => {
                 "|beta| >= alpha or is NaN in normal inverse Gaussian distribution"
             }
@@ -68,6 +73,9 @@ where
 {
     /// Construct a new `NormalInverseGaussian` distribution with the given alpha (tail heaviness) and
     /// beta (asymmetry) parameters.
+    ///
+    /// Note: If subnormal numbers are not supported or are disabled, this sampler may panic or produce
+    /// incorrect output if alpha is close to the maximum finite float value.
     pub fn new(alpha: F, beta: F) -> Result<NormalInverseGaussian<F>, Error> {
         if !(alpha > F::zero()) {
             return Err(Error::AlphaNegativeOrNull);
@@ -76,12 +84,16 @@ where
         if !(beta.abs() < alpha) {
             return Err(Error::AbsoluteBetaNotLessThanAlpha);
         }
-
-        let gamma = (alpha * alpha - beta * beta).sqrt();
-
+        // Note: this calculation method for gamma = sqrt(alpha * alpha - beta * beta)
+        // avoids overflow if alpha is large, ensuring gamma <= alpha, which implies
+        // (assuming IEEE754 with subnormals) mu = 1.0 / gamma >= 1 / F::max_value() > 0.
+        let r = beta / alpha;
+        let gamma = alpha * (F::one() - r * r).sqrt();
         let mu = F::one() / gamma;
-
-        let inverse_gaussian = InverseGaussian::new(mu, F::one()).unwrap();
+        let inverse_gaussian = InverseGaussian::new(mu, F::one()).map_err(|x| match x {
+            InverseGaussianError::MeanNegativeOrNull => Error::AlphaInfinite,
+            InverseGaussianError::ShapeNegativeOrNull => unreachable!(),
+        })?;
 
         Ok(Self {
             beta,
@@ -124,6 +136,7 @@ mod tests {
         assert!(NormalInverseGaussian::new(-1.0, 1.0).is_err());
         assert!(NormalInverseGaussian::new(-1.0, -1.0).is_err());
         assert!(NormalInverseGaussian::new(1.0, 2.0).is_err());
+        assert!(NormalInverseGaussian::new(f64::INFINITY, 2.0).is_err());
         assert!(NormalInverseGaussian::new(2.0, 1.0).is_ok());
     }