04 - Distribuzioni multivariate e covarianza

Lanci due dadi e annoti due numeri: il più grande e il più piccolo. Sono due variabili aleatorie definite sullo stesso esperimento, e non sono indipendenti: se il massimo è 1, il minimo non può che essere 1.

Quanto «si muovono insieme»? Come si descrive un fenomeno con più grandezze, come altezza e peso o rendimenti di due titoli, dove conta il legame fra le grandezze oltre alla distribuzione di ciascuna? La risposta ha tre ingredienti: la legge congiunta, che dà le probabilità delle coppie; la covarianza, che misura il legame lineare; e, nel caso continuo, la normale multivariata, in cui tutto si calcola con l'algebra delle matrici.

Qui li costruisci sui due dadi e su un esempio gaussiano, e metti alla prova l'idea più fraintesa: che covarianza zero significhi indipendenza. Ti servono le variabili discrete e le continue, e un minimo di prodotto fra matrici.

Legge congiunta, marginali, condizionate

La legge congiunta di due variabili X e Y assegna una probabilità a ogni coppia di valori: nel caso discreto è la tabella P(X = x, Y = y). Per i dadi i trentasei esiti sono equiprobabili e la tabella si costruisce contando in quanti modi si ottiene ogni coppia (massimo, minimo): una volta se i due dadi coincidono, due volte se sono diversi, perché l'ordine dei lanci non conta.

La tabella ha tutto: sommando una riga ottieni la marginale del massimo, sommando una colonna quella del minimo, perché la marginale di X è P(X = x) = Σ_y P(X = x, Y = y).

Dividendo la congiunta per la marginale ottieni la condizionata, P(X = x | Y = y) = P(X = x, Y = y)/P(Y = y), cioè la probabilità di X quando sai il valore di Y: è la stessa definizione dell'articolo 01, applicata a una coppia.

Due variabili sono indipendenti se la congiunta è il prodotto delle marginali per ogni coppia di valori. Nel caso continuo le somme diventano integrali e la tabella una densità f(x, y), con probabilità calcolate come volumi, ma le definizioni sono le stesse.

Covarianza, correlazione, matrice di covarianza

La covarianza Cov(X, Y) = E[(X − E[X])(Y − E[Y])] = E[XY] − E[X]E[Y] misura se le due variabili si scostano dalla propria media nella stessa direzione (valore positivo), in direzioni opposte (negativo) o senza un legame lineare (zero).

La covarianza di una variabile con sé stessa è la varianza. Dipende dall'unità di misura, per cui si normalizza nel coefficiente di correlazione ρ = Cov(X, Y)/(σ_Xσ_Y), che sta sempre fra −1 e 1 (conseguenza della disuguaglianza di Cauchy-Schwarz) e vale ±1 solo se una variabile è una funzione lineare dell'altra.

La covarianza entra nella varianza di una combinazione lineare: Var(aX + bY) = a²Var(X) + b²Var(Y) + 2ab Cov(X, Y); l'ultimo termine è quello che nell'articolo precedente non compariva, perché per variabili indipendenti vale zero.

Con più variabili le covarianze si raccolgono nella matrice di covarianza Σ, con le varianze sulla diagonale e Cov(X_i, X_j) nella posizione (i, j); è simmetrica e, poiché w^TΣw = Var(w^TX) ≥ 0 per ogni vettore w di coefficienti, è semi-definita positiva.

Per i dadi il calcolo è finito e si fa con le frazioni esatte: marginali, medie, E[XY], covarianza e correlazione, più la media del massimo condizionata al minimo uguale a 3. Poi un milione di coppie simulate stima la covarianza campionaria, con errore standard calcolato dai prodotti centrati stessi; l'asserzione richiede che la stima cada entro quattro errori standard dal valore esatto. Un dettaglio che semplifica il calcolo: poiché massimo per minimo è uguale al prodotto dei due dadi, E[XY] vale 7/2 · 7/2 = 49/4.

from fractions import Fraction as F
from itertools import product
import numpy as np

