]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | */ | |
92f5a8d4 TL |
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 | |
166 | probability = 0.6; | |
167 | q = quantile(flip, probability); | |
168 | std::cout << "Quantile (flip, " << probability << ") = " << q << std::endl; // Quantile (flip, 0.6) = 5 | |
7c673cae FG |
169 | } |
170 | catch(const std::exception& e) | |
171 | { | |
172 | // | |
173 | /*` | |
174 | [#coinflip_eg_catch] | |
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. | |
178 | ||
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. | |
182 | */ | |
183 | std::cout << | |
184 | "\n""Message from thrown exception was:\n " << e.what() << std::endl; | |
185 | } | |
186 | //] [binomial_coinflip_example1] | |
187 | return 0; | |
188 | } // int main() | |
189 | ||
190 | // Output: | |
191 | ||
192 | //[binomial_coinflip_example_output | |
193 | /*` | |
194 | ||
195 | [pre | |
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). | |
198 | Mode is 5 | |
199 | Standard deviation is 1.581 | |
200 | So about 2/3 will lie within 1 standard deviation and get between 4 and 6 correct. | |
201 | Skewness is 0 | |
202 | Skewness if success_fraction is 0.5 is 0 | |
203 | ||
204 | For 10 coin flips: | |
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 | |
215 | ||
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% | |
228 | ||
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% | |
241 | ] | |
242 | */ | |
243 | //][/binomial_coinflip_example_output] |