Regressione ai minimi quadrati

Esaminiamo per primo il caso più diffuso in applicazioni reali quando il rumore sulle osservazioni è di tipo additivo gaussiano bianco.

Sia pertanto

$$y = f(\mathbf{x}, \boldsymbol\beta) + \varepsilon$$ (3.2)

funzione, in generale non lineare, di alcuni parametri $\boldsymbol\beta$ e di alcuni ingressi $\mathbf {x}$ a cui viene sommato del rumore additivo, gaussiano, a media nulla e varianza $\sigma$. Per poter stimare in maniera robusta i parametri, il numero di campioni in ingresso $\mathbf{x}=\left\{ \mathbf{x}_1 \ldots \mathbf{x}_n \right\}$ deve essere elevato, molto di più dei parametri.

Si può pensare che la funzione dei parametri non sia la stessa per tutti i campioni ma potrebbero essercene di differenti, osservando di fatto quantità differenti, funzione sempre dei medesimi parametri $\boldsymbol\beta$. In tal caso l'equazione (3.2) può venire generalizzata come

$$y_i = f_i(\boldsymbol\beta) + \varepsilon_i$$ (3.3)

avendo sottointeso con il pedice $i$ sia il tipo di funzione sia l'i-esimo campione in ingresso (di fatto un parametro costante della funzione).

Si introduce il vettore $\mathbf{r}$ definito come

$$%r_i = y_i - f_i(\mathbf{x}_i, \boldsymbol\beta) r_i = y_i - f_i(\boldsymbol\beta)$$ (3.4)

contenete il residuo associato all'osservazione i-esima (o alla funzione i-esima). $r_i$ è funzione di $\boldsymbol\beta$ tanto come $f_i$ e ne condivide le derivate (a meno di un segno con questo formalismo).

Per ottenere uno stimatore a massima verosimiglianza, la quantità da minimizzare è la negative log likelihood (sezione 2.8) della funzione (3.2). Nel caso di rumore gaussiano la funzione di verosimiglianza si scrive infatti come

$$\mathcal{L}(r_i \vert \boldsymbol\beta, \sigma) = \frac{1}{\sqrt{2 \pi \sigma_i^2}} e^{ - \frac{r_i^2}{2\sigma_i^2} }$$ (3.5)

nel caso di osservazioni indipendenti. Applicando alla funzione di verosimiglianza la definizione di negative log likelihood si ottiene che nel caso di rumore gaussiano lo stimatore alla massima verosimiglianza è il metodo dei minimi quadrati.

La regressione ai minimi quadrati è una tecnica di ottimizzazione standard per sistemi sovradimensionati che individua i parametri $\boldsymbol\beta=(\beta_1, \ldots, \beta_m)$ di una funzione $f(\mathbf{x},\boldsymbol\beta): \mathbb{R}^{m} \mapsto \mathbb{R}^{n}$ che minimizzano un errore $S$ calcolato come somma del quadrato (Sum Of Squared Error) dei residui $r_i$ su un insieme di $n$ osservazioni $y_1 \ldots y_n$:

$$S(\boldsymbol\beta) = SSE(\boldsymbol\beta) = \mathbf{r} \cd... ...= \sum_{i=1}^{n} { \Vert y_i - f_i(\boldsymbol\beta) \Vert^2 }$$ (3.6)

$S(\boldsymbol\beta)$ è definito come residual sum of squares o alternativamente come expected squared error.

$S:\mathbb{R}^{m} \mapsto \mathbb{R}$ è una funzione che viene analizzata, al variare dei parametri $\boldsymbol\beta \in \mathbb{R}^{m}$, per cercare il suo valor minimo

$$\beta^{+} = \argmin_\beta S(\beta)$$ (3.7)

Per questa ragione viene chiamata funzione obiettivo o funzione costo. Un minimo ottenuto attraverso un procedimento come quello descritto dall'equazione (3.7) viene definito minimo globale.

Un minimo globale è difficile, dal punto di vista prettamente computazionale, da individuare e normalmente si possono sfruttare tecniche per individuare solamente i minimi locali.

