5.7 Regressione (metodo dei minimi quadrati)
In questa sezione introduciamo il problema generale della regressione, concentrandoci in particolare sul caso lineare con errore quadratico, tradizionalmente detto anche metodo dei minimi quadrati (ordinary least squares, OLS, in inglese).
Il problema della regressione, in generale, si può descrivere nel seguente modo: date due variabili aleatorie \(X \in E\), \(Y \in E'\) determinare una funzione \(g: E \to E'\) tale che \(Y\) sia “molto vicina” a \(g(X)\), \[ Y \sim g(X),\] a partire dall’osservazione congiunta di \((X,Y)\) (sottoforma di una o più copie, solitamente indipendenti). La variabile \(X\) è detta predittore (in inglese predictor) o variabile esplicativa (explanatory variable), mentre la variabile \(Y\) è detta risposta (response) oppure esito (outcome). Evitiamo appositamente il linguaggio trazionale di variabile “indipendente” (che sarebbe la \(X\)) e “dipendente” (la \(Y\)) per evitare di confonderlo con il concetto probabilistico di indipendenza.
In molte applicazioni, \(X\) è una variabile vettoriale, ossia \(E=\R^d\) (se \(d >1\) si parla di regressione multipla), mentre a seconda dell’obiettivo che ci si pone, \(Y\) potrebbe anche essere una variabile discreta (noi ci concentreremo al caso in cui sia una varibile continua, eventualmente anch’essa vettoriale).
Esempio 5.7 La regressione si usa anche per problemi di classificazione, in cui ad esempio bisogna “etichettare” i possibili valori di \(X\) per determinare due (o più) classi disgiunte e quindi la variabile di risposta \(g(X)\) è discreta a valori nell’insieme delle possibili etichette.
Dal punto di vista del calcolo delle probabilità, essendo l’incognita \(g\) non nota e di solito non completamente determinata dall’osservazione di \((X,Y)\), si introduce una variabile aleatoria \(G\) a valori nell’insieme delle possibili funzioni da \(E\) in \(E'\). Il problema diventa quindi determinare la legge di \(G\) sulla base dell’informazione a priori \(I\) e dei dati osservati, ossia \((X,Y)\). La regressione generalizza quindi in senso probabilistico il concetto di “curva interpolante”, o più in generale il problema di determinare una funzione il cui grafico passi per determinati punti \((x,y)\). Questa generalizzazione avviene almeno su due fronti: 1. da un lato \(G\) non è una singola funzione ma una densità di probabilità sulle funzioni (ovviamente poi si dovrà scegliere una stima, ad esempio tramite massima verosimiglianza, ossia la funzione \(g\) tale che la probabilità di osservare \(G=g\) sia massima) 2. dall’altro si introduce una ulteriore flessibilità non richiedendo che la curva interpoli esattamente i punti osservati, ma introducendo un certo “residuo” (o errore), definito spesso come la differenza tra \(Y\) e \(G(X)\) ossia \(Y-G(X)\).
Pertanto, almeno in teoria, tutto il problema si riduce, come in altre situazioni, a specificare una densità a priori per la variabile aleatoria \(G\) e usare la formula di Bayes per stimarla dopo le osservazioni \((X,Y)\), o in alternativa usare la stima di massima verosimiglianza. Tuttavia, in pratica, l’insieme delle funzioni da \(E\) in \(E'\) è troppo grande per essere trattato agevolmente (sia numericamente che analiticamente), e perciò si specifica un modello, ossia una opportuna famiglia parametrizzata di funzioni da \(E\) in \(E'\). Questo fatto può anche essere visto come un modo di introdurre una certa informazione a priori sulla struttura della funzione, non nota, ma neppure totalmente arbitraria.
Tecnicamente, per specificare un modello si suppone che la variabile aleatoria \(G\) (a valori nelle funzioni da \(E\) in \(E'\)) sia una variabile composta tramite un “parametro” \(U\) (che possiamo supporre aleatorio, essendo non noto a priori), solitamente a valori in uno spazio vettoriale \(\R^k\) con dimensione \(k\), piccola rispetto alla dimensione dello spazio di tutte le possibili funzioni. Pertanto, per ogni possibile valore \(U=u\) del parametro, è definita una funzione \(g( \cdot ; u): E \to E'\) che ad ogni \(x \in E\) associa \(g(x; u) \in E'\). Per chiarire le idee, vediamo degli esempi fondamentali.
Esempio 5.8 (modelli lineari) L’esempio piu semplice, su cui ci soffermeremo maggiormente, è il caso in cui \(E' = \R^{d'}\) e la funzione \(G\) sia lineare nel parametro \(U\in \R^k\), cioè che l’associazione \[ u \in \R^k \mapsto g( \cdot ; u) \] sia lineare, della forma \[ g(x; u) = \sum_{i=1}^k g_i(x) u_i\] per opportune funzioni (note e fissate a priori) \(g_i: E \to E' = \R^{d'}\). Notiamo che le singole funzioni \(g_i\), \(x \mapsto g_i(x)\), possono essere anche non lineari, come ad esempio \(g(x) = x^2\), anche se in molti casi pure lo sono. Il punto è che si può sempre considerare, al posto della variabile esplicativa \(X\), la variabile congiunta \(X' = (g_i(X))_{i=1}^k\), in modo che quindi si possa scrivere la dipendenza come se fosse lineare, \[ g(x; u) = \sum_{i=1}^k x_i' u_i.\] Ad esempio, il modello \[g(x; (u_1, u_2)) = u_1 x+ u_2 x^2\] è lineare (in \(u\)) ma non in \(x\), tuttavia diventa lineare se lo si pensa in termini della nuova variabile \(X' = (X, X^2)\).
Un altro punto che non viene spesso evidenziato è che anche un modello “affine” ossia \[ g(x; u) = u_0 + \sum_{i=1}^k x_i u_i \] è in realtà lineare, perché lo è nella variabile \(u = (u_0, u_1, \ldots, u_k)\), mentre alle \(g_i = x_i\) va aggiunta la \(g_0 = 1\) costante.
Esempio 5.9 (modello logistico) Un esempio di modello non lineare si ottiene componendo un modello lineare tramite una funzione non lineare, come ad esempio la funzione logistica (detta anche sigmoide per la forma del grafico) \(\ell: \R \to (0,1)\), \[\ell(z) = \frac{1}{1+e^{-z}}.\]
deltaz <- 0.01
z <- seq(-6, 6, by = deltaz)
l_z <- 1/(1 + exp(-z))
plot(z, l_z, type = "l", lwd = 3, col = miei_colori[2],
xlab = "z", ylab = "l(z)")
Figura 5.19: grafico della funzione logistica \(\ell\)
Si ottiene pertanto un modello della forma \[ g(x; u) = \ell\bra{ \sum_{i=1}^k g_i(x) u_i} = \frac{1}{1+\exp\bra{-\sum_{i=1}^k g_i(x) u_i}}. \] per delle opportune funzioni \(g_i(x)\). La regressione basata su tali modelli è detta appunto logistica e viene spesso usata per problemi di classificazione binaria (ossia per partizionare i valori di \(X\) in due classi, corrispondenti a \(\cur{Y\le 1/2}\) e \({Y>1/2}\).
Presentiamo ora il classico metodo dei minimi quadrati con la notazione introdotta sopra: si suppone \(Y \in E' = \R^{d'}\) e \(U \in \R^k\), per cui si può introdurre, come misura dell’errore nell’approssimazione di \(Y\) mediante \(g(X;U)\) la differenza \(Y-g(X;U)\), detta anche residuo. Il metodo consiste quindi, avendo osservato \(X=x\), \(Y=y\), nel determinare un valore del parametro che minimizzi il “residuo quadratico”, ossia \[ u_{\ols} \in \operatorname{arg} \min_{u \in \R^k} |y - g(x; u)|^2,\] Più in generale, di solito si dispone di \(n\) osservazioni indipendenti di coppie di variabili \((X_i,Y_i) =(x_i,y_i)\) per cui si suppone che il parametro \(U\) da stimare sia lo stesso, ossia \(Y_i \sim g(X_i, U)\), allora il metodo consiste nel minimizzare la somma dei residui quadratici: \[ u_{\ols} \in \operatorname{arg} \min_{u \in \R^k}\sum_{i=1}^n |y_i - g(x_i; u)|^2.\] Il metodo indica anche, come stima della “varianza” del residuo tipico \(Y-g(X; u_{\ols})\), la quantità detta anche errore quadratico medio (in inglese mean squared error, MSE) \[ \frac 1 n \sum_{i=1}^n |y_i - g(x_i; u_{\ols})|^2,\] anche se, come per la varianza campionaria, il denominatore \(n\) è solitamente sosituito con il numero \(n-k\) di “gradi di libertà” per ottenerne una stima “non distorta” (unbiased) – questo non cambia molto fintanto che il numero dei parametri \(k\) è molto più piccolo del numero delle osservazioni \(n\), cosa che avviene solitamente.
Esempio 5.10 (caso lineare reale) Consideriamo il caso di \(X\), \(Y\) a valori reali e un modello lineare parametrizzato da \(u=(a,b) \in \R^2\), ossia \[ g(x; u) = ax +b\] Per il metodo dei minimi quadrati, avendo \(n\) osservazioni indipenenti (supponendole tutte con il medesimo parametro \((a,b)\) da stimare) si deve quindi minimizzare la somma dei residui \[ \sum_{i=1}^n (y_i - a x_i -b)^2,\] come funzione di \((a,b)\). Imponendo che le derivate si annullino si trova un semplice sistema (lineare) nelle incognite \(a\), \(b\), che risolto permette di determinare \[ a_{\ols} = \frac{ \sum_{i=1}^n(y_i - \bar{y})(x_i - \bar{x})}{ \sum_{i=1}^n (x_i- \bar{x})^2} = \frac{ \Sigma_{xy}}{\Sigma_{xx}}\] avendo indicato con \(\Sigma\) le covarianze campionarie, e \[ b_{\ols} = \bar{y} - \frac{ \Sigma_{xy}}{\Sigma_{xx}} \bar{x}.\] In particolare, notiamo che il segno di \(a_{\ols}\) coincide con quello della covarianza campionaria \(\Sigma_{xy}\) (essendo la varianza a denominatore sempre positiva). Recuperiamo quindi il significato di positiva (o negativa) correlazione in termini della “concentrazione” della densità della variabile congiunta \((X,Y)\) intorno ad una retta con coefficiente angolare positivo (o negativo). Si può anche esprimere in alternativa \[ a_{\ols} = \rho_{xy} \frac{\sigma_y}{\sigma_x}, \quad b_{\ols} = \bar{y} - \rho_{xy} \frac{ \sigma_y}{\sigma_x} \bar{x}\] usando il coefficiente di correlazione e le deviazioni standard campionarie \[ \sigma_x = \sqrt{ \Sigma_{xx}}, \quad \sigma_y = \sqrt{ \Sigma_{yy}}, \quad \rho_{xy} = \frac{\Sigma_{xy}}{\sigma_x \sigma_y}.\]
Esempio 5.11 In R la regressione su un modello lineare è implementata tramite la funzione lm(). Vediamo un esempio basandoci sulle osservazioni del dataset Iris. Si vuole predire la lunghezza del sepalo (prima colonna) a partire da quella del petalo (terza colonna). Il parametro di intercetta \(b\) è aggiunto automaticamente, non serve specificarlo.
# la funzione y = u_0 + u_1x_1 + ... +
# u_k x_k è specificata introducendo
# una formula del tipo y ~ x1 + x2+ ...
# + xk, dove le xi sono le colonne del
# data frame. Non serve specificare
# l'intercetta perché è introdotta
# automaticamente. (Per formule
# complicate vi sono anche altri modi
# di inserirle)
x <- iris$Petal.Length
y <- iris$Sepal.Length
iris_reg_lin <- lm(y ~ x)
# L'output della funzione è una lista
# contenente diverse informazioni
# utili, tra cui i coefficienti
iris_reg_lin$coefficients## (Intercept) x
## 4.3066034 0.4089223# i residui y_i - g(x_i, u), che
# possiamo plottare con un istogramma
hist(iris_reg_lin$residuals, col = miei_colori[1],
probability = TRUE, xlab = "residui",
ylab = "densità")
# e i valori 'previsti' dal modello con
# i parametri ottenuti, che possiamo
# plottare accanto a quelli osservati
plot(x, y, type = "p", pch = 16, col = miei_colori[1],
xlab = "lunghezza del petalo", ylab = "lunghezza del sepalo")
points(x, iris_reg_lin$fitted.values, pch = 16,
col = miei_colori[2])
legend("bottomright", fill = miei_colori[1:2],
legend = c("osservati", "previsti"))
# possiamo anche aggiungere una linea
# per meglio rappresentare la retta
# interpolante del modello con il
# comando abline()
abline(iris_reg_lin$coefficients, col = miei_colori[2],
lwd = 2)
# introduciamo un data frame con le
# nuove osservazioni di lunghezze di
# petali
x_osservati <- data.frame(x = c(2, 2.2, 2.5,
3))
y_previsti <- predict(iris_reg_lin, x_osservati)
# rappresentiamo le previsioni in un
# nuovo plot
plot(x, y, type = "p", pch = 16, col = miei_colori[1],
xlab = "lunghezza del petalo", ylab = "lunghezza del sepalo")
points(x_osservati$x, y_previsti, pch = 16,
col = miei_colori[2])
legend("bottomright", fill = miei_colori[1:2],
legend = c("osservazioni", "previsioni"))
# notiamo che sono tutti sulla retta di
# regressione
abline(iris_reg_lin$coefficients, col = miei_colori[2],
lwd = 2)
Consideriamo ora un modello lineare più generale (multipla), in cui \(X \in\R^d\), \(U \in \R^{k}\), \(Y \in \R\) e \[ g(x; u) = \sum_{j=1}^k x_j u_j = x \cdot u,\] indicando con \(\cdot\) il prodotto scalare in \(\R^k\) per alleggerire la notazione. Avendo osservato \(X_i= x_i\), \(Y_i =y_i\) per \(i= 1, \ldots, n\), il metodo dei minimi quadrati consiste nel minimizzare la funzione \[ u \mapsto \sum_{i=1}^n (y_i - x_i \cdot u)^2,\] che è una funzione quadratica nelle variabili \(u = (u_j)_{j=1}^k\). Pertanto, imponendo che le derivate parziali si annullino si ottiene un sistema lineare (di \(k\) equazioni in \(k\) incognite) che ammette come soluzione esplicita il vettore \[ u_{\ols} = (x^T x)^{-1} x^Ty,\] dove \(x \in \R^{n \times d}\) è intesa come matrice le cui righe sono le osservazioni \(x_i \in \R^d\) e si suppone che la matrice \(x^Tx\in \R^{d \times d}\) sia invertibile. (tale matrice è gioca un ruolo simile alla matrice delle covarianze campionarie). La previsione è quindi data dalla funzione \[ z \mapsto g(z, u_{\ols}) = z \cdot (x^T x)^{-1} x^Ty.\]
Esempio 5.12 Il comando lm() permette di effettuare regressione lineare multipla in dimensione arbitraria. Ad esempio, possiamo considerare come predittori della lunghezza del sepalo nel dataset Iris tutte le variabili (eccetto la specie).
iris_reg_gen <- lm(Sepal.Length ~ Sepal.Width +
Petal.Length + Petal.Width, data = iris)
# Possiamo avere un 'riassunto' della
# regressione con il comando summary()
summary(iris_reg_gen)##
## Call:
## lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
## data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82816 -0.21989 0.01875 0.19709 0.84570
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.85600 0.25078 7.401 9.85e-12 ***
## Sepal.Width 0.65084 0.06665 9.765 < 2e-16 ***
## Petal.Length 0.70913 0.05672 12.502 < 2e-16 ***
## Petal.Width -0.55648 0.12755 -4.363 2.41e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3145 on 146 degrees of freedom
## Multiple R-squared: 0.8586, Adjusted R-squared: 0.8557
## F-statistic: 295.5 on 3 and 146 DF, p-value: < 2.2e-16Con la funzione summary() usata nell’esempio sopra si leggono anche altre informazioni rilevanti, come la stima della deviazione standard dei residui (detta anche errore standard dei residui, in inglese residual standard error) definita come la radice quadrata della versione “non distorta” dell’errore quadratico medio,
\[ s = \sqrt{ \frac{1}{n-k} \sum_{i=1}^n|y_i - x_i \cdot u_{\ols}|^2}.\]
Inoltre, ciascun parametro stimato, ossia ogni componente del vettore \(u_{OLS}\) è accompagnato da una stima della deviazione standard (visibile nella seconda colonna Std. Error, accanto a quella contenente la stima Estimate), definito per la componente \(j \in \cur{1, \ldots, k}\) come la quantità
\[ s_j = s \sqrt{ (x^Tx)^{-1}_{jj}}.\]
Remark (varianza spiegata). Una quantità spesso utilizzata per valutare l’efficacia della regressione è il coefficiente di determinazione (in inglese coefficient of determination) definito come
\[ R^2 = 1- \frac{ \sum_{i=1}^n |y_i -x_i \cdot u_{\ols}|^2}{ \sum_{i=1}^n |y_i - \bar{y}|^2},\]
dove riconosciamo (a meno di moltiplicare per \(1/n\) numeratore e denominatore) l’errore quadratico medio e la varianza campionaria. Esso è una quantità minore (o uguale) ad \(1\), e misura l’aderenza del modello lineare ai dati osservati, confrontandolo con il caso di un modello “costante” (per cui si otterrebbe come miglior funzione la media campionaria). Più \(R^2\) risulta vicino ad \(1\), migliore è l’aderenza ai dati osservati. Se nel modello lineare è inclusa la funzione costante, come è automatico nella funzione lm(), allora si può mostrare che \(R^2\) è anche non negativo, quindi sempre compreso tra \(0\) ed \(1\). In questo senso si intepreta come la “percentuale” di varianza (dei dati) spiegata dal modello lineare preso in considerazione.
Idealmente si vorrebbe \(R^2\) molto grande, ma bisogna prestare attenzione al fatto che una aderenza eccessiva ai dati, ossia \(R^2\) troppo vicino ad \(1\) potrebbe essere un segnale di overfit, in cui i parametri introdotti sono troppi (ossia la dimensione \(k\) è troppo grande in confronto al numero di osservazioni \(n\)) e la funzione stimata “insegue” le osservazioni senza veramente “imparare” nulla da esse, ossia fornendo previsioni poco efficaci se testato nuove osservazioni della variabile dipendente. Per superare parzialmente questo problema, si usa una versione “aggiustata” del coefficiente (in inglese adjusted \(R^2\)), anch’essa indicata nel comando summary()
Per meglio comprendere il metodo dei minimi quadrati, ne vediamo ora una intepretazione in termini di stima di massima verosimiglianza. Questo, oltre a dare una giustificazione teorica basata sul calcolo delle probabilità permette anche di comprendere meglio alcune ipotesi che ne permettono una applicazione corretta.
Con la notazione sopra, introduciamo un modello probabilistico tale che, rispetto ad una informazione a priori (prima delle osservazioni), \(X \in E\), \(Y \in E = \R^{d'}\), \(U \in \R^k\) siano variabili aleatorie per cui \[ Y = g(X;U) + W,\] dove \(W\) è una variabile (indicante il residuo) con densità gaussiana vettoriale \(\mathcal{N}(0, v Id)\) (dove \(v>0\) è un parametro), e che le tre variabili \(X\), \(U\) e \(W\) siano tra loro indipendenti. Solamente con queste ipotesi, qualsiasi sia la densità a priori di \(X\), si ottiene come funzione di verosimiglianza \[ \begin{split} L(u ; x, y) & = p(X=x, Y=y | U = u) \\ &= p(Y=y| U=u, X=x) p(X=x| U=u)\\ & = p( Y-g(x;u) = y-g(x;u) | U=u, X=x) p(X=x) \\ & = p( W = y-g(x,u) | U=u, X=x) p(X=x)\\ & = \exp\bra{ -\frac 1 {2v} |y-g(x,u)|^2 }\frac{1}{\sqrt{ 2 \pi v}} p(X=x). \end{split}\] Avendo osservato \(X=x\), \(Y=y\), la stima di massima verosimiglianza per \(U\) consiste quindi nel determinare il minimo della quantità \[ u \mapsto |y-g(x,u)|^2, \] ossia il minimo residuo quadratico. Similmente, se invece si dispone di \(n\) variabili \((X_i, W_i)\) tutte indipendenti tra loro (e indipendenti da una variabile \(U\)) tali che, per ogni \(i=1, \ldots, n\) valga \[ Y_i = g(X_i; U) + W_i,\] allora la funzione di verosimiglianza per \(U\) associata alle osservazioni \(X_i= x_i\), \(Y_i = y_i\) diventa, con calcoli analoghi e scrivendo per brevità \(x = (x_i)_{i=1}^n\), \(y= (y_i)_{i=1}^n\), \[ \begin{split} L(u ; x, y) & = p(X=x, Y=y | U = u) \\ & = \exp\bra{ -\frac 1 {2 v} \sum_{i=1}^n |y_i-g(x_i,u)|^2 }\frac{1}{\sqrt{ (2 \pi v)^n}} \prod_{i=1}^n p(X_i=x_i).\end{split}\] E quindi, qualsiasi sia la densità di \(X_i\) (anche non necessariamente la stessa per tutte le osservazioni) si trova che la stima di massima verosimiglianza per \(U\) è tale che la funzione \[ u \mapsto \sum_{i=1}^n |y_i-g(x_i,u)|^2,\] sia minima, quindi la stima \(u_{\ols}\) indicata dal metodo dei minimi quadrati.
Inoltre, poiché nelle espressioni sopra abbiamo mantenuta esplicita la dipendenza dal parametro \(v\) (la varianza del residuo \(W\)), pensando la verosimiglianza anche come funzione di tale parametro, con gli stessi calcoli per la stima del parametro varianza di una gaussiana a partire da osservazioni indipendenti (è in effetti questa la situazione), otteniamo che la stima di massima verosimiglianza per la coppia \((u, v)\) è \[ u_{\mle } = u _{\ols} \in \arg \min_{u \in \R^k} \sum_{i=1}^n |y_i-g(x_i;u)|^2\] e l’errore quadratico medio (nella versione “distorta”) \[ v_{\mle} = MSE = \frac 1 n \sum_{i=1}^n |y_i-g(x_i;u)|^2.\]
Remark. La giustificazione sopra si basa essenzialmente sull’ipotesi (a priori) che i residui siano variabili gaussiane indipendenti con gli stessi parametri (e media nulla). Questa ipotesi non si può giustificare prima di applicare il metodo, tuttavia è possibile, dopo l’applicazione, considerarne la validità, almeno qualitativamente oppure tramite test (si veda la Sezione 5.8).
Questa derivazione del metodo dei minimi quadrati come stima di massima verosimiglianza suggerisce anche un approccio bayesiano, in cui una densità a priori su \(U\) possa racchiudere dell’informazione già nota sul problema. Ad esempio, supponiamo che sia noto a priori che \(U\) non si discosta troppo da un parametro noto \(u_0\), ad esempio con una variabilità dell’ordine di \(\sigma_u>0\) (lungo ciascuna componente) e che le componenti di \(U\) siano indipendenti tra loro. Questo si traduce nell’assumere che la densità a priori per \(U\) sia vettoriale gaussiana \(\mathcal{N}(u_0, \sigma_u^2 Id)\), e quindi la formula di Bayes darebbe, dopo \(n\) osservazioni \(X_i =x_i\), \(Y_i= y_i\) (con le ipotesi sopra di indipendenza) \[ \begin{split}& p(U = u | X_i=x_i, Y_i= y_i, \, \forall i=1, \ldots, n) \\ & \propto \exp\bra{ -\frac 1 {2\sigma_u^2} |u-u_0|^2} L(u; x, y)\\ & \propto \exp\bra{ -\frac 1 {2} \bra{ \frac 1 v \sum_{i=1}^n |y_i-g(x_i;u)|^2 +\frac 1 {\sigma_u^2} |u-u_0|^2}}\end{split}\] e quindi il massimo della densità a posteriori per \(U\) si ottiene minimizzando la funzione \[ u \mapsto \frac 1 v \sum_{i=1}^n |y_i-g(x_i;u)|^2 +\frac 1 {\sigma_u^2} |u-u_0|^2. \] Rispetto al metodo dei minimi quadrati, è stato quindi introdotto un termine di regolarizzazione (o penalizzazione) alla somma dei residui, che diventa rilevante se \(u\) è troppo lontano dal parametro \(u_0\). L’introduzione di questi ed altre funzioni è spesso utile per regolarizzare appunto la soluzione fornita dal semplice metodo dei minimi quadrati (queste tecniche hanno diversi nomi a seconda del tipo di termini che si aggiungono, ad esempio ridge, weight decay, LASSO, ecc.).
L’approccio bayesiano alla regressione si può approfondire analiticamente nel caso di modelli lineari. Supponendo che \(g(x;u) = x\cdot u\), la densità a posteriori per \(U\) diventa \[ \begin{split} p(U = u & | X_i=x_i, Y_i= y_i, \, \forall i=1, \ldots, n) \\ & \propto \exp\bra{ -\frac 1 {2} \bra{ \frac 1 v \sum_{i=1}^n |y_i- x_i \cdot u|^2 +\frac 1 {\sigma_u^2} |u-u_0|^2}},\end{split}\] che è una densità gaussiana vettoriale (essendo un esponenziale di polinomio di secondo grado rispetto alla variabile \(u\)). Precisamente, con calcoli diretti che qui omettiamo, per isolare i termini lineari e quadratici rispetto alla variabile \(u\), si trova che \(U\) ha come nuovi parametri, avendo osservato \(X_i=x_i\), \(Y_i=y_i\) per \(i=1,\ldots, n\), il vettore dei valor medi \[ u_{|X=x,Y=y} = \bra{ x^T x + (v/\sigma_u^2)Id}^{-1} \bra{ x^T y + (v/\sigma_u^2) u_0}\] e la matrice delle covarianze \[ \Sigma_{U|X=x,Y=y} = v \bra{ x^T x + (v/\sigma_u^2)Id}^{-1}.\]
In particolare, la deviazione standard della componente \(j \in \cur{1, \ldots, k}\) del vettore dei parametri \(U\), si ottiene dal termine diagonale della matrice, \[\begin{split} \sigma_{U_j|X=x, Y=y} & = \sqrt{ \Var{U_j| X=x, Y=y}} \\ & = \sqrt{v \bra{x^T x +(v/\sigma_u^2) Id}^{-1}_{jj}}.\end{split}\]
Queste formule sono un po’ più complicate della stima di massima verosimiglianza per \(U\), ma utilizzando anche l’informazione a priori e permettono di meglio quantificare l’incertezza associata alla stima puntuale.
Nelle formule per la varianza e la deviazione standard, il parametro \(v\) (la varianza di \(W\)) è qui trattato come noto, mentre per ottenere un’analisi più precisa dovrebbe essere pure una variabile aleatoria (abbiamo già discusso un problema simile trattando la stima della varianza di una gaussiana dalle osservazioni). Per semplificare parzialmente il metodo bayesiano, tuttavia si può qui sostituire a \(v\) la stima di massima verosimiglianza già trovata, ossia l’errore quadratico medio nella versione “distorta”.
Osserviamo infine che, nel limite \(v<<\sigma_u^2\) (quando l’informazione a priori su \(U\) diventa insignificante perché la densità tende ad essere “uniforme” su tutto \(\R^k\)), dalle formule sopra si recuperano la stima del metodo classico dei minimi quadrati per il modello lineare \[ u_{\ols} = (x^Tx )^{-1} x^T y\] e (avendo posto \(v\) la stima di massima verosimiglianza) gli errori standard dei parametri, per \(j \in \cur{1, \ldots, k}\), \[ \sigma_j = \sqrt{v \bra{x^T x}^{-1}_{jj}}.\] Questo conclude la giustificazione del metodo dal punto di vista del calcolo delle probabilità. L’unico punto non del tutto giustificato, è che abbiamo trovato così le quantità “distorte” invece di quelle “corrette” che si utilizzano comunemente, ma ripetiamo che per la numerosità delle osservazioni, \(n\), molto più grande del numero di parametri \(k\) la differenza non è poi così grande.
Remark (altre funzioni obiettivo). Vi sono certamente altre scelte ragionevoli e a volte preferibili alla funzione quadratica come funzione di costo (o obiettivo, in inglese loss function) da minimizzare. Una scelta utile in alcuni casi è ad esempio il valore assoluto, se \(E' = \R\), per cui si determina invece \[ u_{\operatorname{LAD}} \in \arg \min_{u \in \R^k}\sum_{i=1}^n |y_i - g(x_i; u)|.\] (questo metodo è detto di least absolute deviation in inglese). Per intepretare anche queste varianti del metodo dei minimi quadrati come stime di massima verosimiglianza (oppure per introdurre termini di penalizzazione non quadratici) basta sostituire alle densità gaussiane dei residui \(W_i\) (oppure dei parametri \(U\)) opportune densità, il cui logaritmo sia proporzionale alla funzione di costo. Ad esempio per il valore assoluto, si usa quindi la densità, detta di Laplace, \(p (W = w ) \propto \exp\bra{- \frac{|w-w_0|}{b}}\) dove \(w_0\in \R\) e \(b>0\) sono opportuni parametri.
5.7.1 Esercizi
Esercizio 5.12 Si consideri il dataset ‘mtcars’ e si effettui una regressione lineare con \(X\) data dalla potenza (colonna ‘hp’) e \(Y\) il tempo impegato per percorrere un quarto di miglio da ferma (colonna ‘qsec’). Aggiungere la retta trovata allo scatterplot e verificare che essa aderisca ai punti. Quale dovrebbe essere la potenza prevista dal modello di un auto affinché percorra il quarto di miglio in 10 secondi? Si confronti la previsione quanto effettivamente accade nella realtà delle auto più veloci in produzione (https://en.wikipedia.org/wiki/List_of_fastest_production_cars_by_acceleration)