CUSUM: sequential analysis for sub-clinical data stream drift detection

This page explains how to detect slow, gradual changes in a system—like a mechanical part wearing out bit by bit. While a standard detector looks at one reading at a time and misses small shifts, CUSUM accumulates these tiny anomalies over time and triggers an alert as soon as the total deviation confirms a real structural trend.

A robot joint wears out a little every day. Each single reading still looks normal, so a detector that judges one sample at a time never complains. Together, though, the readings have moved. CUSUM ("cumulative sum") notices this: it keeps a running total of the small surprises, like coins dropped one by one into a piggy bank, and raises an alarm when the bank is full. For the other three detectors and a benchmark on a real robot arm, see Drift detectors.

1. Watch a stream, one sample at a time

CUSUMDriftDetector takes one value per call and tells you when the drift is real.

import numpy as np
from dense_armor.utility.drift.detector import CUSUMDriftDetector

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, 500), rng.normal(1.5, 1, 500)])
det = CUSUMDriftDetector(reference="fixed")
print(next(i for i, v in enumerate(x) if det.update(v).drift_detected))

511

The first 500 samples are normal noise; from sample 500 the mean moves up by 1.5 standard deviations, the kind of step a worn gear or a new payload produces. The detector raises its first alarm at sample 511, eleven samples after the change. update(v) returns the detector itself, and drift_detected is True only on the sample where the alarm fires. The detector needs the streaming extra; the install command is on the bridge page of the API reference.

2. The piggy bank, as a formula

Every new sample adds its surprise and pays a small fee; the total never goes below zero.

z = [0.2, -0.4, 0.1, 1.8, 2.1, 1.9, 2.3]
k, h, s = 0.5, 5.0, 0.0
for t, zt in enumerate(z, start=1):
    s = max(0.0, s + zt - k)
    print(t, round(s, 2), s > h)

1 0.0 False
2 0.0 False
3 0.0 False
4 1.3 False
5 2.9 False
6 4.3 False
7 6.1 True

The rule is

$$S_t = \max(0,\; S_{t-1} + z_t - k), \qquad \text{alarm when } S_t > h,$$

where:

• $z_t$ is how surprising sample $t$ is, in units of the normal noise (0 means "exactly as usual", 2 means "two noise-widths above usual");
• $k$ is the fee paid at every sample (the slack): small surprises are eaten by the fee and the total stays at zero, so ordinary noise never accumulates;
• $h$ is the size of the piggy bank (the threshold): when the total passes it, the alarm fires and the total starts again from zero;
• $S_t$ is the total after sample $t$.

In the run above the first three values are small, so the fee keeps the total at zero. From the fourth value the readings sit about 2 above normal: each one adds roughly $2 - 0.5 = 1.5$, and the total passes $h = 5$ at the seventh sample. A detector like this also watches the other direction (readings that drop) with a second total, which is what "two-sided" means.

3. Where the surprise $z_t$ comes from

The detector measures each value against a window of the recent past, without looking at the future. The centre of the window is its median and its width is $S = 1.4826 \cdot \mathrm{MAD}$ (the median absolute distance from the median; the factor 1.4826 makes $S$ equal to the standard deviation when the noise is normal, see anomaly detection). Then $z_t = (x_t - \text{median}) / S$.

Two ways to choose the window:

• reference="fixed": the median and width are measured once, on the first radius * ref_mult samples, and kept. Best for a sustained step, like the run in step 1.
• reference="adaptive": they are re-measured on a sliding window at every sample. The window follows a new level, so the detector reacts to the leading edge of a drift and then goes quiet.

4. How long will I wait? (before running anything)

The average run length (ARL) is the expected number of samples before an alarm.

from dense_armor.utility.drift.cusum import two_sided_arl

print(round(two_sided_arl(mu=0.0, k=0.5, h=5.0), 1))
print(round(two_sided_arl(mu=1.0, k=0.5, h=5.0), 1))

469.1
10.3

With $k = 0.5$ and $h = 5$: when nothing has changed ($\mu = 0$) a false alarm comes on average every 469 samples; after a shift of one noise-width ($\mu = 1$) the real alarm comes after about 10 samples. These are the classic values of CUSUM tables, a useful sanity check.

The formula, for one side, with drift $\delta = \mu - k$ in noise units, is

