1 // Copyright Paul A. 2007, 2010
2 // Copyright John Maddock 2006
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)
9 // Simple example of computing probabilities and quantiles for
10 // a Bernoulli random variable representing the flipping of a coin.
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
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.
25 //[binomial_coinflip_example1
27 /*`An example of a [@http://en.wikipedia.org/wiki/Bernoulli_process Bernoulli process]
29 A variable in such a sequence may be called a Bernoulli variable.
31 This example shows using the Binomial distribution to predict the probability
32 of heads and tails when throwing a coin.
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).
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).
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:
46 #include <boost/math/distributions/binomial.hpp>
47 using boost::math::binomial;
50 using std::cout; using std::endl; using std::left;
56 cout << "Using Binomial distribution to predict how many heads and tails." << endl;
60 See note [link coinflip_eg_catch with the catch block]
61 about why a try and catch block is always a good idea.
63 First, construct a binomial distribution with parameters success_fraction
64 1/2, and how many flips.
66 const double success_fraction = 0.5; // = 50% = 1/2 for a 'fair' coin.
68 binomial flip(flips, success_fraction);
72 Then some examples of using Binomial moments (and echoing the parameters).
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.
88 Now we show a variety of predictions on the probability of heads:
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;
94 When we want to calculate the probability for a range or values we can sum the PDF's:
96 cout << "Probability of getting 0 or 1 heads is "
97 << pdf(flip, 0) + pdf(flip, 1) << endl; // sum of exactly == probabilities
99 Or we can use the cdf.
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;
106 cout << "Probability of getting 9 or 10 heads is " << 1. - cdf(flip, 8) << endl;
108 is less accurate than using the complement
110 cout << "Probability of getting 9 or 10 heads is " << cdf(complement(flip, 8)) << endl;
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.
117 To get the probability for a range of heads, we can either add the pdfs for each number of heads
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;
123 But this is probably less efficient than using the cdf
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;
129 Certainly for a bigger range like, 3 to 7
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;
137 Finally, print two tables of probability for the /exactly/ and /at least/ a number of heads.
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;
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;
161 The last (0 to 10 heads) must, of course, be 100% probability.
163 double probability = 0.3;
164 double q = quantile(flip, probability);
165 std::cout << "Quantile (flip, " << probability << ") = " << q << std::endl; // Quantile (flip, 0.3) = 3
167 q = quantile(flip, probability);
168 std::cout << "Quantile (flip, " << probability << ") = " << q << std::endl; // Quantile (flip, 0.6) = 5
170 catch(const std::exception& e)
175 It is always essential to include try & catch blocks because
176 default policies are to throw exceptions on arguments that
177 are out of domain or cause errors like numeric-overflow.
179 Lacking try & catch blocks, the program will abort, whereas the
180 message below from the thrown exception will give some helpful
181 clues as to the cause of the problem.
184 "\n""Message from thrown exception was:\n " << e.what() << std::endl;
186 //] [binomial_coinflip_example1]
192 //[binomial_coinflip_example_output
196 Using Binomial distribution to predict how many heads and tails.
197 From 10 one can expect to get on average 5 heads (or tails).
199 Standard deviation is 1.581
200 So about 2/3 will lie within 1 standard deviation and get between 4 and 6 correct.
202 Skewness if success_fraction is 0.5 is 0
205 Probability of getting no heads is 0.0009766
206 Probability of getting at least one head is 0.999
207 Probability of getting 0 or 1 heads is 0.01074
208 Probability of getting 0 or 1 (<= 1) heads is 0.01074
209 Probability of getting 9 or 10 heads is 0.01074
210 Probability of getting 9 or 10 heads is 0.01074
211 Probability of getting 9 or 10 heads is 0.01074
212 Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6562
213 Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6563
214 Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is 0.8906
216 Probability of getting exactly (==) heads
217 0 0.0009766 or 1 in 1024, or 0.09766%
218 1 0.009766 or 1 in 102.4, or 0.9766%
219 2 0.04395 or 1 in 22.76, or 4.395%
220 3 0.1172 or 1 in 8.533, or 11.72%
221 4 0.2051 or 1 in 4.876, or 20.51%
222 5 0.2461 or 1 in 4.063, or 24.61%
223 6 0.2051 or 1 in 4.876, or 20.51%
224 7 0.1172 or 1 in 8.533, or 11.72%
225 8 0.04395 or 1 in 22.76, or 4.395%
226 9 0.009766 or 1 in 102.4, or 0.9766%
227 10 0.0009766 or 1 in 1024, or 0.09766%
229 Probability of getting upto (<=) heads
230 0 0.0009766 or 1 in 1024, or 0.09766%
231 1 0.01074 or 1 in 93.09, or 1.074%
232 2 0.05469 or 1 in 18.29, or 5.469%
233 3 0.1719 or 1 in 5.818, or 17.19%
234 4 0.377 or 1 in 2.653, or 37.7%
235 5 0.623 or 1 in 1.605, or 62.3%
236 6 0.8281 or 1 in 1.208, or 82.81%
237 7 0.9453 or 1 in 1.058, or 94.53%
238 8 0.9893 or 1 in 1.011, or 98.93%
239 9 0.999 or 1 in 1.001, or 99.9%
240 10 1 or 1 in 1, or 100%
243 //][/binomial_coinflip_example_output]