CUSUM Charts for Detecting Small Sustained Shifts

This how-to builds a tabular cumulative sum (CUSUM) chart in Python: the two one-sided sums that accumulate deviations from target so a small, persistent shift in the process mean crosses a decision threshold quickly, while transient noise resets to zero and never accumulates. It is the CUSUM-focused companion to EWMA and CUSUM charts for small-shift detection; where an EWMA smooths, the CUSUM sums, and it gives you two extras a smoothed chart cannot — an estimate of when the shift began and how large it is. We derive the reference value and decision interval from the shift you want to catch, then verify the chart against a fixture whose signal point you can reason about by hand.

Prerequisites

Before you accumulate a single sum, confirm the upstream state:

  • Python 3.9+ with numpy and pandas installed (pip install numpy pandas).
  • A frozen Phase I baseline: an in-control mean $\mu_0$ and standard deviation $\sigma$ from at least 20–30 stable observations with no known assignable causes.
  • Observations are approximately independent — lag-1 autocorrelation below about 0.25. Serial correlation makes the decision interval wrong, the same failure that pinches an I-MR chart.
  • A stated target shift $\delta$ (in σ units) — this sets the reference value $k$.
  • Data arrives as a numeric pandas Series in observation order, with nulls resolved upstream.

Step 1 — Fix the baseline and target shift

The CUSUM is a Phase II monitor, so $\mu_0$ and $\sigma$ are frozen inputs. Record them with the target shift $\delta$ up front; every parameter below is derived from these three numbers, so capturing them immutably keeps the chart auditable.

from dataclasses import dataclass


@dataclass(frozen=True)
class CusumDesign:
    """Frozen CUSUM design: baseline plus target shift and threshold."""
    mu0: float           # in-control mean (Phase I)
    sigma: float         # in-control std dev (Phase I)
    delta: float = 1.0   # target shift to detect, in sigma units
    h: float = 5.0       # decision interval, in sigma units

    def __post_init__(self) -> None:
        if self.sigma <= 0:
            raise ValueError("sigma must be positive.")
        if self.delta <= 0:
            raise ValueError("delta must be positive.")
        if self.h <= 0:
            raise ValueError("h must be positive.")

Step 2 — Derive the reference value k and decision interval H

The reference value $k$ (the slack) is the amount subtracted from each deviation before it accumulates, and it is set to half the shift you want to detect: to catch a $\delta\sigma$ shift, $k = \tfrac{\delta}{2}\sigma$. The decision interval $H$ is the threshold a sum must exceed to signal, conventionally $H = h\sigma$ with $h$ between 4 and 5.

$$k = \frac{\delta}{2}\,\sigma \qquad H = h\,\sigma$$

The canonical pairing $k = 0.5\sigma$ (for a $1\sigma$ target) with $H = 5\sigma$ gives an in-control run length near 370 — matching a Shewhart chart — while detecting a $1\sigma$ shift in about 10 samples. Halving $H$ toward $4\sigma$ speeds detection but shortens the in-control run length, so move it only with an ARL check in hand.

def cusum_params(design: CusumDesign) -> tuple[float, float]:
    """Return (k, H) in the data's original units from the design."""
    k = (design.delta / 2.0) * design.sigma   # slack = half the target shift
    H = design.h * design.sigma               # decision interval
    return k, H

Step 3 — Accumulate the one-sided sums

The tabular CUSUM keeps two sums. The upper sum $C_t^+$ accumulates evidence of an upward shift; the lower sum $C_t^-$ accumulates evidence of a downward shift. Each is floored at zero by a $\max(0,\cdot)$, which is what makes the chart reset after transient noise and respond only to sustained drift:

$$C_t^+ = \max\!\left(0,\ x_t - (\mu_0 + k) + C_{t-1}^+\right), \qquad C_t^- = \max\!\left(0,\ (\mu_0 - k) - x_t + C_{t-1}^-\right)$$

with both sums seeded at zero. The recursion is inherently sequential, so a Python loop over the array is the correct and clearest implementation.

import numpy as np
import pandas as pd


def cusum_sums(values: pd.Series, design: CusumDesign) -> pd.DataFrame:
    """Compute the one-sided CUSUM sums C+ and C- with a run-length counter."""
    x = np.asarray(values, dtype=np.float64)
    if x.size == 0:
        raise ValueError("Input series is empty.")
    if np.isnan(x).any():
        raise ValueError("Resolve missing values upstream before charting.")

    k, _ = cusum_params(design)
    upper_ref = design.mu0 + k
    lower_ref = design.mu0 - k

    c_plus = np.zeros(x.size)
    c_minus = np.zeros(x.size)
    n_plus = np.zeros(x.size, dtype=int)   # samples since C+ last reset
    n_minus = np.zeros(x.size, dtype=int)
    cp = cm = 0.0
    np_run = nm_run = 0
    for t in range(x.size):
        cp = max(0.0, x[t] - upper_ref + cp)
        cm = max(0.0, lower_ref - x[t] + cm)
        np_run = np_run + 1 if cp > 0 else 0
        nm_run = nm_run + 1 if cm > 0 else 0
        c_plus[t], c_minus[t] = cp, cm
        n_plus[t], n_minus[t] = np_run, nm_run

    return pd.DataFrame(
        {"x": x, "c_plus": c_plus, "c_minus": c_minus,
         "n_plus": n_plus, "n_minus": n_minus},
        index=getattr(values, "index", pd.RangeIndex(x.size)),
    )

