boost/math/distributions/hypergeometric.hpp
// Copyright 2008 Gautam Sewani // Copyright 2008 John Maddock // Copyright 2021 Paul A. Bristow // // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt // or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP #define BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP #include <boost/math/distributions/detail/common_error_handling.hpp> #include <boost/math/distributions/complement.hpp> #include <boost/math/distributions/detail/hypergeometric_pdf.hpp> #include <boost/math/distributions/detail/hypergeometric_cdf.hpp> #include <boost/math/distributions/detail/hypergeometric_quantile.hpp> #include <boost/math/special_functions/fpclassify.hpp> #include <cstdint> namespace boost { namespace math { template <class RealType = double, class Policy = policies::policy<> > class hypergeometric_distribution { public: typedef RealType value_type; typedef Policy policy_type; 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. : m_n(n), m_N(N), m_r(r) { static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution"; RealType ret; check_params(function, &ret); } // Accessor functions. std::uint64_t total() const { return m_N; } std::uint64_t defective() const // successes/failures/events { return m_r; } std::uint64_t sample_count()const { return m_n; } bool check_params(const char* function, RealType* result)const { if(m_r > m_N) { *result = boost::math::policies::raise_domain_error<RealType>( function, "Parameter r out of range: must be <= N but got %1%", static_cast<RealType>(m_r), Policy()); return false; } if(m_n > m_N) { *result = boost::math::policies::raise_domain_error<RealType>( function, "Parameter n out of range: must be <= N but got %1%", static_cast<RealType>(m_n), Policy()); return false; } return true; } bool check_x(std::uint64_t x, const char* function, RealType* result)const { 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)))) { *result = boost::math::policies::raise_domain_error<RealType>( function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast<RealType>(x), Policy()); return false; } if(x > (std::min)(m_r, m_n)) { *result = boost::math::policies::raise_domain_error<RealType>( function, "Random variable out of range: must be less than both n and r but got %1%", static_cast<RealType>(x), Policy()); return false; } return true; } private: // Data members: std::uint64_t m_n; // number of items picked or drawn. std::uint64_t m_N; // number of "total" items. std::uint64_t m_r; // number of "defective/successes/failures/events items. }; // class hypergeometric_distribution typedef hypergeometric_distribution<double> hypergeometric; template <class RealType, class Policy> inline const std::pair<std::uint64_t, std::uint64_t> range(const hypergeometric_distribution<RealType, Policy>& dist) { // Range of permissible values for random variable x. #ifdef _MSC_VER # pragma warning(push) # pragma warning(disable:4267) #endif const auto r = dist.defective(); const auto n = dist.sample_count(); const auto N = dist.total(); 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))); const auto u = (std::min)(r, n); return std::make_pair(l, u); #ifdef _MSC_VER # pragma warning(pop) #endif } template <class RealType, class Policy> inline const std::pair<std::uint64_t, std::uint64_t> support(const hypergeometric_distribution<RealType, Policy>& d) { return range(d); } template <class RealType, class Policy> inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x) { static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)"; RealType result = 0; if(!dist.check_params(function, &result)) return result; if(!dist.check_x(x, function, &result)) return result; return boost::math::detail::hypergeometric_pdf<RealType>( x, dist.defective(), dist.sample_count(), dist.total(), Policy()); } template <class RealType, class Policy, class U> inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x) { BOOST_MATH_STD_USING static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)"; RealType r = static_cast<RealType>(x); auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type())); if(u != r) { return boost::math::policies::raise_domain_error<RealType>( function, "Random variable out of range: must be an integer but got %1%", r, Policy()); } return pdf(dist, u); } template <class RealType, class Policy> inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x) { static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)"; RealType result = 0; if(!dist.check_params(function, &result)) return result; if(!dist.check_x(x, function, &result)) return result; return boost::math::detail::hypergeometric_cdf<RealType>( x, dist.defective(), dist.sample_count(), dist.total(), false, Policy()); } template <class RealType, class Policy, class U> inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x) { BOOST_MATH_STD_USING static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)"; RealType r = static_cast<RealType>(x); auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type())); if(u != r) { return boost::math::policies::raise_domain_error<RealType>( function, "Random variable out of range: must be an integer but got %1%", r, Policy()); } return cdf(dist, u); } template <class RealType, class Policy> inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, std::uint64_t>& c) { static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)"; RealType result = 0; if(!c.dist.check_params(function, &result)) return result; if(!c.dist.check_x(c.param, function, &result)) return result; return boost::math::detail::hypergeometric_cdf<RealType>( c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), true, Policy()); } template <class RealType, class Policy, class U> inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, U>& c) { BOOST_MATH_STD_USING static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)"; RealType r = static_cast<RealType>(c.param); auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type())); if(u != r) { return boost::math::policies::raise_domain_error<RealType>( function, "Random variable out of range: must be an integer but got %1%", r, Policy()); } return cdf(complement(c.dist, u)); } template <class RealType, class Policy> inline RealType quantile(const hypergeometric_distribution<RealType, Policy>& dist, const RealType& p) { BOOST_MATH_STD_USING // for ADL of std functions // Checking function argument RealType result = 0; const char* function = "boost::math::quantile(const hypergeometric_distribution<%1%>&, %1%)"; if (false == dist.check_params(function, &result)) return result; if(false == detail::check_probability(function, p, &result, Policy())) return result; return static_cast<RealType>(detail::hypergeometric_quantile(p, RealType(1 - p), dist.defective(), dist.sample_count(), dist.total(), Policy())); } // quantile template <class RealType, class Policy> inline RealType quantile(const complemented2_type<hypergeometric_distribution<RealType, Policy>, RealType>& c) { BOOST_MATH_STD_USING // for ADL of std functions // Checking function argument RealType result = 0; const char* function = "quantile(const complemented2_type<hypergeometric_distribution<%1%>, %1%>&)"; if (false == c.dist.check_params(function, &result)) return result; if (false == detail::check_probability(function, c.param, &result, Policy())) return result; return static_cast<RealType>(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), Policy())); } // quantile // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution template <class RealType, class Policy> inline RealType mean(const hypergeometric_distribution<RealType, Policy>& dist) { return static_cast<RealType>(dist.defective() * dist.sample_count()) / dist.total(); } // RealType mean(const hypergeometric_distribution<RealType, Policy>& dist) template <class RealType, class Policy> inline RealType variance(const hypergeometric_distribution<RealType, Policy>& dist) { RealType r = static_cast<RealType>(dist.defective()); RealType n = static_cast<RealType>(dist.sample_count()); RealType N = static_cast<RealType>(dist.total()); return n * r * (N - r) * (N - n) / (N * N * (N - 1)); } // RealType variance(const hypergeometric_distribution<RealType, Policy>& dist) template <class RealType, class Policy> inline RealType mode(const hypergeometric_distribution<RealType, Policy>& dist) { BOOST_MATH_STD_USING RealType r = static_cast<RealType>(dist.defective()); RealType n = static_cast<RealType>(dist.sample_count()); RealType N = static_cast<RealType>(dist.total()); return floor((r + 1) * (n + 1) / (N + 2)); } template <class RealType, class Policy> inline RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist) { BOOST_MATH_STD_USING RealType r = static_cast<RealType>(dist.defective()); RealType n = static_cast<RealType>(dist.sample_count()); RealType N = static_cast<RealType>(dist.total()); return (N - 2 * r) * sqrt(N - 1) * (N - 2 * n) / (sqrt(n * r * (N - r) * (N - n)) * (N - 2)); } // RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist) template <class RealType, class Policy> inline RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist) { // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution shown as plain text: // mean | (m n)/N // standard deviation | sqrt((m n(N - m) (N - n))/(N - 1))/N // variance | (m n(1 - m/N) (N - n))/((N - 1) N) // skewness | (sqrt(N - 1) (N - 2 m) (N - 2 n))/((N - 2) sqrt(m n(N - m) (N - n))) // 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)) // Kurtosis[HypergeometricDistribution[n, m, N]] RealType m = static_cast<RealType>(dist.defective()); // Failures or success events. (Also symbols K or M are used). RealType n = static_cast<RealType>(dist.sample_count()); // draws or trials. RealType n2 = n * n; // n^2 RealType N = static_cast<RealType>(dist.total()); // Total population from which n draws or trials are made. RealType N2 = N * N; // N^2 // 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)); 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)); // Agrees with kurtosis hypergeometric distribution(50,200,500) kurtosis = 2.96917 // N[kurtosis[hypergeometricdistribution(50,200,500)], 55] 2.969174035736058474901169623721804275002985337280263464 return result; } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist) template <class RealType, class Policy> inline RealType kurtosis(const hypergeometric_distribution<RealType, Policy>& dist) { return kurtosis_excess(dist) + 3; } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist) }} // namespaces // This include must be at the end, *after* the accessors // for this distribution have been defined, in order to // keep compilers that support two-phase lookup happy. #include <boost/math/distributions/detail/derived_accessors.hpp> #endif // include guard