Distanza di Mahalanobis: Guida Pratica con Esercizi Svolti e Codice Python

In sintesi: la distanza di Mahalanobis misura quanto un’osservazione sia distante dal centro di una distribuzione tenendo conto delle scale delle variabili e delle loro correlazioni. È una generalizzazione multivariata dello z-score, ma non è automaticamente una probabilità di anomalia né una prova che un dato sia errato.

La distanza di Mahalanobis è uno degli strumenti più utili per analizzare dati multivariati. Può essere impiegata nel controllo qualità, nel monitoraggio dei sensori, nella data quality, nella finanza e nel rilevamento delle anomalie.

Il suo vantaggio rispetto alla distanza euclidea è che non tratta tutte le direzioni dello spazio come equivalenti. Una variazione lungo una direzione molto variabile può essere normale, mentre una variazione più piccola lungo una direzione stabile può essere statisticamente rilevante.

In questo articolo partiamo dall’intuizione e dalla formula, poi analizziamo tre esercizi progressivi con codice Python:

1. un caso unidimensionale, equivalente allo z-score;
2. un caso bidimensionale con variabili scorrelate;
3. un caso bidimensionale con covarianza non diagonale.

Alla fine discutiamo anche le soglie statistiche, i limiti del metodo e le alternative robuste.

Indice

• Che cos’è la distanza di Mahalanobis
• Distanza euclidea e distanza statistica
• La formula
• Il ruolo della matrice di covarianza
• Quando il quadrato della distanza segue una chi-quadro
• Esercizio 1: il caso unidimensionale
• Esercizio 2: due sensori scorrelati
• Esercizio 3: due rendimenti correlati
• Implementazione numericamente stabile
• Limiti pratici e covarianza robusta
• Conclusioni
• Riferimenti

Che cos’è la distanza di Mahalanobis

Supponiamo di osservare un vettore di dati

[math]\displaystyle x = (x_1, x_2, \ldots, x_p)^\top [/math]

con [math]p[/math] variabili. Vogliamo capire se l’osservazione sia compatibile con una popolazione caratterizzata da:

• un vettore medio [math]\mu[/math];
• una matrice di covarianza [math]\Sigma[/math].

La distanza di Mahalanobis risponde alla domanda:

Quanto è grande lo scostamento dell’osservazione rispetto alla scala e alla struttura di variabilità della distribuzione?

La parola multivariate è essenziale. Il metodo non analizza ogni variabile isolatamente, ma considera simultaneamente le relazioni tra tutte le variabili.

Un esempio intuitivo

Consideriamo altezza e peso. Una persona alta e pesante può essere perfettamente coerente con la popolazione se altezza e peso crescono insieme. Una persona molto alta e insolitamente leggera può invece rappresentare una combinazione più rara, anche se la sua distanza euclidea dal valore medio non è necessariamente maggiore.

La distanza euclidea considera la posizione geometrica del punto. La distanza di Mahalanobis considera anche la forma e l’orientamento della nuvola dei dati.

Distanza euclidea e distanza statistica

La distanza euclidea tra un punto [math]x[/math] e il centro [math]\mu[/math] è:

[math]\displaystyle d_E(x,\mu) = \sqrt{(x-\mu)^\top(x-\mu)}. [/math]

Questa definizione funziona bene quando le variabili hanno scale confrontabili e quando le direzioni dello spazio possono essere considerate equivalenti.

Nei dati reali, queste condizioni spesso non sono soddisfatte.

La distanza di Mahalanobis corregge questi problemi attraverso la matrice di covarianza.

La formula

La distanza di Mahalanobis tra [math]x[/math] e [math]\mu[/math] è definita come:

[math]\displaystyle D_M(x,\mu) = \sqrt{(x-\mu)^\top \Sigma^{-1}(x-\mu)}. [/math]

Il quadrato della distanza è:

[math]\displaystyle D_M^2(x,\mu) = (x-\mu)^\top \Sigma^{-1}(x-\mu). [/math]

I simboli hanno il seguente significato:

