diff --git a/doc/distributions/binomial_example.qbk b/doc/distributions/binomial_example.qbk
new file mode 100644
index 0000000..7f03e51
--- /dev/null
+++ b/doc/distributions/binomial_example.qbk
@@ -0,0 +1,332 @@
+[section:binom_eg Binomial Distribution Examples]
+
+See also the reference documentation for the __binomial_distrib.
+
+[section:binomial_coinflip_example Binomial Coin-Flipping Example]
+
+[import ../../example/binomial_coinflip_example.cpp]
+[binomial_coinflip_example1]
+
+See [@../../example/binomial_coinflip_example.cpp binomial_coinflip_example.cpp]
+for full source code, the program output looks like this:
+
+[binomial_coinflip_example_output]
+
+[endsect] [/section:binomial_coinflip_example Binomial coinflip example]
+
+[section:binomial_quiz_example Binomial Quiz Example]
+
+[import ../../example/binomial_quiz_example.cpp]
+[binomial_quiz_example1]
+[binomial_quiz_example2]
+[discrete_quantile_real]
+
+See [@../../example/binomial_quiz_example.cpp binomial_quiz_example.cpp]
+for full source code and output.
+
+[endsect] [/section:binomial_coinflip_quiz Binomial Coin-Flipping example]
+
+[section:binom_conf Calculating Confidence Limits on the Frequency of Occurrence for a Binomial Distribution]
+
+Imagine you have a process that follows a binomial distribution: for each
+trial conducted, an event either occurs or does it does not, referred
+to as "successes" and "failures".  If, by experiment, you want to measure the
+frequency with which successes occur, the best estimate is given simply
+by /k/ \/ /N/, for /k/ successes out of /N/ trials.  However our confidence in that
+estimate will be shaped by how many trials were conducted, and how many successes
+were observed.  The static member functions 
+`binomial_distribution<>::find_lower_bound_on_p` and
+`binomial_distribution<>::find_upper_bound_on_p` allow you to calculate
+the confidence intervals for your estimate of the occurrence frequency.
+
+The sample program [@../../example/binomial_confidence_limits.cpp 
+binomial_confidence_limits.cpp] illustrates their use.  It begins by defining
+a procedure that will print a table of confidence limits for various degrees
+of certainty:
+
+   #include <iostream>
+   #include <iomanip>
+   #include <boost/math/distributions/binomial.hpp>
+
+   void confidence_limits_on_frequency(unsigned trials, unsigned successes)
+   {
+      //
+      // trials = Total number of trials.
+      // successes = Total number of observed successes.
+      //
+      // Calculate confidence limits for an observed
+      // frequency of occurrence that follows a binomial
+      // distribution.
+      //
+      using namespace std;
+      using namespace boost::math;
+
+      // Print out general info:
+      cout <<
+         "___________________________________________\n"
+         "2-Sided Confidence Limits For Success Ratio\n"
+         "___________________________________________\n\n";
+      cout << setprecision(7);
+      cout << setw(40) << left << "Number of Observations" << "=  " << trials << "\n";
+      cout << setw(40) << left << "Number of successes" << "=  " << successes << "\n";
+      cout << setw(40) << left << "Sample frequency of occurrence" << "=  " << double(successes) / trials << "\n";
+
+The procedure now defines a table of significance levels: these are the 
+probabilities that the true occurrence frequency lies outside the calculated
+interval:
+
+      double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+
+Some pretty printing of the table header follows:
+
+   cout << "\n\n"
+           "_______________________________________________________________________\n"
+           "Confidence        Lower CP       Upper CP       Lower JP       Upper JP\n"
+           " Value (%)        Limit          Limit          Limit          Limit\n"
+           "_______________________________________________________________________\n";
+
+
+And now for the important part - the intervals themselves - for each
+value of /alpha/, we call `find_lower_bound_on_p` and 
+`find_lower_upper_on_p` to obtain lower and upper bounds
+respectively.  Note that since we are calculating a two-sided interval,
+we must divide the value of alpha in two.
+
+Please note that calculating two separate /single sided bounds/, each with risk
+level [alpha][space]is not the same thing as calculating a two sided interval.
+Had we calculate two single-sided intervals each with a risk
+that the true value is outside the interval of [alpha], then:
+
+* The risk that it is less than the lower bound is [alpha].
+
+and
+
+* The risk that it is greater than the upper bound is also [alpha].
+
+So the risk it is outside *upper or lower bound*, is *twice* alpha, and the 
+probability that it is inside the bounds is therefore not nearly as high as 
+one might have thought.  This is why [alpha]/2 must be used in 
+the calculations below.
+
+In contrast, had we been calculating a 
+single-sided interval, for example: ['"Calculate a lower bound so that we are P%
+sure that the true occurrence frequency is greater than some value"]
+then we would *not* have divided by two.
+
+Finally note that `binomial_distribution` provides a choice of two
+methods for the calculation, we print out the results from both 
+methods in this example:
+
+      for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+      {
+         // Confidence value:
+         cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]);
+         // Calculate Clopper Pearson bounds:
+         double l = binomial_distribution<>::find_lower_bound_on_p(
+                        trials, successes, alpha[i]/2);
+         double u = binomial_distribution<>::find_upper_bound_on_p(
+                        trials, successes, alpha[i]/2);
+         // Print Clopper Pearson Limits:
+         cout << fixed << setprecision(5) << setw(15) << right << l;
+         cout << fixed << setprecision(5) << setw(15) << right << u;
+         // Calculate Jeffreys Prior Bounds:
+         l = binomial_distribution<>::find_lower_bound_on_p(
+               trials, successes, alpha[i]/2, 
+               binomial_distribution<>::jeffreys_prior_interval);
+         u = binomial_distribution<>::find_upper_bound_on_p(
+               trials, successes, alpha[i]/2, 
+               binomial_distribution<>::jeffreys_prior_interval);
+         // Print Jeffreys Prior Limits:
+         cout << fixed << setprecision(5) << setw(15) << right << l;
+         cout << fixed << setprecision(5) << setw(15) << right << u << std::endl;
+      }
+      cout << endl;
+   }
+
+And that's all there is to it.  Let's see some sample output for a 2 in 10
+success ratio, first for 20 trials:
+
+[pre'''___________________________________________
+2-Sided Confidence Limits For Success Ratio
+___________________________________________
+
+Number of Observations                  =  20
+Number of successes                     =  4
+Sample frequency of occurrence          =  0.2
+
+
+_______________________________________________________________________
+Confidence        Lower CP       Upper CP       Lower JP       Upper JP
+ Value (%)        Limit          Limit          Limit          Limit
+_______________________________________________________________________
+    50.000        0.12840        0.29588        0.14974        0.26916
+    75.000        0.09775        0.34633        0.11653        0.31861
+    90.000        0.07135        0.40103        0.08734        0.37274
+    95.000        0.05733        0.43661        0.07152        0.40823
+    99.000        0.03576        0.50661        0.04655        0.47859
+    99.900        0.01905        0.58632        0.02634        0.55960
+    99.990        0.01042        0.64997        0.01530        0.62495
+    99.999        0.00577        0.70216        0.00901        0.67897
+''']
+
+As you can see, even at the 95% confidence level the bounds are
+really quite wide (this example is chosen to be easily compared to the one
+in the __handbook
+[@http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
+here]).  Note also that the Clopper-Pearson calculation method (CP above)
+produces quite noticeably more pessimistic estimates than the Jeffreys Prior
+method (JP above).
+
+
+Compare that with the program output for
+2000 trials:
+
+[pre'''___________________________________________
+2-Sided Confidence Limits For Success Ratio
+___________________________________________
+
+Number of Observations                  =  2000
+Number of successes                     =  400
+Sample frequency of occurrence          =  0.2000000
+
+
+_______________________________________________________________________
+Confidence        Lower CP       Upper CP       Lower JP       Upper JP
+ Value (%)        Limit          Limit          Limit          Limit
+_______________________________________________________________________
+    50.000        0.19382        0.20638        0.19406        0.20613
+    75.000        0.18965        0.21072        0.18990        0.21047
+    90.000        0.18537        0.21528        0.18561        0.21503
+    95.000        0.18267        0.21821        0.18291        0.21796
+    99.000        0.17745        0.22400        0.17769        0.22374
+    99.900        0.17150        0.23079        0.17173        0.23053
+    99.990        0.16658        0.23657        0.16681        0.23631
+    99.999        0.16233        0.24169        0.16256        0.24143
+''']
+
+Now even when the confidence level is very high, the limits are really
+quite close to the experimentally calculated value of 0.2.  Furthermore
+the difference between the two calculation methods is now really quite small.
+
+[endsect]
+
+[section:binom_size_eg Estimating Sample Sizes for a Binomial Distribution.]
+
+Imagine you have a critical component that you know will fail in 1 in
+N "uses" (for some suitable definition of "use").  You may want to schedule
+routine replacement of the component so that its chance of failure between
+routine replacements is less than P%.  If the failures follow a binomial
+distribution (each time the component is "used" it either fails or does not)
+then the static member function `binomial_distibution<>::find_maximum_number_of_trials`
+can be used to estimate the maximum number of "uses" of that component for some
+acceptable risk level /alpha/.
+
+The example program 
+[@../../example/binomial_sample_sizes.cpp binomial_sample_sizes.cpp]
+demonstrates its usage.  It centres on a routine that prints out
+a table of maximum sample sizes for various probability thresholds:
+
+   void find_max_sample_size(
+      double p,              // success ratio.
+      unsigned successes)    // Total number of observed successes permitted.
+   {
+
+The routine then declares a table of probability thresholds: these are the
+maximum acceptable probability that /successes/ or fewer events will be
+observed.  In our example, /successes/ will be always zero, since we want
+no component failures, but in other situations non-zero values may well
+make sense.
+
+   double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+
+Much of the rest of the program is pretty-printing, the important part
+is in the calculation of maximum number of permitted trials for each
+value of alpha:
+
+   for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+   {
+      // Confidence value:
+      cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]);
+      // calculate trials:
+      double t = binomial::find_maximum_number_of_trials(
+                     successes, p, alpha[i]);
+      t = floor(t);
+      // Print Trials:
+      cout << fixed << setprecision(5) << setw(15) << right << t << endl;
+   }
+
+Note that since we're
+calculating the maximum number of trials permitted, we'll err on the safe
+side and take the floor of the result.  Had we been calculating the
+/minimum/ number of trials required to observe a certain number of /successes/
+using `find_minimum_number_of_trials` we would have taken the ceiling instead.
+
+We'll finish off by looking at some sample output, firstly for
+a 1 in 1000 chance of component failure with each use:
+
+[pre
+'''________________________
+Maximum Number of Trials
+________________________
+
+Success ratio                           =  0.001
+Maximum Number of "successes" permitted =  0
+
+
+____________________________
+Confidence        Max Number
+ Value (%)        Of Trials
+____________________________
+    50.000            692
+    75.000            287
+    90.000            105
+    95.000             51
+    99.000             10
+    99.900              0
+    99.990              0
+    99.999              0'''
+]
+
+So 51 "uses" of the component would yield a 95% chance that no
+component failures would be observed.
+
+Compare that with a 1 in 1 million chance of component failure:
+
+[pre'''
+________________________
+Maximum Number of Trials
+________________________
+
+Success ratio                           =  0.0000010
+Maximum Number of "successes" permitted =  0
+
+
+____________________________
+Confidence        Max Number
+ Value (%)        Of Trials
+____________________________
+    50.000         693146
+    75.000         287681
+    90.000         105360
+    95.000          51293
+    99.000          10050
+    99.900           1000
+    99.990            100
+    99.999             10'''
+]
+
+In this case, even 1000 uses of the component would still yield a 
+less than 1 in 1000 chance of observing a component failure 
+(i.e. a 99.9% chance of no failure).
+
+[endsect] [/section:binom_size_eg Estimating Sample Sizes for a Binomial Distribution.]
+
+[endsect][/section:binom_eg Binomial Distribution]
+
+[/ 
+  Copyright 2006 John Maddock and Paul A. Bristow.
+  Distributed under 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).
+]
+