Sia pertanto $S(\boldsymbol\beta)$^3.1 differenziabile, ovvero $f$ differenziabile. La condizione necessaria che $\boldsymbol\beta$ sia un minimo è che, in quel punto dello spazio dei parametri, il gradiente di $S(\boldsymbol\beta)$ si annulli, ovvero

$$\frac{\partial S(\boldsymbol\beta)}{\partial \beta_j} = 2 \m... ...i(\boldsymbol\beta)}{\partial \beta_j} = 0 \qquad j=1,\ldots,m$$ (3.8)

Condizione sufficiente che un punto stazionario ( $S'(\boldsymbol\beta)=0$) sia un minimo è che $S''(\boldsymbol\beta)$ (l'hessiana) sia definita positiva. Chiaramente l'esistenza del minimo locale garantisce solo che esiste un intorno $\delta$ di $\boldsymbol\beta$ tale che la funzione $S(\boldsymbol\beta + \delta) \geq S(\boldsymbol\beta)$.

Tutta la discussione affrontata fin ora ha come ipotesi che il rumore sia additivo $\varepsilon$ con varianza costante tra tutti i campioni (homoscedasticity). Nel caso in cui il rumore di misura sia comunque gaussiano additivo a media nulla ma con varianza non costante, ogni singola osservazione $y_i$ è una variabile aleatoria indipendente a cui è associata la varianza $\sigma^{2}_i$. Intuitivamente si capisce che la regressione ottima in questo caso dovrà pesare di più i campioni con varianza bassa mentre dovranno essere pesati meno i campioni con varianza elevata. Per ottenere questo risultato si fa uso di una normalizzazione, simile a quella mostrata in sezione 2.4.1 e diretta conseguenza della likelihood di equazione (3.5), e pertanto non si deve più minimizzare la semplice somma dei residui al quadrato, ma piuttosto la somma pesata dei residui:

$$\chi^{2} = \sum^{n}_{i=1} \frac{ \Vert r_i \Vert^2 } { \sigma_i }$$ (3.9)

La funzione costo, ora somma di una variabile aleatoria di varianza unitaria al quadrato, diventa una distribuzione chi-quadro e per questo motivo viene indicata come $\chi^{2}$. Il minimo di questa funzione costo coincide con quello ottenuto in precedenza dai minimi quadrati quando la varianza è invece costante. La condizione (3.8) per ottenere il minimo si modifica anch'essa di conseguenza:

$$%\sum_{i=1}^{n} \frac{r_i}{\sigma_i} \frac{\partial f(\math... ...i(\boldsymbol\beta)}{\partial \beta_j} = 0 \qquad j=1,\ldots,m$$ (3.10)

Generalizzando ulteriormente questo concetto, quando sull'osservazione è presente del rumore gaussiano con matrice di covarianza nota $\boldsymbol\Sigma$, la Weighted Sum of Squared Error (WSSE) si può scrivere infine come

$$\chi^{2} = \mathbf{r}^{\top} \boldsymbol\Sigma^{-1} \mathbf{r}$$ (3.11)

È da notare che questa formulazione della funzione di costo equivale a quella di equazione (3.6) dove però, invece della distanza euclidea, viene usata la distanza di Mahalanobis (sezione 2.4).

Qualunque Weighted Least Squares può essere ricondotto a un problema non pesato $\boldsymbol\Sigma = I$ premoltiplicando i residui $\mathbf{r}$ (e di conseguenza le derivate) per una matrice $\mathbf{L}^{\top}$ tale che $\boldsymbol\Sigma^{-1} = \mathbf{L} \mathbf{L}^{\top}$, usando per esempio una decomposizione di Cholesky nel caso in cui tale matrice non sia diagonale.

Tutti questi stimatori, che tengono conto della varianza dell'osservazione, coincidono con il negative log likelihood per la variabile $\mathbf{y}$ perturbata da rumore gaussiano di media zero e covarianza $\boldsymbol\Sigma$.

Footnotes

...#tex2html_wrap_inline8764#^3.1

In letteratura la funzione $S$ viene spesso codificata con un fattore $1/2$ di scala per rendere il gradiente di $S$ non viziato dal fattore $2$ ed il segno concorde con $f$ per semplificare la notazione.