Exploring Bias-Variance Trade-off in Machine Learning with Python

Understanding Bias-Variance Trade-off in Machine Learning

The bias-variance trade-off is a fundamental concept in machine learning that affects model performance. It represents the balance between underfitting and overfitting. This slideshow will explore this concept using Python examples and practical applications.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

# Generate sample data
np.random.seed(42)
X = np.linspace(0, 10, 100).reshape(-1, 1)
y = 3 * X + np.sin(X) + np.random.normal(0, 0.5, (100, 1))

# Split data into train and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Plot the data
plt.scatter(X, y, color='blue', label='Data')
plt.xlabel('X')
plt.ylabel('y')
plt.title('Sample Data for Bias-Variance Trade-off')
plt.legend()
plt.show()

Bias: Underfitting the Data

Bias refers to the error introduced by approximating a real-world problem with a simplified model. High bias can lead to underfitting, where the model fails to capture the underlying patterns in the data.

# Fit a linear model (high bias)
model_linear = LinearRegression()
model_linear.fit(X_train, y_train)

# Predict using the linear model
y_pred_linear = model_linear.predict(X)

# Plot the results
plt.scatter(X, y, color='blue', label='Data')
plt.plot(X, y_pred_linear, color='red', label='Linear Model')
plt.xlabel('X')
plt.ylabel('y')
plt.title('High Bias Model (Underfitting)')
plt.legend()
plt.show()

# Calculate Mean Squared Error
mse_linear = np.mean((y - y_pred_linear) ** 2)
print(f"Mean Squared Error (Linear Model): {mse_linear:.4f}")

Variance: Overfitting the Data

Variance refers to the model's sensitivity to small fluctuations in the training data. High variance can lead to overfitting, where the model captures noise in the data rather than the underlying pattern.

# Fit a high-degree polynomial model (high variance)
poly_features = PolynomialFeatures(degree=15)
X_poly = poly_features.fit_transform(X)
X_train_poly = poly_features.transform(X_train)

model_poly = LinearRegression()
model_poly.fit(X_train_poly, y_train)

# Predict using the polynomial model
y_pred_poly = model_poly.predict(X_poly)

# Plot the results
plt.scatter(X, y, color='blue', label='Data')
plt.plot(X, y_pred_poly, color='green', label='Polynomial Model')
plt.xlabel('X')
plt.ylabel('y')
plt.title('High Variance Model (Overfitting)')
plt.legend()
plt.show()

# Calculate Mean Squared Error
mse_poly = np.mean((y - y_pred_poly) ** 2)
print(f"Mean Squared Error (Polynomial Model): {mse_poly:.4f}")

The Trade-off

The bias-variance trade-off involves finding the right balance between a model that is too simple (high bias) and one that is too complex (high variance). The goal is to minimize both bias and variance to achieve optimal model performance.

# Fit models with different polynomial degrees
degrees = range(1, 16)
train_errors = []
test_errors = []

for degree in degrees:
    poly_features = PolynomialFeatures(degree=degree)
    X_poly = poly_features.fit_transform(X)
    X_train_poly = poly_features.transform(X_train)
    X_test_poly = poly_features.transform(X_test)
    
    model = LinearRegression()
    model.fit(X_train_poly, y_train)
    
    train_pred = model.predict(X_train_poly)
    test_pred = model.predict(X_test_poly)
    
    train_errors.append(np.mean((y_train - train_pred) ** 2))
    test_errors.append(np.mean((y_test - test_pred) ** 2))

# Plot the results
plt.plot(degrees, train_errors, label='Training Error')
plt.plot(degrees, test_errors, label='Testing Error')
plt.xlabel('Polynomial Degree')
plt.ylabel('Mean Squared Error')
plt.title('Bias-Variance Trade-off')
plt.legend()
plt.show()

Real-life Example: Housing Price Prediction

Consider a scenario where we want to predict housing prices based on various features. We'll use a simplified dataset to demonstrate the bias-variance trade-off.

from sklearn.datasets import make_regression
from sklearn.model_selection import learning_curve

# Generate a synthetic dataset
X, y = make_regression(n_samples=1000, n_features=1, noise=20, random_state=42)

