]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // boost\math\distributions\bernoulli.hpp |
2 | ||
3 | // Copyright John Maddock 2006. | |
4 | // Copyright Paul A. Bristow 2007. | |
5 | ||
6 | // Use, modification and distribution are subject to the | |
7 | // Boost Software License, Version 1.0. | |
8 | // (See accompanying file LICENSE_1_0.txt | |
9 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | ||
11 | // http://en.wikipedia.org/wiki/bernoulli_distribution | |
12 | // http://mathworld.wolfram.com/BernoulliDistribution.html | |
13 | ||
14 | // bernoulli distribution is the discrete probability distribution of | |
15 | // the number (k) of successes, in a single Bernoulli trials. | |
16 | // It is a version of the binomial distribution when n = 1. | |
17 | ||
18 | // But note that the bernoulli distribution | |
19 | // (like others including the poisson, binomial & negative binomial) | |
20 | // is strictly defined as a discrete function: only integral values of k are envisaged. | |
21 | // However because of the method of calculation using a continuous gamma function, | |
22 | // it is convenient to treat it as if a continous function, | |
23 | // and permit non-integral values of k. | |
24 | // To enforce the strict mathematical model, users should use floor or ceil functions | |
25 | // on k outside this function to ensure that k is integral. | |
26 | ||
27 | #ifndef BOOST_MATH_SPECIAL_BERNOULLI_HPP | |
28 | #define BOOST_MATH_SPECIAL_BERNOULLI_HPP | |
29 | ||
30 | #include <boost/math/distributions/fwd.hpp> | |
31 | #include <boost/math/tools/config.hpp> | |
32 | #include <boost/math/distributions/complement.hpp> // complements | |
33 | #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks | |
34 | #include <boost/math/special_functions/fpclassify.hpp> // isnan. | |
35 | ||
36 | #include <utility> | |
37 | ||
38 | namespace boost | |
39 | { | |
40 | namespace math | |
41 | { | |
42 | namespace bernoulli_detail | |
43 | { | |
44 | // Common error checking routines for bernoulli distribution functions: | |
45 | template <class RealType, class Policy> | |
46 | inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& /* pol */) | |
47 | { | |
48 | if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1)) | |
49 | { | |
50 | *result = policies::raise_domain_error<RealType>( | |
51 | function, | |
52 | "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, Policy()); | |
53 | return false; | |
54 | } | |
55 | return true; | |
56 | } | |
57 | template <class RealType, class Policy> | |
58 | inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& /* pol */, const mpl::true_&) | |
59 | { | |
60 | return check_success_fraction(function, p, result, Policy()); | |
61 | } | |
62 | template <class RealType, class Policy> | |
63 | inline bool check_dist(const char* , const RealType& , RealType* , const Policy& /* pol */, const mpl::false_&) | |
64 | { | |
65 | return true; | |
66 | } | |
67 | template <class RealType, class Policy> | |
68 | inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& /* pol */) | |
69 | { | |
70 | return check_dist(function, p, result, Policy(), typename policies::constructor_error_check<Policy>::type()); | |
71 | } | |
72 | ||
73 | template <class RealType, class Policy> | |
74 | inline bool check_dist_and_k(const char* function, const RealType& p, RealType k, RealType* result, const Policy& pol) | |
75 | { | |
76 | if(check_dist(function, p, result, Policy(), typename policies::method_error_check<Policy>::type()) == false) | |
77 | { | |
78 | return false; | |
79 | } | |
80 | if(!(boost::math::isfinite)(k) || !((k == 0) || (k == 1))) | |
81 | { | |
82 | *result = policies::raise_domain_error<RealType>( | |
83 | function, | |
84 | "Number of successes argument is %1%, but must be 0 or 1 !", k, pol); | |
85 | return false; | |
86 | } | |
87 | return true; | |
88 | } | |
89 | template <class RealType, class Policy> | |
90 | inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& /* pol */) | |
91 | { | |
92 | if((check_dist(function, p, result, Policy(), typename policies::method_error_check<Policy>::type()) && detail::check_probability(function, prob, result, Policy())) == false) | |
93 | { | |
94 | return false; | |
95 | } | |
96 | return true; | |
97 | } | |
98 | } // namespace bernoulli_detail | |
99 | ||
100 | ||
101 | template <class RealType = double, class Policy = policies::policy<> > | |
102 | class bernoulli_distribution | |
103 | { | |
104 | public: | |
105 | typedef RealType value_type; | |
106 | typedef Policy policy_type; | |
107 | ||
108 | bernoulli_distribution(RealType p = 0.5) : m_p(p) | |
109 | { // Default probability = half suits 'fair' coin tossing | |
110 | // where probability of heads == probability of tails. | |
111 | RealType result; // of checks. | |
112 | bernoulli_detail::check_dist( | |
113 | "boost::math::bernoulli_distribution<%1%>::bernoulli_distribution", | |
114 | m_p, | |
115 | &result, Policy()); | |
116 | } // bernoulli_distribution constructor. | |
117 | ||
118 | RealType success_fraction() const | |
119 | { // Probability. | |
120 | return m_p; | |
121 | } | |
122 | ||
123 | private: | |
124 | RealType m_p; // success_fraction | |
125 | }; // template <class RealType> class bernoulli_distribution | |
126 | ||
127 | typedef bernoulli_distribution<double> bernoulli; | |
128 | ||
129 | template <class RealType, class Policy> | |
130 | inline const std::pair<RealType, RealType> range(const bernoulli_distribution<RealType, Policy>& /* dist */) | |
131 | { // Range of permissible values for random variable k = {0, 1}. | |
132 | using boost::math::tools::max_value; | |
133 | return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); | |
134 | } | |
135 | ||
136 | template <class RealType, class Policy> | |
137 | inline const std::pair<RealType, RealType> support(const bernoulli_distribution<RealType, Policy>& /* dist */) | |
138 | { // Range of supported values for random variable k = {0, 1}. | |
139 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
140 | return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); | |
141 | } | |
142 | ||
143 | template <class RealType, class Policy> | |
144 | inline RealType mean(const bernoulli_distribution<RealType, Policy>& dist) | |
145 | { // Mean of bernoulli distribution = p (n = 1). | |
146 | return dist.success_fraction(); | |
147 | } // mean | |
148 | ||
149 | // Rely on dereived_accessors quantile(half) | |
150 | //template <class RealType> | |
151 | //inline RealType median(const bernoulli_distribution<RealType, Policy>& dist) | |
152 | //{ // Median of bernoulli distribution is not defined. | |
153 | // return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); | |
154 | //} // median | |
155 | ||
156 | template <class RealType, class Policy> | |
157 | inline RealType variance(const bernoulli_distribution<RealType, Policy>& dist) | |
158 | { // Variance of bernoulli distribution =p * q. | |
159 | return dist.success_fraction() * (1 - dist.success_fraction()); | |
160 | } // variance | |
161 | ||
162 | template <class RealType, class Policy> | |
163 | RealType pdf(const bernoulli_distribution<RealType, Policy>& dist, const RealType& k) | |
164 | { // Probability Density/Mass Function. | |
165 | BOOST_FPU_EXCEPTION_GUARD | |
166 | // Error check: | |
167 | RealType result = 0; // of checks. | |
168 | if(false == bernoulli_detail::check_dist_and_k( | |
169 | "boost::math::pdf(bernoulli_distribution<%1%>, %1%)", | |
170 | dist.success_fraction(), // 0 to 1 | |
171 | k, // 0 or 1 | |
172 | &result, Policy())) | |
173 | { | |
174 | return result; | |
175 | } | |
176 | // Assume k is integral. | |
177 | if (k == 0) | |
178 | { | |
179 | return 1 - dist.success_fraction(); // 1 - p | |
180 | } | |
181 | else // k == 1 | |
182 | { | |
183 | return dist.success_fraction(); // p | |
184 | } | |
185 | ||
186 | ||
187 | template <class RealType, class Policy> | |
188 | inline RealType cdf(const bernoulli_distribution<RealType, Policy>& dist, const RealType& k) | |
189 | { // Cumulative Distribution Function Bernoulli. | |
190 | RealType p = dist.success_fraction(); | |
191 | // Error check: | |
192 | RealType result = 0; | |
193 | if(false == bernoulli_detail::check_dist_and_k( | |
194 | "boost::math::cdf(bernoulli_distribution<%1%>, %1%)", | |
195 | p, | |
196 | k, | |
197 | &result, Policy())) | |
198 | { | |
199 | return result; | |
200 | } | |
201 | if (k == 0) | |
202 | { | |
203 | return 1 - p; | |
204 | } | |
205 | else | |
206 | { // k == 1 | |
207 | return 1; | |
208 | } | |
209 | } // bernoulli cdf | |
210 | ||
211 | template <class RealType, class Policy> | |
212 | inline RealType cdf(const complemented2_type<bernoulli_distribution<RealType, Policy>, RealType>& c) | |
213 | { // Complemented Cumulative Distribution Function bernoulli. | |
214 | RealType const& k = c.param; | |
215 | bernoulli_distribution<RealType, Policy> const& dist = c.dist; | |
216 | RealType p = dist.success_fraction(); | |
217 | // Error checks: | |
218 | RealType result = 0; | |
219 | if(false == bernoulli_detail::check_dist_and_k( | |
220 | "boost::math::cdf(bernoulli_distribution<%1%>, %1%)", | |
221 | p, | |
222 | k, | |
223 | &result, Policy())) | |
224 | { | |
225 | return result; | |
226 | } | |
227 | if (k == 0) | |
228 | { | |
229 | return p; | |
230 | } | |
231 | else | |
232 | { // k == 1 | |
233 | return 0; | |
234 | } | |
235 | } // bernoulli cdf complement | |
236 | ||
237 | template <class RealType, class Policy> | |
238 | inline RealType quantile(const bernoulli_distribution<RealType, Policy>& dist, const RealType& p) | |
239 | { // Quantile or Percent Point Bernoulli function. | |
240 | // Return the number of expected successes k either 0 or 1. | |
241 | // for a given probability p. | |
242 | ||
243 | RealType result = 0; // of error checks: | |
244 | if(false == bernoulli_detail::check_dist_and_prob( | |
245 | "boost::math::quantile(bernoulli_distribution<%1%>, %1%)", | |
246 | dist.success_fraction(), | |
247 | p, | |
248 | &result, Policy())) | |
249 | { | |
250 | return result; | |
251 | } | |
252 | if (p <= (1 - dist.success_fraction())) | |
253 | { // p <= pdf(dist, 0) == cdf(dist, 0) | |
254 | return 0; | |
255 | } | |
256 | else | |
257 | { | |
258 | return 1; | |
259 | } | |
260 | } // quantile | |
261 | ||
262 | template <class RealType, class Policy> | |
263 | inline RealType quantile(const complemented2_type<bernoulli_distribution<RealType, Policy>, RealType>& c) | |
264 | { // Quantile or Percent Point bernoulli function. | |
265 | // Return the number of expected successes k for a given | |
266 | // complement of the probability q. | |
267 | // | |
268 | // Error checks: | |
269 | RealType q = c.param; | |
270 | const bernoulli_distribution<RealType, Policy>& dist = c.dist; | |
271 | RealType result = 0; | |
272 | if(false == bernoulli_detail::check_dist_and_prob( | |
273 | "boost::math::quantile(bernoulli_distribution<%1%>, %1%)", | |
274 | dist.success_fraction(), | |
275 | q, | |
276 | &result, Policy())) | |
277 | { | |
278 | return result; | |
279 | } | |
280 | ||
281 | if (q <= 1 - dist.success_fraction()) | |
282 | { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0) | |
283 | return 1; | |
284 | } | |
285 | else | |
286 | { | |
287 | return 0; | |
288 | } | |
289 | } // quantile complemented. | |
290 | ||
291 | template <class RealType, class Policy> | |
292 | inline RealType mode(const bernoulli_distribution<RealType, Policy>& dist) | |
293 | { | |
294 | return static_cast<RealType>((dist.success_fraction() <= 0.5) ? 0 : 1); // p = 0.5 can be 0 or 1 | |
295 | } | |
296 | ||
297 | template <class RealType, class Policy> | |
298 | inline RealType skewness(const bernoulli_distribution<RealType, Policy>& dist) | |
299 | { | |
300 | BOOST_MATH_STD_USING; // Aid ADL for sqrt. | |
301 | RealType p = dist.success_fraction(); | |
302 | return (1 - 2 * p) / sqrt(p * (1 - p)); | |
303 | } | |
304 | ||
305 | template <class RealType, class Policy> | |
306 | inline RealType kurtosis_excess(const bernoulli_distribution<RealType, Policy>& dist) | |
307 | { | |
308 | RealType p = dist.success_fraction(); | |
309 | // Note Wolfram says this is kurtosis in text, but gamma2 is the kurtosis excess, | |
310 | // and Wikipedia also says this is the kurtosis excess formula. | |
311 | // return (6 * p * p - 6 * p + 1) / (p * (1 - p)); | |
312 | // But Wolfram kurtosis article gives this simpler formula for kurtosis excess: | |
313 | return 1 / (1 - p) + 1/p -6; | |
314 | } | |
315 | ||
316 | template <class RealType, class Policy> | |
317 | inline RealType kurtosis(const bernoulli_distribution<RealType, Policy>& dist) | |
318 | { | |
319 | RealType p = dist.success_fraction(); | |
320 | return 1 / (1 - p) + 1/p -6 + 3; | |
321 | // Simpler than: | |
322 | // return (6 * p * p - 6 * p + 1) / (p * (1 - p)) + 3; | |
323 | } | |
324 | ||
325 | } // namespace math | |
326 | } // namespace boost | |
327 | ||
328 | // This include must be at the end, *after* the accessors | |
329 | // for this distribution have been defined, in order to | |
330 | // keep compilers that support two-phase lookup happy. | |
331 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
332 | ||
333 | #endif // BOOST_MATH_SPECIAL_BERNOULLI_HPP | |
334 | ||
335 | ||
336 |