Esercizio - OLS, R quadro e aumento polinomiale delle feature

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

Testo (laboratorio LAB2). Si generano dati con np.random.seed(0):

1. Due feature. $n=300$ punti con $x_1,x_2$ uniformi in $[0,50]$ (ore di studio e di social), $y=1{,}5x_1-x_2+50+\varepsilon$ con rumore $\varepsilon\sim N(0;0{,}1^2)$. Calcolare $\hat\beta$ con la formula chiusa dell'OLS, MSE e $R^2$; ripetere con il modello sbagliato $\beta=(0,-1,1)$.
2. Una feature. $n=20$ punti $x\sim U(0,1)$, $y=3x^2-2x+1+\varepsilon$ con $\varepsilon\sim N(0;0{,}1^2)$: adattare una retta e una parabola e confrontare $\mathrm{MSE}$ e $R^2$.

Teoria usata: 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, MSE, $R^2$, espansione di base); le equazioni normali sono il Metodo dei minimi quadratiQuando un sistema AX = b non ha soluzioni si cerca X che rende minima la norma di AX − b: AX è la proiezione ortogonale di b su Im A, e X si trova risolvendo le equazioni normali AᵀA X = Aᵀb; applicazione alla retta di regressione.Metodo dei minimi quadrati →, con 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 → e Operazioni tra matriciLe matrici m×n formano uno spazio vettoriale (somma e prodotto per scalare elemento per elemento); il prodotto righe per colonne corrisponde alla composizione di funzioni lineari, è associativo ma non commutativo; la trasposta scambia righe e colonne e (AB)^T = B^T A^T.Operazioni tra matrici →.

1. Regressione con due feature

Matrice e formula. Si aggiunge la colonna di uni: $X=[\mathbf 1,\,x_1,\,x_2]$ ($300\times3$). La soluzione è $\hat\beta=(X^TX)^{-1}X^Ty$ (np.linalg.inv(X.T @ X) @ X.T @ y). Con il seme del laboratorio si trova $$\hat\beta=(50{,}0080;\ 1{,}5001;\ -0{,}9999),$$ quasi uguale ai parametri veri $(50;\ 1{,}5;\ -1)$ (l'errore è dovuto al rumore). Interpretazione: ogni ora di studio in più aumenta il voto stimato di $1{,}5$ punti, ogni ora di social lo riduce di $1$; con zero ore su entrambi il voto atteso è $50$.

MSE. Con le previsioni $\hat y=X\hat\beta$: $\mathrm{MSE}=\frac1n\sum(y_i-\hat y_i)^2=0{,}0084$ (circa la varianza del rumore, $0{,}01$: l'OLS non può fare meglio del rumore).

$R^2$. $SS_{res}=\sum(y_i-\hat y_i)^2=300\cdot0{,}0084=2{,}52$; $SS_{tot}=\sum(y_i-\bar y)^2=198\,878{,}7$; $R^2=1-\frac{2{,}52}{198\,878{,}7}=0{,}999987$: il modello spiega quasi tutta la variabilità.

Modello sbagliato $\beta=(0,-1,1)$: $\hat y=-x_1+x_2$. $\mathrm{MSE}=6007{,}0$ e $R^2=1-\frac{6007\cdot300}{198\,878{,}7}=-8{,}06$: negativo, cioè molto peggio di prevedere sempre la media di $y$ ($R^2=0$).

2. Retta e parabola

Retta ($X=[\mathbf 1,x]$): $\hat\beta=(0{,}458;\ 1{,}071)$, $\mathrm{MSE}=0{,}0801$, $R^2=0{,}521$. La retta non può seguire una curva: la metà della variabilità resta non spiegata (underfitting).

Parabola (espansione di base, $X=[\mathbf 1,x,x^2]$): $\hat\beta=(0{,}9835;\ -2{,}1337;\ 3{,}2289)$ per $(1,\,x,\,x^2)$, cioè $\hat y=0{,}984-2{,}134x+3{,}229x^2$, vicino al vero $1-2x+3x^2$. $\mathrm{MSE}=0{,}0131$ e $R^2=0{,}922$. Il modello resta lineare nei parametri: si risolve con la stessa formula, con una colonna in più.

Con troppi termini. Con grado $9$ l'MSE sul training scende ancora ($0{,}0061$, $R^2=0{,}964$) ma i coefficienti (standardizzati) diventano enormi, dell'ordine di $10^3$-$10^4$: overfitting. L'esercizio sulla ridge (Esercizio - Ridge regression con standardizzazione e coefficienti in unità originali) mostra come ridurlo.

Codice

import numpy as np
np.random.seed(0); n = 300
x1 = np.random.rand(n) * 50; x2 = np.random.rand(n) * 50          # (queste due righe fissano lo stato del generatore)
Xm = np.random.rand(n, 2) * 50
y = Xm @ np.array([1.5, -1.0]) + 50 + np.random.normal(0, 0.1, n)
X = np.hstack([np.ones((n, 1)), Xm])
beta = np.linalg.inv(X.T @ X) @ X.T @ y                           # [50.008, 1.500, -0.9999]
mse = lambda y, p: np.mean((y - p) ** 2)
r2  = lambda y, p: 1 - np.sum((y - p) ** 2) / np.sum((y - y.mean()) ** 2)

(Nel notebook le prime due righe di rand consumano numeri del generatore e poi Xm ne estrae altri: per riprodurre esattamente i valori bisogna rispettare l'ordine delle chiamate.)