histogram/utility/wilson_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_WILSON_INTERVAL_HPP
0008 #define BOOST_HISTOGRAM_UTILITY_WILSON_INTERVAL_HPP
0009
0010 #include <boost/histogram/fwd.hpp>
0011 #include <boost/histogram/utility/binomial_proportion_interval.hpp>
0012 #include <cmath>
0013 #include <utility>
0014
0015 namespace boost {
0016 namespace histogram {
0017 namespace utility {
0018
0019 /**
0020   Wilson interval.
0021
0022   The Wilson score interval is simple to compute, has good coverage. Intervals are
0023   automatically bounded between 0 and 1 and never empty. The interval is asymmetric.
0024
0025   Wilson, E. B. (1927). "Probable inference, the law of succession, and statistical
0026   inference". Journal of the American Statistical Association. 22 (158): 209-212.
0027   doi:10.1080/01621459.1927.10502953. JSTOR 2276774.
0028
0029   The coverage probability for a random ensemble of fractions is close to the nominal
0030   value. Unlike the Clopper-Pearson interval, the Wilson score interval is not
0031   conservative. For some values of the fractions, the interval undercovers and overcovers
0032   for neighboring values. This is a shared property of all alternatives to the
0033   Clopper-Pearson interval.
0034
0035   The Wilson score intervals is widely recommended for general use in the literature. For
0036   a review of the literature, see R. D. Cousins, K. E. Hymes, J. Tucker, Nucl. Instrum.
0037   Meth. A 612 (2010) 388-398.
0038 */
0039 template <class ValueType>
0040 class wilson_interval : public binomial_proportion_interval<ValueType> {
0041 public:
0042   using value_type = typename wilson_interval::value_type;
0043   using interval_type = typename wilson_interval::interval_type;
0044
0045   /** Construct Wilson interval computer.
0046
0047     @param d Number of standard deviations for the interval. The default value 1
0048     corresponds to a confidence level of 68 %. Both `deviation` and `confidence_level`
0049     objects can be used to initialize the interval.
0050   */
0051   explicit wilson_interval(deviation d = deviation{1.0}) noexcept
0052       : z_{static_cast<value_type>(d)} {}
0053
0054   using binomial_proportion_interval<ValueType>::operator();
0055
0056   /** Compute interval for given number of successes and failures.
0057
0058     @param successes Number of successful trials.
0059     @param failures Number of failed trials.
0060   */
0061   interval_type operator()(value_type successes, value_type failures) const noexcept {
0062     // See https://en.wikipedia.org/wiki/
0063     //   Binomial_proportion_confidence_interval
0064     //   #Wilson_score_interval
0065
0066     // We make sure calculation is done in single precision if value_type is float
0067     // by converting all literals to value_type. Double literals in the equation
0068     // would turn intermediate values to double.
0069     const value_type half{0.5}, quarter{0.25}, zsq{z_ * z_};
0070     const value_type total = successes + failures;
0071     const value_type minv = 1 / (total + zsq);
0072     const value_type t1 = (successes + half * zsq) * minv;
0073     const value_type t2 =
0074         z_ * minv * std::sqrt(successes * failures / total + quarter * zsq);
0075     return {t1 - t2, t1 + t2};
0076   }
0077
0078 private:
0079   value_type z_;
0080 };
0081
0082 } // namespace utility
0083 } // namespace histogram
0084 } // namespace boost
0085
0086 #endif