Illustrating Bias vs. Variance

Deriving Bias and Variance

For this section, recall our model:

\[ y = f(x) + \epsilon \]

This tells us that some of \(y\) can be predicted by the true \(f(x)\), and some is just noise \(\epsilon\). In our simulation from last week, \(f(x) = x^2\), so \(y = x^2 + \epsilon\).

• We want to predict \(y\). We call our prediction \(\hat{y}\).

• Our best guess for \(f(x)\) is \(\hat{f}(x)\), where \(\hat{f}(x)\) is our model. It might be from a linear regression with 1, 2, 9, 15, etc. predictors or interactions of predictors. It might be from a k-nearest-neighbors estimation with k = 4. It might be from a regression tree with cp = .1 and minsplit=2.

• Even when we really nail \(\hat{f}(x)\) (which means \(\hat{f}(x) = f(x)\)), there is still error in our prediction because of \(\epsilon\)

  • \(y \neq \hat{y}\)

So we think of two different measures of error:

\[ EPE = E[(y - \hat{y})^2] = \underbrace{\mathbb{E}_{\mathcal{D}} \left[ \left(f(x) - \hat{f}(x) \right)^2 \right]}_\textrm{reducible error - MSE from imperfect model} + \underbrace{\mathbb{V}_{Y \mid X} \left[ Y \mid X = x \right]}_{\textrm{irreducible error from }\epsilon} \]

And:

\[ MSE(f(x), \hat{f}(x)) = \mathbb{E}_{\mathcal{D}}\left[\left(f(x) - \hat{f}(x)\right)^2\right] \]

Some of you asked about this equation from last time that decomposed our MSE:

\[ \text{MSE}\left(f(x), \hat{f}(x)\right) = \mathbb{E}_{\mathcal{D}} \left[ \left(f(x) - \hat{f}(x) \right)^2 \right] = \underbrace{\left(f(x) - \mathbb{E} \left[ \hat{f}(x) \right] \right)^2}_{\text{bias}^2 \left(\hat{f}(x) \right)} + \underbrace{\mathbb{E} \left[ \left( \hat{f}(x) - \mathbb{E} \left[ \hat{f}(x) \right] \right)^2 \right]}_{\text{var} \left(\hat{f}(x) \right)} \] This can be derived by:

\[ \begin{eqnarray*} \mathbb{E}_{\mathcal{D}} \left[ \left(f(x) - \hat{f}(x) \right)^2 \right] &=& \mathbb{E}_{\mathcal{D}} \left[\left(f(x) - E[\hat{f}(x)] + E[\hat{f}(x)] - \hat{f}(x)\right)^2 \right] \\ &=& \mathbb{E}_{\mathcal{D}} \left[\left(f(x) - E[\hat{f}(x)]\right)^2 + \left(E[\hat{f}(x)] - \hat{f}(x)\right)^2 + 2\left(f(x) - E[\hat{f}(x)]\right) \left(E[\hat{f}(x)] - \hat{f}(x)\right) \right] \\ &=& \mathbb{E}_{\mathcal{D}} \left[\left(f(x) - E[\hat{f}(x)]\right)^2\right] + \mathbb{E}_{\mathcal{D}} \left[\left(E[\hat{f}(x)] - \hat{f}(x)\right)^2\right] + \mathbb{E}_{\mathcal{D}} \left[2\left(f(x) - E[\hat{f}(x)]\right) \left(E[\hat{f}(x)] - \hat{f}(x)\right) \right] \\ &=& \left(f(x) - E[\hat{f}(x)]\right)^2 + Var\left(\hat{f}(x)\right) + 0 \end{eqnarray*} \] Let’s talk about what’s in this equation:

• MSE is the Mean Squared Error between \(f(x)\) and \(\hat{f}(x)\)
  • It does not have the \(\epsilon\) in it
• It is an expectation over all the possible \(\mathcal{D}\) draws of the data we could have
  • Because of this \(f(x)\) and \(E\left[\hat{f}(x)\right]\) can move out of the expectation. This lets us cancel that last term with the “2” in it.

