math/distributions/bernoulli.hpp

0001 // boost\math\distributions\bernoulli.hpp
0002
0003 // Copyright John Maddock 2006.
0004 // Copyright Paul A. Bristow 2007.
0005 // Copyright Matt Borland 2024.
0006
0007 // Use, modification and distribution are subject to the
0008 // Boost Software License, Version 1.0.
0009 // (See accompanying file LICENSE_1_0.txt
0010 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0011
0012 // http://en.wikipedia.org/wiki/bernoulli_distribution
0013 // http://mathworld.wolfram.com/BernoulliDistribution.html
0014
0015 // bernoulli distribution is the discrete probability distribution of
0016 // the number (k) of successes, in a single Bernoulli trials.
0017 // It is a version of the binomial distribution when n = 1.
0018
0019 // But note that the bernoulli distribution
0020 // (like others including the poisson, binomial & negative binomial)
0021 // is strictly defined as a discrete function: only integral values of k are envisaged.
0022 // However because of the method of calculation using a continuous gamma function,
0023 // it is convenient to treat it as if a continuous function,
0024 // and permit non-integral values of k.
0025 // To enforce the strict mathematical model, users should use floor or ceil functions
0026 // on k outside this function to ensure that k is integral.
0027
0028 #ifndef BOOST_MATH_SPECIAL_BERNOULLI_HPP
0029 #define BOOST_MATH_SPECIAL_BERNOULLI_HPP
0030
0031 #include <boost/math/tools/config.hpp>
0032 #include <boost/math/tools/tuple.hpp>
0033 #include <boost/math/tools/type_traits.hpp>
0034 #include <boost/math/distributions/complement.hpp> // complements
0035 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
0036 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
0037 #include <boost/math/policies/policy.hpp>
0038 #include <boost/math/policies/error_handling.hpp>
0039
0040 #ifndef BOOST_MATH_HAS_NVRTC
0041 #include <utility>
0042 #include <boost/math/distributions/fwd.hpp>
0043 #endif
0044
0045 namespace boost
0046 {
0047   namespace math
0048   {
0049     namespace bernoulli_detail
0050     {
0051       // Common error checking routines for bernoulli distribution functions:
0052       template <class RealType, class Policy>
0053       BOOST_MATH_GPU_ENABLED inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& /* pol */)
0054       {
0055         if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1))
0056         {
0057           *result = policies::raise_domain_error<RealType>(
0058             function,
0059             "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, Policy());
0060           return false;
0061         }
0062         return true;
0063       }
0064       template <class RealType, class Policy>
0065       BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& /* pol */, const boost::math::true_type&)
0066       {
0067         return check_success_fraction(function, p, result, Policy());
0068       }
0069       template <class RealType, class Policy>
0070       BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* , const RealType& , RealType* , const Policy& /* pol */, const boost::math::false_type&)
0071       {
0072          return true;
0073       }
0074       template <class RealType, class Policy>
0075       BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& /* pol */)
0076       {
0077          return check_dist(function, p, result, Policy(), typename policies::constructor_error_check<Policy>::type());
0078       }
0079
0080       template <class RealType, class Policy>
0081       BOOST_MATH_GPU_ENABLED inline bool check_dist_and_k(const char* function, const RealType& p, RealType k, RealType* result, const Policy& pol)
0082       {
0083         if(check_dist(function, p, result, Policy(), typename policies::method_error_check<Policy>::type()) == false)
0084         {
0085           return false;
0086         }
0087         if(!(boost::math::isfinite)(k) || !((k == 0) || (k == 1)))
0088         {
0089           *result = policies::raise_domain_error<RealType>(
0090             function,
0091             "Number of successes argument is %1%, but must be 0 or 1 !", k, pol);
0092           return false;
0093         }
0094        return true;
0095       }
0096       template <class RealType, class Policy>
0097       BOOST_MATH_GPU_ENABLED inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& /* pol */)
0098       {
0099         if((check_dist(function, p, result, Policy(), typename policies::method_error_check<Policy>::type()) && detail::check_probability(function, prob, result, Policy())) == false)
0100         {
0101           return false;
0102         }
0103         return true;
0104       }
0105     } // namespace bernoulli_detail
0106
0107
0108     template <class RealType = double, class Policy = policies::policy<> >
0109     class bernoulli_distribution
0110     {
0111     public:
0112       typedef RealType value_type;
0113       typedef Policy policy_type;
0114
0115       BOOST_MATH_GPU_ENABLED bernoulli_distribution(RealType p = 0.5) : m_p(p)
0116       { // Default probability = half suits 'fair' coin tossing
0117         // where probability of heads == probability of tails.
0118         RealType result; // of checks.
0119         bernoulli_detail::check_dist(
0120            "boost::math::bernoulli_distribution<%1%>::bernoulli_distribution",
0121           m_p,
0122           &result, Policy());
0123       } // bernoulli_distribution constructor.
0124
0125       BOOST_MATH_GPU_ENABLED RealType success_fraction() const
0126       { // Probability.
0127         return m_p;
0128       }
0129
0130     private:
0131       RealType m_p; // success_fraction
0132     }; // template <class RealType> class bernoulli_distribution
0133
0134     typedef bernoulli_distribution<double> bernoulli;
0135
0136     #ifdef __cpp_deduction_guides
0137     template <class RealType>
0138     bernoulli_distribution(RealType)->bernoulli_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0139     #endif
0140
0141     template <class RealType, class Policy>
0142     BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const bernoulli_distribution<RealType, Policy>& /* dist */)
0143     { // Range of permissible values for random variable k = {0, 1}.
0144       using boost::math::tools::max_value;
0145       return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
0146     }
0147
0148     template <class RealType, class Policy>
0149     BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const bernoulli_distribution<RealType, Policy>& /* dist */)
0150     { // Range of supported values for random variable k = {0, 1}.
0151       // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0152       return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
0153     }
0154
0155     template <class RealType, class Policy>
0156     BOOST_MATH_GPU_ENABLED inline RealType mean(const bernoulli_distribution<RealType, Policy>& dist)
0157     { // Mean of bernoulli distribution = p (n = 1).
0158       return dist.success_fraction();
0159     } // mean
0160
0161     // Rely on derived_accessors quantile(half)
0162     //template <class RealType>
0163     //inline RealType median(const bernoulli_distribution<RealType, Policy>& dist)
0164     //{ // Median of bernoulli distribution is not defined.
0165     //  return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN());
0166     //} // median
0167
0168     template <class RealType, class Policy>
0169     BOOST_MATH_GPU_ENABLED inline RealType variance(const bernoulli_distribution<RealType, Policy>& dist)
0170     { // Variance of bernoulli distribution =p * q.
0171       return  dist.success_fraction() * (1 - dist.success_fraction());
0172     } // variance
0173
0174     template <class RealType, class Policy>
0175     BOOST_MATH_GPU_ENABLED RealType pdf(const bernoulli_distribution<RealType, Policy>& dist, const RealType& k)
0176     { // Probability Density/Mass Function.
0177       BOOST_FPU_EXCEPTION_GUARD
0178       // Error check:
0179       RealType result = 0; // of checks.
0180       if(false == bernoulli_detail::check_dist_and_k(
0181         "boost::math::pdf(bernoulli_distribution<%1%>, %1%)",
0182         dist.success_fraction(), // 0 to 1
0183         k, // 0 or 1
0184         &result, Policy()))
0185       {
0186         return result;
0187       }
0188       // Assume k is integral.
0189       if (k == 0)
0190       {
0191         return 1 - dist.success_fraction(); // 1 - p
0192       }
0193       else  // k == 1
0194       {
0195         return dist.success_fraction(); // p
0196       }
0197     } // pdf
0198
0199     template <class RealType, class Policy>
0200     BOOST_MATH_GPU_ENABLED inline RealType cdf(const bernoulli_distribution<RealType, Policy>& dist, const RealType& k)
0201     { // Cumulative Distribution Function Bernoulli.
0202       RealType p = dist.success_fraction();
0203       // Error check:
0204       RealType result = 0;
0205       if(false == bernoulli_detail::check_dist_and_k(
0206         "boost::math::cdf(bernoulli_distribution<%1%>, %1%)",
0207         p,
0208         k,
0209         &result, Policy()))
0210       {
0211         return result;
0212       }
0213       if (k == 0)
0214       {
0215         return 1 - p;
0216       }
0217       else
0218       { // k == 1
0219         return 1;
0220       }
0221     } // bernoulli cdf
0222
0223     template <class RealType, class Policy>
0224     BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<bernoulli_distribution<RealType, Policy>, RealType>& c)
0225     { // Complemented Cumulative Distribution Function bernoulli.
0226       RealType const& k = c.param;
0227       bernoulli_distribution<RealType, Policy> const& dist = c.dist;
0228       RealType p = dist.success_fraction();
0229       // Error checks:
0230       RealType result = 0;
0231       if(false == bernoulli_detail::check_dist_and_k(
0232         "boost::math::cdf(bernoulli_distribution<%1%>, %1%)",
0233         p,
0234         k,
0235         &result, Policy()))
0236       {
0237         return result;
0238       }
0239       if (k == 0)
0240       {
0241         return p;
0242       }
0243       else
0244       { // k == 1
0245         return 0;
0246       }
0247     } // bernoulli cdf complement
0248
0249     template <class RealType, class Policy>
0250     BOOST_MATH_GPU_ENABLED inline RealType quantile(const bernoulli_distribution<RealType, Policy>& dist, const RealType& p)
0251     { // Quantile or Percent Point Bernoulli function.
0252       // Return the number of expected successes k either 0 or 1.
0253       // for a given probability p.
0254
0255       RealType result = 0; // of error checks:
0256       if(false == bernoulli_detail::check_dist_and_prob(
0257         "boost::math::quantile(bernoulli_distribution<%1%>, %1%)",
0258         dist.success_fraction(),
0259         p,
0260         &result, Policy()))
0261       {
0262         return result;
0263       }
0264       if (p <= (1 - dist.success_fraction()))
0265       { // p <= pdf(dist, 0) == cdf(dist, 0)
0266         return 0;
0267       }
0268       else
0269       {
0270         return 1;
0271       }
0272     } // quantile
0273
0274     template <class RealType, class Policy>
0275     BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<bernoulli_distribution<RealType, Policy>, RealType>& c)
0276     { // Quantile or Percent Point bernoulli function.
0277       // Return the number of expected successes k for a given
0278       // complement of the probability q.
0279       //
0280       // Error checks:
0281       RealType q = c.param;
0282       const bernoulli_distribution<RealType, Policy>& dist = c.dist;
0283       RealType result = 0;
0284       if(false == bernoulli_detail::check_dist_and_prob(
0285         "boost::math::quantile(bernoulli_distribution<%1%>, %1%)",
0286         dist.success_fraction(),
0287         q,
0288         &result, Policy()))
0289       {
0290         return result;
0291       }
0292
0293       if (q <= 1 - dist.success_fraction())
0294       { // // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
0295         return 1;
0296       }
0297       else
0298       {
0299         return 0;
0300       }
0301     } // quantile complemented.
0302
0303     template <class RealType, class Policy>
0304     BOOST_MATH_GPU_ENABLED inline RealType mode(const bernoulli_distribution<RealType, Policy>& dist)
0305     {
0306       return static_cast<RealType>((dist.success_fraction() <= 0.5) ? 0 : 1); // p = 0.5 can be 0 or 1
0307     }
0308
0309     template <class RealType, class Policy>
0310     BOOST_MATH_GPU_ENABLED inline RealType skewness(const bernoulli_distribution<RealType, Policy>& dist)
0311     {
0312       BOOST_MATH_STD_USING; // Aid ADL for sqrt.
0313       RealType p = dist.success_fraction();
0314       return (1 - 2 * p) / sqrt(p * (1 - p));
0315     }
0316
0317     template <class RealType, class Policy>
0318     BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const bernoulli_distribution<RealType, Policy>& dist)
0319     {
0320       RealType p = dist.success_fraction();
0321       // Note Wolfram says this is kurtosis in text, but gamma2 is the kurtosis excess,
0322       // and Wikipedia also says this is the kurtosis excess formula.
0323       // return (6 * p * p - 6 * p + 1) / (p * (1 - p));
0324       // But Wolfram kurtosis article gives this simpler formula for kurtosis excess:
0325       return 1 / (1 - p) + 1/p -6;
0326     }
0327
0328     template <class RealType, class Policy>
0329     BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const bernoulli_distribution<RealType, Policy>& dist)
0330     {
0331       RealType p = dist.success_fraction();
0332       return 1 / (1 - p) + 1/p -6 + 3;
0333       // Simpler than:
0334       // return (6 * p * p - 6 * p + 1) / (p * (1 - p)) + 3;
0335     }
0336
0337   } // namespace math
0338 } // namespace boost
0339
0340 // This include must be at the end, *after* the accessors
0341 // for this distribution have been defined, in order to
0342 // keep compilers that support two-phase lookup happy.
0343 #include <boost/math/distributions/detail/derived_accessors.hpp>
0344
0345 #endif // BOOST_MATH_SPECIAL_BERNOULLI_HPP
0346
0347
0348