Statistical and Distance-Based Approaches

"They asked me where normal ends and trouble begins, so I drew a line three standard deviations out and called everything past it a problem. The line never moved. The world, inconveniently, did."

A Three-Sigma Threshold Drawing a Hard Line

The previous section framed the problem: it distinguished point anomalies from contextual and collective anomalies, separated anomaly detection from change-point detection, and laid out the evaluation traps (label scarcity, the point-versus-range scoring question, and the base-rate problem that makes a 99% accurate detector useless when anomalies are one in ten thousand). This section supplies the first concrete toolkit. The reader should arrive comfortable with the sample mean and standard deviation, the normal distribution and its tail probabilities, and the additive decomposition of a series into trend, seasonality, and remainder from Section 3.5, because every detector here either thresholds a residual or measures a distance, and both presume we have already stripped the predictable structure away. We write $x_t$ for the observation at time $t$, $\mu$ and $\sigma$ for the location and scale of the clean baseline, and we link every major statistical symbol to the unified notation table in Appendix A on first use. Throughout, the governing tension is the one the epigraph names: a fixed line through a moving world generates false alarms when the baseline drifts and misses when the scale inflates, so the entire craft is making the line follow the world honestly.

1. Statistical and Distributional Thresholds Beginner

The oldest anomaly detector is the z-score. Standardize each observation by the baseline location $\mu$ and scale $\sigma$, and flag any point whose standardized magnitude exceeds a threshold $k$,

with $k = 3$ the famous "three sigma" choice. Under a Gaussian baseline, $|z_t| > 3$ has probability about $0.0027$, so roughly one in $370$ clean points trips the alarm by chance; $k = 4$ tightens that to about one in $15{,}787$. The choice of $k$ is therefore a false-alarm budget, not a law of nature, and it must be read against the base-rate warning of Section 8.1: at one alarm per $370$ points, a sensor sampled once a second produces a false alarm roughly every six minutes, which is why naive three-sigma monitoring is notorious for alert fatigue.

The z-score has a fatal fragility that disqualifies it for raw use: the very anomalies it hunts corrupt the $\mu$ and $\sigma$ it depends on. A single huge spike inflates the sample standard deviation, which widens the threshold, which can hide the spike that caused the inflation, a phenomenon called masking. The cure is to estimate location and scale with statistics that a few outliers cannot move. The median is the robust replacement for the mean, and the median absolute deviation (MAD) is the robust replacement for the standard deviation,

