From f36a3c8cb6c4e7f1f7dd4ffa8e165de97fffd30c Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 7 May 2026 14:23:25 +0100 Subject: [PATCH 1/5] InverseGaussian: Improve quantile behaviour. Provide a better initial guess for the root, and bracket the root to avoid shooting off to infinity. Fixes: https://github.com/boostorg/math/issues/1392 Refs: https://github.com/scipy/scipy/issues/25096 --- .../math/distributions/inverse_gaussian.hpp | 63 +++++++++++-------- test/test_inverse_gaussian.cpp | 7 +++ 2 files changed, 44 insertions(+), 26 deletions(-) diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp index 20d3b6bdd5..5747e0b907 100644 --- a/include/boost/math/distributions/inverse_gaussian.hpp +++ b/include/boost/math/distributions/inverse_gaussian.hpp @@ -52,6 +52,7 @@ #include #include #include // for erf/erfc. +#include // for erf/erfc. #include #include #include @@ -298,7 +299,7 @@ struct inverse_gaussian_quantile_complement_functor namespace detail { template - BOOST_MATH_GPU_ENABLED inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1) + BOOST_MATH_GPU_ENABLED inline RealType guess_ig(RealType p, RealType q, RealType mu, RealType lambda) { // guess at random variate value x for inverse gaussian quantile. BOOST_MATH_STD_USING using boost::math::policies::policy; @@ -323,27 +324,16 @@ namespace detail x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi)); } else - { // phi < 2 so much less symmetrical with long tail, - // so use gamma distribution as an approximation. - using boost::math::gamma_distribution; - - // Define the distribution, using gamma_nooverflow: - using gamma_nooverflow = gamma_distribution; - - gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - - // R qgamma(0.2, 0.5, 1) = 0.0320923 - RealType qg = quantile(complement(g, p)); - x = lambda / (qg * 2); - // - if (x > mu/2) // x > mu /2? - { // x too large for the gamma approximation to work well. - //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807 - RealType q = quantile(g, p); - // x = mu * exp(q * static_cast(0.1)); // Said to improve at high p - // x = mu * x; // Improves at high p? - x = mu * exp(q / sqrt(phi) - 1/(2 * phi)); - } + { + // Construct a gamma distribution with the same mean and standard deviation, and hope it + // get's us in more or less the right ballpark: + inverse_gaussian_distribution ig(mu, lambda); + RealType m = mean(ig); + RealType v = standard_deviation(ig); + RealType k = m * m / v; + RealType theta = v / m; + gamma_distribution n01(k, theta); + x = p < q ? quantile(n01, p) : quantile(complement(n01, q)); } return x; } // guess_ig @@ -379,7 +369,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut return result; // infinity; } - RealType guess = detail::guess_ig(p, dist.mean(), dist.scale()); + RealType guess = detail::guess_ig(p, 1 - p, dist.mean(), dist.scale()); using boost::math::tools::max_value; RealType min = static_cast(0); // Minimum possible value is bottom of range of distribution. @@ -390,8 +380,19 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut int get_digits = policies::digits();// get digits from policy, boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. using boost::math::tools::newton_raphson_iterate; + inverse_gaussian_quantile_functor func(dist, p); + // + // Our guess is often not very good, so lets bracket the root just in case: + // + RealType f0 = std::get<0>(func(guess)); + if (f0 < 0) + tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); + else + tools::detail::bracket_root_towards_min(func, guess, f0, min, max, max_iter); + if((guess < min) || (guess > max)) + guess = (min + max) / 2; result = - newton_raphson_iterate(inverse_gaussian_quantile_functor(dist, p), guess, min, max, get_digits, max_iter); + newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); if (max_iter >= policies::get_max_root_iterations()) { return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE @@ -455,7 +456,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type(); boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); using boost::math::tools::newton_raphson_iterate; - result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor(c.dist, q), guess, min, max, get_digits, max_iter); + inverse_gaussian_quantile_complement_functor func(c.dist, q); + RealType f0 = std::get<0>(func(guess)); + if (f0 > 0) + tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); + else + tools::detail::bracket_root_towards_min(func, guess, f0, min, max, max_iter); + if ((guess < min) || (guess > max)) + guess = (min + max) / 2; + result = + newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); + result = newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); if (max_iter >= policies::get_max_root_iterations()) { return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index 3825da9397..c0904d8278 100644 --- a/test/test_inverse_gaussian.cpp +++ b/test/test_inverse_gaussian.cpp @@ -252,6 +252,13 @@ BOOST_AUTO_TEST_CASE( test_main ) // but better than R that gives up completely! // R Error in SuppDists::qinverse_gaussian(4.87914430108515e-219, 1, 1) : Infinite value in NewtonRoot() + inverse_gaussian w_big(66.99652081); + BOOST_CHECK_CLOSE_FRACTION( + quantile(w_big, 0.97969), 591.567880739988823, 10 * tolfeweps); + BOOST_CHECK_CLOSE_FRACTION( + quantile(complement(w_big, 1 - 0.97969)), 591.567880739988823, 10 * tolfeweps); + + BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.5), static_cast(0.87878257893544476), tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( From 54fd419a106c7a2de9605c644b8a20aca9d2d6f9 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 7 May 2026 17:40:37 +0100 Subject: [PATCH 2/5] InverseGaussian: Correct conceptual errors. --- include/boost/math/distributions/inverse_gaussian.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp index 5747e0b907..795ad23d4a 100644 --- a/include/boost/math/distributions/inverse_gaussian.hpp +++ b/include/boost/math/distributions/inverse_gaussian.hpp @@ -369,7 +369,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut return result; // infinity; } - RealType guess = detail::guess_ig(p, 1 - p, dist.mean(), dist.scale()); + RealType guess = detail::guess_ig(p, RealType(1 - p), dist.mean(), dist.scale()); using boost::math::tools::max_value; RealType min = static_cast(0); // Minimum possible value is bottom of range of distribution. @@ -384,7 +384,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut // // Our guess is often not very good, so lets bracket the root just in case: // - RealType f0 = std::get<0>(func(guess)); + RealType f0 = boost::math::get<0>(func(guess)); if (f0 < 0) tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); else @@ -456,7 +456,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type(); using boost::math::tools::newton_raphson_iterate; inverse_gaussian_quantile_complement_functor func(c.dist, q); - RealType f0 = std::get<0>(func(guess)); + RealType f0 = boost::math::get<0>(func(guess)); if (f0 > 0) tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); else From def77a8a4176b2aab2ed4a6278f718516a83070b Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 7 May 2026 19:30:20 +0100 Subject: [PATCH 3/5] InverseGaussian: Make bracket_root_towards_min/max non-recursive on CUDA. --- include/boost/math/tools/roots.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index b0b0fc246c..7ccd56ae36 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -456,8 +456,13 @@ namespace detail { if (count) { max = guess; + // + // We can't have this extra recursive tidy up step on CUDA: + // +#ifndef BOOST_MATH_ENABLE_CUDA if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count); +#endif } return guess0 - (max + min) / 2; } @@ -520,8 +525,13 @@ namespace detail { if (count) { min = guess; + // + // We can't have this extra recursive tidy up step on CUDA: + // +#ifndef BOOST_MATH_ENABLE_CUDA if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count); +#endif } return guess0 - (max + min) / 2; } From d09d634209c077c45738aea8c3040c3f6231bb96 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 8 May 2026 10:00:54 +0100 Subject: [PATCH 4/5] InverseGaussian: Update SYCL recursion config. --- include/boost/math/tools/roots.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 7ccd56ae36..af9b76d998 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -459,7 +459,7 @@ namespace detail { // // We can't have this extra recursive tidy up step on CUDA: // -#ifndef BOOST_MATH_ENABLE_CUDA +#if !defined(BOOST_MATH_ENABLE_CUDA) && !defined(BOOST_MATH_ENABLE_SYCL) if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count); #endif @@ -528,7 +528,7 @@ namespace detail { // // We can't have this extra recursive tidy up step on CUDA: // -#ifndef BOOST_MATH_ENABLE_CUDA +#if !defined(BOOST_MATH_ENABLE_CUDA) && !defined(BOOST_MATH_ENABLE_SYCL) if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count); #endif From 90530b35e8c71ce99ddf0892a5a51c945602b762 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 9 May 2026 10:32:20 +0100 Subject: [PATCH 5/5] InverseGaussian: adjust test tolerance for SYCL. --- test/test_inverse_gaussian.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index c0904d8278..0fe64223dc 100644 --- a/test/test_inverse_gaussian.cpp +++ b/test/test_inverse_gaussian.cpp @@ -254,7 +254,7 @@ BOOST_AUTO_TEST_CASE( test_main ) inverse_gaussian w_big(66.99652081); BOOST_CHECK_CLOSE_FRACTION( - quantile(w_big, 0.97969), 591.567880739988823, 10 * tolfeweps); + quantile(w_big, 0.97969), 591.567880739988823, 15 * tolfeweps); BOOST_CHECK_CLOSE_FRACTION( quantile(complement(w_big, 1 - 0.97969)), 591.567880739988823, 10 * tolfeweps);