summaryrefslogtreecommitdiffstats
path: root/src/boost/libs/math/example/binomial_coinflip_example.cpp
blob: 5f1826b94a5ff4b5c2f5b2e25f183ce2a192757c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
// 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 up to 10 heads.
    cout << "Probability of getting up to (<=) 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 up to (<=) 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]