mean.hpp 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179
  1. // Copyright 2015-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. #ifndef BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
  7. #define BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
  8. #include <boost/core/nvp.hpp>
  9. #include <boost/histogram/detail/square.hpp>
  10. #include <boost/histogram/fwd.hpp> // for mean<>
  11. #include <boost/throw_exception.hpp>
  12. #include <cassert>
  13. #include <stdexcept>
  14. #include <type_traits>
  15. namespace boost {
  16. namespace histogram {
  17. namespace accumulators {
  18. /** Calculates mean and variance of sample.
  19. Uses Welfords's incremental algorithm to improve the numerical
  20. stability of mean and variance computation.
  21. */
  22. template <class ValueType>
  23. class mean {
  24. public:
  25. using value_type = ValueType;
  26. using const_reference = const value_type&;
  27. mean() = default;
  28. /// Allow implicit conversion from mean<T>.
  29. template <class T>
  30. mean(const mean<T>& o) noexcept
  31. : sum_{o.sum_}, mean_{o.mean_}, sum_of_deltas_squared_{o.sum_of_deltas_squared_} {}
  32. /// Initialize to external count, mean, and variance.
  33. mean(const_reference n, const_reference mean, const_reference variance) noexcept
  34. : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (n - 1)) {}
  35. /// Insert sample x.
  36. void operator()(const_reference x) noexcept {
  37. sum_ += static_cast<value_type>(1);
  38. const auto delta = x - mean_;
  39. mean_ += delta / sum_;
  40. sum_of_deltas_squared_ += delta * (x - mean_);
  41. }
  42. /// Insert sample x with weight w.
  43. void operator()(const weight_type<value_type>& w, const_reference x) noexcept {
  44. sum_ += w.value;
  45. const auto delta = x - mean_;
  46. mean_ += w.value * delta / sum_;
  47. sum_of_deltas_squared_ += w.value * delta * (x - mean_);
  48. }
  49. /// Add another mean accumulator.
  50. mean& operator+=(const mean& rhs) noexcept {
  51. if (rhs.sum_ == 0) return *this;
  52. /*
  53. sum_of_deltas_squared
  54. = sum_i (x_i - mu)^2
  55. = sum_i (x_i - mu)^2 + sum_k (x_k - mu)^2
  56. = sum_i (x_i - mu1 + (mu1 - mu))^2 + sum_k (x_k - mu2 + (mu2 - mu))^2
  57. first part:
  58. sum_i (x_i - mu1 + (mu1 - mu))^2
  59. = sum_i (x_i - mu1)^2 + n1 (mu1 - mu))^2 + 2 (mu1 - mu) sum_i (x_i - mu1)
  60. = sum_i (x_i - mu1)^2 + n1 (mu1 - mu))^2
  61. since sum_i (x_i - mu1) = n1 mu1 - n1 mu1 = 0
  62. Putting it together:
  63. sum_of_deltas_squared
  64. = sum_of_deltas_squared_1 + n1 (mu1 - mu))^2
  65. + sum_of_deltas_squared_2 + n2 (mu2 - mu))^2
  66. */
  67. const auto n1 = sum_;
  68. const auto mu1 = mean_;
  69. const auto n2 = rhs.sum_;
  70. const auto mu2 = rhs.mean_;
  71. sum_ += rhs.sum_;
  72. mean_ = (n1 * mu1 + n2 * mu2) / sum_;
  73. sum_of_deltas_squared_ += rhs.sum_of_deltas_squared_;
  74. sum_of_deltas_squared_ += n1 * detail::square(mean_ - mu1);
  75. sum_of_deltas_squared_ += n2 * detail::square(mean_ - mu2);
  76. return *this;
  77. }
  78. /** Scale by value.
  79. This acts as if all samples were scaled by the value.
  80. */
  81. mean& operator*=(const_reference s) noexcept {
  82. mean_ *= s;
  83. sum_of_deltas_squared_ *= s * s;
  84. return *this;
  85. }
  86. bool operator==(const mean& rhs) const noexcept {
  87. return sum_ == rhs.sum_ && mean_ == rhs.mean_ &&
  88. sum_of_deltas_squared_ == rhs.sum_of_deltas_squared_;
  89. }
  90. bool operator!=(const mean& rhs) const noexcept { return !operator==(rhs); }
  91. /// Return how many samples were accumulated.
  92. const_reference count() const noexcept { return sum_; }
  93. /** Return mean value of accumulated samples.
  94. The result is undefined, if `count() < 1`.
  95. */
  96. const_reference value() const noexcept { return mean_; }
  97. /** Return variance of accumulated samples.
  98. The result is undefined, if `count() < 2`.
  99. */
  100. value_type variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); }
  101. template <class Archive>
  102. void serialize(Archive& ar, unsigned version) {
  103. if (version == 0) {
  104. // read only
  105. std::size_t sum;
  106. ar& make_nvp("sum", sum);
  107. sum_ = static_cast<value_type>(sum);
  108. } else {
  109. ar& make_nvp("sum", sum_);
  110. }
  111. ar& make_nvp("mean", mean_);
  112. ar& make_nvp("sum_of_deltas_squared", sum_of_deltas_squared_);
  113. }
  114. private:
  115. value_type sum_{};
  116. value_type mean_{};
  117. value_type sum_of_deltas_squared_{};
  118. };
  119. } // namespace accumulators
  120. } // namespace histogram
  121. } // namespace boost
  122. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  123. namespace boost {
  124. namespace serialization {
  125. template <class T>
  126. struct version;
  127. // version 1 for boost::histogram::accumulators::mean<T>
  128. template <class T>
  129. struct version<boost::histogram::accumulators::mean<T>> : std::integral_constant<int, 1> {
  130. };
  131. } // namespace serialization
  132. } // namespace boost
  133. namespace std {
  134. template <class T, class U>
  135. /// Specialization for boost::histogram::accumulators::mean.
  136. struct common_type<boost::histogram::accumulators::mean<T>,
  137. boost::histogram::accumulators::mean<U>> {
  138. using type = boost::histogram::accumulators::mean<common_type_t<T, U>>;
  139. };
  140. } // namespace std
  141. #endif
  142. #endif