hypergeometric.hpp 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293
  1. // Copyright 2008 Gautam Sewani
  2. // Copyright 2008 John Maddock
  3. //
  4. // Use, modification and distribution are subject to the
  5. // Boost Software License, Version 1.0.
  6. // (See accompanying file LICENSE_1_0.txt
  7. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  8. #ifndef BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP
  9. #define BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP
  10. #include <boost/math/distributions/detail/common_error_handling.hpp>
  11. #include <boost/math/distributions/complement.hpp>
  12. #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
  13. #include <boost/math/distributions/detail/hypergeometric_cdf.hpp>
  14. #include <boost/math/distributions/detail/hypergeometric_quantile.hpp>
  15. #include <boost/math/special_functions/fpclassify.hpp>
  16. namespace boost { namespace math {
  17. template <class RealType = double, class Policy = policies::policy<> >
  18. class hypergeometric_distribution
  19. {
  20. public:
  21. typedef RealType value_type;
  22. typedef Policy policy_type;
  23. hypergeometric_distribution(unsigned r, unsigned n, unsigned N) // Constructor.
  24. : m_n(n), m_N(N), m_r(r)
  25. {
  26. static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution";
  27. RealType ret;
  28. check_params(function, &ret);
  29. }
  30. // Accessor functions.
  31. unsigned total()const
  32. {
  33. return m_N;
  34. }
  35. unsigned defective()const
  36. {
  37. return m_r;
  38. }
  39. unsigned sample_count()const
  40. {
  41. return m_n;
  42. }
  43. bool check_params(const char* function, RealType* result)const
  44. {
  45. if(m_r > m_N)
  46. {
  47. *result = boost::math::policies::raise_domain_error<RealType>(
  48. function, "Parameter r out of range: must be <= N but got %1%", static_cast<RealType>(m_r), Policy());
  49. return false;
  50. }
  51. if(m_n > m_N)
  52. {
  53. *result = boost::math::policies::raise_domain_error<RealType>(
  54. function, "Parameter n out of range: must be <= N but got %1%", static_cast<RealType>(m_n), Policy());
  55. return false;
  56. }
  57. return true;
  58. }
  59. bool check_x(unsigned x, const char* function, RealType* result)const
  60. {
  61. if(x < static_cast<unsigned>((std::max)(0, (int)(m_n + m_r) - (int)(m_N))))
  62. {
  63. *result = boost::math::policies::raise_domain_error<RealType>(
  64. function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast<RealType>(x), Policy());
  65. return false;
  66. }
  67. if(x > (std::min)(m_r, m_n))
  68. {
  69. *result = boost::math::policies::raise_domain_error<RealType>(
  70. function, "Random variable out of range: must be less than both n and r but got %1%", static_cast<RealType>(x), Policy());
  71. return false;
  72. }
  73. return true;
  74. }
  75. private:
  76. // Data members:
  77. unsigned m_n; // number of items picked
  78. unsigned m_N; // number of "total" items
  79. unsigned m_r; // number of "defective" items
  80. }; // class hypergeometric_distribution
  81. typedef hypergeometric_distribution<double> hypergeometric;
  82. template <class RealType, class Policy>
  83. inline const std::pair<unsigned, unsigned> range(const hypergeometric_distribution<RealType, Policy>& dist)
  84. { // Range of permissible values for random variable x.
  85. #ifdef BOOST_MSVC
  86. # pragma warning(push)
  87. # pragma warning(disable:4267)
  88. #endif
  89. unsigned r = dist.defective();
  90. unsigned n = dist.sample_count();
  91. unsigned N = dist.total();
  92. unsigned l = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
  93. unsigned u = (std::min)(r, n);
  94. return std::pair<unsigned, unsigned>(l, u);
  95. #ifdef BOOST_MSVC
  96. # pragma warning(pop)
  97. #endif
  98. }
  99. template <class RealType, class Policy>
  100. inline const std::pair<unsigned, unsigned> support(const hypergeometric_distribution<RealType, Policy>& d)
  101. {
  102. return range(d);
  103. }
  104. template <class RealType, class Policy>
  105. inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const unsigned& x)
  106. {
  107. static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
  108. RealType result = 0;
  109. if(!dist.check_params(function, &result))
  110. return result;
  111. if(!dist.check_x(x, function, &result))
  112. return result;
  113. return boost::math::detail::hypergeometric_pdf<RealType>(
  114. x, dist.defective(), dist.sample_count(), dist.total(), Policy());
  115. }
  116. template <class RealType, class Policy, class U>
  117. inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)
  118. {
  119. BOOST_MATH_STD_USING
  120. static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
  121. RealType r = static_cast<RealType>(x);
  122. unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
  123. if(u != r)
  124. {
  125. return boost::math::policies::raise_domain_error<RealType>(
  126. function, "Random variable out of range: must be an integer but got %1%", r, Policy());
  127. }
  128. return pdf(dist, u);
  129. }
  130. template <class RealType, class Policy>
  131. inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const unsigned& x)
  132. {
  133. static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
  134. RealType result = 0;
  135. if(!dist.check_params(function, &result))
  136. return result;
  137. if(!dist.check_x(x, function, &result))
  138. return result;
  139. return boost::math::detail::hypergeometric_cdf<RealType>(
  140. x, dist.defective(), dist.sample_count(), dist.total(), false, Policy());
  141. }
  142. template <class RealType, class Policy, class U>
  143. inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)
  144. {
  145. BOOST_MATH_STD_USING
  146. static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
  147. RealType r = static_cast<RealType>(x);
  148. unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
  149. if(u != r)
  150. {
  151. return boost::math::policies::raise_domain_error<RealType>(
  152. function, "Random variable out of range: must be an integer but got %1%", r, Policy());
  153. }
  154. return cdf(dist, u);
  155. }
  156. template <class RealType, class Policy>
  157. inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, unsigned>& c)
  158. {
  159. static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
  160. RealType result = 0;
  161. if(!c.dist.check_params(function, &result))
  162. return result;
  163. if(!c.dist.check_x(c.param, function, &result))
  164. return result;
  165. return boost::math::detail::hypergeometric_cdf<RealType>(
  166. c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), true, Policy());
  167. }
  168. template <class RealType, class Policy, class U>
  169. inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, U>& c)
  170. {
  171. BOOST_MATH_STD_USING
  172. static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
  173. RealType r = static_cast<RealType>(c.param);
  174. unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
  175. if(u != r)
  176. {
  177. return boost::math::policies::raise_domain_error<RealType>(
  178. function, "Random variable out of range: must be an integer but got %1%", r, Policy());
  179. }
  180. return cdf(complement(c.dist, u));
  181. }
  182. template <class RealType, class Policy>
  183. inline RealType quantile(const hypergeometric_distribution<RealType, Policy>& dist, const RealType& p)
  184. {
  185. BOOST_MATH_STD_USING // for ADL of std functions
  186. // Checking function argument
  187. RealType result = 0;
  188. const char* function = "boost::math::quantile(const hypergeometric_distribution<%1%>&, %1%)";
  189. if (false == dist.check_params(function, &result)) return result;
  190. if(false == detail::check_probability(function, p, &result, Policy())) return result;
  191. return static_cast<RealType>(detail::hypergeometric_quantile(p, RealType(1 - p), dist.defective(), dist.sample_count(), dist.total(), Policy()));
  192. } // quantile
  193. template <class RealType, class Policy>
  194. inline RealType quantile(const complemented2_type<hypergeometric_distribution<RealType, Policy>, RealType>& c)
  195. {
  196. BOOST_MATH_STD_USING // for ADL of std functions
  197. // Checking function argument
  198. RealType result = 0;
  199. const char* function = "quantile(const complemented2_type<hypergeometric_distribution<%1%>, %1%>&)";
  200. if (false == c.dist.check_params(function, &result)) return result;
  201. if(false == detail::check_probability(function, c.param, &result, Policy())) return result;
  202. return static_cast<RealType>(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), Policy()));
  203. } // quantile
  204. template <class RealType, class Policy>
  205. inline RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)
  206. {
  207. return static_cast<RealType>(dist.defective() * dist.sample_count()) / dist.total();
  208. } // RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)
  209. template <class RealType, class Policy>
  210. inline RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)
  211. {
  212. RealType r = static_cast<RealType>(dist.defective());
  213. RealType n = static_cast<RealType>(dist.sample_count());
  214. RealType N = static_cast<RealType>(dist.total());
  215. return n * r * (N - r) * (N - n) / (N * N * (N - 1));
  216. } // RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)
  217. template <class RealType, class Policy>
  218. inline RealType mode(const hypergeometric_distribution<RealType, Policy>& dist)
  219. {
  220. BOOST_MATH_STD_USING
  221. RealType r = static_cast<RealType>(dist.defective());
  222. RealType n = static_cast<RealType>(dist.sample_count());
  223. RealType N = static_cast<RealType>(dist.total());
  224. return floor((r + 1) * (n + 1) / (N + 2));
  225. }
  226. template <class RealType, class Policy>
  227. inline RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)
  228. {
  229. BOOST_MATH_STD_USING
  230. RealType r = static_cast<RealType>(dist.defective());
  231. RealType n = static_cast<RealType>(dist.sample_count());
  232. RealType N = static_cast<RealType>(dist.total());
  233. return (N - 2 * r) * sqrt(N - 1) * (N - 2 * n) / (sqrt(n * r * (N - r) * (N - n)) * (N - 2));
  234. } // RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)
  235. template <class RealType, class Policy>
  236. inline RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
  237. {
  238. RealType r = static_cast<RealType>(dist.defective());
  239. RealType n = static_cast<RealType>(dist.sample_count());
  240. RealType N = static_cast<RealType>(dist.total());
  241. RealType t1 = N * N * (N - 1) / (r * (N - 2) * (N - 3) * (N - r));
  242. RealType t2 = (N * (N + 1) - 6 * N * (N - r)) / (n * (N - n))
  243. + 3 * r * (N - r) * (N + 6) / (N * N) - 6;
  244. return t1 * t2;
  245. } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
  246. template <class RealType, class Policy>
  247. inline RealType kurtosis(const hypergeometric_distribution<RealType, Policy>& dist)
  248. {
  249. return kurtosis_excess(dist) + 3;
  250. } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
  251. }} // namespaces
  252. // This include must be at the end, *after* the accessors
  253. // for this distribution have been defined, in order to
  254. // keep compilers that support two-phase lookup happy.
  255. #include <boost/math/distributions/detail/derived_accessors.hpp>
  256. #endif // include guard