Decomposing mean squared error into bias and variance

Question

It is well known that an estimator's MSE can be decomposed into the sum of the variance and the squared bias. I'd like to actually perform this decomposition. Here is some code to set up and train a small model.

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.kernel_ridge import KernelRidge


x = np.linspace(0,10,1001).reshape(-1,1)
X = np.random.uniform(low = 0, high =10, size = 2000).reshape(-1,1)
y = 3*np.sin(X) + 2.5*np.random.normal(size = X.shape)

X_train, X_test, y_train, y_test = train_test_split(X,y,train_size = 0.6)

reg = KernelRidge(kernel='rbf', gamma = 1).fit(X,y)

mse = mean_squared_error(y_test, reg.predict(X_test))

print(mse)

How can I go about computing the squared bias and the variance for this estimator?

Answer 1

The MSE and its components (squared bias and variance) are random variables. Therefore, in order to accurately access these statistics you need to iterate the process many times. For this purpose, we generally use bootstrapping.

The procedure for obtaining bias and variance terms is as follows:

1. Generate multiple training data sets by bootstrapping (e.g. K=200).
2. For each set, train your model. You will end up with K=200 models.
3. For each model, predict the targets for the out-of-bag samples (samples which did not appear in the training sets). Now, for each out-of-bag sample $(x_i,y_i)$ you can obtain K predictions (say $h^1, h^2,.....,h^k$).
4. Average these predictions to obtain $\bar{h}$
5. Then, the bias of your model (for a single sample) will be $(\bar{h}-y_i)$
6. The variance will be $\sum_{k} \frac{(y_i^k-\bar{h})^2}{K-1}$
7. Do this for all out-of-bag samples and average bias and variance values for better estimates.

Answer 2

The key for understanding the previous MSE decomposition is to realise that $\theta$ represents a point variable with a fixed value and $\hat{\theta}$ represents a vector containing samples from the estimator (random variable).

As clarification, this formula does not apply to the case where $\theta$ represents a series of values and $\hat{\theta}$ contains element-wise estimations of these points.

A simple numpy example that verifies the previous formula might be useful for understanding this:

import numpy as np

# theta is the target value (point)
theta = 10
# theta_hat contains estimations of theta (vector of 1000 samples)
theta_hat = np.random.normal(theta,0.2,1000)

mse = np.mean(np.square(theta-theta_hat))

var = np.var(theta_hat)
bias = np.mean(theta_hat) - theta

print("MSE result:", mse)
print("Decomposed result:", var + bias*bias)

Which in my case returns 0.042668... in both cases.