Instrument Control

Rolling Z-Score Anomaly Detection for Sensor Drift

A fixed absolute limit cannot catch a sensor that is slowly walking away from calibration, because the number it should compare against is itself moving. Rolling z-score detection solves that specific problem: it scores each sample against a baseline computed from a sliding window, so a photodiode drifting with lamp age or an RTD creeping with ambient temperature trips a flag the moment it departs from its own recent behaviour. This guide builds an online detector for that case, a companion to Threshold Tuning & Alerting within Data Capture, Validation & Metadata Sync.

Where a Static Limit Fails and a Windowed Score Wins

The scope here is deliberately narrow: detecting slow drift and transient outliers in a single scalar channel where the baseline is non-stationary. A spectrometer’s dark current climbs over a run, a mass-flow controller’s zero-point wanders with manifold temperature, an electrochemical probe fouls gradually — in every case the value five minutes from now is a poor match for a limit chosen at instrument start-up. A static upper = 4.2 mA either fires the moment the baseline shifts under it or, if set loose enough to survive the shift, stops catching anything worth catching.

The windowed z-score reframes the question. Instead of asking “is the value above an absolute line?” it asks “is the value far from what this channel has been doing over the last N samples, measured in units of that window’s own noise?” The band travels with the signal, so genuine drift only trips when the rate of departure outruns the window’s ability to absorb it, and a sudden spike trips immediately because it is many sigma from a baseline that has not yet moved. This detector produces a numeric score and a drift flag; the downstream logic that debounces and holds an alert state — Hysteresis Thresholds to Suppress Alert Flapping and Debouncing Instrument Alarms with Time Windows — is a separate concern layered on top of the score this page computes. Every sample it sees is assumed already integrity-checked, having cleared Checksum & CRC Validation; a corrupted reading fed into the window poisons μ and σ for the entire window length.

The Rolling Z-Score, Welford Updates, and the Robust Variant

For a window W of the last N samples with mean μ_w and standard deviation σ_w, the score for the current sample x_t is its standardised distance from the windowed baseline:

Recomputing μ and σ from scratch every sample is O(N) and needless. Welford’s algorithm maintains the mean and the sum of squared deviations M2 incrementally in O(1) with a single pass and far better numerical stability than the textbook E[x²] − E[x]², which suffers catastrophic cancellation when the variance is small relative to the mean — exactly the regime of a stable sensor. The update on adding a sample is:

A sliding window also has to remove the sample that fell off the tail, which the forward recurrence does not express. The detector below applies the exact inverse update on eviction and periodically recomputes from the live deque to bound the floating-point error that accumulates over millions of add/remove pairs.

The plain z-score has a self-defeating weakness: the very outlier you want to detect inflates both μ and σ of the window it enters, dragging the baseline toward the anomaly and shrinking the score. The robust remedy replaces the mean with the median and σ with the median absolute deviation (MAD), giving the modified z-score, where the constant 0.6745 rescales MAD to be consistent with σ for normally distributed data:

Median and MAD have a breakdown point of 50%: up to half the window can be contaminated before the estimate is destroyed, versus a single bad sample for the mean. The cost is O(N log N) per sample for the median, acceptable for the modest windows (128–1024) typical of instrument channels. Choosing between the two is a real trade-off — the classical z-score is cheaper and more responsive to genuine step changes, the modified z-score is far harder to fool with spiky data. The detector supports both.

Rolling z-score band tracking a drifting signal A horizontal plot of a sensor signal that drifts upward over time. A shaded band bounded by an upper and lower line runs diagonally upward, centred on a dashed mean line that rises with the signal; the band represents plus and minus k sigma computed from a sliding window. The jagged signal stays inside the band while it drifts, so slow drift alone does not trip the detector. A single sharp spike near the middle punches above the upper band and is marked as an outlier where the absolute z-score exceeds k. A translucent rectangle over a slice of the signal marks the current rolling window of N samples from which the mean and sigma are recomputed as it slides right. time / sample index → value μ_w + kσ_w μ_w − kσ_w rolling mean μ_w outlier |z| > k rolling window (N) slides → slow drift stays inside the band — the baseline moves with it

The band is recomputed from the last N samples and rides the drifting mean, so gradual drift alone does not fire; only a departure that outruns the window — a spike, or drift faster than the window absorbs — escapes the ±kσ_w envelope.

A Production RollingZScore Detector