La matrice inversa assegna un peso diverso alle varie direzioni. Le direzioni in cui la distribuzione è molto variabile vengono penalizzate meno. Le direzioni in cui la distribuzione è concentrata vengono penalizzate di più.

Caso di covarianza diagonale

Se le variabili sono scorrelate, la matrice di covarianza è diagonale:

[math]\displaystyle \Sigma = \operatorname{diag}(\sigma_1^2,\sigma_2^2,\ldots,\sigma_p^2). [/math]

In questo caso:

[math]\displaystyle D_M^2 = \sum_{i=1}^{p}\frac{(x_i-\mu_i)^2}{\sigma_i^2}. [/math]

Definendo gli z-score

[math]\displaystyle z_i = \frac{x_i-\mu_i}{\sigma_i}, [/math]

otteniamo:

[math]\displaystyle D_M^2 = \sum_{i=1}^{p}z_i^2. [/math]

La distanza di Mahalanobis è quindi la radice della somma dei quadrati degli z-score quando le variabili sono scorrelate.

Il ruolo della matrice di covarianza

La matrice di covarianza contiene le informazioni sulla variabilità e sulle relazioni tra le variabili:

[math]\displaystyle \Sigma = \begin{pmatrix} \sigma_1^2 & \sigma_{12} & \cdots \\ \sigma_{21} & \sigma_2^2 & \cdots \\ \vdots & \vdots & \ddots \end{pmatrix}. [/math]

Gli elementi sulla diagonale sono le varianze. Gli elementi fuori diagonale sono le covarianze.

Se due variabili sono positivamente correlate, la nuvola dei dati tende a essere orientata verso la direzione [math](1,1)[/math]. Se sono negativamente correlate, tende a essere orientata verso la direzione [math](1,-1)[/math].

Nel caso bidimensionale, le superfici a distanza di Mahalanobis costante sono ellissi. In dimensioni superiori diventano ellissoidi.

La decomposizione spettrale rende esplicita questa geometria:

[math]\displaystyle \Sigma = V\Lambda V^\top, [/math]

dove le colonne di [math]V[/math] sono gli autovettori e [math]\Lambda[/math] contiene gli autovalori.

Se [math]d=x-\mu[/math], allora:

[math]\displaystyle D_M^2 = \sum_{i=1}^{p}\frac{(d^\top v_i)^2}{\lambda_i}. [/math]

Ogni componente dello scarto viene divisa per la varianza lungo la direzione corrispondente. Un autovalore piccolo produce un peso grande; un autovalore grande produce un peso più contenuto.

Quando il quadrato della distanza segue una chi-quadro

Se [math]x[/math] proviene da una distribuzione normale multivariata con media [math]\mu[/math], matrice di covarianza [math]\Sigma[/math] e parametri noti, allora:

[math]\displaystyle D_M^2 \sim \chi^2_p, [/math]

dove [math]p[/math] è il numero di variabili.

Questo risultato permette di confrontare [math]D_M^2[/math] con un quantile della distribuzione chi-quadro. Per esempio, con due variabili, il quantile al 95% è circa:

[math]\displaystyle \chi^2_{2,0.95} \approx 5.991. [/math]

La regola decisionale sarebbe:

[math]\displaystyle D_M^2 > \chi^2_{p,1-\alpha} [/math]

per segnalare un’osservazione estrema al livello [math]\alpha[/math], sotto le ipotesi del modello.

Questa soglia non deve però essere interpretata come una garanzia assoluta. Nella pratica:

• media e covarianza sono spesso stimate dal campione;
• i dati possono non essere gaussiani;
• gli outlier possono distorcere la covarianza;
• la matrice può essere quasi singolare;
• la dipendenza temporale può invalidare l’ipotesi di osservazioni indipendenti.

Per questo è più corretto parlare di segnalazione statistica o osservazione incompatibile con il modello, non di errore certo.

Esercizio 1: il caso unidimensionale

Problema

Un’azienda produce aste con lunghezza nominale di 50 mm. Il processo ha media:

