Simulação do dilema do Bias e da Variância

Todo bom livro de machine learning trata do dilema do bias-variância e utiliza sua decomposição na discussão das técnicas de validação cruzada e ajuste de hiper-parâmetros.

Para além da fundamentação teórica, apresento, neste post, um framework básico para decomposição do bias, variância e do erro irredutível independentementede um modelo preditivo qualquer. Ao final chegaremos aos famosos gráficos de decomposição do bias-variância apresentados nos livros, mas que comumente não possuem explicação quanto sua construção.

Revisão de estatística

Dois conceitos são chave na teoria de probabilidades e no desenvolvimento do dilema do bias-variância: a esperança matemática (valor esperado) e a variância.

A esperança estatística é definida matematicamente como a integral em todo intervalo real do produto da variável aleatória x pela sua probabilidade de ocorrência.

$$E[X] = \int_{-\infty}^{\infty} x f(x) \, dx$$

Tratando-se de uma variável x discreta, a integral assume a forma de um somatório.

$$E[X] = \sum_{i=1}^{\infty} x_i p(x_i)$$

A análise da equação acima nos permite identificar que a esperança de uma variável discreta corresponde a sua média aritmética, uma vez que \(p(x)\) pode ser expresso como a razão entre o número de observações de um determinado \(x_i\) e o número total de observações.

A variância é definida matematicamente como a esperança do quadrado da diferença entre a variável e seu valor esperado.

$$\text{Var}(X) = E\left[(X - E(X))^2\right]$$

A expansão dos termos da equação acima, nos leva a seguinte expressão:

$$Var(X)=E[(X2−2⋅X⋅E(X)+(E(X))2)]$$

E, ao aplicarmos o operador esperança no polinômio acima e utilizando a sua propriedade cumulativa chegamos a seguinte expressão:

$$\text{Var}(X) = E\left[X^2 - 2 \cdot X \cdot E(X) + (E(X))^2\right]$$

Sendo a esperança de uma variável X uma constante, a esperança dela equivalente à própria.

$$\text{Var}(X) = E(X^2) - 2 \cdot (E(X))^2 + (E(X))^2 = E(X^2) - (E(X))^2$$

Para facilitar a decomposição do bias-variância que veremos a seguir, vamos organizar a equação acima da seguinte forma:

$$E(X^2) = \text{Var}(X) + (E(X))^2$$

Decomposição do erro esperado

Dada uma variável \(y^∗\) cuja predição é desejada, iremos assumir a existência de um modelo idealizado \(f(x)\) capaz de explicar integralmente o comportamento desta variável.

$$y^* = f(x)$$

Na prática, \(y^∗\) está sujeita à sua própria variabilidade, de forma que expandimos a equação acima para a forma:

$$y = f(x) + \sigma_y$$

Apesar de desconhecermos \(f(x)\), é possível predizer os valores de y através de modelos matemáticos empíricos baseados em dados, o qual denominaremos \(\hat{f}(x)\). O objetivo de toda modelagem é obter uma função \(\hat{f}(x)\) cujos valores preditos sejam os mais próximos possíveis de y para um dado vetor de preditores \(x=[x_1,x_2,…x_p]\).

Vamos supor que para medir o quão “próximo” são os valores de \(y\) e \(\hat{f}(x)\) utilizaremos uma medida quadrática tal como o erro quadrático, isto é \((y−\hat{f}(x_0))^2\), neste caso, nosso erro esperado. Observe que estamos avaliando a função em um único ponto x0.

É fato que para cada conjunto de dados utilizado obteremos um diferente valor de \(\hat{f}(x_0)\) e, de forma equivalente, um diferente valor de \(y\). Para compensar essa variabilidade e obter um valor representativo podemos trabalhar com a esperança estatística da quantidade quadrática que estamos avaliando.

$$E(\epsilon) = E\left((y - \hat{f}(x_0))^2\right)$$

Expandindo a quantidade acima, chegamos a seguinte expressão.

$$E\left((y - \hat{f}(x_0))^2\right) = E\left(y^2 + (\hat{f}(x_0))^2 - 2 \cdot y \cdot \hat{f}(x_0)\right)$$

$$E\left((y - \hat{f}(x_0))^2\right) = E(y^2) + E((\hat{f}(x_0))^2) - 2 \cdot E(y \cdot \hat{f}(x_0))$$

