RouteMatching

Hampel Filter vs Z-Score for GPS Spike Rejection

The consensus scoring approach in outlier removal for raw telematics streams uses IQR bounds computed once over an entire trip. A rolling detector — recomputing a local baseline as it slides across the series — catches a narrower class of problem: a positional spike that is only anomalous relative to its immediate neighbours, not the whole trip’s distribution. The obvious first implementation reaches for .rolling().mean() and .rolling().std(), because pandas ships both for free. That rolling z-score has a specific, well-documented failure mode against the exact spikes it is meant to catch: the spike inflates the very statistics used to detect it. This page compares that rolling z-score against a rolling Hampel identifier (rolling median + median absolute deviation), which is built from the ground up to resist this failure, and ships a comparison harness so you can validate the difference on your own fleet data before switching detectors in production.

Rolling z-score masking vs rolling Hampel robustness The same raw GPS series is processed by two detectors in parallel. The rolling z-score branch shows the spike inflating mean and standard deviation, which masks the spike. The rolling Hampel branch shows the median and MAD staying stable, so the spike is correctly flagged. Raw GPS series window w, spike at t Rolling Z-Score mean ± n_sigmas·std spike moves both Rolling Hampel median ± t0·MAD spike moves neither Masking inflated std lowers z-score spike may go undetected Robust detection MAD unmoved by spike spike reliably flagged

Compatibility and Configuration Requirements

Requirement Minimum version / value Notes
Python 3.10 Type hints used throughout the module below
pandas 2.0 .rolling(center=True) and .to_numpy() behaviour assumed
numpy 1.24 np.median, vectorised absolute-deviation computation
Input series 1-D numeric, time-ordered Derive from implied speed, per-axis coordinate delta, or heading rate — not raw lat/lon pairs directly
window Odd full-window size recommended Passed here as a half-window w; full window is 2w + 1 samples
n_sigmas / Hampel t0 2.5-4.0 Same parameter name used for both detectors in this module so their thresholds are directly comparable
MAD scale constant 1.4826 Fixed; rescales MAD to be a consistent estimator of standard deviation under normality

Python Module: Rolling Hampel Filter, Rolling Z-Score, and a Comparison Harness

The module below implements both detectors against the same interface and window semantics, plus a compare_detectors function that highlights exactly where the two disagree — the rows worth inspecting manually before trusting either one in production.

import numpy as np
import pandas as pd


MAD_SCALE = 1.4826  # rescales MAD to be a consistent estimator of std under normality


def rolling_hampel_filter(
    series: pd.Series,
    window: int = 7,
    n_sigmas: float = 3.0,
) -> pd.DataFrame:
    """
    Rolling Hampel identifier for a 1-D GPS-derived series.

    Parameters
    ----------
    series : pd.Series
        Time-ordered numeric series, e.g. implied speed (km/h) or a
        per-axis coordinate delta (metres).
    window : int
        Half-window size in samples. The full window spans
        2 * window + 1 points centred on each sample. 5-9 is typical
        for 1 Hz fleet telematics; widen for sparser loggers.
    n_sigmas : float
        Threshold multiplier (commonly written t0) applied to the
        scaled MAD. 3.0 is a standard default; lower toward 2.0 for
        more aggressive rejection, raise toward 4.0-5.0 to only catch
        extreme spikes.

    Returns
    -------
    pd.DataFrame with columns: value, rolling_median, mad, is_outlier
    """
    values = series.to_numpy(dtype=float)
    n = len(values)
    rolling_median = np.full(n, np.nan)
    mad = np.full(n, np.nan)
    is_outlier = np.zeros(n, dtype=bool)

    for i in range(n):
        lo = max(0, i - window)
        hi = min(n, i + window + 1)
        local = values[lo:hi]
        med = np.median(local)
        scaled_mad = MAD_SCALE * np.median(np.abs(local - med))
        rolling_median[i] = med
        mad[i] = scaled_mad
        # Guard against a zero MAD in a perfectly flat local window,
        # which would otherwise flag any tiny deviation as an outlier.
        threshold = n_sigmas * scaled_mad
        is_outlier[i] = scaled_mad > 0 and abs(values[i] - med) > threshold

    return pd.DataFrame({
        "value": values,
        "rolling_median": rolling_median,
        "mad": mad,
        "is_outlier": is_outlier,
    }, index=series.index)