[math]\displaystyle \mu=50.2\text{ mm} [/math]

ed errore standard:

[math]\displaystyle \sigma=0.3\text{ mm}. [/math]

Un’asta misura:

[math]\displaystyle x=51.1\text{ mm}. [/math]

Il pezzo viene scartato quando la distanza di Mahalanobis è almeno 3.

Soluzione

In una dimensione la matrice di covarianza è lo scalare [math]\sigma^2[/math]. Pertanto:

[math]\displaystyle D_M(x)=\frac{|x-\mu|}{\sigma}. [/math]

Lo scarto dalla media è:

[math]\displaystyle x-\mu=51.1-50.2=0.9\text{ mm}. [/math]

La distanza è:

[math]\displaystyle D_M=\frac{0.9}{0.3}=3. [/math]

Il valore coincide con la soglia. Se il criterio aziendale è “scartare quando [math]D_M\ge3[/math]”, l’asta deve essere scartata.

Se la stessa misura fosse osservata su una linea con [math]\sigma=0.1[/math] mm, avremmo:

[math]\displaystyle D_M=\frac{0.9}{0.1}=9. [/math]

La stessa misura è quindi molto più anomala rispetto a un processo più preciso.

Codice Python

from __future__ import annotations


def mahalanobis_1d(x: float, mu: float, sigma: float) -> float:
    """Calcola la distanza di Mahalanobis in una dimensione."""
    if sigma <= 0:
        raise ValueError("sigma deve essere positiva")
    return abs(x - mu) / sigma


mu = 50.2
x = 51.1

for sigma in (0.3, 0.1):
    distance = mahalanobis_1d(x, mu, sigma)
    print(f"sigma={sigma:.1f} -> D_M={distance:.1f}")

sigma=0.3 -> D_M=3.0
sigma=0.1 -> D_M=9.0

In una dimensione la distanza di Mahalanobis non aggiunge informazione rispetto allo z-score assoluto. La sua utilità emerge quando le variabili sono multiple e possono essere correlate.

Esercizio 2: due sensori scorrelati

Problema

Un impianto industriale misura temperatura e pressione. I dati storici sono descritti da:

[math]\displaystyle \mu= \begin{pmatrix} 25.0\\ 1.5 \end{pmatrix} [/math]

ed:

[math]\displaystyle \Sigma= \begin{pmatrix} 4.0 & 0\\ 0 & 0.25 \end{pmatrix}. [/math]

Una nuova lettura è:

[math]\displaystyle x= \begin{pmatrix} 28.0\\ 1.0 \end{pmatrix}. [/math]

Verifichiamo se il punto appartiene alla regione teorica al 95%.

Soluzione

Lo scarto è:

[math]\displaystyle d=x-\mu= \begin{pmatrix} 3.0\\ -0.5 \end{pmatrix}. [/math]

Le deviazioni standard sono:

[math]\displaystyle \sigma_T=2, \qquad \sigma_P=0.5. [/math]

Gli z-score sono quindi:

[math]\displaystyle z_T=\frac{3}{2}=1.5, \qquad z_P=\frac{-0.5}{0.5}=-1. [/math]

Poiché la covarianza è diagonale:

[math]\displaystyle D_M^2=1.5^2+(-1)^2=2.25+1=3.25. [/math]

Di conseguenza:

[math]\displaystyle D_M=\sqrt{3.25}\approx1.803. [/math]

Con due variabili, la soglia teorica al 95% è circa 5.991. Poiché:

[math]\displaystyle 3.25<5.991, [/math]

la lettura è interna alla regione ellittica teorica al 95%, assumendo un modello gaussiano adeguato.

Codice Python

from __future__ import annotations

import numpy as np
from scipy.stats import chi2


