fix: Stabilize inverse beta lower tails

src/distribution/beta.rs

@@ -183,7 +183,7 @@ impl ContinuousCDF<f64, f64> for Beta {
     ///
     /// # Panics
     ///
-    /// If x is not in `[0, 1]`.
+    /// If x is not in `[0, 1]` or the numerical method does not converge.
     ///
     /// # Formula
     ///
@@ -206,7 +206,7 @@ impl ContinuousCDF<f64, f64> for Beta {
     ///
     /// # Returns an error instead of a panic
     ///
-    /// If x is not in `[0, 1]`.
+    /// If x is not in `[0, 1]` or the numerical method does not converge.
     ///
     /// # Formula
     ///
@@ -220,7 +220,8 @@ impl ContinuousCDF<f64, f64> for Beta {
         if !(0.0..=1.0).contains(&x) {
             Err(InverseCdfError::ArgumentOutOfRange)
         } else {
-            Ok(beta::inv_beta_reg(self.shape_a, self.shape_b, x))
+            beta::try_inv_beta_reg(self.shape_a, self.shape_b, x)
+                .map_err(|_| InverseCdfError::ConvergenceFailed)
         }
     }
 }
@@ -714,6 +715,29 @@ mod tests {
         }
     }
 
+    #[test]
+    fn test_try_inverse_cdf_extreme_shape_ratio() {
+        let distribution = Beta::new(1e20, 10.0).unwrap();
+        assert_eq!(distribution.try_inverse_cdf(0.3), Ok(1.0));
+    }
+
+    #[test]
+    fn test_try_inverse_cdf_deep_tail() {
+        let distribution = Beta::new(10.0, 1e10).unwrap();
+        let actual = distribution.try_inverse_cdf(0.3).unwrap();
+        let expected = 8.132928235539348e-10_f64;
+        assert!(actual.to_bits().abs_diff(expected.to_bits()) <= 8);
+    }
+
+    #[test]
+    fn test_try_inverse_cdf_reports_unrepresentable_quantile() {
+        let distribution = Beta::new(1e308, 1e308).unwrap();
+        assert_eq!(
+            distribution.try_inverse_cdf(0.3),
+            Err(InverseCdfError::ConvergenceFailed)
+        );
+    }
+
     #[test]
     fn test_cdf_input_lt_0() {
         let cdf = |arg: f64| move |x: Beta| x.cdf(arg);

src/distribution/mod.rs

@@ -102,13 +102,17 @@ mod ziggurat_tables;
 pub enum InverseCdfError {
     /// The argument `p` is outside the closed interval `[0, 1]`.
     ArgumentOutOfRange,
+
+    /// The numerical method did not converge.
+    ConvergenceFailed,
 }
 
 impl core::fmt::Display for InverseCdfError {
     #[cfg_attr(coverage_nightly, coverage(off))]
     fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result {
         match self {
             InverseCdfError::ArgumentOutOfRange => write!(f, "argument is outside [0, 1]"),
+            InverseCdfError::ConvergenceFailed => write!(f, "computation did not converge"),
         }
     }
 }