diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp index 20d3b6bdd5..795ad23d4a 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, 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. @@ -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 = boost::math::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 = boost::math::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/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index b0b0fc246c..af9b76d998 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: + // +#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 } 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: + // +#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 } return guess0 - (max + min) / 2; } diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index 3825da9397..0fe64223dc 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, 15 * 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(