def mahalanobis(x: np.ndarray, mu: np.ndarray, sigma: np.ndarray) -> float:
    """Calcola D_M usando un sistema lineare, senza invertire Sigma."""
    x = np.asarray(x, dtype=float)
    mu = np.asarray(mu, dtype=float)
    sigma = np.asarray(sigma, dtype=float)

    d = x - mu
    solution = np.linalg.solve(sigma, d)
    squared_distance = d @ solution

    if squared_distance < 0:
        raise ValueError("La matrice di covarianza non è definita positiva")

    return float(np.sqrt(squared_distance))


mu = np.array([25.0, 1.5])
sigma = np.diag([4.0, 0.25])
x = np.array([28.0, 1.0])

D = mahalanobis(x, mu, sigma)
D2 = D**2
threshold = chi2.ppf(0.95, df=2)

print(f"D_M^2 = {D2:.3f}")
print(f"D_M   = {D:.3f}")
print(f"Soglia 95% = {threshold:.3f}")
print(f"Dentro la regione? {D2 <= threshold}")

D_M^2 = 3.250
D_M   = 1.803
Soglia 95% = 5.991
Dentro la regione? True

L’esempio mostra perché una variabile con bassa variabilità può contribuire molto alla distanza: lo scarto viene valutato rispetto alla sua scala naturale.

Esercizio 3: due rendimenti correlati

Problema

Consideriamo i rendimenti giornalieri, espressi in percentuale, di due titoli finanziari. La media è:

[math]\displaystyle \mu= \begin{pmatrix} 0.10\\ 0.05 \end{pmatrix} [/math]

La matrice di covarianza è:

[math]\displaystyle \Sigma= \begin{pmatrix} 4 & 3\\ 3 & 9 \end{pmatrix}. [/math]

La covarianza positiva indica che i due titoli tendono a muoversi nella stessa direzione.

Osserviamo:

[math]\displaystyle x= \begin{pmatrix} 2.10\\ -0.95 \end{pmatrix}. [/math]

Calcoliamo la distanza di Mahalanobis e confrontiamola con la distanza euclidea.

Calcolo diretto

Lo scarto è:

[math]\displaystyle d=x-\mu= \begin{pmatrix} 2\\ -1 \end{pmatrix}. [/math]

Il determinante della matrice è:

[math]\displaystyle \det(\Sigma)=4\cdot9-3^2=27. [/math]

L’inversa è:

[math]\displaystyle \Sigma^{-1}=\frac{1}{27} \begin{pmatrix} 9 & -3\\ -3 & 4 \end{pmatrix}. [/math]

Calcoliamo:

[math]\displaystyle \Sigma^{-1}d =\frac{1}{27} \begin{pmatrix} 21\\ -10 \end{pmatrix}. [/math]

Quindi:

[math]\displaystyle D_M^2=d^\top\Sigma^{-1}d =\frac{52}{27} \approx1.926. [/math]

La distanza è:

[math]\displaystyle D_M\approx1.388. [/math]

La distanza euclidea è invece:

[math]\displaystyle \|d\|_2=\sqrt{2^2+(-1)^2}=\sqrt{5}\approx2.236. [/math]

In questo esempio:

[math]\displaystyle D_M\approx1.388<2.236\approx\|d\|_2. [/math]

Il punto non è estremo secondo la soglia chi-quadro al 95% perché:

[math]\displaystyle D_M^2\approx1.926<5.991. [/math]

Interpretazione spettrale: perché Mahalanobis comprime la distanza?

La correlazione positiva orienta l’ellisse verso la direzione [math](1,1)[/math]. Lo scarto osservato è invece circa [math](2,-1)[/math], quindi è molto vicino alla direzione discordante associata all’asse di minima varianza.

Gli autovalori di [math]\Sigma[/math] sono:

[math]\displaystyle \lambda_1=\frac{13+\sqrt{61}}{2}\approx10.405, [/math]

[math]\displaystyle \lambda_2=\frac{13-\sqrt{61}}{2}\approx2.595. [/math]

La deviazione standard lungo la direzione associata all’autovalore minore è quindi:

[math]\displaystyle \sqrt{\lambda_2}\approx1.611. [/math]

Il vettore [math]d=(2,-1)[/math] non coincide però esattamente con l’autovettore associato a [math]\lambda_2[/math]. Infatti:

