histogram/accumulators/weighted_mean.hpp

0001 // Copyright 2018 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_ACCUMULATORS_WEIGHTED_MEAN_HPP
0008 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_MEAN_HPP
0009
0010 #include <boost/core/nvp.hpp>
0011 #include <boost/histogram/detail/square.hpp>
0012 #include <boost/histogram/fwd.hpp> // for weighted_mean<>
0013 #include <boost/histogram/weight.hpp>
0014 #include <cassert>
0015 #include <type_traits>
0016
0017 namespace boost {
0018 namespace histogram {
0019 namespace accumulators {
0020
0021 /**
0022   Calculates mean and variance of weighted sample.
0023
0024   Uses West's incremental algorithm to improve numerical stability
0025   of mean and variance computation.
0026 */
0027 template <class ValueType>
0028 class weighted_mean {
0029 public:
0030   using value_type = ValueType;
0031   using const_reference = const value_type&;
0032
0033   weighted_mean() = default;
0034
0035   /// Allow implicit conversion from other weighted_means.
0036   template <class T>
0037   weighted_mean(const weighted_mean<T>& o)
0038       : sum_of_weights_{o.sum_of_weights_}
0039       , sum_of_weights_squared_{o.sum_of_weights_squared_}
0040       , weighted_mean_{o.weighted_mean_}
0041       , sum_of_weighted_deltas_squared_{o.sum_of_weighted_deltas_squared_} {}
0042
0043   /// Initialize to external sum of weights, sum of weights squared, mean, and variance.
0044   weighted_mean(const_reference wsum, const_reference wsum2, const_reference mean,
0045                 const_reference variance)
0046       : sum_of_weights_(wsum)
0047       , sum_of_weights_squared_(wsum2)
0048       , weighted_mean_(mean)
0049       , sum_of_weighted_deltas_squared_(
0050             variance * (sum_of_weights_ - sum_of_weights_squared_ / sum_of_weights_)) {}
0051
0052   /// Insert sample x.
0053   void operator()(const_reference x) { operator()(weight(1), x); }
0054
0055   /// Insert sample x with weight w.
0056   void operator()(const weight_type<value_type>& w, const_reference x) {
0057     sum_of_weights_ += w.value;
0058     sum_of_weights_squared_ += w.value * w.value;
0059     const auto delta = x - weighted_mean_;
0060     weighted_mean_ += w.value * delta / sum_of_weights_;
0061     sum_of_weighted_deltas_squared_ += w.value * delta * (x - weighted_mean_);
0062   }
0063
0064   /// Add another weighted_mean.
0065   weighted_mean& operator+=(const weighted_mean& rhs) {
0066     if (rhs.sum_of_weights_ == 0) return *this;
0067
0068     // see mean.hpp for derivation of correct formula
0069
0070     const auto n1 = sum_of_weights_;
0071     const auto mu1 = weighted_mean_;
0072     const auto n2 = rhs.sum_of_weights_;
0073     const auto mu2 = rhs.weighted_mean_;
0074
0075     sum_of_weights_ += rhs.sum_of_weights_;
0076     sum_of_weights_squared_ += rhs.sum_of_weights_squared_;
0077     weighted_mean_ = (n1 * mu1 + n2 * mu2) / sum_of_weights_;
0078
0079     sum_of_weighted_deltas_squared_ += rhs.sum_of_weighted_deltas_squared_;
0080     sum_of_weighted_deltas_squared_ += n1 * detail::square(weighted_mean_ - mu1);
0081     sum_of_weighted_deltas_squared_ += n2 * detail::square(weighted_mean_ - mu2);
0082
0083     return *this;
0084   }
0085
0086   /** Scale by value.
0087
0088    This acts as if all samples were scaled by the value.
0089   */
0090   weighted_mean& operator*=(const_reference s) noexcept {
0091     weighted_mean_ *= s;
0092     sum_of_weighted_deltas_squared_ *= s * s;
0093     return *this;
0094   }
0095
0096   bool operator==(const weighted_mean& rhs) const noexcept {
0097     return sum_of_weights_ == rhs.sum_of_weights_ &&
0098            sum_of_weights_squared_ == rhs.sum_of_weights_squared_ &&
0099            weighted_mean_ == rhs.weighted_mean_ &&
0100            sum_of_weighted_deltas_squared_ == rhs.sum_of_weighted_deltas_squared_;
0101   }
0102
0103   bool operator!=(const weighted_mean& rhs) const noexcept { return !operator==(rhs); }
0104
0105   /// Return sum of weights.
0106   const_reference sum_of_weights() const noexcept { return sum_of_weights_; }
0107
0108   /// Return sum of weights squared (variance of weight distribution).
0109   const_reference sum_of_weights_squared() const noexcept {
0110     return sum_of_weights_squared_;
0111   }
0112
0113   /** Return effective counts.
0114
0115     This corresponds to the equivalent number of unweighted samples that would
0116     have the same variance as this sample. count() should be used to check whether
0117     value() and variance() are defined, see documentation of value() and variance().
0118     count() can be used to compute the variance of the mean by dividing variance()
0119     by count().
0120   */
0121   value_type count() const noexcept {
0122     // see https://en.wikipedia.org/wiki/Effective_sample_size#weighted_samples
0123     return detail::square(sum_of_weights_) / sum_of_weights_squared_;
0124   }
0125
0126   /** Return mean value of accumulated weighted samples.
0127
0128     The result is undefined, if count() == 0.
0129   */
0130   const_reference value() const noexcept { return weighted_mean_; }
0131
0132   /** Return variance of accumulated weighted samples.
0133
0134     The result is undefined, if count() == 0 or count() == 1.
0135   */
0136   value_type variance() const {
0137     // see https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights
0138     return sum_of_weighted_deltas_squared_ /
0139            (sum_of_weights_ - sum_of_weights_squared_ / sum_of_weights_);
0140   }
0141
0142   template <class Archive>
0143   void serialize(Archive& ar, unsigned /* version */) {
0144     ar& make_nvp("sum_of_weights", sum_of_weights_);
0145     ar& make_nvp("sum_of_weights_squared", sum_of_weights_squared_);
0146     ar& make_nvp("weighted_mean", weighted_mean_);
0147     ar& make_nvp("sum_of_weighted_deltas_squared", sum_of_weighted_deltas_squared_);
0148   }
0149
0150 private:
0151   value_type sum_of_weights_{};
0152   value_type sum_of_weights_squared_{};
0153   value_type weighted_mean_{};
0154   value_type sum_of_weighted_deltas_squared_{};
0155 };
0156
0157 } // namespace accumulators
0158 } // namespace histogram
0159 } // namespace boost
0160
0161 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
0162 namespace std {
0163 template <class T, class U>
0164 /// Specialization for boost::histogram::accumulators::weighted_mean.
0165 struct common_type<boost::histogram::accumulators::weighted_mean<T>,
0166                    boost::histogram::accumulators::weighted_mean<U>> {
0167   using type = boost::histogram::accumulators::weighted_mean<common_type_t<T, U>>;
0168 };
0169 } // namespace std
0170 #endif
0171
0172 #endif