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