
BoostCommit : 
From: pbristow_at_[hidden]
Date: 20070814 04:41:15
Author: pbristow
Date: 20070814 04:41:13 EDT (Tue, 14 Aug 2007)
New Revision: 38646
URL: http://svn.boost.org/trac/boost/changeset/38646
Log:
update for new example using policies
Text files modified:
sandbox/math_toolkit/libs/math/example/negative_binomial_construction_examples.cpp  64 +++
sandbox/math_toolkit/libs/math/example/negative_binomial_example1.cpp  399 ++++++++++
sandbox/math_toolkit/libs/math/example/negative_binomial_example2.cpp  225 +++++++++++
sandbox/math_toolkit/libs/math/example/negative_binomial_example3.cpp  182 ++++++++++++++++
4 files changed, 413 insertions(+), 457 deletions()
Modified: sandbox/math_toolkit/libs/math/example/negative_binomial_construction_examples.cpp
==============================================================================
 sandbox/math_toolkit/libs/math/example/negative_binomial_construction_examples.cpp (original)
+++ sandbox/math_toolkit/libs/math/example/negative_binomial_construction_examples.cpp 20070814 04:41:13 EDT (Tue, 14 Aug 2007)
@@ 1,3 +1,6 @@
+NOte now obselete  see deistribution_construction.cpp
+
+
// negative_binomial_example2.cpp
// Copyright Paul A. Bristow 2007.
@@ 7,7 +10,7 @@
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
// Example 2 of using constructing distributions, mainly negative_binomial.
+// Example 2 of using constructing distributions, shown here for negative_binomial.
#include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
using boost::math::negative_binomial_distribution; // default type is double.
@@ 20,51 +23,25 @@
int main()
{
 cout << "Example 2 constructing Distributions (Negative_binomial).";
 #if defined(__FILE__) && defined(__TIMESTAMP__)
 cout << " " << __FILE__ << ' ' << __TIMESTAMP__ << ' '<< _MSC_FULL_VER << "\n";
 #endif
 cout << endl;

// Several examples of constructing distributions, for example, negative binomial:
+ // A negative binomial with 8 successes and a success fraction 0.25, 25% or 1 in 4 is constructed like this:
 // Fundamentally constructed like this:
boost::math::negative_binomial_distribution<double> mydist0(8., 0.25);
// But this is inconveniently long.
+ // The prefix boost::math:: can be avoided by
using boost::math::negative_binomial_distribution;
// Allows convenient reference to negative_binomial_distribution.
 // You can provide the type explicitly thus:
 negative_binomial_distribution<double> mydist1(8., 0.25); // Explicit double.
 negative_binomial_distribution<float> mydist2(8., 0.25); // Explicit float, double arguments > float.
 negative_binomial_distribution<float> mydist3(8, 0.25); // Explicit float, integer & double arguments > float.
 negative_binomial_distribution<float> mydist4(8.F, 0.25F); // Explicit float, float arguments, no conversion.
 negative_binomial_distribution<float> mydist5(8, 1); // Explicit integer, integer arguments > float.
 negative_binomial_distribution<double> mydist6(8., 0.25); // Explicit double.
 negative_binomial_distribution<long double> mydist7(8., 0.25); // Explicit long double.
 // And if you have your own RealType called MyFPType,
 // for example NTL quad_float (128bit floatingpoint), then:
 // negative_binomial_distribution<MyFPType> mydist6(8, 1); // Integer arguments > MyFPType.

 // Note that default constructor arguments are only provided for some distributions.
 // negative_binomial_distribution<> mydist;
 // error C2512 no appropriate default constructor available.
 // Since there are no accessor functions, no default constructor are provided,
 // because it is difficult to chose any sensible default values for this distribution.
 // For other distribution, like the normal distribution,
 // it is obviously very useful to provide
 // defaults for the mean and standard deviation thus:
 // normal_distribution(RealType mean = 0, RealType sd = 1)

negative_binomial_distribution<> mydist9(8., 0.25); // Uses default RealType = double.
// But the name "negative_binomial_distribution" is still inconveniently long,
// so for most distributions, a typedef is provided, for example:
// typedef negative_binomial_distribution<double> negative_binomial; // Reserved name of type double.
+ using boost::math::negative_binomial; // Allows convenient access to the name negative_binomial.
+
// Some examples using the provided typedef:
 using boost::math::negative_binomial; // Convenient access to the name.
// Allows convenient reference to negative_binomial of default type double.
negative_binomial mydist10(5., 0.4); // Both arguments double.
// And automatic conversion takes place, so you can use integers and floats:
@@ 83,14 +60,33 @@
// NOT binomial myb010(1, 0.5); but
binomial_distribution<> myb1(1, 0.5);
+ // You can also provide the type RealType explicitly thus:
+ negative_binomial_distribution<double> mydist1(8., 0.25); // Explicit double.
+ negative_binomial_distribution<float> mydist2(8., 0.25); // Explicit float, double arguments > float.
+ negative_binomial_distribution<float> mydist3(8, 0.25); // Explicit float, integer & double arguments > float.
+ negative_binomial_distribution<float> mydist4(8.F, 0.25F); // Explicit float, float arguments, no conversion.
+ negative_binomial_distribution<float> mydist5(8, 1); // Explicit integer, integer arguments > float.
+ negative_binomial_distribution<double> mydist6(8., 0.25); // Explicit double.
+ negative_binomial_distribution<long double> mydist7(8., 0.25); // Explicit long double.
+ // And if you have your own RealType called MyFPType,
+ // for example NTL quad_float (128bit floatingpoint), then:
+ // negative_binomial_distribution<MyFPType> mydist6(8, 1); // Integer arguments > MyFPType.
+
+ // Note that default constructor arguments are only provided for some distributions.
+ // negative_binomial_distribution<> mydist;
+ // error C2512 no appropriate default constructor available.
+ // Since there are no accessor functions, no default constructor are provided,
+ // because it is difficult to chose any sensible default values for this distribution.
+ // For other distribution, like the normal distribution,
+ // it is obviously very useful to provide
+ // defaults for the mean and standard deviation thus:
+ // normal_distribution(RealType mean = 0, RealType sd = 1)
+
return 0;
} // int main()
/*
Output is:

math_toolkit\libs\math\example\negative_binomial_construction_examples.cpp Wed Aug 1 13:59:34 2007 140050727
*/
Modified: sandbox/math_toolkit/libs/math/example/negative_binomial_example1.cpp
==============================================================================
 sandbox/math_toolkit/libs/math/example/negative_binomial_example1.cpp (original)
+++ sandbox/math_toolkit/libs/math/example/negative_binomial_example1.cpp 20070814 04:41:13 EDT (Tue, 14 Aug 2007)
@@ 9,69 +9,33 @@
// Example 1 of using negative_binomial distribution.
// http://en.wikipedia.org/wiki/Negative_binomial_distribution

// (After a problem by Dr. Diane Evans,
// Professor of mathematics at RoseHulman Institute of Technology)

// Pat is required to sell candy bars to raise money for the 6th grade field trip.
// There are thirty houses in the neighborhood,
// and Pat is not supposed to return home until five candy bars have been sold.
// So the child goes door to door, selling candy bars.
// At each house, there is a 0.4 probability (40%) of selling one candy bar
// and a 0.6 probability (60%) of selling nothing.

