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}

rtc_base/BUILD.gn

@@ -577,6 +577,8 @@ rtc_library("weak_ptr") {
 
 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 @@ if (rtc_include_tests) {
     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",

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

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&section=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_

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