diff options
Diffstat (limited to 'boost/math/distributions/detail/hypergeometric_quantile.hpp')
| -rw-r--r-- | boost/math/distributions/detail/hypergeometric_quantile.hpp | 46 |
1 files changed, 46 insertions, 0 deletions
diff --git a/boost/math/distributions/detail/hypergeometric_quantile.hpp b/boost/math/distributions/detail/hypergeometric_quantile.hpp index a855a4a777..2cdee23c10 100644 --- a/boost/math/distributions/detail/hypergeometric_quantile.hpp +++ b/boost/math/distributions/detail/hypergeometric_quantile.hpp @@ -130,6 +130,29 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned unsigned x = base; result = hypergeometric_pdf<T>(x, r, n, N, pol); T diff = result; + if (diff == 0) + { + ++x; + // We want to skip through x values as fast as we can until we start getting non-zero values, + // otherwise we're just making lots of expensive PDF calls: + T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol) + + boost::math::lgamma(static_cast<T>(r + 1), pol) + + boost::math::lgamma(static_cast<T>(N - n + 1), pol) + + boost::math::lgamma(static_cast<T>(N - r + 1), pol) + - boost::math::lgamma(static_cast<T>(N + 1), pol) + - boost::math::lgamma(static_cast<T>(x + 1), pol) + - boost::math::lgamma(static_cast<T>(n - x + 1), pol) + - boost::math::lgamma(static_cast<T>(r - x + 1), pol) + - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol); + while (log_pdf < tools::log_min_value<T>()) + { + log_pdf += -log(static_cast<T>(x + 1)) + log(static_cast<T>(n - x)) + log(static_cast<T>(r - x)) - log(static_cast<T>(N - n - r + x + 1)); + ++x; + } + // By the time we get here, log_pdf may be fairly inaccurate due to + // roundoff errors, get a fresh PDF calculation before proceding: + diff = hypergeometric_pdf<T>(x, r, n, N, pol); + } while(result < p) { diff = (diff > tools::min_value<T>() * 8) @@ -155,6 +178,29 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned unsigned x = lim; result = 0; T diff = hypergeometric_pdf<T>(x, r, n, N, pol); + if (diff == 0) + { + // We want to skip through x values as fast as we can until we start getting non-zero values, + // otherwise we're just making lots of expensive PDF calls: + --x; + T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol) + + boost::math::lgamma(static_cast<T>(r + 1), pol) + + boost::math::lgamma(static_cast<T>(N - n + 1), pol) + + boost::math::lgamma(static_cast<T>(N - r + 1), pol) + - boost::math::lgamma(static_cast<T>(N + 1), pol) + - boost::math::lgamma(static_cast<T>(x + 1), pol) + - boost::math::lgamma(static_cast<T>(n - x + 1), pol) + - boost::math::lgamma(static_cast<T>(r - x + 1), pol) + - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol); + while (log_pdf < tools::log_min_value<T>()) + { + log_pdf += log(static_cast<T>(x)) - log(static_cast<T>(n - x + 1)) - log(static_cast<T>(r - x + 1)) + log(static_cast<T>(N - n - r + x)); + --x; + } + // By the time we get here, log_pdf may be fairly inaccurate due to + // roundoff errors, get a fresh PDF calculation before proceding: + diff = hypergeometric_pdf<T>(x, r, n, N, pol); + } while(result + diff / 2 < q) { result += diff; |