$$\mathrm{ARL}(\delta) = \frac{e^{-2\delta h'} - 1 + 2\delta h'}{2\delta^2}, \qquad h' = h + 1.166,$$

and for $\delta = 0$ it becomes $h'^2$. This is equation (6) of Reynolds (1975), obtained by treating the running total as a Brownian motion; the $+1.166$ is Siegmund's correction for the last jump that overshoots the threshold. The two sides combine as $1/\mathrm{ARL} = 1/\mathrm{ARL}^+ + 1/\mathrm{ARL}^-$.

5. The same estimate in your sensor's units

detectability_report turns the ARL into numbers for a real deployment: the noise level of your signal and the size of the shift you care about.

from dense_armor.utility.drift.cusum import detectability_report

r = detectability_report(local_noise_scale=0.02, k=0.5, h=20.0, candidate_shift=0.04)
for key, val in r.items():
    print(key, val)

false_alarm_arl 1556958360.0943098
detection_arl 13.888444444444445
shift_in_sigma 2.0

A joint whose noise is 0.02 rad/s and a shift of 0.04 rad/s is a shift of 2 noise-widths. With the default threshold $h = 20$ the alarm comes about 14 samples after it, and a false alarm is expected only once in about 1.6 billion samples. The detector reports the same two numbers for its own settings:

from dense_armor.utility.drift.detector import CUSUMDriftDetector

det = CUSUMDriftDetector(reference="fixed")
print(round(det.expected_detection_delay(1.0), 1))
print(round(det.expected_false_alarm_run(), 1))

40.3
1556958360.1

For a 1-sigma shift it predicts 40.3 samples; measured on a stream with reference="fixed", the delay is 41.

6. A whole recording at once

cusum_detector runs on a recorded series in one call and fires at exactly the same samples as the streaming detector.

import numpy as np
from dense_armor.utility.drift.cusum import cusum_detector
from dense_armor.utility.drift.detector import CUSUMDriftDetector

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(0, 1, 500), rng.normal(1.5, 1, 500)])
flags, s = cusum_detector(x, reference="fixed")
det = CUSUMDriftDetector(reference="fixed")
print(np.flatnonzero(flags)[:4], [i for i, v in enumerate(x) if det.update(v).drift_detected][:4])

[511 524 533 547] [511, 524, 533, 547]

Use the batch function to study a recording, and the streaming detector inside the control loop: the alarms are the same.

7. Results: CUSUM against two classic drift detectors

20 seeded streams of 1,000 N(0, 1) samples per row, mean shift at sample 500. Each cell is recall (fraction of streams where the shift was found) / delay in samples / false alarms per stream:

shift (σ)	CUSUM adaptive	CUSUM fixed	Page–Hinkley	ADWIN
0.5	0.00 / — / 0.00	0.55 / 96.0 / 0.10	1.00 / 81.0 / 0.05	1.00 / 127.8 / 0.00
1.0	0.00 / — / 0.00	0.90 / 42.3 / 0.50	1.00 / 41.2 / 0.10	1.00 / 54.2 / 0.00
1.5	0.15 / 15.0 / 0.05	1.00 / 20.9 / 0.40	1.00 / 42.9 / 0.10	1.00 / 43.0 / 0.00
2.0	0.05 / 6.0 / 0.00	1.00 / 12.6 / 0.35	1.00 / 14.9 / 0.05	1.00 / 41.4 / 0.00
3.0	0.85 / 9.8 / 0.00	1.00 / 7.8 / 0.05	1.00 / 11.2 / 0.05	1.00 / 11.0 / 0.00

"fixed" finds a sustained shift fastest from 1.5σ upward, at the price of more false alarms than Page–Hinkley and ADWIN. "adaptive" is built for something else: its sliding reference catches up with a new level, so it reacts to the leading edge of a drift and misses most small step changes. ADWIN never raised a false alarm here but is the slowest on small shifts.

API reference

cusum

utility/cusum.py

Two-sided CUSUM (cumulative sum) change-point detector -- Page, E.S. (1954), "Continuous Inspection Schemes", Biometrika 41(1-2):100-115. Standard textbook change-point detection, not specific to this project; added to close a real, measured gap: test/test_benchmark_v0_runtime_ behavioral_drift.py found that utility/arbiter.py's classify_segments (instantaneous n_sigmas threshold on a windowed robust z-score) detects a SHARP transient spike perfectly (100% detection, latency 0 at every severity tested) but is much less sensitive to a SLOW, GRADUAL sustained drift (8.9%/24.4%/26.7% detected across the ramp+settling window at mild/moderate/strong severity) -- an instantaneous threshold test is structurally the wrong tool for a shift too small, at any single step, to itself cross the threshold. CUSUM is the standard classical answer: it accumulates small, sustained deviations over time instead of judging each point in isolation, so a drift too gradual to trip an instantaneous threshold still eventually trips the accumulated sum.

