7.3 Anomaly detection

An anomaly is a case unlike the others. Without labels, “unlike” must be defined: far from the centre, in a sparse region, easy to isolate. Each definition finds different anomalies, and anomalies themselves can hide the very spread that should have revealed them.

Code: data and style used in this chapter 
from collections.abc import Callable

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.covariance import EmpiricalCovariance, MinCovDet
from sklearn.ensemble import IsolationForest
from sklearn.metrics import precision_recall_curve
from sklearn.neighbors import LocalOutlierFactor

from statly import style

style.use()
rng = np.random.default_rng(1)

wine = datasets.load_wine(as_frame=True)  # 178 Italian wines, 13 chemical measurements


def score_map(
    ax: Axes,
    points: np.ndarray,
    score_fn: Callable[[np.ndarray], np.ndarray],
    title: str,
    flagged: np.ndarray,
    flip: bool = False,
) -> None:
    """Background shading by anomaly score (darker: more anomalous) with the points on top."""
    xs = np.linspace(points[:, 0].min() - 1, points[:, 0].max() + 1, 150)
    ys = np.linspace(points[:, 1].min() - 1, points[:, 1].max() + 1, 150)
    grid = np.array(np.meshgrid(xs, ys)).reshape(2, -1).T
    scores = (
        -score_fn(grid) if flip else score_fn(grid)
    )  # some scorers give higher = more normal
    ax.contourf(xs, ys, scores.reshape(150, 150), levels=12, cmap="Reds", alpha=0.6)
    ax.scatter(*points[~flagged].T, s=6, color=style.INK, alpha=0.6)
    ax.scatter(*points[flagged].T, s=40, facecolors="none", edgecolors=style.BLUE, lw=1.5)
    ax.set(title=title, xticks=[], yticks=[])

7.3.1 One variable: how far is far?

The simplest rule flags values more than three standard deviations from the mean. Its weakness is that the anomalies themselves inflate the mean and the standard deviation they are measured against (masking): a group of extreme values can hide itself. Robust versions measure distance from the median, in units of the median absolute deviation (MAD), which a few extreme values barely move:

zi=xi−x‾sversuszirobust=xi−median1.4826⋅MAD z_i = \frac{x_i - \bar x}{s} \quad\text{versus}\quad z_i^{\text{robust}} = \frac{x_i - \text{median}}{1.4826 \cdot \text{MAD}}

(The factor 1.4826 makes the MAD match the standard deviation for normal data.)

Code 
values = np.r_[rng.normal(0, 1, 200), rng.normal(4.6, 0.3, 35)]
mean, sd = values.mean(), values.std()
median = np.median(values)
mad = 1.4826 * np.median(np.abs(values - median))

fig, axes = style.subplots(2, height=3, width=10)
for ax, (centre, spread, k, color, title) in zip(
    axes,
    (
        (mean, sd, 3, style.GREY, "mean ± 3 sd"),
        (median, mad, 3.5, style.BLUE, "median ± 3.5 robust sd"),
    ),
):
    flagged = np.abs(values - centre) > k * spread
    style.hist(ax, values, bins=50, color=style.GREY)
    ax.axvspan(centre - k * spread, centre + k * spread, color=color, alpha=0.15)
    ax.scatter(
        values[flagged],
        np.zeros(flagged.sum()),
        color=style.CORAL,
        s=25,
        zorder=3,
        clip_on=False,
    )
    style.note(ax, f"{flagged.sum()} flagged", x=0.03, y=0.95, color=style.INK)
    ax.set(title=title, xlabel="measurement", ylabel="readings")
fig.tight_layout()
plt.show()
Figure 7.3.1: 200 normal measurements with 35 faulty readings clustered around 4.6. Left: the faulty readings pull the mean right and inflate the standard deviation, so the ±3 sd band (grey) stretches to include them, and none is flagged. Right: the median and MAD hardly move, and the ±3.5 robust band (blue) flags nearly all of them.

7.3.2 Several variables: unusual combinations

A case can be ordinary on every measurement and unusual in their combination: a tall person with a child’s weight. The Mahalanobis distance measures distance from the centre in units that follow the shape of the data cloud (Chapter 3.9):

di=(xi−x‾)⊤S−1(xi−x‾) d_i = \sqrt{(x_i - \bar x)^\top S^{-1} (x_i - \bar x)}

With classical estimates of the centre and covariance, the same masking happens in many dimensions. The minimum covariance determinant estimator fits the centre and covariance to the most tightly packed half of the data, so a cluster of anomalies cannot stretch it.

Code 
normal = rng.multivariate_normal([0, 0], [[1, 0.9], [0.9, 1]], 300)
odd = rng.multivariate_normal([1.2, -1.2], [[0.05, 0], [0, 0.05]], 40)
points = np.r_[normal, odd]
is_odd = np.r_[np.zeros(300, bool), np.ones(40, bool)]

