math/distributions/hypergeometric.hpp

0001 // Copyright 2008 Gautam Sewani
0002 // Copyright 2008 John Maddock
0003 // Copyright 2021 Paul A. Bristow
0004 //
0005 // Use, modification and distribution are subject to the
0006 // Boost Software License, Version 1.0.
0007 // (See accompanying file LICENSE_1_0.txt
0008 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0009
0010 #ifndef BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP
0011 #define BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP
0012
0013 #include <boost/math/distributions/detail/common_error_handling.hpp>
0014 #include <boost/math/distributions/complement.hpp>
0015 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
0016 #include <boost/math/distributions/detail/hypergeometric_cdf.hpp>
0017 #include <boost/math/distributions/detail/hypergeometric_quantile.hpp>
0018 #include <boost/math/special_functions/fpclassify.hpp>
0019 #include <cstdint>
0020
0021 namespace boost { namespace math {
0022
0023    template <class RealType = double, class Policy = policies::policy<> >
0024    class hypergeometric_distribution
0025    {
0026    public:
0027       typedef RealType value_type;
0028       typedef Policy policy_type;
0029
0030       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.
0031          : m_n(n), m_N(N), m_r(r)
0032       {
0033          static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution";
0034          RealType ret;
0035          check_params(function, &ret);
0036       }
0037       // Accessor functions.
0038       std::uint64_t total() const
0039       {
0040          return m_N;
0041       }
0042
0043       std::uint64_t defective() const // successes/failures/events
0044       {
0045          return m_r;
0046       }
0047
0048       std::uint64_t sample_count()const
0049       {
0050          return m_n;
0051       }
0052
0053       bool check_params(const char* function, RealType* result)const
0054       {
0055          if(m_r > m_N)
0056          {
0057             *result = boost::math::policies::raise_domain_error<RealType>(
0058                function, "Parameter r out of range: must be <= N but got %1%", static_cast<RealType>(m_r), Policy());
0059             return false;
0060          }
0061          if(m_n > m_N)
0062          {
0063             *result = boost::math::policies::raise_domain_error<RealType>(
0064                function, "Parameter n out of range: must be <= N but got %1%", static_cast<RealType>(m_n), Policy());
0065             return false;
0066          }
0067          return true;
0068       }
0069       bool check_x(std::uint64_t x, const char* function, RealType* result)const
0070       {
0071          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))))
0072          {
0073             *result = boost::math::policies::raise_domain_error<RealType>(
0074                function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast<RealType>(x), Policy());
0075             return false;
0076          }
0077          if(x > (std::min)(m_r, m_n))
0078          {
0079             *result = boost::math::policies::raise_domain_error<RealType>(
0080                function, "Random variable out of range: must be less than both n and r but got %1%", static_cast<RealType>(x), Policy());
0081             return false;
0082          }
0083          return true;
0084       }
0085
0086    private:
0087       // Data members:
0088       std::uint64_t m_n;  // number of items picked or drawn.
0089       std::uint64_t m_N; // number of "total" items.
0090       std::uint64_t m_r; // number of "defective/successes/failures/events items.
0091
0092    }; // class hypergeometric_distribution
0093
0094    typedef hypergeometric_distribution<double> hypergeometric;
0095
0096    template <class RealType, class Policy>
0097    inline const std::pair<std::uint64_t, std::uint64_t> range(const hypergeometric_distribution<RealType, Policy>& dist)
0098    { // Range of permissible values for random variable x.
0099 #ifdef _MSC_VER
0100 #  pragma warning(push)
0101 #  pragma warning(disable:4267)
0102 #endif
0103       const auto r = dist.defective();
0104       const auto n = dist.sample_count();
0105       const auto N = dist.total();
0106       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)));
0107       const auto u = (std::min)(r, n);
0108       return std::make_pair(l, u);
0109 #ifdef _MSC_VER
0110 #  pragma warning(pop)
0111 #endif
0112    }
0113
0114    template <class RealType, class Policy>
0115    inline const std::pair<std::uint64_t, std::uint64_t> support(const hypergeometric_distribution<RealType, Policy>& d)
0116    {
0117       return range(d);
0118    }
0119
0120    template <class RealType, class Policy>
0121    inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x)
0122    {
0123       static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
0124       RealType result = 0;
0125       if(!dist.check_params(function, &result))
0126          return result;
0127       if(!dist.check_x(x, function, &result))
0128          return result;
0129
0130       return boost::math::detail::hypergeometric_pdf<RealType>(
0131          x, dist.defective(), dist.sample_count(), dist.total(), Policy());
0132    }
0133
0134    template <class RealType, class Policy, class U>
0135    inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)
0136    {
0137       BOOST_MATH_STD_USING
0138       static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
0139       RealType r = static_cast<RealType>(x);
0140       auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
0141       if(u != r)
0142       {
0143          return boost::math::policies::raise_domain_error<RealType>(
0144             function, "Random variable out of range: must be an integer but got %1%", r, Policy());
0145       }
0146       return pdf(dist, u);
0147    }
0148
0149    template <class RealType, class Policy>
0150    inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const std::uint64_t& x)
0151    {
0152       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
0153       RealType result = 0;
0154       if(!dist.check_params(function, &result))
0155          return result;
0156       if(!dist.check_x(x, function, &result))
0157          return result;
0158
0159       return boost::math::detail::hypergeometric_cdf<RealType>(
0160          x, dist.defective(), dist.sample_count(), dist.total(), false, Policy());
0161    }
0162
0163    template <class RealType, class Policy, class U>
0164    inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)
0165    {
0166       BOOST_MATH_STD_USING
0167       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
0168       RealType r = static_cast<RealType>(x);
0169       auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
0170       if(u != r)
0171       {
0172          return boost::math::policies::raise_domain_error<RealType>(
0173             function, "Random variable out of range: must be an integer but got %1%", r, Policy());
0174       }
0175       return cdf(dist, u);
0176    }
0177
0178    template <class RealType, class Policy>
0179    inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, std::uint64_t>& c)
0180    {
0181       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
0182       RealType result = 0;
0183       if(!c.dist.check_params(function, &result))
0184          return result;
0185       if(!c.dist.check_x(c.param, function, &result))
0186          return result;
0187
0188       return boost::math::detail::hypergeometric_cdf<RealType>(
0189          c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), true, Policy());
0190    }
0191
0192    template <class RealType, class Policy, class U>
0193    inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, U>& c)
0194    {
0195       BOOST_MATH_STD_USING
0196       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
0197       RealType r = static_cast<RealType>(c.param);
0198       auto u = static_cast<std::uint64_t>(lltrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
0199       if(u != r)
0200       {
0201          return boost::math::policies::raise_domain_error<RealType>(
0202             function, "Random variable out of range: must be an integer but got %1%", r, Policy());
0203       }
0204       return cdf(complement(c.dist, u));
0205    }
0206
0207    template <class RealType, class Policy>
0208    inline RealType quantile(const hypergeometric_distribution<RealType, Policy>& dist, const RealType& p)
0209    {
0210       BOOST_MATH_STD_USING // for ADL of std functions
0211
0212       // Checking function argument
0213       RealType result = 0;
0214       const char* function = "boost::math::quantile(const hypergeometric_distribution<%1%>&, %1%)";
0215       if (false == dist.check_params(function, &result))
0216          return result;
0217
0218       if(false == detail::check_probability(function, p, &result, Policy()))
0219          return result;
0220
0221       return static_cast<RealType>(detail::hypergeometric_quantile(p, RealType(1 - p), dist.defective(), dist.sample_count(), dist.total(), Policy()));
0222    } // quantile
0223
0224    template <class RealType, class Policy>
0225    inline RealType quantile(const complemented2_type<hypergeometric_distribution<RealType, Policy>, RealType>& c)
0226    {
0227       BOOST_MATH_STD_USING // for ADL of std functions
0228
0229       // Checking function argument
0230       RealType result = 0;
0231       const char* function = "quantile(const complemented2_type<hypergeometric_distribution<%1%>, %1%>&)";
0232       if (false == c.dist.check_params(function, &result))
0233          return result;
0234       if (false == detail::check_probability(function, c.param, &result, Policy()))
0235          return result;
0236
0237       return static_cast<RealType>(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), Policy()));
0238    } // quantile
0239
0240    // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution
0241
0242    template <class RealType, class Policy>
0243    inline RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)
0244    {
0245       return static_cast<RealType>(dist.defective() * dist.sample_count()) / dist.total();
0246    } // RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)
0247
0248    template <class RealType, class Policy>
0249    inline RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)
0250    {
0251       RealType r = static_cast<RealType>(dist.defective());
0252       RealType n = static_cast<RealType>(dist.sample_count());
0253       RealType N = static_cast<RealType>(dist.total());
0254       return n * r  * (N - r) * (N - n) / (N * N * (N - 1));
0255    } // RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)
0256
0257    template <class RealType, class Policy>
0258    inline RealType mode(const hypergeometric_distribution<RealType, Policy>& dist)
0259    {
0260       BOOST_MATH_STD_USING
0261       RealType r = static_cast<RealType>(dist.defective());
0262       RealType n = static_cast<RealType>(dist.sample_count());
0263       RealType N = static_cast<RealType>(dist.total());
0264       return floor((r + 1) * (n + 1) / (N + 2));
0265    }
0266
0267    template <class RealType, class Policy>
0268    inline RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)
0269    {
0270       BOOST_MATH_STD_USING
0271       RealType r = static_cast<RealType>(dist.defective());
0272       RealType n = static_cast<RealType>(dist.sample_count());
0273       RealType N = static_cast<RealType>(dist.total());
0274       return (N - 2 * r) * sqrt(N - 1) * (N - 2 * n) / (sqrt(n * r * (N - r) * (N - n)) * (N - 2));
0275    } // RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)
0276
0277    template <class RealType, class Policy>
0278    inline RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
0279    {
0280       // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution shown as plain text:
0281       //  mean | (m n)/N
0282       //  standard deviation | sqrt((m n(N - m) (N - n))/(N - 1))/N
0283       //  variance | (m n(1 - m/N) (N - n))/((N - 1) N)
0284       //  skewness | (sqrt(N - 1) (N - 2 m) (N - 2 n))/((N - 2) sqrt(m n(N - m) (N - n)))
0285       //  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))
0286      // Kurtosis[HypergeometricDistribution[n, m, N]]
0287       RealType m = static_cast<RealType>(dist.defective()); // Failures or success events. (Also symbols K or M are used).
0288       RealType n = static_cast<RealType>(dist.sample_count()); // draws or trials.
0289       RealType n2 = n * n; // n^2
0290       RealType N = static_cast<RealType>(dist.total()); // Total population from which n draws or trials are made.
0291       RealType N2 = N * N; // N^2
0292       // 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));
0293       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));
0294       // Agrees with kurtosis hypergeometric distribution(50,200,500) kurtosis = 2.96917
0295       // N[kurtosis[hypergeometricdistribution(50,200,500)], 55]  2.969174035736058474901169623721804275002985337280263464
0296       return result;
0297    } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
0298
0299    template <class RealType, class Policy>
0300    inline RealType kurtosis(const hypergeometric_distribution<RealType, Policy>& dist)
0301    {
0302       return kurtosis_excess(dist) + 3;
0303    } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
0304 }} // namespaces
0305
0306 // This include must be at the end, *after* the accessors
0307 // for this distribution have been defined, in order to
0308 // keep compilers that support two-phase lookup happy.
0309 #include <boost/math/distributions/detail/derived_accessors.hpp>
0310
0311 #endif // include guard