1 // Copyright 2018 Hans Dembinski 2 // 3 // Distributed under the Boost Software License, version 1.0. 4 // (See accompanying file LICENSE_1_0.txt 5 // or copy at http://www.boost.org/LICENSE_1_0.txt) 6 7 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_MEAN_HPP 8 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_MEAN_HPP 9 10 #include <boost/core/nvp.hpp> 11 #include <boost/histogram/fwd.hpp> // for weighted_mean<> 12 #include <boost/histogram/weight.hpp> 13 #include <cassert> 14 #include <type_traits> 15 16 namespace boost { 17 namespace histogram { 18 namespace accumulators { 19 20 /** 21 Calculates mean and variance of weighted sample. 22 23 Uses West's incremental algorithm to improve numerical stability 24 of mean and variance computation. 25 */ 26 template <class ValueType> 27 class weighted_mean { 28 public: 29 using value_type = ValueType; 30 using const_reference = const value_type&; 31 32 weighted_mean() = default; 33 34 /// Allow implicit conversion from other weighted_means 35 template <class T> weighted_mean(const weighted_mean<T> & o)36 weighted_mean(const weighted_mean<T>& o) 37 : sum_of_weights_{o.sum_of_weights_} 38 , sum_of_weights_squared_{o.sum_of_weights_squared_} 39 , weighted_mean_{o.weighted_mean_} 40 , sum_of_weighted_deltas_squared_{o.sum_of_weighted_deltas_squared_} {} 41 42 /// Initialize to external sum of weights, sum of weights squared, mean, and variance weighted_mean(const_reference wsum,const_reference wsum2,const_reference mean,const_reference variance)43 weighted_mean(const_reference wsum, const_reference wsum2, const_reference mean, 44 const_reference variance) 45 : sum_of_weights_(wsum) 46 , sum_of_weights_squared_(wsum2) 47 , weighted_mean_(mean) 48 , sum_of_weighted_deltas_squared_( 49 variance * (sum_of_weights_ - sum_of_weights_squared_ / sum_of_weights_)) {} 50 51 /// Insert sample x operator ()(const_reference x)52 void operator()(const_reference x) { operator()(weight(1), x); } 53 54 /// Insert sample x with weight w operator ()(const weight_type<value_type> & w,const_reference x)55 void operator()(const weight_type<value_type>& w, const_reference x) { 56 sum_of_weights_ += w.value; 57 sum_of_weights_squared_ += w.value * w.value; 58 const auto delta = x - weighted_mean_; 59 weighted_mean_ += w.value * delta / sum_of_weights_; 60 sum_of_weighted_deltas_squared_ += w.value * delta * (x - weighted_mean_); 61 } 62 63 /// Add another weighted_mean operator +=(const weighted_mean & rhs)64 weighted_mean& operator+=(const weighted_mean& rhs) { 65 if (sum_of_weights_ != 0 || rhs.sum_of_weights_ != 0) { 66 const auto tmp = 67 weighted_mean_ * sum_of_weights_ + rhs.weighted_mean_ * rhs.sum_of_weights_; 68 sum_of_weights_ += rhs.sum_of_weights_; 69 sum_of_weights_squared_ += rhs.sum_of_weights_squared_; 70 weighted_mean_ = tmp / sum_of_weights_; 71 } 72 sum_of_weighted_deltas_squared_ += rhs.sum_of_weighted_deltas_squared_; 73 return *this; 74 } 75 76 /** Scale by value 77 78 This acts as if all samples were scaled by the value. 79 */ operator *=(const_reference s)80 weighted_mean& operator*=(const_reference s) { 81 weighted_mean_ *= s; 82 sum_of_weighted_deltas_squared_ *= s * s; 83 return *this; 84 } 85 operator ==(const weighted_mean & rhs) const86 bool operator==(const weighted_mean& rhs) const noexcept { 87 return sum_of_weights_ == rhs.sum_of_weights_ && 88 sum_of_weights_squared_ == rhs.sum_of_weights_squared_ && 89 weighted_mean_ == rhs.weighted_mean_ && 90 sum_of_weighted_deltas_squared_ == rhs.sum_of_weighted_deltas_squared_; 91 } 92 operator !=(const weighted_mean & rhs) const93 bool operator!=(const weighted_mean& rhs) const noexcept { return !operator==(rhs); } 94 95 /// Return sum of weights sum_of_weights() const96 const_reference sum_of_weights() const noexcept { return sum_of_weights_; } 97 98 /// Return sum of weights squared (variance of weight distribution) sum_of_weights_squared() const99 const_reference sum_of_weights_squared() const noexcept { 100 return sum_of_weights_squared_; 101 } 102 103 /// Return mean of accumulated weighted samples value() const104 const_reference value() const noexcept { return weighted_mean_; } 105 106 /// Return variance of accumulated weighted samples variance() const107 value_type variance() const { 108 return sum_of_weighted_deltas_squared_ / 109 (sum_of_weights_ - sum_of_weights_squared_ / sum_of_weights_); 110 } 111 112 template <class Archive> serialize(Archive & ar,unsigned)113 void serialize(Archive& ar, unsigned /* version */) { 114 ar& make_nvp("sum_of_weights", sum_of_weights_); 115 ar& make_nvp("sum_of_weights_squared", sum_of_weights_squared_); 116 ar& make_nvp("weighted_mean", weighted_mean_); 117 ar& make_nvp("sum_of_weighted_deltas_squared", sum_of_weighted_deltas_squared_); 118 } 119 120 private: 121 value_type sum_of_weights_{}; 122 value_type sum_of_weights_squared_{}; 123 value_type weighted_mean_{}; 124 value_type sum_of_weighted_deltas_squared_{}; 125 }; 126 127 } // namespace accumulators 128 } // namespace histogram 129 } // namespace boost 130 131 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED 132 namespace std { 133 template <class T, class U> 134 /// Specialization for boost::histogram::accumulators::weighted_mean. 135 struct common_type<boost::histogram::accumulators::weighted_mean<T>, 136 boost::histogram::accumulators::weighted_mean<U>> { 137 using type = boost::histogram::accumulators::weighted_mean<common_type_t<T, U>>; 138 }; 139 } // namespace std 140 #endif 141 142 #endif 143