ceph/src/boost/libs/math/example/binomial_coinflip_example.cpp

   1 // Copyright Paul A. 2007, 2010
   2 // Copyright John Maddock 2006
   3
   4 // Use, modification and distribution are subject to the
   5 // Boost Software License, Version 1.0.
   6 // (See accompanying file LICENSE_1_0.txt
   7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 // Simple example of computing probabilities and quantiles for
  10 // a Bernoulli random variable representing the flipping of a coin.
  11
  12 // http://mathworld.wolfram.com/CoinTossing.html
  13 // http://en.wikipedia.org/wiki/Bernoulli_trial
  14 // Weisstein, Eric W. "Dice." From MathWorld--A Wolfram Web Resource.
  15 // http://mathworld.wolfram.com/Dice.html
  16 // http://en.wikipedia.org/wiki/Bernoulli_distribution
  17 // http://mathworld.wolfram.com/BernoulliDistribution.html
  18 //
  19 // An idealized coin consists of a circular disk of zero thickness which,
  20 // when thrown in the air and allowed to fall, will rest with either side face up
  21 // ("heads" H or "tails" T) with equal probability. A coin is therefore a two-sided die.
  22 // Despite slight differences between the sides and nonzero thickness of actual coins,
  23 // the distribution of their tosses makes a good approximation to a p==1/2 Bernoulli distribution.
  24
  25 //[binomial_coinflip_example1
  26
  27 /*`An example of a [@http://en.wikipedia.org/wiki/Bernoulli_process Bernoulli process]
  28 is coin flipping.
  29 A variable in such a sequence may be called a Bernoulli variable.
  30
  31 This example shows using the Binomial distribution to predict the probability
  32 of heads and tails when throwing a coin.
  33
  34 The number of correct answers (say heads),
  35 X, is distributed as a binomial random variable
  36 with binomial distribution parameters number of trials (flips) n = 10 and probability (success_fraction) of getting a head p = 0.5 (a 'fair' coin).
  37
  38 (Our coin is assumed fair, but we could easily change the success_fraction parameter p
  39 from 0.5 to some other value to simulate an unfair coin,
  40 say 0.6 for one with chewing gum on the tail,
  41 so it is more likely to fall tails down and heads up).
  42
  43 First we need some includes and using statements to be able to use the binomial distribution, some std input and output, and get started:
  44 */
  45
  46 #include <boost/math/distributions/binomial.hpp>
  47   using boost::math::binomial;
  48
  49 #include <iostream>
  50   using std::cout;  using std::endl;  using std::left;
  51 #include <iomanip>
  52   using std::setw;
  53
  54 int main()
  55 {
  56   cout << "Using Binomial distribution to predict how many heads and tails." << endl;
  57   try
  58   {
  59 /*`
  60 See note [link coinflip_eg_catch with the catch block]
  61 about why a try and catch block is always a good idea.
  62
  63 First, construct a binomial distribution with parameters success_fraction
  64 1/2, and how many flips.
  65 */
  66     const double success_fraction = 0.5; // = 50% = 1/2 for a 'fair' coin.
  67     int flips = 10;
  68     binomial flip(flips, success_fraction);
  69
  70     cout.precision(4);
  71 /*`
  72  Then some examples of using Binomial moments (and echoing the parameters).
  73 */
  74     cout << "From " << flips << " one can expect to get on average "
  75       << mean(flip) << " heads (or tails)." << endl;
  76     cout << "Mode is " << mode(flip) << endl;
  77     cout << "Standard deviation is " << standard_deviation(flip) << endl;
  78     cout << "So about 2/3 will lie within 1 standard deviation and get between "
  79       <<  ceil(mean(flip) - standard_deviation(flip))  << " and "
  80       << floor(mean(flip) + standard_deviation(flip)) << " correct." << endl;
  81     cout << "Skewness is " << skewness(flip) << endl;
  82     // Skewness of binomial distributions is only zero (symmetrical)
  83     // if success_fraction is exactly one half,
  84     // for example, when flipping 'fair' coins.
  85     cout << "Skewness if success_fraction is " << flip.success_fraction()
  86       << " is " << skewness(flip) << endl << endl; // Expect zero for a 'fair' coin.
  87 /*`
  88 Now we show a variety of predictions on the probability of heads:
  89 */
  90     cout << "For " << flip.trials() << " coin flips: " << endl;
  91     cout << "Probability of getting no heads is " << pdf(flip, 0) << endl;
  92     cout << "Probability of getting at least one head is " << 1. - pdf(flip, 0) << endl;
  93 /*`
  94 When we want to calculate the probability for a range or values we can sum the PDF's:
  95 */
  96     cout << "Probability of getting 0 or 1 heads is "
  97       << pdf(flip, 0) + pdf(flip, 1) << endl; // sum of exactly == probabilities
  98 /*`
  99 Or we can use the cdf.
 100 */
 101     cout << "Probability of getting 0 or 1 (<= 1) heads is " << cdf(flip, 1) << endl;
 102     cout << "Probability of getting 9 or 10 heads is " << pdf(flip, 9) + pdf(flip, 10) << endl;
 103 /*`
 104 Note that using
 105 */
 106     cout << "Probability of getting 9 or 10 heads is " << 1. - cdf(flip, 8) << endl;
 107 /*`
 108 is less accurate than using the complement
 109 */
 110     cout << "Probability of getting 9 or 10 heads is " << cdf(complement(flip, 8)) << endl;
 111 /*`
 112 Since the subtraction may involve
 113 [@http://docs.sun.com/source/806-3568/ncg_goldberg.html cancellation error],
 114 where as `cdf(complement(flip, 8))`
 115 does not use such a subtraction internally, and so does not exhibit the problem.
 116
 117 To get the probability for a range of heads, we can either add the pdfs for each number of heads
 118 */
 119     cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "
 120       //  P(X == 4) + P(X == 5) + P(X == 6)
 121       << pdf(flip, 4) + pdf(flip, 5) + pdf(flip, 6) << endl;
 122 /*`
 123 But this is probably less efficient than using the cdf
 124 */
 125     cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "
 126       // P(X <= 6) - P(X <= 3) == P(X < 4)
 127       << cdf(flip, 6) - cdf(flip, 3) << endl;
 128 /*`
 129 Certainly for a bigger range like, 3 to 7
 130 */
 131     cout << "Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is "
 132       // P(X <= 7) - P(X <= 2) == P(X < 3)
 133       << cdf(flip, 7) - cdf(flip, 2) << endl;
 134     cout << endl;
 135
 136 /*`
 137 Finally, print two tables of probability for the /exactly/ and /at least/ a number of heads.
 138 */
 139     // Print a table of probability for the exactly a number of heads.
 140     cout << "Probability of getting exactly (==) heads" << endl;
 141     for (int successes = 0; successes <= flips; successes++)
 142     { // Say success means getting a head (or equally success means getting a tail).
 143       double probability = pdf(flip, successes);
 144       cout << left << setw(2) << successes << "     " << setw(10)
 145         << probability << " or 1 in " << 1. / probability
 146         << ", or " << probability * 100. << "%" << endl;
 147     } // for i
 148     cout << endl;
 149
 150     // Tabulate the probability of getting between zero heads and 0 upto 10 heads.
 151     cout << "Probability of getting upto (<=) heads" << endl;
 152     for (int successes = 0; successes <= flips; successes++)
 153     { // Say success means getting a head
 154       // (equally success could mean getting a tail).
 155       double probability = cdf(flip, successes); // P(X <= heads)
 156       cout << setw(2) << successes << "        " << setw(10) << left
 157         << probability << " or 1 in " << 1. / probability << ", or "
 158         << probability * 100. << "%"<< endl;
 159     } // for i
 160 /*`
 161 The last (0 to 10 heads) must, of course, be 100% probability.
 162 */
 163   }
 164   catch(const std::exception& e)
 165   {
 166     //
 167     /*`
 168     [#coinflip_eg_catch]
 169     It is always essential to include try & catch blocks because
 170     default policies are to throw exceptions on arguments that
 171     are out of domain or cause errors like numeric-overflow.
 172
 173     Lacking try & catch blocks, the program will abort, whereas the
 174     message below from the thrown exception will give some helpful
 175     clues as to the cause of the problem.
 176     */
 177     std::cout <<
 178       "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
 179   }
 180 //] [binomial_coinflip_example1]
 181   return 0;
 182 } // int main()
 183
 184 // Output:
 185
 186 //[binomial_coinflip_example_output
 187 /*`
 188
 189 [pre
 190 Using Binomial distribution to predict how many heads and tails.
 191 From 10 one can expect to get on average 5 heads (or tails).
 192 Mode is 5
 193 Standard deviation is 1.581
 194 So about 2/3 will lie within 1 standard deviation and get between 4 and 6 correct.
 195 Skewness is 0
 196 Skewness if success_fraction is 0.5 is 0
 197
 198 For 10 coin flips:
 199 Probability of getting no heads is 0.0009766
 200 Probability of getting at least one head is 0.999
 201 Probability of getting 0 or 1 heads is 0.01074
 202 Probability of getting 0 or 1 (<= 1) heads is 0.01074
 203 Probability of getting 9 or 10 heads is 0.01074
 204 Probability of getting 9 or 10 heads is 0.01074
 205 Probability of getting 9 or 10 heads is 0.01074
 206 Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6562
 207 Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6563
 208 Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is 0.8906
 209
 210 Probability of getting exactly (==) heads
 211 0      0.0009766  or 1 in 1024, or 0.09766%
 212 1      0.009766   or 1 in 102.4, or 0.9766%
 213 2      0.04395    or 1 in 22.76, or 4.395%
 214 3      0.1172     or 1 in 8.533, or 11.72%
 215 4      0.2051     or 1 in 4.876, or 20.51%
 216 5      0.2461     or 1 in 4.063, or 24.61%
 217 6      0.2051     or 1 in 4.876, or 20.51%
 218 7      0.1172     or 1 in 8.533, or 11.72%
 219 8      0.04395    or 1 in 22.76, or 4.395%
 220 9      0.009766   or 1 in 102.4, or 0.9766%
 221 10     0.0009766  or 1 in 1024, or 0.09766%
 222
 223 Probability of getting upto (<=) heads
 224 0         0.0009766  or 1 in 1024, or 0.09766%
 225 1         0.01074    or 1 in 93.09, or 1.074%
 226 2         0.05469    or 1 in 18.29, or 5.469%
 227 3         0.1719     or 1 in 5.818, or 17.19%
 228 4         0.377      or 1 in 2.653, or 37.7%
 229 5         0.623      or 1 in 1.605, or 62.3%
 230 6         0.8281     or 1 in 1.208, or 82.81%
 231 7         0.9453     or 1 in 1.058, or 94.53%
 232 8         0.9893     or 1 in 1.011, or 98.93%
 233 9         0.999      or 1 in 1.001, or 99.9%
 234 10        1          or 1 in 1, or 100%
 235 ]
 236 */
 237 //][/binomial_coinflip_example_output]