# Due dadi: X = massimo, Y = minimo. Spazio di probabilita': 36 coppie equiprobabili.
joint = {}
for d1, d2 in product(range(1, 7), repeat=2):
    chiave = (max(d1, d2), min(d1, d2))
    joint[chiave] = joint.get(chiave, 0) + F(1, 36)
assert sum(joint.values()) == 1

pX = {x: sum(p for (a, _), p in joint.items() if a == x) for x in range(1, 7)}   # marginale del massimo
pY = {y: sum(p for (_, b), p in joint.items() if b == y) for y in range(1, 7)}   # marginale del minimo
EX = sum(x * p for x, p in pX.items()); EY = sum(y * p for y, p in pY.items())
EXY = sum(x * y * p for (x, y), p in joint.items())
cov = EXY - EX * EY
varX = sum(x * x * p for x, p in pX.items()) - EX**2
varY = sum(y * y * p for y, p in pY.items()) - EY**2
rho = float(cov) / (float(varX) * float(varY)) ** 0.5
print("E[max] =", EX, " E[min] =", EY, " E[max*min] =", EXY)
print("Cov =", cov, "=", round(float(cov), 4), " Var(max) = Var(min) =", varX, varY == varX, " rho =", round(rho, 4))
print("P(max=1, min=6) =", joint.get((1, 6), 0), " ma P(max=1) P(min=6) =", pX[1] * pY[6], " -> dipendenti")
cond = sum(x * p for (x, y), p in joint.items() if y == 3) / pY[3]
print("E[max | min=3] =", cond, "=", round(float(cond), 4), "contro E[max] =", round(float(EX), 4))

N = 10**6
rng = np.random.default_rng(31)
d = rng.integers(1, 7, size=(N, 2))
x, y = d.max(axis=1), d.min(axis=1)
prodotto = (x - x.mean()) * (y - y.mean())
cov_s, se = prodotto.mean(), prodotto.std() / N**0.5
print(f"Monte Carlo: Cov = {cov_s:.4f} (esatta {float(cov):.4f}, errore standard {se:.4f}), rho = {np.corrcoef(x, y)[0, 1]:.4f}")
assert abs(cov_s - float(cov)) < 4 * se

# Output:
# E[max] = 161/36  E[min] = 91/36  E[max*min] = 49/4
# Cov = 1225/1296 = 0.9452  Var(max) = Var(min) = 2555/1296 True  rho = 0.4795
# P(max=1, min=6) = 0  ma P(max=1) P(min=6) = 1/1296  -> dipendenti
# E[max | min=3] = 33/7 = 4.7143 contro E[max] = 4.4722
# Monte Carlo: Cov = 0.9430 (esatta 0.9452, errore standard 0.0019), rho = 0.4786

La covarianza esatta è 1225/1296, circa 0,9452, e la correlazione 0,4795: massimo e minimo crescono insieme in modo moderato, perché un massimo alto rende possibile un minimo alto. La riga sull'indipendenza mostra la dipendenza nel modo più netto: la coppia (massimo 1, minimo 6) ha probabilità 0, mentre il prodotto delle marginali varrebbe 1/1296, quindi non possono essere indipendenti.

La media condizionata conferma il legame: sapendo che il minimo è 3, il massimo atteso sale da 4,4722 a 33/7 = 4,7143. La simulazione trova 0,9430 con errore standard 0,0019, cioè uno scarto di 1,2 errori standard dal valore esatto, e una correlazione di 0,4786.

Incorrelate non vuol dire indipendenti

Se X e Y sono indipendenti, E[XY] = E[X]E[Y] e la covarianza è zero. Il viceversa è falso, e il controesempio più semplice è una variabile e una sua funzione: X uniforme su {−2, −1, 0, 1, 2} e Y = X².

Qui Y è determinata da X, eppure Cov(X, Y) = E[X³] − E[X]E[X²] = 0, perché i valori di X si bilanciano attorno a zero e il legame (a parabola) non è lineare. L'esercizio 23 della sezione 6.2 di Grinstead e Snell chiede esattamente questo: dimostrare che l'indipendenza implica covarianza nulla e trovare un esempio in cui covarianza nulla e dipendenza coesistono. La covarianza vede solo la componente lineare del legame.

