Residuals and thresholds

"I never flag the world for being strange. I flag it for being stranger than my model expected, and only after I have measured exactly how surprised I am."

A Model-Based AI Agent

This section builds directly on the state-estimation machinery of Part III. The cleanest residual you will ever meet is the Kalman innovation from Chapter 9, and we lean on the tail probabilities and Gaussian reasoning from Chapter 4. It is the model-based cousin of the control-chart view in Section 12.2: there we charted a stationary statistic, here we chart the error of a predictor that can track a moving normal.

What a residual is, and why it is the right substrate

Let a predictor produce \(\hat{y}_t\) for the quantity you actually observe, \(y_t\). The residual is simply

The predictor can be anything: a one-step forecast from an ARIMA or exponential-smoothing model, the reconstruction of an autoencoder, a physics model of the plant, or a Bayesian filter running in parallel with the real sensor. The point of subtracting is decorrelation from context. A model that already knows the daily temperature cycle, the current motor speed, or the patient's activity level folds all of that structure into \(\hat{y}_t\), so \(r_t\) no longer carries it. What remains is the part the model could not anticipate, which is exactly the part worth alarming on. This is why detecting on residuals turns hard contextual anomalies (Section 12.1) into ordinary point anomalies on a well-behaved signal: the model has already done the contextualizing.

The gold-standard residual is the filter innovation \(\nu_t = z_t - H\hat{x}_{t|t-1}\), the measurement minus its prediction. Its beauty is that a correctly tuned Kalman filter tells you not only the innovation but its covariance \(S_t = HP_{t|t-1}H^\top + R\), so you know precisely how large an innovation should look when nothing is wrong. That paired uncertainty is what makes the whitening step below principled rather than a guess.

Whitening: turning a residual into a comparable surprise

A residual of \(+3\) is large for a temperature channel that normally errs by \(0.1\,^\circ\)C and trivial for one that swings by \(50\). To compare across channels and to attach a probability to a residual, you must normalize by its expected scale. For a scalar residual with healthy-state mean \(\mu_r\) and standard deviation \(\sigma_r\), the standardized residual (a running z-score) is

For a vector residual you cannot divide componentwise, because the components are usually correlated. The correct generalization weights by the inverse covariance, giving the squared Mahalanobis distance, which for a filter innovation is the normalized innovation squared (NIS):

The payoff is distributional. If the healthy residual is Gaussian, \(s_t\) is a standard normal and \(d_t^2\) follows a chi-square distribution with degrees of freedom equal to the measurement dimension. That is what lets you translate "how big is this residual" into "how improbable is this residual under normal," which is the currency a threshold actually spends. A crucial practical caution: estimate \(\mu_r\), \(\sigma_r\), and the covariance from a window that contains only healthy data. If a fault leaks into the statistics you normalize by, the detector calibrates itself to think the fault is normal, a leakage failure mode we return to in Chapter 5.

Setting the threshold so the false-alarm rate is a decision

A threshold \(\tau\) fires an alarm whenever \(|s_t| > \tau\) (or \(d_t^2 > \tau\)). The amateur move is to pick \(\tau = 3\) because "three sigma" sounds rigorous. The professional move is to start from the false-alarm rate you can afford and solve backwards for \(\tau\). If \(s_t\) is standard normal under normal operation, the per-sample false-alarm probability is

where \(\Phi\) is the standard-normal CDF. Invert it: \(\tau = \Phi^{-1}(1 - \alpha/2)\). Three sigma is \(\alpha \approx 0.0027\), which sounds tiny until you remember rate. At 100 Hz that is roughly one false alarm every four seconds per channel, a pager that never stops. The threshold is not really about sigmas; it is about how many false alarms per hour your operators will tolerate, a tradeoff we make explicit in Section 12.7.

Two robustness upgrades matter in the field. First, real residual tails are almost never exactly Gaussian, so a Gaussian \(\tau\) can be badly optimistic. When false alarms are expensive, set the threshold from an empirical high quantile of the healthy residuals, or model the tail directly with extreme value theory, rather than trusting the normal approximation out where it is weakest. Second, estimate the scale robustly. A single spike in the calibration window inflates \(\sigma_r\) and desensitizes the detector; using the median absolute deviation, \(\hat{\sigma} = 1.4826 \cdot \mathrm{MAD}\), keeps the scale honest in the presence of the very outliers you are hunting.

Adaptive thresholds and residual smoothing

Two refinements separate a toy detector from a deployable one. A single large residual is often just sensor noise, whereas a fault produces a small bias that persists. Thresholding a smoothed residual, an exponentially weighted moving average, trades a little detection latency for a large drop in false alarms because it integrates evidence across samples instead of gambling on one. Second, normal behavior drifts: sensors age, seasons turn, a machine is re-commissioned. A fixed threshold slowly goes stale, so the scale estimate should be refreshed on a rolling window of recent healthy data, with the important guard that you never fold flagged samples back into the calibration statistics. The code below assembles these pieces: a one-step predictor, a robust standardized residual, an EWMA smoother, and a false-alarm-driven threshold.

import numpy as np

def robust_scale(x):                 # MAD-based sigma, spike-resistant
    med = np.median(x)
    return 1.4826 * np.median(np.abs(x - med)) + 1e-9

def residual_detector(y, warmup=200, alpha=1e-3, lam=0.2):
    yhat = np.empty_like(y)          # one-step "predict last value" model
    yhat[0] = y[0]
    yhat[1:] = y[:-1]                # replace with any real forecaster
    r = y - yhat                     # the residual signal

    mu, sigma = np.median(r[:warmup]), robust_scale(r[:warmup])
    tau = np.sqrt(2) * erfinv_safe(1 - alpha)   # Gaussian false-alarm level

    s_ewma, alarms = 0.0, []
    for t in range(len(y)):
        s = (r[t] - mu) / sigma      # standardized residual
        s_ewma = lam * s + (1 - lam) * s_ewma    # smooth the evidence
        alarms.append(abs(s_ewma) > tau)
    return np.array(alarms), tau

def erfinv_safe(p):                  # tau from a target false-alarm rate
    from scipy.special import erfinv
    return erfinv(p)

rng = np.random.default_rng(0)
y = np.r_[rng.normal(0, 1, 600), rng.normal(0, 1, 400) + 2.5]  # step fault
alarms, tau = residual_detector(y)
print(f"threshold={tau:.2f}  first alarm at t={np.argmax(alarms)}")

The detector above finds the step change because the EWMA lets a small persistent bias accumulate past a threshold that a single sample would never cross. That accumulation idea, turned into an optimal sequential test, is exactly the CUSUM procedure of the next section.

What's Next

In Section 12.4, we stop treating each smoothed residual as an independent test and ask a sharper question: given a running stream of residuals, what is the statistically optimal way to declare that their distribution has shifted? That is the change-point problem, and its two classic answers, the cumulative-sum (CUSUM) chart and Bayesian online change-point detection, turn the ad hoc EWMA of this section into a principled sequential decision.