Esercizio - Regressione ai minimi quadrati su quattro punti

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

Testo. Dati i punti $(x_i,y_i)$: $(0,1)$, $(1,3)$, $(2,4)$, $(3,8)$.

1. Trovare la retta $\hat y=\beta_0+\beta_1x$ che minimizza l'errore quadratico con le equazioni normali.
2. Calcolare residui, $\mathrm{MSE}$, $\mathrm{RMSE}$ e $R^2$, e verificare che i residui sono ortogonali alle colonne di $X$.
3. Prevedere $y$ per $x=5$ e confrontare con la regressione di un modello quadratico.

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 →; stesso problema del 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 → (Esercizio 87 · retta dei minimi quadrati per quattro punti è il caso analogo con dati diversi); prodotti tra matrici in 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 →, inversa in 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 →.

1. Equazioni normali

La matrice dei dati con la colonna di uni e il vettore delle uscite sono $$X=\begin{bmatrix}1&0\\1&1\\1&2\\1&3\end{bmatrix},\qquad y=\begin{bmatrix}1\\3\\4\\8\end{bmatrix}.$$ Si calcolano le due matrici delle equazioni normali $X^TX\beta=X^Ty$: $$X^TX=\begin{bmatrix}4&0+1+2+3\\0+1+2+3&0+1+4+9\end{bmatrix}=\begin{bmatrix}4&6\\6&14\end{bmatrix},\qquad X^Ty=\begin{bmatrix}1+3+4+8\\0\cdot1+1\cdot3+2\cdot4+3\cdot8\end{bmatrix}=\begin{bmatrix}16\\35\end{bmatrix}.$$ Il determinante è $4\cdot14-6\cdot6=20$ e l'inversa $\frac1{20}\begin{bmatrix}14&-6\\-6&4\end{bmatrix}$. Allora $$\hat\beta=\frac1{20}\begin{bmatrix}14\cdot16-6\cdot35\\-6\cdot16+4\cdot35\end{bmatrix}=\frac1{20}\begin{bmatrix}224-210\\-96+140\end{bmatrix}=\frac1{20}\begin{bmatrix}14\\44\end{bmatrix}=\begin{bmatrix}0{,}7\\2{,}2\end{bmatrix}.$$ La retta è $\hat y=0{,}7+2{,}2x$.

Controllo con le formule scalari. $\bar x=1{,}5$, $\bar y=4$; $\sum(x_i-\bar x)(y_i-\bar y)=(-1{,}5)(-3)+(-0{,}5)(-1)+0{,}5\cdot0+1{,}5\cdot4=4{,}5+0{,}5+0+6=11$ e $\sum(x_i-\bar x)^2=2{,}25+0{,}25+0{,}25+2{,}25=5$: $\beta_1=11/5=2{,}2$, $\beta_0=4-2{,}2\cdot1{,}5=0{,}7$ ✓.

2. Residui e qualità

Previsioni: $\hat y=[0{,}7;\ 2{,}9;\ 5{,}1;\ 7{,}3]$. Residui $y-\hat y=[0{,}3;\ 0{,}1;\ -1{,}1;\ 0{,}7]$ (somma $0$). Quadrati: $0{,}09+0{,}01+1{,}21+0{,}49=1{,}8$. Quindi $\mathrm{MSE}=1{,}8/4=0{,}45$ e $\mathrm{RMSE}=\sqrt{0{,}45}=0{,}671$.

$SS_{tot}=\sum(y_i-\bar y)^2=(-3)^2+(-1)^2+0^2+4^2=9+1+0+16=26$ e $R^2=1-\frac{1{,}8}{26}=0{,}931$: la retta spiega il $93\%$ della variabilità.

Ortogonalità. $X^T(y-\hat y)$: prima componente $0{,}3+0{,}1-1{,}1+0{,}7=0$; seconda $0\cdot0{,}3+1\cdot0{,}1+2\cdot(-1{,}1)+3\cdot0{,}7=0{,}1-2{,}2+2{,}1=0$. I residui sono ortogonali alle colonne di $X$: $\hat y$ è la proiezione ortogonale di $y$ sul piano generato da $(1,1,1,1)$ e $(0,1,2,3)$ (Complemento ortogonale e proiezioni ortogonaliL'ortogonale U⊥ di un sottospazio è un sottospazio di dimensione n − dim U, e R^n = U ⊕ U⊥; ogni vettore si scompone in proiezione su U più componente ortogonale; la proiezione è il punto di U più vicino e si calcola con un sistema o con la matrice di proiezione A(AᵀA)⁻¹Aᵀ.Complemento ortogonale e proiezioni ortogonali →).

3. Previsione e modello quadratico

Per $x=5$: $\hat y=0{,}7+2{,}2\cdot5=11{,}7$.

Con il termine quadratico $X=[\mathbf1,x,x^2]$ la matrice $X$ è $4\times3$. Le equazioni normali hanno $X^TX=\begin{bmatrix}4&6&14\\6&14&36\\14&36&98\end{bmatrix}$ e $X^Ty=\begin{bmatrix}16\\35\\91\end{bmatrix}$ (l'ultima componente è $0\cdot1+1\cdot3+4\cdot4+9\cdot8=91$) e danno $\hat\beta=(1{,}2;\ 0{,}7;\ 0{,}5)$, cioè $\hat y=1{,}2+0{,}7x+0{,}5x^2$ (risolto con Python). I residui sono $[-0{,}2;\ 0{,}6;\ -0{,}6;\ 0{,}2]$, $\mathrm{MSE}=0{,}2$ e $R^2=1-\frac{0{,}8}{26}=0{,}969$: migliore della retta ($0{,}931$) sul training, per costruzione (la parabola contiene la retta come caso particolare). Con 4 punti e 3 parametri l'adattamento è però quasi forzato (con 4 parametri, un polinomio cubico, passerebbe per tutti i punti: overfitting). Per $x=5$ la parabola dà $1{,}2+0{,}7\cdot5+0{,}5\cdot25=17{,}2$, contro $11{,}7$ della retta: le due estrapolazioni sono molto diverse, e con 4 dati non c'è modo di sapere quale sia migliore senza un test set (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 →).

Verifica

import numpy as np
x = np.array([0,1,2,3.]); y = np.array([1,3,4,8.])
X = np.c_[np.ones(4), x]; b = np.linalg.solve(X.T @ X, X.T @ y)       # [0.7, 2.2]
r = y - X @ b; print(b, r, (r**2).mean(), X.T @ r)                    # residui [0.3 0.1 -1.1 0.7]; MSE 0.45; ~0
X2 = np.c_[np.ones(4), x, x**2]; b2 = np.linalg.solve(X2.T @ X2, X2.T @ y); print(b2)   # [1.2 0.7 0.5]