// What is the probability mass/density function for selling the last (fifth) candy bar at the nth house?

// The Negative Binomial(r, p) distribution describes the probability of k failures
// and r successes in k+r Bernoulli(p) trials with success on the last trial.
// Selling five candy bars means getting five successes, so successes r = 5.
// The total number of trials n (in this case, houses) this takes is therefore
// = sucesses + failures or k + r = k + 5.
// The random variable we are interested in is the number of houses k
// that must be visited to sell five candy bars,
// so we substitute k = n 5 into a NegBin(5, 0.4) mass/density function
// and obtain the following mass/density function of the distribution of houses (for n >= 5):
// Obviously, the best case is that Pat makes sales on all the first five houses.

// What is the probability that Pat finishes ON the tenth house?

// f(10) = 0.1003290624, or about 1 in 10

// What is the probability that Pat finishes ON OR BEFORE reaching the eighth house?

// To finish on or before the eighth house,
// Pat must finish at the fifth, sixth, seventh, or eighth house.
// Sum those probabilities:

 // f(5) = 0.01024
 // f(6) = 0.03072
 // f(7) = 0.055296
 // f(8) = 0.0774144
 // sum {j=5 to 8} f(j) = 0.17367

// What is the probability that Pat exhausts all 30 houses in the neighborhood,
// and still doesn't sell the required 5 candy bars?

// 1  sum{j=5 to 30} f(j) = 1  incomplete beta (p = 0.4)(5, 305+1) =~ 1  0.99849 = 0.00151 = 0.15%.

// see also http://www.math.uah.edu/stat/bernoulli/Introduction.xhtml
// http://www.codecogs.com/pages/catgen.php?category=stats/dists/discrete
+//[negative_binomial_eg1_1
+/*`
+First we need to #define macros to control the error and discrete handling policies.
+For this simple example, we want to avoid throwing an exception (the default policy) and just return infinity.
+We want to treat the distribution as if it was continuous, so we choose a policy of real,
+rather than the default policy integer_outside.
+*/
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
// Avoid throwing an exception (the default policy) and just return infinity.
+#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+
+/*`After that (vital to #include *after* #defines), we need some includes
+to provide easy access to the negative binomial distribution,
+and some std library iostream, of course.
+*/
#include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
 using boost::math::negative_binomial_distribution;
+#include <boost/math/distributions/negative_binomial.hpp>
+// for negative_binomial_distribution
using boost::math::negative_binomial; // typedef provides default type is double.
using ::boost::math::cdf;
 using ::boost::math::pdf; // Probability of negative_binomial.
+ using ::boost::math::pdf; // Probability mass of negative_binomial.
using ::boost::math::quantile;
#include <iostream>
 using std::cout;
 using std::endl;
 using std::noshowpoint;
