Chapter 2 — Supervised Learning¶

import numpy as np
from sklearn.linear_model import Ridge, Lasso, LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline


def fit_linear_models(
    X_train: np.ndarray,
    y_train: np.ndarray,
    alpha: float = 1.0,
) -> dict[str, Pipeline]:
    """Fit OLS, Ridge, and Lasso regression pipelines."""
    models: dict[str, Pipeline] = {
        "ols": Pipeline([("scaler", StandardScaler()), ("model", LinearRegression())]),
        "ridge": Pipeline([("scaler", StandardScaler()), ("model", Ridge(alpha=alpha))]),
        "lasso": Pipeline([("scaler", StandardScaler()), ("model", Lasso(alpha=alpha, max_iter=10_000))]),
    }
    for name, pipe in models.items():
        pipe.fit(X_train, y_train)
        print(f"{name:>5}: fitted")
    return models

import numpy as np


class LogisticRegressionScratch:
    """Binary logistic regression implemented from first principles."""

    def __init__(self, learning_rate: float = 0.01, n_epochs: int = 1000) -> None:
        self.lr = learning_rate
        self.n_epochs = n_epochs
        self.weights_: np.ndarray | None = None
        self.bias_: float = 0.0
        self.loss_history_: list[float] = []

    @staticmethod
    def _sigmoid(z: np.ndarray) -> np.ndarray:
        return 1.0 / (1.0 + np.exp(-np.clip(z, -500, 500)))

    @staticmethod
    def _bce_loss(y_true: np.ndarray, y_pred: np.ndarray) -> float:
        eps = 1e-15
        y_pred = np.clip(y_pred, eps, 1 - eps)
        return -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))

    def fit(self, X: np.ndarray, y: np.ndarray) -> "LogisticRegressionScratch":
        n_samples, n_features = X.shape
        self.weights_ = np.zeros(n_features)
        self.bias_ = 0.0

        for epoch in range(self.n_epochs):
            z = X @ self.weights_ + self.bias_
            y_hat = self._sigmoid(z)
            error = y_hat - y

            dw = (X.T @ error) / n_samples
            db = error.mean()

            self.weights_ -= self.lr * dw
            self.bias_ -= self.lr * db

            if epoch % 100 == 0:
                loss = self._bce_loss(y, y_hat)
                self.loss_history_.append(loss)

        return self

    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        return self._sigmoid(X @ self.weights_ + self.bias_)

    def predict(self, X: np.ndarray, threshold: float = 0.5) -> np.ndarray:
        return (self.predict_proba(X) >= threshold).astype(int)

from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

lr_pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("clf", LogisticRegression(C=1.0, solver="lbfgs", max_iter=1000, random_state=42)),
])
lr_pipeline.fit(X_train, y_train)

from sklearn.tree import DecisionTreeClassifier, export_text
from sklearn.model_selection import cross_val_score
import numpy as np

# Find optimal alpha via cross-validation
dt = DecisionTreeClassifier(random_state=42)
path = dt.cost_complexity_pruning_path(X_train, y_train)
alphas = path.ccp_alphas[::5]   # sample every 5th alpha to save time

cv_scores = []
for alpha in alphas:
    tree = DecisionTreeClassifier(ccp_alpha=alpha, random_state=42)
    scores = cross_val_score(tree, X_train, y_train, cv=5, scoring="f1_macro")
    cv_scores.append(scores.mean())

best_alpha = alphas[np.argmax(cv_scores)]
best_tree = DecisionTreeClassifier(ccp_alpha=best_alpha, random_state=42)
best_tree.fit(X_train, y_train)
print(export_text(best_tree, feature_names=feature_names, max_depth=4))

from sklearn.ensemble import RandomForestClassifier

rf = RandomForestClassifier(
    n_estimators=300,
    max_features="sqrt",
    oob_score=True,        # enable OOB estimate
    random_state=42,
    n_jobs=-1,
)
rf.fit(X_train, y_train)
print(f"OOB accuracy: {rf.oob_score_:.4f}")

from sklearn.inspection import permutation_importance
import pandas as pd

result = permutation_importance(rf, X_val, y_val, n_repeats=10, random_state=42, n_jobs=-1)
importance_df = pd.DataFrame({
    "feature": feature_names,
    "importance_mean": result.importances_mean,
    "importance_std": result.importances_std,
}).sort_values("importance_mean", ascending=False)
print(importance_df.head(10).to_string(index=False))

import xgboost as xgb
import lightgbm as lgb

# XGBoost
xgb_model = xgb.XGBClassifier(
    n_estimators=500,
    max_depth=6,
    learning_rate=0.05,
    subsample=0.8,
    colsample_bytree=0.8,
    eval_metric="logloss",
    early_stopping_rounds=50,
    random_state=42,
)
xgb_model.fit(X_train, y_train, eval_set=[(X_val, y_val)], verbose=False)

# LightGBM
lgb_model = lgb.LGBMClassifier(
    n_estimators=500,
    num_leaves=63,
    learning_rate=0.05,
    subsample=0.8,
    colsample_bytree=0.8,
    random_state=42,
)
lgb_model.fit(X_train, y_train, eval_set=[(X_val, y_val)], callbacks=[lgb.early_stopping(50, verbose=False)])

from sklearn.svm import SVC
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

svm_pipeline = Pipeline([
    ("scaler", StandardScaler()),   # SVMs require feature scaling
    ("svm", SVC(kernel="rbf", C=1.0, gamma="scale", probability=True, random_state=42)),
])
svm_pipeline.fit(X_train, y_train)

from sklearn.model_selection import GridSearchCV

param_grid = {
    "clf__C": [0.01, 0.1, 1.0, 10.0],
    "clf__gamma": ["scale", 0.001, 0.01, 0.1],
}
grid_search = GridSearchCV(
    svm_pipeline, param_grid,
    cv=5, scoring="f1_macro", n_jobs=-1, verbose=2,
)
grid_search.fit(X_train, y_train)
print(f"Best params: {grid_search.best_params_}")
print(f"Best CV F1: {grid_search.best_score_:.4f}")

import optuna
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.model_selection import cross_val_score
import numpy as np

optuna.logging.set_verbosity(optuna.logging.WARNING)


def objective(trial: optuna.Trial) -> float:
    """Optuna objective: returns mean CV F1 score to maximise."""
    params = {
        "n_estimators": trial.suggest_int("n_estimators", 100, 600, step=50),
        "max_depth": trial.suggest_int("max_depth", 2, 8),
        "learning_rate": trial.suggest_float("learning_rate", 1e-3, 0.3, log=True),
        "subsample": trial.suggest_float("subsample", 0.5, 1.0),
        "min_samples_split": trial.suggest_int("min_samples_split", 2, 20),
    }
    model = GradientBoostingClassifier(**params, random_state=42)
    scores = cross_val_score(
        model, X_train, y_train,
        cv=5, scoring="f1_macro", n_jobs=-1,
    )
    return float(scores.mean())


study = optuna.create_study(direction="maximize", sampler=optuna.samplers.TPESampler(seed=42))
study.optimize(objective, n_trials=100, timeout=600, show_progress_bar=True)

print(f"Best trial value: {study.best_value:.4f}")
print(f"Best params:      {study.best_params}")