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