Utilizando a transformação quanto a variância apresentada na seção 1, podemos expandir os termos com a esperança da variável ao quadrado para a forma:

$$E\left((y - \hat{f}(x_0))^2\right) = \text{Var}(y) + (E(y))^2 + \text{Var}(\hat{f}(x_0)) + (E(\hat{f}(x_0)))^2 - 2 \cdot E(y \cdot \hat{f}(x_0))$$

Rearranjamos os termos para fins de praticidade dos próximos passos:

$$E\left((y - \hat{f}(x_0))^2\right) = \text{Var}(y) + \text{Var}(\hat{f}(x_0)) + \left[(E(y))^2 + (E(\hat{f}(x_0)))^2 - 2 \cdot E(y \cdot \hat{f}(x_0))\right]$$

A expressão acima pode ter seu último termo simplificado pelo conceito de produto notável para a forma:

$$E\left((y - \hat{f}(x_0))^2\right) = \text{Var}(y) + \text{Var}(\hat{f}(x_0)) + \left[E(y) - E(\hat{f}(x_0))\right]^2$$

Dado que \(y=f(x)+σ_y\), a esperança de \(y - E(y)\) - é igual ao próprio \(f(x)\) visto que se espera que a esperança, o valor médio, de \(σ_y - E(σ_y)\) - seja igual a zero. Por outro lado, a variância de \(y - Var(y)\) é igual ao próprio \(σ_y\), uma vez que se espera que \(f(x)\) tenha variância nula por ser o modelo idealizado. Desta forma, o erro esperado toma a seguinte expressão:

$$E\left((y - \hat{f}(x_0))^2\right) = \sigma_y + \text{Var}(\hat{f}(x_0)) + \left[f(x) - E(\hat{f}(x_0))\right]^2$$

A expressão acima é o que a literatura indica como a decomposição em bias-variância do erro esperado.

$$E\left((y - \hat{f}(x_0))^2\right) = \text{Erro irredutível} + \text{Variância} + (\text{Bias})^2$$

Interpretação dos termos da decomposição

