1. Optimism

L’ottimismo (Opt) è definito come la differenza attesa tra il Mean Squared Error del test set e quello del training set.
Si tratta di una quantità informativa dell’over-ottimismo dovuto alla stima sul training set (che avrà sempre errore inferiore rispetto al test) e permette di calcolare il forecast error (valore atteso del test error) come somma tra l’errore nel training set e l’ottimismo.

\[\begin{aligned} \mathrm{Opt} &= \mathbb{E}(\mathrm{MSE}_{\mathrm{Te}}) - \mathbb{E}(\mathrm{MSE}_{\mathrm{Tr}}) \end{aligned}\]

\[\begin{aligned} \widehat{\mathrm{ErrF}} &= \mathrm{MSE}_{\mathrm{Tr}} + \mathrm{OptF}\\ \end{aligned}\]

In generale si può dimostrare che l’Opt è pari alla media del doppio della covarianza tra il processo stocastico Y ed i suoi valori fittati:
\[\mathrm{OptF} = \frac{2}{n}\sum_{i=1}^{n}\mathbb{C}\mathrm{ov}(y_i,\hat{f}(x_i))\]

Per modelli lineari si dimostra che l’Opt è pari a \(\frac{2\sigma^2 p}{n}\), però è poco realistico assumere che la vera varianza del processo stocastico è nota. Per questo motivo si rimpiazza la quantità incognita (\(\sigma^2\)) con la sua stima campionaria:

\[\hat{\sigma}^2= \frac{\mathrm{RSS}}{n-p}= \frac{\sum_{i=1}^{n}( y_i - \hat{f}(x_i))^2}{n-p}\]

Lo stimatore del forecast error per modelli lineari, con l’Opt stimato attraverso la varianza campionaria, si chiama CP di Mallow e nel caso di Fixed-X setting è pari a:

\[\mathrm{Cp} = \mathrm{MSE}_{\mathrm{Tr}} + \frac{2 \hat{\sigma}^2 p}{n} = \frac{\mathrm{RSS}}{n}(1+ \frac{2p}{n-p})\] N.B. \(\mathrm{p}\) è il numero di parametri da stimare nel modello lineare, intercetta inclusa

# dgp : Y = f(x) + error
# caso: sigma noto 

ftrue <- c(0.4342,0.4780,0.5072,0.5258,0.5369,0.5426,0.5447,0.5444,0.5425,0.5397,0.5364,0.5329,0.5294,0.5260,0.5229,0.5200,0.5174,0.5151,0.5131,0.5113,0.5097,0.5083,0.5071,0.5061,0.5052,0.5044,0.5037,0.5032,0.5027,0.5023)

x <- seq(0.5, 3, length=30)

n <- length(x)
sigma <- 0.01

set.seed(123)

y <- ftrue + rnorm(n, 0, sigma)

fun <- function(d) {
   fit <- lm(y ~ poly(x, degree=d))
}  
  
# calcolo dell'Opt per d=1,2,...,15
ds <- 1:15
ps <- ds + 1
fits <- lapply(ds, fun)

MSEs <- unlist(lapply(fits, deviance))/n
RSSs <- n*MSEs

Opts <- 2*sigma^2*ps/n

Vars <- sigma^2*ps/n

Bias2 <- sapply(ds, function(d) 
  mean( ( ftrue - fitted(lm(ftrue ~ poly(x, degree=d))) )^2 )
  )

# calcolo del Forecast Error per d=1,2,..,15
forecast_errors <- Bias2 + Vars + sigma^2

# stima del Forecast Error per d=1,2,...,15
ErrF <- MSEs + Opts

# plot degli errori
plot(ps, MSEs, type="b", xlab="p", ylab="forecast_errors")
lines(ps, ErrF, type="b", col=2)
lines(ps, forecast_errors, type="b", col=3)
legend("topright", c("MSEs Training","MSEs_Tr + Optimism", "Forecast Errors"), col=c(1,2,3), pch=19)

# stima del Mallow's CP
# caso: sigma INCOGNITO 

hat_sigma2 <- n*MSEs/(n - ps)
hat_Opts <- 2*hat_sigma2*ps/n
CP <- MSEs + hat_Opts 

plot(ps, MSEs, type="b", xlab="p", ylab="forecast errors")
lines(ps, CP, type="b", col=4)
lines(ps, forecast_errors, type="b", col=2)
legend("topright", c("MSes Training", "Forecast Errors", "Mallow's CP"), col=c(1,2,4), pch=19)

2. Information Criteria

Si tratta di altri criteri (alternativi al forecast error e al Mallow’s CP) per decidere il miglior modello tra i vari addestrati. Ci muoviamo quindi ancora nell’ambito model selection e presentiamo i seguenti criteri informativi:

  1. AIC : Akaike Information Criteria

\[\mathrm{AIC} = -2 \cdot \mathrm{loglikelihood}(\hat{\beta},\hat{\sigma}^2) + 2p\]

  1. BIC : Bayesian Information Criteria

