246 lines · plain
1// Copyright 2008 John Maddock2//3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP9#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP10 11#include <boost/math/policies/error_handling.hpp>12#include <boost/math/distributions/detail/hypergeometric_pdf.hpp>13 14namespace boost{ namespace math{ namespace detail{15 16template <class T>17inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)18{19 if((p < cum * fudge_factor) && (x != lbound))20 {21 BOOST_MATH_INSTRUMENT_VARIABLE(x-1);22 return --x;23 }24 return x;25}26 27template <class T>28inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)29{30 if((cum < p * fudge_factor) && (x != ubound))31 {32 BOOST_MATH_INSTRUMENT_VARIABLE(x+1);33 return ++x;34 }35 return x;36}37 38template <class T>39inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)40{41 if(p >= 0.5)42 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());43 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());44}45 46template <class T>47inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)48{49 if(p >= 0.5)50 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());51 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());52}53 54template <class T>55inline std::uint64_t round_x_from_p(std::uint64_t x, T /*p*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)56{57 return x;58}59 60template <class T>61inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)62{63 if((q * fudge_factor > cum) && (x != lbound))64 {65 BOOST_MATH_INSTRUMENT_VARIABLE(x-1);66 return --x;67 }68 return x;69}70 71template <class T>72inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)73{74 if((q < cum * fudge_factor) && (x != ubound))75 {76 BOOST_MATH_INSTRUMENT_VARIABLE(x+1);77 return ++x;78 }79 return x;80}81 82template <class T>83inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)84{85 if(q < 0.5)86 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());87 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());88}89 90template <class T>91inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)92{93 if(q >= 0.5)94 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());95 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());96}97 98template <class T>99inline std::uint64_t round_x_from_q(std::uint64_t x, T /*q*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)100{101 return x;102}103 104template <class T, class Policy>105std::uint64_t hypergeometric_quantile_imp(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy& pol)106{107#ifdef _MSC_VER108# pragma warning(push)109# pragma warning(disable:4267)110#endif111 typedef typename Policy::discrete_quantile_type discrete_quantile_type;112 BOOST_MATH_STD_USING113 BOOST_FPU_EXCEPTION_GUARD114 T result;115 T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);116 std::uint64_t base = static_cast<std::uint64_t>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));117 std::uint64_t lim = (std::min)(r, n);118 119 BOOST_MATH_INSTRUMENT_VARIABLE(p);120 BOOST_MATH_INSTRUMENT_VARIABLE(q);121 BOOST_MATH_INSTRUMENT_VARIABLE(r);122 BOOST_MATH_INSTRUMENT_VARIABLE(n);123 BOOST_MATH_INSTRUMENT_VARIABLE(N);124 BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);125 BOOST_MATH_INSTRUMENT_VARIABLE(base);126 BOOST_MATH_INSTRUMENT_VARIABLE(lim);127 128 if(p <= 0.5)129 {130 std::uint64_t x = base;131 result = hypergeometric_pdf<T>(x, r, n, N, pol);132 T diff = result;133 if (diff == 0)134 {135 ++x;136 // We want to skip through x values as fast as we can until we start getting non-zero values,137 // otherwise we're just making lots of expensive PDF calls:138 T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)139 + boost::math::lgamma(static_cast<T>(r + 1), pol)140 + boost::math::lgamma(static_cast<T>(N - n + 1), pol)141 + boost::math::lgamma(static_cast<T>(N - r + 1), pol)142 - boost::math::lgamma(static_cast<T>(N + 1), pol)143 - boost::math::lgamma(static_cast<T>(x + 1), pol)144 - boost::math::lgamma(static_cast<T>(n - x + 1), pol)145 - boost::math::lgamma(static_cast<T>(r - x + 1), pol)146 - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);147 while (log_pdf < tools::log_min_value<T>())148 {149 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));150 ++x;151 }152 // By the time we get here, log_pdf may be fairly inaccurate due to153 // roundoff errors, get a fresh PDF calculation before proceeding:154 diff = hypergeometric_pdf<T>(x, r, n, N, pol);155 }156 while(result < p)157 {158 diff = (diff > tools::min_value<T>() * 8)159 ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))160 : hypergeometric_pdf<T>(x + 1, r, n, N, pol);161 if(result + diff / 2 > p)162 break;163 ++x;164 result += diff;165#ifdef BOOST_MATH_INSTRUMENT166 if(diff != 0)167 {168 BOOST_MATH_INSTRUMENT_VARIABLE(x);169 BOOST_MATH_INSTRUMENT_VARIABLE(diff);170 BOOST_MATH_INSTRUMENT_VARIABLE(result);171 }172#endif173 }174 return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());175 }176 else177 {178 std::uint64_t x = lim;179 result = 0;180 T diff = hypergeometric_pdf<T>(x, r, n, N, pol);181 if (diff == 0)182 {183 // We want to skip through x values as fast as we can until we start getting non-zero values,184 // otherwise we're just making lots of expensive PDF calls:185 --x;186 T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)187 + boost::math::lgamma(static_cast<T>(r + 1), pol)188 + boost::math::lgamma(static_cast<T>(N - n + 1), pol)189 + boost::math::lgamma(static_cast<T>(N - r + 1), pol)190 - boost::math::lgamma(static_cast<T>(N + 1), pol)191 - boost::math::lgamma(static_cast<T>(x + 1), pol)192 - boost::math::lgamma(static_cast<T>(n - x + 1), pol)193 - boost::math::lgamma(static_cast<T>(r - x + 1), pol)194 - boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);195 while (log_pdf < tools::log_min_value<T>())196 {197 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));198 --x;199 }200 // By the time we get here, log_pdf may be fairly inaccurate due to201 // roundoff errors, get a fresh PDF calculation before proceeding:202 diff = hypergeometric_pdf<T>(x, r, n, N, pol);203 }204 while(result + diff / 2 < q)205 {206 result += diff;207 diff = (diff > tools::min_value<T>() * 8)208 ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))209 : hypergeometric_pdf<T>(x - 1, r, n, N, pol);210 --x;211#ifdef BOOST_MATH_INSTRUMENT212 if(diff != 0)213 {214 BOOST_MATH_INSTRUMENT_VARIABLE(x);215 BOOST_MATH_INSTRUMENT_VARIABLE(diff);216 BOOST_MATH_INSTRUMENT_VARIABLE(result);217 }218#endif219 }220 return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());221 }222#ifdef _MSC_VER223# pragma warning(pop)224#endif225}226 227template <class T, class Policy>228inline std::uint64_t hypergeometric_quantile(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)229{230 BOOST_FPU_EXCEPTION_GUARD231 typedef typename tools::promote_args<T>::type result_type;232 typedef typename policies::evaluation<result_type, Policy>::type value_type;233 typedef typename policies::normalise<234 Policy,235 policies::promote_float<false>,236 policies::promote_double<false>,237 policies::assert_undefined<> >::type forwarding_policy;238 239 return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());240}241 242}}} // namespaces243 244#endif245 246