HONEST DEVIATION FROM THE CITED ALGORITHM: Page's original CUSUM (and the textbook treatment since) accumulates deviation against a FIXED reference mean, decided once from an in-control baseline and never updated. This implementation instead recomputes med/scale from a SLIDING causal window every step (same convention as arbiter.py), so the "reference" itself drifts along with the data over time. This is a deliberate choice (matching classify_segments' own adaptive-window convention, and needed for a streaming setting where there is no fixed known-good baseline to pin against) but it is NOT literally Page's fixed-reference scheme -- call this an adaptive-reference CUSUM, not an unqualified claim to be running the textbook algorithm unmodified. The practical effect, measured in test/test_benchmark_v0_runtime_behavioral_ drift.py: it detects the LEADING EDGE of a sustained drift quickly (3-8 points into a 15-point ramp), then stops accumulating once its own reference window has caught up to the new level -- it will NOT keep accumulating forever against a genuinely fixed target the way Page's original scheme does.

This is a NEW, SEPARATE detector, not a replacement for classify_segments and not a retuning of its shipped thresholds (radius=10, ref_mult=3, n_sigmas=3.0, spike_run_max=2 in arbiter.py are untouched by this file, and remain exactly the preregistered defaults the benchmark declared it would never adjust after seeing results) -- classify_segments stays the right tool for transient spikes/glitches (see the benchmark's condition 3), this is the complementary tool for condition 2's slow-drift blind spot. A production system would run both channels and flag on either.

Reference window and robust standardization deliberately mirror utility/arbiter.py's own convention (causal window, radiusref_mult span, median/1.4826MAD center-scale) so a point's "z" here is in the same robust-sigma units classify_segments already reports, rather than introducing a second, incompatible notion of scale.

cusum_detector

cusum_detector(
    x: ndarray,
    radius: int = 10,
    ref_mult: int = 3,
    k: float = 0.5,
    h: float = 20.0,
    eps: float = 1e-09,
    reference: str = "adaptive",
) -> Tuple[np.ndarray, np.ndarray]

Two-sided CUSUM on the same causal robust-z-score classify_segments uses, so a sustained drift too gradual to cross classify_segments' OWN instantaneous n_sigmas threshold still accumulates and eventually trips this one.

S+[i] = max(0, S+[i-1] + (z[i] - k)) -- accumulates upward drift S-[i] = min(0, S-[i-1] + (z[i] + k)) -- accumulates downward drift flags when |S+[i]| > h or |S-[i]| > h; both sums reset to 0 the instant they flag, so a second, later shift can be detected too (standard CUSUM reset behavior, not specific to this implementation).

k (the "slack"/reference value, in robust-sigma units) is the classical convention: half the shift size you want to detect -- Page's original scheme and the textbook tuning since (e.g. Hawkins & Olwell, "Cumulative Sum Charts and Charting for Quality Improvement", 1998). h (the decision threshold) is NOT the commonly-cited h~4-5: checked directly (2026-09-03, prompted by a parallel investigation porting this same algorithm to a streaming detector) before shipping that value -- h=5.0 gives a 100% ("adaptive") / 85% ("fixed") stream-level false-alarm rate on a 1000-sample purely stable N(0,1) series (200 trials, radius=10, ref_mult=3, k=0.5), because its average run length under no-change is only ~19-38 samples, nowhere near 1000. h=20.0 (this default) empirically gives 3.5% ("adaptive") / 15.5% ("fixed") over the same 1000-sample horizon. "fixed" stays structurally more exposed than "adaptive" at any shared h -- its target never updates, so an unlucky warmup-window estimate is never self-corrected the way "adaptive"'s sliding reference recovers from one; not a tuning gap, a direct consequence of what "fixed" is for (see its own docstring below). "Textbook k=0.5/h=5.0" is a real, commonly-cited tuning for short/industrial monitoring horizons; it was not previously re-validated here against the longer streams typical of software/ML drift monitoring. test/test_cusum.py's existing assertions (large, unambiguous 10-sigma steps; a <10% POINT-level noise flag rate, which h=5.0 already satisfied despite its much higher STREAM-level false-alarm rate -- the same per-point-vs-per-stream distinction that motivated this fix) all still hold at h=20.0, re-verified directly, not assumed.

Parameters

x : np.ndarray 1D signal. radius, ref_mult : int Causal reference window span = radiusref_mult, same convention as protect.arbiter.classify_segments. k : float CUSUM slack, in robust-sigma units. h : float Decision threshold, in accumulated robust-sigma units. eps : float Degenerate-scale guard, same convention as classify_segments. reference : {"adaptive", "fixed"}, default "adaptive" "adaptive" (this module's own default, see the module docstring for the honest disclosure of how this differs from Page's original scheme): med/scale recomputed every step from a sliding causal window, so the reference itself drifts along with the data -- catches the LEADING EDGE of a sustained drift, then stops accumulating once the window has caught up to the new level. "fixed": med/scale computed ONCE from x[0:span] (the first radiusref_mult points, assumed in-control/baseline) and never updated afterward -- this IS Page's original scheme. Points before index span have no valid fixed reference yet and are skipped (flagged=False, cusum=0), same warmup convention as "adaptive"'s own causal window. Unlike "adaptive", a genuinely sustained shift will keep accumulating against this fixed target indefinitely if never reset by a flag, rather than fading once a sliding window adapts to the new level.

Returns

(flagged, cusum) flagged : bool array, True where either accumulator crossed h. cusum : float array, signed accumulated statistic (S+ where positive-going, S- where negative-going) at each point, for inspection/plotting.

Source code in dense_armor/utility/drift/cusum.py

def cusum_detector(
    x: np.ndarray, radius: int = 10, ref_mult: int = 3,
    k: float = 0.5, h: float = 20.0, eps: float = 1e-9,
    reference: str = "adaptive",
) -> Tuple[np.ndarray, np.ndarray]:
    """Two-sided CUSUM on the same causal robust-z-score classify_segments
    uses, so a sustained drift too gradual to cross classify_segments'
    OWN instantaneous n_sigmas threshold still accumulates and eventually
    trips this one.

    S+[i] = max(0, S+[i-1] + (z[i] - k))   -- accumulates upward drift
    S-[i] = min(0, S-[i-1] + (z[i] + k))   -- accumulates downward drift
    flags when |S+[i]| > h or |S-[i]| > h; both sums reset to 0 the
    instant they flag, so a second, later shift can be detected too
    (standard CUSUM reset behavior, not specific to this implementation).

    k (the "slack"/reference value, in robust-sigma units) is the
    classical convention: half the shift size you want to detect --
    Page's original scheme and the textbook tuning since (e.g. Hawkins &
    Olwell, "Cumulative Sum Charts and Charting for Quality Improvement",
    1998). h (the decision threshold) is NOT the commonly-cited h~4-5:
    checked directly (2026-09-03, prompted by a parallel investigation
    porting this same algorithm to a streaming detector) before shipping that
    value -- h=5.0 gives a 100% ("adaptive") / 85% ("fixed") stream-level
    false-alarm rate on a 1000-sample purely stable N(0,1) series (200
    trials, radius=10, ref_mult=3, k=0.5), because its average run length
    under no-change is only ~19-38 samples, nowhere near 1000. h=20.0
    (this default) empirically gives 3.5% ("adaptive") / 15.5% ("fixed")
    over the same 1000-sample horizon. "fixed" stays structurally more
    exposed than "adaptive" at any shared h -- its target never updates,
    so an unlucky warmup-window estimate is never self-corrected the way
    "adaptive"'s sliding reference recovers from one; not a tuning gap,
    a direct consequence of what "fixed" is for (see its own docstring
    below). "Textbook k=0.5/h=5.0" is a real, commonly-cited tuning for
    short/industrial monitoring horizons; it was not previously
    re-validated here against the longer streams typical of software/ML
    drift monitoring. test/test_cusum.py's existing assertions (large,
    unambiguous 10-sigma steps; a <10% POINT-level noise flag rate, which
    h=5.0 already satisfied despite its much higher STREAM-level
    false-alarm rate -- the same per-point-vs-per-stream distinction that
    motivated this fix) all still hold at h=20.0, re-verified directly,
    not assumed.

    Parameters
    ----------
    x : np.ndarray
        1D signal.
    radius, ref_mult : int
        Causal reference window span = radius*ref_mult, same convention
        as protect.arbiter.classify_segments.
    k : float
        CUSUM slack, in robust-sigma units.
    h : float
        Decision threshold, in accumulated robust-sigma units.
    eps : float
        Degenerate-scale guard, same convention as classify_segments.
    reference : {"adaptive", "fixed"}, default "adaptive"
        "adaptive" (this module's own default, see the module docstring
        for the honest disclosure of how this differs from Page's
        original scheme): med/scale recomputed every step from a sliding
        causal window, so the reference itself drifts along with the
        data -- catches the LEADING EDGE of a sustained drift, then
        stops accumulating once the window has caught up to the new
        level.
        "fixed": med/scale computed ONCE from x[0:span] (the first
        radius*ref_mult points, assumed in-control/baseline) and never
        updated afterward -- this IS Page's original scheme. Points
        before index `span` have no valid fixed reference yet and are
        skipped (flagged=False, cusum=0), same warmup convention as
        "adaptive"'s own causal window. Unlike "adaptive", a genuinely
        sustained shift will keep accumulating against this fixed target
        indefinitely if never reset by a flag, rather than fading once a
        sliding window adapts to the new level.

    Returns
    -------
    (flagged, cusum)
        flagged : bool array, True where either accumulator crossed h.
        cusum   : float array, signed accumulated statistic (S+ where
                  positive-going, S- where negative-going) at each point,
                  for inspection/plotting.
    """
    if reference not in ("adaptive", "fixed"):
        raise ValueError(f"reference must be 'adaptive' or 'fixed', got {reference!r}")
    x = np.asarray(x, dtype=np.float64).ravel()
    n = x.size
    flagged = np.zeros(n, dtype=bool)
    cusum = np.zeros(n, dtype=np.float64)
    s_pos, s_neg = 0.0, 0.0
    span = radius * ref_mult

    fixed_med, fixed_scale = None, None
    if reference == "fixed" and n > span:
        init_window = x[:span]
        if np.all(np.isfinite(init_window)):
            fixed_med, fixed_scale = _robust_center_scale(init_window)

    for i in range(n):
        if reference == "fixed":
            if i < span or fixed_med is None or not np.isfinite(fixed_scale) or fixed_scale < eps:
                continue
            med, scale = fixed_med, fixed_scale
            if not np.isfinite(x[i]):
                continue
            z = (x[i] - med) / scale
            s_pos = max(0.0, s_pos + z - k)
            s_neg = min(0.0, s_neg + z + k)
            cusum[i] = s_pos if s_pos >= abs(s_neg) else s_neg
            if s_pos > h:
                flagged[i] = True
                s_pos = 0.0
            elif s_neg < -h:
                flagged[i] = True
                s_neg = 0.0
            continue

        if not np.isfinite(x[i]):
            # Explicit, not incidental: without this guard, a NaN x[i] would
            # flow into z (nan) and then into max(0.0, s_pos+z-k)/min(0.0,
            # s_neg+z+k) -- Python's max/min happen to return the FIRST
            # argument when compared against NaN (max(0.0, nan) == 0.0,
            # min(0.0, nan) == 0.0, verified directly), which silently
            # clamps the accumulator to 0 only because `0.0` is written
            # first in this file's own calls below -- swap the argument
            # order in a future edit and NaN would propagate instead.
            # Holding state (skip, don't touch s_pos/s_neg) is the explicit,
            # order-independent version of that same "don't corrupt the
            # accumulator" intent.
            continue
        w_ref = _window_causal(x, i, span)
        if w_ref.size < 4:
            continue
        med, scale = _robust_center_scale(w_ref)
        if not np.isfinite(scale) or scale < eps:
            # not np.isfinite(scale): a NaN/Inf point earlier still inside
            # this causal window (span points back) would otherwise make
            # med/scale themselves NaN -- same explicit-not-incidental
            # reasoning as above, this time for a contaminated window
            # rather than the current point.
            continue
        z = (x[i] - med) / scale
        s_pos = max(0.0, s_pos + z - k)
        s_neg = min(0.0, s_neg + z + k)
        cusum[i] = s_pos if s_pos >= abs(s_neg) else s_neg
        if s_pos > h:
            flagged[i] = True
            s_pos = 0.0
        elif s_neg < -h:
            flagged[i] = True
            s_neg = 0.0

    return flagged, cusum

one_sided_arl

one_sided_arl(
    delta: float, h: float, corrected: bool = True
) -> float

Wald/Brownian-motion approximation to the Average Run Length (ARL) of a one-sided CUSUM with standardized drift delta (true mean minus reference value k) and decision boundary h.

THEORY: Page, E.S. (1954), "Continuous Inspection Schemes", Biometrika 41(1-2):100-115 (this module's own citation) introduced CUSUM. The closed-form Brownian-motion/Wald-type ARL approximation used here follows Reynolds, M.R. Jr. (1975), "Approximations to the Average Run Length in Cumulative Sum Control Charts", Technometrics 17(1):65-71, with the boundary/continuity correction refined by Siegmund, D. (1985), Sequential Analysis: Tests and Confidence Intervals, Theorem 10.16, to ~1.166 (Reynolds' own empirical value was ~1.2). Verified directly against Monte Carlo simulation of the exact idealized process before use (see Dense-Evolution-Discovery's scripts/cusum_detectability_theory/validate_arl_theory.py): the corrected formula matches to ~1%, the uncorrected one underestimates by ~20% at realistic h -- Reynolds' own 1975 finding, reproduced.

Parameters

delta : float Standardized drift (true mean - k). delta=0 is the in-control (null) case for that branch. h : float CUSUM decision boundary, standardized (robust-sigma) units, same as this module's own h parameter. corrected : bool, default True Apply Siegmund's boundary correction (h -> h + 1.166) -- leave True unless deliberately reproducing the plain, uncorrected textbook Wald formula.

Returns

float Approximate average run length, in samples.

Source code in dense_armor/utility/drift/cusum.py

def one_sided_arl(delta: float, h: float, corrected: bool = True) -> float:
    """Wald/Brownian-motion approximation to the Average Run Length (ARL)
    of a one-sided CUSUM with standardized drift `delta` (true mean minus
    reference value k) and decision boundary `h`.

    THEORY: Page, E.S. (1954), "Continuous Inspection Schemes", Biometrika
    41(1-2):100-115 (this module's own citation) introduced CUSUM. The
    closed-form Brownian-motion/Wald-type ARL approximation used here
    follows Reynolds, M.R. Jr. (1975), "Approximations to the Average Run
    Length in Cumulative Sum Control Charts", Technometrics 17(1):65-71,
    with the boundary/continuity correction refined by Siegmund, D.
    (1985), Sequential Analysis: Tests and Confidence Intervals, Theorem
    10.16, to ~1.166 (Reynolds' own empirical value was ~1.2). Verified
    directly against Monte Carlo simulation of the exact idealized
    process before use (see Dense-Evolution-Discovery's
    scripts/cusum_detectability_theory/validate_arl_theory.py): the
    corrected formula matches to ~1%, the uncorrected one underestimates
    by ~20% at realistic h -- Reynolds' own 1975 finding, reproduced.

    Parameters
    ----------
    delta : float
        Standardized drift (true mean - k). delta=0 is the in-control
        (null) case for that branch.
    h : float
        CUSUM decision boundary, standardized (robust-sigma) units, same
        as this module's own `h` parameter.
    corrected : bool, default True
        Apply Siegmund's boundary correction (h -> h + 1.166) -- leave
        True unless deliberately reproducing the plain, uncorrected
        textbook Wald formula.

    Returns
    -------
    float
        Approximate average run length, in samples.
    """
    hh = h + BOUNDARY_CORRECTION if corrected else h
    if abs(delta) < 1e-12:
        return hh ** 2
    z = -2 * delta * hh
    if z > 700.0:
        return float("inf")
    return (np.exp(z) - 1 + 2 * delta * hh) / (2 * delta ** 2)

two_sided_arl

two_sided_arl(
    mu: float, k: float, h: float, corrected: bool = True
) -> float

ARL of this module's symmetric two-sided CUSUM (same k, h on both the S+/S- accumulators, matching cusum_detector's own convention) under a true standardized mean shift mu. Combines the two one-sided branches via 1/ARL = 1/ARL+ + 1/ARL- (standard two-sided-from- one-sided combination rule, e.g. Reynolds 1975 eq. 13).

Source code in dense_armor/utility/drift/cusum.py

def two_sided_arl(mu: float, k: float, h: float, corrected: bool = True) -> float:
    """ARL of this module's symmetric two-sided CUSUM (same k, h on both
    the S+/S- accumulators, matching `cusum_detector`'s own convention)
    under a true standardized mean shift `mu`. Combines the two one-sided
    branches via 1/ARL = 1/ARL+ + 1/ARL- (standard two-sided-from-
    one-sided combination rule, e.g. Reynolds 1975 eq. 13)."""
    arl_pos = one_sided_arl(mu - k, h, corrected=corrected)
    arl_neg = one_sided_arl(-mu - k, h, corrected=corrected)
    return 1.0 / (1.0 / arl_pos + 1.0 / arl_neg)

detectability_report

detectability_report(
    local_noise_scale: float,
    k: float,
    h: float,
    candidate_shift: float,
) -> dict

Convert the standardized-units ARL theory above into REAL units for a specific cusum_detector deployment -- a pre-flight estimate of how long detection/false-alarms should take, computable BEFORE running a benchmark, from a detector's real local noise level and a candidate shift size.

VALIDATED, honestly, on TWO independent real physical domains (not just synthetic Monte Carlo), promoted from Dense-Evolution-Discovery after both checks (see scripts/cusum_detectability_theory/ there): - Real lidar (Sydney Urban Objects, 7 real independent points): in all 7 cases, real detection latency was LOWER than the predicted mean ARL -- a consistent one-directional bias (real threshold crossings happen faster than the idealized iid-Gaussian model predicts), not scatter. - Real accelerometer (UCI HAR, 5 real independent points): a genuinely MIXED result at moderate SNR (2/5 points faster than predicted, 3/5 slower) -- NOT the same one-directional bias lidar showed. Documented as found, not forced to match. Also surfaced the extreme-SNR floor issue below: this domain's real local noise was quiet enough that the raw formula predicted a fractional, physically meaningless ARL.

HONEST TAKEAWAY: this is a pre-flight estimate under classical iid-Gaussian assumptions, not an oracle for real, non-Gaussian sensor data -- use it to reason about detectability before running a benchmark, then measure the real false-positive/detection rate for the real deployment, don't trust the number alone. The DIRECTION of the theory-vs-reality gap is not guaranteed to be consistent across physically different domains (confirmed: lidar biased one way, accelerometer gave a mixed result) -- do not assume a single correction factor transfers between domains.

Parameters

local_noise_scale : float The detector's own real local noise scale (e.g. a causal window's median/MAD*1.4826, the same quantity cusum_detector computes internally) -- real units (e.g. meters, g, seconds). k, h : float The CUSUM's own k, h parameters, in STANDARDIZED (robust-sigma) units, matching cusum_detector's own convention. candidate_shift : float A real-unit shift magnitude to evaluate detectability for (e.g. "a +10m persistent offset").