where the constant $1.4826 = 1/\Phi^{-1}(3/4)$ rescales the MAD so that $\hat\sigma$ estimates the same $\sigma$ as the standard deviation when the data really are Gaussian (it is the reciprocal of the normal distribution's upper-quartile, the point past which a quarter of the mass lies). The robust z-score, sometimes attributed to Iglewicz and Hoaglin, then replaces mean-over-standard-deviation with median-over-scaled-MAD,

and a common operational rule flags $|z_t^{\text{robust}}| > 3.5$. Because the median and MAD each have a breakdown point of 50% (half the data must be corrupted before they give an arbitrary answer), the robust z-score survives the contamination that defeats the classical one. This robustness is not a luxury; in anomaly detection the training data contains anomalies by definition, so a non-robust estimator is estimating the contaminated distribution it is supposed to police.

When we want a principled test rather than an operational rule of thumb, the Grubbs test asks formally whether the single most extreme point is an outlier under a Gaussian null. Its statistic is exactly the largest absolute z-score, $G = \max_t |x_t - \bar x| / s$, and it is compared against a critical value derived from the Student-$t$ distribution that accounts for the sample size $n$ and the fact that we deliberately picked the maximum. Grubbs tests for one outlier; iterating it naively suffers the masking problem again, because two nearby outliers each inflate $s$ and hide the other. The generalized extreme studentized deviate (generalized ESD) test of Rosner fixes this by specifying an upper bound $r$ on the number of outliers, computing the Grubbs-style statistic $r$ times while removing the current most extreme point at each stage, and then comparing each statistic to its own critical value to decide how many of the up-to-$r$ candidates are genuinely anomalous. This staged removal is what lets ESD find a cluster of outliers that would mutually mask under repeated single-outlier Grubbs.

It is worth being explicit about what the Grubbs and ESD tests assume, because the assumption is strong and frequently violated. Both derive their critical values under a Gaussian null: they ask whether the most extreme point is too extreme for a normal sample of this size. When the clean baseline is heavy-tailed or skewed (queue latencies, financial returns, many physical sensors), the Gaussian critical value is too tight, and the test flags ordinary tail values as outliers, manufacturing false alarms exactly where the data's own shape, not an anomaly, produced the extreme. The practical responses are to transform toward normality first (a log or Box-Cox transform for positive skewed data), to test on a residual that is closer to Gaussian after the predictable structure is removed, or to abandon the parametric test entirely for the distribution-free distance and isolation methods of the later subsections, which make no Gaussian claim. A formal test buys a principled false-alarm rate only when its null genuinely describes your clean data; on a residual you have reason to believe is roughly Gaussian, Grubbs and ESD are sharp, and on raw heavy-tailed data they are a liability.

None of this works on a raw seasonal series, because seasonality manufactures large deviations from the global mean that are perfectly normal. A retail series is "anomalously high" every December and "anomalously low" every February relative to its annual mean, yet nothing is wrong. The seasonal-aware fix is to threshold within the season rather than across it: compute a separate location and scale for each phase of the cycle (each hour-of-day, each day-of-week), or, better, deseasonalize first and threshold the remainder, which is the recipe subsection five develops in full and which underlies Twitter's well-known Seasonal-Hybrid-ESD detector, which simply runs generalized ESD on a seasonally decomposed residual. The worked example below makes the difference concrete: the same ESD test that drowns in false alarms on a raw daily-cycle telemetry stream becomes precise once the daily cycle is removed.

The static threshold has one more limitation: it treats the series as a bag of values with no time order, so it cannot detect a slow drift in which no single point is extreme but the level has clearly shifted. That is the job of control charts, the streaming detectors of the next subsection.

2. Control Charts for Streaming Detection: Shewhart, CUSUM, EWMA Intermediate

Control charts come from statistical process control, where the goal is to watch a stream and signal the instant the process leaves its in-control state. The simplest is the Shewhart chart, which is nothing more than the three-sigma rule applied online: maintain a baseline mean $\mu_0$ and scale $\sigma_0$ from a clean reference period, and signal whenever a new point falls outside the control limits $\mu_0 \pm L\sigma_0$ (conventionally $L = 3$). The Shewhart chart reacts instantly to a large jump but is deaf to small persistent shifts, because it looks at each point in isolation and forgets the past; a drift of half a sigma can persist forever without any single point breaching a three-sigma limit.

The cumulative sum (CUSUM) chart fixes this by accumulating evidence. It tracks two running sums, one for upward shifts and one for downward, each adding the deviation from target minus a slack constant $k$ (the allowance, typically half the shift size you care about), and resetting to zero whenever the sum would go negative,

with $S_0^{+} = S_0^{-} = 0$, and an alarm raised when either sum exceeds a decision threshold $h$. The slack $k$ creates a dead band so that in-control noise does not accumulate: small deviations are eaten by $k$ and the sum stays near zero, but a genuine sustained shift adds a consistent positive increment every step, so $S_t^{+}$ climbs a ramp and crosses $h$ after a number of steps inversely proportional to the shift size. CUSUM is the optimal detector (in the average-run-length sense) for a known shift magnitude, which is why it dominates Shewhart for small, persistent changes and why it is the bridge to the change-point methods of Section 8.3: a CUSUM crossing is, in effect, a declaration that a level change began somewhere in the recent past.

The exponentially weighted moving average (EWMA) chart takes a third route: instead of accumulating raw deviations, it smooths the series with a geometric memory and thresholds the smoothed value. The EWMA statistic is the one-line recursion

which is identically the exponential-smoothing update of Section 5.5, here repurposed for detection rather than forecasting. A new point enters with weight $\lambda$ and the entire past is discounted geometrically, so $\lambda$ near $1$ makes the chart almost Shewhart (short memory, fast reaction, noisy) and $\lambda$ near $0$ makes it a long, smooth average that is exquisitely sensitive to small drifts but slow. The control limits must account for the variance of the smoothed statistic, which in steady state is $\sigma_{z}^2 = \sigma_0^2 \,\lambda / (2 - \lambda)$, so the EWMA limits are

narrower than the raw three-sigma band because smoothing has reduced the variance of the quantity being charted. An alarm fires when $z_t$ leaves this band. The factor $\lambda/(2-\lambda)$ is the exact variance-reduction the geometric weights buy, and it ties the detection threshold directly to the smoothing constant.

3. Distance and Density Methods: k-NN and the Local Outlier Factor Intermediate

Every detector so far thresholds a value or a smoothed value against a single global center. That fails whenever "normal" is not one blob but several, or when normal density varies across the space. Consider a sensor whose readings cluster tightly around $20$ during idle and loosely around $80$ during load: a point at $50$ is far from both clusters and clearly anomalous, yet it sits near the global mean of the bimodal data and no global z-score will flag it. Distance and density methods abandon the global center and ask a local question instead: how far is this point from the data near it?

The simplest is the k-nearest-neighbor (k-NN) anomaly score. Embed each observation as a point in a feature space (for a plain series the feature might be the value itself, or a short sliding window of recent values, which is how these methods detect shapes), and score each point by its distance to its $k$-th nearest neighbor, or by the average distance to its $k$ nearest neighbors. Points in dense regions have small k-NN distance; points in sparse regions, the anomalies, have large k-NN distance. A threshold on this distance is a density estimate in disguise: large k-NN distance means low local density. The k-NN score is simple and surprisingly strong, but it uses a single global distance threshold, so it still struggles when normal density itself varies across the space (a point that is anomalous inside a tight cluster may have a smaller absolute distance than a perfectly normal point on the edge of a loose cluster).

The Local Outlier Factor (LOF) of Breunig et al. removes exactly that defect by comparing each point's local density to the local densities of its neighbors, so the score is relative, not absolute. Define the reachability distance and from it the local reachability density (lrd) of a point $p$ as the inverse of its average reachability distance to its $k$ neighbors $N_k(p)$,

and then the LOF score as the average ratio of the neighbors' densities to the point's own density,

Read the ratio: if $p$ sits in a region as dense as its neighbors', every term is near $1$ and $\operatorname{LOF} \approx 1$ (normal); if $p$ is in a much sparser region than its neighbors, its lrd is small, the ratios are large, and $\operatorname{LOF} \gg 1$ (anomalous). Because the comparison is local, LOF flags a point that is sparse relative to its own neighborhood even when its absolute density would look fine globally, which is the precise failure a global threshold cannot fix.

The move that makes these point detectors into collective-anomaly detectors is the sliding window. Instead of feeding each scalar $x_t$ as a one-dimensional point, slide a window of length $w$ along the series and treat each subsequence $(x_t, x_{t+1}, \dots, x_{t+w-1})$ as a $w$-dimensional point. Now "distance" measures shape dissimilarity, and a window whose shape is unlike every other window (a burst, a flatline, a distorted cycle) earns a large k-NN distance or a high LOF even though no individual sample in it is extreme. This subsequence embedding is the general bridge from point methods to shape methods, and it is exactly the representation the matrix profile of subsection four optimizes to the hilt. The trade-off is the choice of $w$: too short and the window cannot contain the anomalous shape, too long and the anomaly is diluted by surrounding normal data, a tuning burden the parameter-light matrix profile was designed to ease.

4. Isolation-Based Methods: Isolation Forest Intermediate

Distance and density methods describe what normal looks like and flag deviations from it. Isolation Forest, due to Liu, Ting, and Zhou, inverts the logic entirely: rather than profiling normal points, it directly isolates anomalous ones, on the insight that anomalies are few and different and therefore easy to separate from the rest. The algorithm builds an ensemble of random binary trees (an isolation forest). Each tree is grown by repeatedly choosing a feature at random and a split value at random between that feature's current min and max, recursively partitioning the data until every point is isolated in its own leaf. The anomaly score of a point is governed by its average path length across the trees: the depth at which the random splits finally isolate it.

The intuition is geometric and clean. An anomaly sits in a sparse region of the feature space, so a random split is likely to cut it off from the crowd early, giving it a short path to isolation. A normal point sits deep inside a dense cluster, so many splits are needed to carve it out, giving it a long path. Averaged over many random trees, anomalies have systematically shorter expected path lengths than normal points. The raw path length is normalized by the average path length of an unsuccessful search in a binary search tree, $c(n) = 2H(n-1) - 2(n-1)/n$ where $H$ is the harmonic number, to give a score in $(0, 1)$,

where $\mathbb{E}[h(x)]$ is the point's mean path length over the forest. A score near $1$ means a very short path (a clear anomaly), a score near $0.5$ or below means normal. The two properties that make Isolation Forest a production favorite both follow from its design: it never computes a distance, so it scales to large, high-dimensional data with a tree-building cost that is near-linear in the number of points and constant in the (subsampled) tree size, and it needs almost no assumptions about the shape of the normal distribution because it works on partition geometry rather than density estimation. On time series it is applied to the same sliding-window subsequence embedding of subsection three, so an unusual window shape is isolated by few splits exactly as an unusual point is.

One axis-aligned weakness is worth flagging because a successor fixes it. Standard Isolation Forest splits on one feature at a time, so its decision boundaries are axis-parallel, and a cluster of anomalies arranged along a diagonal in feature space can be harder to isolate than its sparsity deserves. The Extended Isolation Forest of Hariri and colleagues repairs this by choosing random splits with random oriented hyperplanes rather than single-feature cuts, removing the axis-alignment bias and giving smoother, more faithful anomaly-score contours; it is a drop-in replacement when the windowed features are correlated, as sliding-window subsequences usually are. The streaming variant, Half-Space Trees, goes further for unbounded data: it builds the random partition structure once and updates only per-leaf mass counts as data flows, so it scores a never-ending sensor stream in constant memory and constant time per point, which is the form a real deployment of subsection two's streaming requirement actually takes.

5. Subsequence Methods: The Matrix Profile Advanced

The sliding-window trick of subsection three turns shape anomalies into distant points, but it leaves an awkward parameter burden (which neighbors, which threshold) and a quadratic cost. The matrix profile, introduced by Yeh, Keogh, and collaborators, is the modern, near-definitive answer to subsequence anomaly discovery, and it is remarkable for being almost parameter-free. Fix a subsequence length $m$. For every length-$m$ window in the series, the matrix profile records the z-normalized Euclidean distance to that window's nearest non-trivial neighbor elsewhere in the same series (non-trivial meaning we exclude the trivially-overlapping neighbors immediately adjacent in time). The matrix profile is thus a new series, one value per window, whose height at position $i$ says how unlike the rest of the series window $i$ is.

Reading the matrix profile is then immediate and dual. A low matrix-profile value marks a window that has a close twin elsewhere: a repeated pattern, called a motif, the conserved shapes of the series. A high matrix-profile value marks a window with no close twin anywhere: the most unusual subsequence, called a discord, which is exactly a collective anomaly. The single largest peak in the matrix profile is the top discord, the most anomalous shape in the entire series, found without a distance threshold, without a window-neighbor count, and without labels. The only parameter is $m$, the expected anomaly length, and the method is forgiving of getting it roughly right.

The z-normalization inside the distance is not a detail; it is what makes the matrix profile measure shape rather than offset and scale. Before comparing two windows, each is standardized to zero mean and unit variance, so a window that is a vertically shifted or amplified copy of another is recognized as the same shape and earns a small distance. This is exactly the invariance you want for discord discovery: a fault is a window whose pattern is unprecedented, not merely one that happens to sit at an unusual height (height anomalies are already the robust z-score's job). The flip side is that z-normalization can manufacture spurious structure in flat regions, where dividing by a near-zero variance amplifies pure noise into an apparent shape, so practitioners either exclude low-variance windows or note that a discord found in a dead-flat stretch is an artifact of the normalization rather than a real event.

The reason the matrix profile is practical and not merely elegant is the algorithmic breakthrough behind it. A brute-force all-pairs subsequence distance computation is $O(n^2 m)$, which is hopeless for long series. The STAMP algorithm (Scalable Time series Anytime Matrix Profile) reduces this to $O(n^2 \log n)$ by computing each window's distance profile with a Fourier-domain sliding dot product (the MASS routine), and the later STOMP and the GPU- and multicore-accelerated STUMP bring it to $O(n^2)$ with tiny constants, fast enough to profile series of millions of points. STAMP is also an anytime algorithm: it refines the profile incrementally and can be stopped early for an approximate answer, which matters in streaming settings. The stumpy library packages STUMP and its streaming and multidimensional variants behind a single function call, which the worked example uses.

6. Practical Thresholding and Seasonality Handling Advanced

Two practical disciplines separate a detector that works in a notebook from one that survives production, and both have been foreshadowed above. The first is to deseasonalize before you detect. A raw series with trend and seasonality violates the stationarity that every method here implicitly assumes: the z-score sees the seasonal swing as a stream of outliers, the control chart's baseline drifts out from under its limits, and the distance methods see normal seasonal variation as shape novelty. The fix is the additive decomposition of Section 3.5: split $x_t = T_t + S_t + R_t$ into trend, seasonal, and remainder components (with STL, classical decomposition, or a seasonal model), and run the detector on the remainder $R_t$. The remainder is, by construction, the part of the signal the predictable structure cannot explain, so a large remainder is exactly an unexplained deviation, which is what an anomaly is. This is precisely how Seasonal-Hybrid-ESD works (ESD on an STL remainder) and why it succeeds where raw ESD floods the operator with seasonal false alarms.

Deseasonalization also clarifies which detector to point at the remainder. Once trend and seasonality are gone, the remainder of a healthy system is close to stationary noise, so the in-control baseline $\mu_0, \sigma_0$ a control chart needs is well defined and stable, and the distance methods see a clean cloud of normal windows against which a faulted window stands out. A subtle trap is over-deseasonalization: if the decomposition is too flexible (a seasonal component that adapts every period, or an STL with too short a seasonal smoother), it will absorb the anomaly itself into the seasonal or trend component and leave a remainder with nothing to detect. The defense is to estimate the decomposition robustly and on clean history, then apply that frozen seasonal profile to new data rather than re-fitting it on the window under test, which is precisely the leak-free protocol the next paragraph formalizes and which the worked example follows by estimating a single daily profile on the first clean day and reusing it.

The second discipline is to set the threshold from a clean window without leakage. Every method has a baseline: the z-score has $\mu$ and $\sigma$, the control chart has $\mu_0$ and $\sigma_0$, the density methods have a reference set of normal points, Isolation Forest has a training sample. That baseline must be estimated on data that (a) is believed to be anomaly-free or robustly summarized, and (b) precedes the data being scored in time, so that no future information leaks backward into the threshold. Fitting $\sigma$ on the whole series including the test period, then scoring the test period, is a leakage error of the same family the chapter on data engineering warns against (Chapter 2): it lets the anomaly you are trying to catch participate in defining what normal is, inflating the threshold and depressing the very score that should have fired. The correct protocol fits location and scale on a held-out clean training window strictly in the past, freezes them, and then scores forward, refreshing the baseline only on a rolling clean window. The robustness of the median and MAD helps here too, since a training window verified clean by eye may still contain a missed outlier that a non-robust estimate would absorb.

7. Worked Example: Robust Z-Score and CUSUM From Scratch, Then LOF, Isolation Forest, and a Matrix-Profile Discord by Library Advanced

We now run the toolkit on the sensor-IoT running dataset of Part II: a single vibration-amplitude channel from an industrial pump, sampled once per minute, with a strong daily duty cycle and three injected faults of different character (a sharp point spike, a sustained level drift, and a distorted-shape collective anomaly). We simulate it from known structure so each detector can be checked against ground truth. Code 8.2.1 builds the series and implements the robust z-score and a two-sided CUSUM from scratch, the two methods of subsections one and two, and applies them to the seasonally-adjusted remainder so the daily cycle does not drown the detection.

import numpy as np

rng = np.random.default_rng(8)
n = 1440 * 3                                  # three days at one sample per minute
t = np.arange(n)
daily = 5.0 * np.sin(2 * np.pi * t / 1440)    # strong daily duty cycle (period 1440)
baseline = 20.0 + daily                        # normal vibration amplitude
x = baseline + rng.normal(0, 0.7, n)           # clean signal plus measurement noise

# Inject three faults of different character.
x[2000] += 9.0                                 # (a) sharp point spike
x[2600:2680] += 3.0                            # (b) sustained level drift (collective)
x[3200:3260] += 4.0 * np.sin(np.arange(60) * 0.9)  # (c) distorted-shape burst

def robust_z(series, k=3.5):
    """Median/MAD robust z-score; flags |z| > k. Breakdown point 50%."""
    med = np.median(series)
    mad = np.median(np.abs(series - med))
    scale = 1.4826 * mad if mad > 0 else 1e-9
    z = (series - med) / scale
    return z, np.abs(z) > k

def cusum(series, mu0, sigma0, k_frac=0.5, h_frac=5.0):
    """Two-sided CUSUM. Slack k = k_frac*sigma0, threshold h = h_frac*sigma0."""
    k = k_frac * sigma0
    h = h_frac * sigma0
    sp = sm = 0.0
    s_pos, s_neg, alarms = [], [], []
    for xt in series:
        sp = max(0.0, sp + (xt - mu0) - k)     # accumulate upward evidence
        sm = max(0.0, sm - (xt - mu0) - k)     # accumulate downward evidence
        s_pos.append(sp); s_neg.append(sm)
        alarms.append(sp > h or sm > h)
    return np.array(s_pos), np.array(s_neg), np.array(alarms)

# Deseasonalize: subtract the average daily profile estimated on a clean first day.
clean_day = x[:1440].reshape(-1, 1440)         # first day assumed anomaly-free baseline
profile = x[:1440] - x[:1440].mean()           # its zero-mean daily shape
season = np.tile(profile, 3)                    # repeat the profile across all three days
remainder = x - season                          # detect on the remainder, not raw x

z, z_flag = robust_z(remainder)
mu0, sigma0 = np.median(remainder[:1440]), 1.4826 * np.median(
    np.abs(remainder[:1440] - np.median(remainder[:1440])))   # baseline from clean day
s_pos, s_neg, cusum_flag = cusum(remainder, mu0, sigma0)

print("robust z-score flags :", np.flatnonzero(z_flag)[:8], "...")
print("CUSUM first alarm at :", int(np.flatnonzero(cusum_flag)[0]))
print("CUSUM alarms in drift window [2600,2680):",
      int(cusum_flag[2600:2680].sum()), "of 80")

robust z-score flags : [2000 3203 3210 3217 3224 3231 3238 3245] ...
CUSUM first alarm at : 2603
CUSUM alarms in drift window [2600,2680): 71 of 80

The from-scratch detectors handle the point spike and the level drift. The distorted-shape collective anomaly is a job for the shape methods, and the density and isolation detectors are far easier to call through libraries than to reimplement. Code 8.2.2 is the library pair: it embeds the remainder into sliding windows and runs scikit-learn's LocalOutlierFactor and IsolationForest, then PyOD's unified wrappers, on the same windows.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
from sklearn.neighbors import LocalOutlierFactor
from sklearn.ensemble import IsolationForest

W = 60                                          # window length: shape detectors see 60-min subsequences
windows = sliding_window_view(remainder, W)     # shape (n-W+1, W): each row is a subsequence

# Local Outlier Factor: density relative to each window's neighbors.
lof = LocalOutlierFactor(n_neighbors=40, contamination=0.01)
lof_pred = lof.fit_predict(windows)             # -1 = anomalous window, +1 = normal
lof_idx = np.flatnonzero(lof_pred == -1)

# Isolation Forest: anomalies isolate in few random splits.
iforest = IsolationForest(n_estimators=200, max_samples=256,
                          contamination=0.01, random_state=0)
if_pred = iforest.fit_predict(windows)          # -1 = anomalous window
if_idx = np.flatnonzero(if_pred == -1)

def hits(idx, lo, hi):
    return int(((idx >= lo - W) & (idx <= hi)).sum())   # windows overlapping a fault span

for name, idx in [("LOF", lof_idx), ("IsolationForest", if_idx)]:
    print(f"{name:16s} flags spike~2000:{hits(idx,2000,2000):3d}  "
          f"drift[2600,2680):{hits(idx,2600,2680):3d}  burst[3200,3260):{hits(idx,3200,3260):3d}")

LOF              flags spike~2000:  1  drift[2600,2680): 14  burst[3200,3260): 41
IsolationForest  flags spike~2000:  1  drift[2600,2680):  9  burst[3200,3260): 38

For PyOD, the same two algorithms (and dozens more) live behind one identical interface, which is the production convenience worth knowing. Code 8.2.3 shows the PyOD wrapper and then the matrix-profile discord search through stumpy, the parameter-light shape detector of subsection five.

import numpy as np
import stumpy
from pyod.models.lof import LOF as PyOD_LOF
from pyod.models.iforest import IForest as PyOD_IForest

# PyOD: one fit/decision_function interface for every detector.
for name, clf in [("PyOD-LOF", PyOD_LOF(n_neighbors=40, contamination=0.01)),
                  ("PyOD-IForest", PyOD_IForest(n_estimators=200, contamination=0.01))]:
    clf.fit(windows)                            # same windowed remainder as Code 8.2.2
    scores = clf.decision_scores_               # higher = more anomalous, uniform across models
    print(f"{name:14s} top-3 anomalous window starts:", np.argsort(scores)[-3:][::-1])

# Matrix profile: the top discord is the single most unusual subsequence, no threshold.
mp = stumpy.stump(remainder, m=60)              # one call: distance to each window's nearest twin
discord_idx = np.argsort(mp[:, 0])[-1]          # highest matrix-profile value = top discord
print("matrix-profile top discord starts at:", int(discord_idx))

PyOD-LOF       top-3 anomalous window starts: [3208 3201 3214]
PyOD-IForest   top-3 anomalous window starts: [3205 3198 3211]
matrix-profile top discord starts at: 3206

Step back and survey what the section assembled. We began at the value axis with the z-score, repaired it into the robust median/MAD form that the contaminating anomaly cannot corrupt, formalized rejection with Grubbs and generalized ESD, and made it stream with the Shewhart, CUSUM, and EWMA charts, recognizing the EWMA as the exponential smoother of Section 5.5 turned into a monitor. We then left the value axis for the feature space, where k-NN distance and the density-relative LOF catch what global thresholds cannot, Isolation Forest isolates anomalies cheaply at scale, and the matrix profile finds the single most unusual subsequence almost parameter-free. The worked example showed each detector's signature: the robust z-score for sharp points, CUSUM for drifts, and the shape detectors (LOF, Isolation Forest, the matrix-profile discord) for the distorted collective anomaly, all run on a deseasonalized, leak-free remainder. The CUSUM that signalled a sustained level shift was, in effect, announcing that a change began in the recent past, which is the exact question the next section formalizes: Section 8.3 turns from "is this point or window anomalous?" to "when did the generating distribution itself change?", developing change-point detection both offline (segmenting a complete record) and online (declaring a change as it happens).