histogram/utility/clopper_pearson_interval.hpp

0001 // Copyright 2022 Jay Gohil, Hans Dembinski
0002 //
0003 // Distributed under the Boost Software License, version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_UTILITY_CLOPPER_PEARSON_INTERVAL_HPP
0008 #define BOOST_HISTOGRAM_UTILITY_CLOPPER_PEARSON_INTERVAL_HPP
0009
0010 #include <boost/histogram/fwd.hpp>
0011 #include <boost/histogram/utility/binomial_proportion_interval.hpp>
0012 #include <boost/math/distributions/beta.hpp>
0013 #include <cmath>
0014
0015 namespace boost {
0016 namespace histogram {
0017 namespace utility {
0018
0019 /**
0020   Clopper-Pearson interval.
0021
0022   This is the classic frequentist interval obtained with the Neyman construction.
0023   It is therefore often called the 'exact' interval. It is guaranteed to have at least the
0024   requested confidence level for all values of the fraction.
0025
0026   The interval is wider than others that produce coverage closer to the expected
0027   confidence level over a random ensemble of factions. The Clopper-Pearson interval
0028   essentially always overcovers for such a random ensemble, which is undesirable in
0029   practice. The Clopper-Pearson interval is recommended when it is important to be
0030   conservative, but the Wilson interval should be preferred in most applications.
0031
0032   C. Clopper, E.S. Pearson (1934), Biometrika 26 (4): 404-413.
0033   doi:10.1093/biomet/26.4.404.
0034 */
0035 template <class ValueType>
0036 class clopper_pearson_interval : public binomial_proportion_interval<ValueType> {
0037 public:
0038   using value_type = typename clopper_pearson_interval::value_type;
0039   using interval_type = typename clopper_pearson_interval::interval_type;
0040
0041   /** Construct Clopper-Pearson interval computer.
0042
0043     @param cl Confidence level for the interval. The default value produces a
0044     confidence level of 68 % equivalent to one standard deviation. Both `deviation` and
0045     `confidence_level` objects can be used to initialize the interval.
0046   */
0047   explicit clopper_pearson_interval(confidence_level cl = deviation{1}) noexcept
0048       : alpha_half_{static_cast<value_type>(0.5 - 0.5 * static_cast<double>(cl))} {}
0049
0050   using binomial_proportion_interval<ValueType>::operator();
0051
0052   /** Compute interval for given number of successes and failures.
0053
0054     @param successes Number of successful trials.
0055     @param failures Number of failed trials.
0056   */
0057   interval_type operator()(value_type successes, value_type failures) const noexcept {
0058     // analytical solution when successes or failures are zero
0059     // T. Mans (2014), Electronic Journal of Statistics. 8 (1): 817-840.
0060     // arXiv:1303.1288. doi:10.1214/14-EJS909.
0061     const value_type total = successes + failures;
0062     if (successes == 0) return {0, 1 - std::pow(alpha_half_, 1 / total)};
0063     if (failures == 0) return {std::pow(alpha_half_, 1 / total), 1};
0064
0065     // Source:
0066     // https://en.wikipedia.org/wiki/
0067     //   Binomial_proportion_confidence_interval#Clopper%E2%80%93Pearson_interval
0068     math::beta_distribution<value_type> beta_a(successes, failures + 1);
0069     const value_type a = math::quantile(beta_a, alpha_half_);
0070     math::beta_distribution<value_type> beta_b(successes + 1, failures);
0071     const value_type b = math::quantile(beta_b, 1 - alpha_half_);
0072     return {a, b};
0073   }
0074
0075 private:
0076   value_type alpha_half_;
0077 };
0078
0079 } // namespace utility
0080 } // namespace histogram
0081 } // namespace boost
0082
0083 #endif