Il secondo caso riguarda la normale, perché si sente spesso dire che per le gaussiane incorrelato significa indipendente. È vero per variabili congiuntamente normali, non per due variabili che sono ciascuna normale. Il programma costruisce X normale standard e Y = S·X, con S = ±1 con probabilità 1/2 indipendente da X. Poiché la normale standard è simmetrica, anche Y è normale standard; ma Y ha sempre lo stesso modulo di X. Il programma stampa la correlazione, le statistiche di Y, la correlazione dei quadrati e due probabilità condizionate.

from fractions import Fraction as F
import numpy as np

# X uniforme su {-2,-1,0,1,2}, Y = X^2: Y e' funzione di X, ma la covarianza e' nulla.
valori = [-2, -1, 0, 1, 2]
p = F(1, 5)
EX = sum(v * p for v in valori)
EY = sum(v * v * p for v in valori)
EXY = sum(v * v * v * p for v in valori)
print("Cov(X, X^2) =", EXY - EX * EY, " con E[X] =", EX, "E[X^2] =", EY)
P_Y4 = sum(p for v in valori if v * v == 4)
print("P(Y=4) =", P_Y4, " ma P(Y=4 | X=2) = 1 -> dipendenti")

# Caso continuo e caso gaussiano: due normali standard, incorrelate, ma non indipendenti
N = 10**6
rng = np.random.default_rng(37)
x = rng.standard_normal(N)
segno = rng.choice([-1.0, 1.0], size=N)       # +/-1 con probabilita' 1/2, indipendente da x
y = segno * x                                   # anche y e' normale standard (la normale e' simmetrica)
print(f"corr(X, Y) = {np.corrcoef(x, y)[0, 1]:+.4f} (errore standard circa {1 / N**0.5:.4f})")
print(f"media e deviazione di Y: {y.mean():+.4f}, {y.std():.4f}   ->  marginale normale standard")
print(f"corr(X^2, Y^2) = {np.corrcoef(x**2, y**2)[0, 1]:.4f}   (Y^2 = X^2 sempre)")
print(f"P(|X|<0.5) = {np.mean(np.abs(x) < .5):.4f}, P(|X|<0.5 | |Y|>2) = {np.mean(np.abs(x[np.abs(y) > 2]) < .5):.4f}")

# Output:
# Cov(X, X^2) = 0  con E[X] = 0 E[X^2] = 2
# P(Y=4) = 2/5  ma P(Y=4 | X=2) = 1 -> dipendenti
# corr(X, Y) = -0.0001 (errore standard circa 0.0010)
# media e deviazione di Y: -0.0016, 0.9987   ->  marginale normale standard
# corr(X^2, Y^2) = 1.0000   (Y^2 = X^2 sempre)
# P(|X|<0.5) = 0.3841, P(|X|<0.5 | |Y|>2) = 0.0000

