histogram/utility/jeffreys_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_JEFFREYS_INTERVAL_HPP
0008 #define BOOST_HISTOGRAM_UTILITY_JEFFREYS_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   Jeffreys interval.
0021
0022   This is the Bayesian credible interval with a Jeffreys prior. Although it has a
0023   Bayesian derivation, it has good coverage. The interval boundaries are close to the
0024   Wilson interval. A special property of this interval is that it is equal-tailed; the
0025   probability of the true value to be above or below the interval is approximately equal.
0026
0027   To avoid coverage probability tending to zero when the fraction approaches 0 or 1,
0028   this implementation uses a modification described in section 4.1.2 of the
0029   paper by L.D. Brown, T.T. Cai, A. DasGupta, Statistical Science 16 (2001) 101-133,
0030   doi:10.1214/ss/1009213286.
0031 */
0032 template <class ValueType>
0033 class jeffreys_interval : public binomial_proportion_interval<ValueType> {
0034 public:
0035   using value_type = typename jeffreys_interval::value_type;
0036   using interval_type = typename jeffreys_interval::interval_type;
0037
0038   /** Construct Jeffreys interval computer.
0039
0040     @param cl Confidence level for the interval. The default value produces a
0041     confidence level of 68 % equivalent to one standard deviation. Both `deviation` and
0042     `confidence_level` objects can be used to initialize the interval.
0043   */
0044   explicit jeffreys_interval(confidence_level cl = deviation{1}) noexcept
0045       : alpha_half_{static_cast<value_type>(0.5 - 0.5 * static_cast<double>(cl))} {}
0046
0047   using binomial_proportion_interval<ValueType>::operator();
0048
0049   /** Compute interval for given number of successes and failures.
0050
0051     @param successes Number of successful trials.
0052     @param failures Number of failed trials.
0053   */
0054   interval_type operator()(value_type successes, value_type failures) const noexcept {
0055     // See L.D. Brown, T.T. Cai, A. DasGupta, Statistical Science 16 (2001) 101-133,
0056     // doi:10.1214/ss/1009213286, section 4.1.2.
0057     const value_type half{0.5};
0058     const value_type total = successes + failures;
0059
0060     // if successes or failures are 0, modified interval is equal to Clopper-Pearson
0061     if (successes == 0) return {0, 1 - std::pow(alpha_half_, 1 / total)};
0062     if (failures == 0) return {std::pow(alpha_half_, 1 / total), 1};
0063
0064     math::beta_distribution<value_type> beta(successes + half, failures + half);
0065     const value_type a = successes == 1 ? 0 : math::quantile(beta, alpha_half_);
0066     const value_type b = failures == 1 ? 1 : math::quantile(beta, 1 - alpha_half_);
0067     return {a, b};
0068   }
0069
0070 private:
0071   value_type alpha_half_;
0072 };
0073
0074 } // namespace utility
0075 } // namespace histogram
0076 } // namespace boost
0077
0078 #endif