Skip to content

InverseGaussian: Improve quantile behaviour. #1395

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jzmaddock wants to merge 5 commits into
base: develop
Choose a base branch
from
+54 −26
Open

InverseGaussian: Improve quantile behaviour. #1395

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f36a3c8
InverseGaussian: Improve quantile behaviour.
jzmaddock May 7, 2026
54fd419
InverseGaussian: Correct conceptual errors.
jzmaddock May 7, 2026
def77a8
InverseGaussian: Make bracket_root_towards_min/max non-recursive on C…
jzmaddock May 7, 2026
d09d634
InverseGaussian: Update SYCL recursion config.
jzmaddock May 8, 2026
c12fe11
Merge branch 'develop' into issue1392
jzmaddock May 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 37 additions & 26 deletions include/boost/math/distributions/inverse_gaussian.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@
#include <boost/math/tools/config.hpp>
#include <boost/math/tools/cstdint.hpp>
#include <boost/math/special_functions/erf.hpp> // for erf/erfc.
#include <boost/math/special_functions/pow.hpp> // for erf/erfc.
#include <boost/math/distributions/complement.hpp>
#include <boost/math/distributions/detail/common_error_handling.hpp>
#include <boost/math/distributions/normal.hpp>
Expand Down Expand Up @@ -298,7 +299,7 @@ struct inverse_gaussian_quantile_complement_functor
namespace detail
{
template <class RealType>
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;
Expand All @@ -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<RealType, no_overthrow_policy>;

gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(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<RealType>(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<RealType> ig(mu, lambda);
RealType m = mean(ig);
RealType v = standard_deviation(ig);
RealType k = m * m / v;
RealType theta = v / m;
gamma_distribution<RealType> n01(k, theta);
x = p < q ? quantile(n01, p) : quantile(complement(n01, q));
}
return x;
} // guess_ig
Expand Down Expand Up @@ -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<RealType>(0); // Minimum possible value is bottom of range of distribution.
Expand All @@ -390,8 +380,19 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut
int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
using boost::math::tools::newton_raphson_iterate;
inverse_gaussian_quantile_functor<RealType, Policy> 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<RealType, Policy>(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<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
Expand Down Expand Up @@ -455,7 +456,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<inverse
if(false == detail::check_probability(function, q, &result, Policy()))
return result;

RealType guess = detail::guess_ig(q, mean, scale);
RealType guess = detail::guess_ig(RealType(1 - q), q, mean, scale);
// Complement.
using boost::math::tools::max_value;

Expand All @@ -466,7 +467,17 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<inverse
int get_digits = policies::digits<RealType, Policy>();
boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
using boost::math::tools::newton_raphson_iterate;
result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, max_iter);
inverse_gaussian_quantile_complement_functor<RealType, Policy> 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<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE
Expand Down
10 changes: 10 additions & 0 deletions include/boost/math/tools/roots.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand Down
7 changes: 7 additions & 0 deletions test/test_inverse_gaussian.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>(0.87878257893544476), tolfeweps); // pdf
BOOST_CHECK_CLOSE_FRACTION(
Expand Down
Loading