[ML] Improve robustness to very low variance data (#232)

Author: tveasey

include/maths/CPeriodicityHypothesisTests.h

@@ -208,7 +208,7 @@ class MATHS_EXPORT CPeriodicityHypothesisTests {
 
     //! \brief A collection of statistics used during testing.
     struct STestStats {
-        STestStats();
+        explicit STestStats(double meanMagnitude);
         //! Set the various test thresholds.
         void setThresholds(double vt, double at, double Rt);
         //! Check if the null hypothesis is good enough to not need an
@@ -233,6 +233,8 @@ class MATHS_EXPORT CPeriodicityHypothesisTests {
         double s_NonEmptyBuckets;
         //! The average number of measurements per bucket value.
         double s_MeasurementsPerBucket;
+        //! The mean magnitude of the bucket values.
+        double s_MeanMagnitude;
         //! The null hypothesis periodic components.
         CPeriodicityHypothesisTestsResult s_H0;
         //! The variance estimate of H0.

lib/maths/CPeriodicityHypothesisTests.cc

@@ -1154,8 +1154,16 @@ CPeriodicityHypothesisTests::best(const TNestedHypothesesVec& hypotheses) const
     THypothesisSummaryVec summaries;
     summaries.reserve(hypotheses.size());
 
+    double meanMagnitude{CBasicStatistics::mean(std::accumulate(
+        m_BucketValues.begin(), m_BucketValues.end(), TMeanAccumulator{},
+        [](TMeanAccumulator partial, const TFloatMeanAccumulator& value) {
+            partial.add(std::fabs(CBasicStatistics::mean(value)),
+                        CBasicStatistics::count(value));
+            return partial;
+        }))};
+
     for (const auto& hypothesis : hypotheses) {
-        STestStats stats;
+        STestStats stats{meanMagnitude};
         stats.s_TrendSegments = static_cast<double>(hypothesis.trendSegments());
         CPeriodicityHypothesisTestsResult resultForHypothesis{hypothesis.test(stats)};
         if (stats.s_NonEmptyBuckets > stats.s_DF0) {
@@ -2195,13 +2203,14 @@ bool CPeriodicityHypothesisTests::testAmplitude(const TTimeTimePr2Vec& window,
 const double CPeriodicityHypothesisTests::ACCURATE_TEST_POPULATED_FRACTION{0.9};
 const double CPeriodicityHypothesisTests::MINIMUM_COEFFICIENT_OF_VARIATION{1e-4};
 
-CPeriodicityHypothesisTests::STestStats::STestStats()
+CPeriodicityHypothesisTests::STestStats::STestStats(double meanMagnitude)
     : s_TrendSegments(1.0), s_HasPeriod(false), s_HasPartition(false),
       s_VarianceThreshold(COMPONENT_SIGNIFICANT_VARIANCE_REDUCTION[E_HighThreshold]),
       s_AmplitudeThreshold(SEASONAL_SIGNIFICANT_AMPLITUDE[E_HighThreshold]),
       s_AutocorrelationThreshold(SEASONAL_SIGNIFICANT_AUTOCORRELATION[E_HighThreshold]),
       s_Range(0.0), s_NonEmptyBuckets(0.0), s_MeasurementsPerBucket(0.0),
-      s_V0(0.0), s_R0(0.0), s_DF0(0.0), s_StartOfPartition(0) {
+      s_MeanMagnitude(meanMagnitude), s_V0(0.0), s_R0(0.0), s_DF0(0.0),
+      s_StartOfPartition(0) {
 }
 
 void CPeriodicityHypothesisTests::STestStats::setThresholds(double vt, double at, double Rt) {
@@ -2211,16 +2220,7 @@ void CPeriodicityHypothesisTests::STestStats::setThresholds(double vt, double at
 }
 
 bool CPeriodicityHypothesisTests::STestStats::nullHypothesisGoodEnough() const {
-    TMeanAccumulator mean;
-    for (const auto& t : s_T0) {
-        mean += std::accumulate(t.begin(), t.end(), TMeanAccumulator(),
-                                [](TMeanAccumulator m, double x) {
-                                    m.add(std::fabs(x));
-                                    return m;
-                                });
-    }
-    return std::sqrt(s_V0) <=
-           MINIMUM_COEFFICIENT_OF_VARIATION * CBasicStatistics::mean(mean);
+    return std::sqrt(s_V0) <= MINIMUM_COEFFICIENT_OF_VARIATION * s_MeanMagnitude;
 }
 
 CPeriodicityHypothesisTests::CNestedHypotheses::CNestedHypotheses(TTestFunc test)

lib/maths/CTimeSeriesDecompositionDetail.cc

@@ -1215,14 +1215,14 @@ void CTimeSeriesDecompositionDetail::CComponents::handle(const SAddValue& messag
         for (std::size_t i = 1u; i <= m; ++i) {
             CSeasonalComponent* component{seasonalComponents[i - 1]};
             CComponentErrors* error_{seasonalErrors[i - 1]};
-            double varianceIncrease{variances[i] / variance / expectedVarianceIncrease};
+            double varianceIncrease{variance == 0.0 ? 1.0 : variances[i] / variance / expectedVarianceIncrease};
             component->add(time, values[i], weight);
             error_->add(referenceError, error, predictions[i - 1], varianceIncrease, weight);
         }
         for (std::size_t i = m + 1; i <= m + n; ++i) {
             CCalendarComponent* component{calendarComponents[i - m - 1]};
             CComponentErrors* error_{calendarErrors[i - m - 1]};
-            double varianceIncrease{variances[i] / variance / expectedVarianceIncrease};
+            double varianceIncrease{variance == 0.0 ? 1.0 : variances[i] / variance / expectedVarianceIncrease};
             component->add(time, values[i], weight);
             error_->add(referenceError, error, predictions[i - 1], varianceIncrease, weight);
         }