Step 4 — Flag signals and interpret them

A signal fires the first time either sum exceeds $H$. Two diagnostics fall straight out of the table: the side that signaled tells you the shift direction, and the run-length counter $N^+$ or $N^-$ (samples since that sum last reset) estimates how many observations back the shift began.

def cusum_chart(values: pd.Series, design: CusumDesign) -> pd.DataFrame:
    """Assemble the CUSUM sums with signal flags and direction."""
    _, H = cusum_params(design)
    df = cusum_sums(values, design)
    df["signal_up"] = df["c_plus"] > H
    df["signal_down"] = df["c_minus"] > H
    df["signal"] = df["signal_up"] | df["signal_down"]
    return df


def estimate_shift_onset(df: pd.DataFrame) -> int | None:
    """Positional index where the shift likely began, from the run-length counter."""
    hits = df.index[df["signal"]]
    if len(hits) == 0:
        return None
    first = df.index.get_loc(hits[0])
    run = df["n_plus"].iloc[first] if df["signal_up"].iloc[first] else df["n_minus"].iloc[first]
    return int(first - run + 1)

Once flagged, route the event into the shared out-of-control rule detection layer so a CUSUM signal is de-duplicated against any Shewhart run rules watching the same stream, rather than raising a second independent alarm for one drift.

Step 5 (optional) — The V-mask equivalent

The tabular form above is what automated pipelines implement, but the classic graphical CUSUM uses a V-mask: a wedge laid over the plotted cumulative sum, with its lead distance $d$ and half-angle $\theta$ chosen so the mask arms are equivalent to the tabular $k$ and $H$. It signals when any past point falls outside the arms. The two are mathematically equivalent for the two-sided case; prefer the tabular version in code because its entire state is two scalars, which stream and serialize trivially, whereas the V-mask needs the full plotted history to re-evaluate on each point.

Verification

Prove the chart on a fixture whose behavior you can reason about: a dead-on series must never signal, a sustained shift must signal on the correct side within a handful of samples, and a lone spike must be absorbed by the reset.

import numpy as np
import pandas as pd


def test_cusum_flags_sustained_upward_shift():
    design = CusumDesign(mu0=100.0, sigma=2.0, delta=1.0, h=5.0)
    x = np.concatenate([np.full(20, 100.0), np.full(20, 102.0)])  # +1 sigma
    df = cusum_chart(pd.Series(x), design)

    assert not df["signal"].iloc[:20].any()          # silent in control
    assert df["signal_up"].any() and not df["signal_down"].any()
    first = int(df["signal"].values.argmax())
    assert 20 <= first <= 32
    assert estimate_shift_onset(df) >= 19            # points near true onset


def test_single_spike_does_not_accumulate():
    design = CusumDesign(mu0=100.0, sigma=2.0)
    x = np.full(30, 100.0)
    x[10] = 108.0                                     # one transient outlier
    df = cusum_chart(pd.Series(x), design)
    assert not df["signal"].any()                     # reset absorbs the spike


def test_deadon_series_is_silent():
    design = CusumDesign(mu0=0.0, sigma=1.0)
    df = cusum_chart(pd.Series(np.zeros(50)), design)
    assert (df["c_plus"] == 0).all() and (df["c_minus"] == 0).all()

The spike test is the one that catches sign errors: because each deviation is offset by $k$ and floored at zero, a single point 4σ above target adds one term and then decays, never reaching $H = 5\sigma$. If that test signals, the sign or the reference offset in the recursion is wrong.

Root-Cause Table

Symptom Cause Fix
Slow to signal, or misses the target shift $k$ not matched to the shift, or $H$ too large Set $k = \tfrac{\delta}{2}\sigma$; use $H = 4$–$5\sigma$
False alarms on stable data $H$ too small, or observations autocorrelated Raise $H$ toward $5\sigma$; check lag-1 ACF and model the series
A lone outlier raises an alarm Sign error or missing $k$ offset in the recursion Verify the $\max(0,\cdot)$ reset and the $\mu_0 \pm k$ references
Both sums grow together on stable data Reference offset omitted (accumulating raw deviations) Subtract $k$ (upper) and add $k$ (lower) before flooring at zero
Onset estimate is nonsense Run-length counter not reset when a sum returns to zero Reset $N^+$/$N^-$ to zero whenever the corresponding sum is zero
Sums never signal after a real shift Baseline $\mu_0$ recomputed on shifted live data Freeze $\mu_0,\sigma$ in Phase I; recalibrate only after a verified change

For chart selection across every variable and attribute chart, see SPC Fundamentals & Control Chart Taxonomy.