The detector maintains a bounded deque, updates the Welford statistics on every add and eviction, and exposes update(x) -> float returning the signed score alongside a drift flag raised only when the score stays beyond the threshold for a configurable run of samples — the criterion that separates sustained drift from a lone spike. A single-sample excursion sets outlier but not drift; a persistent one-sided departure sets both. NumPy is used only for the robust median/MAD path and degrades to pure Python when unavailable.

from __future__ import annotations

import math
from collections import deque
from dataclasses import dataclass
from enum import Enum
from typing import Deque, Optional

try:
    import numpy as np
except ImportError:  # numpy optional; only the robust path needs it
    np = None  # type: ignore[assignment]


class Method(str, Enum):
    ZSCORE = "zscore"        # Welford mean/std, responsive to step changes
    MODIFIED = "modified"    # median/MAD, robust to contamination


@dataclass(frozen=True)
class Score:
    z: float
    outlier: bool
    drift: bool
    mean: float
    scale: float  # sigma (zscore) or MAD-derived scale (modified)


class RollingZScore:
    """Online rolling-window anomaly detector for non-stationary sensor channels.

    Feed one scalar sample per call to :meth:`update`. Returns a signed score in
    units of the window's own dispersion. ``outlier`` flags a single-sample
    excursion; ``drift`` flags a departure sustained for ``drift_run`` samples.
    """

    def __init__(
        self,
        window: int = 256,
        threshold: float = 3.5,
        method: Method = Method.ZSCORE,
        drift_run: int = 8,
        min_scale: float = 1e-9,
        recompute_every: int = 4096,
    ) -> None:
        if window < 8:
            raise ValueError("window must be >= 8 for a meaningful baseline")
        if method is Method.MODIFIED and np is None:
            raise RuntimeError("Method.MODIFIED requires numpy for median/MAD")
        self._buf: Deque[float] = deque(maxlen=window)
        self._threshold = threshold
        self._method = method
        self._drift_run = drift_run
        self._min_scale = min_scale
        self._recompute_every = recompute_every
        self._n = 0
        self._mean = 0.0
        self._m2 = 0.0        # running sum of squared deviations (Welford)
        self._since_recompute = 0
        self._run_sign = 0     # +1/-1 direction of the current exceedance run
        self._run_len = 0

    def update(self, x: float) -> Score:
        """Ingest one sample and return its rolling anomaly score."""
        if not math.isfinite(x):
            raise ValueError(f"non-finite sample rejected: {x!r}")

        # Score against the CURRENT window, before x is admitted, so a spike is
        # measured against an uncontaminated baseline.
        warm = len(self._buf) >= max(8, self._buf.maxlen // 4)
        if self._method is Method.ZSCORE:
            mean, scale = self._mean, self._welford_std()
        else:
            mean, scale = self._median_mad()

        scale = max(scale, self._min_scale)
        z = (x - mean) / scale if warm else 0.0
        outlier = warm and abs(z) >= self._threshold

        # Drift = a run of same-signed exceedances; a lone spike does not qualify.
        sign = (z > 0) - (z < 0)
        if outlier and sign == self._run_sign:
            self._run_len += 1
        elif outlier:
            self._run_sign, self._run_len = sign, 1
        else:
            self._run_sign, self._run_len = 0, 0
        drift = self._run_len >= self._drift_run

        self._admit(x)
        return Score(z=z, outlier=outlier, drift=drift, mean=mean, scale=scale)

    # --- window maintenance -------------------------------------------------
    def _admit(self, x: float) -> None:
        if len(self._buf) == self._buf.maxlen:
            self._remove(self._buf[0])  # evicted by the append below
        self._buf.append(x)
        self._add(x)
        self._since_recompute += 1
        if self._since_recompute >= self._recompute_every:
            self._recompute()

    def _add(self, x: float) -> None:
        self._n += 1
        delta = x - self._mean
        self._mean += delta / self._n
        self._m2 += delta * (x - self._mean)

    def _remove(self, x: float) -> None:
        if self._n <= 1:
            self._n, self._mean, self._m2 = 0, 0.0, 0.0
            return
        prev_mean = self._mean
        self._mean = (self._mean * self._n - x) / (self._n - 1)
        self._m2 -= (x - prev_mean) * (x - self._mean)
        self._n -= 1
        if self._m2 < 0.0:  # guard tiny negative from rounding
            self._m2 = 0.0

    def _recompute(self) -> None:
        """Rebuild Welford state from the live buffer to bound rounding drift."""
        self._n, self._mean, self._m2 = 0, 0.0, 0.0
        for v in self._buf:
            self._add(v)
        self._since_recompute = 0

    def _welford_std(self) -> float:
        return math.sqrt(self._m2 / (self._n - 1)) if self._n > 1 else 0.0

    def _median_mad(self) -> tuple[float, float]:
        arr = np.fromiter(self._buf, dtype=float)
        med = float(np.median(arr))
        mad = float(np.median(np.abs(arr - med)))
        return med, mad / 0.6745  # rescale MAD to a sigma-equivalent scale

Scoring the sample before admitting it to the window is the load-bearing detail: it measures each reading against a baseline that has not yet been moved by the reading itself, which is what lets a single spike register a large score instead of half-canceling itself.

Validating Detection Latency and False-Positive Rate

Confirm behaviour against a synthetic signal with known ground truth before trusting the detector on a rig: a stationary baseline, a slow linear ramp (drift), and injected spikes (outliers). Two numbers matter — detection latency (samples from drift onset to the drift flag) and false-positive rate on the clean baseline.

import random


def test_ramp_and_spikes() -> None:
    rng = random.Random(42)
    det = RollingZScore(window=256, threshold=3.5, drift_run=8)

    # 400 stationary samples: measure the false-positive rate.
    false_pos = 0
    for _ in range(400):
        s = det.update(rng.gauss(10.0, 0.20))
        false_pos += s.outlier
    assert false_pos <= 4, f"FP rate too high: {false_pos}/400"

    # Inject a spike: should flag outlier, but NOT drift (single sample).
    spike = det.update(10.0 + 8 * 0.20)
    assert spike.outlier and not spike.drift

    # Slow ramp of +0.02/sample: drift must latch within a bounded latency.
    latency: int | None = None
    for i in range(300):
        s = det.update(10.0 + 0.02 * i + rng.gauss(0.0, 0.20))
        if s.drift and latency is None:
            latency = i
    assert latency is not None and latency < 60, f"drift latency {latency}"

On a live instrument the same two indicators apply. Log the score stream and overlay it on the raw channel: a healthy detector shows |z| hovering under 1 during steady operation and rising monotonically as a fault develops. Reconcile flagged events against the instrument’s own diagnostic register where one exists — many DAQ cards expose a saturation or over-range bit that should coincide with your outlier flags. A steady trickle of outliers on a channel you know to be quiet means the threshold is too tight or the window too short; zero flags across a shift where calibration records show drift means the reverse.

Failure Modes Specific to Windowed Z-Score Detection

Variance collapse on a constant input. A sensor sitting on a perfectly flat value drives σ_w (or MAD) to zero, and the next (x − μ)/σ divides by zero, producing inf or nan scores that flag every subsequent sample. The min_scale floor is the guard, but it must be sized to the channel’s quantisation step, not left at a token 1e-9: too small and the first sub-LSB wobble reads as hundreds of sigma. Diagnose by watching for scores that saturate immediately after a period of a stuck reading.

The window swallows the anomaly and masks it. With the plain z-score, a large spike admitted into the window inflates σ_w so much that the next real excursion scores low and slips through — the anomaly raised the bar for detecting the next one. This is exactly why the detector scores before admitting and why the modified z-score exists; on channels prone to bursts, prefer Method.MODIFIED, whose median/MAD baseline is barely moved by a handful of contaminating samples.

Non-stationary noise, not just non-stationary mean. The z-score assumes the dispersion is locally stable even as the mean drifts. When the noise amplitude itself changes — a heater entering turbulent flow, an amplifier warming up — a window sized for the quiet regime floods with false positives once the noise grows. A window sized for the noisy regime goes deaf in the quiet one. Size the window to the noisiest operating regime you must tolerate, and treat a sustained rise in the baseline σ_w as its own signal worth exporting to Threshold Tuning & Alerting.

Quantised sensor gives zero MAD. A coarsely quantised ADC channel that dwells on one or two codes yields a median absolute deviation of exactly zero even when the mean is not constant, because more than half the window shares the modal code. MAD then collapses like σ did above. Detect the condition explicitly — if MAD == 0 over a full window, fall back to the classical σ path or widen the window until at least three distinct codes are present, and record the fallback so a downstream reviewer knows the robust estimator was not in play.

← Back to Threshold Tuning & Alerting

References