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