+ using std::cout;
+ using std::endl;
+ using std::noshowpoint;
using std::fixed;
using std::right;
using std::left;
@@ 79,44 +43,57 @@
using std::setprecision;
using std::setw;
+#include <limits>
+ using std::numeric_limits;
+
int main()
{
 cout << "Example 1 using the Negative Binomial Distribution.";
+ cout << "Example 1 using the Negative Binomial Distribution.";
#if defined(__FILE__) && defined(__TIMESTAMP__)
 cout << " " << __FILE__ << ' ' << __TIMESTAMP__ << ' '<< _MSC_FULL_VER << "\n";
+ cout << " " << __FILE__ << ' ' << __TIMESTAMP__ << "\n";
#endif
 cout << endl;
 cout.precision(5); // NB INF shows wrongly with < 5 !
+ cout << endl;
+ cout.precision(5);
+ // None of the values calculated have a useful accuracy as great this, but
+ // INF shows wrongly with < 5 !
// https://connect.microsoft.com/VisualStudio/feedback/ViewFeedback.aspx?FeedbackID=240227

try
{
double sales_quota = 5; // Pat's sales quota  successes (r).
 double success_fraction = 0.4; // success_fraction (p)  so fail_fraction is 0.6.
+ double success_fraction = 0.4; // success_fraction (p)  so failure_fraction is 0.6.
negative_binomial nb(sales_quota, success_fraction); // double by default.
int all_houses = 30; // The number of houses on the estate.
 cout <<"Selling candy bars  an example of using the negative binomial distribution. "
 << "\n""An example by Dr. Diane Evans,"
+ cout <<"Selling candy bars  using the negative binomial distribution."
+ << "\n""by Dr. Diane Evans,"
"\n""Professor of Mathematics at RoseHulman Institute of Technology,"
<< "\n""see http://en.wikipedia.org/wiki/Negative_binomial_distribution""\n"
<< endl;
cout << "Pat has a sales per house success rate of " << success_fraction
 << ".""\n""Therefore he would, on average, sell " << nb.success_fraction() * 100 << " bars after trying 100 houses." << endl;
+ << ".""\n""Therefore he would, on average, sell " << nb.success_fraction() * 100
+ << " bars after trying 100 houses." << endl;
 cout << "With a success rate of " << nb.success_fraction() << ", he might expect, on average,""\n"
 " to need to visit about " << success_fraction * all_houses << " houses in order to sell all " << nb.successes() << " candy bars. " << endl;

 // To finish on or before the 8th house, Pat must finish at the 5th, 6th, 7th or 8th house.
 // (Obviously he could not finish on fewer than 5 houses because he must sell 5 candy bars.
 // so the 5th house is the first that he could possibly finish on).
 // The probability that he will finish on EXACTLY on any house is the Probability Density Function (pdf).
 cout << "Probability that Pat finishes on the " << sales_quota << "th house is " << "f(5) = " << pdf(nb, nb.successes()) << endl;
 cout << "Probability that Pat finishes on the 6th house is " << pdf(nb, 6  sales_quota) << endl;
 cout << "Probability that Pat finishes on the 7th house is " << pdf(nb, 7  sales_quota) << endl;
 cout << "Probability that Pat finishes on the 8th house is " << pdf(nb, 8  sales_quota) << endl;
+ cout << "With a success rate of " << nb.success_fraction()
+ << ", he might expect, on average,""\n"
+ " to need to visit about " << success_fraction * all_houses
+ << " houses in order to sell all " << nb.successes() << " candy bars. " << endl;
+/*`
+To finish on or before the 8th house, Pat must finish at the 5th, 6th, 7th or 8th house.
+(Obviously he could not finish on fewer than 5 houses because he must sell 5 candy bars.
+So the 5th house is the first that he could possibly finish on).
+The probability that he will finish on EXACTLY on any house
+is the Probability Density Function (pdf).
+*/
+ cout << "Probability that Pat finishes on the " << sales_quota << "th house is "
+ << "f(5) = " << pdf(nb, nb.successes()) << endl;
+ cout << "Probability that Pat finishes on the 6th house is "
+ << pdf(nb, 6  sales_quota) << endl;
+ cout << "Probability that Pat finishes on the 7th house is "
+ << pdf(nb, 7  sales_quota) << endl;
+ cout << "Probability that Pat finishes on the 8th house is "
+ << pdf(nb, 8  sales_quota) << endl;
 // The sum of the probabilities for these houses is the Cumulative Distribution Function (cdf).
+ // Sum of the probabilities for these houses is the Cumulative Distribution Function (cdf).
cout << "Probability that Pat finishes on or before the 8th house is sum "
"\n" << "pdf(sales_quota) + pdf(6) + pdf(7) + pdf(8) = "
// Sum each of the mass/density probabilities for houses sales_quota = 5, 6, 7, & 8.
@@ 126,28 +103,35 @@
+ pdf(nb, 8  sales_quota) // 3
<< endl;
 // Or using the negative binomial **cumulative** distribution function (cdf instead sum of the pdfs):
+ // Or using the negative binomial *cumulative* distribution function
+ // (cdf instead sum of the pdfs):
cout << "\n""Probability of selling his quota of " << sales_quota
<< " candy bars""\n""on or before the " << 8 << "th house is "
<< cdf(nb, 8  sales_quota) << endl;
 cout << "\n""Probability that Pat finishes exactly on the 10th house is " << pdf(nb, 10  sales_quota) << endl;
+ cout << "\n""Probability that Pat finishes exactly on the 10th house is "
+ << pdf(nb, 10  sales_quota) << endl;
cout << "\n""Probability of selling his quota of " << sales_quota
<< " candy bars""\n""on or before the " << 10 << "th house is "
<< cdf(nb, 10  sales_quota) << endl;
 cout << "Probability that Pat finishes on the 11th house is " << pdf(nb, 11  sales_quota) << endl;
+ cout << "Probability that Pat finishes on the 11th house is "
+ << pdf(nb, 11  sales_quota) << endl;
cout << "\n""Probability of selling his quota of " << sales_quota
<< " candy bars""\n""on or before the " << 11 << "th house is "
<< cdf(nb, 11  sales_quota) << endl;
 cout << "Probability that Pat finishes on the 12th house is " << pdf(nb, 12  sales_quota) << endl;
+ cout << "Probability that Pat finishes on the 12th house is "
+ << pdf(nb, 12  sales_quota) << endl;
cout << "\n""Probability of selling his quota of " << sales_quota
<< " candy bars""\n""on or before the " << 12 << "th house is "
<< cdf(nb, 12  sales_quota) << endl;
 // Finally consider the risk of Pat not setting his quota of 5 bars even after visiting all the houses.
 // Calculate the probability that he would sell on the (nonexistent) lastplus1 house.
+/*`
+Finally consider the risk of Pat not selling his quota of 5 bars
+even after visiting all the houses.
+Calculate the probability that he would sell on the (nonexistent) lastplus1 house.
+*/
cout << "\n""Probability of selling his quota of " << sales_quota
<< " candy bars""\n""on or before the " << all_houses + 1 << "th house is "
<< cdf(nb, all_houses + 1  sales_quota) << endl;
@@ 161,11 +145,15 @@
// Probability of meeting sales quota on or before 8th house is 0.174
cout << "If the confidence of meeting sales quota is " << p
<< ", then the finishing house is " << quantile(nb, p) + sales_quota << endl;
 // Also try wanting absolute certainty that all 5 will be sold
 // which implies an infinite number of sales.
 // (Of course, there are only 30 houses on the estate, so he can't even be certain of selling his quota.
+/*`
+Also try demanding absolute certainty that all 5 will be sold,
+which implies an infinite number of trials.
+(Of course, there are only 30 houses on the estate,
+so he can't even be certain of selling his quota).
+*/
cout << "If the confidence of meeting sales quota is " << 1.
 << ", then the finishing house is " << quantile(nb, 1) + sales_quota << endl; // 1.#INF == infinity.
+ << ", then the finishing house is " << quantile(nb, 1) + sales_quota << endl;
+ // 1.#INF == infinity.
cout << "If the confidence of meeting sales quota is " << 0.
<< ", then the finishing house is " << quantile(nb, 0.) + sales_quota << endl;
@@ 176,14 +164,15 @@
cout << "If the confidence of meeting sales quota is " << 1  0.00151 // 30 th
<< ", then the finishing house is " << quantile(nb, 1  0.00151) + sales_quota << endl;
 // If the opposite is true, we don't want to assume any confidence, then
 // this is tantamount to assuming that all the first sales_quota will be successful sales.
+/*
+If the opposite is true, we don't want to assume any confidence, then this is tantamount
+to assuming that all the first sales_quota trials will be successful sales.
+*/
cout << "If confidence of meeting quota is zero (we assume all houses are successful sales)"
", then finishing house is " << sales_quota << endl;

double ps[] = {0., 0.001, 0.01, 0.05, 0.1, 0.5, 0.9, 0.95, 0.99, 0.999, 1.};
 // Confidence as fraction = 1alpha, as percent = 100 * (1alpha[i])
+ // Confidence as fraction = 1alpha, as percent = 100 * (1alpha[i]) %
for (int i = 0; i < sizeof(ps)/sizeof(ps[0]); i++)
{
cout << "If confidence of meeting quota is " << ps[i]
@@ 199,190 +188,24 @@
cout << left << setw(6) << cdf(nb, i  sales_quota) << " " << setw(3) << i<< endl;
}
cout << endl;
 // Define a table of significance levels:
 double alpha[] = { 1., 0.99999, 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 , 0.000000001};

 // Now turn the question on its head as ask how many houses for a given confidence.

 cout << "\n"
 "Confidence (%) Minimum houses for " << sales_quota << " sales" << endl;
 for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
 { // Confidence values %:
 cout << fixed << setprecision(3) << setw(10) << right << 100 * (1alpha[i]) << " "
 // find_minimum_number_of_trials
 << setw(6) << right
 << (int)ceil(negative_binomial::find_minimum_number_of_trials(sales_quota, success_fraction, alpha[i]))
 << endl;
 }
 cout << endl;

 /*
 Notes:
Confidence (%) Minimum houses for 5 sales
 0.000 5 only if confidence is really zero can we assume that all visits will be successful.
 0.001 6 So even a tiny increase means we must allow for one failure to sell.
 50.000 10

 But above using quantile(nb, 0.5) + sales_quota

" If the confidence of meeting sales quota is 0.5, then the finishing house is 12 "

Probability of selling his quota of 5 candy bars
on or before the 10th house is 0.3669

Probability that Pat finishes on the 11th house is 0.10033
Probability of selling his quota of 5 candy bars
on or before the 11th house is 0.46723

Probability that Pat finishes on the 12th house is 0.094596
Probability of selling his quota of 5 candy bars
on or before the 12th house is 0.56182

So is this formula wrong / over optimistic?

 75.000 11
 90.000 13
 95.000 15
 99.000 18
 99.900 21
 99.990 25
 99.999 28

 */


}
catch(const std::exception& e)
 {
 std::cout <<
 "\n""Message from thrown exception was:\n " << e.what() << std::endl;
 }

+ { // Since we have set an overflow policy of ignore_error,
+ // an exception should never be thrown.
+ std::cout << "\n""Message from thrown exception was:\n " << e.what() << std::endl;
+ }
return 0;
} // int main()
+//] [/ negative_binomial_eg1_1]
+
/*
Output is:
Example 1 using the Negative Binomial Distribution. ..\..\..\..\..\..\Boostsandbox\math_toolkit\libs\math\example\negative_binomial_example1.cpp Wed Aug 1 14:12:58 2007 140050727
Selling candy bars  an example of using the negative binomial distribution.
An example by Dr. Diane Evans,
Professor of Mathematics at RoseHulman Institute of Technology,
see http://en.wikipedia.org/wiki/Negative_binomial_distribution
Pat has a sales per house success rate of 0.4.
Therefore he would, on average, sell 40 bars after trying 100 houses.
With a success rate of 0.4, he might expect, on average,
 to need to visit about 12 houses in order to sell all 5 candy bars.
Probability that Pat finishes on the 5th house is f(5) = 0.10033
Probability that Pat finishes on the 6th house is 0.03072
Probability that Pat finishes on the 7th house is 0.055296
Probability that Pat finishes on the 8th house is 0.077414
Probability that Pat finishes on or before the 8th house is sum
pdf(sales_quota) + pdf(6) + pdf(7) + pdf(8) = 0.17367
Probability of selling his quota of 5 candy bars
on or before the 8th house is 0.17367
Probability that Pat finishes exactly on the 10th house is 0.10033
Probability of selling his quota of 5 candy bars
on or before the 10th house is 0.3669
Probability that Pat finishes on the 11th house is 0.10033
Probability of selling his quota of 5 candy bars
on or before the 11th house is 0.46723
Probability that Pat finishes on the 12th house is 0.094596
Probability of selling his quota of 5 candy bars
on or before the 12th house is 0.56182
Probability of selling his quota of 5 candy bars
on or before the 31th house is 0.99897
Probability of failing to sell his quota of 5 candy bars
even after visiting all 30 houses is 0.0010314
Probability of meeting sales quota on or before 8th house is 0.17367
If the confidence of meeting sales quota is 0.17367, then the finishing house is 8
Message from thrown exception was:
 Error in function boost::math::quantile(const negative_binomial_distribution<double>&, double): Probability argument is 1, which implies infinite failures !

After adding

#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error

to get the previous overflow policy

Compiling...
negative_binomial_example1.cpp
Linking...
Autorun "i:\boost0605031300\libs\math\test\Math_test\debug\negative_binomial_example1.exe"
Autorun "i:\boost0605031300\libs\math\test\Math_test\debug\negative_binomial_example1.exe"
Example 1 using the Negative Binomial Distribution. ..\..\..\..\..\..\boostsandbox\math_toolkit\libs\math\example\negative_binomial_example1.cpp Mon Aug 6 12:57:09 2007 140050727
Selling candy bars  an example of using the negative binomial distribution.
An example by Dr. Diane Evans,
Professor of Mathematics at RoseHulman Institute of Technology,
see http://en.wikipedia.org/wiki/Negative_binomial_distribution
Pat has a sales per house success rate of 0.4.
Therefore he would, on average, sell 40 bars after trying 100 houses.
With a success rate of 0.4, he might expect, on average,
 to need to visit about 12 houses in order to sell all 5 candy bars.
Probability that Pat finishes on the 5th house is f(5) = 0.10033
Probability that Pat finishes on the 6th house is 0.03072
Probability that Pat finishes on the 7th house is 0.055296
Probability that Pat finishes on the 8th house is 0.077414
Probability that Pat finishes on or before the 8th house is sum
pdf(sales_quota) + pdf(6) + pdf(7) + pdf(8) = 0.17367
Probability of selling his quota of 5 candy bars
on or before the 8th house is 0.17367
Probability that Pat finishes exactly on the 10th house is 0.10033
Probability of selling his quota of 5 candy bars
on or before the 10th house is 0.3669
Probability that Pat finishes on the 11th house is 0.10033
Probability of selling his quota of 5 candy bars
on or before the 11th house is 0.46723
Probability that Pat finishes on the 12th house is 0.094596
Probability of selling his quota of 5 candy bars
on or before the 12th house is 0.56182
Probability of selling his quota of 5 candy bars
on or before the 31th house is 0.99897
Probability of failing to sell his quota of 5 candy bars
even after visiting all 30 houses is 0.0010314
Probability of meeting sales quota on or before 8th house is 0.17367
If the confidence of meeting sales quota is 0.17367, then the finishing house is 8
If the confidence of meeting sales quota is 1, then the finishing house is 1.#INF
If the confidence of meeting sales quota is 0, then the finishing house is 5
If the confidence of meeting sales quota is 0.5, then the finishing house is 12
If the confidence of meeting sales quota is 0.99849, then the finishing house is 31
If confidence of meeting quota is zero (we assume all houses are successful sales), then finishing house is 5
If confidence of meeting quota is 0, then finishing house is 5
If confidence of meeting quota is 0.001, then finishing house is 5
If confidence of meeting quota is 0.01, then finishing house is 5
If confidence of meeting quota is 0.05, then finishing house is 6
If confidence of meeting quota is 0.1, then finishing house is 7
If confidence of meeting quota is 0.5, then finishing house is 12
If confidence of meeting quota is 0.9, then finishing house is 18
If confidence of meeting quota is 0.95, then finishing house is 21
If confidence of meeting quota is 0.99, then finishing house is 25
If confidence of meeting quota is 0.999, then finishing house is 32
If confidence of meeting quota is 1, then finishing house is 1.#INF
If we demand a confidence of meeting sales quota of unity, then we can never be certain of selling 5 bars, so the finishing house is infinite!
Confidence (%) Minimum houses for 5 sales
 0.000 5
 0.001 6
 50.000 10
 75.000 11
 90.000 13
 95.000 15
 99.000 18
 99.900 21
 99.990 25
 99.999 28


 Diagnostic output+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

 Build started: Project: negative_binomial_example1, Configuration: Debug Win32 
Compiling...
negative_binomial_example1.cpp
Linking...
Autorun "i:\boost0605031300\libs\math\test\Math_test\debug\negative_binomial_example1.exe"
Example 1 using the Negative Binomial Distribution. ..\..\..\..\..\..\boostsandbox\math_toolkit\libs\math\example\negative_binomial_example1.cpp Mon Aug 6 15:06:46 2007 140050727
Selling candy bars  an example of using the negative binomial distribution.
An example by Dr. Diane Evans,
+Example 1 using the Negative Binomial Distribution. ..\..\..\..\..\..\boostsandbox\math_toolkit\libs\math\example\negative_binomial_example1.cpp Sun Aug 12 22:28:11 2007
+Selling candy bars  using the negative binomial distribution.
+by Dr. Diane Evans,
Professor of Mathematics at RoseHulman Institute of Technology,
see http://en.wikipedia.org/wiki/Negative_binomial_distribution
Pat has a sales per house success rate of 0.4.
@@ 414,36 +237,38 @@
If the confidence of meeting sales quota is 0.17367, then the finishing house is 8
If the confidence of meeting sales quota is 1, then the finishing house is 1.#INF
If the confidence of meeting sales quota is 0, then the finishing house is 5
If the confidence of meeting sales quota is 0.5, then the finishing house is 12
If the confidence of meeting sales quota is 0.99849, then the finishing house is 31
+If the confidence of meeting sales quota is 0.5, then the finishing house is 11.337
+If the confidence of meeting sales quota is 0.99849, then the finishing house is 30
If confidence of meeting quota is zero (we assume all houses are successful sales), then finishing house is 5
If confidence of meeting quota is 0, then finishing house is 5
If confidence of meeting quota is 0.001, then finishing house is 5
If confidence of meeting quota is 0.01, then finishing house is 5
If confidence of meeting quota is 0.05, then finishing house is 6
If confidence of meeting quota is 0.1, then finishing house is 7
+If confidence of meeting quota is 0.05, then finishing house is 7
+If confidence of meeting quota is 0.1, then finishing house is 8
If confidence of meeting quota is 0.5, then finishing house is 12
If confidence of meeting quota is 0.9, then finishing house is 18
If confidence of meeting quota is 0.95, then finishing house is 21
If confidence of meeting quota is 0.99, then finishing house is 25
If confidence of meeting quota is 0.999, then finishing house is 32
If confidence of meeting quota is 1, then finishing house is 1.#INF
If we demand a confidence of meeting sales quota of unity, then we can never be certain of selling 5 bars, so the finishing house is infinite!
+If we demand a confidence of meeting sales quota of unity,
+then we can never be certain of selling 5 bars,
+so the finishing house is infinite!
Probability (%) House for (last) 5 th sale.
0.01024 5
0.04096 6
0.096256 7
+0.096256 7
0.17367 8
0.26657 9
0.3669 10
+0.3669 10
0.46723 11
0.56182 12 <<<<<<<<<<<< 0.5 should at least 12?
+0.56182 12
0.64696 13
0.72074 14
0.78272 15
0.83343 16
0.874 17
0.90583 18 <<<<<<<<<<<<<<<<< 0.9 should be at least 18
+0.874 17
+0.90583 18
0.93039 19
0.94905 20
0.96304 21
@@ 454,30 +279,12 @@
0.99337 26
0.99539 27
0.99681 28
0.9978 29
+0.9978 29
0.99849 30
0.99897 31
0.9993 32
+0.9993 32
0.99952 33
0.99968 34

Confidence (%) Minimum houses for 5 sales
 0.000 5
 0.001 6
 50.000 10
 75.000 11
 90.000 13
 95.000 15
 99.000 18
 99.900 21
 99.990 25
 99.999 28
 100.000 40
Build Time 0:03
Build log was saved at "file://i:\boost0605031300\libs\math\test\Math_test\negative_binomial_example1\Debug\BuildLog.htm"
negative_binomial_example1  0 error(s), 0 warning(s)
========== Build: 1 succeeded, 0 failed, 0 uptodate, 0 skipped ==========

+0.99968 34
*/
Modified: sandbox/math_toolkit/libs/math/example/negative_binomial_example2.cpp
==============================================================================
 sandbox/math_toolkit/libs/math/example/negative_binomial_example2.cpp (original)
