From ccfb864daf3ef682afbe238247e3db637550ac67 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Thu, 30 Nov 2006 10:47:26 +0000 Subject: [PATCH] Added Jeffreys prior method to the binomial confidence limits. Fixed a few typos in the neg binomial. [SVN r3482] --- doc/distributions/binomial.qbk | 74 ++++++++++++++- doc/distributions/binomial_example.qbk | 91 ++++++++++-------- .../negative_binomial_example.qbk | 2 +- example/binomial_confidence_limits.cpp | 92 ++++++++++--------- example/neg_binomial_sample_sizes.cpp | 2 +- include/boost/math/distributions/binomial.hpp | 25 +++-- test/test_binomial.cpp | 32 +++++++ 7 files changed, 223 insertions(+), 95 deletions(-) diff --git a/doc/distributions/binomial.qbk b/doc/distributions/binomial.qbk index 05e166233..9f3fa36c3 100644 --- a/doc/distributions/binomial.qbk +++ b/doc/distributions/binomial.qbk @@ -15,6 +15,11 @@ class binomial_distribution { public: + typedef RealType value_type; + + static const ``['unspecified-type]`` cloppper_pearson_exact_interval; + static const ``['unspecified-type]`` jeffreys_prior_interval; + // construct: binomial_distribution(RealType n, RealType p); @@ -26,11 +31,13 @@ static RealType estimate_lower_bound_on_p( RealType trials, RealType successes, - RealType probability); + RealType probability, + ``['unspecified-type]`` method = clopper_pearson_exact_interval); static RealType estimate_upper_bound_on_p( RealType trials, RealType successes, - RealType probability); + RealType probability, + ``['unspecified-type]`` method = clopper_pearson_exact_interval); // estimate min/max number of trials: static RealType estimate_number_of_trials( @@ -91,7 +98,8 @@ Returns the parameter /n/ from which this distribution was constructed. static RealType estimate_lower_bound_on_p( RealType trials, RealType successes, - RealType alpha); + RealType alpha, + ``['unspecified-type]`` method = clopper_pearson_exact_interval); Returns a lower bound on the success fraction: @@ -100,6 +108,8 @@ Returns a lower bound on the success fraction: [[successes][The number of successes that occurred.]] [[alpha][The largest acceptable probability that the true value of the success fraction is [*less than] the value returned.]] +[[method][An optional parameter that specifies the method to be used + to compute the interval (See below).]] ] For example, if you observe /k/ successes from /n/ trials the @@ -111,11 +121,61 @@ want to be 95% sure that the true value is [*greater than] some value, n, k, 0.05); [link math_toolkit.dist.stat_tut.weg.binom_eg.binom_conf See worked example.] - + +There are currently two possible values available for the /method/ +optional parameter: /clopper_pearson_exact_interval/ +or /jeffreys_prior_interval/. These constants are both members of +class template `binomial_distribution`, so usage is for example: + + p = binomial_distribution::estimate_lower_bound_on_p( + n, k, 0.05, binomial_distribution::jeffreys_prior_interval); + +The default method if this parameter is not specified is the Clopper Pearson +"exact" interval. This produces an interval that guarantees at least +`100(1-alpha)%` coverage, but which is known to be overly conservative, +sometimes producing intervals with much greater than the requested coverage. + +The alternative calculation method produces a non-informative +Jeffreys Prior interval. It produces `100(1-alpha)%` coverage only +['in the average case], though is typically very close to the requested +coverage level. It is one of the main methods of calculation recomended +in the review by Brown, Cai and DasGupta. + +Please note that the "textbook" calculation method using +a normal approximation (the Wald interval) is deliberately +not provided: it is known to produce consistently poor results, +even when the sample size is surprisingly large. +Refer to Brown, Cai and DasGupta for a full explanation. Many other methods +of calculation are available, and may be more appropriate for specific +situations. Unfortunately there appears to be no consensus amoungst +stataticians as to which is "best": refer to the discussion at the end of +Brown, Cai and DasGupta for examples. + +The two methods provided here were chosen principally because they +can be used for both one and two sided intervals. +See also: + +Lawrence D. Brown, T. Tony Cai and Anirban DasGupta (2001), +Interval Estimation for a Binomial Proportion, +Statistical Science, Vol. 16, No. 2, 101-133. + +T. Tony Cai (2005), +One-sided confidence intervals in discrete distributions, +Journal of Statistical Planning and Inference 131, 63-88. + +Agresti, A. and Coull, B. A. (1998). Approximate is better than +"exact" for interval estimation of binomial proportions. Amer. +Statist. 52 119-126. + +Clopper, C. J. and Pearson, E. S. (1934). The use of confidence +or fiducial limits illustrated in the case of the binomial. +Biometrika 26 404-413. + static RealType estimate_upper_bound_on_p( RealType trials, RealType successes, - RealType alpha); + RealType alpha, + ``['unspecified-type]`` method = clopper_pearson_exact_interval); Returns an upper bound on the success fraction: @@ -124,6 +184,10 @@ Returns an upper bound on the success fraction: [[successes][The number of successes that occurred.]] [[alpha][The largest acceptable probability that the true value of the success fraction is [*greater than] the value returned.]] +[[method][An optional parameter that specifies the method to be used + to compute the interval. Refer to the documentation for + `estimate_upper_bound_on_p` above for the meaning of the + method options.]] ] For example, if you observe /k/ successes from /n/ trials the diff --git a/doc/distributions/binomial_example.qbk b/doc/distributions/binomial_example.qbk index eda917cde..63008f307 100644 --- a/doc/distributions/binomial_example.qbk +++ b/doc/distributions/binomial_example.qbk @@ -55,11 +55,11 @@ interval: Some pretty printing of the table header follows: - cout << "\n\n" - "___________________________________________\n" - "Confidence Lower Upper\n" - " Value (%) Limit Limit\n" - "___________________________________________\n"; + cout << "\n\n" + "_______________________________________________________________________\n" + "Confidence Lower CP Upper CP Lower JP Upper JP\n" + " Value (%) Limit Limit Limit Limit\n" + "_______________________________________________________________________\n"; And now for the important part - the intervals themselves - for each @@ -89,16 +89,32 @@ single-sided interval, for example: ['"Calculate a lower bound so that we are P% sure that the true occurrence frequency is greater than some value"] then we would *not* have divided by two. +Finally note that `binomial_distribution` provides a choice of two +methods for the calculation, we print out the results from both +methods in this example: + for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i) { // Confidence value: cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]); - // Calculate bounds of interval: - double l = binomial_distribution<>::estimate_lower_bound_on_p(trials, successes, alpha[i]/2); - double u = binomial_distribution<>::estimate_upper_bound_on_p(trials, successes, alpha[i]/2); - // Print Limits: + // Calculate Clopper Pearson bounds: + double l = binomial_distribution<>::estimate_lower_bound_on_p( + trials, successes, alpha[i]/2); + double u = binomial_distribution<>::estimate_upper_bound_on_p( + trials, successes, alpha[i]/2); + // Print Clopper Pearson Limits: cout << fixed << setprecision(5) << setw(15) << right << l; - cout << fixed << setprecision(5) << setw(15) << right << u << endl; + cout << fixed << setprecision(5) << setw(15) << right << u; + // Calculate Jeffreys Prior Bounds: + l = binomial_distribution<>::estimate_lower_bound_on_p( + trials, successes, alpha[i]/2, + binomial_distribution<>::jeffreys_prior_interval); + u = binomial_distribution<>::estimate_upper_bound_on_p( + trials, successes, alpha[i]/2, + binomial_distribution<>::jeffreys_prior_interval); + // Print Jeffreys Prior Limits: + cout << fixed << setprecision(5) << setw(15) << right << l; + cout << fixed << setprecision(5) << setw(15) << right << u << std::endl; } cout << endl; } @@ -115,25 +131,27 @@ Number of successes = 4 Sample frequency of occurrence = 0.2 -___________________________________________ -Confidence Lower Upper - Value (%) Limit Limit -___________________________________________ - 50.000 0.12840 0.29588 - 75.000 0.09775 0.34633 - 90.000 0.07135 0.40103 - 95.000 0.05733 0.43661 - 99.000 0.03576 0.50661 - 99.900 0.01905 0.58632 - 99.990 0.01042 0.64997 - 99.999 0.00577 0.70216 +_______________________________________________________________________ +Confidence Lower CP Upper CP Lower JP Upper JP + Value (%) Limit Limit Limit Limit +_______________________________________________________________________ + 50.000 0.12840 0.29588 0.14974 0.26916 + 75.000 0.09775 0.34633 0.11653 0.31861 + 90.000 0.07135 0.40103 0.08734 0.37274 + 95.000 0.05733 0.43661 0.07152 0.40823 + 99.000 0.03576 0.50661 0.04655 0.47859 + 99.900 0.01905 0.58632 0.02634 0.55960 + 99.990 0.01042 0.64997 0.01530 0.62495 + 99.999 0.00577 0.70216 0.00901 0.67897 '''] As you can see, even at the 95% confidence level the bounds are really quite wide (this example is chosen to be easily compared to the one in the __handbook [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm -here]). +here]). Note also that the Clopper-Pearson calculation method (CP above) +produces quite noticably more pessimistic estimates than the Jeffreys Prior +method (JP above). Compare that with the program output for @@ -148,22 +166,23 @@ Number of successes = 400 Sample frequency of occurrence = 0.2000000 -___________________________________________ -Confidence Lower Upper - Value (%) Limit Limit -___________________________________________ - 50.000 0.19382 0.20638 - 75.000 0.18965 0.21072 - 90.000 0.18537 0.21528 - 95.000 0.18267 0.21821 - 99.000 0.17745 0.22400 - 99.900 0.17150 0.23079 - 99.990 0.16658 0.23657 - 99.999 0.16233 0.24169 +_______________________________________________________________________ +Confidence Lower CP Upper CP Lower JP Upper JP + Value (%) Limit Limit Limit Limit +_______________________________________________________________________ + 50.000 0.19382 0.20638 0.19406 0.20613 + 75.000 0.18965 0.21072 0.18990 0.21047 + 90.000 0.18537 0.21528 0.18561 0.21503 + 95.000 0.18267 0.21821 0.18291 0.21796 + 99.000 0.17745 0.22400 0.17769 0.22374 + 99.900 0.17150 0.23079 0.17173 0.23053 + 99.990 0.16658 0.23657 0.16681 0.23631 + 99.999 0.16233 0.24169 0.16256 0.24143 '''] Now even when the confidence level is very high, the limits are really -quite close to the experimentally calculated value of 0.2. +quite close to the experimentally calculated value of 0.2. Furthermore +the difference between the two calculation methods is now really quite small. [endsect] diff --git a/doc/distributions/negative_binomial_example.qbk b/doc/distributions/negative_binomial_example.qbk index 2e7a58ad6..5f57b65cf 100644 --- a/doc/distributions/negative_binomial_example.qbk +++ b/doc/distributions/negative_binomial_example.qbk @@ -200,7 +200,7 @@ Note that since we're calculating the /minimum/ number of trials required, we'll err on the safe side and take the ceiling of the result. Had we been calculating the /maximum/ number of trials permitted to observe less than a certain -number of /failures/ then we would have taken the ceiling instead. We +number of /failures/ then we would have taken the floor instead. We would also have wrapped the arguments to `estimate_number_of_trials` in a call to complement like this: diff --git a/example/binomial_confidence_limits.cpp b/example/binomial_confidence_limits.cpp index 29a645cdc..5ddc6d72a 100644 --- a/example/binomial_confidence_limits.cpp +++ b/example/binomial_confidence_limits.cpp @@ -44,10 +44,10 @@ void confidence_limits_on_frequency(unsigned trials, unsigned successes) // Print table header: // cout << "\n\n" - "___________________________________________\n" - "Confidence Lower Upper\n" - " Value (%) Limit Limit\n" - "___________________________________________\n"; + "_______________________________________________________________________\n" + "Confidence Lower CP Upper CP Lower JP Upper JP\n" + " Value (%) Limit Limit Limit Limit\n" + "_______________________________________________________________________\n"; // // Now print out the data for the table rows. // @@ -55,12 +55,18 @@ void confidence_limits_on_frequency(unsigned trials, unsigned successes) { // Confidence value: cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]); - // calculate bounds: + // Calculate Clopper Pearson bounds: double l = binomial_distribution<>::estimate_lower_bound_on_p(trials, successes, alpha[i]/2); double u = binomial_distribution<>::estimate_upper_bound_on_p(trials, successes, alpha[i]/2); - // Print Limits: + // Print Clopper Pearson Limits: cout << fixed << setprecision(5) << setw(15) << right << l; - cout << fixed << setprecision(5) << setw(15) << right << u << endl; + cout << fixed << setprecision(5) << setw(15) << right << u; + // Calculate Jeffreys Prior Bounds: + l = binomial_distribution<>::estimate_lower_bound_on_p(trials, successes, alpha[i]/2, binomial_distribution<>::jeffreys_prior_interval); + u = binomial_distribution<>::estimate_upper_bound_on_p(trials, successes, alpha[i]/2, binomial_distribution<>::jeffreys_prior_interval); + // Print Jeffreys Prior Limits: + cout << fixed << setprecision(5) << setw(15) << right << l; + cout << fixed << setprecision(5) << setw(15) << right << u << std::endl; } cout << endl; } // void confidence_limits_on_frequency() @@ -90,18 +96,18 @@ Number of successes = 4 Sample frequency of occurrence = 0.2 -___________________________________________ -Confidence Lower Upper - Value (%) Limit Limit -___________________________________________ - 50.000 0.12840 0.29588 - 75.000 0.09775 0.34633 - 90.000 0.07135 0.40103 - 95.000 0.05733 0.43661 - 99.000 0.03576 0.50661 - 99.900 0.01905 0.58632 - 99.990 0.01042 0.64997 - 99.999 0.00577 0.70216 +_______________________________________________________________________ +Confidence Lower CP Upper CP Lower JP Upper JP + Value (%) Limit Limit Limit Limit +_______________________________________________________________________ + 50.000 0.12840 0.29588 0.14974 0.26916 + 75.000 0.09775 0.34633 0.11653 0.31861 + 90.000 0.07135 0.40103 0.08734 0.37274 + 95.000 0.05733 0.43661 0.07152 0.40823 + 99.000 0.03576 0.50661 0.04655 0.47859 + 99.900 0.01905 0.58632 0.02634 0.55960 + 99.990 0.01042 0.64997 0.01530 0.62495 + 99.999 0.00577 0.70216 0.00901 0.67897 ___________________________________________ 2-Sided Confidence Limits For Success Ratio @@ -112,18 +118,18 @@ Number of successes = 40 Sample frequency of occurrence = 0.2000000 -___________________________________________ -Confidence Lower Upper - Value (%) Limit Limit -___________________________________________ - 50.000 0.17949 0.22259 - 75.000 0.16701 0.23693 - 90.000 0.15455 0.25225 - 95.000 0.14689 0.26223 - 99.000 0.13257 0.28218 - 99.900 0.11703 0.30601 - 99.990 0.10489 0.32652 - 99.999 0.09492 0.34485 +_______________________________________________________________________ +Confidence Lower CP Upper CP Lower JP Upper JP + Value (%) Limit Limit Limit Limit +_______________________________________________________________________ + 50.000 0.17949 0.22259 0.18190 0.22001 + 75.000 0.16701 0.23693 0.16934 0.23429 + 90.000 0.15455 0.25225 0.15681 0.24956 + 95.000 0.14689 0.26223 0.14910 0.25951 + 99.000 0.13257 0.28218 0.13468 0.27940 + 99.900 0.11703 0.30601 0.11902 0.30318 + 99.990 0.10489 0.32652 0.10677 0.32366 + 99.999 0.09492 0.34485 0.09670 0.34197 ___________________________________________ 2-Sided Confidence Limits For Success Ratio @@ -134,18 +140,18 @@ Number of successes = 400 Sample frequency of occurrence = 0.2000000 -___________________________________________ -Confidence Lower Upper - Value (%) Limit Limit -___________________________________________ - 50.000 0.19382 0.20638 - 75.000 0.18965 0.21072 - 90.000 0.18537 0.21528 - 95.000 0.18267 0.21821 - 99.000 0.17745 0.22400 - 99.900 0.17150 0.23079 - 99.990 0.16658 0.23657 - 99.999 0.16233 0.24169 +_______________________________________________________________________ +Confidence Lower CP Upper CP Lower JP Upper JP + Value (%) Limit Limit Limit Limit +_______________________________________________________________________ + 50.000 0.19382 0.20638 0.19406 0.20613 + 75.000 0.18965 0.21072 0.18990 0.21047 + 90.000 0.18537 0.21528 0.18561 0.21503 + 95.000 0.18267 0.21821 0.18291 0.21796 + 99.000 0.17745 0.22400 0.17769 0.22374 + 99.900 0.17150 0.23079 0.17173 0.23053 + 99.990 0.16658 0.23657 0.16681 0.23631 + 99.999 0.16233 0.24169 0.16256 0.24143 */ diff --git a/example/neg_binomial_sample_sizes.cpp b/example/neg_binomial_sample_sizes.cpp index 80d272a6d..ff5ccc056 100644 --- a/example/neg_binomial_sample_sizes.cpp +++ b/example/neg_binomial_sample_sizes.cpp @@ -34,7 +34,7 @@ void estimate_number_of_trials(double failures, double p) // p = success ratio. // successes = required number of successes. // - // Calculate how many trials we can have to ensure the + // Calculate how many trials we need to ensure the // required number of failures DOES exceed "failures". // Define a table of significance levels: diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp index 2ab8f841d..01e0452ba 100644 --- a/include/boost/math/distributions/binomial.hpp +++ b/include/boost/math/distributions/binomial.hpp @@ -188,18 +188,23 @@ namespace boost { // Total number of trials. return m_n; } + + enum interval_type{ + clopper_pearson_exact_interval, + jeffreys_prior_interval + }; + // // Estimation of the success fraction parameter. // The best estimate is actually simply successes/trials, // these functions are used // to obtain confidence intervals for the success fraction. - // Note these functions return the Clopper-Pearson interval - // which is known to be overly conservative. // static RealType estimate_lower_bound_on_p( RealType trials, RealType successes, - RealType probability) + RealType probability, + interval_type t = clopper_pearson_exact_interval) { // Error checks: RealType result; @@ -213,15 +218,17 @@ namespace boost if(successes == 0) return 0; - // NOTE!!! This formual uses "successes" not "successes+1" as usual - // to get the lower bound, + // NOTE!!! The Clopper Pearson formula uses "successes" not + // "successes+1" as usual to get the lower bound, // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm - return ibeta_inv(successes, trials - successes + 1, probability); + return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability) + : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability); } static RealType estimate_upper_bound_on_p( RealType trials, RealType successes, - RealType probability) + RealType probability, + interval_type t = clopper_pearson_exact_interval) { // Error checks: RealType result; @@ -235,9 +242,9 @@ namespace boost if(trials == successes) return 1; - return ibetac_inv(successes + 1, trials - successes, probability); + return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability) + : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability); } - // Estimate number of trials parameter: // // "How many trials do I need to be P% sure of seeing k events?" diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp index d002ccd68..45f118f0d 100644 --- a/test/test_binomial.cpp +++ b/test/test_binomial.cpp @@ -100,6 +100,7 @@ void test_spot( if(Q < P) { + // Default method (Clopper Pearson) BOOST_CHECK( binomial_distribution::estimate_lower_bound_on_p( N, k, Q) @@ -114,9 +115,25 @@ void test_spot( binomial_distribution::estimate_upper_bound_on_p( N, k, Q)) ); + // Bayes Method (Jeffreys Prior) + BOOST_CHECK( + binomial_distribution::estimate_lower_bound_on_p( + N, k, Q, binomial_distribution::jeffreys_prior_interval) + <= + binomial_distribution::estimate_upper_bound_on_p( + N, k, Q, binomial_distribution::jeffreys_prior_interval) + ); + BOOST_CHECK(( + binomial_distribution::estimate_lower_bound_on_p( + N, k, Q, binomial_distribution::jeffreys_prior_interval) + <= k/N) && (k/N <= + binomial_distribution::estimate_upper_bound_on_p( + N, k, Q, binomial_distribution::jeffreys_prior_interval)) + ); } else { + // Default method (Clopper Pearson) BOOST_CHECK( binomial_distribution::estimate_lower_bound_on_p( N, k, P) @@ -131,6 +148,21 @@ void test_spot( binomial_distribution::estimate_upper_bound_on_p( N, k, P)) ); + // Bayes Method (Jeffreys Prior) + BOOST_CHECK( + binomial_distribution::estimate_lower_bound_on_p( + N, k, P, binomial_distribution::jeffreys_prior_interval) + <= + binomial_distribution::estimate_upper_bound_on_p( + N, k, P, binomial_distribution::jeffreys_prior_interval) + ); + BOOST_CHECK( + (binomial_distribution::estimate_lower_bound_on_p( + N, k, P, binomial_distribution::jeffreys_prior_interval) + <= k / N) && (k/N <= + binomial_distribution::estimate_upper_bound_on_p( + N, k, P, binomial_distribution::jeffreys_prior_interval)) + ); } } //