Con i cinque valori la covarianza è esattamente 0, ma P(Y = 4) = 2/5 mentre P(Y = 4 | X = 2) = 1. Nel caso gaussiano la correlazione simulata fra X e Y è −0,0001, compatibile con zero (l'errore standard è circa 0,001), e Y ha media −0,0016 e deviazione 0,9987, quindi marginale normale standard; ma la correlazione fra i quadrati è 1,0000 e la probabilità che |X| sia sotto 0,5 passa da 0,3841 a 0 se sai che |Y| supera 2.

La coppia (X, Y) non è normale bidimensionale: la sua distribuzione vive su due rette, non su un'ellisse. Conseguenza pratica: un coefficiente di correlazione vicino a zero non dimostra che due grandezze siano indipendenti, e per cercare dipendenze non lineari servono altri strumenti, come l'informazione mutua o i test basati su kernel, che qui non tratto.

La normale multivariata

Un vettore X con d componenti è normale multivariato, X ~ N(μ, Σ), se ha densità proporzionale a exp(−(x − μ)^TΣ⁻¹(x − μ)/2), dove μ è il vettore delle medie e Σ la matrice di covarianza, qui supposta definita positiva (invertibile). La costante di normalizzazione è (2π)^−d/2 det(Σ)^−1/2.

La quantità (x − μ)^TΣ⁻¹(x − μ) è il quadrato della distanza di Mahalanobis da μ: una distanza che tiene conto delle scale e dell'orientamento della nuvola di punti. Se Z ha componenti normali standard indipendenti e Σ = LL^T (fattorizzazione di Cholesky, con L triangolare inferiore), allora X = μ + LZ ha esattamente covarianza Σ: è il modo standard di generare campioni.

L'esempio ha due variabili con deviazioni standard 2 e 1 e correlazione 0,9 (covarianza 1,8), quindi una nuvola molto allungata. Il programma genera 200.000 punti con Cholesky e verifica media e covarianza campionarie.

Poi diagonalizza Σ con numpy.linalg.eigh, adatta alle matrici simmetriche, per trovare le direzioni principali; verifica la varianza di una combinazione lineare con w^TΣw; controlla che il quadrato della distanza di Mahalanobis abbia la distribuzione chi-quadro con 2 gradi di libertà (quindi il 95% dei punti sta sotto il quantile 95% di quella distribuzione); e calcola la regressione di X₁ su X₂, il cui coefficiente è Σ₁₂/Σ₂₂ e la cui varianza residua è Σ₁₁ − Σ₁₂²/Σ₂₂.

import numpy as np
from scipy.stats import chi2

mu = np.array([1.0, 2.0])
Sigma = np.array([[4.0, 1.8],
                  [1.8, 1.0]])                  # sigma1 = 2, sigma2 = 1, rho = 1.8/(2*1) = 0.9
L = np.linalg.cholesky(Sigma)                   # Sigma = L L^T
rng = np.random.default_rng(41)
N = 200_000
Z = rng.standard_normal((N, 2))                  # componenti standard indipendenti
X = mu + Z @ L.T                                 # X = mu + L Z, riga per riga

print("media campionaria:", X.mean(axis=0).round(3), " (vera", mu, ")")
S = np.cov(X, rowvar=False)
print("covarianza campionaria:\n", S.round(3))
assert np.allclose(S, Sigma, atol=0.05) and np.allclose(X.mean(axis=0), mu, atol=0.02)

autoval, autovet = np.linalg.eigh(Sigma)         # eigh: matrici simmetriche, autovalori crescenti
print("autovalori:", autoval.round(4), " quota di varianza della 1a componente principale:", round(autoval[-1] / autoval.sum(), 4))
print("direzione della 1a componente:", np.abs(autovet[:, -1]).round(4))

w = np.array([1.0, -1.0])                        # combinazione X1 - X2
print(f"Var(X1-X2): formula w'Sigma w = {w @ Sigma @ w:.3f}, simulazione {np.var(X @ w):.3f}")

# Distanza di Mahalanobis: d^2 = (x-mu)' Sigma^-1 (x-mu) segue un chi-quadro con 2 gradi di liberta'
D = X - mu
d2 = np.einsum('ij,ij->i', D @ np.linalg.inv(Sigma), D)
soglia = chi2.ppf(0.95, df=2)
oss = np.mean(d2 <= soglia)
tol = 3 * (0.95 * 0.05 / N) ** 0.5
print(f"quantile 95% del chi-quadro(2) = {soglia:.4f}; quota con d^2 <= soglia: {oss:.4f} (atteso 0.95 +/- {tol:.4f})")
assert abs(oss - 0.95) < tol
dist_euclidea = np.sum(D**2, axis=1)
print(f"con la distanza euclidea la stessa soglia lascerebbe dentro il {np.mean(dist_euclidea <= soglia):.4f}")

# Regressione di X1 su X2: coefficiente Sigma12/Sigma22 e varianza residua Sigma11 - Sigma12^2/Sigma22
beta = Sigma[0, 1] / Sigma[1, 1]
res = X[:, 0] - mu[0] - beta * (X[:, 1] - mu[1])
print(f"beta = {beta:.2f}; varianza residua: formula {Sigma[0,0] - Sigma[0,1]**2 / Sigma[1,1]:.3f}, simulazione {res.var():.3f}; corr(residuo, X2) = {np.corrcoef(res, X[:, 1])[0, 1]:+.4f}")

# Output:
# media campionaria: [0.991 1.996]  (vera [1. 2.] )
# covarianza campionaria:
#  [[4.01  1.805]
#  [1.805 1.002]]
# autovalori: [0.1569 4.8431]  quota di varianza della 1a componente principale: 0.9686
# direzione della 1a componente: [0.9056 0.4242]
# Var(X1-X2): formula w'Sigma w = 1.400, simulazione 1.401
# quantile 95% del chi-quadro(2) = 5.9915; quota con d^2 <= soglia: 0.9498 (atteso 0.95 +/- 0.0015)
# con la distanza euclidea la stessa soglia lascerebbe dentro il 0.7272
# beta = 1.80; varianza residua: formula 0.760, simulazione 0.757; corr(residuo, X2) = +0.0017

Gli autovalori sono 4,8431 e 0,1569, e la prima componente principale raccoglie il 96,86% della varianza totale, con direzione (0,9056, 0,4242): la nuvola è quasi una retta, e proiettare su quell'asse perde il 3% dell'informazione. La varianza di X₁ − X₂ è 1,400 dalla formula w^TΣw e 1,401 dalla simulazione.

La soglia del chi-quadro è 5,9915 e la quota di punti entro quella distanza di Mahalanobis è 0,9498, dentro la tolleranza 0,95 ± 0,0015; con la distanza euclidea, che ignora la correlazione, la stessa soglia lascerebbe dentro solo il 72,72% dei punti, e per questo la distanza di Mahalanobis è la scelta giusta per riconoscere i valori anomali di una nuvola correlata.

La regressione dà coefficiente 1,80 e varianza residua 0,757 contro 0,760 teorici, con correlazione fra residuo e X₂ di +0,0017: nel caso gaussiano il residuo è incorrelato (e anzi indipendente) dalla variabile esplicativa, ed è la proprietà che regge la regressione lineare con errori normali.

Che cosa resta fuori

Questo articolo non dimostra le proprietà della normale multivariata (marginali e condizionate gaussiane, la caratterizzazione per combinazioni lineari), che si ricavano con la funzione generatrice dei momenti. Non tratta le copule, che separano le marginali dalla struttura di dipendenza, né l'informazione mutua, né la stima di Σ da pochi dati in dimensione alta, dove la covarianza campionaria diventa singolare o rumorosa.

Quest'ultimo punto è materia di inferenza. Per sapere quanto la media di molte osservazioni dipenda dal legame fra di esse, serve il teorema centrale del limite, che nella versione per vettori ha come limite proprio una normale multivariata.

Domande frequenti

Correlazione e causalità sono la stessa cosa?

No. La correlazione misura un legame lineare osservato fra due grandezze e non dice in che direzione agisce, né se c'è una causa comune. Questo articolo non entra nell'inferenza causale: stabilisce solo che cosa misura il coefficiente.

Perché la matrice di covarianza deve essere semi-definita positiva?

Perché la varianza di qualsiasi combinazione lineare delle variabili, w^TΣw, non può essere negativa. Se vale zero per qualche w non nullo, la matrice non è invertibile e il vettore vive in un sottospazio, senza densità nel senso usuale.

Quando si usa la distanza di Mahalanobis invece della euclidea?

Quando le componenti hanno scale diverse o sono correlate: normalizza per la matrice di covarianza e rende confrontabili le distanze. Nel programma la stessa soglia contiene il 95% dei punti con Mahalanobis e il 72,72% con la distanza euclidea.

Che cosa sai fare adesso

Sai leggere una legge congiunta, ricavarne marginali e condizionate, calcolare covarianza e correlazione con le frazioni esatte e interpretarle per quello che sono, un legame lineare. Sai che covarianza zero non basta per concludere indipendenza, e sai produrre il controesempio, anche gaussiano. Sai generare, diagonalizzare e interrogare una normale multivariata con poche righe di algebra lineare, e perché la distanza di Mahalanobis è il metro giusto in una nuvola correlata. Nel prossimo articolo passi dalle distribuzioni al comportamento delle medie di molte osservazioni: la legge dei grandi numeri e il teorema centrale del limite.

Nella stessa serie trovi anche Analisi dati esplorativa (EDA) con Python. L'elenco completo è nella pagina della serie, Probabilità e Statistica.