Esercizio - Statistica T2 di Hotelling e soglia chi quadro

Questa pagina non ha ancora la versione ripasso: qui sotto c'è il testo completo.

Testo (laboratorio LAB6, anomaly detection).

1. Per una variabile con media $0$ e varianza $1$ (caso univariato) calcolare $T^2$ per $x=5$ e confrontare con la soglia $\chi^2_{1,0{,}95}=3{,}841$.
2. Per due variabili con media $(0,0)$, varianze $1$ e covarianza $0{,}8$ calcolare $T^2$ per $(1,-1)$, $(2,2)$ e $(1{,}5;\,1{,}5)$ con soglia $\chi^2_{2,0{,}95}=5{,}991$; spiegare perché $(1,-1)$ è anomalo e $(2,2)$ no.
3. Implementare la classe HotellingT2 e applicarla al dataset cow (1857 punti nel piano a forma di mucca): contare le anomalie con $\alpha=0{,}05$ e $\alpha=0{,}01$ e commentare i limiti del metodo.

Teoria usata: Anomaly detectionUn'anomalia (outlier) è un'osservazione che si discosta tanto dalle altre da far pensare che sia generata da un meccanismo diverso. Rilevarle serve come pulizia dei dati (solo se sono errori o rumore, non per migliorare artificialmente le metriche), come obiettivo finale (frodi, guasti, cybersicurezza) e per il monitoraggio di un modello in produzione. Metodi semplici: box plot (oltre $1{,}5,\mathrm{IQR}$), carte di controllo univariate ($\mu\pm3\sigma$) e la statistica multivariata di Hotelling $T^2=(x-\bar x)^TS^{-1}(x-\bar x)$ con soglia $\chi^2_{p,1-\alpha}$, valida per dati gaussiani e unimodali. Metodi non supervisionati multivariati danno un anomaly score: l'isolation forest isola ogni punto con split casuali (le anomalie hanno cammini corti) e calcola $s(x,n)=2^{-E(h(x))/c(n)}$, con soglia scelta dalla contaminazione. Senza etichette si valuta con esperti, eventi noti o anomalie sintetiche. Approfondimento: non nel programma di Telecomunicazioni.Anomaly detection → (carte di controllo, $T^2$ di Hotelling); Vettori gaussianiX = (X₁, ..., Xₙ) è un vettore gaussiano N(m, Σ) se ogni combinazione lineare a·X è gaussiana (equivalentemente X = m + AZ con Z gaussiane standard indipendenti); se Σ è invertibile ha densità exp(−½(x−m)ᵀΣ⁻¹(x−m)) / √((2π)ⁿ det Σ). Proprietà chiave: AX + b ~ N(Am + b, AΣAᵀ), le marginali sono gaussiane e componenti non correlate sono indipendenti.Vettori gaussiani → e Distribuzione gammaΓ(α, λ) ha densità λ^α x^(α−1) e^(−λx) / Γ(α) per x > 0, dove Γ(α) = ∫ x^(α−1) e^(−x) dx è la funzione Gamma (Γ(n) = (n − 1)!, Γ(1/2) = √π); media α/λ, varianza α/λ²; Γ(1, λ) = Exp(λ), e la somma di n esponenziali Exp(λ) indipendenti è Γ(n, λ), il tempo d'attesa dell'n-esimo evento.Distribuzione gamma → per la distribuzione chi-quadro; la matrice inversa è in Matrice inversaL'inversa di una matrice quadrata A è la matrice A⁻¹ con A A⁻¹ = A⁻¹ A = I; esiste se e solo se rango(A) = n e si calcola con Gauss-Jordan riducendo (A | I) fino a (I | A⁻¹).Matrice inversa →; l'ellisse è una forma quadratica (Teorema spettrale e forme quadraticheUna funzione lineare è simmetrica se f(v)·w = v·f(w); in una base ortonormale ha matrice simmetrica. Teorema spettrale: f è simmetrica se e solo se esiste una base ortonormale di autovettori, cioè A simmetrica ⇔ PᵀAP diagonale con P ortogonale. Applicato alle forme quadratiche, permette di scriverle come somma di quadrati con gli autovalori come coefficienti.Teorema spettrale e forme quadratiche →).

1. Caso univariato

$T^2=\dfrac{(x-\bar x)^2}{s^2}$. Con un campione di 100 valori $N(0,1)$ (seme 0) e $x=5$: $\bar x\approx0{,}06$ e $s^2\approx1{,}0$ danno $T^2=24{,}03\gg3{,}841$: il punto è un'anomalia. Il valore $3{,}841=1{,}96^2$ è il quadrato del quantile $97{,}5\%$ della normale: in una dimensione la soglia $\chi^2_{1,0{,}95}$ equivale a $|z|>1{,}96$.

