[Boost-commit] svn:boost r38683 - sandbox/math_toolkit/libs/math/example

From: pbristow_at_[hidden]
Date: 2007-08-15 11:03:12

Author: pbristow
Date: 2007-08-15 11:03:05 EDT (Wed, 15 Aug 2007)
New Revision: 38683
URL: http://svn.boost.org/trac/boost/changeset/38683

Log:
New example referenced from binomial_example.qbk
Added:
   sandbox/math_toolkit/libs/math/example/binomial_coinflip_example.cpp   (contents, props changed)

Added: sandbox/math_toolkit/libs/math/example/binomial_coinflip_example.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/example/binomial_coinflip_example.cpp	2007-08-15 11:03:05 EDT (Wed, 15 Aug 2007)
@@ -0,0 +1,214 @@
+// Copyright Paul A. 2007
+// 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 to some includes and using statements to be able to use the binomial distribution, 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 try'n'catch is 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); // Might be able to calculate an appropriate precision from how many flips?
+/*`
+ 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; // 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;
+/*`
+But this may be inaccurate, so better to use either of
+*/
+    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;
+/*`
+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 << setw(2) << successes << "     " << left << setw(10) 
+        << probability << " or 1 in " << 1. / probability
+        << ", or " << probability * 100. << "%" << endl;
+    } // for i
+    cout << endl;
+
+    // Tabulate the probability of getting between zero and 0 to 10 heads.
+    cout << "Probability of getting upto (<=) heads" << endl;
+    for (int successes = 0; successes <= flips; successes++) 
+    { // Say success means getting a head (or equally success means 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)
+  { // Always useful to include try & catch blocks because
+    // default policies are to throw exceptions on arguments that cause
+    // errors like underflow, overflow. 
+    // Lacking try & catch blocks, the program will abort without a message below,
+    // which may give some helpful clues as to the cause of the exception.
+    std::cout <<
+      "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
+  }
+//] [binomial_coinflip_example1]
+  return 0;
+} // int main()
+
+/*
+
+Output:
+
+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.6563
+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%
+
+*/
+