The bias-variance decomposition can be expressed as:
\begin{align} \newcommand{\Var}{{\rm Var}} E[(y_0 - \hat{f}(x_0)) ^ 2] &= (E[\hat{f}(x_0)] - f(x_0)) ^ 2 + E[(\hat{f}(x_0) - E[\hat{f}(x_0)]) ^ 2] + \sigma ^ 2\\ &= [{\rm Bias}(\hat{f}(x_0))] ^ 2 + \Var(\hat{f}(x_0)) + \Var(\varepsilon) \end{align}
I tried to verify the expression using a simulated experiment, but the result seems to suggest that the left hand side of the expression is not equal to the right hand side. Here's the R code to reproduce the experiment:
library(tidyverse)
set.seed(1)
coefs <- rnorm(4) # generates four parameters of the underlying linear model (see below)
training_sets <- lapply(1:100, function(...) { # 100 training sets, each with 100 observations
tibble(
X1 = rnorm(100),
X2 = rnorm(100),
X3 = rnorm(100),
fX = coefs[1] * X1 + coefs[2] * X2 + coefs[3] * X3 + coefs[4], # a simple linear regression model with three predictors and an intercept
Y = fX + rnorm(100) # adds irreducible error
)
})
fX_estimates <- training_sets %>% # train a linear model on each training set
lapply(function(training_set_i) {
lm(Y ~ X1 + X2 + X3 + 1, # the last term (i.e., '1') represents the intercept
data = training_set_i)
})
test_set <- tibble( # generate a test set with one observation, from the same population as the training sets
X1 = rnorm(1),
X2 = rnorm(1),
X3 = rnorm(1),
fX = coefs[1] * X1 + coefs[2] * X2 + coefs[3] * X3 + coefs[4],
Y = fX + rnorm(1) # adds irreducible error
)
y0 <- test_set$Y
fx0 <- test_set$fX
fx0_estimates <- fX_estimates %>% # estimates of the outcome based on models trained with different training sets
sapply(function(fX_estimate_i) {
predict(fX_estimate_i, newdata = test_set)
})
lhs <- mean((y0 - fx0_estimates) ^ 2) # value of the left hand side of the expression
rhs <- (mean(fx0_estimates) - fx0) ^ 2 + mean((fx0_estimates - mean(fx0_estimates)) ^ 2) + 1 ^ 2 # value of the right hand side of the expression
Here're some of my guesses that may explain the results:
- Some misunderstanding of the bias-variance decomposition, and hence the wrong math expressions and the wrong experiment.
- Some kind of bias was introduced by the way I generated the data and the underlying
f(X)
. - The inequality was introduced by random noise. However, increasing the sample size is not helpful. Additionally, when I ran the code with different random seeds for 100 time, $E[(y_0 - \hat{f}(x_0)) ^ 2]$ has a mean of 0.66 and a standard deviation of 0.78, while $(E[\hat{f}(x_0)] - f(x_0)) ^ 2 + E[(\hat{f}(x_0) - E[\hat{f}(x_0)]) ^ 2] + \sigma ^ 2$ has a mean of 1.04 and a standard deviation of 0.03. They are still not equal, and it seems that $E[(y_0 - \hat{f}(x_0)) ^ 2]$ is much more variable than $(E[\hat{f}(x_0)] - f(x_0)) ^ 2 + E[(\hat{f}(x_0) - E[\hat{f}(x_0)]) ^ 2] + \sigma ^ 2$.
So, what seems to be the problem?