Esercizio - PCA da zero sul dataset California Housing

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

Testo (laboratorio LAB1, parte 2). Con il dataset California Housing (20 640 osservazioni, 8 variabili; si toglie MedInc e restano 7 variabili: HouseAge, AveRooms, AveBedrms, Population, AveOccup, Latitude, Longitude):

1. implementare la PCA con il metodo della covarianza usando solo NumPy;
2. trovare il numero $k$ di componenti che spiega almeno il $95\%$ della varianza e disegnare lo scree plot;
3. proiettare i dati sulle prime due componenti;
4. ripetere senza standardizzare e confrontare.

Teoria usata: Analisi delle componenti principali (PCA)Con $p>3$ variabili non si può disegnare il dataset. La PCA (analisi delle componenti principali) lo proietta su pochi assi ortogonali, le componenti principali: dopo aver centrato (e di solito standardizzato) i dati, le direzioni sono gli autovettori della matrice di covarianza $S=\frac1{n-1}X_c^TX_c$ ordinati per autovalore $\lambda_1\ge\lambda_2\ge\dots$; $\lambda_k$ è la varianza lungo la componente $k$ e $\lambda_k/\sum\lambda_j$ la frazione spiegata (scree plot). Trovare la retta che minimizza le distanze dai punti equivale a massimizzare la varianza delle proiezioni (Pitagora). È lineare e conserva la struttura globale, non quella locale; t-SNE e UMAP sono alternative non lineari solo per visualizzare. Approfondimento: non nel programma di Telecomunicazioni.Analisi delle componenti principali (PCA) →; per la matrice di covarianza Covarianza e coefficiente di correlazioneCov(X, Y) = E[(X − E X)(Y − E Y)] = E[XY] − E[X]E[Y] misura quanto X e Y variano insieme; è bilineare, Cov(X, X) = Var(X), Var(X + Y) = Var X + Var Y + 2Cov(X, Y); ρ = Cov / (σ_X σ_Y) sta in [−1, 1] e vale ±1 solo per legami lineari. Indipendenti ⇒ non correlate, ma non viceversa (tranne per i vettori gaussiani).Covarianza e coefficiente di correlazione →; per autovalori e autovettori Autovalori e autovettoriUn autovettore è un vettore non nullo che una funzione lineare manda in un suo multiplo; si trovano gli autovalori come radici del polinomio caratteristico det(A − λI) e gli autovettori come nucleo di A − λI. Matrici simili hanno gli stessi autovalori.Autovalori e autovettori → e 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 →; la standardizzazione è in Statistica per il machine learningI dati di un problema ML si organizzano nella matrice di progetto $X$ ($n$ osservazioni, $p$ variabili). La statistica serve a capirli, ripulirli e prepararli: i momenti (media $\mu$, varianza $\sigma^2$, asimmetria, curtosi), i quartili con lo scarto interquartile $\mathrm{IQR}=Q_3-Q_1$ (all'esame senza interpolazione), la moda per i dati categorici. Con queste quantità si imputano i dati mancanti (media o mediana), si eliminano le variabili costanti e si standardizza con lo z-score $z=(x-\mu)/\sigma$, usando sempre media e deviazione standard del solo training set.Statistica per il machine learning →.

Svolgimento

Passo 1: standardizzare. Per ogni colonna si sottrae la media e si divide per la deviazione standard (calcolata con $n$): Z = (X - X.mean(0)) / X.std(0). Così ogni colonna ha media $0$ e varianza $1$ e nessuna variabile domina per via della scala (Population ha deviazione standard $1132$, AveBedrms $0{,}47$).

Passo 2: covarianza. $S=\frac1{n-1}Z^TZ$, matrice $7\times7$ simmetrica con diagonale $\approx1$ (esattamente $n/(n-1)$ perché la standardizzazione usa $n$). In NumPy: np.cov(Z, rowvar=False) (rowvar=False perché le variabili sono sulle colonne).

Passo 3: autovalori e autovettori di $S$: lam, V = np.linalg.eig(S) (o eigh, che per matrici simmetriche dà autovettori ortonormali già ordinati in modo crescente). Gli autovalori sono

$$2{,}0224;\ 1{,}8333;\ 1{,}2575;\ 1{,}0031;\ 0{,}6698;\ 0{,}1478;\ 0{,}0665,$$

con somma $7{,}0000$ (la varianza totale è $p=7$, perché ogni variabile ha varianza $1$).

Passo 4: ordinare per autovalore decrescente (np.argsort(lam)[::-1]), riordinando anche le colonne di $V$.

Passo 5: varianza spiegata. Frazioni $\lambda_k/7$: $0{,}2889;\ 0{,}2619;\ 0{,}1796;\ 0{,}1433;\ 0{,}0957;\ 0{,}0211;\ 0{,}0095$. Cumulate: $0{,}2889;\ 0{,}5508;\ 0{,}7304;\ 0{,}8737;\ \mathbf{0{,}9694};\ 0{,}9905;\ 1$. Il primo valore $\ge0{,}95$ è il quinto: $k=5$ (con np.argmax(cumulata >= 0.95) + 1). Lo scree plot è il grafico a barre delle frazioni.

Passo 6: proiezione. $Z_2=Z\,V_{[:,\,0:2]}$ ($20640\times2$). Le due colonne hanno varianza $2{,}0224$ e $1{,}8333$ (gli autovalori) e correlazione $0$: le componenti sono incorrelate. Le prime due componenti spiegano il $55{,}1\%$ della varianza: il grafico 2D è una rappresentazione grossolana dei dati.

Che cosa sono le componenti. I loading (autovettori, a meno del segno) delle prime due:

PC1 contrappone latitudine e longitudine (la posizione lungo la costa: Latitude e Longitude hanno correlazione $-0{,}92$, vedi Correlazione e visualizzazione dei datiLa correlazione di Pearson $r=\sum(X_i-\bar X)(Y_i-\bar Y)/\big(\sqrt{\sum(X_i-\bar X)^2}\sqrt{\sum(Y_i-\bar Y)^2}\big)\in[-1,1]$ misura la relazione lineare tra due variabili (covarianza divisa per le deviazioni standard); correlazione non implica causalità. Serve a capire quali variabili contano per il target e a eliminare quelle quasi duplicate (|r| molto alto). Gli indicatori di sintesi non bastano (quartetto di Anscombe, Datasaurus): vanno affiancati ai grafici: istogramma, KDE, box plot, violin plot, heatmap di correlazione, scatter plot e matrice di scatter plot.Correlazione e visualizzazione dei dati →); PC2 è dominata da stanze e camere da letto (correlazione $0{,}85$).

Senza standardizzazione

Si omette il passo 1 e si fa la covarianza dei dati grezzi. La varianza della colonna Population ($1132^2\approx1{,}28\cdot10^6$) è circa $2\cdot10^5$ volte quella di AveRooms ($2{,}47^2=6{,}1$): la PC1 è praticamente l'asse di Population (loading $\approx-1{,}0$ su Population e $\le0{,}004$ in modulo sulle altre) e spiega il $99{,}98\%$ della varianza; PC2 e successive, quasi nulla. Il grafico proiettato è quasi un segmento: la PCA non ha informazione sulle altre variabili. Con la standardizzazione la proiezione mostra una struttura molto più ricca. Il segno degli autovettori è arbitrario: due grafici specchiati sono entrambi corretti.

Codice

import numpy as np
Z = (X - X.mean(0)) / X.std(0)                      # 1. standardizza
S = np.cov(Z, rowvar=False)                         # 2. covarianza 7x7
lam, V = np.linalg.eigh(S)                          # 3. eigh: simmetrica
order = np.argsort(lam)[::-1]; lam, V = lam[order], V[:, order]     # 4. ordine
frac = lam / lam.sum(); k = np.argmax(np.cumsum(frac) >= 0.95) + 1  # 5 -> k = 5
scores = Z @ V[:, :2]                               # 6. proiezione 2D
# senza standardizzare: S = np.cov(X, rowvar=False); frac[0] = 0.9998

Esercizio aggiuntivo del laboratorio

Se si includono tutte e 8 le variabili (con MedInc) gli autovalori standardizzati sono $2{,}027;\ 1{,}881;\ 1{,}271;\ 1{,}031;\ 1{,}003;\ 0{,}659;\ 0{,}082;\ 0{,}046$ (somma $8$) e per il $95\%$ servono $k=6$ componenti (cumulate $0{,}253;\ 0{,}489;\ 0{,}647;\ 0{,}776;\ 0{,}902;\ 0{,}984;\ \dots$).