312 lines · plain
1// Copyright 2008 Gautam Sewani2// Copyright 2008 John Maddock3// Copyright 2021 Paul A. Bristow4//5// Use, modification and distribution are subject to the6// Boost Software License, Version 1.0.7// (See accompanying file LICENSE_1_0.txt8// or copy at http://www.boost.org/LICENSE_1_0.txt)9 10#ifndef BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP11#define BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP12 13#include <boost/math/distributions/detail/common_error_handling.hpp>14#include <boost/math/distributions/complement.hpp>15#include <boost/math/distributions/detail/hypergeometric_pdf.hpp>16#include <boost/math/distributions/detail/hypergeometric_cdf.hpp>17#include <boost/math/distributions/detail/hypergeometric_quantile.hpp>18#include <boost/math/special_functions/fpclassify.hpp>19#include <cstdint>20 21namespace boost { namespace math {22 23 template <class RealType = double, class Policy = policies::policy<> >24 class hypergeometric_distribution25 {26 public:27 typedef RealType value_type;28 typedef Policy policy_type;29 30 hypergeometric_distribution(std::uint64_t r, std::uint64_t n, std::uint64_t N) // Constructor. r=defective/failures/success, n=trials/draws, N=total population.31 : m_n(n), m_N(N), m_r(r)32 {33 static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution";34 RealType ret;35 check_params(function, &ret);36 }37 // Accessor functions.38 std::uint64_t total() const39 {40 return m_N;41 }42 43 std::uint64_t defective() const // successes/failures/events44 {45 return m_r;46 }47 48 std::uint64_t sample_count()const49 {50 return m_n;51 }52 53 bool check_params(const char* function, RealType* result)const54 {55 if(m_r > m_N)56 {57 *result = boost::math::policies::raise_domain_error<RealType>(58 function, "Parameter r out of range: must be <= N but got %1%", static_cast<RealType>(m_r), Policy());59 return false;60 }61 if(m_n > m_N)62 {63 *result = boost::math::policies::raise_domain_error<RealType>(64 function, "Parameter n out of range: must be <= N but got %1%", static_cast<RealType>(m_n), Policy());65 return false;66 }67 return true;68 }69 bool check_x(std::uint64_t x, const char* function, RealType* result)const70 {71 if(x < static_cast<std::uint64_t>((std::max)(INT64_C(0), static_cast<std::int64_t>(m_n + m_r) - static_cast<std::int64_t>(m_N))))72 {73 *result = boost::math::policies::raise_domain_error<RealType>(74 function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast<RealType>(x), Policy());75 return false;76 }77 if(x > (std::min)(m_r, m_n))78 {79 *result = boost::math::policies::raise_domain_error<RealType>(80 function, "Random variable out of range: must be less than both n and r but got %1%", static_cast<RealType>(x), Policy());81 return false;82 }83 return true;84 }85 86 private:87 // Data members:88 std::uint64_t m_n; // number of items picked or drawn.89 std::uint64_t m_N; // number of "total" items.90 std::uint64_t m_r; // number of "defective/successes/failures/events items.91 92 }; // class hypergeometric_distribution93 94 typedef hypergeometric_distribution<double> hypergeometric;95 96 template <class RealType, class Policy>97 inline const std::pair<std::uint64_t, std::uint64_t> range(const hypergeometric_distribution<RealType, Policy>& dist)98 { // Range of permissible values for random variable x.99#ifdef _MSC_VER100# pragma warning(push)101# pragma warning(disable:4267)102#endif103 const auto r = dist.defective();104 const auto n = dist.sample_count();105 const auto N = dist.total();106 const auto l = static_cast<std::uint64_t>((std::max)(INT64_C(0), static_cast<std::int64_t>(n + r) - static_cast<std::int64_t>(N)));107 const auto u = (std::min)(r, n);108 return std::make_pair(l, u);109#ifdef _MSC_VER110# pragma warning(pop)111#endif112 }113 114 template <class RealType, class Policy>115 inline const std::pair<std::uint64_t, std::uint64_t> support(const hypergeometric_distribution<RealType, Policy>& d)116 {117 return range(d);118 }119 120 template <class RealType, class Policy>121 inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x)122 {123 static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";124 RealType result = 0;125 if(!dist.check_params(function, &result))126 return result;127 if(!dist.check_x(x, function, &result))128 return result;129 130 return boost::math::detail::hypergeometric_pdf<RealType>(131 x, dist.defective(), dist.sample_count(), dist.total(), Policy());132 }133 134 template <class RealType, class Policy, class U>135 inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)136 {137 BOOST_MATH_STD_USING138 static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";139 RealType r = static_cast<RealType>(x);140 auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));141 if(u != r)142 {143 return boost::math::policies::raise_domain_error<RealType>(144 function, "Random variable out of range: must be an integer but got %1%", r, Policy());145 }146 return pdf(dist, u);147 }148 149 template <class RealType, class Policy>150 inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x)151 {152 static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";153 RealType result = 0;154 if(!dist.check_params(function, &result))155 return result;156 if(!dist.check_x(x, function, &result))157 return result;158 159 return boost::math::detail::hypergeometric_cdf<RealType>(160 x, dist.defective(), dist.sample_count(), dist.total(), false, Policy());161 }162 163 template <class RealType, class Policy, class U>164 inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)165 {166 BOOST_MATH_STD_USING167 static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";168 RealType r = static_cast<RealType>(x);169 auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));170 if(u != r)171 {172 return boost::math::policies::raise_domain_error<RealType>(173 function, "Random variable out of range: must be an integer but got %1%", r, Policy());174 }175 return cdf(dist, u);176 }177 178 template <class RealType, class Policy>179 inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, std::uint64_t>& c)180 {181 static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";182 RealType result = 0;183 if(!c.dist.check_params(function, &result))184 return result;185 if(!c.dist.check_x(c.param, function, &result))186 return result;187 188 return boost::math::detail::hypergeometric_cdf<RealType>(189 c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), true, Policy());190 }191 192 template <class RealType, class Policy, class U>193 inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, U>& c)194 {195 BOOST_MATH_STD_USING196 static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";197 RealType r = static_cast<RealType>(c.param);198 auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));199 if(u != r)200 {201 return boost::math::policies::raise_domain_error<RealType>(202 function, "Random variable out of range: must be an integer but got %1%", r, Policy());203 }204 return cdf(complement(c.dist, u));205 }206 207 template <class RealType, class Policy>208 inline RealType quantile(const hypergeometric_distribution<RealType, Policy>& dist, const RealType& p)209 {210 BOOST_MATH_STD_USING // for ADL of std functions211 212 // Checking function argument213 RealType result = 0;214 const char* function = "boost::math::quantile(const hypergeometric_distribution<%1%>&, %1%)";215 if (false == dist.check_params(function, &result)) 216 return result;217 218 if(false == detail::check_probability(function, p, &result, Policy())) 219 return result;220 221 return static_cast<RealType>(detail::hypergeometric_quantile(p, RealType(1 - p), dist.defective(), dist.sample_count(), dist.total(), Policy()));222 } // quantile223 224 template <class RealType, class Policy>225 inline RealType quantile(const complemented2_type<hypergeometric_distribution<RealType, Policy>, RealType>& c)226 {227 BOOST_MATH_STD_USING // for ADL of std functions228 229 // Checking function argument230 RealType result = 0;231 const char* function = "quantile(const complemented2_type<hypergeometric_distribution<%1%>, %1%>&)";232 if (false == c.dist.check_params(function, &result)) 233 return result;234 if (false == detail::check_probability(function, c.param, &result, Policy()))235 return result;236 237 return static_cast<RealType>(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), Policy()));238 } // quantile239 240 // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution241 242 template <class RealType, class Policy>243 inline RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)244 {245 return static_cast<RealType>(dist.defective() * dist.sample_count()) / dist.total();246 } // RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)247 248 template <class RealType, class Policy>249 inline RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)250 {251 RealType r = static_cast<RealType>(dist.defective());252 RealType n = static_cast<RealType>(dist.sample_count());253 RealType N = static_cast<RealType>(dist.total());254 return n * r * (N - r) * (N - n) / (N * N * (N - 1));255 } // RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)256 257 template <class RealType, class Policy>258 inline RealType mode(const hypergeometric_distribution<RealType, Policy>& dist)259 {260 BOOST_MATH_STD_USING261 RealType r = static_cast<RealType>(dist.defective());262 RealType n = static_cast<RealType>(dist.sample_count());263 RealType N = static_cast<RealType>(dist.total());264 return floor((r + 1) * (n + 1) / (N + 2));265 }266 267 template <class RealType, class Policy>268 inline RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)269 {270 BOOST_MATH_STD_USING271 RealType r = static_cast<RealType>(dist.defective());272 RealType n = static_cast<RealType>(dist.sample_count());273 RealType N = static_cast<RealType>(dist.total());274 return (N - 2 * r) * sqrt(N - 1) * (N - 2 * n) / (sqrt(n * r * (N - r) * (N - n)) * (N - 2));275 } // RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)276 277 template <class RealType, class Policy>278 inline RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)279 {280 // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution shown as plain text:281 // mean | (m n)/N282 // standard deviation | sqrt((m n(N - m) (N - n))/(N - 1))/N283 // variance | (m n(1 - m/N) (N - n))/((N - 1) N)284 // skewness | (sqrt(N - 1) (N - 2 m) (N - 2 n))/((N - 2) sqrt(m n(N - m) (N - n)))285 // kurtosis | ((N - 1) N^2 ((3 m(N - m) (n^2 (-N) + (n - 2) N^2 + 6 n(N - n)))/N^2 - 6 n(N - n) + N(N + 1)))/(m n(N - 3) (N - 2) (N - m) (N - n))286 // Kurtosis[HypergeometricDistribution[n, m, N]]287 RealType m = static_cast<RealType>(dist.defective()); // Failures or success events. (Also symbols K or M are used).288 RealType n = static_cast<RealType>(dist.sample_count()); // draws or trials.289 RealType n2 = n * n; // n^2290 RealType N = static_cast<RealType>(dist.total()); // Total population from which n draws or trials are made.291 RealType N2 = N * N; // N^2292 // result = ((N - 1) N^2 ((3 m(N - m) (n^2 (-N) + (n - 2) N^2 + 6 n(N - n)))/N^2 - 6 n(N - n) + N(N + 1)))/(m n(N - 3) (N - 2) (N - m) (N - n));293 RealType result = ((N-1)*N2*((3*m*(N-m)*(n2*(-N)+(n-2)*N2+6*n*(N-n)))/N2-6*n*(N-n)+N*(N+1)))/(m*n*(N-3)*(N-2)*(N-m)*(N-n));294 // Agrees with kurtosis hypergeometric distribution(50,200,500) kurtosis = 2.96917295 // N[kurtosis[hypergeometricdistribution(50,200,500)], 55] 2.969174035736058474901169623721804275002985337280263464296 return result;297 } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)298 299 template <class RealType, class Policy>300 inline RealType kurtosis(const hypergeometric_distribution<RealType, Policy>& dist)301 {302 return kurtosis_excess(dist) + 3;303 } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)304}} // namespaces305 306// This include must be at the end, *after* the accessors307// for this distribution have been defined, in order to308// keep compilers that support two-phase lookup happy.309#include <boost/math/distributions/detail/derived_accessors.hpp>310 311#endif // include guard312