
BoostCommit : 
From: pbristow_at_[hidden]
Date: 20070815 11:03:12
Author: pbristow
Date: 20070815 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 20070815 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 MathWorldA 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 twosided 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%
+
+*/
+
BoostCommit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk