350 lines · plain
1// boost\math\distributions\bernoulli.hpp2 3// Copyright John Maddock 2006.4// Copyright Paul A. Bristow 2007.5// Copyright Matt Borland 2024.6 7// Use, modification and distribution are subject to the8// Boost Software License, Version 1.0.9// (See accompanying file LICENSE_1_0.txt10// or copy at http://www.boost.org/LICENSE_1_0.txt)11 12// http://en.wikipedia.org/wiki/bernoulli_distribution13// http://mathworld.wolfram.com/BernoulliDistribution.html14 15// bernoulli distribution is the discrete probability distribution of16// the number (k) of successes, in a single Bernoulli trials.17// It is a version of the binomial distribution when n = 1.18 19// But note that the bernoulli distribution20// (like others including the poisson, binomial & negative binomial)21// is strictly defined as a discrete function: only integral values of k are envisaged.22// However because of the method of calculation using a continuous gamma function,23// it is convenient to treat it as if a continuous function,24// and permit non-integral values of k.25// To enforce the strict mathematical model, users should use floor or ceil functions26// on k outside this function to ensure that k is integral.27 28#ifndef BOOST_MATH_SPECIAL_BERNOULLI_HPP29#define BOOST_MATH_SPECIAL_BERNOULLI_HPP30 31#include <boost/math/tools/config.hpp>32#include <boost/math/tools/tuple.hpp>33#include <boost/math/tools/type_traits.hpp>34#include <boost/math/tools/promotion.hpp>35#include <boost/math/distributions/complement.hpp> // complements36#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks37#include <boost/math/special_functions/fpclassify.hpp> // isnan.38#include <boost/math/policies/policy.hpp>39#include <boost/math/policies/error_handling.hpp>40 41#ifndef BOOST_MATH_HAS_NVRTC42#include <utility>43#include <boost/math/distributions/fwd.hpp>44#endif45 46namespace boost47{48 namespace math49 {50 namespace bernoulli_detail51 {52 // Common error checking routines for bernoulli distribution functions:53 template <class RealType, class Policy>54 BOOST_MATH_GPU_ENABLED inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& /* pol */)55 {56 if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1))57 {58 *result = policies::raise_domain_error<RealType>(59 function,60 "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, Policy());61 return false;62 }63 return true;64 }65 template <class RealType, class Policy>66 BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& /* pol */, const boost::math::true_type&)67 {68 return check_success_fraction(function, p, result, Policy());69 }70 template <class RealType, class Policy>71 BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* , const RealType& , RealType* , const Policy& /* pol */, const boost::math::false_type&)72 {73 return true;74 }75 template <class RealType, class Policy>76 BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& /* pol */)77 {78 return check_dist(function, p, result, Policy(), typename policies::constructor_error_check<Policy>::type());79 }80 81 template <class RealType, class Policy>82 BOOST_MATH_GPU_ENABLED inline bool check_dist_and_k(const char* function, const RealType& p, RealType k, RealType* result, const Policy& pol)83 {84 if(check_dist(function, p, result, Policy(), typename policies::method_error_check<Policy>::type()) == false)85 {86 return false;87 }88 if(!(boost::math::isfinite)(k) || !((k == 0) || (k == 1)))89 {90 *result = policies::raise_domain_error<RealType>(91 function,92 "Number of successes argument is %1%, but must be 0 or 1 !", k, pol);93 return false;94 }95 return true;96 }97 template <class RealType, class Policy>98 BOOST_MATH_GPU_ENABLED inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& /* pol */)99 {100 if((check_dist(function, p, result, Policy(), typename policies::method_error_check<Policy>::type()) && detail::check_probability(function, prob, result, Policy())) == false)101 {102 return false;103 }104 return true;105 }106 } // namespace bernoulli_detail107 108 109 template <class RealType = double, class Policy = policies::policy<> >110 class bernoulli_distribution111 {112 public:113 typedef RealType value_type;114 typedef Policy policy_type;115 116 BOOST_MATH_GPU_ENABLED bernoulli_distribution(RealType p = 0.5) : m_p(p)117 { // Default probability = half suits 'fair' coin tossing118 // where probability of heads == probability of tails.119 RealType result; // of checks.120 bernoulli_detail::check_dist(121 "boost::math::bernoulli_distribution<%1%>::bernoulli_distribution",122 m_p,123 &result, Policy());124 } // bernoulli_distribution constructor.125 126 BOOST_MATH_GPU_ENABLED RealType success_fraction() const127 { // Probability.128 return m_p;129 }130 131 private:132 RealType m_p; // success_fraction133 }; // template <class RealType> class bernoulli_distribution134 135 typedef bernoulli_distribution<double> bernoulli;136 137 #ifdef __cpp_deduction_guides138 template <class RealType>139 bernoulli_distribution(RealType)->bernoulli_distribution<typename boost::math::tools::promote_args<RealType>::type>;140 #endif141 142 template <class RealType, class Policy>143 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const bernoulli_distribution<RealType, Policy>& /* dist */)144 { // Range of permissible values for random variable k = {0, 1}.145 using boost::math::tools::max_value;146 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));147 }148 149 template <class RealType, class Policy>150 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const bernoulli_distribution<RealType, Policy>& /* dist */)151 { // Range of supported values for random variable k = {0, 1}.152 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.153 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));154 }155 156 template <class RealType, class Policy>157 BOOST_MATH_GPU_ENABLED inline RealType mean(const bernoulli_distribution<RealType, Policy>& dist)158 { // Mean of bernoulli distribution = p (n = 1).159 return dist.success_fraction();160 } // mean161 162 // Rely on derived_accessors quantile(half)163 //template <class RealType>164 //inline RealType median(const bernoulli_distribution<RealType, Policy>& dist)165 //{ // Median of bernoulli distribution is not defined.166 // return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN());167 //} // median168 169 template <class RealType, class Policy>170 BOOST_MATH_GPU_ENABLED inline RealType variance(const bernoulli_distribution<RealType, Policy>& dist)171 { // Variance of bernoulli distribution =p * q.172 return dist.success_fraction() * (1 - dist.success_fraction());173 } // variance174 175 template <class RealType, class Policy>176 BOOST_MATH_GPU_ENABLED RealType pdf(const bernoulli_distribution<RealType, Policy>& dist, const RealType& k)177 { // Probability Density/Mass Function.178 BOOST_FPU_EXCEPTION_GUARD179 // Error check:180 RealType result = 0; // of checks.181 if(false == bernoulli_detail::check_dist_and_k(182 "boost::math::pdf(bernoulli_distribution<%1%>, %1%)",183 dist.success_fraction(), // 0 to 1184 k, // 0 or 1185 &result, Policy()))186 {187 return result;188 }189 // Assume k is integral.190 if (k == 0)191 {192 return 1 - dist.success_fraction(); // 1 - p193 }194 else // k == 1195 {196 return dist.success_fraction(); // p197 }198 } // pdf199 200 template <class RealType, class Policy>201 BOOST_MATH_GPU_ENABLED inline RealType cdf(const bernoulli_distribution<RealType, Policy>& dist, const RealType& k)202 { // Cumulative Distribution Function Bernoulli.203 RealType p = dist.success_fraction();204 // Error check:205 RealType result = 0;206 if(false == bernoulli_detail::check_dist_and_k(207 "boost::math::cdf(bernoulli_distribution<%1%>, %1%)",208 p,209 k,210 &result, Policy()))211 {212 return result;213 }214 if (k == 0)215 {216 return 1 - p;217 }218 else219 { // k == 1220 return 1;221 }222 } // bernoulli cdf223 224 template <class RealType, class Policy>225 BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<bernoulli_distribution<RealType, Policy>, RealType>& c)226 { // Complemented Cumulative Distribution Function bernoulli.227 RealType const& k = c.param;228 bernoulli_distribution<RealType, Policy> const& dist = c.dist;229 RealType p = dist.success_fraction();230 // Error checks:231 RealType result = 0;232 if(false == bernoulli_detail::check_dist_and_k(233 "boost::math::cdf(bernoulli_distribution<%1%>, %1%)",234 p,235 k,236 &result, Policy()))237 {238 return result;239 }240 if (k == 0)241 {242 return p;243 }244 else245 { // k == 1246 return 0;247 }248 } // bernoulli cdf complement249 250 template <class RealType, class Policy>251 BOOST_MATH_GPU_ENABLED inline RealType quantile(const bernoulli_distribution<RealType, Policy>& dist, const RealType& p)252 { // Quantile or Percent Point Bernoulli function.253 // Return the number of expected successes k either 0 or 1.254 // for a given probability p.255 256 RealType result = 0; // of error checks:257 if(false == bernoulli_detail::check_dist_and_prob(258 "boost::math::quantile(bernoulli_distribution<%1%>, %1%)",259 dist.success_fraction(),260 p,261 &result, Policy()))262 {263 return result;264 }265 if (p <= (1 - dist.success_fraction()))266 { // p <= pdf(dist, 0) == cdf(dist, 0)267 return 0;268 }269 else270 {271 return 1;272 }273 } // quantile274 275 template <class RealType, class Policy>276 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<bernoulli_distribution<RealType, Policy>, RealType>& c)277 { // Quantile or Percent Point bernoulli function.278 // Return the number of expected successes k for a given279 // complement of the probability q.280 //281 // Error checks:282 RealType q = c.param;283 const bernoulli_distribution<RealType, Policy>& dist = c.dist;284 RealType result = 0;285 if(false == bernoulli_detail::check_dist_and_prob(286 "boost::math::quantile(bernoulli_distribution<%1%>, %1%)",287 dist.success_fraction(),288 q,289 &result, Policy()))290 {291 return result;292 }293 294 if (q <= 1 - dist.success_fraction())295 { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0)296 return 1;297 }298 else299 {300 return 0;301 }302 } // quantile complemented.303 304 template <class RealType, class Policy>305 BOOST_MATH_GPU_ENABLED inline RealType mode(const bernoulli_distribution<RealType, Policy>& dist)306 {307 return static_cast<RealType>((dist.success_fraction() <= 0.5) ? 0 : 1); // p = 0.5 can be 0 or 1308 }309 310 template <class RealType, class Policy>311 BOOST_MATH_GPU_ENABLED inline RealType skewness(const bernoulli_distribution<RealType, Policy>& dist)312 {313 BOOST_MATH_STD_USING; // Aid ADL for sqrt.314 RealType p = dist.success_fraction();315 return (1 - 2 * p) / sqrt(p * (1 - p));316 }317 318 template <class RealType, class Policy>319 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const bernoulli_distribution<RealType, Policy>& dist)320 {321 RealType p = dist.success_fraction();322 // Note Wolfram says this is kurtosis in text, but gamma2 is the kurtosis excess,323 // and Wikipedia also says this is the kurtosis excess formula.324 // return (6 * p * p - 6 * p + 1) / (p * (1 - p));325 // But Wolfram kurtosis article gives this simpler formula for kurtosis excess:326 return 1 / (1 - p) + 1/p -6;327 }328 329 template <class RealType, class Policy>330 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const bernoulli_distribution<RealType, Policy>& dist)331 {332 RealType p = dist.success_fraction();333 return 1 / (1 - p) + 1/p -6 + 3;334 // Simpler than:335 // return (6 * p * p - 6 * p + 1) / (p * (1 - p)) + 3;336 }337 338 } // namespace math339} // namespace boost340 341// This include must be at the end, *after* the accessors342// for this distribution have been defined, in order to343// keep compilers that support two-phase lookup happy.344#include <boost/math/distributions/detail/derived_accessors.hpp>345 346#endif // BOOST_MATH_SPECIAL_BERNOULLI_HPP347 348 349 350