Repository: alan-turing-institute/sktime
Pull request #3338 (2022-08-25 07:42:08)

[ENH] STRAY anomaly detection
#### Reference Issues/PRs
Issue #3337

#### What does this implement/fix? Explain your changes.
Implements the novel anomaly detection algorithm STRAY (Search TRace AnomalY), converted from R code.

An initial implementation, assuming this will be refactored once we have more of an idea about annotators. This is an interesting example to reference because a lot of logic in the function `_find_HDoutliers` deals with missing data (i.e. removes all missing data, then re-indexes back to the original data). Using pandas would make this indexing conversion easier, I think?
In any case, it highlights an open question: how should we handle this in general.

Similar to Hidalgo, all the logic is done in `fit`, so that `predict` just returns a specific reference to `self`.

Happy to hear any feedback!

@lmmentel @Lovkush-A @fkiraly @miraep8 @aiwalter 

KatieBuc: @lmmentel @Lovkush-A @fkiraly @miraep8 @aiwalter this is ready for review!

As suggested, for experiment, this anomaly detection algorithm inherits from `BaseTransformer`. The only thing I find clunky is that we can't write to self in `_transform`, hence the warning. A suggestion would be to have instead an `annotate` function that can write to self (instead of `predict` or `transform`). It would be nice to be able to bundle the logic in `annotate` and have public variables updated to self. Does that make any sense? What does everyone else think?



KatieBuc: Hey @lmmentel @Lovkush-A @fkiraly - something to flag about this normalization transformation pipeline.

