1. Forecast Error

\[Y_i = f(x_i)+ \varepsilon_i, \quad i=1,\ldots,n\]

Il setting è detto Fixed-X quando sia nel training che nel test set i predictors sono trattati come non-stocastici, i.e. le realizzazioni di \(X\) sono uguali in entrambi gli insiemi di dati.
Abbiamo quindi un campione \((x_1,y_1),\ldots,(x_n,y_n)\) iid per il training set ed un campione \((x_1,y^*_1),\ldots,(x_n,y^*_n)\) iid per il test set.

In notazione matriciale, \(y\) è detto vettore di risposte e la matrice \(X\) che ha come prima colonna il vettore unitario è detta design matrix.
Il numero di colonne di \(X\) è pari al numero di basis functions (\(p\)) + 1 per includere l’intercetta nella regressione. Se abbiamo una sola variabile esplicativa e vogliamo fittare una regressione cubica, avremo 3+1=4 colonne nella design matrix.

\[\underset{n\times 1}{\mathbf{y}} = \left[ \begin{array}{c} y_1 \\ \cdots\\ y_i \\ \cdots\\ y_n \\ \end{array}\right] \qquad \underset{n\times p}{\mathbf{X}} = \left[ \begin{array}{cccccc} x_{1}^\mathsf{T} \\ x_{2}^\mathsf{T} \\ \cdots \\ x_{i}^\mathsf{T} \\ \cdots\\ x_{n}^\mathsf{T}\\ \end{array}\right] = \left[ \begin{array}{cccccc} x_{11} & x_{12} & \cdots & x_{1j} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2j} & \cdots & x_{2p} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ x_{i1} & x_{i2} & \cdots & x_{ij}& \cdots & x_{ip} \\ \cdots & \cdots & \cdots & \cdots & \cdots\\ x_{n1} & x_{n2} & \cdots & x_{nj} & \cdots & x_{np}\\ \end{array}\right]\]


Assumiamo che tra \(Y\) ed \(X\) esista una relazione (non nota) descritta dalla legge \(f()\): questo si chiama segnale e si contrappone al rumore (noise), che invece è uno scostamento rispetto alla relazione da segnale. L’obiettivo nel capire la relazione è trovare la funzione \(f()\), non interpolare tutti i punti, altrimenti incorreremmo in un errore di scarsa generalizzazione (over-ottimismo).

Più formalmente, vogliamo minimizzare nel test set - rispetto alla stima del segnale - l’errore di previsione (forecast error) definito come media degli scarti al quadrato:
\[\mathrm{ErrF} = \mathbb{E}(\mathrm{MSE}_{\mathrm{Te}})=\mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}( Y^*_i - \hat{f}(x_i))^2\right]= \frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[( Y^*_i - \hat{f}(x_i))^2\right]\] In altre parole il forecast error è il valore atteso del mean squared error calcolato nel test set, laddove consideriamo come loss function l’errore quadratico medio. In generale è il valore atteso della funzione di perdita calcolata nel test set, data la stima del segnale \(\hat{f}\).
Si noti che la stima del segnale è basata solo sul training set.

Una buona stima sarà quindi quella funzione che si avvicina molto alla vera funzione di regressione (ignota) \(f\), e non solo per caso, ma in modo sistematico, i.e. anche dato un nuovo insieme di training. Sulla base di questa considerazione possiamo individuare almeno 3 tipi di errori:

  1. Errore Irriducibile
    dato dal noise, non si può eliminare nemmeno se conoscessimo la vera funzione di regressione.

  2. Bias della Stima
    dato dalla differenza tra la vera \(f\) e la stima \(\hat{f}\).

  3. Varianza della Stima
    al variare del training set, cambierà (in modo sostanziale) la funzione stimata \(\hat{f}\)? O sarà più o meno stabile (bassa varianza) ?

2. Bias-Variance Decomposition

Focalizziamo l’attenzione sull’errore Riducibile, i.e. gli errori dati dalla stima distorta e ad alta varianza (punti 2 e 3 di cui sopra).
Si dimostra che esiste un trade-off tra i due tipi di errore: si può ottenere una stima non distorta, ma questa sarà molto sensibile al training set e cambiando dati la funzione stimata potrebbe cambiare in modo sostanziale (alta varianza).
Viceversa, si può stimare una funzione robusta a nuovi insiemi di dati, ma questa potrebbe essere sensibilmente diversa dalla vera funzione di regressione (ad es. stimo una retta ma la vera relazione è quadratica).