def rolling_modified_zscore(
    series: pd.Series,
    window: int = 7,
    n_sigmas: float = 3.0,
) -> pd.DataFrame:
    """
    Rolling z-score outlier detector using a centred rolling mean and
    standard deviation, for direct comparison against the Hampel
    filter above at the same window size and threshold multiplier.

    Parameters
    ----------
    series : pd.Series
        Same input contract as rolling_hampel_filter.
    window : int
        Half-window size; converted internally to the same
        2 * window + 1 full window used by the Hampel filter.
    n_sigmas : float
        Number of standard deviations from the rolling mean beyond
        which a point is flagged. Same parameter name as the Hampel
        filter's n_sigmas so both detectors can be run with identical
        strictness for comparison.

    Returns
    -------
    pd.DataFrame with columns: value, rolling_mean, rolling_std, is_outlier
    """
    full_window = 2 * window + 1
    rolling_mean = series.rolling(full_window, center=True, min_periods=1).mean()
    rolling_std = series.rolling(full_window, center=True, min_periods=1).std(ddof=0)

    deviation = (series - rolling_mean).abs()
    threshold = n_sigmas * rolling_std
    is_outlier = (rolling_std > 0) & (deviation > threshold)

    return pd.DataFrame({
        "value": series.to_numpy(dtype=float),
        "rolling_mean": rolling_mean.to_numpy(),
        "rolling_std": rolling_std.to_numpy(),
        "is_outlier": is_outlier.to_numpy(),
    }, index=series.index)


def compare_detectors(
    series: pd.Series,
    window: int = 7,
    n_sigmas: float = 3.0,
) -> pd.DataFrame:
    """
    Run both detectors over the same series, window, and threshold and
    return a side-by-side comparison, including the disagreement rows
    that matter most when validating a migration from z-score to
    Hampel on real fleet data.

    Returns
    -------
    pd.DataFrame with columns: value, hampel_flag, zscore_flag, agree,
    hampel_only, zscore_only
    """
    hampel = rolling_hampel_filter(series, window=window, n_sigmas=n_sigmas)
    zscore = rolling_modified_zscore(series, window=window, n_sigmas=n_sigmas)

    comparison = pd.DataFrame({
        "value": series.to_numpy(dtype=float),
        "hampel_flag": hampel["is_outlier"].to_numpy(),
        "zscore_flag": zscore["is_outlier"].to_numpy(),
    }, index=series.index)
    comparison["agree"] = comparison["hampel_flag"] == comparison["zscore_flag"]
    comparison["hampel_only"] = comparison["hampel_flag"] & ~comparison["zscore_flag"]
    comparison["zscore_only"] = comparison["zscore_flag"] & ~comparison["hampel_flag"]

    return comparison


# ---------------------------------------------------------------------------
# Usage example
# ---------------------------------------------------------------------------
# implied_speed = df["implied_speed_kmh"]  # from the Haversine spatial check
#
# result = compare_detectors(implied_speed, window=7, n_sigmas=3.0)
# print(result["hampel_only"].sum(), "spikes only the Hampel filter caught")
# print(result["zscore_only"].sum(), "points only the z-score flagged")
#
# clean = df.loc[~result["hampel_flag"]].reset_index(drop=True)

Key parameter notes:

  • Both detectors share the window and n_sigmas argument names deliberately, so a fleet migrating from the z-score to the Hampel filter can run compare_detectors at their existing production threshold before changing anything downstream.
  • The Hampel filter’s inner loop is a plain Python for loop over np.median calls, which is clear to audit but not the fastest option at fleet scale — see the computational cost pitfall below.
  • rolling_modified_zscore uses ddof=0 (population standard deviation) rather than the pandas default ddof=1, to stay consistent with the Hampel filter’s population-style MAD computation at the same window size.

Comparison: Breakdown Point, Masking, and Computational Cost

Property Rolling Hampel (median + MAD) Rolling Z-Score (mean + std)
Breakdown point ~50% — the median tolerates up to half the window being corrupted before it moves 0% — a single extreme value can shift the mean and inflate the standard deviation arbitrarily
Masking effect Resistant — the spike is measured against a baseline it cannot itself move Present — the spike inflates its own std, which can push its z-score back under the threshold
Swamping effect (false positives) Low — normal points near a spike are still compared against a stable median Moderate — an inflated std widens the accepted range for every point in the window, but can also flag borderline-normal points when a spike pulls the mean sideways
Computational cost per point O(w log w) for the sort inside each local np.median call O(1) amortised with a running-sum formulation, or O(w) with the naive .rolling() implementation above
pandas built-in support None — requires the manual loop shown above, or a vectorised/Numba rewrite for scale Native via .rolling().mean() and .rolling().std()
Best suited to Streams where spikes can cluster (urban canyon multipath lasting several consecutive fixes) Pre-cleaned streams where remaining outliers are already known to be isolated single points

