1 /* 2 * Copyright (c) 2019 The WebRTC project authors. All Rights Reserved. 3 * 4 * Use of this source code is governed by a BSD-style license 5 * that can be found in the LICENSE file in the root of the source 6 * tree. An additional intellectual property rights grant can be found 7 * in the file PATENTS. All contributing project authors may 8 * be found in the AUTHORS file in the root of the source tree. 9 */ 10 11 #ifndef RTC_BASE_NUMERICS_RUNNING_STATISTICS_H_ 12 #define RTC_BASE_NUMERICS_RUNNING_STATISTICS_H_ 13 14 #include <algorithm> 15 #include <cmath> 16 #include <limits> 17 18 #include "absl/types/optional.h" 19 #include "rtc_base/checks.h" 20 #include "rtc_base/numerics/math_utils.h" 21 22 namespace webrtc { 23 24 // tl;dr: Robust and efficient online computation of statistics, 25 // using Welford's method for variance. [1] 26 // 27 // This should be your go-to class if you ever need to compute 28 // min, max, mean, variance and standard deviation. 29 // If you need to get percentiles, please use webrtc::SamplesStatsCounter. 30 // 31 // Please note RemoveSample() won't affect min and max. 32 // If you want a full-fledged moving window over N last samples, 33 // please use webrtc::RollingAccumulator. 34 // 35 // The measures return absl::nullopt if no samples were fed (Size() == 0), 36 // otherwise the returned optional is guaranteed to contain a value. 37 // 38 // [1] 39 // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm 40 41 // The type T is a scalar which must be convertible to double. 42 // Rationale: we often need greater precision for measures 43 // than for the samples themselves. 44 template <typename T> 45 class RunningStatistics { 46 public: 47 // Update stats //////////////////////////////////////////// 48 49 // Add a value participating in the statistics in O(1) time. AddSample(T sample)50 void AddSample(T sample) { 51 max_ = std::max(max_, sample); 52 min_ = std::min(min_, sample); 53 ++size_; 54 // Welford's incremental update. 55 const double delta = sample - mean_; 56 mean_ += delta / size_; 57 const double delta2 = sample - mean_; 58 cumul_ += delta * delta2; 59 } 60 61 // Remove a previously added value in O(1) time. 62 // Nb: This doesn't affect min or max. 63 // Calling RemoveSample when Size()==0 is incorrect. RemoveSample(T sample)64 void RemoveSample(T sample) { 65 RTC_DCHECK_GT(Size(), 0); 66 // In production, just saturate at 0. 67 if (Size() == 0) { 68 return; 69 } 70 // Since samples order doesn't matter, this is the 71 // exact reciprocal of Welford's incremental update. 72 --size_; 73 const double delta = sample - mean_; 74 mean_ -= delta / size_; 75 const double delta2 = sample - mean_; 76 cumul_ -= delta * delta2; 77 } 78 79 // Merge other stats, as if samples were added one by one, but in O(1). MergeStatistics(const RunningStatistics<T> & other)80 void MergeStatistics(const RunningStatistics<T>& other) { 81 if (other.size_ == 0) { 82 return; 83 } 84 max_ = std::max(max_, other.max_); 85 min_ = std::min(min_, other.min_); 86 const int64_t new_size = size_ + other.size_; 87 const double new_mean = 88 (mean_ * size_ + other.mean_ * other.size_) / new_size; 89 // Each cumulant must be corrected. 90 // * from: sum((x_i - mean_)²) 91 // * to: sum((x_i - new_mean)²) 92 auto delta = [new_mean](const RunningStatistics<T>& stats) { 93 return stats.size_ * (new_mean * (new_mean - 2 * stats.mean_) + 94 stats.mean_ * stats.mean_); 95 }; 96 cumul_ = cumul_ + delta(*this) + other.cumul_ + delta(other); 97 mean_ = new_mean; 98 size_ = new_size; 99 } 100 101 // Get Measures //////////////////////////////////////////// 102 103 // Returns number of samples involved via AddSample() or MergeStatistics(), 104 // minus number of times RemoveSample() was called. Size()105 int64_t Size() const { return size_; } 106 107 // Returns minimum among all seen samples, in O(1) time. 108 // This isn't affected by RemoveSample(). GetMin()109 absl::optional<T> GetMin() const { 110 if (size_ == 0) { 111 return absl::nullopt; 112 } 113 return min_; 114 } 115 116 // Returns maximum among all seen samples, in O(1) time. 117 // This isn't affected by RemoveSample(). GetMax()118 absl::optional<T> GetMax() const { 119 if (size_ == 0) { 120 return absl::nullopt; 121 } 122 return max_; 123 } 124 125 // Returns mean in O(1) time. GetMean()126 absl::optional<double> GetMean() const { 127 if (size_ == 0) { 128 return absl::nullopt; 129 } 130 return mean_; 131 } 132 133 // Returns unbiased sample variance in O(1) time. GetVariance()134 absl::optional<double> GetVariance() const { 135 if (size_ == 0) { 136 return absl::nullopt; 137 } 138 return cumul_ / size_; 139 } 140 141 // Returns unbiased standard deviation in O(1) time. GetStandardDeviation()142 absl::optional<double> GetStandardDeviation() const { 143 if (size_ == 0) { 144 return absl::nullopt; 145 } 146 return std::sqrt(*GetVariance()); 147 } 148 149 private: 150 int64_t size_ = 0; // Samples seen. 151 T min_ = infinity_or_max<T>(); 152 T max_ = minus_infinity_or_min<T>(); 153 double mean_ = 0; 154 double cumul_ = 0; // Variance * size_, sometimes noted m2. 155 }; 156 157 } // namespace webrtc 158 159 #endif // RTC_BASE_NUMERICS_RUNNING_STATISTICS_H_ 160