hypergeometric_cdf.hpp 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  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. #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
  8. #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
  9. #include <boost/math/policies/error_handling.hpp>
  10. #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
  11. namespace boost{ namespace math{ namespace detail{
  12. template <class T, class Policy>
  13. T hypergeometric_cdf_imp(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy& pol)
  14. {
  15. #ifdef BOOST_MSVC
  16. # pragma warning(push)
  17. # pragma warning(disable:4267)
  18. #endif
  19. BOOST_MATH_STD_USING
  20. T result = 0;
  21. T mode = floor(T(r + 1) * T(n + 1) / (N + 2));
  22. if(x < mode)
  23. {
  24. result = hypergeometric_pdf<T>(x, r, n, N, pol);
  25. T diff = result;
  26. unsigned lower_limit = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
  27. while(diff > (invert ? T(1) : result) * tools::epsilon<T>())
  28. {
  29. diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));
  30. result += diff;
  31. BOOST_MATH_INSTRUMENT_VARIABLE(x);
  32. BOOST_MATH_INSTRUMENT_VARIABLE(diff);
  33. BOOST_MATH_INSTRUMENT_VARIABLE(result);
  34. if(x == lower_limit)
  35. break;
  36. --x;
  37. }
  38. }
  39. else
  40. {
  41. invert = !invert;
  42. unsigned upper_limit = (std::min)(r, n);
  43. if(x != upper_limit)
  44. {
  45. ++x;
  46. result = hypergeometric_pdf<T>(x, r, n, N, pol);
  47. T diff = result;
  48. while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>()))
  49. {
  50. diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r));
  51. result += diff;
  52. ++x;
  53. BOOST_MATH_INSTRUMENT_VARIABLE(x);
  54. BOOST_MATH_INSTRUMENT_VARIABLE(diff);
  55. BOOST_MATH_INSTRUMENT_VARIABLE(result);
  56. }
  57. }
  58. }
  59. if(invert)
  60. result = 1 - result;
  61. return result;
  62. #ifdef BOOST_MSVC
  63. # pragma warning(pop)
  64. #endif
  65. }
  66. template <class T, class Policy>
  67. inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy&)
  68. {
  69. BOOST_FPU_EXCEPTION_GUARD
  70. typedef typename tools::promote_args<T>::type result_type;
  71. typedef typename policies::evaluation<result_type, Policy>::type value_type;
  72. typedef typename policies::normalise<
  73. Policy,
  74. policies::promote_float<false>,
  75. policies::promote_double<false>,
  76. policies::discrete_quantile<>,
  77. policies::assert_undefined<> >::type forwarding_policy;
  78. value_type result;
  79. result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, forwarding_policy());
  80. if(result > 1)
  81. {
  82. result = 1;
  83. }
  84. if(result < 0)
  85. {
  86. result = 0;
  87. }
  88. return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_cdf<%1%>(%1%,%1%,%1%,%1%)");
  89. }
  90. }}} // namespaces
  91. #endif