Statistics Meets Machine Learning

"Machine learning is statistics minus any checking of models and assumptions." — Richard Breiman (provocatively)

The reality is more nuanced: statistics and machine learning share deep mathematical foundations while differing in emphasis, scale, and philosophy.

The Learning Problem

Empirical Risk Minimization

Bias-Variance Tradeoff

Vapnik-Chervonenkis (VC) Dimension

Rademacher Complexity

Model Selection: AIC, BIC, and Cross-Validation

Akaike Information Criterion (AIC)

Bayesian Information Criterion (BIC)

Cross-Validation

Regularization Theory

Ensemble Methods: A Statistical Perspective

Bagging (Bootstrap Aggregating)

Random Forests

Boosting as Gradient Descent

Connections: Statistics ↔ Machine Learning

Python Implementation

import numpy as np
from sklearn.datasets import make_classification, make_regression
from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error, r2_score
import warnings
warnings.filterwarnings('ignore')

# --- Bias-Variance Decomposition ---
def bias_variance_decomposition(X, y, model_class, n_bootstrap=200):
    """Empirically decompose prediction error into bias^2, variance, and noise."""
    n = len(y)
    predictions = np.zeros((n_bootstrap, n))

    for b in range(n_bootstrap):
        idx = np.random.choice(n, n, replace=True)
        model = model_class()
        model.fit(X[idx], y[idx])
        predictions[b] = model.predict(X)

    mean_pred = predictions.mean(axis=0)
    bias_sq = np.mean((mean_pred - y) ** 2)
    variance = np.mean(predictions.var(axis=0))
    noise = np.var(y - mean_pred)

    return bias_sq, variance, noise

np.random.seed(42)
X, y = make_regression(n_samples=500, n_features=20, noise=10, random_state=42)

print("=== Bias-Variance Decomposition ===")
for name, model_cls in [("Decision Tree (depth=1)", lambda: DecisionTreeRegressor(max_depth=1)),
                          ("Decision Tree (depth=10)", lambda: DecisionTreeRegressor(max_depth=10)),
                          ("Ridge (alpha=1)", lambda: Ridge(alpha=1)),
                          ("Random Forest", lambda: RandomForestRegressor(n_estimators=100, random_state=42))]:
    b2, v, n = bias_variance_decomposition(X, y, model_cls)
    print(f"{name:30s}: Bias²={b2:8.1f}, Var={v:8.1f}, Noise={n:8.1f}")

# --- AIC / BIC for Linear Models ---
def compute_aic_bic(model, X, y):
    n = len(y)
    k = X.shape[1] + 1  # +1 for intercept
    y_pred = model.predict(X)
    rss = np.sum((y - y_pred) ** 2)
    sigma2 = rss / n
    log_likelihood = -n / 2 * (np.log(2 * np.pi * sigma2) + 1)
    aic = -2 * log_likelihood + 2 * k
    bic = -2 * log_likelihood + k * np.log(n)
    return aic, bic

from sklearn.preprocessing import PolynomialFeatures

X_poly1 = PolynomialFeatures(1).fit_transform(X[:, :5])
X_poly2 = PolynomialFeatures(2).fit_transform(X[:, :5])
X_poly3 = PolynomialFeatures(3).fit_transform(X[:, :5])

print("\n=== AIC/BIC for Model Selection ===")
for name, Xfeat in [("Linear", X_poly1), ("Quadratic", X_poly2), ("Cubic", X_poly3)]:
    model = Ridge(alpha=0.01).fit(Xfeat, y)
    aic, bic = compute_aic_bic(model, Xfeat, y)
    r2 = r2_score(y, model.predict(Xfeat))
    print(f"{name:12s}: AIC={aic:10.1f}, BIC={bic:10.1f}, R²={r2:.4f}")

# --- Cross-Validation Comparison ---
print("\n=== Cross-Validation MSE ===")
models = {
    "Ridge (a=1)": Ridge(alpha=1),
    "Lasso (a=0.1)": Lasso(alpha=0.1),
    "ElasticNet": ElasticNet(alpha=0.1, l1_ratio=0.5),
    "SVR (RBF)": SVR(kernel='rbf', C=10),
    "Random Forest": RandomForestRegressor(n_estimators=100, random_state=42),
    "Gradient Boosting": GradientBoostingRegressor(n_estimators=100, random_state=42),
}

kf = KFold(n_splits=5, shuffle=True, random_state=42)
for name, model in models.items():
    scores = cross_val_score(model, X, y, cv=kf, scoring='neg_mean_squared_error')
    mse = -scores.mean()
    se = scores.std() / np.sqrt(5)
    print(f"{name:25s}: MSE={mse:8.2f} (SE={se:.2f})")

# --- Ensemble Variance Analysis ---
print("\n=== Bagging Variance Reduction ===")
np.random.seed(42)
n_trees_list = [1, 5, 10, 25, 50, 100]
for n_trees in n_trees_list:
    oob_preds = []
    for seed in range(50):
        rf = RandomForestRegressor(n_estimators=n_trees, random_state=seed, oob_score=True)
        rf.fit(X, y)
        oob_preds.append(rf.oob_prediction_)
    oob_preds = np.array(oob_preds)
    avg_mse = np.mean([mean_squared_error(y, p) for p in oob_preds])
    avg_var = np.mean([np.var(p) for p in oob_preds])
    print(f"  {n_trees:3d} trees: Avg OOB MSE={avg_mse:.2f}, Variance={avg_var:.2f}")

Key Takeaways

Next Steps

• -> High-Dimensional Statistics
• -> Machine Learning Career Guide
• -> Meta-Analysis