From db0eb53b441846065e25ae8987443fc926ff1ce3 Mon Sep 17 00:00:00 2001 From: Przemyslaw Olszewski Date: Wed, 26 Aug 2026 15:16:09 +0200 Subject: [PATCH 1/5] fix: Stabilize inverse beta tails --- src/function/beta.rs | 282 ++++++++++-------------------- src/function/beta/inverse.rs | 207 ++++++++++++++++++++++ src/function/beta/large_params.rs | 2 +- 3 files changed, 305 insertions(+), 186 deletions(-) create mode 100644 src/function/beta/inverse.rs diff --git a/src/function/beta.rs b/src/function/beta.rs index 033a0afe..eeb5fc06 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -3,6 +3,7 @@ //! //! This module sets the default precision more tightly than crate defaults for `DEFAULT_EPS` +mod inverse; mod large_params; mod temme; @@ -336,192 +337,9 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { } /// Computes the inverse of the regularized incomplete beta function -// This code is based on the implementation in the ["special"][1] crate, -// which in turn is based on a [C implementation][2] by John Burkardt. The -// original algorithm was published in Applied Statistics and is known as -// [Algorithm AS 64][3] and [Algorithm AS 109][4]. -// -// [1]: https://docs.rs/special/0.8.1/ -// [2]: http://people.sc.fsu.edu/~jburkardt/c_src/asa109/asa109.html -// [3]: http://www.jstor.org/stable/2346798 -// [4]: http://www.jstor.org/stable/2346887 -// -// > Copyright 2014–2019 The special Developers -// > -// > Permission is hereby granted, free of charge, to any person obtaining a copy of -// > this software and associated documentation files (the "Software"), to deal in -// > the Software without restriction, including without limitation the rights to -// > use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of -// > the Software, and to permit persons to whom the Software is furnished to do so, -// > subject to the following conditions: -// > -// > The above copyright notice and this permission notice shall be included in all -// > copies or substantial portions of the Software. -// > -// > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS -// > FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR -// > COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER -// > IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN -// > CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -pub fn inv_beta_reg(mut a: f64, mut b: f64, mut x: f64) -> f64 { - // Algorithm AS 64 - // http://www.jstor.org/stable/2346798 - // - // An approximation x₀ to x if found from (cf. Scheffé and Tukey, 1944) - // - // 1 + x₀ 4p + 2q - 2 - // ------ = ----------- - // 1 - x₀ χ²(α) - // - // where χ²(α) is the upper α point of the χ² distribution with 2q - // degrees of freedom and is obtained from Wilson and Hilferty's - // approximation (cf. Wilson and Hilferty, 1931) - // - // χ²(α) = 2q (1 - 1 / (9q) + y(α) sqrt(1 / (9q)))^3, - // - // y(α) being Hastings' approximation (cf. Hastings, 1955) for the upper - // α point of the standard normal distribution. If χ²(α) < 0, then - // - // x₀ = 1 - ((1 - α)q B(p, q))^(1 / q). - // - // Again if (4p + 2q - 2) / χ²(α) does not exceed 1, x₀ is obtained from - // - // x₀ = (αp B(p, q))^(1 / p). - // - // The final solution is obtained by the Newton–Raphson method from the - // relation - // - // f(x[i - 1]) - // x[i] = x[i - 1] - ------------ - // f'(x[i - 1]) - // - // where - // - // f(x) = I(x, p, q) - α. - let ln_beta = ln_beta(a, b); - - // Remark AS R83 - // http://www.jstor.org/stable/2347779 - const SAE: i32 = -30; - const FPU: f64 = 1e-30; // 10^SAE - - debug_assert!((0.0..=1.0).contains(&x) && a > 0.0 && b > 0.0); - - if x == 0.0 { - return 0.0; - } - if x == 1.0 { - return 1.0; - } - - let mut p; - let mut q; - - let flip = 0.5 < x; - if flip { - p = a; - a = b; - b = p; - x = 1.0 - x; - } - - p = (-(x * x).ln()).sqrt(); - q = p - (2.30753 + 0.27061 * p) / (1.0 + (0.99229 + 0.04481 * p) * p); - - if 1.0 < a && 1.0 < b { - // Remark AS R19 and Algorithm AS 109 - // http://www.jstor.org/stable/2346887 - // - // For a and b > 1, the approximation given by Carter (1947), which - // improves the Fisher–Cochran formula, is generally better. For - // other values of a and b en empirical investigation has shown that - // the approximation given in AS 64 is adequate. - let r = (q * q - 3.0) / 6.0; - let s = 1.0 / (2.0 * a - 1.0); - let t = 1.0 / (2.0 * b - 1.0); - let h = 2.0 / (s + t); - let w = q * (h + r).sqrt() / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h)); - p = a / (a + b * (2.0 * w).exp()); - } else { - let mut t = 1.0 / (9.0 * b); - t = 2.0 * b * (1.0 - t + q * t.sqrt()).powf(3.0); - if t <= 0.0 { - p = 1.0 - ((((1.0 - x) * b).ln() + ln_beta) / b).exp(); - } else { - t = 2.0 * (2.0 * a + b - 1.0) / t; - if t <= 1.0 { - p = (((x * a).ln() + ln_beta) / a).exp(); - } else { - p = 1.0 - 2.0 / (t + 1.0); - } - } - } - - p = p.clamp(0.0001, 0.9999); - - // Remark AS R83 - // http://www.jstor.org/stable/2347779 - let e = (-5.0 / a / a - 1.0 / x.powf(0.2) - 13.0) as i32; - let acu = if e > SAE { f64::powi(10.0, e) } else { FPU }; - - let mut pnext; - let mut qprev = 0.0; - let mut sq = 1.0; - let mut prev = 1.0; - - 'outer: loop { - // Remark AS R19 and Algorithm AS 109 - // http://www.jstor.org/stable/2346887 - q = beta_reg(a, b, p); - q = (q - x) * (ln_beta + (1.0 - a) * p.ln() + (1.0 - b) * (1.0 - p).ln()).exp(); - - // Remark AS R83 - // http://www.jstor.org/stable/2347779 - if q * qprev <= 0.0 { - prev = if sq > FPU { sq } else { FPU }; - } - - // Remark AS R19 and Algorithm AS 109 - // http://www.jstor.org/stable/2346887 - let mut g = 1.0; - loop { - loop { - let adj = g * q; - sq = adj * adj; - - if sq < prev { - pnext = p - adj; - if (0.0..=1.0).contains(&pnext) { - break; - } - } - g /= 3.0; - } - - if prev <= acu || q * q <= acu { - p = pnext; - break 'outer; - } - - if pnext != 0.0 && pnext != 1.0 { - break; - } - - g /= 3.0; - } - - if pnext == p { - break; - } - - p = pnext; - qprev = q; - } - - if flip { 1.0 - p } else { p } +pub fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { + inverse::inv_beta_reg(a, b, probability) } - #[cfg(test)] mod tests { use super::*; @@ -1081,6 +899,100 @@ mod tests { assert!(checked_beta_reg(1.0, 1.0, 2.0).is_err()); } + #[test] + fn test_inv_beta_reg_extreme_probability_does_not_panic() { + let actual = inv_beta_reg(200.0, 2.0, 1e-165).to_bits(); + assert!(actual.abs_diff(0x3fc2_aa4f_7f31_6421) <= 1); + } + + #[test] + fn test_inv_beta_reg_extreme_probability_terminates() { + let actual = inv_beta_reg(200.0, 2.0, 1e-60).to_bits(); + assert!(actual.abs_diff(0x3fdf_5753_caf6_9652) <= 1); + } + + #[test] + fn test_inv_beta_reg_small_shape_lower_tail_is_monotone() { + let cases = [ + (1e-33, 0.0), + (1e-32, f64::from_bits(2)), + (1e-31, 1.215703604971242e-313), + (1e-30, 1.2157036049544172e-303), + (1e-20, 1.2157036049544e-203), + (1e-10, 1.2157036049543856e-103), + (1e-4, 1.2157036049543764e-43), + (1e-2, 1.215703604954373e-23), + ]; + let mut previous = 0.0; + + for (probability, expected) in cases { + let actual = inv_beta_reg(0.1, 500.0, probability); + if expected == 0.0 { + assert_eq!(actual, expected); + continue; + } + if expected < f64::MIN_POSITIVE { + assert!(actual.to_bits().abs_diff(expected.to_bits()) <= 1); + previous = actual; + continue; + } + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-12, + "probability={probability:?}, actual={actual:?}, expected={expected:?}, relative_error={relative_error:?}" + ); + assert!(actual > previous); + previous = actual; + } + } + + #[test] + fn test_inv_beta_reg_regular_shape_lower_tail() { + let cases = [ + (1e-300, 7.053456158585983e-153), + (1e-40, 7.053456158585983e-23), + (1e-30, 7.053456158585999e-18), + (1e-20, 7.053456158916007e-13), + ]; + + for (probability, expected) in cases { + let actual = inv_beta_reg(2.0, 200.0, probability); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-12, + "probability={probability:?}, actual={actual:?}, expected={expected:?}, relative_error={relative_error:?}" + ); + } + } + + #[test] + fn test_inv_beta_reg_is_monotone_and_round_trips() { + let probabilities = [1e-12, 1e-6, 0.01, 0.25, 0.5, 0.75, 0.99, 1.0 - 1e-6]; + for (a, b) in [ + (0.1, 0.1), + (0.1, 500.0), + (2.0, 200.0), + (200.0, 2.0), + (100.0, 200.0), + (1e8, 2e8), + ] { + let mut previous = 0.0; + for probability in probabilities { + let quantile = inv_beta_reg(a, b, probability); + if quantile > 0.0 && quantile < 1.0 { + let recovered = beta_reg(a, b, quantile); + let tolerance = 5e-11 * probability.max(1.0 - probability); + assert!( + (recovered - probability).abs() <= tolerance, + "a={a:?}, b={b:?}, probability={probability:?}, quantile={quantile:?}, recovered={recovered:?}" + ); + } + assert!(quantile >= previous); + previous = quantile; + } + } + } + #[test] fn test_error_is_sync_send() { fn assert_sync_send() {} diff --git a/src/function/beta/inverse.rs b/src/function/beta/inverse.rs new file mode 100644 index 00000000..74934e42 --- /dev/null +++ b/src/function/beta/inverse.rs @@ -0,0 +1,207 @@ +use super::{ + BetaFuncError, beta_continued_fraction, checked_beta_reg, + large_params::{self, LogPrefactor}, + ln_beta, +}; +use crate::function::gamma; +#[cfg(not(feature = "std"))] +use num_traits::Float as _; + +const MAX_ITERATIONS: usize = 128; + +#[derive(Clone, Copy)] +struct Point { + log_x: f64, + x: f64, + log_cdf: f64, +} + +fn log_regularized_beta(a: f64, b: f64, x: f64, log_beta: f64) -> Result { + if x == 0.0 { + return Ok(f64::NEG_INFINITY); + } + if x == 1.0 { + return Ok(0.0); + } + + let symmetry_split = (a + 1.0) / (a + b + 2.0); + if x >= symmetry_split { + return checked_beta_reg(a, b, x).map(f64::ln); + } + + let fraction = beta_continued_fraction(a, b, x)?; + let log_prefactor = match large_params::log_prefactor(a, b, x) { + Some(LogPrefactor::Value(parts)) => parts.0 + parts.1, + Some(LogPrefactor::Underflow) => f64::NEG_INFINITY, + None => a * x.ln() + b * (-x).ln_1p() - log_beta, + }; + Ok(log_prefactor + (fraction / a).ln()) +} + +fn closer_point(lower: Point, upper: Point, log_probability: f64) -> f64 { + let lower_error = (lower.log_cdf - log_probability).abs(); + let upper_error = (upper.log_cdf - log_probability).abs(); + if lower_error <= upper_error { + lower.x + } else { + upper.x + } +} + +fn closest_representable( + a: f64, + b: f64, + log_beta: f64, + lower: Point, + upper: Point, + log_probability: f64, +) -> Result { + let mut best = lower.x; + let mut best_error = (lower.log_cdf - log_probability).abs(); + for bits in lower.x.to_bits() + 1..=upper.x.to_bits() { + let x = f64::from_bits(bits); + let error = (log_regularized_beta(a, b, x, log_beta)? - log_probability).abs(); + if error < best_error { + best = x; + best_error = error; + } + } + Ok(best) +} + +fn inverse_log_beta(a: f64, b: f64) -> f64 { + if b == 1.0 { + return -a.ln(); + } + if b == 2.0 { + return -a.ln() - (a + 1.0).ln(); + } + if a == 1.0 { + return -b.ln(); + } + if a == 2.0 { + return -b.ln() - (b + 1.0).ln(); + } + if a.min(b) >= 10.0 { + let scale = a.max(b); + let scaled_a = a / scale; + let scaled_b = b / scale; + let x = scaled_a / (scaled_a + scaled_b); + if let Some(LogPrefactor::Value(prefactor)) = large_params::log_prefactor(a, b, x) { + return a * x.ln() + b * (-x).ln_1p() - prefactor.0 - prefactor.1; + } + } + let smaller = a.min(b); + let larger = a.max(b); + if larger >= 10.0 { + let ratio = smaller / larger; + let gamma_difference = -smaller * larger.ln() - (larger + smaller - 0.5) * ratio.ln_1p() + + smaller + + large_params::stirling_correction(larger) + - large_params::stirling_correction(larger + smaller); + return gamma::ln_gamma(smaller) + gamma_difference; + } + ln_beta(a, b) +} + +fn solve_lower_tail(a: f64, b: f64, probability: f64) -> f64 { + let log_probability = probability.ln(); + let log_beta = inverse_log_beta(a, b); + let smallest = f64::from_bits(1); + let smallest_log = smallest.ln(); + let smallest_log_cdf = log_regularized_beta(a, b, smallest, log_beta) + .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + + if smallest_log_cdf > log_probability { + return 0.0; + } + if smallest_log_cdf == log_probability { + return smallest; + } + + let mut lower = Point { + log_x: smallest_log, + x: smallest, + log_cdf: smallest_log_cdf, + }; + let mut upper = Point { + log_x: 0.0, + x: 1.0, + log_cdf: 0.0, + }; + let estimate = (log_probability + a.ln() + log_beta) / a; + let mut log_x = estimate.clamp(lower.log_x, upper.log_x); + if log_x == lower.log_x || log_x == upper.log_x { + log_x = lower.log_x + 0.5 * (upper.log_x - lower.log_x); + } + + for _ in 0..MAX_ITERATIONS { + let x = log_x.exp(); + let log_cdf = log_regularized_beta(a, b, x, log_beta) + .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + let current = Point { log_x, x, log_cdf }; + + if log_cdf < log_probability { + lower = current; + } else if log_cdf > log_probability { + upper = current; + } else { + return x; + } + + if upper.x.to_bits() - lower.x.to_bits() <= 8 { + return closest_representable(a, b, log_beta, lower, upper, log_probability) + .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + } + let log_tolerance = 32.0 * f64::EPSILON * log_probability.abs().max(1.0); + if (lower.log_cdf - log_probability).abs() <= log_tolerance + && (upper.log_cdf - log_probability).abs() <= log_tolerance + { + return closer_point(lower, upper, log_probability); + } + + let log_density = (a - 1.0) * log_x + (b - 1.0) * (-x).ln_1p() - log_beta; + let slope = (log_x + log_density - log_cdf).exp(); + let newton = log_x - (log_cdf - log_probability) / slope; + let midpoint = lower.log_x + 0.5 * (upper.log_x - lower.log_x); + let mut next = if newton.is_finite() && newton > lower.log_x && newton < upper.log_x { + newton + } else { + midpoint + }; + + if next.exp() == x { + next = midpoint; + } + log_x = next; + } + + panic!( + "inv_beta_reg did not converge for a={a}, b={b}, probability={probability}, lower={:?}, upper={:?}, lower_log_cdf={:?}, upper_log_cdf={:?}", + lower.x, upper.x, lower.log_cdf, upper.log_cdf + ) +} + +pub(super) fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { + assert!(a.is_finite() && a > 0.0, "a must be finite and positive"); + assert!(b.is_finite() && b > 0.0, "b must be finite and positive"); + assert!( + probability.is_finite() && (0.0..=1.0).contains(&probability), + "probability must be finite and in [0, 1]" + ); + + if probability == 0.0 { + return 0.0; + } + if probability == 1.0 { + return 1.0; + } + if probability == 0.5 && a == b { + return 0.5; + } + if probability <= 0.5 { + return solve_lower_tail(a, b, probability); + } + + 1.0 - solve_lower_tail(b, a, 1.0 - probability) +} diff --git a/src/function/beta/large_params.rs b/src/function/beta/large_params.rs index 71d94313..4b96a025 100644 --- a/src/function/beta/large_params.rs +++ b/src/function/beta/large_params.rs @@ -28,7 +28,7 @@ pub(super) fn log_ratio(a: f64, b: f64, x: f64) -> (f64, (f64, f64)) { (residual.0 + residual.1, value) } -fn stirling_correction(value: f64) -> f64 { +pub(super) fn stirling_correction(value: f64) -> f64 { let inverse = 1.0 / value; let inverse_squared = inverse * inverse; let mut series: f64 = 7.0 / 1_092.0; From 5046dc9ce699c7ae72f7a8c605ddb660a2b09560 Mon Sep 17 00:00:00 2001 From: Przemyslaw Olszewski Date: Wed, 26 Aug 2026 16:51:13 +0200 Subject: [PATCH 2/5] fix: Propagate inverse beta failures --- src/distribution/beta.rs | 20 ++++++++++++++-- src/distribution/mod.rs | 4 ++++ src/function/beta.rs | 15 +++++++++++- src/function/beta/inverse.rs | 44 +++++++++++++++++------------------- 4 files changed, 57 insertions(+), 26 deletions(-) diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs index a5b4d36b..7cd83727 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -206,7 +206,7 @@ impl ContinuousCDF 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 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,21 @@ 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_reports_convergence_failure() { + let distribution = Beta::new(10.0, 1e10).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); diff --git a/src/distribution/mod.rs b/src/distribution/mod.rs index c31b53d7..807c23a5 100644 --- a/src/distribution/mod.rs +++ b/src/distribution/mod.rs @@ -102,6 +102,9 @@ 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 { @@ -109,6 +112,7 @@ impl core::fmt::Display for InverseCdfError { 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"), } } } diff --git a/src/function/beta.rs b/src/function/beta.rs index eeb5fc06..48992df1 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -336,7 +336,15 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { Ok(beta_reg_from_fraction(log_prefactor, fraction, a)) } -/// Computes the inverse of the regularized incomplete beta function +pub(crate) fn try_inv_beta_reg(a: f64, b: f64, probability: f64) -> Result { + inverse::try_inv_beta_reg(a, b, probability) +} + +/// Computes the inverse of the regularized incomplete beta function. +/// +/// # Panics +/// +/// If the parameters are invalid or the numerical method does not converge. pub fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { inverse::inv_beta_reg(a, b, probability) } @@ -911,6 +919,11 @@ mod tests { assert!(actual.abs_diff(0x3fdf_5753_caf6_9652) <= 1); } + #[test] + fn test_inv_beta_reg_extreme_shape_ratio() { + assert_eq!(inv_beta_reg(1e20, 10.0, 0.3), 1.0); + } + #[test] fn test_inv_beta_reg_small_shape_lower_tail_is_monotone() { let cases = [ diff --git a/src/function/beta/inverse.rs b/src/function/beta/inverse.rs index 74934e42..f6aa45e1 100644 --- a/src/function/beta/inverse.rs +++ b/src/function/beta/inverse.rs @@ -1,5 +1,5 @@ use super::{ - BetaFuncError, beta_continued_fraction, checked_beta_reg, + BetaFuncError, beta_continued_fraction, beta_reg_use_complement, checked_beta_reg, large_params::{self, LogPrefactor}, ln_beta, }; @@ -24,8 +24,7 @@ fn log_regularized_beta(a: f64, b: f64, x: f64, log_beta: f64) -> Result= symmetry_split { + if beta_reg_use_complement(a, b, x) { return checked_beta_reg(a, b, x).map(f64::ln); } @@ -104,19 +103,18 @@ fn inverse_log_beta(a: f64, b: f64) -> f64 { ln_beta(a, b) } -fn solve_lower_tail(a: f64, b: f64, probability: f64) -> f64 { +fn solve_lower_tail(a: f64, b: f64, probability: f64) -> Result { let log_probability = probability.ln(); let log_beta = inverse_log_beta(a, b); let smallest = f64::from_bits(1); let smallest_log = smallest.ln(); - let smallest_log_cdf = log_regularized_beta(a, b, smallest, log_beta) - .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + let smallest_log_cdf = log_regularized_beta(a, b, smallest, log_beta)?; if smallest_log_cdf > log_probability { - return 0.0; + return Ok(0.0); } if smallest_log_cdf == log_probability { - return smallest; + return Ok(smallest); } let mut lower = Point { @@ -137,8 +135,7 @@ fn solve_lower_tail(a: f64, b: f64, probability: f64) -> f64 { for _ in 0..MAX_ITERATIONS { let x = log_x.exp(); - let log_cdf = log_regularized_beta(a, b, x, log_beta) - .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + let log_cdf = log_regularized_beta(a, b, x, log_beta)?; let current = Point { log_x, x, log_cdf }; if log_cdf < log_probability { @@ -146,18 +143,17 @@ fn solve_lower_tail(a: f64, b: f64, probability: f64) -> f64 { } else if log_cdf > log_probability { upper = current; } else { - return x; + return Ok(x); } if upper.x.to_bits() - lower.x.to_bits() <= 8 { - return closest_representable(a, b, log_beta, lower, upper, log_probability) - .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + return closest_representable(a, b, log_beta, lower, upper, log_probability); } let log_tolerance = 32.0 * f64::EPSILON * log_probability.abs().max(1.0); if (lower.log_cdf - log_probability).abs() <= log_tolerance && (upper.log_cdf - log_probability).abs() <= log_tolerance { - return closer_point(lower, upper, log_probability); + return Ok(closer_point(lower, upper, log_probability)); } let log_density = (a - 1.0) * log_x + (b - 1.0) * (-x).ln_1p() - log_beta; @@ -176,13 +172,10 @@ fn solve_lower_tail(a: f64, b: f64, probability: f64) -> f64 { log_x = next; } - panic!( - "inv_beta_reg did not converge for a={a}, b={b}, probability={probability}, lower={:?}, upper={:?}, lower_log_cdf={:?}, upper_log_cdf={:?}", - lower.x, upper.x, lower.log_cdf, upper.log_cdf - ) + Err(BetaFuncError::ConvergenceFailed) } -pub(super) fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { +pub(super) fn try_inv_beta_reg(a: f64, b: f64, probability: f64) -> Result { assert!(a.is_finite() && a > 0.0, "a must be finite and positive"); assert!(b.is_finite() && b > 0.0, "b must be finite and positive"); assert!( @@ -191,17 +184,22 @@ pub(super) fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { ); if probability == 0.0 { - return 0.0; + return Ok(0.0); } if probability == 1.0 { - return 1.0; + return Ok(1.0); } if probability == 0.5 && a == b { - return 0.5; + return Ok(0.5); } if probability <= 0.5 { return solve_lower_tail(a, b, probability); } - 1.0 - solve_lower_tail(b, a, 1.0 - probability) + solve_lower_tail(b, a, 1.0 - probability).map(|value| 1.0 - value) +} + +pub(super) fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { + try_inv_beta_reg(a, b, probability) + .unwrap_or_else(|error| panic!("inv_beta_reg failed: {error}")) } From 94ffb0fd40a6940cf214b1fda74d603c9af2e631 Mon Sep 17 00:00:00 2001 From: Przemyslaw Olszewski Date: Wed, 26 Aug 2026 17:05:27 +0200 Subject: [PATCH 3/5] fix: Prevent inverse beta stalls --- src/distribution/beta.rs | 10 +++++++++- src/function/beta/inverse.rs | 24 ++++++++++++++++++------ 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs index 7cd83727..2640f6be 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -722,8 +722,16 @@ mod tests { } #[test] - fn test_try_inverse_cdf_reports_convergence_failure() { + 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) diff --git a/src/function/beta/inverse.rs b/src/function/beta/inverse.rs index f6aa45e1..0427aff4 100644 --- a/src/function/beta/inverse.rs +++ b/src/function/beta/inverse.rs @@ -47,6 +47,12 @@ fn closer_point(lower: Point, upper: Point, log_probability: f64) -> f64 { } } +fn representable_midpoint(lower: f64, upper: f64) -> Option { + let lower_bits = lower.to_bits(); + let distance = upper.to_bits() - lower_bits; + (distance > 1).then(|| f64::from_bits(lower_bits + distance / 2)) +} + fn closest_representable( a: f64, b: f64, @@ -132,9 +138,10 @@ fn solve_lower_tail(a: f64, b: f64, probability: f64) -> Result Result lower.log_x && newton < upper.log_x { + let next = if newton.is_finite() && newton > lower.log_x && newton < upper.log_x { newton } else { midpoint }; - if next.exp() == x { - next = midpoint; - } - log_x = next; + let next_x = next.exp(); + x = if next_x <= lower.x || next_x >= upper.x { + let Some(value) = representable_midpoint(lower.x, upper.x) else { + return closest_representable(a, b, log_beta, lower, upper, log_probability); + }; + value + } else { + next_x + }; } Err(BetaFuncError::ConvergenceFailed) From 49eeeb6b03d3b29fc3ade0f6efb59e11e8bed6b7 Mon Sep 17 00:00:00 2001 From: Przemyslaw Olszewski Date: Wed, 26 Aug 2026 17:11:59 +0200 Subject: [PATCH 4/5] fix: Validate fallible inverse beta --- src/distribution/beta.rs | 2 +- src/function/beta.rs | 16 ++++++++++++++++ src/function/beta/inverse.rs | 15 +++++++++------ 3 files changed, 26 insertions(+), 7 deletions(-) diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs index 2640f6be..6c81b7c6 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -183,7 +183,7 @@ impl ContinuousCDF 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 /// diff --git a/src/function/beta.rs b/src/function/beta.rs index 48992df1..353c38dd 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -924,6 +924,22 @@ mod tests { assert_eq!(inv_beta_reg(1e20, 10.0, 0.3), 1.0); } + #[test] + fn test_try_inv_beta_reg_rejects_invalid_arguments() { + assert_eq!( + try_inv_beta_reg(0.0, 1.0, 0.5), + Err(BetaFuncError::ANotGreaterThanZero) + ); + assert_eq!( + try_inv_beta_reg(1.0, f64::INFINITY, 0.5), + Err(BetaFuncError::BNotGreaterThanZero) + ); + assert_eq!( + try_inv_beta_reg(1.0, 1.0, f64::NAN), + Err(BetaFuncError::XOutOfRange) + ); + } + #[test] fn test_inv_beta_reg_small_shape_lower_tail_is_monotone() { let cases = [ diff --git a/src/function/beta/inverse.rs b/src/function/beta/inverse.rs index 0427aff4..e31554e5 100644 --- a/src/function/beta/inverse.rs +++ b/src/function/beta/inverse.rs @@ -188,12 +188,15 @@ fn solve_lower_tail(a: f64, b: f64, probability: f64) -> Result Result { - assert!(a.is_finite() && a > 0.0, "a must be finite and positive"); - assert!(b.is_finite() && b > 0.0, "b must be finite and positive"); - assert!( - probability.is_finite() && (0.0..=1.0).contains(&probability), - "probability must be finite and in [0, 1]" - ); + if !a.is_finite() || a <= 0.0 { + return Err(BetaFuncError::ANotGreaterThanZero); + } + if !b.is_finite() || b <= 0.0 { + return Err(BetaFuncError::BNotGreaterThanZero); + } + if !probability.is_finite() || !(0.0..=1.0).contains(&probability) { + return Err(BetaFuncError::XOutOfRange); + } if probability == 0.0 { return Ok(0.0); From e1a19628d10b40d75307283a421351538c57a906 Mon Sep 17 00:00:00 2001 From: Przemyslaw Olszewski Date: Wed, 26 Aug 2026 17:40:49 +0200 Subject: [PATCH 5/5] test: Tighten inverse beta round trips --- src/function/beta.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/function/beta.rs b/src/function/beta.rs index 353c38dd..8d79a350 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -1010,7 +1010,7 @@ mod tests { let quantile = inv_beta_reg(a, b, probability); if quantile > 0.0 && quantile < 1.0 { let recovered = beta_reg(a, b, quantile); - let tolerance = 5e-11 * probability.max(1.0 - probability); + let tolerance = 5e-11 * probability.min(1.0 - probability); assert!( (recovered - probability).abs() <= tolerance, "a={a:?}, b={b:?}, probability={probability:?}, quantile={quantile:?}, recovered={recovered:?}"