Bias-variance decomposition

Examples

Error decomposition

Take a data generator whose mean is \(f\) and whose noise has standard deviation \(\sigma\), and a test function \(\hat f\) that is not \(f\). Estimate the expected error at one input by sampling.

import numpy as np

rng = np.random.default_rng(3)

def f(x):
    return 0.3 * np.sin(10 * x) + 0.7 * x

def conditional_generator(x, n=50, sigma=0.1):
    return f(x) + sigma * rng.standard_normal(n)

def expected_error(f_hat, x, sigma=0.1):
    return np.mean((conditional_generator(x, 10**6, sigma) - f_hat(x)) ** 2)

def f_hat(x):
    return 0.1 + x

print("expected error at 0.1:  ", round(float(expected_error(f_hat, 0.1)), 5))
print("reducible error at 0.1: ", round(float((f(0.1) - f_hat(0.1)) ** 2), 5))
print("irreducible error:      ", round(float(expected_error(f, 0.1)), 5))

expected error at 0.1:   0.02501
reducible error at 0.1:  0.01499
irreducible error:       0.00998

The third number is an estimate of \(\sigma^2 = 0.1^2 = 0.01\), and the first is the sum of the other two. That holds at every input, and for every noise level.

Above, samples from the generator with \(f\) and \(\hat f(x) = 0.1 + x\). Below, the expected error of \(\hat f\) estimated from \(10^5\) samples at each input, the reducible error computed from \(f\) directly, and the irreducible error. σ_noise is the noise level of the generator.

The reducible error, and with it the expected error, changes with \(x\): it is near zero where the two curves cross and largest where they are furthest apart. The irreducible error is flat, by construction of the generator, and moving σ_noise lifts the whole expected error curve by exactly that amount.

Bias and variance

The decomposition above holds for any fixed \(\hat f\). Now let the training set vary too, and fit a polynomial of the same degree to each one.

Top: 100 training sets from the same generator, six of them drawn with their fits, together with \(f\) and the average of all 100 fits. Middle: at the input x_test, the six predictions on the left and test data on the right, with the bias, the variance and the irreducible error as measuring bars. Bottom: the three terms at every input, their sum, and the expected error measured directly from predictions and test samples.

The bottom plot is the whole story in one picture: the dashed sum lies on the measured expected error everywhere. At degree 1 the fits are nearly the same line, so the variance is invisible and the bias carries the error — at \(x = 0.45\), bias\(^2\) is 0.153 and the variance is 0.0002. By degree 6 the bias at that input has gone, 0.00003, and the variance has taken over at 0.0013.

Where flexibility is actually paid for is at the two ends of the range, where every input has data on one side only. The variance there is 0.14 at degree 6 and leaves the top of the plot by degree 10, while in the middle it is still well under the irreducible error.

Both terms can be measured directly: fit the same model to 100 training sets, and the spread of the fits at an input is the variance while the distance from their average to \(f\) is the bias.

from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

gen = np.random.default_rng(123)                      # the panel's data sets
xgrid = np.round(np.arange(0, 1.001, 0.01), 2)
x = np.sort(gen.choice(xgrid, 50)).reshape(-1, 1)
noise = gen.standard_normal((100, 50))
grid = xgrid.reshape(-1, 1)

middle, edge = 45, 100                                # x = 0.45 and x = 1

for degree in [1, 4, 6, 10]:
    fits = np.vstack([
        make_pipeline(PolynomialFeatures(degree, include_bias=False),
                      StandardScaler(), LinearRegression())
        .fit(x, f(x).ravel() + 0.1 * z)
        .predict(grid)
        for z in noise])
    bias2 = (fits.mean(axis=0) - f(grid).ravel()) ** 2
    var = fits.var(axis=0)
    print(f"degree {degree:2d}   at x=0.45: bias² {bias2[middle]:.5f} "
          f"var {var[middle]:.5f}   at x=1: var {var[edge]:.4f}")

1

Standardizing the powers before the fit costs nothing and keeps a degree 10 design matrix from losing most of its digits.

degree  1   at x=0.45: bias² 0.15302 var 0.00020   at x=1: var 0.0010
degree  4   at x=0.45: bias² 0.00210 var 0.00076   at x=1: var 0.0158
degree  6   at x=0.45: bias² 0.00003 var 0.00114   at x=1: var 0.1407
degree 10   at x=0.45: bias² 0.00000 var 0.00124   at x=1: var 5.1189