libs/math/example/binomial_coinflip_example.cpp
// Copyright Paul A. 2007, 2010 // Copyright John Maddock 2006 // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt // or copy at http://www.boost.org/LICENSE_1_0.txt) // Simple example of computing probabilities and quantiles for // a Bernoulli random variable representing the flipping of a coin. // http://mathworld.wolfram.com/CoinTossing.html // http://en.wikipedia.org/wiki/Bernoulli_trial // Weisstein, Eric W. "Dice." From MathWorld--A Wolfram Web Resource. // http://mathworld.wolfram.com/Dice.html // http://en.wikipedia.org/wiki/Bernoulli_distribution // http://mathworld.wolfram.com/BernoulliDistribution.html // // An idealized coin consists of a circular disk of zero thickness which, // when thrown in the air and allowed to fall, will rest with either side face up // ("heads" H or "tails" T) with equal probability. A coin is therefore a two-sided die. // Despite slight differences between the sides and nonzero thickness of actual coins, // the distribution of their tosses makes a good approximation to a p==1/2 Bernoulli distribution. //[binomial_coinflip_example1 /*`An example of a [@http://en.wikipedia.org/wiki/Bernoulli_process Bernoulli process] is coin flipping. A variable in such a sequence may be called a Bernoulli variable. This example shows using the Binomial distribution to predict the probability of heads and tails when throwing a coin. The number of correct answers (say heads), X, is distributed as a binomial random variable with binomial distribution parameters number of trials (flips) n = 10 and probability (success_fraction) of getting a head p = 0.5 (a 'fair' coin). (Our coin is assumed fair, but we could easily change the success_fraction parameter p from 0.5 to some other value to simulate an unfair coin, say 0.6 for one with chewing gum on the tail, so it is more likely to fall tails down and heads up). First we need some includes and using statements to be able to use the binomial distribution, some std input and output, and get started: */ #include <boost/math/distributions/binomial.hpp> using boost::math::binomial; #include <iostream> using std::cout; using std::endl; using std::left; #include <iomanip> using std::setw; int main() { cout << "Using Binomial distribution to predict how many heads and tails." << endl; try { /*` See note [link coinflip_eg_catch with the catch block] about why a try and catch block is always a good idea. First, construct a binomial distribution with parameters success_fraction 1/2, and how many flips. */ const double success_fraction = 0.5; // = 50% = 1/2 for a 'fair' coin. int flips = 10; binomial flip(flips, success_fraction); cout.precision(4); /*` Then some examples of using Binomial moments (and echoing the parameters). */ cout << "From " << flips << " one can expect to get on average " << mean(flip) << " heads (or tails)." << endl; cout << "Mode is " << mode(flip) << endl; cout << "Standard deviation is " << standard_deviation(flip) << endl; cout << "So about 2/3 will lie within 1 standard deviation and get between " << ceil(mean(flip) - standard_deviation(flip)) << " and " << floor(mean(flip) + standard_deviation(flip)) << " correct." << endl; cout << "Skewness is " << skewness(flip) << endl; // Skewness of binomial distributions is only zero (symmetrical) // if success_fraction is exactly one half, // for example, when flipping 'fair' coins. cout << "Skewness if success_fraction is " << flip.success_fraction() << " is " << skewness(flip) << endl << endl; // Expect zero for a 'fair' coin. /*` Now we show a variety of predictions on the probability of heads: */ cout << "For " << flip.trials() << " coin flips: " << endl; cout << "Probability of getting no heads is " << pdf(flip, 0) << endl; cout << "Probability of getting at least one head is " << 1. - pdf(flip, 0) << endl; /*` When we want to calculate the probability for a range or values we can sum the PDF's: */ cout << "Probability of getting 0 or 1 heads is " << pdf(flip, 0) + pdf(flip, 1) << endl; // sum of exactly == probabilities /*` Or we can use the cdf. */ cout << "Probability of getting 0 or 1 (<= 1) heads is " << cdf(flip, 1) << endl; cout << "Probability of getting 9 or 10 heads is " << pdf(flip, 9) + pdf(flip, 10) << endl; /*` Note that using */ cout << "Probability of getting 9 or 10 heads is " << 1. - cdf(flip, 8) << endl; /*` is less accurate than using the complement */ cout << "Probability of getting 9 or 10 heads is " << cdf(complement(flip, 8)) << endl; /*` Since the subtraction may involve [@http://docs.sun.com/source/806-3568/ncg_goldberg.html cancellation error], where as `cdf(complement(flip, 8))` does not use such a subtraction internally, and so does not exhibit the problem. To get the probability for a range of heads, we can either add the pdfs for each number of heads */ cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is " // P(X == 4) + P(X == 5) + P(X == 6) << pdf(flip, 4) + pdf(flip, 5) + pdf(flip, 6) << endl; /*` But this is probably less efficient than using the cdf */ cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is " // P(X <= 6) - P(X <= 3) == P(X < 4) << cdf(flip, 6) - cdf(flip, 3) << endl; /*` Certainly for a bigger range like, 3 to 7 */ cout << "Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is " // P(X <= 7) - P(X <= 2) == P(X < 3) << cdf(flip, 7) - cdf(flip, 2) << endl; cout << endl; /*` Finally, print two tables of probability for the /exactly/ and /at least/ a number of heads. */ // Print a table of probability for the exactly a number of heads. cout << "Probability of getting exactly (==) heads" << endl; for (int successes = 0; successes <= flips; successes++) { // Say success means getting a head (or equally success means getting a tail). double probability = pdf(flip, successes); cout << left << setw(2) << successes << " " << setw(10) << probability << " or 1 in " << 1. / probability << ", or " << probability * 100. << "%" << endl; } // for i cout << endl; // Tabulate the probability of getting between zero heads and 0 upto 10 heads. cout << "Probability of getting upto (<=) heads" << endl; for (int successes = 0; successes <= flips; successes++) { // Say success means getting a head // (equally success could mean getting a tail). double probability = cdf(flip, successes); // P(X <= heads) cout << setw(2) << successes << " " << setw(10) << left << probability << " or 1 in " << 1. / probability << ", or " << probability * 100. << "%"<< endl; } // for i /*` The last (0 to 10 heads) must, of course, be 100% probability. */ } catch(const std::exception& e) { // /*` [#coinflip_eg_catch] It is always essential to include try & catch blocks because default policies are to throw exceptions on arguments that are out of domain or cause errors like numeric-overflow. Lacking try & catch blocks, the program will abort, whereas the message below from the thrown exception will give some helpful clues as to the cause of the problem. */ std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl; } //] [binomial_coinflip_example1] return 0; } // int main() // Output: //[binomial_coinflip_example_output /*` [pre Using Binomial distribution to predict how many heads and tails. From 10 one can expect to get on average 5 heads (or tails). Mode is 5 Standard deviation is 1.581 So about 2/3 will lie within 1 standard deviation and get between 4 and 6 correct. Skewness is 0 Skewness if success_fraction is 0.5 is 0 For 10 coin flips: Probability of getting no heads is 0.0009766 Probability of getting at least one head is 0.999 Probability of getting 0 or 1 heads is 0.01074 Probability of getting 0 or 1 (<= 1) heads is 0.01074 Probability of getting 9 or 10 heads is 0.01074 Probability of getting 9 or 10 heads is 0.01074 Probability of getting 9 or 10 heads is 0.01074 Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6562 Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6563 Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is 0.8906 Probability of getting exactly (==) heads 0 0.0009766 or 1 in 1024, or 0.09766% 1 0.009766 or 1 in 102.4, or 0.9766% 2 0.04395 or 1 in 22.76, or 4.395% 3 0.1172 or 1 in 8.533, or 11.72% 4 0.2051 or 1 in 4.876, or 20.51% 5 0.2461 or 1 in 4.063, or 24.61% 6 0.2051 or 1 in 4.876, or 20.51% 7 0.1172 or 1 in 8.533, or 11.72% 8 0.04395 or 1 in 22.76, or 4.395% 9 0.009766 or 1 in 102.4, or 0.9766% 10 0.0009766 or 1 in 1024, or 0.09766% Probability of getting upto (<=) heads 0 0.0009766 or 1 in 1024, or 0.09766% 1 0.01074 or 1 in 93.09, or 1.074% 2 0.05469 or 1 in 18.29, or 5.469% 3 0.1719 or 1 in 5.818, or 17.19% 4 0.377 or 1 in 2.653, or 37.7% 5 0.623 or 1 in 1.605, or 62.3% 6 0.8281 or 1 in 1.208, or 82.81% 7 0.9453 or 1 in 1.058, or 94.53% 8 0.9893 or 1 in 1.011, or 98.93% 9 0.999 or 1 in 1.001, or 99.9% 10 1 or 1 in 1, or 100% ] */ //][/binomial_coinflip_example_output]