123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243 |
- // 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.
- */
- double probability = 0.3;
- double q = quantile(flip, probability);
- std::cout << "Quantile (flip, " << probability << ") = " << q << std::endl; // Quantile (flip, 0.3) = 3
- probability = 0.6;
- q = quantile(flip, probability);
- std::cout << "Quantile (flip, " << probability << ") = " << q << std::endl; // Quantile (flip, 0.6) = 5
- }
- 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]
|