fig, axes = style.subplots(2, height=3.6, width=10)
for ax, estimator, title in (
    (axes[0], EmpiricalCovariance(), "classical"),
    (axes[1], MinCovDet(random_state=0), "robust"),
):
    fit = estimator.fit(points)
    distance = fit.mahalanobis(points)
    flagged = np.zeros(len(points), bool)
    flagged[np.argsort(distance)[-40:]] = True
    score_map(ax, points, fit.mahalanobis, title, flagged)
    ax.scatter(*odd.T, s=10, color=style.CORAL, zorder=3)
    style.note(
        ax,
        f"anomalies among top 40: {(flagged & is_odd).sum()}",
        x=0.03,
        y=0.95,
        color=style.INK,
    )
fig.tight_layout()
plt.show()
Figure 7.3.2: Two strongly correlated measurements for 300 normal cases and 40 anomalies (red) that are unremarkable on each measurement alone but break the correlation. Left: contours of the classical Mahalanobis distance, stretched by the anomalies, so that some normal cases at the ends of the cloud rank above them. Right: the robust version (minimum covariance determinant) follows the normal cases closely, and nearly all anomalies lie far outside its contours. Blue circles: the 40 most distant points by each method.

7.3.3 Isolation forests

An isolation forest turns the question around: instead of describing normal data, it measures how easy a case is to separate from the rest. It builds many random trees, each splitting on a random measurement at a random value, until every case sits alone in a leaf. Anomalies, being few and different, are isolated after few splits; normal cases, deep in the crowd, need many. The anomaly score is based on the average depth at which a case is isolated.

Code 
def isolate(
    points: np.ndarray, target: int, local: np.random.Generator, box: list[float]
) -> list[tuple]:
    """Random cuts until `target` is alone; returns the cuts that contained it."""
    cuts, members = [], np.arange(len(points))
    lo, hi = np.array(box[:2], float), np.array(box[2:], float)
    while len(members) > 1:
        axis = local.integers(2)
        value = local.uniform(points[members, axis].min(), points[members, axis].max())
        cuts.append((axis, value, lo.copy(), hi.copy()))
        side = points[members, axis] < value
        keep = side if points[target, axis] < value else ~side
        if points[target, axis] < value:
            hi[axis] = value
        else:
            lo[axis] = value
        members = members[keep]
    return cuts


def draw_cuts(ax: Axes, cuts: list[tuple], color: str) -> None:
    """Draw each cut as a line segment across the region it split."""
    for axis, value, lo, hi in cuts:
        if axis == 0:
            ax.plot([value, value], [lo[1], hi[1]], color=color, lw=1)
        else:
            ax.plot([lo[0], hi[0]], [value, value], color=color, lw=1)


box = [points[:, 0].min(), points[:, 1].min(), points[:, 0].max(), points[:, 1].max()]
edge_case = int(np.argmax(np.abs(points[:, 0]) + np.abs(points[:, 1])))
middle_case = int(np.argmin(np.linalg.norm(points, axis=1)))
forest = IsolationForest(n_estimators=200, random_state=0).fit(points)
iso_score = -forest.score_samples(points)
flagged = np.zeros(len(points), bool)
flagged[np.argsort(iso_score)[-40:]] = True

fig, axes = style.subplots(2, height=3.6, width=10)
axes[0].scatter(*points.T, s=5, color=style.GREY)
for case, color in ((edge_case, style.CORAL), (middle_case, style.BLUE)):
    cuts = isolate(points, case, np.random.default_rng(3), box)
    draw_cuts(axes[0], cuts, color)
    axes[0].scatter(*points[case], s=60, color=color, zorder=3)
    axes[0].text(*points[case] + 0.15, f"{len(cuts)} cuts", color=color, fontsize=9)
axes[0].set(title="One isolation tree", xticks=[], yticks=[])
score_map(axes[1], points, forest.score_samples, "Isolation forest", flagged, flip=True)
fig.tight_layout()
plt.show()
Figure 7.3.3: Left: one random isolation tree, drawn as its cuts of the plane. The isolated point at the edge (red) is cut off after a few cuts; a point in the middle of the crowd (blue) needs many more. Right: isolation forest scores (200 trees) for the data of Figure 2.3.3.2; darker regions are easier to isolate. Blue circles: the 40 highest-scoring cases.

7.3.4 Local density: the local outlier factor

Some anomalies are only unusual relative to their neighbourhood: a point slightly outside a tight cluster can be more suspicious than a point in a loose cluster at the same distance. The local outlier factor (LOF) compares the density around each case with the density around its k nearest neighbours. A ratio near 1 means the case fits its neighbourhood; well above 1, it sits in a sparser spot than its neighbours.

Code 
tight = rng.normal([0, 0], 0.25, (150, 2))
loose = rng.normal([5, 4], 1.4, (150, 2))
near_tight = np.array([[1.2, 0.2], [-0.4, 1.3], [0.9, -1.0]])
mixed = np.r_[tight, loose, near_tight]
lof = LocalOutlierFactor(n_neighbors=20, novelty=True).fit(mixed)
iso = IsolationForest(n_estimators=200, random_state=0).fit(mixed)

