mean.hpp 2.5 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091
  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/histogram/fwd.hpp>
  9. #include <cstddef>
  10. #include <type_traits>
  11. namespace boost {
  12. namespace histogram {
  13. namespace accumulators {
  14. /** Calculates mean and variance of sample.
  15. Uses Welfords's incremental algorithm to improve the numerical
  16. stability of mean and variance computation.
  17. */
  18. template <class RealType>
  19. class mean {
  20. public:
  21. mean() = default;
  22. mean(const std::size_t n, const RealType& mean, const RealType& variance)
  23. : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (sum_ - 1)) {}
  24. void operator()(const RealType& x) {
  25. sum_ += 1;
  26. const auto delta = x - mean_;
  27. mean_ += delta / sum_;
  28. sum_of_deltas_squared_ += delta * (x - mean_);
  29. }
  30. template <class T>
  31. mean& operator+=(const mean<T>& rhs) {
  32. const auto tmp = mean_ * sum_ + static_cast<RealType>(rhs.mean_ * rhs.sum_);
  33. sum_ += rhs.sum_;
  34. mean_ = tmp / sum_;
  35. sum_of_deltas_squared_ += static_cast<RealType>(rhs.sum_of_deltas_squared_);
  36. return *this;
  37. }
  38. mean& operator*=(const RealType& s) {
  39. mean_ *= s;
  40. sum_of_deltas_squared_ *= s * s;
  41. return *this;
  42. }
  43. template <class T>
  44. bool operator==(const mean<T>& rhs) const noexcept {
  45. return sum_ == rhs.sum_ && mean_ == rhs.mean_ &&
  46. sum_of_deltas_squared_ == rhs.sum_of_deltas_squared_;
  47. }
  48. template <class T>
  49. bool operator!=(const mean<T>& rhs) const noexcept {
  50. return !operator==(rhs);
  51. }
  52. std::size_t count() const noexcept { return sum_; }
  53. const RealType& value() const noexcept { return mean_; }
  54. RealType variance() const { return sum_of_deltas_squared_ / (sum_ - 1); }
  55. template <class Archive>
  56. void serialize(Archive&, unsigned /* version */);
  57. private:
  58. std::size_t sum_ = 0;
  59. RealType mean_ = 0, sum_of_deltas_squared_ = 0;
  60. };
  61. } // namespace accumulators
  62. } // namespace histogram
  63. } // namespace boost
  64. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  65. namespace std {
  66. template <class T, class U>
  67. /// Specialization for boost::histogram::accumulators::mean.
  68. struct common_type<boost::histogram::accumulators::mean<T>,
  69. boost::histogram::accumulators::mean<U>> {
  70. using type = boost::histogram::accumulators::mean<common_type_t<T, U>>;
  71. };
  72. } // namespace std
  73. #endif
  74. #endif