diff --git a/rtc_base/BUILD.gn b/rtc_base/BUILD.gn index 0aee6d1165..a7f0c9ee75 100644 --- a/rtc_base/BUILD.gn +++ b/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", 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 0000000000..18242bd5f9 --- /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 + +#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(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::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 0000000000..a72aa271ef --- /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 +#include +#include +#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::infinity(); + // This is the ratio between variance of the estimate and variance of samples. + double estimator_variance_ = 1; + absl::optional 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 0000000000..53b094e10e --- /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 + +#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::infinity(), average.GetVariance()); + EXPECT_EQ(std::numeric_limits::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::infinity(), average.GetVariance()); + EXPECT_EQ(std::numeric_limits::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