Binomial Coin-Flipping Example

| Home | Libraries | People | FAQ | More

###### Binomial Coin-Flipping Example

An example of a 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 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 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)
    {
      //

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;
    }

See binomial_coinflip_example.cpp for full source code, the program output looks like this:

    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%