Returns

dict with: false_alarm_arl : expected samples between false alarms (mu=0), floored at 1.0. detection_arl : expected samples to detect candidate_shift, floored at 1.0. shift_in_sigma : candidate_shift / local_noise_scale.

FLOOR AT 1.0: at extreme shift_in_sigma (real local noise tiny relative to candidate_shift -- the real accelerometer check above hit this at >1000 sigma), the raw Wald/Siegmund formula returns a fractional ARL below 1, which has no physical meaning (a detector cannot flag in under one real sample). Both ARLs are floored at 1.0 here; treat any raw value that needed flooring as "near-instant detection", not a precise number -- the formula is an asymptotic approximation, least trustworthy exactly in this regime.

Source code in dense_armor/utility/drift/cusum.py

def detectability_report(local_noise_scale: float, k: float, h: float, candidate_shift: float) -> dict:
    """Convert the standardized-units ARL theory above into REAL units
    for a specific `cusum_detector` deployment -- a pre-flight estimate
    of how long detection/false-alarms should take, computable BEFORE
    running a benchmark, from a detector's real local noise level and a
    candidate shift size.

    VALIDATED, honestly, on TWO independent real physical domains (not
    just synthetic Monte Carlo), promoted from Dense-Evolution-Discovery
    after both checks (see scripts/cusum_detectability_theory/ there):
    - Real lidar (Sydney Urban Objects, 7 real independent points): in
      all 7 cases, real detection latency was LOWER than the predicted
      mean ARL -- a consistent one-directional bias (real threshold
      crossings happen faster than the idealized iid-Gaussian model
      predicts), not scatter.
    - Real accelerometer (UCI HAR, 5 real independent points): a
      genuinely MIXED result at moderate SNR (2/5 points faster than
      predicted, 3/5 slower) -- NOT the same one-directional bias lidar
      showed. Documented as found, not forced to match. Also surfaced
      the extreme-SNR floor issue below: this domain's real local noise
      was quiet enough that the raw formula predicted a fractional,
      physically meaningless ARL.

    HONEST TAKEAWAY: this is a pre-flight estimate under classical
    iid-Gaussian assumptions, not an oracle for real, non-Gaussian
    sensor data -- use it to reason about detectability before running a
    benchmark, then measure the real false-positive/detection rate for
    the real deployment, don't trust the number alone. The DIRECTION of
    the theory-vs-reality gap is not guaranteed to be consistent across
    physically different domains (confirmed: lidar biased one way,
    accelerometer gave a mixed result) -- do not assume a single
    correction factor transfers between domains.

    Parameters
    ----------
    local_noise_scale : float
        The detector's own real local noise scale (e.g. a causal
        window's median/MAD*1.4826, the same quantity `cusum_detector`
        computes internally) -- real units (e.g. meters, g, seconds).
    k, h : float
        The CUSUM's own k, h parameters, in STANDARDIZED (robust-sigma)
        units, matching `cusum_detector`'s own convention.
    candidate_shift : float
        A real-unit shift magnitude to evaluate detectability for (e.g.
        "a +10m persistent offset").

    Returns
    -------
    dict with:
        false_alarm_arl : expected samples between false alarms (mu=0),
                           floored at 1.0.
        detection_arl    : expected samples to detect `candidate_shift`,
                            floored at 1.0.
        shift_in_sigma   : candidate_shift / local_noise_scale.

    FLOOR AT 1.0: at extreme shift_in_sigma (real local noise tiny
    relative to candidate_shift -- the real accelerometer check above
    hit this at >1000 sigma), the raw Wald/Siegmund formula returns a
    fractional ARL below 1, which has no physical meaning (a detector
    cannot flag in under one real sample). Both ARLs are floored at 1.0
    here; treat any raw value that needed flooring as "near-instant
    detection", not a precise number -- the formula is an asymptotic
    approximation, least trustworthy exactly in this regime.
    """
    mu_std = candidate_shift / local_noise_scale
    return dict(
        false_alarm_arl=max(1.0, two_sided_arl(0.0, k, h)),
        detection_arl=max(1.0, two_sided_arl(mu_std, k, h)),
        shift_in_sigma=mu_std,
    )