# Split the data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Function to plot learning curves
def plot_learning_curve(estimator, title, X, y, ylim=None, cv=None, n_jobs=None, train_sizes=np.linspace(.1, 1.0, 5)):
    plt.figure()
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel("Training examples")
    plt.ylabel("Score")
    train_sizes, train_scores, test_scores = learning_curve(
        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes)
    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)
    plt.grid()
    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1, color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std, alpha=0.1, color="g")
    plt.plot(train_sizes, train_scores_mean, 'o-', color="r", label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g", label="Cross-validation score")
    plt.legend(loc="best")
    return plt

# Plot learning curve for linear regression
plot_learning_curve(LinearRegression(), "Learning Curve for Linear Regression", X, y, ylim=(0.7, 1.01), cv=5)
plt.show()

Interpreting the Learning Curve

The learning curve shows how the model's performance changes as we increase the number of training examples. The gap between training and cross-validation scores represents the variance, while the difference between the maximum possible score and the training score represents the bias.

# Function to create polynomial features
def create_poly_features(X, degree):
    return PolynomialFeatures(degree=degree).fit_transform(X)

# Plot learning curves for different polynomial degrees
degrees = [1, 3, 5, 10]
for degree in degrees:
    X_poly = create_poly_features(X, degree)
    plot_learning_curve(LinearRegression(), f"Learning Curve (Degree {degree})", X_poly, y, ylim=(0.7, 1.01), cv=5)
    plt.show()

Balancing Bias and Variance

To find the optimal model complexity, we can use techniques like cross-validation to estimate the model's performance on unseen data. Let's implement k-fold cross-validation to find the best polynomial degree.

from sklearn.model_selection import cross_val_score

# Perform cross-validation for different polynomial degrees
degrees = range(1, 16)
cv_scores = []

for degree in degrees:
    X_poly = create_poly_features(X, degree)
    model = LinearRegression()
    scores = cross_val_score(model, X_poly, y, cv=5, scoring='neg_mean_squared_error')
    cv_scores.append(-scores.mean())

# Plot the results
plt.plot(degrees, cv_scores, marker='o')
plt.xlabel('Polynomial Degree')
plt.ylabel('Mean Squared Error')
plt.title('Cross-Validation Scores for Different Polynomial Degrees')
plt.show()

# Find the best degree
best_degree = degrees[np.argmin(cv_scores)]
print(f"Best polynomial degree: {best_degree}")

Real-life Example: Image Classification

In image classification tasks, the bias-variance trade-off is crucial. Let's consider a simplified example using the MNIST dataset.

from sklearn.datasets import load_digits
from sklearn.svm import SVC
from sklearn.model_selection import validation_curve

# Load the digits dataset
digits = load_digits()
X, y = digits.data, digits.target

# Define parameter range for SVM's C parameter
param_range = np.logspace(-6, 6, 13)

# Compute validation curve
train_scores, test_scores = validation_curve(
    SVC(), X, y, param_name="C", param_range=param_range, cv=5, scoring="accuracy", n_jobs=-1)

# Calculate mean and standard deviation
train_scores_mean = np.mean(train_scores, axis=1)
train_scores_std = np.std(train_scores, axis=1)
test_scores_mean = np.mean(test_scores, axis=1)
test_scores_std = np.std(test_scores, axis=1)

# Plot validation curve
plt.title("Validation Curve with SVM")
plt.xlabel("C")
plt.ylabel("Score")
plt.semilogx(param_range, train_scores_mean, label="Training score", color="r")
plt.fill_between(param_range, train_scores_mean - train_scores_std,
                 train_scores_mean + train_scores_std, alpha=0.2, color="r")
plt.semilogx(param_range, test_scores_mean, label="Cross-validation score", color="g")
plt.fill_between(param_range, test_scores_mean - test_scores_std,
                 test_scores_mean + test_scores_std, alpha=0.2, color="g")
plt.legend(loc="best")
plt.show()

Regularization: A Tool to Control Variance

Regularization is a technique used to reduce model complexity and prevent overfitting. Let's explore how L2 regularization (Ridge regression) affects the bias-variance trade-off.

from sklearn.linear_model import Ridge

# Generate sample data
np.random.seed(42)
X = np.linspace(0, 10, 100).reshape(-1, 1)
y = 3 * X + np.sin(X) + np.random.normal(0, 0.5, (100, 1))

# Create polynomial features
poly_features = PolynomialFeatures(degree=10)
X_poly = poly_features.fit_transform(X)

# Train models with different regularization strengths
alphas = [0, 0.1, 1, 10]
plt.figure(figsize=(12, 8))

for i, alpha in enumerate(alphas):
    model = Ridge(alpha=alpha)
    model.fit(X_poly, y)
    y_pred = model.predict(X_poly)
    
    plt.subplot(2, 2, i+1)
    plt.scatter(X, y, color='blue', label='Data')
    plt.plot(X, y_pred, color='red', label='Prediction')
    plt.title(f'Ridge Regression (alpha={alpha})')
    plt.xlabel('X')
    plt.ylabel('y')
    plt.legend()

plt.tight_layout()
plt.show()

Feature Selection and Dimensionality Reduction

Feature selection and dimensionality reduction techniques can help manage the bias-variance trade-off by reducing the number of features and focusing on the most important ones.

from sklearn.decomposition import PCA
from sklearn.feature_selection import SelectKBest, f_regression

# Generate sample data with irrelevant features
np.random.seed(42)
X = np.random.randn(200, 20)
y = 3 * X[:, 0] + 2 * X[:, 1] + np.random.randn(200)

# Perform PCA
pca = PCA(n_components=5)
X_pca = pca.fit_transform(X)

# Perform feature selection
selector = SelectKBest(f_regression, k=5)
X_selected = selector.fit_transform(X, y)

# Train models and evaluate
models = {
    'Full': LinearRegression().fit(X, y),
    'PCA': LinearRegression().fit(X_pca, y),
    'Selected': LinearRegression().fit(X_selected, y)
}

for name, model in models.items():
    if name == 'PCA':
        score = model.score(pca.transform(X), y)
    elif name == 'Selected':
        score = model.score(selector.transform(X), y)
    else:
        score = model.score(X, y)
    print(f"{name} model R-squared: {score:.4f}")

Ensemble Methods: Combining Models

Ensemble methods, such as Random Forests and Gradient Boosting, can help balance bias and variance by combining multiple models.

from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

# Generate sample data
np.random.seed(42)
X = np.random.rand(1000, 5)
y = 3 * X[:, 0] + 2 * X[:, 1] - X[:, 2] + 0.5 * X[:, 3] + np.random.randn(1000)

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Train models
models = {
    'Linear Regression': LinearRegression(),
    'Random Forest': RandomForestRegressor(n_estimators=100, random_state=42),
    'Gradient Boosting': GradientBoostingRegressor(n_estimators=100, random_state=42)
}

for name, model in models.items():
    model.fit(X_train, y_train)
    train_score = model.score(X_train, y_train)
    test_score = model.score(X_test, y_test)
    print(f"{name}:")
    print(f"  Train R-squared: {train_score:.4f}")
    print(f"  Test R-squared: {test_score:.4f}")

Cross-Validation: Estimating Model Performance

Cross-validation is a crucial technique for assessing model performance and selecting the best model while considering the bias-variance trade-off.

from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.tree import DecisionTreeRegressor

# Generate sample data
np.random.seed(42)
X = np.random.rand(1000, 5)
y = 3 * X[:, 0] + 2 * X[:, 1] - X[:, 2] + 0.5 * X[:, 3] + np.random.randn(1000)

# Define models
models = {
    'Linear Regression': LinearRegression(),
    'Ridge Regression': Ridge(alpha=1.0),
    'Decision Tree': DecisionTreeRegressor(max_depth=5)
}

# Perform cross-validation
cv = KFold(n_splits=5, shuffle=True, random_state=42)
for name, model in models.items():
    scores = cross_val_score(model, X, y, cv=cv, scoring='neg_mean_squared_error')
    mse_scores = -scores
    print(f"{name}:")
    print(f"  Mean MSE: {mse_scores.mean():.4f}")
    print(f"  Std MSE: {mse_scores.std():.4f}")

# Plot cross-validation results
plt.figure(figsize=(10, 6))
plt.boxplot([cross_val_score(model, X, y, cv=cv, scoring='neg_mean_squared_error') for model in models.values()], labels=models.keys())
plt.title('Cross-Validation Results')
plt.ylabel('Negative Mean Squared Error')
plt.show()

Hyperparameter Tuning: Finding the Sweet Spot

Hyperparameter tuning is essential for finding the optimal model complexity that balances bias and variance. We'll use GridSearchCV to find the best hyperparameters for a decision tree.

from sklearn.model_selection import GridSearchCV

# Generate sample data
np.random.seed(42)
X = np.random.rand(1000, 5)
y = 3 * X[:, 0] + 2 * X[:, 1] - X[:, 2] + 0.5 * X[:, 3] + np.random.randn(1000)

# Define the model and parameter grid
tree = DecisionTreeRegressor(random_state=42)
param_grid = {
    'max_depth': [3, 5, 7, 9],
    'min_samples_split': [2, 5, 10],
    'min_samples_leaf': [1, 2, 4]
}

# Perform grid search
grid_search = GridSearchCV(tree, param_grid, cv=5, scoring='neg_mean_squared_error')
grid_search.fit(X, y)

# Print results
print("Best parameters:", grid_search.best_params_)
print("Best MSE:", -grid_search.best_score_)

# Plot the results
results = grid_search.cv_results_
plt.figure(figsize=(12, 6))
for i, param in enumerate(['max_depth', 'min_samples_split', 'min_samples_leaf']):
    plt.subplot(1, 3, i+1)
    param_values = results[f'param_{param}'].data
    mse_scores = -results['mean_test_score']
    plt.plot(param_values, mse_scores)
    plt.title(f'MSE vs {param}')
    plt.xlabel(param)
    plt.ylabel('Mean Squared Error')
plt.tight_layout()
plt.show()

Bias-Variance Decomposition

Understanding the components of prediction error can help in diagnosing and addressing bias-variance trade-off issues. Let's implement a simple bias-variance decomposition.

def bias_variance_decomposition(model, X, y, test_size=0.3, n_iterations=100):
    mse_list, bias_list, var_list = [], [], []
    
    for _ in range(n_iterations):
        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size)
        model.fit(X_train, y_train)
        y_pred = model.predict(X_test)
        
        mse = np.mean((y_test - y_pred) ** 2)
        bias = np.mean((y_test - np.mean(y_pred)) ** 2)
        var = np.var(y_pred)
        
        mse_list.append(mse)
        bias_list.append(bias)
        var_list.append(var)
    
    return np.mean(mse_list), np.mean(bias_list), np.mean(var_list)