O bias - \([f(x)−E(\hat{f}(x_0)]\) - corresponde ao grau de proximidade do modelo proposto ao modelo ideal - \(f(x)\). Observe, contudo, que não tratamos diretamente do modelo proposto \(\hat{f}(x_0)\), mas da esperança (média) de modelos propostos - \(E(\hat{f}(x_0))\).

E o que é a esperança de um modelo?

O termo \(E(\hat{f}(x_0))\) corresponde ao valor médio de diferentes modelos \(f(x)\) aplicados no ponto \(x_0\).

Por esta razão que o bias-variância é simulado.

Na prática teremos um único banco de dados, obteremos um determinado \(\hat{f}(x_0)\) e teremos uma estimativa do erro esperado.

Não é possível decompor o erro esperado em termos do bias e da variância com um único modelo construído. Por esta razão, simulamos diferentes bancos de dados para a construção de diferentes modelos e a obtenção de diferentes estimativas pontuais.

A diferença entre o valor do modelo perfeito em um ponto \(x_0\) e o valor médio de diferentes modelos \(\hat{f}(x)\) no mesmo ponto é o que matematicamente denomina-se como bias.

E a variância?

A variância segue o mesmo princípio, sua estimativa passa pela construção de diferentes modelos vindo de diferentes bancos de dados, a aplicação destes modelos em um ponto \(x_0\) e então a estimativa da variância para cada valor de \(\hat{f}(x_0)\).

O último termo de nossa equação - erro irredutível - é a própria variabilidade de \(y\) e nos fornece uma importante conclusão: Nenhum modelo pode ter erro esperado menor que a variabilidade da variável predita.

Tá… e o onde está o tradeoff?

Certo! Até agora, verificamos como decompor o erro esperado nas suas parcelas de bias, variância e erro irredutível. Para enxergar o dilema (tradeoff) existente, passemos para um exemplo intuitivo!

Noção intuitiva do dilema do bias-variância

Digamos que possuímos o conjunto de dados apresentado na figura a seguir. Nossa variável de interesse \(y\) é uma função linear da preditora \(x\).

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
from pydataset import data

np.random.seed(42)
cars = data('cars')

filtered = cars[cars['dist'] < 80]
grouped = filtered.groupby('speed', as_index=False)['dist'].mean()
sampled = grouped.sample(n=10, random_state=42)

plt.figure(figsize=(10, 6))
sns.scatterplot(data=sampled, x='speed', y='dist')
plt.xlabel("x", fontsize=12, fontweight='bold', family='serif')
plt.ylabel("y", fontsize=12, fontweight='bold', family='serif')
plt.xticks(fontsize=12, fontweight='bold', family='serif')
plt.yticks(fontsize=12, fontweight='bold', family='serif')
plt.tight_layout()
plt.show()

Podemos a partir do conjunto de dados acima propor diferentes modelos polinomiais, digamos que ajustaremos os modelos de primeira ordem (linear), terceira ordem e quinta ordem e nona ordem. O perfil de cada um destes modelos é apresentado a seguir.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from pydataset import data

data_seed = 42
np.random.seed(data_seed)
sns.set(style="whitegrid")

df = data("cars")
df = df[df['dist'] < 80]
df_grouped = df.groupby('speed', as_index=False)['dist'].mean()
sampled = df_grouped.sample(n=10, random_state=data_seed).sort_values(by="speed")

def plot_poly(sampled_data, degree):
    X = sampled_data['speed'].values.reshape(-1, 1)
    y = sampled_data['dist'].values

    poly = PolynomialFeatures(degree)
    X_poly = poly.fit_transform(X)

    model = LinearRegression()
    model.fit(X_poly, y)

    x_range = np.linspace(0, 25, 200).reshape(-1, 1)
    x_range_poly = poly.transform(x_range)
    y_pred = model.predict(x_range_poly)

    plt.figure(figsize=(6, 4))
    plt.scatter(sampled_data['speed'], sampled_data['dist'], color='blue')
    plt.plot(x_range, y_pred, color='red')
    plt.xlim(0, 25)
    plt.ylim(0, 70)
    plt.xlabel("x", fontsize=12, fontweight='bold', family='serif')
    plt.ylabel("y", fontsize=12, fontweight='bold', family='serif')
    plt.xticks(fontsize=12, fontweight='bold', family='serif')
    plt.yticks(fontsize=12, fontweight='bold', family='serif')
    plt.title(f'Degree {degree}', fontsize=14, fontweight='bold', family='serif')
    plt.tight_layout()
    return plt

plots = [plot_poly(sampled, deg) for deg in [1, 3, 5, 7]]

fig, axs = plt.subplots(2, 2, figsize=(12, 8))
degrees = [1, 3, 5, 7]

for ax, degree in zip(axs.ravel(), degrees):
    poly = PolynomialFeatures(degree)
    X = sampled['speed'].values.reshape(-1, 1)
    y = sampled['dist'].values
    X_poly = poly.fit_transform(X)

    model = LinearRegression()
    model.fit(X_poly, y)

    x_range = np.linspace(0, 25, 200).reshape(-1, 1)
    y_pred = model.predict(poly.transform(x_range))

    ax.scatter(sampled['speed'], sampled['dist'], color='blue')
    ax.plot(x_range, y_pred, color='red')
    ax.set_xlim(0, 25)
    ax.set_ylim(0, 70)
    ax.set_title(f'Degree {degree}', fontsize=12, fontweight='bold', family='serif')
    ax.set_xlabel("x", fontsize=10, fontweight='bold', family='serif')
    ax.set_ylabel("y", fontsize=10, fontweight='bold', family='serif')
    ax.tick_params(labelsize=10)

plt.tight_layout()
plt.show()

Analisando os gráficos acima, vemos que conforme aumentamos o grau do polinômio, ou seja, a complexidade do modelo, obtemos uma curva mais aderente aos dados - dizemos, desta forma, que conforme aumentamos a complexidade do modelo, reduzimos o bias associado ao modelo.

Esta observação se torna mais clara, se ao invés de poucos pontos, povoássemos todo o espaço da variável preditora. Neste caso, apesar do último modelo ser de grau 7, seu comportamento seria equivalente ao de uma reta, portanto, se aproximando do modelo ideal linear.

import numpy as np
import pandas as pd
from pydataset import data
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

np.random.seed(42)

cars = data('cars')

dc = (
    cars[cars['dist'] < 80]
    .groupby('speed', as_index=False)
    .agg({'dist': 'mean'})
    .sample(n=10, random_state=42)
)

model = LinearRegression()
model.fit(dc[['speed']], dc['dist'])

x_seq = np.arange(0, 25.01, 0.01).reshape(-1, 1)
y_pred = model.predict(x_seq)

y_noisy = y_pred + np.random.normal(0, 0.2, size=x_seq.shape[0])

df_final = pd.DataFrame({
    'speed': x_seq.flatten(),
    'dist': y_noisy
}).query('dist > 0')

sns.set(style="whitegrid")
plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_final, x='speed', y='dist', s=10)