+++ sandbox/math_toolkit/libs/math/example/negative_binomial_example2.cpp 20070814 04:41:13 EDT (Tue, 14 Aug 2007)
@@ 7,8 +7,7 @@
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
// Simple examples demonstrating use of the Negative Binomial Distribution.
// (See other examples for practical applications).
+// Simple example demonstrating use of the Negative Binomial Distribution.
#include <boost/math/distributions/negative_binomial.hpp>
using boost::math::negative_binomial_distribution;
@@ 33,63 +32,67 @@
int main()
{
 cout << "negative_binomial distribution  simple example 2" << endl;

 // Construct distribution: 8 successes (r), 0.25 success fraction = 25% or 1 in 4 successes.
 // negative_binomial_distribution<double> my8dist(8, 0.25);
 negative_binomial my8dist(8, 0.25); // Shorter method using typedef.

 cout.precision(17); // max_digits10 precision
 cout << "mean(my8dist) = " << mean(my8dist) << endl; // 24
 cout << "my8dist.successes() = " << my8dist.successes() << endl;
+ cout << "Negative_binomial distribution  simple example 2" << endl;
+ // Construct a negative binomial distribution with:
+ // 8 successes (r), success fraction (p) 0.25 = 25% or 1 in 4 successes.
+ negative_binomial mynbdist(8, 0.25); // Shorter method using typedef.
+
+ // Display (to check) properties of the distribution just constructed.
+ cout << "mean(mynbdist) = " << mean(mynbdist) << endl; // 24
+ cout << "mynbdist.successes() = " << mynbdist.successes() << endl; // 8
// r th successful trial, after k failures, is r + k th trial.
 cout << "my8dist.success_fraction() = " << my8dist.success_fraction() << endl;
+ cout << "mynbdist.success_fraction() = " << mynbdist.success_fraction() << endl;
// success_fraction = failures/successes or k/r = 0.25 or 25%.
 cout << "my8dist.percent success = " << my8dist.success_fraction() * 100 << "%" << endl;
 cout << "cdf(my8dist, 2.) = " << cdf(my8dist, 2.) << endl; // 0.000415802001953125
 cout << "cdf(my8dist, 8.) = " << cdf(my8dist, 8.) << endl; // 0.027129956288263202
 cout << "cdf(complement(my8dist, 8.)) = " << cdf(complement(my8dist, 8.)) << endl; // 0.9728700437117368
+ cout << "mynbdist.percent success = " << mynbdist.success_fraction() * 100 << "%" << endl;
+ // Show as % too.
+ // Show some cumulative distribution function values for failures k = 2 and 8
+ cout << "cdf(mynbdist, 2.) = " << cdf(mynbdist, 2.) << endl; // 0.000415802001953125
+ cout << "cdf(mynbdist, 8.) = " << cdf(mynbdist, 8.) << endl; // 0.027129956288263202
+ cout << "cdf(complement(mynbdist, 8.)) = " << cdf(complement(mynbdist, 8.)) << endl; // 0.9728700437117368
// Check that cdf plus its complement is unity.
 cout << "cdf + complement = " << cdf(my8dist, 8.) + cdf(complement(my8dist, 8.)) << endl; // 1
+ cout << "cdf + complement = " << cdf(mynbdist, 8.) + cdf(complement(mynbdist, 8.)) << endl; // 1
// Note: No complement for pdf!
// Compare cdf with sum of pdfs.
 double sum = 0.; // of pdfs
 int k = 20;
+ double sum = 0.; // Calculate the sum of all the pdfs,
+ int k = 20; // for 20 failures
for(signed i = 0; i <= k; ++i)
{
 sum += pdf(my8dist, double(i));
+ sum += pdf(mynbdist, double(i));
}
 // Compare with
 double cdf8 = cdf(my8dist, static_cast<double>(k));
 double diff = sum  cdf8;

 cout << setprecision(17) << "Sum pdfs = " << sum << ' ' // 0.40025683281803698
 << ", cdf = " << cdf(my8dist, static_cast<double>(k)) // cdf = 0.40025683281803687
 << ", difference = "
 << diff/ std::numeric_limits<double>::epsilon() // difference = 0.50000000000000000
+ // Compare with the cdf
+ double cdf8 = cdf(mynbdist, static_cast<double>(k));
+ double diff = sum  cdf8; // Expect the diference to be very small.
+ cout << setprecision(17) << "Sum pdfs = " << sum << ' ' // sum = 0.40025683281803698
+ << ", cdf = " << cdf(mynbdist, static_cast<double>(k)) // cdf = 0.40025683281803687
+ << ", difference = " // difference = 0.50000000000000000
+ << setprecision(1) << diff/ (std::numeric_limits<double>::epsilon() * sum)
<< " in epsilon units." << endl;
 // Note: Use boost::math::tools::epsilon rather than std::numeric_limits to cover
 // RealTypes that do not specialize numeric_limits.
+ // Note: Use boost::math::tools::epsilon rather than std::numeric_limits
+ // to cover RealTypes that do not specialize numeric_limits.
+
+//[neg_binomial_example2
 // Print a list of values that can be used to plot
+ // Print a table of values that can be used to plot
// using Excel, or some other superior graphical display tool.
 cout << setprecision(17) << showpoint << endl; // max_digits10 precision, including trailing zeros.
 int maxk = static_cast<int>(2. * my8dist.successes() / my8dist.success_fraction());
 // This maxk shows most of the range of interest, probability about 0.0001 to 0.9999.
 cout << " k pdf cdf""\n" << endl;
+ cout.precision(17); // Use max_digits10 precision, the maximum available for a reference table.
+ cout << showpoint << endl; // include trailing zeros.
+ // This is a maximum possible precision for the type (here double) to suit a reference table.
+ int maxk = static_cast<int>(2. * mynbdist.successes() / mynbdist.success_fraction());
+ // This maxk shows most of the range of interest, probability about 0.0001 to 0.999.
+ cout << "\n"" k pdf cdf""\n" << endl;
for (int k = 0; k < maxk; k++)
{
cout << right << setprecision(17) << showpoint
<< right << setw(3) << k << ", "
 << left << setw(25) << pdf(my8dist, static_cast<double>(k)) << ", "
 << left << setw(25) << cdf(my8dist, static_cast<double>(k))
+ << left << setw(25) << pdf(mynbdist, static_cast<double>(k))
+ << left << setw(25) << cdf(mynbdist, static_cast<double>(k))
<< endl;
}
cout << endl;

+//] [/ neg_binomial_example2]
return 0;
} // int main()
@@ 98,80 +101,82 @@
Output is:
negative_binomial distribution  simple example 2
mean(my8dist) = 24
my8dist.successes() = 8
my8dist.success_fraction() = 0.25
my8dist.percent success = 25%
cdf(my8dist, 2.) = 0.000415802001953125
cdf(my8dist, 8.) = 0.027129956288263202
cdf(complement(my8dist, 8.)) = 0.9728700437117368
+mean(mynbdist) = 24
+mynbdist.successes() = 8
+mynbdist.success_fraction() = 0.25
+mynbdist.percent success = 25%
+cdf(mynbdist, 2.) = 0.000415802001953125
+cdf(mynbdist, 8.) = 0.027129956288263202
+cdf(complement(mynbdist, 8.)) = 0.9728700437117368
cdf + complement = 1
Sum pdfs = 0.40025683281803692 , cdf = 0.40025683281803687, difference = 0.25 in epsilon units.
 k pdf cdf
 0, 1.5258789062500000e005 , 1.5258789062500003e005
 1, 9.1552734375000000e005 , 0.00010681152343750000
 2, 0.00030899047851562522 , 0.00041580200195312500
 3, 0.00077247619628906272 , 0.0011882781982421875
 4, 0.0015932321548461918 , 0.0027815103530883789
 5, 0.0028678178787231476 , 0.0056493282318115234
 6, 0.0046602040529251142 , 0.010309532284736633
 7, 0.0069903060793876605 , 0.017299838364124298
 8, 0.0098301179241389001 , 0.027129956288263202
 9, 0.013106823898851871 , 0.040236780187115073
 10, 0.016711200471036140 , 0.056947980658151209
 11, 0.020509200578089786 , 0.077457181236241013
 12, 0.024354675686481652 , 0.10181185692272265
 13, 0.028101548869017230 , 0.12991340579173993
 14, 0.031614242477644432 , 0.16152764826938440
 15, 0.034775666725408917 , 0.19630331499479325
 16, 0.037492515688331451 , 0.23379583068312471
 17, 0.039697957787645101 , 0.27349378847076977
 18, 0.041352039362130305 , 0.31484582783290005
 19, 0.042440250924291580 , 0.35728607875719176
 20, 0.042970754060845245 , 0.40025683281803687
 21, 0.042970754060845225 , 0.44322758687888220
 22, 0.042482450037426581 , 0.48571003691630876
 23, 0.041558918514873783 , 0.52726895543118257
 24, 0.040260202311284021 , 0.56752915774246648
 25, 0.038649794218832620 , 0.60617895196129912
 26, 0.036791631035234917 , 0.64297058299653398
 27, 0.034747651533277427 , 0.67771823452981139
 28, 0.032575923312447595 , 0.71029415784225891
 29, 0.030329307911589130 , 0.74062346575384819
 30, 0.028054609818219924 , 0.76867807557206813
 31, 0.025792141284492545 , 0.79447021685656061
 32, 0.023575629142856460 , 0.81804584599941710
 33, 0.021432390129869489 , 0.83947823612928651
 34, 0.019383705779220189 , 0.85886194190850684
 35, 0.017445335201298231 , 0.87630727710980494
 36, 0.015628112784496322 , 0.89193538989430121
 37, 0.013938587078064250 , 0.90587397697236549
 38, 0.012379666154859701 , 0.91825364312722524
 39, 0.010951243136991251 , 0.92920488626421649
 40, 0.0096507830144735539 , 0.93885566927869002
 41, 0.0084738582566109364 , 0.94732952753530097
 42, 0.0074146259745345548 , 0.95474415350983555
 43, 0.0064662435824429246 , 0.96121039709227851
 44, 0.0056212231142827853 , 0.96683162020656122
 45, 0.0048717266990450708 , 0.97170334690560634
 46, 0.0042098073105878630 , 0.97591315421619418
 47, 0.0036275999165703964 , 0.97954075413276465
 48, 0.0031174686783026818 , 0.98265822281106729
 49, 0.0026721160099737302 , 0.98533033882104104
 50, 0.0022846591885275322 , 0.98761499800956853
 51, 0.0019486798960970148 , 0.98956367790566557
 52, 0.0016582516423517923 , 0.99122192954801736
 53, 0.0014079495076571762 , 0.99262987905567457
 54, 0.0011928461106539983 , 0.99382272516632852
 55, 0.0010084971662802015 , 0.99483122233260868
 56, 0.00085091948404891532 , 0.99568214181665760
 57, 0.00071656377604119542 , 0.99639870559269883
 58, 0.00060228420831048650 , 0.99700098980100937
 59, 0.00050530624256557675 , 0.99750629604357488
 60, 0.00042319397814867202 , 0.99792949002172360
 61, 0.00035381791615708398 , 0.99828330793788067
 62, 0.00029532382517950324 , 0.99857863176306016
 63, 0.00024610318764958566 , 0.99882473495070978