fig, axes = style.subplots(2, height=3.6, width=10)
for ax, scores, title in (
    (axes[0], -iso.score_samples(mixed), "Isolation forest"),
    (axes[1], -lof.score_samples(mixed), "Local outlier factor"),
):
    top = np.zeros(len(mixed), bool)
    top[np.argsort(scores)[-10:]] = True
    ax.scatter(*mixed[~top].T, s=6, color=style.INK, alpha=0.6)
    ax.scatter(*near_tight.T, s=20, color=style.CORAL, zorder=3)
    ax.scatter(*mixed[top].T, s=60, facecolors="none", edgecolors=style.BLUE, lw=1.5)
    ax.set(title=title, xticks=[], yticks=[])
fig.tight_layout()
plt.show()
Figure 7.3.4: A tight cluster and a loose cluster, with three points just outside the tight one (red). Left: the isolation forest (like any global distance) flags mostly the edges of the loose cluster. Right: the local outlier factor (20 neighbours) compares each point with its own neighbourhood and flags the three points beside the tight cluster first. Blue circles: the 10 highest scores by each method.

7.3.5 How good is a detector?

If some anomalies are known (fraud cases, faulty parts), a detector can be judged like a classifier on a very rare class (Chapter 6.1): with a precision-recall curve, or by the share of real anomalies among the top k flagged cases, which is what an investigator with time for k cases cares about.

Code 
scaled = ((wine.data - wine.data.mean()) / wine.data.std()).to_numpy()
cultivar = wine.target.to_numpy()
flavanoids = list(wine.data.columns).index("flavanoids")


def flipped() -> tuple[np.ndarray, np.ndarray]:
    """10 wines with the sign of their flavanoid value flipped."""
    data = scaled.copy()
    victims = rng.choice(np.flatnonzero(np.abs(scaled[:, flavanoids]) > 1), 10, replace=False)
    data[victims, flavanoids] *= -1
    return data, np.isin(np.arange(len(data)), victims)


def exaggerated() -> tuple[np.ndarray, np.ndarray]:
    """10 wines pushed 2.5 times as far from their own cultivar's average profile."""
    data = scaled.copy()
    victims = rng.choice(len(scaled), 10, replace=False)
    for i in victims:
        centre = scaled[cultivar == cultivar[i]].mean(axis=0)
        data[i] = centre + 2.5 * (scaled[i] - centre)
    return data, np.isin(np.arange(len(data)), victims)


def detector_scores(data: np.ndarray) -> dict[str, np.ndarray]:
    """Anomaly scores (higher = more anomalous) from the three detectors."""
    return {
        "robust Mahalanobis": MinCovDet(random_state=0).fit(data).mahalanobis(data),
        "isolation forest": -IsolationForest(n_estimators=300, random_state=0)
        .fit(data)
        .score_samples(data),
        "local outlier factor": -LocalOutlierFactor(n_neighbors=20)
        .fit(data)
        .negative_outlier_factor_,
    }


fig, axes = style.subplots(2, height=3.3, width=10.5)
for ax, (data, truth), title in zip(
    axes, (flipped(), exaggerated()), ("Broken combinations", "Exaggerated wines")
):
    for (name, score), color in zip(
        detector_scores(data).items(), (style.BLUE, style.GREEN, style.GOLD)
    ):
        precision, recall, _ = precision_recall_curve(truth, score)
        ax.step(recall, precision, where="post", color=color, lw=2, label=name)
    ax.axhline(truth.mean(), color=style.INK, ls="--", lw=1)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(
        title=title,
        xlabel="recall (share of anomalies found)",
        ylabel="precision",
        ylim=(0, 1.02),
    )
axes[0].legend(fontsize=8)
fig.tight_layout()
plt.show()
Figure 7.3.5: Three detectors on the 178 wines with 10 hidden anomalies of two kinds, as precision-recall curves (dashed: chance). Left: wines whose flavanoid value was flipped in sign, an ordinary value on its own that breaks its strong link with the other phenol measurements. The robust Mahalanobis distance, which judges combinations against one ellipse, finds them; the isolation forest and LOF barely do. Right: wines pushed two and a half times further from the average profile of their own cultivar, unusual in several measurements at once. Now the isolation forest, which does not assume a single elliptical cloud (the wines form three clusters), does best, and the Mahalanobis distance has lost its lead.

Anomaly detection in time series looks for values far from what the series’ own pattern predicts (residuals from a seasonal model, Chapter 9.1) or for points where the pattern itself changes (change points, Chapter 9.3).

7.3.6 In one picture

  • Anomalies inflate the spread they are judged against (masking); use median and MAD, and robust covariance (minimum covariance determinant).
  • Mahalanobis distance finds unusual combinations of ordinary values.
  • Isolation forests score how few random cuts isolate a case.
  • The local outlier factor compares each case with its own neighbourhood.
  • Judge detectors on known anomalies with precision-recall or precision among the top k; the best method depends on what the anomalies look like.

Next: Part 8, deep learning, networks that learn their own measurements.