"""Sequential change detection: CUSUM, GLR, and their run lengths. A fixed-window detector asks "is the signal in this block?". A sequential change detector asks a different question, one sample at a time and without end: "has the process I am watching stopped being the process it was?". There is no window to average over and no final answer, only two competing costs -- how long the detector takes to notice a real change, and how often it cries wolf when nothing has happened. Parameterisation, stated up front per the arc's standing rule ------------------------------------------------------------- Everything here works on the **standardized** observation z[n] = (y[n] - mu0) / sigma so a change of size *delta* means a shift of ``delta`` standard deviations in the mean of ``z``. ``delta`` is therefore a per-sample amplitude SNR, NOT an energy-to-noise ratio like the ``enr`` of ``topics/detection-theory/detection.py`` (for an N-sample block those differ by a factor N), and NOT a power ratio. A CUSUM designed for "delta = 1" is designed for a one-sigma shift. The two design numbers follow the same convention: k reference value ("slack"), in sigma units. k = delta/2 makes the CUSUM the exact sequential likelihood-ratio test for a shift of size delta. h decision threshold, in the same units as the accumulated statistic sum(z - k). Both are dimensionless. A threshold quoted in log-likelihood units is ``delta * h``, a factor this module never applies silently: see ``llr_increment`` for the conversion and the module's tests for the pin. Run length, not false-alarm probability --------------------------------------- A sequential detector run forever alarms with probability one, so there is no Pfa to quote. The performance pair is instead ARL0 mean samples to a false alarm when nothing changes ARL1 mean samples to detection after a change of size delta (the detection *delay*, measured from the change point) and the run length itself is close to geometric under H0, so ARL0 is a mean and not a guarantee: see ``run_length_quantile``. Three routes to the ARL are provided deliberately, because a run-length formula quoted from memory is exactly as treacherous as a Cramer-Rao bound quoted from memory: arl_siegmund closed form, a diffusion approximation with an overshoot correction; cheap, and accurate to ~1% only while k is small arl_markov Brook & Evans quadrature of the integral equation; near-exact, with Richardson extrapolation on by default arl_mc direct simulation of the recursion itself The tests check all three against each other, and the topic page plots the region where the closed form stops being trustworthy. References ---------- Page (1954) introduced the CUSUM; Lorden (1971) and Moustakides (1986) established its minimax optimality; Siegmund (1985) is the source of the ARL approximation used here; Brook & Evans (1972) the Markov-chain route; Basseville & Nikiforov (1993) the standard book-length treatment, including the windowed GLR of the second half of this module. """ import numpy as np from scipy.stats import chi2, norm # ---------------------------------------------------------------------- # Standardization and the likelihood-ratio increment # ---------------------------------------------------------------------- def standardize(y, mu0: float, sigma: float) -> np.ndarray: """Map observations onto the z = (y - mu0)/sigma scale this module uses. Parameters ---------- y : array_like Observations. mu0 : float In-control (pre-change) mean. sigma : float In-control standard deviation, > 0. Returns ------- ndarray Standardized observations. """ y = np.asarray(y, dtype=float) if sigma <= 0: raise ValueError("sigma must be positive") return (y - mu0) / sigma def llr_increment(z, delta: float) -> np.ndarray: """Log-likelihood-ratio increment for a mean shift of *delta* sigma. For standardized data z ~ N(0, 1) under H0 and N(delta, 1) under H1, s[n] = ln f1(z[n]) / f0(z[n]) = delta * (z[n] - delta/2) which is ``delta`` times the CUSUM increment ``z - k`` at the matched reference value ``k = delta/2``. That factor is the whole difference between a statistic in log-likelihood units and one in sigma units, and it is the easiest way to ship a threshold that is wrong by a constant. Parameters ---------- z : array_like Standardized observations. delta : float Shift size in sigma units. Returns ------- ndarray Log-likelihood-ratio increments (nats). """ z = np.asarray(z, dtype=float) return delta * (z - delta / 2.0) # ---------------------------------------------------------------------- # CUSUM # ---------------------------------------------------------------------- def cusum_stat(z, k: float, two_sided: bool = False): """Page's CUSUM statistic, without alarms or resets. The upper recursion is S[n] = max(0, S[n-1] + z[n] - k), S[-1] = 0 which is *exactly* the running maximum over all candidate change times of the accumulated log-likelihood ratio, computed in one add and one compare per sample. That identity, not the accumulation, is what makes the CUSUM cheap: the optimization over "when did it change" has a closed-form recursion. Parameters ---------- z : array_like Standardized observations. k : float Reference value in sigma units, >= 0. Use ``delta/2`` to match a shift of size delta. two_sided : bool, default False Also return the mirrored recursion that watches for a downward shift. Returns ------- ndarray or (ndarray, ndarray) The upper statistic, or ``(upper, lower)`` when *two_sided*. The lower statistic is returned as a non-negative quantity, so both are compared against the same positive threshold. """ z = np.asarray(z, dtype=float) if k < 0: raise ValueError("k must be non-negative") n = len(z) hi = np.empty(n) s = 0.0 for i in range(n): s = max(0.0, s + z[i] - k) hi[i] = s if not two_sided: return hi lo = np.empty(n) s = 0.0 for i in range(n): s = max(0.0, s - z[i] - k) lo[i] = s return hi, lo def cusum_detect(z, k: float, h: float, two_sided: bool = False, reset: bool = True): """Run the CUSUM to alarms. Parameters ---------- z : array_like Standardized observations. k, h : float Reference value and decision threshold, both in sigma units. two_sided : bool, default False Watch for shifts in both directions. Note that this roughly halves ARL0 at a fixed h, because two arms each get a chance to fire; ``arl_*`` functions describe ONE arm. reset : bool, default True Zero the statistic after each alarm, so a long record yields a sequence of alarms rather than one. Returns ------- alarms : ndarray of int Sample indices at which the statistic crossed h. stat : ndarray The statistic actually run (the pointwise maximum of the two arms when *two_sided*), including any resets. """ z = np.asarray(z, dtype=float) if h <= 0: raise ValueError("h must be positive") n = len(z) stat = np.empty(n) alarms = [] hi = lo = 0.0 for i in range(n): hi = max(0.0, hi + z[i] - k) if two_sided: lo = max(0.0, lo - z[i] - k) cur = max(hi, lo) else: cur = hi stat[i] = cur if cur >= h: alarms.append(i) if reset: hi = lo = 0.0 return np.asarray(alarms, dtype=int), stat # ---------------------------------------------------------------------- # CUSUM run length: three independent routes # ---------------------------------------------------------------------- def arl_siegmund(k: float, h: float, delta: float = 0.0) -> float: """Siegmund's ARL approximation for the one-sided CUSUM. With increments ``z - k`` of mean ``D = delta - k`` and unit variance, treating the statistic as Brownian motion with a reflecting barrier at 0 and an absorbing one at h gives ARL = ( exp(-2 D b) + 2 D b - 1 ) / (2 D^2), b = h + 1.166 The constant 1.166 corrects for the *overshoot*: a discrete-time random walk does not stop exactly at the threshold, it jumps past it, which is worth a little more than one extra unit of threshold. Dropping the correction (b = h) is not a small sin: at thresholds designed for ARL0 = 1000 it underestimates the run length by 46% at k = 0.25, 69% at k = 0.5 and 90% at k = 1.0 (measured against ``arl_markov``; the topic page plots it). Set ``delta = 0`` for ARL0 (false alarms) and ``delta = 2k`` for ARL1 at the shift the detector was designed for. This is a *diffusion* approximation: it is excellent while the per-sample increments are small compared with the threshold, i.e. while k is small, and it loses accuracy as k grows (measured at thresholds designed for ARL0 = 1000: 0.02% error at k = 0.1, 0.1% at k = 0.25, 0.8% at k = 0.5, 5.9% at k = 1.0, 23% at k = 1.5). Small k is the regime that matters -- a CUSUM is for changes too small to see in one sample -- but the limit is real, and ``arl_markov`` costs milliseconds. Parameters ---------- k : float Reference value, > 0. h : float Decision threshold, > 0. delta : float, default 0.0 True mean shift, in sigma units. Returns ------- float Mean run length in samples. """ if k <= 0 or h <= 0: raise ValueError("k and h must be positive") d = delta - k b = h + 1.166 if abs(d) < 1e-9: # D -> 0 limit of the expression return float(b * b) return float((np.exp(-2 * d * b) + 2 * d * b - 1) / (2 * d * d)) def arl_markov(k: float, h: float, delta: float = 0.0, m: int = 200, richardson: bool = True) -> float: """Brook & Evans Markov-chain ARL for the one-sided CUSUM. The statistic lives on [0, h); discretize it into *m* cells, build the transition matrix of the resulting Markov chain, and solve ``(I - P) N = 1`` for the mean number of steps to absorption. The zero-start ARL is then one step from the exact value 0 followed by N. Convergence of the raw quadrature is O(1/m), which is slow enough to matter: at k = 0.5, h = 5 the raw values are 925.0 (m = 400) and 927.9 (m = 800) against a limit of 930.9, so even 800 cells are still 0.3% low. With *richardson* the routine evaluates at m and 2m and extrapolates ``2 A(2m) - A(m)``, which converges as O(1/m^2) for about three times the cost: at the default m = 200 that is 0.003% on a typical design (k = 0.5, h = 5) and 0.12% on the widest threshold the tests cover (k = 0.1, h = 20, where the cells are coarsest relative to the unit-variance increments). Parameters ---------- k, h : float Reference value and threshold, both > 0. delta : float, default 0.0 True mean shift in sigma units. m : int, default 200 Number of cells. With *richardson* on, 200 is already within 0.1% of the m = 400 answer on every case the tests cover. richardson : bool, default True Extrapolate from m and 2m cells. Returns ------- float Mean run length in samples. """ if k <= 0 or h <= 0: raise ValueError("k and h must be positive") if m < 8: raise ValueError("m must be at least 8") def _raw(mm): w = h / mm edges = np.arange(mm + 1) * w centres = (np.arange(mm) + 0.5) * w # From value c the next value is max(0, c + z - k), z ~ N(delta, 1), # so before the reflection it is Normal(c - k + delta, 1). loc = centres - k + delta p = (norm.cdf(edges[1:][None, :] - loc[:, None]) - norm.cdf(edges[:-1][None, :] - loc[:, None])) # Cell 0 also absorbs everything the max(0, .) reflects to zero. p[:, 0] = norm.cdf(edges[1] - loc) n_steps = np.linalg.solve(np.eye(mm) - p, np.ones(mm)) loc0 = delta - k p0 = (norm.cdf(edges[1:] - loc0) - norm.cdf(edges[:-1] - loc0)) p0[0] = norm.cdf(edges[1] - loc0) return 1.0 + float(p0 @ n_steps) if not richardson: return _raw(m) return 2.0 * _raw(2 * m) - _raw(m) def arl_mc(k: float, h: float, delta: float = 0.0, trials: int = 2000, rng=None, max_samples: int = 2_000_000) -> tuple: """Direct simulation of the CUSUM run length. Ground truth for the two approximations above: it runs the actual recursion on actual pseudorandom noise until it alarms. Parameters ---------- k, h : float Reference value and threshold. delta : float, default 0.0 True mean shift in sigma units. trials : int, default 2000 Independent runs. rng : numpy.random.Generator, optional max_samples : int Safety stop; runs still alive are counted at this length, which biases the estimate DOWNWARD, so keep it well above the ARL. Returns ------- (float, float) Mean run length and its standard error. """ rng = np.random.default_rng() if rng is None else rng s = np.zeros(trials) length = np.zeros(trials) alive = np.ones(trials, dtype=bool) n = 0 while alive.any() and n < max_samples: block = min(4096, max_samples - n) z = rng.standard_normal((block, int(alive.sum()))) + delta sa = s[alive] la = length[alive] still = np.ones(sa.shape, dtype=bool) for j in range(block): sa = np.maximum(0.0, sa + z[j] - k) la = np.where(still, la + 1, la) still &= sa < h sa = np.where(still, sa, 0.0) if not still.any(): break s[alive] = sa length[alive] = la idx = np.where(alive)[0] alive[idx[~still]] = False n += block return float(length.mean()), float(length.std(ddof=1) / np.sqrt(trials)) def run_length_quantile(arl0: float, q: float) -> float: """Quantile of the run length, treating it as geometric. The CUSUM's H0 run length is very close to geometric (memoryless): measured at k = 0.5, h = 5 over 40,000 runs, the standard deviation equals the mean to within 1% and the empirical 10th, 50th, 90th and 99th percentiles match ``-ARL0 ln(1-q)`` to about 1%. So an ARL0 of 1000 samples does NOT mean "safe for 1000 samples": it means a 10% chance of a false alarm within the first 105. Parameters ---------- arl0 : float Mean run length under H0. q : float Probability level in (0, 1). Returns ------- float The q-quantile of the run length, in samples. """ if not 0.0 < q < 1.0: raise ValueError("q must be in (0, 1)") return float(-arl0 * np.log1p(-q)) def design_cusum(delta: float, arl0: float, m: int = 200) -> tuple: """Pick (k, h) for a target shift and a target false-alarm rate. ``k = delta/2`` is the likelihood-ratio choice; h is then solved from ``arl_markov`` by bisection, not from the closed form, so the design is right even where Siegmund's approximation is not. Parameters ---------- delta : float Shift size to be detected quickly, in sigma units, > 0. arl0 : float Target mean samples between false alarms, > 1. m : int, default 200 Cell count handed to ``arl_markov``. Returns ------- (float, float) ``(k, h)``. """ if delta <= 0: raise ValueError("delta must be positive") if arl0 <= 1: raise ValueError("arl0 must exceed 1") k = delta / 2.0 lo, hi = 1e-3, 1.0 while arl_markov(k, hi, 0.0, m) < arl0: hi *= 2.0 if hi > 1e4: raise ValueError("no threshold reaches that ARL0") for _ in range(40): mid = 0.5 * (lo + hi) if arl_markov(k, mid, 0.0, m) < arl0: lo = mid else: hi = mid return k, 0.5 * (lo + hi) # ---------------------------------------------------------------------- # Windowed GLR: the post-change level is unknown # ---------------------------------------------------------------------- def glr_stat(z, window: int) -> np.ndarray: """Windowed generalised-likelihood-ratio change statistic. The CUSUM needs to be told how big a change to expect. When the post-change mean is unknown, estimate it by maximum likelihood for every candidate change time and keep the best: G[n] = max_{1<=j<=W} ( sum_{i=n-j+1..n} z[i] )^2 / (2 j) The inner quantity is the log-likelihood ratio with the ML estimate ``mu_hat = (1/j) sum z`` plugged in, so G is a two-sided statistic: it responds to shifts of either sign, at any magnitude, with no design shift to get wrong. The price is arithmetic -- O(W) per sample against the CUSUM's O(1) -- and a threshold with no usable closed form (see ``glr_threshold_mc``). Parameters ---------- z : array_like Standardized observations. window : int Maximum candidate change age W, >= 1. Changes older than W samples are no longer searched for, which bounds both the cost and the memory. Returns ------- ndarray G[n], same length as z. The first W-1 entries search only the candidates that exist yet, so they are a genuine warm-up transient, not padding. """ z = np.asarray(z, dtype=float) w = int(window) if w < 1: raise ValueError("window must be at least 1") n = len(z) if n == 0: return np.zeros(0) cs = np.concatenate(([0.0], np.cumsum(z))) pad = np.concatenate((np.full(w - 1, np.nan), cs)) view = np.lib.stride_tricks.sliding_window_view(pad, w)[:n] # view[i] = [cs[i-w+1], ..., cs[i]] -> j = w, w-1, ..., 1 j = np.arange(w, 0, -1, dtype=float) sums = cs[1:n + 1, None] - view with np.errstate(invalid="ignore"): g = np.nanmax(sums * sums / (2.0 * j), axis=1) return g def glr_detect(z, window: int, h: float, reset: bool = True): """Run the windowed GLR to alarms, restarting after each one. Parameters ---------- z : array_like Standardized observations. window : int Candidate-age window W. h : float Threshold on G, in nats. reset : bool, default True Discard all history after an alarm, which is what a monitoring system does when it has acted on the change. Returns ------- alarms : ndarray of int Indices into z at which G crossed h. stat : ndarray The statistic as run, including restarts. """ z = np.asarray(z, dtype=float) if h <= 0: raise ValueError("h must be positive") n = len(z) stat = np.zeros(n) alarms = [] start = 0 while start < n: g = glr_stat(z[start:], window) hit = np.nonzero(g >= h)[0] if len(hit) == 0 or not reset: stat[start:] = g if len(hit): alarms.extend((hit + start).tolist()) break i = int(hit[0]) stat[start:start + i + 1] = g[:i + 1] alarms.append(start + i) start += i + 1 return np.asarray(alarms, dtype=int), stat def glr_threshold_naive(p_sample: float) -> float: """The threshold a single-candidate calculation suggests -- and the trap. For ONE candidate change time the statistic ``(sum z)^2 / (2j)`` is exactly ``chi2_1 / 2`` under H0, whatever j, so P(G > h) = 2 Q(sqrt(2h)) = chi2_1.sf(2h) inviting the conclusion that ``h = chi2_1.isf(p)/2`` buys a false alarm every ``1/p`` samples. It does not. The deployed statistic is a maximum over W correlated candidates evaluated at every sample, and the measured run length is shorter by a factor that saturates near 3.6 as W grows (measured: 3.2x at W = 10, 3.7x at W = 50, 3.6x at W = 200, at a nominal 1e-3). Use ``glr_threshold_mc`` for anything that must hold a false-alarm budget; this function exists to make the size of the trap measurable, and is named accordingly. Parameters ---------- p_sample : float Per-sample, per-candidate tail probability, in (0, 1). Returns ------- float Threshold on G. """ if not 0.0 < p_sample < 1.0: raise ValueError("p_sample must be in (0, 1)") return float(chi2.isf(p_sample, 1) / 2.0) def glr_arl0_mc(h: float, window: int, trials: int = 400, rng=None, max_samples: int = 200_000) -> tuple: """Simulated mean samples to a false alarm for the windowed GLR. There is no closed form worth trusting here, so this is the primary tool rather than a check on one. Parameters ---------- h : float Threshold on G. window : int Candidate-age window W. trials : int, default 400 Independent runs. rng : numpy.random.Generator, optional max_samples : int Safety stop per run; runs still alive are counted at this length, biasing the estimate downward. Returns ------- (float, float) Mean run length and its standard error. """ rng = np.random.default_rng() if rng is None else rng w = int(window) lengths = np.empty(trials) for t in range(trials): done = 0 tail = np.zeros(0) found = None while found is None and done < max_samples: block = 4096 z = np.concatenate((tail, rng.standard_normal(block))) g = glr_stat(z, w)[len(tail):] hit = np.nonzero(g >= h)[0] if len(hit): found = done + int(hit[0]) + 1 else: done += block tail = z[-(w - 1):] if w > 1 else np.zeros(0) lengths[t] = found if found is not None else max_samples return float(lengths.mean()), float(lengths.std(ddof=1) / np.sqrt(trials)) def glr_threshold_mc(arl0_target: float, window: int, trials: int = 200, rng=None, bracket=(1.0, 30.0), iters: int = 12) -> float: """Calibrate the GLR threshold by simulating the statistic as deployed. Bisection on ``glr_arl0_mc``. This is the whole lesson of ``glr_threshold_naive`` applied: calibrate the statistic the code computes, on the data the code sees, rather than the textbook marginal of one of its ingredients. The result carries Monte Carlo error of roughly ``1/sqrt(trials)`` in the ARL, which the log-scale of the run length converts into a much smaller error in h; a few hundred trials is usually plenty for a design number. Parameters ---------- arl0_target : float Desired mean samples between false alarms. window : int Candidate-age window W. trials : int, default 200 Runs per bisection step. rng : numpy.random.Generator, optional bracket : (float, float) Initial threshold bracket. iters : int, default 12 Bisection steps. Returns ------- float Threshold h. """ rng = np.random.default_rng() if rng is None else rng lo, hi = bracket for _ in range(iters): mid = 0.5 * (lo + hi) arl, _ = glr_arl0_mc(mid, window, trials=trials, rng=rng, max_samples=int(50 * arl0_target) + 10_000) if arl < arl0_target: lo = mid else: hi = mid return 0.5 * (lo + hi) def detection_delay_mc(detector, delta: float, trials: int = 2000, burn_in: int = 200, rng=None, max_samples: int = 4_000) -> tuple: """Mean delay from a change to the alarm, with the detector warmed up. A detection delay measured from a cold start flatters the detector: the statistic begins at zero, which is where H0 keeps it anyway. Real changes arrive at a detector that has been running for a while, so this routine feeds *burn_in* samples of in-control data first and only then applies the shift. Parameters ---------- detector : callable ``detector(z) -> index of first alarm, or None``. delta : float Shift applied after the burn-in, in sigma units. trials : int, default 2000 burn_in : int, default 200 In-control samples before the change. rng : numpy.random.Generator, optional max_samples : int, default 4000 Post-change record length. A run that has not alarmed by then is dropped, which truncates the tail of the delay distribution and biases the mean downward exactly as censoring would, so keep it far above the expected delay. The routine raises if more than 1% of runs need dropping, rather than quietly returning an optimistic number. Returns ------- (float, float) Mean delay in samples and its standard error. Runs that alarm during the burn-in are excluded; they are false alarms, not detections. """ rng = np.random.default_rng() if rng is None else rng delays = [] never = 0 for _ in range(trials): pre = rng.standard_normal(burn_in) post = rng.standard_normal(max_samples) + delta z = np.concatenate((pre, post)) idx = detector(z) if idx is None: never += 1 # ran out of record: see below continue if idx < burn_in: continue # false alarm before the change delays.append(idx - burn_in + 1) if never > 0.01 * trials: raise RuntimeError( f"{never}/{trials} runs never alarmed within max_samples=" f"{max_samples}; the reported mean would be biased low") d = np.asarray(delays, dtype=float) if len(d) < 2: raise RuntimeError("too few clean detections; raise trials or h") return float(d.mean()), float(d.std(ddof=1) / np.sqrt(len(d)))