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> check_success_fraction(const char * function,const RealType & p,RealType * result,const 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> check_dist(const char * function,const RealType & p,RealType * result,const Policy &,const mpl::true_ &)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> check_dist(const char *,const RealType &,RealType *,const Policy &,const mpl::false_ &)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> check_dist(const char * function,const RealType & p,RealType * result,const 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> check_dist_and_k(const char * function,const RealType & p,RealType k,RealType * result,const Policy & pol)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> check_dist_and_prob(const char * function,RealType p,RealType prob,RealType * result,const 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 bernoulli_distribution(RealType p=0.5)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 success_fraction() const118 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> range(const bernoulli_distribution<RealType,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> support(const bernoulli_distribution<RealType,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> mean(const bernoulli_distribution<RealType,Policy> & dist)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> variance(const bernoulli_distribution<RealType,Policy> & dist)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 } // pdf 186 187 template <class RealType, class Policy> cdf(const bernoulli_distribution<RealType,Policy> & dist,const RealType & k)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> cdf(const complemented2_type<bernoulli_distribution<RealType,Policy>,RealType> & c)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> quantile(const bernoulli_distribution<RealType,Policy> & dist,const RealType & p)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> quantile(const complemented2_type<bernoulli_distribution<RealType,Policy>,RealType> & c)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> mode(const bernoulli_distribution<RealType,Policy> & dist)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> skewness(const bernoulli_distribution<RealType,Policy> & dist)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> kurtosis_excess(const bernoulli_distribution<RealType,Policy> & dist)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> kurtosis(const bernoulli_distribution<RealType,Policy> & dist)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 337