1. In the original implementation of the algorithm, the NA's were dropped before normalization. Since one NA value deemed the whole row invalid. However, applying a column-wise transform for data with NA results in different normalization values i.e.
![image](https://user-images.githubusercontent.com/34525024/188807912-4338594a-03ba-4fa6-a440-5d129a8e843f.png)
where the top are the computed values within the algorithm, the bottom are when the transformation is applied _first_, then na removed.

A simple solution would be to add the tag "handles-missing-data: False" (does this mean if na's are passed in they will be automatically filtered out?) and what does this mean for the output? 

For example, 
with the input data
```
    X = np.array(
        [
            [-1.20706575, -0.57473996],
            [0.27742924, -0.54663186],
            [1.08444118, np.nan],
            [-2.3456977, -0.89003783],
            [0.42912469, -0.4771927],
            [0.50605589, -0.99838644],
        ]
    )
```

The anomaly boolean flags and scores are:
`y_expected = np.array([0, 0, 0, 1, 1, 1])`
`    y_scores_expected = np.array(
        [0.5233413, 0.5233413, np.nan, 0.7248368, 0.6035040, 0.8704687]
    )
`
according to the original index (with na)

2. Another thing to note is that
```
    >>> scaler = MinMaxScaler()
    >>> X = scaler.fit_transform(X)
```
throws an error when the input is 1D, hence `X = np.array(X).reshape(-1, 1)` is needed.

3. Also, in another example (DOBIN), the scaler changes type from `pd.DataFrame` to `np.Array`
```
    >>> _, X = load_uschange()
    >>> scaler = MinMaxScaler()
    >>> X = scaler.fit_transform(X)
    >>> model = DOBIN()
    >>> X_outlier = model.fit_transform(pd.DataFrame(X))
```
I think 2,3 make a case for interfacing, not sure?

and 1 is something for us to consider when thinking about anomaly outputs with missing data.
lmmentel: @KatieBuc thanks for elaborated comments. Here's some of my thoughts.

## Ad 1.

Does not seem like an algorithmic problem, rather a convention. You can decide what data your implementation expects and how to deal with it. To be honest I'm not seeing the issue here since it seems that for tests you are taking input and output from the R implementation (which was more of a pipeline and not a modular estimator) and trying to reproduce it here. Why not simply test the estimator on some other data? The rest of the issues are probably best left as docs comments/example on what are the options for preparing data, i.e.
- filter out `NaN`s -> normalize -> STRAY
- normalize -> filter out `NaN`s -> STRAY
- normalize -> STRAY

I'm not sure what your recommendation would be but the options above could be examples on how STRAY can be used not necessarily tests since their scope is larger that the estimator.

## Ad 2.

Yup, I think all sklearn X's are assumed to be 2D i.e. tabular data.

## Ad 3.

It seems your example already fixed the issue of type conversion right? This seems like another good example for the docs on "how to use STRAY with the existing sktime/sklearn preprocessors".

Patch:
diff --git a/sktime/annotation/stray.py b/sktime/annotation/stray.py
new file mode 100644
index 000000000..5d5d83d14
--- /dev/null
+++ b/sktime/annotation/stray.py
@@ -0,0 +1,266 @@
+# -*- coding: utf-8 -*-
+"""Tests for STRAY (Search TRace AnomalY) outlier estimator."""
+
+import warnings
+from typing import Dict
+
+import numpy as np
+import numpy.typing as npt
+from sklearn.neighbors import NearestNeighbors
+
+from sktime.transformations.base import BaseTransformer
+
+__author__ = ["KatieBuc"]
+__all__ = ["STRAY"]
+
+
+class STRAY(BaseTransformer):
+    """STRAY: robust anomaly detection in data streams with concept drift.
+
+    This is based on STRAY (Search TRace AnomalY) _[1], which is a modification
+    of HDoutliers _[2]. HDoutliers is a powerful algorithm for the detection of
+    anomalous observations in a dataset, which has (among other advantages) the
+    ability to detect clusters of outliers in multi-dimensional data without
+    requiring a model of the typical behavior of the system. However, it suffers
+    from some limitations that affect its accuracy. STRAY is an extension of
+    HDoutliers that uses extreme value theory for the anomolous threshold
+    calculation, to deal with data streams that exhibit non-stationary behavior.
+
+    Parameters
+    ----------
+    alpha : float, optional (default=0.01)
+        Threshold for determining the cutoff for outliers. Observations are
+        considered outliers if they fall in the (1 - alpha) tail of
+        the distribution of the nearest-neighbor distances between exemplars.
+    k : int, optional (default=10)
+        Number of neighbours considered.
+    knn_algorithm : str {"auto", "ball_tree", "kd_tree", "brute"}, optional
+        (default="brute")
+        Algorithm used to compute the nearest neighbors, from
+        sklearn.neighbors.NearestNeighbors
+    p : float, optional (default=0.5)
+        Proportion of possible candidates for outliers. This defines the starting point
+        for the bottom up searching algorithm.
+    size_threshold : int, optional (default=50)
+        Sample size to calculate an emperical threshold.
+    outlier_tail : str {"min", "max"}, optional (default="max")
+        Direction of the outlier tail.
+
+    Attributes
+    ----------
+    score_ : pd.Series
+        Outlier score of each data point in X.
+    y_ : pd.Series
+        Outlier boolean flag for each data point in X.
+
+    References
+    ----------
+    .. [1] Talagala, Priyanga Dilini, Rob J. Hyndman, and Kate Smith-Miles.
+    "Anomaly detection in high-dimensional data." Journal of Computational
+    and Graphical Statistics 30.2 (2021): 360-374.
+    .. [2] Wilkinson, Leland. "Visualizing big data outliers through
+    distributed aggregation." IEEE transactions on visualization and
+    computer graphics 24.1 (2017): 256-266.
+
+    Examples
+    --------
+    >>> from sktime.annotation.stray import STRAY
+    >>> from sktime.datasets import load_airline
+    >>> from sklearn.preprocessing import MinMaxScaler
+    >>> import numpy as np
+    >>> X = load_airline().to_frame().to_numpy()
+    >>> scaler = MinMaxScaler()
+    >>> X = scaler.fit_transform(X)
+    >>> model = STRAY(k=3)
+    >>> y = model.fit_transform(X)
+    >>> y[:5]
+    array([False, False, False, False, False])
+    """
+
+    _tags = {
+        "handles-missing-data": True,
+        "X_inner_mtype": "np.ndarray",
+        "fit_is_empty": False,
+        "skip-inverse-transform": True,
+    }
+
+    def __init__(
+        self,
+        alpha: float = 0.01,
+        k: int = 10,
+        knn_algorithm: str = "brute",
+        p: float = 0.5,
+        size_threshold: int = 50,
+        outlier_tail: str = "max",
+    ):
+        self.alpha = alpha
+        self.k = k
+        self.knn_algorithm = knn_algorithm
+        self.p = p
+        self.size_threshold = size_threshold
+        self.outlier_tail = outlier_tail
+        super(STRAY, self).__init__()
+
+    def _find_threshold(self, outlier_score: npt.ArrayLike, n: int) -> npt.ArrayLike:
+        """Find Outlier Threshold.
+
+        Parameters
+        ----------
+        outlier_score : np.ArrayLike
+            The outlier scores determined by k nearast neighbours distance
+        n : int
+            The number of rows remaining in X when NA's are removed.
+
+        Returns
+        -------
+        array of indices of the observations determined to be outliers.
+        """
+        if self.outlier_tail == "min":
+            outlier_score = -outlier_score
+
+        order = np.argsort(outlier_score)
+        gaps = np.append(0, np.diff(outlier_score[order]))
+        n4 = int(max(min(self.size_threshold, np.floor(n / 4)), 2))
+
+        J = np.array(range(2, n4 + 1))
+        start = int(max(np.floor(n * (1 - self.p)), 1))
+
+        ghat = [
+            0.0 if i < start else sum((J / (n4 - 1)) * gaps[i - J + 1])
+            for i in range(n)
+        ]
+
+        log_alpha = np.log(1 / self.alpha)
+        bound = np.Inf
+
+        for i in range(start, n):
+            if gaps[i] > log_alpha * ghat[i]:
+                bound = outlier_score[order][i - 1]
+                break
+
+        return np.where(outlier_score > bound)[0]
+
+    def _find_outliers_kNN(self, X: npt.ArrayLike, n: int) -> Dict:
+        """Find outliers using kNN distance with maximum gap.
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+        n : int
+            The number of rows remaining in X when NA's are removed.
+
+        Returns
+        -------
+        dict of index of outliers and the outlier scores
+        """
+        if len(X.shape) == 1:
+            nbrs = NearestNeighbors(n_neighbors=self.k + 1).fit(X.reshape(-1, 1))
+            distances, _ = nbrs.kneighbors(X.reshape(-1, 1))
+        else:
+            nbrs = NearestNeighbors(
+                n_neighbors=n if self.k >= n else self.k + 1,
+                algorithm=self.knn_algorithm,
+            ).fit(X)
+            distances, _ = nbrs.kneighbors(X)
+
+        if self.k == 1:
+            d = distances[:, 1]
+        else:
+            diff = np.apply_along_axis(np.diff, 1, distances)
+            d = distances[range(n), np.apply_along_axis(np.argmax, 1, diff) + 1]
+
+        out_index = self._find_threshold(d, n)
+        return {"idx_outliers": out_index, "out_scores": d}
+
+    def _find_outliers(self, X):
+        """Detect Anomalies in High Dimensional Data.
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+
+        Returns
+        -------
+        dict of anomalies and their corresponding scores
+        """
+        r = np.shape(X)[0]
+        idx_dropna = np.array([i for i in range(r) if not np.isnan(X[i]).any()])
+        X_dropna = X[
+            idx_dropna,
+        ]
+
+        n = np.shape(X_dropna)[0]
+        outliers = self._find_outliers_kNN(X_dropna, n)
+
+        # adjusted back to length r, for missing data
+        slice_ = [True if i in outliers["idx_outliers"] else False for i in range(n)]
+        idx_outliers = idx_dropna[slice_]  # index values from 1:r
+        outlier_bool = np.array([1 if i in idx_outliers else 0 for i in range(r)])
+
+        list_scores = outliers["out_scores"].tolist()
+        outlier_scores = np.array(
+            [list_scores.pop(0) if i in idx_dropna else np.nan for i in range(r)]
+        )
+
+        return {
+            "outlier_scores": outlier_scores,
+            "outlier_bool": outlier_bool,
+        }
+
+    def _fit(self, X: npt.ArrayLike, y: npt.ArrayLike = None):
+        """Find outliers using STRAY (Search TRace AnomalY).
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+        y : pd.Series, optional
+            Not used for this unsupervsed method.
+
+        Returns
+        -------
+        self :
+            Reference to self.
+        """
+        # remember X for transform
+        self._X = X
+
+        info_dict = self._find_outliers(X)
+        self.score_ = info_dict["outlier_scores"]
+        self.y_ = info_dict["outlier_bool"]
+
+        return self
+
+    def _transform(self, X: npt.ArrayLike, y: npt.ArrayLike = None) -> npt.ArrayLike:
+        """Return anomaly detection.
+
+        Parameters
+        ----------
+        X : np.ArrayLike
+            Data for anomaly detection (time series).
+
+        Returns
+        -------
+        y_ : np.ArrayLike
+            Anomaly detection, boolean.
+        """
+        # fit again if data is different to fit, but don't store anything
+        if not np.allclose(X, self._X, equal_nan=True):
+            new_obj = STRAY(
+                alpha=self.alpha,
+                k=self.k,
+                knn_algorithm=self.knn_algorithm,
+                p=self.p,
+                size_threshold=self.size_threshold,
+                outlier_tail=self.outlier_tail,
+            ).fit(X)
+            warnings.warn(
+                "Warning: Input data X differs from that given to fit(). "
+                "Refitting with new input data, not storing updated public class "
+                "attributes. For this, explicitly use fit(X) or fit_transform(X)."
+            )
+            return new_obj.y_.astype(bool)
+
+        return self.y_.astype(bool)


Test patch:
diff --git a/sktime/annotation/tests/test_stray.py b/sktime/annotation/tests/test_stray.py
new file mode 100644
index 000000000..bcb6127ce
--- /dev/null
+++ b/sktime/annotation/tests/test_stray.py
@@ -0,0 +1,289 @@
+# -*- coding: utf-8 -*-
+"""Tests for STRAY outlier estimator."""
+
+__author__ = ["KatieBuc"]
+
+import numpy as np
+from sklearn.preprocessing import MinMaxScaler, RobustScaler
+
+from sktime.annotation.stray import STRAY
+
+
+def test_default_1D():
+    """Test with default parameters and 1D input array."""
+    X = np.array(
+        [
+            -7.207066,
+            -5.722571,
+            -4.915559,
+            -8.345698,
+            -5.570875,
+            -5.493944,
+            -6.574740,
+            -6.546632,
+            -6.564452,
+            -6.890038,
+            -6.477193,
+            -6.998386,
+            -6.776254,
+            -5.935541,
+            -5.040506,
+            0.000000,
+            5.889715,
+            5.488990,
+            5.088805,
+            5.162828,
+            8.415835,
+            6.134088,
+            5.509314,
+            5.559452,
+            6.459589,
+            5.306280,
+            4.551795,
+            6.574756,
+            4.976344,
+            5.984862,
+            5.064051,
+            7.102298,
+        ]
+    )
+
+    y_expected = np.array(
+        [
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            True,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+        ]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X.reshape(-1, 1))
+    model = STRAY()
+    y_actual = model.fit_transform(X)
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_default_2D():
+    """Test with default parameters and 2D input array."""
+    X = np.array(
+        [
+            [-1.207, -0.776, -0.694],
+            [0.277, 0.064, -1.448],
+            [1.084, 0.959, 0.575],
+            [-2.346, -0.11, -1.024],
+            [0.429, -0.511, -0.015],
+            [0.506, -0.911, -0.936],
+            [-0.575, -0.837, 1.102],
+            [-0.547, 2.416, -0.476],
+            [-0.564, 0.134, -0.709],
+            [-0.89, -0.491, -0.501],
+            [-0.477, -0.441, -1.629],
+            [-0.998, 0.46, -1.168],
+            [10.0, 12.0, 10.0],
+            [3.0, 7.0, 10.0],
+        ]
+    )
+
+    y_expected = np.array(
+        [
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            False,
+            True,
+            True,
+        ]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY()
+    y_actual = model.fit_transform(X)
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_1D_score_with_na():
+    """Test score with 1D input array with missing values."""
+    X = np.array(
+        [
+            np.nan,
+            -5.72257076,
+            -4.91555882,
+            -8.3456977,
+            -5.57087531,
+            0.0,
+            6.50605589,
+            5.42526004,
+            5.45336814,
+            5.435548,
+            5.10996217,
+        ]
+    )
+
+    y_scores_expected = np.array(
+        [
+            np.nan,
+            0.17662069,
+            0.23095851,
+            0.17662069,
+            0.18683466,
+            0.33097498,
+            0.07087969,
+            0.02122967,
+            0.02312225,
+            0.02192238,
+            0.02122967,
+        ]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X.reshape(-1, 1))
+    model = STRAY(k=3)
+    fitted_model = model.fit(X)
+    y_scores_actual = fitted_model.score_
+    assert np.allclose(y_scores_actual, y_scores_expected, equal_nan=True)
+
+
+def test_1D_bool_with_na():
+    """Test anomaly detection with 1D input array with missing values."""
+    X = np.array(
+        [
+            np.nan,
+            -5.72257076,
+            -4.91555882,
+            -8.3456977,
+            -5.57087531,
+            0.0,
+            6.50605589,
+            5.42526004,
+            5.45336814,
+            5.435548,
+            5.10996217,
+        ]
+    )
+
+    y_expected = np.array([0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0])
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X.reshape(-1, 1))
+    model = STRAY(k=3)
+    fitted_model = model.fit(X)
+    y_actual = fitted_model.y_
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_2D_score_with_na():
+    """Test score with 2D input array with missing values."""
+    X = np.array(
+        [
+            [-1.20706575, -0.57473996],
+            [0.27742924, -0.54663186],
+            [1.08444118, np.nan],
+            [-2.3456977, -0.89003783],
+            [0.42912469, -0.4771927],
+            [0.50605589, -0.99838644],
+        ]
+    )
+
+    y_scores_expected = np.array(
+        [0.43612712, 0.43612712, np.nan, 0.69004259, 0.51240831, 0.85697804]
+    )
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY(k=2, size_threshold=4)
+    fitted_model = model.fit(X)
+    y_scores_actual = fitted_model.score_
+    assert np.allclose(y_scores_actual, y_scores_expected, equal_nan=True)
+
+
+def test_2D_bool_with_na():
+    """Test anomaly detection with 2D input array with missing values."""
+    X = np.array(
+        [
+            [-1.20706575, -0.57473996],
+            [0.27742924, -0.54663186],
+            [1.08444118, np.nan],
+            [-2.3456977, -0.89003783],
+            [0.42912469, -0.4771927],
+            [0.50605589, -0.99838644],
+        ]
+    )
+
+    y_expected = np.array([0, 0, 0, 1, 1, 1])
+
+    scaler = MinMaxScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY(k=2, size_threshold=4)
+    fitted_model = model.fit(X)
+    y_actual = fitted_model.y_
+    assert np.allclose(y_actual, y_expected)
+
+
+def test_2D_score_with_standardize():
+    """Test score with 2D input array and median/IQR normalization."""
+    X = np.array(
+        [
+            [-1.20706575, -0.57473996],
+            [0.27742924, -0.54663186],
+            [1.08444118, -0.5644520],
+            [-2.3456977, -0.89003783],
+            [0.42912469, -0.4771927],
+            [0.50605589, -0.99838644],
+        ]
+    )
+
+    y_scores_expected = np.array(
+        [
+            1.1274565,
+            0.6139288,
+            0.5982989,
+            1.4866554,
+            0.5982989,
+            1.7245212,
+        ]
+    )
+
+    scaler = RobustScaler()
+    X = scaler.fit_transform(X)
+    model = STRAY(k=2, size_threshold=4)
+    fitted_model = model.fit(X)
+    y_scores_actual = fitted_model.score_
+    assert np.allclose(y_scores_actual, y_scores_expected)