Robust filters (standalone detectors)

You already have a recorded series: a joint velocity at 100 Hz, a temperature log, a sensor reading. A few samples in it are wrong, and you want to know which ones and to see the series with them replaced. That is what this page does.

The library has four classic batch filters for this job: Hampel, Tukey fences, Chauvenet and sigma clipping. Each one takes the whole series, compares every sample with its neighbours, and returns two things: the series with the flagged samples replaced by a local estimate, and the list of positions it flagged. The four differ in how they decide what "strange" means, so they often disagree on borderline samples.

pressure_valve settles the disagreement: it runs all four, weighs them by how confident each one is, and returns one cleaned series, one list of flagged positions, and one pressure score per sample. The sections below show each method in turn, then the combination.

These filters look at samples after the current one (the window is centred on it). For a series that is arriving now, where looking ahead is not allowed, use the causal versions on Streaming.

1. Hampel filter

The batch Hampel filter uses a window centred on each point (it sees the future) and replaces outliers with the window median.

import numpy as np
from dense_armor.utility.anomaly.robust_filters import hampel_filter

rng = np.random.default_rng(42)
v = 0.5 + 0.1 * rng.standard_normal(3000)
v[1500] = 2.5
clean, idx = hampel_filter(v, radius=15, n_sigmas=3.0)
print(1500 in idx, round(float(clean[1500]), 3), len(idx))

True 0.498 35

The stream is a joint velocity at 100 Hz: a level of 0.5, Gaussian noise with standard deviation 0.1, and one collision at sample 1500 that drives the reading to 2.5.

The Hampel filter flags the collision and replaces it with the local median, 0.498. It also flags 34 ordinary samples: with a 31-point window and a threshold of 3 in the Hampel rule, some noise values fall beyond it. The total is 35 flagged positions.

Every function on this page returns two values, clean and idx. clean is the input series where every unflagged sample keeps its value and every flagged sample is replaced by the local median. idx is the list of flagged positions.

The rule is |x_k - m_k| > t * S_k, where m_k is the median of the window centred on sample k, S_k = 1.4826 * MAD is the robust estimate of the standard deviation of that window, and t is the threshold (Pearson et al., 2016, eqs. 3 and 4; step by step on Streaming).

2. The other three

tukey_fences, chauvenet_criterion and sigma_clip implement the other three classic rules, with the same (clean, idx) return as hampel_filter.

import numpy as np
from dense_armor.utility.anomaly.robust_filters import (
    hampel_filter, tukey_fences, chauvenet_criterion, sigma_clip)

rng = np.random.default_rng(42)
v = 0.5 + 0.1 * rng.standard_normal(3000)
v[1500] = 2.5
for f in (hampel_filter, tukey_fences, chauvenet_criterion, sigma_clip):
    clean, idx = f(v, radius=15)
    print(f.__name__, 1500 in idx, len(idx))

hampel_filter True 35
tukey_fences True 65
chauvenet_criterion True 48
sigma_clip True 24

All four catch the collision; they differ in how many ordinary samples they also flag. Tukey uses the quartiles (outside $[Q_1 - 1.5\,\mathrm{IQR},\, Q_3 + 1.5\,\mathrm{IQR}]$), Chauvenet rejects when the expected number of samples that far out, $N \cdot P$, is below 0.5, sigma clipping removes values beyond n_sigmas standard deviations and repeats on what remains.

3. pressure_valve: four estimates combined

Each of the four methods (Hampel, Tukey fences, Chauvenet, sigma clipping) returns two numbers for every sample. The first is a centre c_k: the method's estimate of the true value at that sample. The second is a scale s_k: the method's estimate of its own uncertainty. Both come from the window around the sample, and the pair differs per method:

method	centre c_k	scale s_k
Hampel	median	1.4826 * MAD
Tukey fences	median	IQR / 1.349
Chauvenet	mean	standard deviation
sigma clipping	clipped mean	standard deviation

The four scales can differ a lot: a method whose window already contains an outlier reports a large scale, a method whose window is clean reports a small one.

pressure_valve weights each centre by the inverse of its variance, so a method that is sure of itself counts more than one that is not. The weight of method k is

$$w_k = \frac{1/s_k^2}{\sum_j 1/s_j^2},$$

the minimum-variance (inverse-variance) combination of independent estimates. No weight is chosen by hand.

import numpy as np

s = np.array([0.10, 0.12, 0.10, 0.30])
w = (1 / s**2) / np.sum(1 / s**2)
print(np.round(w, 3), round(float(1 / np.sqrt(np.sum(1 / s**2))), 4))

[0.356 0.248 0.356 0.04 ] 0.0597

A method with a small scale gets a large weight; a method with a large scale gets a small one. The combined estimate is more precise than any single method, and its scale is

$$s = \frac{1}{\sqrt{\sum_j 1/s_j^2}}.$$

The pressure of a new sample x is how far it lies from the combined centre, measured in units of the combined scale:

$$p = \frac{\left|x - \sum_k w_k c_k\right|}{s}.$$

In the code above the four scales are 0.10, 0.12, 0.10 and 0.30. The method with scale 0.30 gets weight 0.04; the two methods with scale 0.10 get weight 0.36 each. The combined scale is 0.0597, smaller than any of the four. The four weights are the solution of the calculation above, not numbers chosen by hand. This is the minimum-variance combination of independent estimates.

4. pressure_valve on the joint signal

import numpy as np
from dense_armor.utility.anomaly.robust_filters import pressure_valve