+//[neg_binomial_example2_1
+ k pdf cdf
+ 0, 1.5258789062500000e005 1.5258789062500003e005
+ 1, 9.1552734375000000e005 0.00010681152343750000
+ 2, 0.00030899047851562522 0.00041580200195312500
+ 3, 0.00077247619628906272 0.0011882781982421875
+ 4, 0.0015932321548461918 0.0027815103530883789
+ 5, 0.0028678178787231476 0.0056493282318115234
+ 6, 0.0046602040529251142 0.010309532284736633
+ 7, 0.0069903060793876605 0.017299838364124298
+ 8, 0.0098301179241389001 0.027129956288263202
+ 9, 0.013106823898851871 0.040236780187115073
+ 10, 0.016711200471036140 0.056947980658151209
+ 11, 0.020509200578089786 0.077457181236241013
+ 12, 0.024354675686481652 0.10181185692272265
+ 13, 0.028101548869017230 0.12991340579173993
+ 14, 0.031614242477644432 0.16152764826938440
+ 15, 0.034775666725408917 0.19630331499479325
+ 16, 0.037492515688331451 0.23379583068312471
+ 17, 0.039697957787645101 0.27349378847076977
+ 18, 0.041352039362130305 0.31484582783290005
+ 19, 0.042440250924291580 0.35728607875719176
+ 20, 0.042970754060845245 0.40025683281803687
+ 21, 0.042970754060845225 0.44322758687888220
+ 22, 0.042482450037426581 0.48571003691630876
+ 23, 0.041558918514873783 0.52726895543118257
+ 24, 0.040260202311284021 0.56752915774246648
+ 25, 0.038649794218832620 0.60617895196129912
+ 26, 0.036791631035234917 0.64297058299653398
+ 27, 0.034747651533277427 0.67771823452981139
+ 28, 0.032575923312447595 0.71029415784225891
+ 29, 0.030329307911589130 0.74062346575384819
+ 30, 0.028054609818219924 0.76867807557206813
+ 31, 0.025792141284492545 0.79447021685656061
+ 32, 0.023575629142856460 0.81804584599941710
+ 33, 0.021432390129869489 0.83947823612928651
+ 34, 0.019383705779220189 0.85886194190850684
+ 35, 0.017445335201298231 0.87630727710980494
+ 36, 0.015628112784496322 0.89193538989430121
+ 37, 0.013938587078064250 0.90587397697236549
+ 38, 0.012379666154859701 0.91825364312722524
+ 39, 0.010951243136991251 0.92920488626421649
+ 40, 0.0096507830144735539 0.93885566927869002
+ 41, 0.0084738582566109364 0.94732952753530097
+ 42, 0.0074146259745345548 0.95474415350983555
+ 43, 0.0064662435824429246 0.96121039709227851
+ 44, 0.0056212231142827853 0.96683162020656122
+ 45, 0.0048717266990450708 0.97170334690560634
+ 46, 0.0042098073105878630 0.97591315421619418
+ 47, 0.0036275999165703964 0.97954075413276465
+ 48, 0.0031174686783026818 0.98265822281106729
+ 49, 0.0026721160099737302 0.98533033882104104
+ 50, 0.0022846591885275322 0.98761499800956853
+ 51, 0.0019486798960970148 0.98956367790566557
+ 52, 0.0016582516423517923 0.99122192954801736
+ 53, 0.0014079495076571762 0.99262987905567457
+ 54, 0.0011928461106539983 0.99382272516632852
+ 55, 0.0010084971662802015 0.99483122233260868
+ 56, 0.00085091948404891532 0.99568214181665760
+ 57, 0.00071656377604119542 0.99639870559269883
+ 58, 0.00060228420831048650 0.99700098980100937
+ 59, 0.00050530624256557675 0.99750629604357488
+ 60, 0.00042319397814867202 0.99792949002172360
+ 61, 0.00035381791615708398 0.99828330793788067
+ 62, 0.00029532382517950324 0.99857863176306016
+ 63, 0.00024610318764958566 0.99882473495070978
+//] [neg_binomial_example2_1 end of Quickbook]
*/
Modified: sandbox/math_toolkit/libs/math/example/negative_binomial_example3.cpp
==============================================================================
 sandbox/math_toolkit/libs/math/example/negative_binomial_example3.cpp (original)
