|
|
Boost-Commit : |
From: pbristow_at_[hidden]
Date: 2007-08-14 04:41:15
Author: pbristow
Date: 2007-08-14 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 2007-08-14 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 (128-bit floating-point), 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 (128-bit floating-point), 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 2007-08-14 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 Rose-Hulman 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, 30-5+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 Rose-Hulman 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 (non-existent) last-plus-1 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 (non-existent) last-plus-1 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 = 1-alpha, as percent = 100 * (1-alpha[i])
+ // Confidence as fraction = 1-alpha, as percent = 100 * (1-alpha[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 * (1-alpha[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. ..\..\..\..\..\..\Boost-sandbox\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 Rose-Hulman 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:\boost-06-05-03-1300\libs\math\test\Math_test\debug\negative_binomial_example1.exe"
-Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\negative_binomial_example1.exe"
-Example 1 using the Negative Binomial Distribution. ..\..\..\..\..\..\boost-sandbox\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 Rose-Hulman 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:\boost-06-05-03-1300\libs\math\test\Math_test\debug\negative_binomial_example1.exe"
-Example 1 using the Negative Binomial Distribution. ..\..\..\..\..\..\boost-sandbox\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. ..\..\..\..\..\..\boost-sandbox\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 Rose-Hulman 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:\boost-06-05-03-1300\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 up-to-date, 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 2007-08-14 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.5258789062500000e-005 , 1.5258789062500003e-005
- 1, 9.1552734375000000e-005 , 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.5258789062500000e-005 1.5258789062500003e-005
+ 1, 9.1552734375000000e-005 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 2007-08-14 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, 12-10) << endl; // 0.013428
+ cout << "P(X == k) " << pdf(nb, 12-10) << 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, 12-11)) << 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 on-by-one randomly and found 6th defective at 30th inspection.
+ // 2-sided 95% confidence 0.0771355 to 0.357748
+ double alpha = 0.05; // Note divide by 2 if two-sided 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 Clopper-Pearson 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; // r-1 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 . ..\..\..\..\..\..\boost-sandbox\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
*/
-
-
-
Boost-Commit 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