rng = np.random.default_rng(42)
v = 0.5 + 0.1 * rng.standard_normal(3000)
v[1500] = 2.5
clean, idx, p, thr = pressure_valve(v, radius=15)
print(idx, round(float(p[1500]), 1), round(float(np.median(thr)), 2))

[949, 1500, 2707] 42.7 8.46

Of the 3,000 samples, the filter flags three: the collision at 1500 (pressure 42.7) and two noise values at 949 and 2707.

The threshold is not the same at every sample. At each sample it depends on how different the local window looks from a wider reference window of size radius * ref_mult. The difference is measured by the Jensen–Shannon divergence, a number in [0, 1]: near 0 when the two windows have the same shape (steady noise), near 1 when they have different shapes (a real change of level is under way). The threshold is

$$\text{threshold}i = \text{soglia_pressione} \cdot (1 + k{\text{molla}} \cdot \mathrm{JSD}_i),$$

where soglia_pressione (default 8.0) and k_molla are the two parameters of the code.

When the two windows look the same, JSD is near 0 and the threshold stays at its base value, 8.0. When they differ, the threshold rises and the filter tolerates a larger deviation, so a genuine change of level is not flagged as an anomaly. On the 3,000 samples above the median threshold is 8.46.

API reference

robust_filters

robust_filters.py -- quattro rilevatori di anomalie "a basso costo", nessuna meccanica pesante (nessun modello dinamico, nessuno stato, nessuna differenziazione automatica): statistiche classiche, note e usate da decenni in telemetria/astronomia, ognuna calcolata su una finestra locale centrata (come healing_filter.py in questo stesso pacchetto -- non causale, pensate per pulizia offline/batch di una serie gia' registrata, non per il ciclo real-time di core/hybrid_engine.py).

Ogni funzione ha la stessa firma (x, radius=..., **soglia) -> (pulito, anomalie): pulito -- array 1D, stessa shape di x: valore grezzo dove non anomalo, mediana locale dove anomalo. anomalie -- lista di indici (int) segnalati come anomali.

I quattro metodi, in ordine di introduzione storica: - Chauvenet (1863): soglia di rigetto basata sulla probabilita' attesa che un campione della finestra locale (size N) si discosti tanto per puro caso -- se il numero atteso di tali eventi (N * P) e' < 0.5, il punto e' un outlier statisticamente "troppo improbabile per il campione a disposizione". - Tukey's fences / IQR (anni '70): nessuna assunzione sulla distribuzione, solo quartili -- oltre Q1-1.5IQR o Q3+1.5IQR e' outlier. - Hampel filter: mediana + MAD scalato (1.4826*MAD, costante di consistenza per una gaussiana) come soglia -- la generalizzazione robusta di "3 sigma" che non usa mai una media/std sensibile agli outlier che sta cercando. - Sigma-clipping iterativo: standard in astronomia (stacking di immagini, fotometria) -- ricalcola media/std della finestra escludendo i punti oltre N sigma, ripete finche' stabile o max_iters, poi giudica il punto centrale contro le statistiche "ripulite".

chauvenet_criterion

chauvenet_criterion(
    x: ndarray, radius: int = 10
) -> Tuple[np.ndarray, List[int]]

Criterio di Chauvenet (1863) su finestra locale di raggio radius.

Per il punto i, con la finestra come campione (media mu, std sigma, size N): d = |x[i] - mu| / sigma; P = erfc(d / sqrt(2)) (probabilita' a due code di un tale scostamento sotto ipotesi gaussiana). Se N*P < 0.5 -- il numero atteso di scostamenti pari a questo, ripetendo l'esperimento, e' inferiore a mezzo campione -- il punto e' rigettato come outlier.

Nota: usa mean/std (non robuste), come nella formulazione originale -- su una finestra con PIU' di un outlier, mean/std sono gia' distorte dagli altri outlier prima ancora di giudicare il punto i (limite noto e documentato del criterio originale, non specifico di questa implementazione; per finestre con outlier multipli preferire hampel_filter o tukey_fences, che usano statistiche robuste).

Source code in dense_armor/utility/anomaly/robust_filters.py

def chauvenet_criterion(x: np.ndarray, radius: int = 10) -> Tuple[np.ndarray, List[int]]:
    """Criterio di Chauvenet (1863) su finestra locale di raggio `radius`.

    Per il punto i, con la finestra come campione (media mu, std sigma, size
    N): d = |x[i] - mu| / sigma; P = erfc(d / sqrt(2)) (probabilita' a due code
    di un tale scostamento sotto ipotesi gaussiana). Se N*P < 0.5 -- il numero
    atteso di scostamenti pari a questo, ripetendo l'esperimento, e' inferiore
    a mezzo campione -- il punto e' rigettato come outlier.

    Nota: usa mean/std (non robuste), come nella formulazione originale --
    su una finestra con PIU' di un outlier, mean/std sono gia' distorte dagli
    altri outlier prima ancora di giudicare il punto i (limite noto e
    documentato del criterio originale, non specifico di questa implementazione;
    per finestre con outlier multipli preferire hampel_filter o tukey_fences,
    che usano statistiche robuste)."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    out = np.copy(x)
    anomalie = []
    for i in range(n):
        w = _window(x, i, radius)
        if w.size < 3:
            continue
        mu, sigma = float(np.mean(w)), float(np.std(w))
        if sigma < 1e-12:
            continue
        d = abs(x[i] - mu) / sigma
        p = math.erfc(d / math.sqrt(2))
        if w.size * p < 0.5:
            anomalie.append(i)
            out[i] = float(np.median(w))
    return out, anomalie

tukey_fences

tukey_fences(
    x: ndarray, radius: int = 10, k: float = 1.5
) -> Tuple[np.ndarray, List[int]]

Tukey's fences / IQR su finestra locale di raggio radius.

Q1/Q3 = 25/75-esimo percentile della finestra, IQR = Q3-Q1. Il punto i e' outlier se fuori da [Q1 - kIQR, Q3 + kIQR] (k=1.5 e' lo standard, k=3.0 la variante "extreme outliers"). Nessuna assunzione sulla distribuzione -- l'unico dei quattro a non presupporre normalita'.

Source code in dense_armor/utility/anomaly/robust_filters.py

def tukey_fences(x: np.ndarray, radius: int = 10, k: float = 1.5) -> Tuple[np.ndarray, List[int]]:
    """Tukey's fences / IQR su finestra locale di raggio `radius`.

    Q1/Q3 = 25/75-esimo percentile della finestra, IQR = Q3-Q1. Il punto i e'
    outlier se fuori da [Q1 - k*IQR, Q3 + k*IQR] (k=1.5 e' lo standard, k=3.0
    la variante "extreme outliers"). Nessuna assunzione sulla distribuzione --
    l'unico dei quattro a non presupporre normalita'."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    out = np.copy(x)
    anomalie = []
    for i in range(n):
        w = _window(x, i, radius)
        if w.size < 4:
            continue
        q1, q3 = np.percentile(w, [25, 75])
        iqr = q3 - q1
        lo_fence, hi_fence = q1 - k * iqr, q3 + k * iqr
        if x[i] < lo_fence or x[i] > hi_fence:
            anomalie.append(i)
            out[i] = float(np.median(w))
    return out, anomalie

hampel_filter

hampel_filter(
    x: ndarray, radius: int = 10, n_sigmas: float = 3.0
) -> Tuple[np.ndarray, List[int]]

Hampel filter su finestra locale di raggio radius.

mediana e MAD (Median Absolute Deviation) della finestra, scalato per 1.4826 (costante di consistenza: rende scaled_MAD uno stimatore non distorto della deviazione standard SE il dato fosse gaussiano, ma senza mai calcolare una media/std vera -- ogni singola statistica usata qui e' gia' robusta, a differenza di Chauvenet). Outlier se |x[i] - mediana| > n_sigmas * scaled_MAD.

Source code in dense_armor/utility/anomaly/robust_filters.py

def hampel_filter(x: np.ndarray, radius: int = 10, n_sigmas: float = 3.0) -> Tuple[np.ndarray, List[int]]:
    """Hampel filter su finestra locale di raggio `radius`.

    mediana e MAD (Median Absolute Deviation) della finestra, scalato per
    1.4826 (costante di consistenza: rende scaled_MAD uno stimatore
    non distorto della deviazione standard SE il dato fosse gaussiano,
    ma senza mai calcolare una media/std vera -- ogni singola statistica
    usata qui e' gia' robusta, a differenza di Chauvenet). Outlier se
    |x[i] - mediana| > n_sigmas * scaled_MAD."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    out = np.copy(x)
    anomalie = []
    for i in range(n):
        w = _window(x, i, radius)
        if w.size < 3:
            continue
        med = float(np.median(w))
        mad = float(np.median(np.abs(w - med)))
        scaled_mad = 1.4826 * mad
        if scaled_mad < 1e-12:
            if abs(x[i] - med) > 1e-9:
                anomalie.append(i)
                out[i] = med
            continue
        if abs(x[i] - med) > n_sigmas * scaled_mad:
            anomalie.append(i)
            out[i] = med
    return out, anomalie

sigma_clip

sigma_clip(
    x: ndarray,
    radius: int = 10,
    n_sigmas: float = 3.0,
    max_iters: int = 5,
) -> Tuple[np.ndarray, List[int]]

Sigma-clipping iterativo (standard in astronomia: stacking, fotometria) su finestra locale di raggio radius.

Sulla finestra (esclude il punto i stesso, cosi' un vero outlier al centro non puo' auto-proteggersi gonfiando la propria std): calcola media/std, scarta i punti oltre n_sigmas, ricalcola su quelli rimasti, ripete fino a stabilita' o max_iters. Il punto i e' outlier se fuori da n_sigmas dalla media/std FINALI (ripulite dagli altri outlier della finestra, a differenza di chauvenet_criterion che usa mean/std grezze).

Source code in dense_armor/utility/anomaly/robust_filters.py

def sigma_clip(x: np.ndarray, radius: int = 10, n_sigmas: float = 3.0, max_iters: int = 5) -> Tuple[np.ndarray, List[int]]:
    """Sigma-clipping iterativo (standard in astronomia: stacking, fotometria)
    su finestra locale di raggio `radius`.

    Sulla finestra (esclude il punto i stesso, cosi' un vero outlier al
    centro non puo' auto-proteggersi gonfiando la propria std): calcola
    media/std, scarta i punti oltre n_sigmas, ricalcola su quelli rimasti,
    ripete fino a stabilita' o max_iters. Il punto i e' outlier se fuori da
    n_sigmas dalla media/std FINALI (ripulite dagli altri outlier della
    finestra, a differenza di chauvenet_criterion che usa mean/std grezze)."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    out = np.copy(x)
    anomalie = []
    for i in range(n):
        lo, hi = max(0, i - radius), min(n, i + radius + 1)
        neighbor_idx = [j for j in range(lo, hi) if j != i]
        if len(neighbor_idx) < 3:
            continue
        w = x[neighbor_idx].astype(float)
        mask = np.ones(w.shape, dtype=bool)
        for _ in range(max_iters):
            subset = w[mask]
            if subset.size < 2:
                break
            mu, sigma = float(np.mean(subset)), float(np.std(subset))
            if sigma < 1e-12:
                break
            new_mask = np.abs(w - mu) <= n_sigmas * sigma
            if np.array_equal(new_mask, mask):
                break
            mask = new_mask
        subset = w[mask]
        if subset.size < 1:
            continue
        mu, sigma = float(np.mean(subset)), float(np.std(subset))
        if sigma < 1e-12:
            if abs(x[i] - mu) > 1e-9:
                anomalie.append(i)
                out[i] = float(np.median(subset))
            continue
        if abs(x[i] - mu) > n_sigmas * sigma:
            anomalie.append(i)
            out[i] = float(np.median(subset))
    return out, anomalie

pressure_valve

pressure_valve(
    x: ndarray,
    radius: int = 10,
    soglia_pressione: float = 8.0,
    n_sigmas: float = 3.0,
    ref_mult: int = 3,
    k_molla: float = 3.0,
) -> Tuple[np.ndarray, List[int], np.ndarray, np.ndarray]

Orchestratore automatico dei quattro rilevatori sopra: non un voto (quanti metodi segnalano il punto), una combinazione a MINIMA VARIANZA vincolata -- lo stimatore BLUE (Best Linear Unbiased Estimator) classico della statistica, derivato con un moltiplicatore di Lagrange.

Ogni metodo produce una coppia (centro, scala) locale: Chauvenet usa media/std (non robuste, la formulazione originale -- la debolezza e' voluta, vedi sotto), sigma-clipping usa media/std ripulite iterativamente dagli outlier della finestra, Hampel usa mediana/MAD scalato (1.4826*MAD), Tukey usa mediana/IQR scalato (IQR/1.349, la stessa costante che rende IQR uno stimatore di sigma per una gaussiana). Quattro stime indipendenti dello stesso "vero" centro locale, ognuna con la propria incertezza (scala).

Combinarle non con una media semplice ma cercando i pesi w_k che minimizzano la VARIANZA della combinazione pesata Sum(w_k * centro_k), vincolati a Sum(w_k) = 1 (combinazione non distorta) -- un problema di ottimizzazione vincolata, risolto con un moltiplicatore di Lagrange lambda: L(w, lambda) = Sum(w_k^2 * scala_k^2) - lambda * (Sum(w_k) - 1) dL/dw_k = 0 => w_k = lambda / (2 * scala_k^2) dal vincolo => w_k = (1/scala_k^2) / Sum_j(1/scala_j^2) Il risultato (peso inversamente proporzionale al QUADRATO della propria incertezza) e' lo stimatore a minima varianza tra tutte le combinazioni lineari non distorte di stime indipendenti (Gauss-Markov) -- non e' un peso scelto a mano, e' quello che il vincolo impone. Un metodo la cui scala si gonfia (es. Chauvenet quando la finestra contiene gia' un outlier, la sua std non e' robusta) viene automaticamente pesato meno, senza bisogno di scartarlo esplicitamente.

centro_combinato = Sum(w_k * centro_k) scala_combinata = 1 / sqrt(Sum(1/scala_k^2)) (la stessa derivazione da' anche la varianza della combinazione stessa -- SEMPRE piu' stretta di qualunque scala_k singola, e' la proprieta' che rende utile combinare piu' stime indipendenti) pressione = |x[i] - centro_combinato| / scala_combinata

Una singola soglia finale (soglia_pressione) decide -- trigger binario, come il resto di questo ecosistema (core/hybrid_engine.py, dense_evolution/healing.py): mai un terzo stato intermedio, mai un blend tra grezzo e pulito. Non e' piu' un valore "in sigma" nel senso classico (scala_combinata e' piu' stretta di ogni scala_k, quindi la stessa deviazione reale produce una pressione piu' alta che con un singolo metodo) -- soglia_pressione=8.0 e' il default calibrato direttamente (vedi test/test_robust_filters.py): su rumore gaussiano puro (N=300, seed fisso) i percentili della pressione sono 50%=1.4, 95%=4.6, 99%=6.0, 99.9%=7.6 -- 8.0 lascia 1/300 falsi positivi mantenendo comunque un vero outlier (pressione~88) e una coppia di outlier ravvicinati (pressione~40) ben sopra soglia con ampio margine. n_sigmas resta il parametro del sotto-passo sigma-clipping interno (quante sigma oltre cui un vicino viene escluso nell'iterazione, non la soglia finale).

LA MOLLA (JSD): la soglia stessa non e' fissa. Ad ogni punto si confronta la finestra locale (raggio radius, la stessa usata dai 4 metodi) con una finestra di riferimento piu' ampia (raggio radius * ref_mult, default ref_mult=3 -- la stessa costante di healing_filter.py in questo pacchetto) via la divergenza di Jensen-Shannon (_jensen_shannon, 0..1). Se le due distribuzioni sono simili (rumore stazionario, nessuna vera transizione in corso), JSD~0 e la soglia resta soglia_pressione. Se sono molto diverse (una transizione di regime vera e' in corso: la finestra locale, piu' vicina al nuovo livello, non somiglia piu' alla finestra di riferimento, ancora dominata dal vecchio), JSD sale e la soglia si allarga proporzionalmente: soglia_effettiva[i] = soglia_pressione * (1 + k_molla * JSD[i]) "la molla cede" -- lo scudo diventa meno nervoso proprio nei punti dove un cambiamento genuino e' piu' plausibile, senza mai smettere di giudicare in binario (anomalia sopra soglia_effettiva, nient'altro).

_jensen_shannon usa gia' n_bins adattivo e smoothing di Laplace (vedi il suo docstring) apposta per questo uso: su finestre piccole, una JSD "ingenua" (bin fissi, epsilon quasi-zero) risultava PIU' alta su rumore stazionario puro che vicino a una vera transizione -- la molla si sarebbe attivata ovunque, non "solo dove serve". Con la versione corretta, verificato su rumore stazionario a bassa ampiezza vs un vero gradino 1.0->5.0 (radius=10, k_molla=3.0): soglia_effettiva media ~8.3 sul rumore stazionario, ~8.7-8.8 in prossimita' della transizione (piu' alta, direzione giusta) -- comunque un segnale rumoroso per natura (non separa mai perfettamente le due code), ma la soglia puo' solo ALLARGARSI (JSD >= 0), mai restringersi sotto soglia_pressione: anche nel caso peggiore (molla che si attiva "a sproposito" su rumore), il risultato e' uno scudo leggermente piu' permissivo, mai piu' nervoso.

Nessuna configurazione richiesta dal chiamante oltre i default: pensato per un contesto dove chi legge il risultato non deve scegliere pesi o soglie per metodo, solo leggere anomalie/pressione. Verificato (test/test_robust_filters.py): pochi falsi positivi su rumore gaussiano puro, rileva un vero outlier con pressione ben sopra soglia, non tocca un gradino genuino sostenuto.

LIMITE NOTO, ORA RISOLTO -- mascheramento reciproco tra outlier ravvicinati su una baseline quasi piatta (test/test_pressure_valve_vs_ orca_experiment.py, scenario "impulsi alternati" +/-500 su fondo costante 120, zero rumore reale): la finestra locale usata per ogni stima ESCLUDE gia' il punto i stesso (fix applicato qui sotto), ma se altri outlier della stessa esplosione restano nella finestra (qui: raggio 10, 3 picchi entro 2 passi l'uno dall'altro), mediana/IQR e mediana/MAD restano comunque ESATTAMENTE 0 -- non e' un bug della loro implementazione, e' matematicamente corretto: con 18/20 punti della finestra identici, MAD e IQR (breakdown ~50%/~25%) non vedono niente di anomalo. Un M-estimator di Huber (IRLS, vedi Sung et al. 2019, arXiv:1912.04982, usato li' per spettroscopia a 2 qubit con outlier sperimentali) e' stato provato come alternativa e converge ANCH'ESSO a scala 0, per lo stesso motivo di fondo (la sua scala finale e' a sua volta una MAD dei residui). Letto direttamente il paper (non assunto dalla formula): la loro tecnica robusta funziona perche' la scala sigma^2=var(O_tau) viene stimata da MISURE RIPETUTE sullo STESSO punto (rumore di shot, M ripetizioni proiettive), un asse dati ortogonale alla finestra temporale -- non disponibile qui (una singola serie scalare, nessuna ripetizione per punto), quindi la loro soluzione non e' portabile cosi' com'e'.

Provata anche una scala di fallback dalla finestra di riferimento piu' ampia (radius*ref_mult): VERIFICATA NON FUNZIONARE, e per un motivo strutturale, non un errore di implementazione -- allargare la finestra riduce sempre la FRAZIONE di outlier rispetto al totale (3 punti su 21 locali = 14%, 3 su 61 di riferimento = 5%), quindi non puo' mai far superare a IQR/MAD (stimatori a breakdown frazionario fisso) la soglia di rottura che li fa collassare a zero: verificato numericamente, entrambe le scale restavano 0.0 anche sulla finestra di riferimento.

LA VERA CORREZIONE: quando IQR o MAD collassano davvero a zero (finestra locale genuinamente omogenea) ma il punto i si scosta comunque dalla mediana locale, quello NON e' "nessuna informazione" da scartare -- e' la prova piu' forte possibile (deviazione reale contro dispersione locale nulla). Il bug era scartare questa evidenza insieme allo stimatore degenere. Ora: scala-locale-nulla + deviazione reale => pressione forzata a infinito (anomalia certa), bypassando la combinazione pesata per quel punto. Verificato: sullo scenario reale "impulsi alternati" (+/-500 su fondo 120 senza rumore), prima del fix pressione ai 3 picchi = 3.16/5.19/3.16 (sotto soglia 8.0, zero anomalie, RMSE=86.6); dopo il fix, tutti e 3 i picchi rilevati, RMSE=0.0. Nessuna regressione: 223/223 test della suite passano, incluso il tasso di falsi positivi su rumore gaussiano puro (invariato) e tutti gli altri 6 scenari del confronto a tre vie.

Ritorna (pulito, anomalie, pressione, soglia_effettiva): pressione e' la deviazione combinata continua per ogni punto, `soglia_effettiva' la soglia realmente applicata in quel punto (>= soglia_pressione, si vede dove/quanto la molla ha ceduto).

Source code in dense_armor/utility/anomaly/robust_filters.py

def pressure_valve(
    x: np.ndarray, radius: int = 10, soglia_pressione: float = 8.0, n_sigmas: float = 3.0,
    ref_mult: int = 3, k_molla: float = 3.0,
) -> Tuple[np.ndarray, List[int], np.ndarray, np.ndarray]:
    """Orchestratore automatico dei quattro rilevatori sopra: non un voto
    (quanti metodi segnalano il punto), una combinazione a MINIMA VARIANZA
    vincolata -- lo stimatore BLUE (Best Linear Unbiased Estimator) classico
    della statistica, derivato con un moltiplicatore di Lagrange.

    Ogni metodo produce una coppia (centro, scala) locale: Chauvenet usa
    media/std (non robuste, la formulazione originale -- la debolezza e'
    voluta, vedi sotto), sigma-clipping usa media/std ripulite
    iterativamente dagli outlier della finestra, Hampel usa mediana/MAD
    scalato (1.4826*MAD), Tukey usa mediana/IQR scalato (IQR/1.349, la
    stessa costante che rende IQR uno stimatore di sigma per una gaussiana).
    Quattro stime indipendenti dello stesso "vero" centro locale, ognuna con
    la propria incertezza (scala).

    Combinarle non con una media semplice ma cercando i pesi w_k che
    minimizzano la VARIANZA della combinazione pesata Sum(w_k * centro_k),
    vincolati a Sum(w_k) = 1 (combinazione non distorta) -- un problema di
    ottimizzazione vincolata, risolto con un moltiplicatore di Lagrange
    lambda:
        L(w, lambda) = Sum(w_k^2 * scala_k^2) - lambda * (Sum(w_k) - 1)
        dL/dw_k = 0  =>  w_k = lambda / (2 * scala_k^2)
        dal vincolo  =>  w_k = (1/scala_k^2) / Sum_j(1/scala_j^2)
    Il risultato (peso inversamente proporzionale al QUADRATO della propria
    incertezza) e' lo stimatore a minima varianza tra tutte le combinazioni
    lineari non distorte di stime indipendenti (Gauss-Markov) -- non e' un
    peso scelto a mano, e' quello che il vincolo impone. Un metodo la cui
    scala si gonfia (es. Chauvenet quando la finestra contiene gia' un
    outlier, la sua std non e' robusta) viene automaticamente pesato meno,
    senza bisogno di scartarlo esplicitamente.

    centro_combinato = Sum(w_k * centro_k)
    scala_combinata  = 1 / sqrt(Sum(1/scala_k^2))   (la stessa derivazione
                        da' anche la varianza della combinazione stessa --
                        SEMPRE piu' stretta di qualunque scala_k singola,
                        e' la proprieta' che rende utile combinare piu'
                        stime indipendenti)
    pressione        = |x[i] - centro_combinato| / scala_combinata

    Una singola soglia finale (soglia_pressione) decide -- trigger binario,
    come il resto di questo ecosistema (core/hybrid_engine.py,
    dense_evolution/healing.py): mai un terzo stato intermedio, mai un
    blend tra grezzo e pulito. Non e' piu' un valore "in sigma" nel senso
    classico (scala_combinata e' piu' stretta di ogni scala_k, quindi la
    stessa deviazione reale produce una pressione piu' alta che con un
    singolo metodo) -- soglia_pressione=8.0 e' il default calibrato
    direttamente (vedi test/test_robust_filters.py): su rumore gaussiano
    puro (N=300, seed fisso) i percentili della pressione sono
    50%=1.4, 95%=4.6, 99%=6.0, 99.9%=7.6 -- 8.0 lascia 1/300 falsi
    positivi mantenendo comunque un vero outlier (pressione~88) e una
    coppia di outlier ravvicinati (pressione~40) ben sopra soglia con
    ampio margine. `n_sigmas` resta il parametro del sotto-passo
    sigma-clipping interno (quante sigma oltre cui un vicino viene escluso
    nell'iterazione, non la soglia finale).

    LA MOLLA (JSD): la soglia stessa non e' fissa. Ad ogni punto si confronta
    la finestra locale (raggio `radius`, la stessa usata dai 4 metodi) con
    una finestra di riferimento piu' ampia (raggio `radius * ref_mult`,
    default ref_mult=3 -- la stessa costante di healing_filter.py in questo
    pacchetto) via la divergenza di Jensen-Shannon (_jensen_shannon, 0..1).
    Se le due distribuzioni sono simili (rumore stazionario, nessuna vera
    transizione in corso), JSD~0 e la soglia resta soglia_pressione. Se
    sono molto diverse (una transizione di regime vera e' in corso: la
    finestra locale, piu' vicina al nuovo livello, non somiglia piu' alla
    finestra di riferimento, ancora dominata dal vecchio), JSD sale e la
    soglia si allarga proporzionalmente:
        soglia_effettiva[i] = soglia_pressione * (1 + k_molla * JSD[i])
    "la molla cede" -- lo scudo diventa meno nervoso proprio nei punti dove
    un cambiamento genuino e' piu' plausibile, senza mai smettere di
    giudicare in binario (anomalia sopra soglia_effettiva, nient'altro).

    _jensen_shannon usa gia' n_bins adattivo e smoothing di Laplace (vedi il
    suo docstring) apposta per questo uso: su finestre piccole, una JSD
    "ingenua" (bin fissi, epsilon quasi-zero) risultava PIU' alta su rumore
    stazionario puro che vicino a una vera transizione -- la molla si sarebbe
    attivata ovunque, non "solo dove serve". Con la versione corretta,
    verificato su rumore stazionario a bassa ampiezza vs un vero gradino
    1.0->5.0 (radius=10, k_molla=3.0): soglia_effettiva media ~8.3 sul
    rumore stazionario, ~8.7-8.8 in prossimita' della transizione (piu'
    alta, direzione giusta) -- comunque un segnale rumoroso per natura (non
    separa mai perfettamente le due code), ma la soglia puo' solo ALLARGARSI
    (JSD >= 0), mai restringersi sotto soglia_pressione: anche nel caso
    peggiore (molla che si attiva "a sproposito" su rumore), il risultato e'
    uno scudo leggermente piu' permissivo, mai piu' nervoso.

    Nessuna configurazione richiesta dal chiamante oltre i default: pensato
    per un contesto dove chi legge il risultato non deve scegliere pesi o
    soglie per metodo, solo leggere `anomalie`/`pressione`. Verificato
    (test/test_robust_filters.py): pochi falsi positivi su rumore
    gaussiano puro, rileva un vero outlier con pressione ben sopra soglia,
    non tocca un gradino genuino sostenuto.

    LIMITE NOTO, ORA RISOLTO -- mascheramento reciproco tra outlier
    ravvicinati su una baseline quasi piatta (test/test_pressure_valve_vs_
    orca_experiment.py, scenario "impulsi alternati" +/-500 su fondo
    costante 120, zero rumore reale): la finestra locale usata per ogni
    stima ESCLUDE gia' il punto i stesso (fix applicato qui sotto), ma se
    altri outlier della stessa esplosione restano nella finestra (qui:
    raggio 10, 3 picchi entro 2 passi l'uno dall'altro), mediana/IQR e
    mediana/MAD restano comunque ESATTAMENTE 0 -- non e' un bug della loro
    implementazione, e' matematicamente corretto: con 18/20 punti della
    finestra identici, MAD e IQR (breakdown ~50%/~25%) non vedono niente
    di anomalo. Un M-estimator di Huber (IRLS, vedi Sung et al. 2019,
    arXiv:1912.04982, usato li' per spettroscopia a 2 qubit con outlier
    sperimentali) e' stato provato come alternativa e converge ANCH'ESSO a
    scala 0, per lo stesso motivo di fondo (la sua scala finale e' a sua
    volta una MAD dei residui). Letto direttamente il paper (non assunto
    dalla formula): la loro tecnica robusta funziona perche' la scala
    sigma^2=var(O_tau) viene stimata da MISURE RIPETUTE sullo STESSO punto
    (rumore di shot, M ripetizioni proiettive), un asse dati ortogonale
    alla finestra temporale -- non disponibile qui (una singola serie
    scalare, nessuna ripetizione per punto), quindi la loro soluzione non
    e' portabile cosi' com'e'.

    Provata anche una scala di fallback dalla finestra di riferimento piu'
    ampia (radius*ref_mult): VERIFICATA NON FUNZIONARE, e per un motivo
    strutturale, non un errore di implementazione -- allargare la finestra
    riduce sempre la FRAZIONE di outlier rispetto al totale (3 punti su 21
    locali = 14%, 3 su 61 di riferimento = 5%), quindi non puo' mai far
    superare a IQR/MAD (stimatori a breakdown frazionario fisso) la soglia
    di rottura che li fa collassare a zero: verificato numericamente,
    entrambe le scale restavano 0.0 anche sulla finestra di riferimento.

    LA VERA CORREZIONE: quando IQR o MAD collassano davvero a zero
    (finestra locale genuinamente omogenea) ma il punto i si scosta
    comunque dalla mediana locale, quello NON e' "nessuna informazione" da
    scartare -- e' la prova piu' forte possibile (deviazione reale contro
    dispersione locale nulla). Il bug era scartare questa evidenza insieme
    allo stimatore degenere. Ora: scala-locale-nulla + deviazione reale =>
    pressione forzata a infinito (anomalia certa), bypassando la
    combinazione pesata per quel punto. Verificato: sullo scenario reale
    "impulsi alternati" (+/-500 su fondo 120 senza rumore), prima del fix
    pressione ai 3 picchi = 3.16/5.19/3.16 (sotto soglia 8.0, zero
    anomalie, RMSE=86.6); dopo il fix, tutti e 3 i picchi rilevati,
    RMSE=0.0. Nessuna regressione: 223/223 test della suite passano,
    incluso il tasso di falsi positivi su rumore gaussiano puro (invariato)
    e tutti gli altri 6 scenari del confronto a tre vie.

    Ritorna (pulito, anomalie, pressione, soglia_effettiva): `pressione` e'
    la deviazione combinata continua per ogni punto, `soglia_effettiva' la
    soglia realmente applicata in quel punto (>= soglia_pressione, si vede
    dove/quanto la molla ha ceduto)."""
    x = np.asarray(x, dtype=float)
    n = len(x)
    pressione = np.zeros(n)
    soglia_effettiva = np.full(n, soglia_pressione)

    for i in range(n):
        w_self = _window(x, i, radius)
        if w_self.size < 4:
            continue
        # Punto i ESCLUSO dalla propria finestra (come gia' faceva solo il
        # sotto-passo sigma-clipping sotto): includerlo permetteva a un
        # outlier estremo di gonfiare la propria stessa scala di
        # riferimento (mascheramento) -- verificato su impulsi isolati
        # (+/-500 su fondo 120): con la finestra self-inclusive, pressione
        # massima 4.6 (sotto qualunque soglia, zero anomalie rilevate);
        # con la finestra esclusiva, l'outlier non contamina piu' la scala
        # con cui viene giudicato.
        w = np.delete(w_self, min(i - max(0, i - radius), w_self.size - 1))
        w_riferimento = _window(x, i, radius * ref_mult)

        centri, scale = [], []

        mu, sigma = float(np.mean(w)), float(np.std(w))
        if sigma > 1e-12:
            centri.append(mu)
            scale.append(sigma)

        med = float(np.median(w))
        q1, q3 = np.percentile(w, [25, 75])
        iqr_scale = (q3 - q1) / 1.349
        if iqr_scale > 1e-12:
            centri.append(med)
            scale.append(iqr_scale)

        mad = float(np.median(np.abs(w - med)))
        scaled_mad = 1.4826 * mad
        if scaled_mad > 1e-12:
            centri.append(med)
            scale.append(scaled_mad)

        # Scala-zero-e'-evidenza-certa: se IQR o MAD collassano davvero a
        # zero (finestra locale genuinamente omogenea, >=25%/>=50% dei
        # punti identici) MA il punto i si scosta comunque dal centro, non
        # e' "nessuna informazione" da scartare -- e' la prova piu' forte
        # possibile (deviazione reale rispetto a dispersione locale nulla).
        # Il bug del "LIMITE NOTO" originale era proprio scartare questa
        # evidenza insieme allo stimatore degenere. Non tentato: una scala
        # di fallback dalla finestra di riferimento piu' ampia (provato e
        # verificato NON funzionare -- vedi sotto: allargare la finestra
        # riduce sempre la frazione di outlier, non puo' mai far rivivere
        # una statistica a breakdown frazionario). Verificato sullo
        # scenario reale "impulsi alternati" (test/test_pressure_valve_vs_
        # orca_experiment.py, +/-500 su fondo 120 senza rumore): prima di
        # questo fix, pressione ai 3 picchi = 3.16/5.19/3.16 (sotto soglia
        # 8.0, zero anomalie, RMSE=86.6); con questo fix, tutti e 3 i
        # picchi rilevati, RMSE=0.0.
        scala_locale_nulla = (iqr_scale <= 1e-12) or (scaled_mad <= 1e-12)
        if scala_locale_nulla and abs(x[i] - med) > 1e-9:
            pressione[i] = float("inf")
            continue

        neighbor_idx = [j for j in range(max(0, i - radius), min(n, i + radius + 1)) if j != i]
        if len(neighbor_idx) >= 3:
            w_neighbors = x[neighbor_idx].astype(float)
            mask = np.ones(w_neighbors.shape, dtype=bool)
            for _ in range(5):
                subset = w_neighbors[mask]
                if subset.size < 2:
                    break
                mu_s, sigma_s = float(np.mean(subset)), float(np.std(subset))
                if sigma_s < 1e-12:
                    break
                new_mask = np.abs(w_neighbors - mu_s) <= n_sigmas * sigma_s
                if np.array_equal(new_mask, mask):
                    break
                mask = new_mask
            subset = w_neighbors[mask]
            if subset.size >= 1:
                mu_s, sigma_s = float(np.mean(subset)), float(np.std(subset))
                if sigma_s > 1e-12:
                    centri.append(mu_s)
                    scale.append(sigma_s)

        if not centri:
            continue

        centri_arr = np.array(centri)
        scale_arr = np.array(scale)
        pesi = (1.0 / scale_arr ** 2)
        pesi /= pesi.sum()  # w_k = (1/scala_k^2) / Sum_j(1/scala_j^2), dal moltiplicatore di Lagrange

        centro_combinato = float(np.sum(pesi * centri_arr))
        scala_combinata = float(1.0 / np.sqrt(np.sum(1.0 / scale_arr ** 2)))

        if scala_combinata > 1e-12:
            pressione[i] = abs(x[i] - centro_combinato) / scala_combinata

        if w_riferimento.size >= 4 and w.size >= 4:
            jsd = _jensen_shannon(w, w_riferimento)
            soglia_effettiva[i] = soglia_pressione * (1.0 + k_molla * jsd)

    out = np.copy(x).astype(float)
    anomalie = []
    for i in range(n):
        if pressione[i] > soglia_effettiva[i]:
            anomalie.append(i)
            out[i] = float(np.median(_window(x, i, radius)))
    return out, anomalie, pressione, soglia_effettiva

---

Details

• Why the base threshold is 8.0. On pure stationary noise (N = 300, fixed seed) the pressure percentiles are 50% = 1.4, 95% = 4.6, 99% = 6.0, 99.9% = 7.6. A base threshold of 8.0 sits above the 99.9th percentile, so on average fewer than one sample in 300 crosses it. A true outlier reaches a pressure well above the base threshold: the collision in the example above, for instance, has pressure 42.7.

• Why Chauvenet keeps its non-robust scale. The Chauvenet rule is used here in its original 1863 formulation, with the mean and the standard deviation of the window. When the window contains an outlier, that scale inflates, and the minimum-variance weights automatically lower Chauvenet's influence. Replacing Chauvenet's mean and standard deviation with a robust centre and scale would change the rule; keeping the original and letting the weighting handle the outlier is the choice made here.

• Jensen–Shannon on small windows. _jensen_shannon uses adaptive bins and Laplace smoothing. A naive version with fixed bins gave a higher divergence on stationary noise than on a real transition, the opposite of what the threshold formula needs: the threshold would rise on noise and would not rise near a transition. Adaptive bins fix this by sizing the histogram to the range the window actually covers.

• Sources. Pearson, R. K. et al. (2016), "Generalized Hampel filters", EURASIP J. Adv. Signal Process., eqs. 3–4. Chauvenet, W. (1863), A Manual of Spherical and Practical Astronomy, vol. 2. Tukey, J. W. (1977), Exploratory Data Analysis.