To check out this repository please hg clone the following URL, or open the URL using EasyMercurial or your preferred Mercurial client.
The primary repository for this project is hosted at https://github.com/sonic-visualiser/sv-dependency-builds .
This repository is a read-only copy which is updated automatically every hour.
root / any / include / boost / math / distributions / detail / hypergeometric_quantile.hpp @ 160:cff480c41f97
History | View | Annotate | Download (9.8 KB)
| 1 |
// Copyright 2008 John Maddock
|
|---|---|
| 2 |
//
|
| 3 |
// Use, modification and distribution are subject to the
|
| 4 |
// Boost Software License, Version 1.0.
|
| 5 |
// (See accompanying file LICENSE_1_0.txt
|
| 6 |
// or copy at http://www.boost.org/LICENSE_1_0.txt)
|
| 7 |
|
| 8 |
#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
|
| 9 |
#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
|
| 10 |
|
| 11 |
#include <boost/math/policies/error_handling.hpp> |
| 12 |
#include <boost/math/distributions/detail/hypergeometric_pdf.hpp> |
| 13 |
|
| 14 |
namespace boost{ namespace math{ namespace detail{ |
| 15 |
|
| 16 |
template <class T> |
| 17 |
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned /*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 |
|
| 27 |
template <class T> |
| 28 |
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned /*lbound*/, unsigned 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 |
|
| 38 |
template <class T> |
| 39 |
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned 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 |
|
| 46 |
template <class T> |
| 47 |
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned 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 |
|
| 54 |
template <class T> |
| 55 |
inline unsigned round_x_from_p(unsigned x, T /*p*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&) |
| 56 |
{
|
| 57 |
return x;
|
| 58 |
} |
| 59 |
|
| 60 |
template <class T> |
| 61 |
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned /*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 |
|
| 71 |
template <class T> |
| 72 |
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned /*lbound*/, unsigned 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 |
|
| 82 |
template <class T> |
| 83 |
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned 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 |
|
| 90 |
template <class T> |
| 91 |
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned 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 |
|
| 98 |
template <class T> |
| 99 |
inline unsigned round_x_from_q(unsigned x, T /*q*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&) |
| 100 |
{
|
| 101 |
return x;
|
| 102 |
} |
| 103 |
|
| 104 |
template <class T, class Policy> |
| 105 |
unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol) |
| 106 |
{
|
| 107 |
#ifdef BOOST_MSVC
|
| 108 |
# pragma warning(push)
|
| 109 |
# pragma warning(disable:4267) |
| 110 |
#endif
|
| 111 |
typedef typename Policy::discrete_quantile_type discrete_quantile_type; |
| 112 |
BOOST_MATH_STD_USING |
| 113 |
BOOST_FPU_EXCEPTION_GUARD |
| 114 |
T result; |
| 115 |
T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N); |
| 116 |
unsigned base = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N))); |
| 117 |
unsigned 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 |
unsigned 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 to
|
| 153 |
// roundoff errors, get a fresh PDF calculation before proceding:
|
| 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_INSTRUMENT
|
| 166 |
if(diff != 0) |
| 167 |
{
|
| 168 |
BOOST_MATH_INSTRUMENT_VARIABLE(x); |
| 169 |
BOOST_MATH_INSTRUMENT_VARIABLE(diff); |
| 170 |
BOOST_MATH_INSTRUMENT_VARIABLE(result); |
| 171 |
} |
| 172 |
#endif
|
| 173 |
} |
| 174 |
return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
|
| 175 |
} |
| 176 |
else
|
| 177 |
{
|
| 178 |
unsigned 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 to
|
| 201 |
// roundoff errors, get a fresh PDF calculation before proceding:
|
| 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_INSTRUMENT
|
| 212 |
if(diff != 0) |
| 213 |
{
|
| 214 |
BOOST_MATH_INSTRUMENT_VARIABLE(x); |
| 215 |
BOOST_MATH_INSTRUMENT_VARIABLE(diff); |
| 216 |
BOOST_MATH_INSTRUMENT_VARIABLE(result); |
| 217 |
} |
| 218 |
#endif
|
| 219 |
} |
| 220 |
return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
|
| 221 |
} |
| 222 |
#ifdef BOOST_MSVC
|
| 223 |
# pragma warning(pop)
|
| 224 |
#endif
|
| 225 |
} |
| 226 |
|
| 227 |
template <class T, class Policy> |
| 228 |
inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&) |
| 229 |
{
|
| 230 |
BOOST_FPU_EXCEPTION_GUARD |
| 231 |
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 |
}}} // namespaces
|
| 243 |
|
| 244 |
#endif
|
| 245 |
|