The main takeaway is that, even given the error, \(\epsilon\), we still have additional error coming from our inability to perfectly get \(\hat{f}(x) = f(x)\).

Begin Content 10 code here:

We will illustrate these decompositions, most importantly the bias-variance tradeoff, through simulation. Suppose we would like to train a model to learn the true regression function function \(f(x) = x^2\).

f = function(x) {
  x ^ 2
}

To carry out a concrete simulation example, we need to fully specify the data generating process. We do so with the following R code.

get_sim_data = function(f, sample_size = NN) {
  x = runif(n = sample_size, min = 0, max = 1)
  eps = rnorm(n = sample_size, mean = 0, sd = SD.of.Bayes.Error)
  y = f(x) + eps
  tibble(x, y)
}

To completely specify the data generating process, we have made more model assumptions than simply \(\mathbb{E}[Y \mid X = x] = x^2\) and \(\mathbb{V}[Y \mid X = x] = \sigma ^ 2\). In particular,

• The \(x_i\) in \(\mathcal{D}\) are sampled from a uniform distribution over \([0, 1]\).
• The \(x_i\) and \(\epsilon\) are independent.
• The \(y_i\) in \(\mathcal{D}\) are sampled from the conditional normal distribution.

Using this setup, we will generate datasets, \(\mathcal{D}\), with a sample size \(NN\) and fit four models.

\[ \begin{aligned} \texttt{predict(fit0, x)} &= \hat{f}_0(x) = \hat{\beta}_0\\ \texttt{predict(fit1, x)} &= \hat{f}_1(x) = \hat{\beta}_0 + \hat{\beta}_1 x \\ \texttt{predict(fit2, x)} &= \hat{f}_2(x) = \hat{\beta}_0 + \hat{\beta}_1 x + \hat{\beta}_2 x^2 \\ \texttt{predict(fit9, x)} &= \hat{f}_9(x) = \hat{\beta}_0 + \hat{\beta}_1 x + \hat{\beta}_2 x^2 + \ldots + \hat{\beta}_9 x^9 \end{aligned} \]

To get a sense of the data and these four models, we generate one simulated dataset, and fit the four models.

set.seed(1)
sim_data = get_sim_data(f)

fit_0 = lm(y ~ 1,                   data = sim_data)
fit_1 = lm(y ~ poly(x, degree = 1), data = sim_data)
fit_2 = lm(y ~ poly(x, degree = 2), data = sim_data)
fit_9 = lm(y ~ poly(x, degree = 9), data = sim_data)

Plotting these four trained models, we see that the zero predictor model does very poorly. The first degree model is reasonable, but we can see that the second degree model fits much better. The ninth degree model seem rather wild.

…

We will now complete a simulation study to understand the relationship between the bias, variance, and mean squared error for the estimates for \(f(x)\) given by these four models at the point \(x = 0.90\). We use simulation to complete this task, as performing the analytical calculations would prove to be rather tedious and difficult.

set.seed(1)
n_sims = 250
n_models = 4
x = data.frame(x = 0.90) # fixed point at which we make predictions
predictions = matrix(0, nrow = n_sims, ncol = n_models)

for (sim in 1:n_sims) {

  # simulate new, random, training data
  # this is the only random portion of the bias, var, and mse calculations
  # this allows us to calculate the expectation over D
  sim_data = get_sim_data(f, sample_size = NN)

  # fit models
  fit_0 = lm(y ~ 1,                   data = sim_data)
  fit_1 = lm(y ~ poly(x, degree = 1), data = sim_data)
  fit_2 = lm(y ~ poly(x, degree = 2), data = sim_data)
  fit_9 = lm(y ~ poly(x, degree = 9), data = sim_data)

  # get predictions
  predictions[sim, 1] = predict(fit_0, x)
  predictions[sim, 2] = predict(fit_1, x)
  predictions[sim, 3] = predict(fit_2, x)
  predictions[sim, 4] = predict(fit_9, x)
}

Compile all of the results: