Detecting Outliers with Normalized Residuals

Once a least-squares control-network adjustment has converged, the next obligation is to prove that no single observation is a blunder corrupting the whole solution — the disciplined way to do that is Baarda data snooping on normalized residuals, a technique that sits directly under least-squares adjustment for control networks within the broader Algorithmic Math & Geodetic Workflows reference. A raw residual of five millimetres means nothing on its own; it only becomes actionable once it is scaled by how well the network geometry could actually detect an error at that observation.

Figure — data snooping rejects at most one observation per pass, then re-adjusts.

What Normalizing the Residual Actually Does #

Least-squares distributes an error across many observations, so a gross blunder in one measurement leaks into the residuals of its well-connected neighbours. The residual vector $\mathbf{v}$ therefore does not point cleanly at the culprit, and its entries do not share a common scale — a residual at a high-redundancy observation is far more diagnostic than one at a barely-checked observation. Baarda’s insight was to divide each residual by its own standard deviation so that, under the null hypothesis of no blunder, every scaled residual follows the same standard normal distribution. The normalized (studentized) residual is

Figure — remove one outlier per pass; removing several at once is how good observations get deleted.

where $v_i$ is the $i$-th residual, $\sigma_0$ is the a-priori reference standard deviation (unit weight), and $q_{v_i v_i}$ is the $i$-th diagonal entry of the cofactor matrix of the residuals, $\mathbf{Q}_{vv} = \mathbf{Q}_{\ell\ell} - \mathbf{A}\,\mathbf{N}^{-1}\mathbf{A}^{\mathsf T}$. That diagonal is the local redundancy number carried in cofactor form: it shrinks toward zero where the network cannot check an observation, which is exactly why raw residuals mislead. Each $w_i$ is then compared against a critical value $c$ — commonly $3.29$ for a per-observation significance of $\alpha = 0.001$ — and the observation with the largest $|w_i|$ above $c$ is flagged as the most likely single blunder. This is the same variance-scaling discipline that underpins error distribution modeling in Python, where residuals are propagated through the full covariance rather than judged in raw metres.

Complete Runnable Implementation #

The routine below takes the residual vector, the residual cofactor matrix $\mathbf{Q}_{vv}$, and the reference standard deviation, then returns each normalized residual and a single-blunder flag following the one-at-a-time data-snooping rule. It never mutates its inputs and clamps near-zero cofactor diagonals so an unchecked observation cannot masquerade as a huge normalized residual.

from __future__ import annotations

import logging
from dataclasses import dataclass

import numpy as np

logger = logging.getLogger("data_snooping")

# Baarda critical values for common per-observation significance levels.
BAARDA_CRITICAL: dict[float, float] = {0.05: 1.96, 0.01: 2.58, 0.001: 3.29}


@dataclass(frozen=True)
class SnoopingResult:
    normalized_residuals: np.ndarray  # w_i, one per observation
    flagged_index: int                # index of the worst blunder, -1 if none
    flagged_w: float                  # |w| at that index (0.0 if none)
    critical_value: float             # c used for the test


def detect_outlier(
    residuals: np.ndarray,
    q_vv: np.ndarray,
    sigma0: float = 1.0,
    alpha: float = 0.001,
    zero_floor: float = 1e-12,
) -> SnoopingResult:
    """Flag the single most likely blunder by Baarda data snooping.

    Computes the normalized residual w_i = v_i / (sigma0 * sqrt(q_vv,i))
    for every observation and returns the index of the largest |w_i| that
    exceeds the critical value, or -1 when the network is clean.

    Args:
        residuals: length-n residual vector v from the adjustment (metres).
        q_vv: n-by-n residual cofactor matrix Q_vv, or its length-n diagonal.
        sigma0: a-priori reference standard deviation of unit weight.
        alpha: per-observation significance level (0.05, 0.01, or 0.001).
        zero_floor: cofactor-diagonal values at or below this are treated as
            uncheckable (redundancy ~ 0) and excluded from the test.

    Returns:
        A SnoopingResult with the full w vector and the worst-blunder flag.
    """
    v = np.asarray(residuals, dtype=np.float64).ravel()
    q = np.asarray(q_vv, dtype=np.float64)
    q_diag = np.diag(q) if q.ndim == 2 else q.ravel()
    if q_diag.shape != v.shape:
        raise ValueError(
            f"residuals ({v.shape}) and q_vv diagonal ({q_diag.shape}) disagree"
        )
    if sigma0 <= 0.0:
        raise ValueError("sigma0 must be positive")

    checkable = q_diag > zero_floor
    w = np.zeros_like(v)
    denom = sigma0 * np.sqrt(np.where(checkable, q_diag, 1.0))
    w[checkable] = v[checkable] / denom[checkable]

    critical = BAARDA_CRITICAL.get(alpha)
    if critical is None:
        raise ValueError(f"no tabulated critical value for alpha={alpha}")

    abs_w = np.abs(w)
    worst = int(np.argmax(abs_w)) if abs_w.size else -1
    if worst < 0 or abs_w[worst] <= critical:
        return SnoopingResult(w, -1, 0.0, critical)
    logger.info("blunder at obs %d: |w|=%.3f > c=%.2f", worst, abs_w[worst], critical)
    return SnoopingResult(w, worst, float(abs_w[worst]), critical)

Parameter Reference #

Worked Example #

A five-observation levelling loop with one deliberately corrupted measurement makes the scaling visible: the blunder does not have the largest raw residual, yet it has the largest normalized residual.

import numpy as np

v = np.array([0.004, -0.003, 0.021, 0.005, -0.002])          # metres
q_vv = np.array([0.62, 0.55, 0.71, 0.10, 0.48])              # cofactor diagonal
res = detect_outlier(v, q_vv, sigma0=0.005, alpha=0.001)
print(np.round(res.normalized_residuals, 2))
print(res.flagged_index, round(res.flagged_w, 2))
# [ 1.02 -0.81  4.98  3.16 -0.58]
# 2 4.98

Validation Check #

Gate the flag before it is used to reject an observation, and re-adjust after every single removal rather than dropping several at once.

Figure — the critical value depends on how many observations you are screening.

assert res.flagged_index == 2, "expected the corrupted observation to be flagged"
assert res.flagged_w > res.critical_value, "flagged w must exceed the critical value"
assert np.isfinite(res.normalized_residuals).all(), "non-finite normalized residual"

Common Mistakes #

Frequently Asked Questions #

Why normalise by the residual’s standard deviation rather than the observation’s? #

Because the adjustment redistributes error: a well-checked observation shows most of its error in its residual, while an observation with little redundancy shows almost none, however wrong it is. Dividing by the observation’s own standard deviation ignores that and systematically under-flags exactly the observations least protected by geometry. The residual covariance matrix carries the right denominator.

Can I remove all the observations that exceed the critical value at once? #

No, and doing so is how good observations get deleted. The presence of one large blunder distorts every other residual, so several observations can exceed the threshold purely because of one bad one. Remove the single largest, re-adjust, and test again. The sequence usually collapses quickly, and the ones that were only guilty by association come back inside the threshold on the second pass.

What significance level should the test use? #

Something small, and corrected for the number of observations being tested. Screening a thousand observations at a nominal 0.1 per cent will flag roughly one good observation per run by chance alone, which trains everyone to ignore flags. Dividing the level by the number of tests — a Bonferroni correction — raises the critical value to about 4.4 for a thousand observations and keeps the flags meaningful.

Related #

• Least-squares adjustment for control networks — the parent adjustment that produces the residuals and cofactor tested here.
• Error distribution modeling in Python — how the surviving residuals are propagated through the full covariance.
• Python least-squares adjustment for control points — the affine fit whose residuals are screened the same way.
• Algorithmic Math & Geodetic Workflows — the pipeline that adjustment, snooping, and covariance propagation feed.