+++ sandbox/math_toolkit/libs/math/example/negative_binomial_example3.cpp 20070814 04:41:13 EDT (Tue, 14 Aug 2007)
@@ 9,38 +9,186 @@
// Example 3 of using constructing distributions, mainly negative_binomial.
+//[neg_binomial_example3
+/*`
+First we need some includes to access the negative binomial distribution
+(and some basic std output of course).
+*/
+
#include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
 using boost::math::negative_binomial_distribution; // default type is double.
using boost::math::negative_binomial; // typedef provides default type is double.
#include <boost/math/distributions/binomial.hpp> // for negative_binomial_distribution
#include <iostream>
using std::cout;
using std::endl;
+ using std::setw;
+#include <limits>
+ using std::numeric_limits;
+
+/*`
+A number of the application examples from K. Krishnamoorthy, Handbook of Statistical Distributions with Applications,
+ISBN 1 58488 635 8, page 100... are implemented using the Math Toolkit library.
int main()
{
 cout << "Example 3 constructing Distributions (Negative_binomial).";
 #if defined(__FILE__) && defined(__TIMESTAMP__)
 cout << " " << __FILE__ << ' ' << __TIMESTAMP__ << ' '<< _MSC_FULL_VER << "\n";
 #endif
 cout << endl;
+A sample example is shown here:
+*/
 // TODO!