Execution and Tuning Guidelines

  • window — 5-9 (half-window) is typical for 1 Hz commercial telematics. Widen it for sparse loggers (1 ping per 30-60 s) so the local baseline still has enough points to compute a meaningful median and MAD; a half-window of 2-3 on sparse data makes both detectors unstable regardless of which one you choose.
  • n_sigmas / Hampel t0 — 3.0 is a reasonable default for both detectors when run through compare_detectors at the same value. Lower toward 2.0 if you need aggressive rejection ahead of a sensitive downstream stage such as Kalman filtering; raise toward 4.0-5.0 if legitimate hard-braking or sharp-turn events are being flagged as spikes.
  • MAD scale (1.4826) — leave this fixed unless you have a specific statistical reason to change it. It is a normality-derived constant, not a tuning knob; changing it decouples the Hampel n_sigmas from the same physical meaning as the z-score’s n_sigmas, which breaks the like-for-like comparison this harness is built for.
  • Choosing between detectors — run compare_detectors on a representative sample of trips first. If hampel_only counts dominate and correspond visually to clustered spikes on a map tile, migrate the production pipeline to the Hampel filter. If the two detectors agree on almost every row, the cheaper rolling z-score is defensible, particularly at very high data volumes where the O(w log w) Hampel cost matters.
  • Feeding the result forward — treat is_outlier exactly like the flags in the base outlier removal pipeline: either drop flagged rows outright or fold the flag into that pipeline’s consensus score alongside the acceleration, heading-rate, and spatial-jump checks rather than filtering on this signal in isolation.

Common Pitfalls

Trusting the z-score on a stream with clustered spikes

Cause: Urban canyon multipath reflections often persist for two or three consecutive GPS fixes, not just one. When two adjacent points in the rolling window are both spikes, they inflate the mean and standard deviation together, and can mask each other even more severely than a single isolated spike would.

Symptom: zscore_only and shared-agreement counts look reasonable on a spot-check of isolated spikes, but a full-trip comparison against rolling_hampel_filter shows a run of consecutive points the z-score missed entirely inside a known urban-canyon segment.

Fix: Run compare_detectors specifically on trip segments known to pass through dense urban corridors before deciding a rolling z-score is sufficient. If clustering is common in your fleet’s operating environment, default to the Hampel filter rather than the z-score.

Zero MAD in a flat local window over-triggers or silently passes everything

Cause: A vehicle stationary for longer than the window duration produces a local window where every value is identical (or nearly so after rounding), making the local MAD exactly zero. Without a guard, dividing by a zero MAD either raises a runtime warning or, if implemented carelessly as a ratio rather than a threshold comparison, flags every non-identical value — including harmless floating-point noise — as an outlier.

Symptom: A burst of is_outlier = True flags appears the moment a vehicle starts moving again after being parked, even though the actual movement is legitimate.

Fix: The scaled_mad > 0 guard in rolling_hampel_filter above prevents the false trigger by treating a zero-MAD window as “insufficient variation to detect an outlier” rather than as an implicit zero threshold. Keep this guard if you rewrite the loop for performance.

Computational cost blows up on fleet-scale batches

Cause: The reference rolling_hampel_filter implementation above runs a plain Python loop with an np.median call per row. At a few thousand rows per trip this is fine; across tens of millions of rows for a full fleet-day batch, the per-row Python overhead dominates and the job can take hours where the equivalent .rolling().std() z-score finishes in seconds.

Symptom: The outlier-detection stage becomes the slowest step in the batch pipeline, disproportionate to its conceptual simplicity, once the fleet scales past a few hundred vehicles.

Fix: Replace the Python loop with pandas’ .rolling().median() for the median term (which is implemented in compiled code) combined with a vectorised MAD computed via .rolling().apply() using a Numba-accelerated function, or move the computation to a compiled bottleneck library. Validate that the optimised version produces is_outlier flags identical to the reference loop on a held-out sample before replacing it in production, since a subtly different window-edge convention can silently change results near the start and end of each trip.

Up: Outlier Removal in Raw Telematics Streams | GPS Data Preprocessing & Cleaning Fundamentals

Related

Up to section Outlier Removal in Raw Telematics Streams Related Automating Outlier Detection in High-Frequency Telematics Data