Esercizio - Cross-validation k-fold e Monte Carlo per OLS e ridge sul dataset Advertising

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

Testo (esercizio finale del laboratorio LAB2). Con il dataset Advertising (200 mercati; ingressi TV, radio, newspaper = budget pubblicitari, uscita sales = vendite):

• costruire la matrice $X$ e il vettore $y$ come array NumPy;
• scegliere un metodo di cross-validation, addestrare e valutare sia l'OLS sia la ridge regression;
• mostrare le prestazioni (MSE e $R^2$) con box plot e commentare.

Teoria usata: Overfitting, ridge regression e cross-validationUna buona prestazione sul training non basta: serve stimare quella su dati nuovi. La cross-validation (K-fold: $k$ parti, ciascuna a turno come test, errore medio; Monte Carlo: $k$ divisioni casuali con quota di test $q$; leave-one-out se $k=n$) evita di dipendere da una sola divisione casuale. L'errore atteso si scompone in $\text{bias}^2+\text{varianza}+\sigma^2$: i modelli semplici fanno underfitting (bias alto), quelli complessi overfitting (varianza alta). La regolarizzazione aggiunge alla perdita una penalità: la ridge regression minimizza $|y-X\beta|^2+\lambda\sum_{j\ge1}\beta_j^2$ e ha soluzione $\beta=(X^TX+\lambda\tilde I)^{-1}X^Ty$ (l'intercetta non si penalizza, le feature si standardizzano): riduce i coefficienti, rende l'inversa stabile con feature collineari, e $\lambda$ è un iperparametro scelto con la validazione (cross-validation annidata per non contaminare il test).Overfitting, ridge regression e cross-validation → (K-fold e Monte Carlo, ridge); Regressione lineareNell'apprendimento supervisionato si impara una funzione $F(x)$ dagli esempi $(x,y)$: regressione se $y$ è continua, classificazione se è categorica. Il modello lineare è $F_\beta(x)=\beta_0+\beta_1x_1+\dots+\beta_px_p=X\beta$ (con una colonna di uni per $\beta_0$) e i parametri si scelgono minimizzando l'errore quadratico medio $\mathrm{MSE}=\frac1n\sum_i(y_i-F_\beta(x_i))^2$, funzione convessa dei parametri. Annullando il gradiente di $J(\beta)=|y-X\beta|^2$ si ottengono le equazioni normali $X^TX\beta=X^Ty$ e la soluzione dei minimi quadrati ordinari $\beta=(X^TX)^{-1}X^Ty$. Il coefficiente di determinazione $R^2=1-SS_{res}/SS_{tot}$ misura la qualità del fit (0 = come la media, negativo = peggio della media). Un modello va valutato su un test set mai usato per addestrare: l'errore sul training è ottimistico e un polinomio di grado alto lo azzera senza generalizzare.Regressione lineare → (OLS, $R^2$); 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 → (standardizzazione dentro ogni fold); 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 → (matrice di correlazione).

Procedura

Per ogni fold (K-fold, $K=5$, 40 campioni ciascuno).

1. Si mescolano gli indici (np.random.permutation) e si dividono in 5 blocchi da 40.
2. Il blocco $i$ è il test, gli altri quattro (160 campioni) il training.
3. Si standardizzano le feature con media e deviazione standard del solo training e si applicano le stesse al test.
4. Si addestrano OLS ($\hat\beta=(X^TX)^{-1}X^Ty$) e ridge ($\hat\beta=(X^TX+\lambda\tilde I)^{-1}X^Ty$ con $\lambda=10$ e intercetta non penalizzata).
5. Si calcolano MSE e $R^2$ sul fold di test.

Poi si riassumono i 5 valori con un box plot (e la media).

Monte Carlo CV. Si ripete 100 volte: divisione casuale $80\%/20\%$ (160 training, 40 test), con gli stessi passi 3-5.

Risultati (seme 0)

K-fold, MSE per fold (ordine dei fold casuale):

$R^2$ medio: OLS $0{,}891$, ridge $0{,}889$. Monte Carlo ($100$ ripetizioni): MSE medio OLS $3{,}00$ (deviazione standard tra le ripetizioni $1{,}07$), ridge $3{,}07$ ($1{,}13$).

