ceph/src/boost/libs/math/include/boost/math/distributions/bernoulli.hpp

   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     } // pdf
 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