+//] [/ neg_binomial_example3]
 return 0;
+int main()
+{
+ cout << "Example 3 Negative_binomial, Krishnamoorthy applications.";
+#if defined(__FILE__) && defined(__TIMESTAMP__)
+ cout << " " << __FILE__ << ' ' << __TIMESTAMP__ << "\n";
+#endif
+ cout << endl;
+ {
+//[neg_binomial_example3_1
+ // K. Krishnamoorthy, Handbook of Statistical Distributions with Applications,
+ // ISBN 1 58488 635 8, page 100, example 7.3.1
+ // r successes = 20, k failures = 18, and success probability (fraction) = 0.6
+ negative_binomial nbk(20, 0.6); // successes, success fraction.
+ cout.precision(6); // default precision.
+ cout << "probability of < 18 = cdf(18) = " << cdf(nbk, 18) << endl; // = 0.862419
+ cout << "probability of >= 18 = 1  cdf(17) = " << 1  cdf(nbk, 17) << endl; // = 0.181983
+ // But this may suffer from inaccuracy, so mcuh better to use the complement.
+ cout << "probability of >= 18 = cdf(complement(nbk, 17)) = "
+ << cdf(complement(nbk, 17)) << endl; // = 0.181983
+ cout << "probability of == 18 = pdf(18) = " << pdf(nbk, 18) << endl; // 0.044402
+ cout << "probability of exactly 18 = pdf(18) = " << pdf(nbk, 18) << endl; // 0.044402
+//] [/negative_binomial_example3_1 Quickbook end]
+ }
+ { // Example 7.3.2 to find the success probability when r = 4, k = 5, P(X <= k) = 0.56 and
+ // calculates success probability = 0.417137.
+ negative_binomial nb(4, 0.417137);
+ cout << cdf(nb, 5) << endl; // expect P(X <= k) = 0.56 got 0.560001
+ // Now check the inverse by calculating the k failures.
+ cout << quantile(nb, 0.56) << endl; // expect to get k = 5. and got 5.000000
+ // P(X <= k) = 0.56001
+ // P(X >= k) = 0.55406
+ // P(X == k) = 0.114061
+
+ // mean 5.58918, sd 3.66045, mode = 4, cv = 0.654918, skew 1.03665, Mean dev 2.86877, kurtosis 4.57463
+ cout << "Mean = " << mean(nb) // 5.589176
+ << ", sd = " << standard_deviation(nb)
+ << ", Coefficient of variation (sd/mean) = " << standard_deviation(nb) / mean(nb)
+ << ", mode = " << mode(nb) << ", skew " << skewness(nb)
+ << ", mean deviation = " << "??" // perhaps todo?
+ << ", kurt = " << kurtosis(nb) << endl;
+ }
+ { // 7.3.3 Coin flipping. What are chances that 10th head will occur at the 12th flip.
+ negative_binomial nb(10, 0.5);
+
+ cout << "Probability that 10th head will occur at the 12th flip is " << pdf(nb, 1210) << endl; // 0.013428
+ cout << "P(X == k) " << pdf(nb, 1210) << endl; // 0.013428
+
+ cout << "Probability that 10th head will occur before the 12th flip is " << cdf(nb, 1) << endl; //
+ cout << "P(X < k) is " << cdf(nb, 1) << endl; // 1 = 12  11, // 0.005859
+
+ cout << "Probability that 10th head will occur on or before the 12th flip is " << cdf(nb, 2) << endl; //
+ cout << "P(X <= k) is " << cdf(nb, 2) << endl; // P(X <= k) = 0.01928711
+
+ cout << "Probability that 10th head will occur on or after the 12th flip is " << 1  cdf(nb, 1) << endl; //
+ cout << "P(X >= k) is " << cdf(complement(nb, 1211)) << endl; // P(X >= k) = 0.994140625
+
+ cout << "Probability that 10th head will occur after the 12th flip is " << 1  cdf(nb, 2) << endl; //
+ cout << "P(X > k) is " << cdf(complement(nb, 2)) << endl; // P(X > k) 0.980713
+ /*
+ Probability that 10th head will occur at the 12th flip is 0.013428
+ P(X == k) 0.013428
+ Probability that 10th head will occur before the 12th flip is 0.005859
+ P(X < k) is 0.005859
+ Probability that 10th head will occur on or before the 12th flip is 0.019287
+ P(X <= k) is 0.019287
+ Probability that 10th head will occur on or after the 12th flip is 0.994141
+ P(X >= k) is 0.994141
+ Probability that 10th head will occur after the 12th flip is 0.980713
+ P(X > k) is 0.980713
+ */
+ }
+ { // Example 7.3.4, Sample up to 30 items from each lot.
+ // Chances of rejecting the lot if it indeed contains 15% defectives?
+ // stop at the 3rd defective, or go on to sample all 30.
+ negative_binomial nb(3, 0.15); // 3rd defective.
+ // probability of observing 27 or less OK items to get the 3rd defective.
+ // P(X <= k) = 0.8485994
+ cout << "Chance of rejecting the lot is " << cdf(nb, 27) << " if the lot actually contains 15% defectives." << endl; // 0.848599
+ }
+ { // K. Krishnamoorthy, ISBN 1 58488 635 8, page 101, section 7.46 confidence Intervals for the Proportion
+ // 7.6.1 Shipment inspected onbyone randomly and found 6th defective at 30th inspection.
+ // 2sided 95% confidence 0.0771355 to 0.357748
+ double alpha = 0.05; // Note divide by 2 if twosided test.
+ double successes = 6;
+ double failures = 24;
+ double trials = successes + failures;
+ double l = negative_binomial::find_lower_bound_on_p(trials, successes, alpha/2); // 0.077136
+ double u = negative_binomial::find_upper_bound_on_p(trials, successes, alpha/2); // 0.357748
+ // These find_ use the ClopperPearson approach.
+ cout << (1  alpha) * 100 << "% confidence interval low " << l << ", up " << u << endl;
+ cout << "So we conclude that the true percentage of defective items is between "
+ << static_cast<int>(l * 100) << " and " << setw(2) << u * 100 << "%." << endl;
+
+ cout << "The rounding style for a double type is: "
+ << numeric_limits<double>::round_style << endl;
+ }
+ { // K. Krishnamoorthy, ISBN 1 58488 635 8, page 102, section 7.4
+ // Suppose required k + r trials to get the rth success.
+ double r = 5; double k = 25; // Example 7.5.1 values
+ double rm1 = r 1; // r1 is the penultimate success.
+ double p = rm1/ (rm1 + k);
+ // A point estimate of the actual proportion of defective items.
+ // 0.137931
+ // 'True' estimate, if this was all the items ever available is successes/failures = r/k = 5/25 = 0.2
+ // so the point estimate 'how we are doing from the info so far' is rather less at 0.14.
+ cout << "Uniformly minimum variance unbiased estimator of success probability is " << p << endl;
+ double v = 0.5 * p * p * (1  p) * (k + k + 2  p) / (k * (k  p + 2));
+ cout << "Variance of estimator is " << v << endl; // 0.000633
+ cout << "Standard deviation of estimator is " << sqrt(v) << endl; // 0.025  seems small?
+
+ // Section 7.5 testing the true success probability.
+ // 7.5.1 A lot is inspected randomly. 5th defective found at 30th inspection.
+ // so 25 are OK and 5 are duds.
+ // Using StatCalc r = 5, k = 25, so r + r = 30, value for p0 0.3
+ negative_binomial nb(5, 0.2); // 1/3rd are defective?
+ cout << cdf(nb, 25) << endl; // 0.812924
+ cout << quantile(nb, 0.3) << endl; // 7.344587
+ // TODO?
+ }
+ return 0;
} // int main()
/*

Output is:
+/*
+Example 3 Negative_binomial . ..\..\..\..\..\..\boostsandbox\math_toolkit\libs\math\example\negative_binomial_example3.cpp Mon Aug 13 14:32:28 2007 140050727
+probability of < 18 = cdf(18) = 0.862419
+probability of >= 18 = 1  cdf(17) = 0.181983
+probability of >= 18 = cdf(complement(nbk, 17)) = 0.181983
+probability of == 18 = pdf(18) = 0.0444024
+probability of exactly 18 = pdf(18) = 0.0444024
+0.560001
+5
+Mean = 5.58918, sd = 3.66045, Coefficient of variation (sd/mean) = 0.654918, mode = 4, skew 1.03665, mean deviation = ??, kurt = 4.57463
+Probability that 10th head will occur at the 12th flip is 0.0134277
+P(X == k) 0.0134277
+Probability that 10th head will occur before the 12th flip is 0.00585938
+P(X < k) is 0.00585938
+Probability that 10th head will occur on or before the 12th flip is 0.0192871
+P(X <= k) is 0.0192871
+Probability that 10th head will occur on or after the 12th flip is 0.994141
+P(X >= k) is 0.994141
+Probability that 10th head will occur after the 12th flip is 0.980713
+P(X > k) is 0.980713
+Chance of rejecting the lot is 0.848599 if the lot actually contains 15% defectives.
+95% confidence interval low 0.0771355, up 0.357748
+So we conclude that the true percentage of defective items is between 7 and 35.7748%.
+The rounding style for a double type is: 1
+Uniformly minimum variance unbiased estimator of success probability is 0.137931
+Variance of estimator is 0.000633295
+Standard deviation of estimator is 0.0251654
+0.744767
+13
*/



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