# Generate sample data
np.random.seed(42)
X = np.random.rand(1000, 1)
y = 3 * X.squeeze() + np.random.randn(1000) * 0.5

# Compare models
models = {
    'Linear Regression': LinearRegression(),
    'Decision Tree (max_depth=2)': DecisionTreeRegressor(max_depth=2),
    'Decision Tree (max_depth=10)': DecisionTreeRegressor(max_depth=10)
}

for name, model in models.items():
    mse, bias, var = bias_variance_decomposition(model, X, y)
    print(f"{name}:")
    print(f"  MSE: {mse:.4f}")
    print(f"  Bias^2: {bias:.4f}")
    print(f"  Variance: {var:.4f}")

# Plot results
results = [bias_variance_decomposition(model, X, y) for model in models.values()]
mse, bias, var = zip(*results)

plt.figure(figsize=(10, 6))
x = range(len(models))
width = 0.25
plt.bar(x, bias, width, label='Bias^2', color='b', alpha=0.7)
plt.bar([i + width for i in x], var, width, label='Variance', color='g', alpha=0.7)
plt.bar([i + 2 * width for i in x], mse, width, label='MSE', color='r', alpha=0.7)
plt.xticks([i + width for i in x], models.keys(), rotation=45, ha='right')
plt.ylabel('Error')
plt.title('Bias-Variance Decomposition')
plt.legend()
plt.tight_layout()
plt.show()

Additional Resources

For a deeper understanding of the bias-variance trade-off and related concepts in machine learning, consider exploring the following resources:

1. "Understanding the Bias-Variance Tradeoff" by Scott Fortmann-Roe https://arxiv.org/abs/1910.09457
2. "A Unified Approach to Interpreting Model Predictions" by Lundberg and Lee https://arxiv.org/abs/1705.07874
3. "Regularization for Deep Learning: A Taxonomy" by Kukačka et al. https://arxiv.org/abs/1710.10686
4. "An Overview of Multi-Task Learning in Deep Neural Networks" by Ruder https://arxiv.org/abs/1706.05098

These papers provide in-depth analyses and discussions on various aspects of model complexity, interpretability, and performance optimization in machine learning.