\[\mathrm{BIC} = -2 \cdot \mathrm{loglikelihood}(\hat{\beta},\hat{\sigma}^2) + \log(n)p\]

Poiché log(n) > 2 per ogni n>7, il BIC assume valori più grandi dell’AIC e tende a selezionare modelli più parsimoniosi in termini di numero di predictors utilizzati rispetto all’AIC.

Solo per i modelli lineari si dimostra che l’AIC è proporzionale al Mean Squared Error del training set: \[-2 \cdot \mathrm{loglikelihood}(\hat{\beta},\hat{\sigma}^2) = n \log(\mathrm{MSE}_{\mathrm{Tr}}) \] Di conseguenza, è proporzionale anche al Mallow’s CP: quindi con i modelli lineari i due criteri diventano equivalenti, perché verranno minimizzati in corrispondenza dello stesso modello.

Si conclude la lezione con un confronto grafico di AIC e BIC per i dati già presentati.

AICs <- unlist(lapply(fits, AIC))
BICs <- unlist(lapply(fits, BIC))

plot(ps, AICs, type="b", xlab="p", ylab="AIC,BIC,CP")
lines(ps, BICs, type="b", col=2)

legend("topright", c("Akaike IC","Bayesian IC"), col=c(1,2), pch=19)

Best model secondo entrambi i criteri: p=7 -> polinomio di grado d=6
N.B. viene selezionato lo stesso modello anche secondo il Mallow’s CP.

3. Random-X setting

Fin’ora abbiamo affrontato il caso Fixed-X in cui i regressori X erano variabili deterministiche.
Adesso affrontiamo il caso più realistico di variabili esplicative (sia nel training che nel test) stocastiche e vediamo come cambiano i risultati già visti.

La coppia di v.aleatorie \((X,Y)\) ha una distribuzione congiunta incognita: l’obiettivo rimane trovare una stima della relazione tra Y ed X ‘fatta bene’ (cioè con basso bias e bassa varianza, i.e. basso errore riducibile).
Condizionatamente ad \(X=x\) (realizzazioni delle variabili esplicative) la variabile target Y è pari alla funzione di regressione \(f(X)\) + l’errore irriducibile \(\varepsilon\) (noise).

\[f(x) = \mathbb{E}(Y|X=x)\] Il problema di stima si può formulare in termini di ottimizzazione: l’obiettivo è trovare una funzione \(g()\) che sia ‘vicina’ alla vera funzione \(f()\), dove ‘vicino’ significa che lo scarto quadratico è minimo. Il problema di minimizzazione in modo formale è il seguente:
\[f = \underset{g}{\arg\min} \,\, \mathbb{E}_{X,Y}[(Y - g(X))^2] \]

La soluzione è il previsore ottimo: \[f(X) = \mathbb{E}(Y|X)\]

Si dimostra che questa è la soluzione del minimization problem prendendo il condizionamento a \(X=x\) ed applicando la legge del valore atteso iterato:
\[\mathbb{E}_{X,Y}\{[Y - g(X)]^2\} = \mathbb{E}_X \mathbb{E}_{Y|X=x} \{[Y - g(x)]^2 | X=x\}\]

Il problema si ‘riduce’ quindi ad una minimizzazione punto-a-punto:
\[f(x) =\underset{c}{\arg\min} \,\, \mathbb{E}_{Y|X=x} \{[Y - c ]^2 | X=x\}\]

Aggiungo e sottraggo la soluzione e riscrivo la funzione obiettivo sfruttando le proprietà dell’operatore valore atteso \(\mathbb{E}\). Si dimsotra così facilmente che la soluzione è la media condizionata e che questa è la miglior previsione per \(Y\) dato \(X=x\).

\[f(x) = \mathbb{E}(Y|X=x)\]

Vediamo quindi in questo Random-X setting come cambiano le formule per calcolare l’Ottimismo e l’Errore di Previsione.

Se assumiamo una distribuzione congiunta Gaussiana, allora il previsore ottimo coincide con il previsore lineare e l’ottimismo diventa

\[Opt.Random = Opt.Fixed + \frac{\sigma^2p}{n}\left( \frac{p+1}{n-p-1}\right) = \frac{\sigma^2p}{n}\left(2 + \frac{p+1}{n-p-1}\right) \]

Si può quindi stimare il Forecast Error sommando all’Opt.Random l’MSE nel training set.

Se la vera varianza del processo non è nota (come in molte real-world applications) allora sostituiamo \(\sigma^2\) con la sua stima \(\hat\sigma^2=RSS/(n-p)\) ed otteniamo il Mallow’s CP in Random-X setting:

\[\mathrm{RCp} = \mathrm{Cp} + \frac{\hat{\sigma}^2p}{n}\left( \frac{p+1}{n-p-1}\right) = \frac{\mathrm{RSS}(n-1)}{(n-p)(n-p-1)}\]