poly_model = make_pipeline(PolynomialFeatures(7), LinearRegression())
poly_model.fit(df_final[['speed']], df_final['dist'])
y_poly = poly_model.predict(x_seq)

plt.plot(x_seq, y_poly, color='red')
plt.xlabel('x', fontsize=12, fontweight='bold', family='serif')
plt.ylabel('y', fontsize=12, fontweight='bold', family='serif')
plt.xticks(fontsize=12, fontweight='bold', family='serif')
plt.yticks(fontsize=12, fontweight='bold', family='serif')
plt.xlim(0, 25)
plt.ylim(0, 70)
plt.title('Polynomial Fit (Degree 7)', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

Experimentemos, agora, gerar outros dados vindo do mesmo banco de dados cars e observar o comportamento dos modelos polinomiais com set.seed(43).

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from pydataset import data

data_seed = 43
np.random.seed(data_seed)
sns.set(style="whitegrid")

df = data("cars")
df = df[df['dist'] < 80]
df_grouped = df.groupby('speed', as_index=False)['dist'].mean()
sampled = df_grouped.sample(n=10, random_state=data_seed).sort_values(by="speed")

def plot_poly(sampled_data, degree):
    X = sampled_data['speed'].values.reshape(-1, 1)
    y = sampled_data['dist'].values

    poly = PolynomialFeatures(degree)
    X_poly = poly.fit_transform(X)

    model = LinearRegression()
    model.fit(X_poly, y)

    x_range = np.linspace(0, 25, 200).reshape(-1, 1)
    x_range_poly = poly.transform(x_range)
    y_pred = model.predict(x_range_poly)

    plt.figure(figsize=(6, 4))
    plt.scatter(sampled_data['speed'], sampled_data['dist'], color='blue')
    plt.plot(x_range, y_pred, color='red')
    plt.xlim(0, 25)
    plt.ylim(0, 70)
    plt.xlabel("x", fontsize=12, fontweight='bold', family='serif')
    plt.ylabel("y", fontsize=12, fontweight='bold', family='serif')
    plt.xticks(fontsize=12, fontweight='bold', family='serif')
    plt.yticks(fontsize=12, fontweight='bold', family='serif')
    plt.title(f'Degree {degree}', fontsize=14, fontweight='bold', family='serif')
    plt.tight_layout()
    return plt

plots = [plot_poly(sampled, deg) for deg in [1, 3, 5, 7]]

fig, axs = plt.subplots(2, 2, figsize=(12, 8))
degrees = [1, 3, 5, 7]

for ax, degree in zip(axs.ravel(), degrees):
    poly = PolynomialFeatures(degree)
    X = sampled['speed'].values.reshape(-1, 1)
    y = sampled['dist'].values
    X_poly = poly.fit_transform(X)

    model = LinearRegression()
    model.fit(X_poly, y)

    x_range = np.linspace(0, 25, 200).reshape(-1, 1)
    y_pred = model.predict(poly.transform(x_range))

    ax.scatter(sampled['speed'], sampled['dist'], color='blue')
    ax.plot(x_range, y_pred, color='red')
    ax.set_xlim(0, 25)
    ax.set_ylim(0, 70)
    ax.set_title(f'Degree {degree}', fontsize=12, fontweight='bold', family='serif')
    ax.set_xlabel("x", fontsize=10, fontweight='bold', family='serif')
    ax.set_ylabel("y", fontsize=10, fontweight='bold', family='serif')
    ax.tick_params(labelsize=10)

plt.tight_layout()
plt.show()

Compare o gráfico acima com o gerado anteriormente.

Observe que conforme se aumenta o grau do polinômio (complexidade do modelo) menos o modelo ajustado é similar ao equivalente na primeira figura - dizemos, neste caso, que a variância entre os modelo aumenta conforme a complexidade.

Com base nesta abordagem intuitiva fica claro que o bias reduz com a complexidade do modelo e o inverso ocorre com a variância. O erro irredutível é constante, e o erro esperado corresponde a soma destes três últimos termos.

Poderíamos repetir esse processo infinitamente, obtendo em cada uma das vezes diferentes ajustes para os parâmetros do modelo e valores diferentes para as predições em um mesmo data point. E, é exatamente isso que fazemos para realizar simulações que nos entregam o formato das curvas de bias e variância.

A partir do código a seguir, realizamos a simulação de 500 observações cujo modelo ideal é uma função quadrática, e adicionamos um ruido normal aleatório nos valores ideais da sua variável predita; ajustamos os parâmetros de modelos lineares de 0 a 8 parâmetros polinomiais e com esses valores calculamos algo que na prática não seria possível: a variância da estimação dos valores preditos e o bias dos modelos.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(42)

fun = lambda x: x**2

x_ref = np.linspace(0, 1, 11)
y_ref = fun(x_ref)
ref_df = pd.DataFrame({'x': x_ref, 'y': y_ref})

def data_sim(n=50):
    x = np.random.uniform(0, 1, n)
    y = fun(x) + np.random.normal(0, 0.2, n)
    return pd.DataFrame({'x': x, 'y': y})

results = []

for i in range(500):
    sim = data_sim()
    eps = np.random.normal(0, 0.2, len(x_ref))
    y_obs = y_ref + eps

    for degree in range(9):
        poly = PolynomialFeatures(degree)
        X_poly = poly.fit_transform(sim[['x']])
        model = LinearRegression().fit(X_poly, sim['y'])

        X_pred = poly.transform(ref_df[['x']])
        pred = model.predict(X_pred)

        for x_val, y_true, y_hat in zip(x_ref, y_obs, pred):
            results.append({
                'model': degree,
                'x': x_val,
                'y': y_true,
                'pred': y_hat
            })

df = pd.DataFrame(results)

var_df = df.groupby(['model', 'x'])['pred'].var().reset_index(name='varfx')
var_summary = var_df.groupby('model')['varfx'].mean().reset_index(name='varf')

mean_pred = df.groupby(['model', 'x'])['pred'].mean().reset_index(name='meanf')
mean_pred['bias2'] = (mean_pred['meanf'] - fun(mean_pred['x']))**2
bias_summary = mean_pred.groupby('model')['bias2'].mean().reset_index(name='bias2f')

df['error'] = (df['pred'] - df['y'])**2
mse_summary = df.groupby('model')['error'].mean().reset_index(name='meanmse')

summary = (
    bias_summary
    .merge(var_summary, on='model')
    .merge(mse_summary, on='model')
)
summary['eps'] = summary['meanmse'] - summary['bias2f'] - summary['varf']

summary_melted = summary.melt(id_vars='model', 
                               value_vars=['bias2f', 'varf', 'eps', 'meanmse'])

plt.figure(figsize=(10, 6))
for var in summary_melted['variable'].unique():
    subset = summary_melted[summary_melted['variable'] == var]
    plt.plot(subset['model'], subset['value'], label=var)
    plt.scatter(subset['model'], subset['value'])

plt.xlabel("Polynomial Degree", fontsize=12, fontweight='bold', family='serif')
plt.ylabel("Value", fontsize=12, fontweight='bold', family='serif')
plt.title("Bias-Variance Tradeoff", fontsize=14, fontweight='bold', family='serif')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

O gráfico acima resume o resultado da simulação. Como é possível verificar existe um ponto de mínimo na curva de essa médio quadrático esperado (meansmse), esse erro nunca pode ser zero (pelo ruido aleatório que adicionamos e que existirá em qualquer medição na vida real) e que se aproxima no seu limite no valor do ruído (eps). Também, conforme o grau do polinômio se aumenta, observamos uma redução do bias (bias2f), mas que por consequência eleva a variância das predições (varf).

Desta forma, conseguimos enxergar que em machine learning sempre estamos lutando por uma otimização que nunca será zero, o erro sempre existirá, cabe a nos minimizarmos ao máximo.

Até uma próxima!