Per dimostrare la Bias-Variance Decomposition consideriamo la differenza attesa tra i veri valori futuri ed i valori previsti utilizzando la funzione stimata con il training set e calcolata nei valori di X (fissati):
\[\mathbb{E}[(Y^*_{i} - \hat{Y}^*_i)^2] = \mathbb{E}[(f(x_i) + \varepsilon^*_i - \hat{f}(x_i))^2] = \underbrace{\mathbb{E}[(f(x_i) - \hat{f}(x_i))^2]}_{\mathrm{Reducible}} + \underbrace{\mathbb{V}\mathrm{ar}(\varepsilon^*_i)}_{\mathrm{Irreducible}}\]

Dove la seconda uguaglianza è giustificata dal fatto che il noise è indipendente da X (e sue trasformazioni) ed ha media 0.
Guardando solo all’errore Reducible possiamo decomporre questo valore atteso applicando nuovamente il quadrato dentro l’operatore dopo aver aggiunto e sottratto la quantità \(\mathbb{E}\hat{f}(x_i)\).
Sommando su tutte le n-osservazioni si ottiene la decomposizione del forecast error come somma di errore irriducibile, bias al quadrato e varianza (di \(\hat{f}\)):
\[\mathrm{ErrF} = \sigma^2 + \underbrace{\frac{1}{n}\sum_{i=1}^{n}(\mathbb{E}\hat{f}(x_i) - f(x_i) )^2}_{\mathrm{Bias}^2} + \underbrace{\frac{1}{n}\sum_{i=1}^{n}\mathbb{V}\mathrm{ar}(\hat{f}(x_i))}_{\mathrm{Variance}}\]

3. Bias-Variance trade-off

In questa sezione finale viene fatta una simulazione per calcolare le due componenti dell’errore riducibile ed ottimizzare la scelta di \(\hat{f}\) in modo da avere una varianza ridotta, seppur al prezzo di una stima distorta. Il miglior modello, tra i diversi che verranno stimati per diverse scelte di grado polinomiale, sarà quello con errore riducibile (somma delle due componenti, calcolate in modo separato) minore.

L’approccio teorico è il seguente: si ripete B-volte la procedura di stima e si applica ad ogni riga (da 1 a n) il valore atteso e la varianza per ottenere \(\mathbb{E}\hat{f}(x_i)\) e \(\mathbb{V}\mathrm{ar}(\hat{f}(x_i))\) per ogni \(\quad i=1,\ldots,n\)

Ripeto B-volte tutte le procedure di stima per le diverse scelte di d (grado del polinomio).

# setting
sigma <- 0.01
x <- seq(0.5,3,length=30)
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)
n <- length(x)

# true regression function
plot(x, ftrue, type="l", col="blue")

# simulation polynomial regression order d >= 1
sim <- function(d) {
  y <- ftrue + rnorm(n, 0, sigma)
  fit <- lm(y ~ poly(x, degree=d)) 
  yhat <- fitted(fit)
}

# call B-times the cubic regression simulation
B <- 100 
d <- 3
set.seed(123)
yhats <- replicate(B, sim(d))
Ehatf <- apply(yhats, 1, mean)
Bias2 <- (Ehatf - ftrue)^2
Vhatf <- apply(yhats, 1, var)

# barplot Bias2 + Variance
barplot(Bias2+Vhatf, ylab="Bias2 + Var", names.arg=round(x,1), ylim=c(0,0.00012))
barplot(Vhatf,add=T, col=1, names.arg=" ")
legend("topright", c("Bias2","Var"), col=c("gray",1), pch=c(19,19))

# forecast error con d=3
errF <- mean(Bias2 + Vhatf + sigma)
errF
[1] 0.01010881
# stimo per d = 1,2,3,...,10 e scelgo dstar con forecast error minimo
ds <- 1:10
forecast_errors <- c()
for (i in 1:10) {
  ys <- replicate(B, sim(ds[i]))
  Ehats <- apply(ys, 1, mean)
  Biass <- (Ehats - ftrue)^2
  Vhats <- apply(ys, 1, var)
  forecast_errors[i] <- mean(Biass+Vhats+sigma)
}



dstar <- which.min(forecast_errors)
dstar
[1] 5
y <- ftrue + rnorm(n, mean=0, sd=sigma)
best_model <- lm(y ~ poly(x, degree=dstar))
plot(x, ftrue, col="blue", type="l")
points(x, y)
lines(x, fitted(best_model), col="red")