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