Add class for ExponentialMovingAverage
Bug: webrtc:11120
Change-Id: I210671e00276546e9d63b148385263cb1256e2b0
Reviewed-on: https://webrtc-review.googlesource.com/c/src/+/160307
Reviewed-by: Harald Alvestrand <hta@webrtc.org>
Reviewed-by: Niels Moller <nisse@webrtc.org>
Commit-Queue: Jonas Oreland <jonaso@webrtc.org>
Cr-Commit-Position: refs/heads/master@{#29901}
diff --git a/rtc_base/BUILD.gn b/rtc_base/BUILD.gn
index 0aee6d1..a7f0c9e 100644
--- a/rtc_base/BUILD.gn
+++ b/rtc_base/BUILD.gn
@@ -577,6 +577,8 @@
rtc_library("rtc_numerics") {
sources = [
+ "numerics/event_based_exponential_moving_average.cc",
+ "numerics/event_based_exponential_moving_average.h",
"numerics/exp_filter.cc",
"numerics/exp_filter.h",
"numerics/math_utils.h",
@@ -1288,6 +1290,7 @@
testonly = true
sources = [
+ "numerics/event_based_exponential_moving_average_unittest.cc",
"numerics/exp_filter_unittest.cc",
"numerics/moving_average_unittest.cc",
"numerics/moving_median_filter_unittest.cc",
diff --git a/rtc_base/numerics/event_based_exponential_moving_average.cc b/rtc_base/numerics/event_based_exponential_moving_average.cc
new file mode 100644
index 0000000..18242bd
--- /dev/null
+++ b/rtc_base/numerics/event_based_exponential_moving_average.cc
@@ -0,0 +1,66 @@
+/*
+ * Copyright 2019 The WebRTC Project Authors. All rights reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "rtc_base/numerics/event_based_exponential_moving_average.h"
+
+#include <cmath>
+
+#include "rtc_base/checks.h"
+
+namespace {
+
+// For a normal distributed value, the 95% double sided confidence interval is
+// is 1.96 * stddev.
+constexpr double ninetyfive_percent_confidence = 1.96;
+
+} // namespace
+
+namespace rtc {
+
+// |half_time| specifies how much weight will be given to old samples,
+// a sample gets exponentially less weight so that it's 50%
+// after |half_time| time units has passed.
+EventBasedExponentialMovingAverage::EventBasedExponentialMovingAverage(
+ int half_time)
+ : tau_(static_cast<double>(half_time) / log(2)) {}
+
+void EventBasedExponentialMovingAverage::AddSample(int64_t now, int sample) {
+ if (!last_observation_timestamp_.has_value()) {
+ value_ = sample;
+ } else {
+ RTC_DCHECK(now > *last_observation_timestamp_);
+ // Variance gets computed after second sample.
+ int64_t age = now - *last_observation_timestamp_;
+ double e = exp(-age / tau_);
+ double alpha = e / (1 + e);
+ double one_minus_alpha = 1 - alpha;
+ double sample_diff = sample - value_;
+ value_ = one_minus_alpha * value_ + alpha * sample;
+ estimator_variance_ =
+ (one_minus_alpha * one_minus_alpha) * estimator_variance_ +
+ (alpha * alpha);
+ if (sample_variance_ == std::numeric_limits<double>::infinity()) {
+ // First variance.
+ sample_variance_ = sample_diff * sample_diff;
+ } else {
+ double new_variance = one_minus_alpha * sample_variance_ +
+ alpha * sample_diff * sample_diff;
+ sample_variance_ = new_variance;
+ }
+ }
+ last_observation_timestamp_ = now;
+}
+
+double EventBasedExponentialMovingAverage::GetConfidenceInterval() const {
+ return ninetyfive_percent_confidence *
+ sqrt(sample_variance_ * estimator_variance_);
+}
+
+} // namespace rtc
diff --git a/rtc_base/numerics/event_based_exponential_moving_average.h b/rtc_base/numerics/event_based_exponential_moving_average.h
new file mode 100644
index 0000000..a72aa27
--- /dev/null
+++ b/rtc_base/numerics/event_based_exponential_moving_average.h
@@ -0,0 +1,63 @@
+/*
+ * Copyright 2019 The WebRTC Project Authors. All rights reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef RTC_BASE_NUMERICS_EVENT_BASED_EXPONENTIAL_MOVING_AVERAGE_H_
+#define RTC_BASE_NUMERICS_EVENT_BASED_EXPONENTIAL_MOVING_AVERAGE_H_
+
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include "absl/types/optional.h"
+
+namespace rtc {
+
+/**
+ * This class implements exponential moving average for time series
+ * estimating both value, variance and variance of estimator based on
+ * https://en.wikipedia.org/w/index.php?title=Moving_average§ion=9#Application_to_measuring_computer_performance
+ * with the additions from nisse@ added to
+ * https://en.wikipedia.org/wiki/Talk:Moving_average.
+ *
+ * A sample gets exponentially less weight so that it's 50%
+ * after |half_time| time units.
+ */
+class EventBasedExponentialMovingAverage {
+ public:
+ // |half_time| specifies how much weight will be given to old samples,
+ // see example above.
+ explicit EventBasedExponentialMovingAverage(int half_time);
+
+ void AddSample(int64_t now, int value);
+
+ double GetAverage() const { return value_; }
+ double GetVariance() const { return sample_variance_; }
+
+ // Compute 95% confidence interval assuming that
+ // - variance of samples are normal distributed.
+ // - variance of estimator is normal distributed.
+ //
+ // The returned values specifies the distance from the average,
+ // i.e if X = GetAverage(), m = GetConfidenceInterval()
+ // then a there is 95% likelihood that the observed variables is inside
+ // [ X +/- m ].
+ double GetConfidenceInterval() const;
+
+ private:
+ const double tau_;
+ double value_ = std::nan("uninit");
+ double sample_variance_ = std::numeric_limits<double>::infinity();
+ // This is the ratio between variance of the estimate and variance of samples.
+ double estimator_variance_ = 1;
+ absl::optional<int64_t> last_observation_timestamp_;
+};
+
+} // namespace rtc
+
+#endif // RTC_BASE_NUMERICS_EVENT_BASED_EXPONENTIAL_MOVING_AVERAGE_H_
diff --git a/rtc_base/numerics/event_based_exponential_moving_average_unittest.cc b/rtc_base/numerics/event_based_exponential_moving_average_unittest.cc
new file mode 100644
index 0000000..53b094e
--- /dev/null
+++ b/rtc_base/numerics/event_based_exponential_moving_average_unittest.cc
@@ -0,0 +1,168 @@
+/*
+ * Copyright 2019 The WebRTC Project Authors. All rights reserved.
+ *
+ * Use of this source code is governed by a BSD-style license
+ * that can be found in the LICENSE file in the root of the source
+ * tree. An additional intellectual property rights grant can be found
+ * in the file PATENTS. All contributing project authors may
+ * be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include "rtc_base/numerics/event_based_exponential_moving_average.h"
+
+#include <cmath>
+
+#include "test/gtest.h"
+
+namespace {
+
+constexpr int kHalfTime = 500;
+constexpr double kError = 0.1;
+
+} // namespace
+
+namespace rtc {
+
+TEST(EventBasedExponentialMovingAverageTest, NoValue) {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+
+ EXPECT_TRUE(std::isnan(average.GetAverage()));
+ EXPECT_EQ(std::numeric_limits<double>::infinity(), average.GetVariance());
+ EXPECT_EQ(std::numeric_limits<double>::infinity(),
+ average.GetConfidenceInterval());
+}
+
+TEST(EventBasedExponentialMovingAverageTest, FirstValue) {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+
+ int64_t time = 23;
+ constexpr int value = 1000;
+ average.AddSample(time, value);
+ EXPECT_NEAR(value, average.GetAverage(), kError);
+ EXPECT_EQ(std::numeric_limits<double>::infinity(), average.GetVariance());
+ EXPECT_EQ(std::numeric_limits<double>::infinity(),
+ average.GetConfidenceInterval());
+}
+
+TEST(EventBasedExponentialMovingAverageTest, Half) {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+
+ int64_t time = 23;
+ constexpr int value = 1000;
+ average.AddSample(time, value);
+ average.AddSample(time + kHalfTime, 0);
+ EXPECT_NEAR(666.7, average.GetAverage(), kError);
+ EXPECT_NEAR(1000000, average.GetVariance(), kError);
+ EXPECT_NEAR(1460.9, average.GetConfidenceInterval(), kError);
+}
+
+TEST(EventBasedExponentialMovingAverageTest, Same) {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+
+ int64_t time = 23;
+ constexpr int value = 1000;
+ average.AddSample(time, value);
+ average.AddSample(time + kHalfTime, value);
+ EXPECT_NEAR(value, average.GetAverage(), kError);
+ EXPECT_NEAR(0, average.GetVariance(), kError);
+ EXPECT_NEAR(0, average.GetConfidenceInterval(), kError);
+}
+
+TEST(EventBasedExponentialMovingAverageTest, Almost100) {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+
+ int64_t time = 23;
+ constexpr int value = 100;
+ average.AddSample(time + 0 * kHalfTime, value - 10);
+ average.AddSample(time + 1 * kHalfTime, value + 10);
+ average.AddSample(time + 2 * kHalfTime, value - 15);
+ average.AddSample(time + 3 * kHalfTime, value + 15);
+ EXPECT_NEAR(100.2, average.GetAverage(), kError);
+ EXPECT_NEAR(372.6, average.GetVariance(), kError);
+ EXPECT_NEAR(19.7, average.GetConfidenceInterval(), kError); // 100 +/- 20
+
+ average.AddSample(time + 4 * kHalfTime, value);
+ average.AddSample(time + 5 * kHalfTime, value);
+ average.AddSample(time + 6 * kHalfTime, value);
+ average.AddSample(time + 7 * kHalfTime, value);
+ EXPECT_NEAR(100.0, average.GetAverage(), kError);
+ EXPECT_NEAR(73.6, average.GetVariance(), kError);
+ EXPECT_NEAR(7.6, average.GetConfidenceInterval(), kError); // 100 +/- 7
+}
+
+// Test that getting a value at X and another at X+1
+// is almost the same as getting another at X and a value at X+1.
+TEST(EventBasedExponentialMovingAverageTest, SameTime) {
+ int64_t time = 23;
+ constexpr int value = 100;
+
+ {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+ average.AddSample(time + 0, value);
+ average.AddSample(time + 1, 0);
+ EXPECT_NEAR(50, average.GetAverage(), kError);
+ EXPECT_NEAR(10000, average.GetVariance(), kError);
+ EXPECT_NEAR(138.6, average.GetConfidenceInterval(),
+ kError); // 50 +/- 138.6
+ }
+
+ {
+ EventBasedExponentialMovingAverage average(kHalfTime);
+ average.AddSample(time + 0, 0);
+ average.AddSample(time + 1, 100);
+ EXPECT_NEAR(50, average.GetAverage(), kError);
+ EXPECT_NEAR(10000, average.GetVariance(), kError);
+ EXPECT_NEAR(138.6, average.GetConfidenceInterval(),
+ kError); // 50 +/- 138.6
+ }
+}
+
+// This test shows behavior of estimator with a half_time of 100.
+// It is unclear if these set of observations are representative
+// of any real world scenarios.
+TEST(EventBasedExponentialMovingAverageTest, NonUniformSamplesHalftime100) {
+ int64_t time = 23;
+ constexpr int value = 100;
+
+ {
+ // The observations at 100 and 101, are significantly close in
+ // time that the estimator returns approx. the average.
+ EventBasedExponentialMovingAverage average(100);
+ average.AddSample(time + 0, value);
+ average.AddSample(time + 100, value);
+ average.AddSample(time + 101, 0);
+ EXPECT_NEAR(50.2, average.GetAverage(), kError);
+ EXPECT_NEAR(86.2, average.GetConfidenceInterval(), kError); // 50 +/- 86
+ }
+
+ {
+ EventBasedExponentialMovingAverage average(100);
+ average.AddSample(time + 0, value);
+ average.AddSample(time + 1, value);
+ average.AddSample(time + 100, 0);
+ EXPECT_NEAR(66.5, average.GetAverage(), kError);
+ EXPECT_NEAR(65.4, average.GetConfidenceInterval(), kError); // 66 +/- 65
+ }
+
+ {
+ EventBasedExponentialMovingAverage average(100);
+ for (int i = 0; i < 10; i++) {
+ average.AddSample(time + i, value);
+ }
+ average.AddSample(time + 100, 0);
+ EXPECT_NEAR(65.3, average.GetAverage(), kError);
+ EXPECT_NEAR(59.1, average.GetConfidenceInterval(), kError); // 55 +/- 59
+ }
+
+ {
+ EventBasedExponentialMovingAverage average(100);
+ average.AddSample(time + 0, 100);
+ for (int i = 90; i <= 100; i++) {
+ average.AddSample(time + i, 0);
+ }
+ EXPECT_NEAR(0.05, average.GetAverage(), kError);
+ EXPECT_NEAR(4.9, average.GetConfidenceInterval(), kError); // 0 +/- 5
+ }
+}
+
+} // namespace rtc