Commento

• Le prestazioni dipendono molto dal fold: l'MSE dell'OLS va da $1{,}46$ a $4{,}40$ (un fattore tre): con una sola divisione train/test si potrebbe pensare di avere un MSE di $1{,}5$ o di $4{,}4$ a seconda della fortuna. La media sui fold e la distribuzione del box plot danno una stima onesta; l'MCCV, con più ripetizioni, mostra lo stesso ordine di grandezza ($3{,}0$).
• La ridge non aiuta qui. Con $\lambda=10$ l'MSE è leggermente peggiore ($3{,}00$ contro $2{,}92$). Motivo: la ridge serve quando le feature sono collineari o troppo numerose rispetto ai dati. Qui sono 3 e quasi incorrelate: le correlazioni sono $0{,}055$ (TV-radio), $0{,}057$ (TV-newspaper), $0{,}354$ (radio-newspaper). Non c'è multicollinearità e l'OLS ha poca varianza. Inoltre la ridge introduce un bias.
• Coefficienti standardizzati (su tutti i dati, OLS): $\beta_0=14{,}02$, TV $3{,}92$, radio $2{,}79$, newspaper $-0{,}02$. Il budget sui giornali non ha effetto (coefficiente $\approx0$). Con la ridge $\lambda=10$: $3{,}74;\ 2{,}65;\ 0{,}04$; con $\lambda=100$: $2{,}64;\ 1{,}84;\ 0{,}26$ (si riducono tutti, newspaper compreso, ma non si annulla: per azzerarlo servirebbe il LASSO, LASSO e discesa del gradienteIl LASSO è la regressione regolarizzata con penalità $L_1$: minimizza $\sum_i(y_i-x_i^T\beta)^2+\lambda\sum_{j\ge1}|\beta_j|$. A differenza della ridge, porta alcuni coefficienti esattamente a zero (soluzione sparsa, selezione delle feature): geometricamente le curve di livello dell'errore toccano il vincolo $\sum|\beta_j|\le s$ (un rombo) in uno spigolo. Non ha formula chiusa, quindi si minimizza con la discesa del gradiente $W\leftarrow W-\eta,\nabla J(W)$, usando il subgradiente $\operatorname{sign}(\beta_j)$ per il valore assoluto (nel punto 0 qualunque valore in $[-1,1]$). Il passo $\eta$ è critico: per l'errore quadratico converge se $\eta<1/\mu_{\max}(X^TX)$; si può usare $\eta_t=\eta_0/(1+\gamma t)$. L'Elastic Net combina le penalità $L_1$ e $L_2$ con $\lambda_1=\alpha\lambda$, $\lambda_2=(1-\alpha)\lambda$.LASSO e discesa del gradiente →). La ridge ha un effetto sull'MSE solo per $\lambda$ molto grande.
• Per il confronto si usano le stesse partizioni per i due modelli (confronto appaiato), e ogni modello è valutato sui dati che non ha visto.

Codice

import numpy as np
A = np.loadtxt("Advertising.csv", delimiter=",", skiprows=1, usecols=(1, 2, 3, 4))
X, y = A[:, :3], A[:, 3]

def fit_ridge(X, y, lam):                          # lam = 0 -> OLS
    X = np.hstack([np.ones((len(X), 1)), X]); I = np.eye(X.shape[1]); I[0, 0] = 0
    return np.linalg.solve(X.T @ X + lam * I, X.T @ y)

idx = np.random.RandomState(0).permutation(len(X))
mse = {0: [], 10: []}
for fold in np.array_split(idx, 5):
    tr = np.setdiff1d(idx, fold)
    mu, sd = X[tr].mean(0), X[tr].std(0)           # statistiche del solo training
    Xtr, Xte = (X[tr] - mu) / sd, (X[fold] - mu) / sd
    for lam in mse:
        b = fit_ridge(Xtr, y[tr], lam)
        mse[lam].append(np.mean((y[fold] - (b[0] + Xte @ b[1:])) ** 2))
# box plot: plt.boxplot([mse[0], mse[10]], tick_labels=["OLS", "Ridge"])

La scelta del seed cambia i singoli valori ma non la conclusione.