histogram/utility/binomial_proportion_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_BINOMIAL_PROPORTION_INTERVAL_HPP
0008 #define BOOST_HISTOGRAM_UTILITY_BINOMIAL_PROPORTION_INTERVAL_HPP
0009
0010 #include <boost/histogram/detail/normal.hpp>
0011 #include <boost/histogram/fwd.hpp>
0012 #include <boost/throw_exception.hpp>
0013 #include <cmath>
0014 #include <stdexcept>
0015 #include <type_traits>
0016
0017 namespace boost {
0018 namespace histogram {
0019 namespace utility {
0020
0021 /**
0022   Common base class for interval calculators.
0023 */
0024 template <class ValueType>
0025 class binomial_proportion_interval {
0026   static_assert(std::is_floating_point<ValueType>::value,
0027                 "Value must be a floating point!");
0028
0029 public:
0030   using value_type = ValueType;
0031   using interval_type = std::pair<value_type, value_type>;
0032
0033   /** Compute interval for given number of successes and failures.
0034
0035     @param successes Number of successful trials.
0036     @param failures Number of failed trials.
0037   */
0038   virtual interval_type operator()(value_type successes,
0039                                    value_type failures) const noexcept = 0;
0040
0041   /** Compute interval for a fraction accumulator.
0042
0043     @param fraction Fraction accumulator.
0044   */
0045   template <class T>
0046   interval_type operator()(const accumulators::fraction<T>& fraction) const noexcept {
0047     return operator()(fraction.successes(), fraction.failures());
0048   }
0049 };
0050
0051 class deviation;
0052 class confidence_level;
0053
0054 /** Confidence level in units of deviations for intervals.
0055
0056   Intervals become wider as the deviation value increases. The standard deviation
0057   corresponds to a value of 1 and corresponds to 68.3 % confidence level. The conversion
0058   between confidence level and deviations is based on a two-sided interval on the normal
0059   distribution.
0060  */
0061 class deviation {
0062 public:
0063   /// constructor from units of standard deviations
0064   explicit deviation(double d) : d_{d} {
0065     if (d <= 0)
0066       BOOST_THROW_EXCEPTION(std::invalid_argument("scaling factor must be positive"));
0067   }
0068
0069   /// explicit conversion to units of standard deviations
0070   template <class T, class = std::enable_if_t<std::is_floating_point<T>::value>>
0071   explicit operator T() const noexcept {
0072     return static_cast<T>(d_);
0073   }
0074
0075   /// implicit conversion to confidence level
0076   operator confidence_level() const noexcept; // need to implement confidence_level first
0077
0078   friend deviation operator*(deviation d, double z) noexcept {
0079     return deviation(d.d_ * z);
0080   }
0081   friend deviation operator*(double z, deviation d) noexcept { return d * z; }
0082   friend bool operator==(deviation a, deviation b) noexcept { return a.d_ == b.d_; }
0083   friend bool operator!=(deviation a, deviation b) noexcept { return !(a == b); }
0084
0085 private:
0086   double d_;
0087 };
0088
0089 /** Confidence level for intervals.
0090
0091   Intervals become wider as the deviation value increases.
0092  */
0093 class confidence_level {
0094 public:
0095   /// constructor from confidence level (a probability)
0096   explicit confidence_level(double cl) : cl_{cl} {
0097     if (cl <= 0 || cl >= 1)
0098       BOOST_THROW_EXCEPTION(std::invalid_argument("0 < cl < 1 is required"));
0099   }
0100
0101   /// explicit conversion to numerical confidence level
0102   template <class T, class = std::enable_if_t<std::is_floating_point<T>::value>>
0103   explicit operator T() const noexcept {
0104     return static_cast<T>(cl_);
0105   }
0106
0107   /// implicit conversion to units of standard deviation
0108   operator deviation() const noexcept {
0109     return deviation{detail::normal_ppf(std::fma(0.5, cl_, 0.5))};
0110   }
0111
0112   friend bool operator==(confidence_level a, confidence_level b) noexcept {
0113     return a.cl_ == b.cl_;
0114   }
0115   friend bool operator!=(confidence_level a, confidence_level b) noexcept {
0116     return !(a == b);
0117   }
0118
0119 private:
0120   double cl_;
0121 };
0122
0123 inline deviation::operator confidence_level() const noexcept {
0124   // solve normal cdf(z) - cdf(-z) = 2 (cdf(z) - 0.5)
0125   return confidence_level{std::fma(2.0, detail::normal_cdf(d_), -1.0)};
0126 }
0127
0128 } // namespace utility
0129 } // namespace histogram
0130 } // namespace boost
0131
0132 #endif