detector

CUSUM drift detector for the streaming interface.

CUSUMDriftDetector

CUSUMDriftDetector(
    radius: int = 10,
    ref_mult: int = 3,
    k: float = 0.5,
    h: float = 20.0,
    two_sided: bool = True,
    eps: float = 1e-09,
    reference: str = "adaptive",
)

Bases: DriftDetector

CUSUM drift detector matching the batch cusum_detector behaviour.

Parameters

radius, ref_mult Causal window span = radius * ref_mult, same convention as cusum_detector. k CUSUM slack, in robust-sigma units. h Decision threshold, in accumulated robust-sigma units. two_sided If True, both upward and downward shifts are monitored. If False, only upward shifts are monitored. eps Degenerate-scale guard, same convention as cusum_detector. reference "adaptive" recomputes median and scale from a sliding causal window every step, matching cusum_detector(reference="adaptive"). The reference drifts with the data, so a sustained shift is caught at its leading edge, then stops accumulating once the window catches up. The detection floor for a sustained step is roughly 2.2 robust sigmas at h=20. "fixed" computes median and scale once from the first span usable (finite, non-flat) values and never updates, matching cusum_detector(reference="fixed") and Page's original scheme. Non-finite samples do not count towards span; the fixed buffer keeps only the last span finite values, so a long flat start does not make the buffer grow without bound, and the reference is retried on every new sample until a usable one is found.