[math]\displaystyle \Sigma d= \begin{pmatrix} 4 & 3\\ 3 & 9 \end{pmatrix} \begin{pmatrix} 2\\ -1 \end{pmatrix} = \begin{pmatrix} 5\\ -3 \end{pmatrix}, [/math]

e il risultato non è proporzionale a [math](2,-1)[/math]. Perciò non è rigoroso scrivere che il rapporto tra le due distanze sia esattamente [math]1/\sqrt{\lambda_2}[/math].

Il rapporto effettivo è:

[math]\displaystyle \frac{D_M}{\|d\|_2} =\frac{1.3878}{2.2361} \approx0.6207. [/math]

Il valore è molto vicino a [math]1/\sqrt{\lambda_2}\approx0.6207[/math] perché lo scarto è quasi allineato con la direzione di minima varianza. La decomposizione spettrale permette quindi di capire sia quanto il punto sia lontano dal centro sia in quali direzioni si stia allontanando.

La distanza viene dunque compressa rispetto alla distanza euclidea, ma non perché il vettore coincida esattamente con un autovettore.

La correlazione può comprimere o amplificare la distanza

Manteniamo le deviazioni standard [math]\sigma_1=2[/math] e [math]\sigma_2=3[/math], ma scriviamo la matrice come:

[math]\displaystyle \Sigma(\rho)= \begin{pmatrix} 4 & 6\rho\\ 6\rho & 9 \end{pmatrix}. [/math]

È importante distinguere due domande diverse. La prima riguarda la geometria complessiva della distribuzione. Il valore di correlazione per cui l’autovalore minore diventa 1 soddisfa:

[math]\displaystyle \lambda_2(\rho)=1. [/math]

La soluzione è:

[math]\displaystyle \rho^*=\frac{2}{3}\approx0.667. [/math]

Questo valore descrive il cambiamento della varianza lungo la direzione più stretta dell’ellisse. Non è però automaticamente la soglia alla quale la distanza di Mahalanobis coincide con la distanza euclidea per ogni vettore.

La seconda domanda riguarda lo specifico scarto [math]d=(2,-1)[/math]. Per trovare quando le due distanze coincidono bisogna imporre direttamente:

[math]\displaystyle D_M^2=d^\top d=5. [/math]

Per la matrice

[math]\displaystyle \Sigma(\rho)= \begin{pmatrix} 4 & 6\rho\\ 6\rho & 9 \end{pmatrix} [/math]

si ottiene:

[math]\displaystyle D_M^2=\frac{2(3\rho+5)}{9(1-\rho^2)}. [/math]

Ponendo questa quantità uguale a 5, la soluzione positiva è:

[math]\displaystyle \rho^*=\frac{4\sqrt{11}-1}{15}\approx0.818. [/math]

Per il particolare scarto considerato, quindi, Mahalanobis è inferiore alla distanza euclidea fino a circa [math]\rho=0.818[/math], coincide con essa intorno a questo valore e diventa superiore per correlazioni maggiori.

La distinzione è fondamentale: la soglia [math]0.667[/math] riguarda un autovalore della matrice e descrive la geometria complessiva della distribuzione; la soglia [math]0.818[/math] riguarda invece il confronto tra due distanze calcolate sul vettore specifico [math]d=(2,-1)[/math].

Perciò:

Il valore [math]0.667[/math] resta corretto per l’analisi dell’autovalore minore. Non deve però essere usato come soglia di equivalenza tra le due distanze per il punto specifico. Il valore corretto in quel caso è circa [math]0.818[/math].

Codice Python

from __future__ import annotations

import numpy as np
from scipy.stats import chi2


def mahalanobis_squared(
    x: np.ndarray,
    mu: np.ndarray,
    sigma: np.ndarray,
) -> float:
    """Restituisce il quadrato della distanza di Mahalanobis."""
    d = np.asarray(x, dtype=float) - np.asarray(mu, dtype=float)
    sigma = np.asarray(sigma, dtype=float)
    return float(d @ np.linalg.solve(sigma, d))


