mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 16:32:10 +00:00
209 lines
10 KiB
C++
209 lines
10 KiB
C++
// negative_binomial_example3.cpp
|
|
|
|
// Copyright Paul A. Bristow 2007.
|
|
|
|
// Use, modification and distribution are subject to the
|
|
// Boost Software License, Version 1.0.
|
|
// (See accompanying file LICENSE_1_0.txt
|
|
// or copy at http://www.boost.org/LICENSE_1_0.txt)
|
|
|
|
// 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; // typedef provides default type is double.
|
|
|
|
#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.
|
|
|
|
A sample example is shown here:
|
|
*/
|
|
|
|
//] [/ neg_binomial_example3]
|
|
|
|
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, P(X <= 18) == cdf(18) = " << cdf(nbk, 18) << endl; // = 0.862419
|
|
cout << "Probability of > 18, P(X < 18) == 1 - cdf(18) = " << 1 - cdf(nbk, 18) << endl; // = 0.137581
|
|
// But this may suffer from inaccuracy, so much better to use the complement.
|
|
cout << "Probability of > 18, P(X < 18) == 1 - cdf(18) = " << cdf(complement(nbk, 18)) << endl; // 0.137581
|
|
// And, of course, the sum of probability of X <= 18 and X > 18 really should be unity!
|
|
BOOST_ASSERT(cdf(nbk, 18) + cdf(complement(nbk, 18))); // = 0.862419 + 0.137581 == 1
|
|
cout << "Probability of < 18 == <= 17 == cdf(17) = " << cdf(nbk, 17) << endl; //
|
|
cout << "Probability of >= 18 == > 17 == 1 - cdf(17) = " << 1 - cdf(nbk, 17) << endl; // 0.181983
|
|
// But this may suffer from inaccuracy, so much better to use the complement.
|
|
cout << "Probability of >= 18 == > 17 == cdf(complement(nbk, 17)) = "
|
|
<< cdf(complement(nbk, 17)) << endl; // = 0.181983
|
|
cout << "Probability of exactly == 18, P(X =18) = pdf(18) = " << pdf(nbk, 18) << endl; // 0.044402
|
|
cout << endl;
|
|
//] [/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 << "P(X <= k) = " << cdf(nb, 5) << endl; // expect P(X <= k) = 0.56 got 0.560001
|
|
// Now check the inverse by calculating the k failures.
|
|
cout << " k = " << 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
|
|
// Compute moments:
|
|
// 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 & 7.5
|
|
// 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 phat 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.138.
|
|
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 smallish?
|
|
// Would expect point estimate 0.14 + a couple of sd = 0.25 * 2 = 0.05 = 1.9 ~= 0.2?
|
|
// So perhaps near enough?
|
|
|
|
// 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 << "P(X <= 25) = " << cdf(nb, 25) << endl; // 0.744767
|
|
cout << "quantile(nb, 0.3) = " << quantile(nb, 0.3) << endl; // 13
|
|
// Don't understand what the p-value calculation does.
|
|
// TODO????
|
|
}
|
|
return 0;
|
|
} // int main()
|
|
|
|
|
|
/*
|
|
|
|
Output is:
|
|
|
|
Example 3 Negative_binomial, Krishnamoorthy applications. ..\..\..\..\..\..\boost-sandbox\math_toolkit\libs\math\example\negative_binomial_example3.cpp Wed Aug 15 11:08:14 2007
|
|
Probability of <= 18, P(X <= 18) == cdf(18) = 0.862419
|
|
Probability of > 18, P(X < 18) == 1 - cdf(18) = 0.137581
|
|
Probability of > 18, P(X < 18) == 1 - cdf(18) = 0.137581
|
|
Probability of < 18 == <= 17 == cdf(17) = 0.818017
|
|
Probability of >= 18 == > 17 == 1 - cdf(17) = 0.181983
|
|
Probability of >= 18 == > 17 == cdf(complement(nbk, 17)) = 0.181983
|
|
Probability of exactly == 18, P(X =18) = pdf(18) = 0.0444024
|
|
P(X <= k) = 0.560001
|
|
k = 5
|
|
Mean = 5.58918, sd = 3.66045, Coefficient of variation (sd/mean) = 0.654918, mode = 4, skew 1.03665, 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
|
|
P(X <= 25) = 0.744767
|
|
13
|
|
|
|
|
|
*/
|
|
|
|
|
|
|