Examples

from dense_armor.utility.drift.detector import CUSUMDriftDetector import numpy as np det = CUSUMDriftDetector(reference="fixed") rng = np.random.default_rng(0) stream = list(rng.normal(0, 1, 500)) + list(rng.normal(1, 1, 500)) for i, x in enumerate(stream): ... _ = det.update(x) ... if det.drift_detected: ... print(i) ... break 516

Source code in dense_armor/utility/drift/detector.py

def __init__(self, radius: int = 10, ref_mult: int = 3, k: float = 0.5,
             h: float = 20.0, two_sided: bool = True, eps: float = 1e-9,
             reference: str = "adaptive") -> None:
    super().__init__()
    if reference not in ("adaptive", "fixed"):
        raise ValueError(
            f"reference must be 'adaptive' or 'fixed', got {reference!r}"
        )
    self.radius = radius
    self.ref_mult = ref_mult
    self.k = k
    self.h = h
    self.two_sided = two_sided
    self.eps = eps
    self.reference = reference
    self._span = radius * ref_mult
    self._buffer: deque[float] = deque(maxlen=self._span)
    self._fixed_buffer: deque[float] = deque(maxlen=self._span)
    self._s_pos = 0.0
    self._s_neg = 0.0
    self._fixed_ready = False
    self._fixed_med = 0.0
    self._fixed_scale = 1.0
    self.n_missing_ = 0