mu = np.array([0.10, 0.05])
sigma = np.array([[4.0, 3.0], [3.0, 9.0]])
x = np.array([2.10, -0.95])

D2 = mahalanobis_squared(x, mu, sigma)
D = np.sqrt(D2)
euclidean = np.linalg.norm(x - mu)
threshold = chi2.ppf(0.95, df=2)

print(f"D_M^2      = {D2:.4f}")
print(f"D_M        = {D:.4f}")
print(f"D euclidea  = {euclidean:.4f}")
print(f"Soglia 95% = {threshold:.4f}")
print(f"Segnalazione? {D2 > threshold}")

D_M^2      = 1.9259
D_M        = 1.3878
D euclidea  = 2.2361
Soglia 95% = 5.9915
Segnalazione? False

Implementazione numericamente stabile

Dal punto di vista numerico, è preferibile evitare l’inversione esplicita della matrice:

np.linalg.inv(sigma) @ d

È meglio risolvere il sistema lineare:

np.linalg.solve(sigma, d)

Questa scelta è generalmente più stabile e più efficiente.

Quando la matrice è simmetrica e definita positiva, una fattorizzazione di Cholesky può essere ancora più adatta. Se la matrice è quasi singolare, è necessario valutare una regolarizzazione o uno stimatore di covarianza più robusto.

Limiti pratici e covarianza robusta

La distanza di Mahalanobis classica presuppone che la media e la covarianza descrivano adeguatamente il nucleo della distribuzione. Se il dataset contiene molti outlier, la media può spostarsi e la covarianza può gonfiarsi. Il risultato è paradossale: gli outlier possono rendere se stessi meno riconoscibili.

In questi casi si possono considerare:

• covarianza robusta basata su Minimum Covariance Determinant;
• MinCovDet di scikit-learn;
• EllipticEnvelope per il rilevamento di osservazioni anomale;
• shrinkage covariance, come Ledoit-Wolf, in presenza di molte variabili;
• trasformazioni o modelli non gaussiani quando la distribuzione è fortemente asimmetrica.

La documentazione di scikit-learn distingue inoltre tra:

• outlier detection, quando il dataset di addestramento può contenere osservazioni anomale;
• novelty detection, quando si dispone di dati storici considerati rappresentativi e si devono valutare nuove osservazioni.

La scelta del metodo deve dipendere dal processo generativo dei dati, non soltanto dalla disponibilità di una formula.

Una procedura pratica

Per applicare Mahalanobis a un problema reale, è utile seguire questi passaggi:

1. definire chiaramente le variabili e il periodo di riferimento;
2. controllare valori mancanti, errori di misura e scale incompatibili;
3. stimare media e covarianza sul dataset appropriato;
4. verificare il condizionamento della matrice di covarianza;
5. scegliere, se necessario, una stima robusta o regolarizzata;
6. calcolare [math]D_M^2[/math], non soltanto [math]D_M[/math], quando si usano soglie chi-quadro;
7. verificare le ipotesi di normalità o utilizzare soglie empiriche;
8. interpretare l’alert nel contesto del processo;
9. validare i falsi positivi e i falsi negativi;
10. documentare la soglia e la versione dei dati usati per stimare il modello.

Una distanza elevata deve generare un controllo, non una decisione automatica priva di contesto.

Conclusioni

La distanza di Mahalanobis è una distanza adattata alla struttura statistica dei dati.

Rispetto alla distanza euclidea:

• tiene conto delle diverse scale;
• incorpora le correlazioni;
• pesa ogni direzione secondo la sua variabilità;
• permette di costruire regioni ellittiche o ellissoidali;
• può essere utilizzata per individuare osservazioni incompatibili con un modello.

In una dimensione coincide con lo z-score assoluto. Con variabili scorrelate è la radice della somma dei quadrati degli z-score. Con variabili correlate, invece, la matrice di covarianza modifica anche l’orientamento della geometria.