Chapter 31: Deep Learning for Tabular Data

import warnings

import numpy as np
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler

from kaggle_companion.activities._common import clean

SEED = 31
REPLICATES = 6          # independent constructed datasets; every estimate is their mean
N_DEV, N_SMALL_DEV, N_NEW = 600, 60, 5000
GRID = np.linspace(0, 1, 21)   # candidate blend weights on the neural model: 0 is the GBM alone, 1 the network alone


def truth(X):
    """Smooth terms (easy for a network) plus two sharp steps and an interaction (easy for trees)."""
    smooth = np.sin(1.5 * X[:, 0]) + 0.6 * X[:, 1] ** 2 + 0.5 * X[:, 2] + 0.4 * np.tanh(2 * X[:, 3])
    steps = 1.2 * ((X[:, 4] > 0.3) & (X[:, 5] > -0.2)) + 0.8 * (X[:, 6] > 0.8)
    return smooth + steps


def generate(n, rng):
    X = rng.normal(size=(n, 8))
    return X, truth(X) + rng.normal(0, 0.5, n)


def mse(pred, y):
    return float(np.mean((pred - y) ** 2))


def best_weight(p_net, p_gbm, y):
    """The weight on the network that minimizes squared error on these rows."""
    return float(GRID[int(np.argmin([mse(w * p_net + (1 - w) * p_gbm, y) for w in GRID]))])


def credible_trees(X, y, X_dev, y_dev):
    """A credible tree baseline: tree size and the number of boosting rounds are chosen on the development rows."""
    best = None
    for leaves in (4, 8):
        trial = HistGradientBoostingRegressor(max_iter=400, learning_rate=0.08, max_leaf_nodes=leaves, min_samples_leaf=5,
                                              random_state=0).fit(X, y)
        errors = [mse(p, y_dev) for p in trial.staged_predict(X_dev)]     # development error after each boosting round
        rounds = int(np.argmin(errors)) + 1
        if best is None or errors[rounds - 1] < best[0]:
            best = (errors[rounds - 1], leaves, rounds)
    return HistGradientBoostingRegressor(max_iter=best[2], learning_rate=0.08, max_leaf_nodes=best[1], min_samples_leaf=5,
                                         random_state=0).fit(X, y)


def one_replicate(train_rows, seed):
    rng = np.random.default_rng(seed)
    X, y = generate(train_rows, rng)
    X_dev, y_dev = generate(N_DEV, rng)         # development rows: tune the trees and choose the blend weight here
    X_new, y_new = generate(N_NEW, rng)         # untouched rows: score here only
    gbm = credible_trees(X, y, X_dev, y_dev)
    scaler, mean, sd = StandardScaler().fit(X), y.mean(), y.std()      # the network needs scaled inputs and target
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")        # a small network may stop at the iteration cap; that is part of the trial
        net = MLPRegressor(hidden_layer_sizes=(32,), alpha=0.1, learning_rate_init=0.003, max_iter=300,
                           random_state=seed).fit(scaler.transform(X), (y - mean) / sd)
    predict_net = lambda Z: net.predict(scaler.transform(Z)) * sd + mean
    g_dev, n_dev, g_new, n_new = gbm.predict(X_dev), predict_net(X_dev), gbm.predict(X_new), predict_net(X_new)

    w = best_weight(n_dev, g_dev, y_dev)                                           # weight from 600 development rows
    w_small = best_weight(n_dev[:N_SMALL_DEV], g_dev[:N_SMALL_DEV], y_dev[:N_SMALL_DEV])   # weight from only 60 rows
    blend = lambda weight: weight * n_new + (1 - weight) * g_new
    return {"gbm": mse(g_new, y_new), "net": mse(n_new, y_new), "equal": mse(blend(0.5), y_new),
            "dev_weight": mse(blend(w), y_new), "small_dev_weight": mse(blend(w_small), y_new),
            "weight": w, "small_weight": w_small,
            "residual_corr": float(np.corrcoef(g_new - y_new, n_new - y_new)[0, 1])}


def run(train_rows):
    reps = [one_replicate(train_rows, SEED * 100 + i) for i in range(REPLICATES)]
    mean = lambda key: float(np.mean([r[key] for r in reps]))
    keys = ("gbm", "net", "equal", "dev_weight", "small_dev_weight")
    return clean({
        "train_rows": train_rows,
        **{k: mean(k) for k in keys},
        "weight": mean("weight"), "small_weight": mean("small_weight"),
        "small_weight_range": [min(r["small_weight"] for r in reps), max(r["small_weight"] for r in reps)],
        "residual_corr": mean("residual_corr"),
        "per_replicate_gain": {k: [r["gbm"] - r[k] for r in reps] for k in ("equal", "dev_weight", "small_dev_weight")},
        "replicates": REPLICATES, "development_rows": N_DEV, "small_development_rows": N_SMALL_DEV, "new_rows": N_NEW,
    })

import pandas as pd
from kaggle_companion.activities.ch31 import explain

VALUES = [150, 400, 1000, 3000]
results = {value: run(value) for value in VALUES}
table = pd.DataFrame({str(v): explain(results[v], v)['metrics'] for v in VALUES})
print(table.to_string())

import matplotlib.pyplot as plt
from kaggle_companion.activities.ch31 import draw, NCOLS, HEIGHT

DEFAULT = 400
fig, axes = plt.subplots(1, NCOLS, figsize=(6.6 if NCOLS == 1 else 10.4, HEIGHT), layout='constrained')
draw(axes, results[DEFAULT], DEFAULT)
plt.show()