expected_false_alarm_run

expected_false_alarm_run() -> float

Expected number of samples between false alarms under no change.

Source code in dense_armor/utility/drift/detector.py

def expected_false_alarm_run(self) -> float:
    """Expected number of samples between false alarms under no change."""
    return two_sided_arl(0.0, self.k, self.h)

expected_detection_delay

expected_detection_delay(shift: float) -> float

Expected number of samples to detect a standardized mean shift shift.

Source code in dense_armor/utility/drift/detector.py

def expected_detection_delay(self, shift: float) -> float:
    """Expected number of samples to detect a standardized mean shift `shift`."""
    return two_sided_arl(shift, self.k, self.h)

Details

• Real-world check of the ARL estimate. detectability_report was validated on two real physical domains before being promoted from Dense-Evolution-Discovery. On real lidar (Sydney Urban Objects, 7 independent points) the real detection was always faster than the predicted mean ARL, a consistent bias in one direction. On a real accelerometer (UCI HAR, 5 points) the result was mixed: 2 points faster than predicted, 3 slower. The estimate assumes independent Gaussian noise; use it to reason before a benchmark, then measure the real rates.
• Very quiet sensors. When the local noise is extremely small the raw formula can predict a fractional ARL, which has no physical meaning; the docstring of detectability_report explains the floor applied.
• Relation to the Arbiter. The Arbiter judges each point against an per-point threshold; CUSUM is its complement for slow, sustained changes.
• Sources. Reynolds, M. R. (1975), "Approximations to the Average Run Length in Cumulative Sum Control Charts", Technometrics 17(1), equation (6). Siegmund, D. (1985), Sequential Analysis, for the 1.166 boundary correction. Page, E. S. (1954), "Continuous inspection schemes", Biometrika 41, for the CUSUM rule.