2. Caso bidimensionale

Con $S=\begin{pmatrix}1&0{,}8\\0{,}8&1\end{pmatrix}$: determinante $1-0{,}64=0{,}36$ e $S^{-1}=\frac1{0{,}36}\begin{pmatrix}1&-0{,}8\\-0{,}8&1\end{pmatrix}$. Per un punto $x=(a,b)$: $T^2=x^TS^{-1}x=\frac{a^2-1{,}6\,ab+b^2}{0{,}36}$.

Con la distanza euclidea normale $(1,-1)$ ha $\lVert x\rVert^2=2$ e $(2,2)$ ha $8$: sembrerebbe il contrario. La differenza è nella correlazione: le due variabili si muovono insieme ($\rho=0{,}8$), quindi $(2,2)$ è nella direzione di maggiore variabilità (autovalore $1{,}8$ dell'autovettore $(1,1)/\sqrt2$, il punto sta a $\frac{2\sqrt2}{\sqrt{1{,}8}}=2{,}1$ deviazioni standard) mentre $(1,-1)$ è nella direzione di minore variabilità (autovalore $0{,}2$: $\frac{\sqrt2}{\sqrt{0{,}2}}=3{,}2$ deviazioni standard). La regione $T^2\le5{,}991$ è un'ellisse allungata lungo la diagonale (vedi il grafico in Anomaly detectionUn'anomalia (outlier) è un'osservazione che si discosta tanto dalle altre da far pensare che sia generata da un meccanismo diverso. Rilevarle serve come pulizia dei dati (solo se sono errori o rumore, non per migliorare artificialmente le metriche), come obiettivo finale (frodi, guasti, cybersicurezza) e per il monitoraggio di un modello in produzione. Metodi semplici: box plot (oltre $1{,}5,\mathrm{IQR}$), carte di controllo univariate ($\mu\pm3\sigma$) e la statistica multivariata di Hotelling $T^2=(x-\bar x)^TS^{-1}(x-\bar x)$ con soglia $\chi^2_{p,1-\alpha}$, valida per dati gaussiani e unimodali. Metodi non supervisionati multivariati danno un anomaly score: l'isolation forest isola ogni punto con split casuali (le anomalie hanno cammini corti) e calcola $s(x,n)=2^{-E(h(x))/c(n)}$, con soglia scelta dalla contaminazione. Senza etichette si valuta con esperti, eventi noti o anomalie sintetiche. Approfondimento: non nel programma di Telecomunicazioni.Anomaly detection →).

3. Classe HotellingT2 e dataset cow

class HotellingT2:
    def __init__(self, alpha): self.alpha = alpha
    def fit(self, X):
        self.mean = X.mean(axis=0); self.inv_cov = np.linalg.inv(np.cov(X, rowvar=False))
    def t2(self, X):
        d = X - self.mean; return np.einsum("ij,jk,ik->i", d, self.inv_cov, d)
    def predict(self, X):
        return self.t2(X) > chi2.ppf(1 - self.alpha, df=X.shape[1])      # True = outlier

Sul dataset cow: media $(-4{,}40;\ 5{,}65)$, covarianza $\begin{pmatrix}19{,}32&-1{,}72\\-1{,}72&5{,}38\end{pmatrix}$ (correlazione $-0{,}17$). Con $\alpha=0{,}05$ la soglia è $\chi^2_{2,0{,}95}=5{,}991$ e il metodo segnala 46 punti su $1857$ ($2{,}5\%$); con $\alpha=0{,}01$ (soglia $9{,}21$) segnala 0 punti, perché il massimo $T^2$ nel dataset vale $6{,}75$.

Commento: i limiti. (1) Il metodo adatta una gaussiana (un'unica ellisse) ai dati: una nube a forma di mucca non è né gaussiana né unimodale. L'ellisse copre molte zone vuote tra le zampe e la testa: un punto lì sarebbe «normale» per $T^2$ ma è certamente anomalo rispetto alla forma. (2) Per dati gaussiani ci si aspetterebbe il $5\%$ di punti fuori soglia con $\alpha=0{,}05$; qui solo il $2{,}5\%$: i dati sono limitati e con code più leggere. (3) Servono metodi che non assumono una forma, come l'isolation forest (Esercizio - Isolation forest da zero e anomalie sul dataset delle abitazioni).

Verifica

import numpy as np
from scipy.stats import chi2
S = np.array([[1, .8], [.8, 1]]); Si = np.linalg.inv(S)
for p in ([1, -1], [2, 2], [1.5, 1.5]):
    p = np.array(p); print(p, p @ Si @ p)             # 10.0, 4.44, 2.5
print(chi2.ppf([0.95, 0.99], 2))                      # [5.991 9.210]