[ML] Improve robustness to outliers after detecting changes in time series

docs/CHANGELOG.asciidoc

@@ -36,6 +36,8 @@
   and {ml-pull}2238[#2238].)
 * Upgrade zlib to version 1.2.12 on Windows. (See {ml-pull}2253[#2253].)
 * Address root cause for actuals equals typical equals zero anomalies. (See {ml-pull}2270[#2270].)
+* Better handling of outliers in update immediately after detecting changes in time
+  series. (See {ml-pull}2280[#2280].)
 
 === Bug Fixes
 

include/maths/time_series/CTimeSeriesDecomposition.h

@@ -196,7 +196,7 @@ class MATHS_TIME_SERIES_EXPORT EMPTY_BASE_OPT CTimeSeriesDecomposition
     double countWeight(core_t::TTime time) const override;
 
     //! Get the derate to apply to the Winsorisation weight at \p time.
-    double winsorisationDerate(core_t::TTime time) const override;
+    double winsorisationDerate(core_t::TTime time, double error) const override;
 
     //! Get the prediction residuals in a recent time window.
     //!

include/maths/time_series/CTimeSeriesDecompositionDetail.h

@@ -249,8 +249,9 @@ class MATHS_TIME_SERIES_EXPORT CTimeSeriesDecompositionDetail
         //! Get the count weight to apply to samples.
         double countWeight(core_t::TTime time) const;
 
-        //! Get the derate to apply to the Winsorisation weight.
-        double winsorisationDerate(core_t::TTime time) const;
+        //! Get the derate to apply to the Winsorization weight for a prediction error
+        //! of size \p error.
+        double winsorisationDerate(core_t::TTime time, double error) const;
 
         //! Age the test to account for the interval \p end - \p start elapsed time.
         void propagateForwards(core_t::TTime start, core_t::TTime end);
@@ -337,13 +338,13 @@ class MATHS_TIME_SERIES_EXPORT CTimeSeriesDecompositionDetail
         TMeanVarAccumulator m_ResidualMoments;
 
         //! The proportion of recent values with significantly prediction error.
-        double m_LargeErrorFraction = 0.0;
+        double m_LargeErrorFraction{0.0};
 
         //! The total adjustment applied to the count weight.
-        double m_TotalCountWeightAdjustment = 0.0;
+        double m_TotalCountWeightAdjustment{0.0};
 
         //! The minimum permitted total adjustment applied to the count weight.
-        double m_MinimumTotalCountWeightAdjustment = 0.0;
+        double m_MinimumTotalCountWeightAdjustment{0.0};
 
         //! The last test time.
         core_t::TTime m_LastTestTime;
@@ -357,6 +358,9 @@ class MATHS_TIME_SERIES_EXPORT CTimeSeriesDecompositionDetail
         //! The last change which was made, if it hasn't been committed, in a form
         //! which can be undone.
         TChangePointUPtr m_UndoableLastChange;
+
+        //! The derate to apply to the Winsorization immediately after the last change point.
+        CWinsorizationDerate m_LastChangeWinsorizationDerate;
     };
 
     //! \brief Scans through increasingly low frequencies looking for significant

include/maths/time_series/CTimeSeriesDecompositionInterface.h

@@ -181,7 +181,7 @@ class MATHS_TIME_SERIES_EXPORT CTimeSeriesDecompositionInterface
     virtual double countWeight(core_t::TTime time) const = 0;
 
     //! Get the derate to apply to the Winsorisation weight at \p time.
-    virtual double winsorisationDerate(core_t::TTime time) const = 0;
+    virtual double winsorisationDerate(core_t::TTime time, double derate) const = 0;
 
     //! Get the prediction residuals in a recent time window.
     //!

include/maths/time_series/CTimeSeriesDecompositionStub.h

@@ -93,7 +93,7 @@ class MATHS_TIME_SERIES_EXPORT CTimeSeriesDecompositionStub
     double countWeight(core_t::TTime time) const override;
 
     //! Returns 0.0.
-    double winsorisationDerate(core_t::TTime time) const override;
+    double winsorisationDerate(core_t::TTime time, double error) const override;
 
     //! Returns an empty vector.
     TFloatMeanAccumulatorVec residuals(bool isNonNegative) const override;

include/maths/time_series/CTimeSeriesTestForChange.h

@@ -21,6 +21,7 @@
 #include <maths/time_series/ImportExport.h>
 
 #include <functional>
+#include <limits>
 #include <memory>
 #include <tuple>
 #include <vector>
@@ -37,6 +38,20 @@ class CSeasonalComponent;
 class CTimeSeriesDecomposition;
 class CTrendComponent;
 
+//! \brief Determines the Winsorization derate which should apply after a change.
+class MATHS_TIME_SERIES_EXPORT CWinsorizationDerate {
+public:
+    CWinsorizationDerate() = default;
+    explicit CWinsorizationDerate(double magnitude);
+    double value(double error) const;
+    bool acceptRestoreTraverser(core::CStateRestoreTraverser& traverser);
+    void acceptPersistInserter(core::CStatePersistInserter& inserter) const;
+    std::uint64_t checksum(std::uint64_t seed = 0) const;
+
+private:
+    double m_Magnitude{std::numeric_limits<double>::max()};
+};
+
 //! \brief Represents a sudden change to a time series model.
 class MATHS_TIME_SERIES_EXPORT CChangePoint {
 public:
@@ -51,6 +66,9 @@ class MATHS_TIME_SERIES_EXPORT CChangePoint {
     virtual ~CChangePoint();
 
     virtual TChangePointUPtr undoable() const = 0;
+    virtual CWinsorizationDerate winsorizationDerate(core_t::TTime startTime,
+                                                     core_t::TTime endTime,
+                                                     const TPredictor& predictor) const = 0;
     virtual bool largeEnough(double threshold) const = 0;
     virtual bool longEnough(core_t::TTime time, core_t::TTime minimumDuration) const = 0;
     virtual bool apply(CTimeSeriesDecomposition&) const { return false; }
@@ -70,7 +88,7 @@ class MATHS_TIME_SERIES_EXPORT CChangePoint {
     TFloatMeanAccumulatorVec& residuals() { return m_Residuals; }
     const TFloatMeanAccumulatorVec& residuals() const { return m_Residuals; }
 
-private:
+protected:
     using TMeanAccumulator = common::CBasicStatistics::SSampleMean<double>::TAccumulator;
 
 private:
@@ -79,15 +97,15 @@ class MATHS_TIME_SERIES_EXPORT CChangePoint {
     }
 
 private:
-    core_t::TTime m_Time = 0;
-    double m_SignificantPValue = 0.0;
+    core_t::TTime m_Time{0};
+    double m_SignificantPValue{0.0};
     TFloatMeanAccumulatorVec m_Residuals;
     TMeanAccumulator m_Mse;
     TMeanAccumulator m_UndoneMse;
 };
 
 //! \brief Represents a level shift of a time series.
-class MATHS_TIME_SERIES_EXPORT CLevelShift : public CChangePoint {
+class MATHS_TIME_SERIES_EXPORT CLevelShift final : public CChangePoint {
 public:
     using TDoubleVec = std::vector<double>;
     using TSizeVec = std::vector<std::size_t>;
@@ -107,6 +125,10 @@ class MATHS_TIME_SERIES_EXPORT CLevelShift : public CChangePoint {
                 double significantPValue);
 
     TChangePointUPtr undoable() const override;
+    CWinsorizationDerate
+    winsorizationDerate(core_t::TTime, core_t::TTime, const TPredictor&) const override {
+        return CWinsorizationDerate{m_Shift};
+    }
     bool largeEnough(double threshold) const override;
     bool longEnough(core_t::TTime time, core_t::TTime minimumDuration) const override;
     bool apply(CTrendComponent& component) const override;
@@ -116,16 +138,16 @@ class MATHS_TIME_SERIES_EXPORT CLevelShift : public CChangePoint {
     std::uint64_t checksum(std::uint64_t seed = 0) const override;
 
 private:
-    double m_Shift = 0.0;
-    core_t::TTime m_ValuesStartTime = 0;
-    core_t::TTime m_BucketLength = 0;
+    double m_Shift{0.0};
+    core_t::TTime m_ValuesStartTime{0};
+    core_t::TTime m_BucketLength{0};
     TFloatMeanAccumulatorVec m_Values;
     TSizeVec m_Segments;
     TDoubleVec m_Shifts;
 };
 
 //! \brief Represents a linear scale of a time series.
-class MATHS_TIME_SERIES_EXPORT CScale : public CChangePoint {
+class MATHS_TIME_SERIES_EXPORT CScale final : public CChangePoint {
 public:
     static const std::string TYPE;
 
@@ -138,6 +160,10 @@ class MATHS_TIME_SERIES_EXPORT CScale : public CChangePoint {
            double significantPValue);
 
     TChangePointUPtr undoable() const override;
+    CWinsorizationDerate
+    winsorizationDerate(core_t::TTime, core_t::TTime, const TPredictor&) const override {
+        return CWinsorizationDerate{m_Magnitude};
+    }
     bool largeEnough(double threshold) const override;
     bool longEnough(core_t::TTime time, core_t::TTime minimumDuration) const override;
     bool apply(CTrendComponent& component) const override;
@@ -149,13 +175,13 @@ class MATHS_TIME_SERIES_EXPORT CScale : public CChangePoint {
     std::uint64_t checksum(std::uint64_t seed = 0) const override;
 
 private:
-    double m_Scale = 1.0;
-    double m_Magnitude = 0.0;
-    double m_MinimumDurationScale = 1.0;
+    double m_Scale{1.0};
+    double m_Magnitude{0.0};
+    double m_MinimumDurationScale{1.0};
 };
 
 //! \brief Represents a time shift of a time series.
-class MATHS_TIME_SERIES_EXPORT CTimeShift : public CChangePoint {
+class MATHS_TIME_SERIES_EXPORT CTimeShift final : public CChangePoint {
 public:
     static const std::string TYPE;
 
@@ -168,6 +194,9 @@ class MATHS_TIME_SERIES_EXPORT CTimeShift : public CChangePoint {
     CTimeShift(core_t::TTime time, core_t::TTime shift, double significantPValue);
 
     TChangePointUPtr undoable() const override;
+    CWinsorizationDerate winsorizationDerate(core_t::TTime startTime,
+                                             core_t::TTime endTime,
+                                             const TPredictor& predictor) const override;
     bool largeEnough(double) const override { return m_Shift != 0; }
     bool longEnough(core_t::TTime time, core_t::TTime minimumDuration) const override;
     bool apply(CTimeSeriesDecomposition& decomposition) const override;
@@ -180,7 +209,7 @@ class MATHS_TIME_SERIES_EXPORT CTimeShift : public CChangePoint {
     double undonePredict(const TPredictor& predictor, core_t::TTime time) const override;
 
 private:
-    core_t::TTime m_Shift = 0;
+    core_t::TTime m_Shift{0};
 };
 
 //! \brief Manages persist and restore of an undoable change point.
@@ -292,9 +321,9 @@ class MATHS_TIME_SERIES_EXPORT CTimeSeriesTestForChange {
             : s_ResidualVariance{residualVariance}, s_TruncatedResidualVariance{truncatedResidualVariance},
               s_NumberParameters{numberParameters}, s_ChangePoint{std::move(changePoint)} {}
 
-        double s_ResidualVariance = 0.0;
-        double s_TruncatedResidualVariance = 0.0;
-        double s_NumberParameters = 0.0;
+        double s_ResidualVariance{0.0};
+        double s_TruncatedResidualVariance{0.0};
+        double s_NumberParameters{0.0};
         TChangePointUPtr s_ChangePoint;
     };
 

lib/api/CAnomalyJob.cc

@@ -98,7 +98,7 @@ const std::string INTERIM_BUCKET_CORRECTOR_TAG("k");
 //! boundary.  Model snapshots generated in 8.x will not be loadable by 7.x, and
 //! when 7.x is end-of-life we'll be able to remove all the 7.x state backwards
 //! compatibility code.)
-const std::string MODEL_SNAPSHOT_MIN_VERSION("8.0.0");
+const std::string MODEL_SNAPSHOT_MIN_VERSION("8.3.0");
 
 //! Persist state as JSON with meaningful tag names.
 class CReadableJsonStatePersistInserter : public core::CJsonStatePersistInserter {

lib/maths/time_series/CTimeSeriesDecomposition.cc

@@ -336,7 +336,9 @@ CTimeSeriesDecomposition::value(core_t::TTime time, double confidence, bool isNo
 CTimeSeriesDecomposition::TFilteredPredictor
 CTimeSeriesDecomposition::predictor(int components) const {
 
-    auto trend_ = m_Components.trend().predictor();
+    auto trend_ = (((components & E_TrendForced) != 0) || ((components & E_Trend) != 0))
+                      ? m_Components.trend().predictor()
+                      : [](core_t::TTime) { return 0.0; };
 
     return [ components, trend = std::move(trend_),
              this ](core_t::TTime time, const TBoolVec& removedSeasonalMask) {
@@ -547,8 +549,8 @@ double CTimeSeriesDecomposition::countWeight(core_t::TTime time) const {
     return m_ChangePointTest.countWeight(time);
 }
 
-double CTimeSeriesDecomposition::winsorisationDerate(core_t::TTime time) const {
-    return m_ChangePointTest.winsorisationDerate(time);
+double CTimeSeriesDecomposition::winsorisationDerate(core_t::TTime time, double error) const {
+    return m_ChangePointTest.winsorisationDerate(time, error);
 }
 
 CTimeSeriesDecomposition::TFloatMeanAccumulatorVec

lib/maths/time_series/CTimeSeriesDecompositionDetail.cc

@@ -317,6 +317,8 @@ const core::TPersistenceTag LAST_CHANGE_POINT_TIME_7_11_TAG{"i", "last_change_po
 const core::TPersistenceTag LAST_CANDIDATE_CHANGE_POINT_TIME_7_11_TAG{
     "j", "last_candidate_change_point_time"};
 const core::TPersistenceTag LAST_CHANGE_POINT_7_11_TAG{"k", "last_change_point"};
+// Version 8.3
+const core::TPersistenceTag WINSORIZATION_DERATE_8_3_TAG{"l", "winsorization_derate"};
 
 // Seasonality Test Tags
 // Version 7.9
@@ -505,7 +507,8 @@ CTimeSeriesDecompositionDetail::CChangePointTest::CChangePointTest(const CChange
       m_TotalCountWeightAdjustment{other.m_TotalCountWeightAdjustment},
       m_MinimumTotalCountWeightAdjustment{other.m_MinimumTotalCountWeightAdjustment},
       m_LastTestTime{other.m_LastTestTime}, m_LastChangePointTime{other.m_LastChangePointTime},
-      m_LastCandidateChangePointTime{other.m_LastCandidateChangePointTime} {
+      m_LastCandidateChangePointTime{other.m_LastCandidateChangePointTime},
+      m_LastChangeWinsorizationDerate{other.m_LastChangeWinsorizationDerate} {
 
     if (isForForecast) {
         this->apply(CD_DISABLE);
@@ -539,6 +542,10 @@ bool CTimeSeriesDecompositionDetail::CChangePointTest::acceptRestoreTraverser(
         ](auto& traverser_) {
             return serializer(m_UndoableLastChange, traverser_);
         }))
+        RESTORE(WINSORIZATION_DERATE_8_3_TAG,
+                traverser.traverseSubLevel([this](core::CStateRestoreTraverser& traverser_) {
+                    return m_LastChangeWinsorizationDerate.acceptRestoreTraverser(traverser_);
+                }))
     } while (traverser.next());
     return true;
 }
@@ -567,6 +574,9 @@ void CTimeSeriesDecompositionDetail::CChangePointTest::acceptPersistInserter(
             this, serializer = CUndoableChangePointStateSerializer{}
         ](auto& inserter_) { serializer(*m_UndoableLastChange, inserter_); });
     }
+    inserter.insertLevel(WINSORIZATION_DERATE_8_3_TAG, [this](auto& inserter_) {
+        return m_LastChangeWinsorizationDerate.acceptPersistInserter(inserter_);
+    });
 }
 
 void CTimeSeriesDecompositionDetail::CChangePointTest::swap(CChangePointTest& other) {
@@ -583,6 +593,7 @@ void CTimeSeriesDecompositionDetail::CChangePointTest::swap(CChangePointTest& ot
     std::swap(m_LastChangePointTime, other.m_LastChangePointTime);
     std::swap(m_LastCandidateChangePointTime, other.m_LastCandidateChangePointTime);
     std::swap(m_UndoableLastChange, other.m_UndoableLastChange);
+    std::swap(m_LastChangeWinsorizationDerate, other.m_LastChangeWinsorizationDerate);
 }
 
 void CTimeSeriesDecompositionDetail::CChangePointTest::handle(const SAddValue& message) {
@@ -622,7 +633,7 @@ void CTimeSeriesDecompositionDetail::CChangePointTest::handle(const SAddValue& m
 }
 
 void CTimeSeriesDecompositionDetail::CChangePointTest::handle(const SDetectedSeasonal& message) {
-    if (m_Window.size() > 0) {
+    if (m_Window.empty() == false) {
         m_Window.assign(m_Window.size(), TFloatMeanAccumulator{});
     }
     m_ResidualMoments = TMeanVarAccumulator{};
@@ -645,10 +656,12 @@ double CTimeSeriesDecompositionDetail::CChangePointTest::countWeight(core_t::TTi
     return 1.0 + std::min(1.0, -m_TotalCountWeightAdjustment);
 }
 
-double CTimeSeriesDecompositionDetail::CChangePointTest::winsorisationDerate(core_t::TTime time) const {
+double CTimeSeriesDecompositionDetail::CChangePointTest::winsorisationDerate(core_t::TTime time,
+                                                                             double error) const {
     return std::max(1.0 - static_cast<double>(time - m_LastChangePointTime) /
                               static_cast<double>(3 * DAY),
-                    0.0);
+                    0.0) *
+           m_LastChangeWinsorizationDerate.value(error);
 }
 
 void CTimeSeriesDecompositionDetail::CChangePointTest::propagateForwards(core_t::TTime start,
@@ -671,7 +684,8 @@ std::uint64_t CTimeSeriesDecompositionDetail::CChangePointTest::checksum(std::ui
     seed = common::CChecksum::calculate(seed, m_LastTestTime);
     seed = common::CChecksum::calculate(seed, m_LastChangePointTime);
     seed = common::CChecksum::calculate(seed, m_LastCandidateChangePointTime);
-    return common::CChecksum::calculate(seed, m_UndoableLastChange);
+    seed = common::CChecksum::calculate(seed, m_UndoableLastChange);
+    return common::CChecksum::calculate(seed, m_LastChangeWinsorizationDerate);
 }
 
 void CTimeSeriesDecompositionDetail::CChangePointTest::debugMemoryUsage(
@@ -752,7 +766,7 @@ void CTimeSeriesDecompositionDetail::CChangePointTest::testForChange(const SAddV
     core_t::TTime time{message.s_Time};
     core_t::TTime lastTime{message.s_LastTime};
     core_t::TTime timeShift{message.s_TimeShift};
-    bool seasonal{message.s_Decomposition->seasonalComponents().size() > 0};
+    bool seasonal{message.s_Decomposition->seasonalComponents().empty() == false};
     const auto& makePredictor = message.s_MakePredictor;
     CTimeSeriesDecomposition& decomposition{*message.s_Decomposition};
 
@@ -799,6 +813,8 @@ void CTimeSeriesDecompositionDetail::CChangePointTest::testForChange(const SAddV
         m_LastCandidateChangePointTime = std::min(
             m_LastCandidateChangePointTime, time - this->maximumIntervalToDetectChange());
         m_UndoableLastChange = change->undoable();
+        m_LastChangeWinsorizationDerate =
+            change->winsorizationDerate(bucketsStartTime, time, predictor);
         this->mediator()->forward(SDetectedChangePoint{time, lastTime, std::move(change)});
     } else if (change != nullptr) {
         m_LastCandidateChangePointTime = change->time();

lib/maths/time_series/CTimeSeriesDecompositionStub.cc

@@ -100,7 +100,8 @@ double CTimeSeriesDecompositionStub::countWeight(core_t::TTime /*time*/) const {
     return 1.0;
 }
 
-double CTimeSeriesDecompositionStub::winsorisationDerate(core_t::TTime /*time*/) const {
+double CTimeSeriesDecompositionStub::winsorisationDerate(core_t::TTime /*time*/,
+                                                         double /*error*/) const {
     return 0.0;
 }
 

lib/maths/time_series/CTimeSeriesModel.cc

@@ -1177,7 +1177,8 @@ void CUnivariateTimeSeriesModel::countWeights(core_t::TTime time,
     maths_t::setSeasonalVarianceScale(seasonalWeight[0], weights);
     double winsorisationWeight{winsorisation::weight(
         *m_ResidualModel, weights,
-        std::max(winsorisationDerate, m_TrendModel->winsorisationDerate(time)), sample)};
+        std::max(winsorisationDerate, m_TrendModel->winsorisationDerate(time, sample)),
+        sample)};
 
     double changeWeight{m_TrendModel->countWeight(time)};
     trendCountWeight /= countVarianceScale;
@@ -2572,7 +2573,8 @@ void CMultivariateTimeSeriesModel::countWeights(core_t::TTime time,
         maths_t::setSeasonalVarianceScale(seasonalWeight[d], weights);
         double winsorisationWeight{winsorisation::weight(
             *conditional(*m_ResidualModel, d, sample), weights,
-            std::max(winsorisationDerate, m_TrendModel[d]->winsorisationDerate(time)),
+            std::max(winsorisationDerate,
+                     m_TrendModel[d]->winsorisationDerate(time, sample[d])),
             sample[d])};
         residualCountWeights[d] *= changeWeight;
         trendWinsorisationWeight[d] = winsorisationWeight * changeWeight;