Chapter 8: The Five-Stage Feature Engineering Process

import numpy as np
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import KFold

from kaggle_companion.activities._common import clean

SEED = 8
REPLICATES = 8        # independent constructed datasets; every score is their mean
N_TRAIN, N_NEW = 6000, 8000


def generate(n, sigma, rng):
    """Five features; the target is multiplicative: y = exp(f(x) + sigma * noise). exp(f) is its median."""
    X = rng.normal(size=(n, 5))
    f = 3.0 + 0.8 * X[:, 0] + 0.5 * X[:, 1] - 0.4 * X[:, 2] + 0.3 * X[:, 0] * X[:, 1]
    y = np.exp(f + sigma * rng.normal(size=n))
    return X, y, f


def model():
    return HistGradientBoostingRegressor(max_iter=100, learning_rate=0.08, max_depth=3, min_samples_leaf=30, random_state=0)


def rmse(y, p):
    return float(np.sqrt(np.mean((y - p) ** 2)))


def rmsle(y, p):
    return float(np.sqrt(np.mean((np.log1p(y) - np.log1p(np.maximum(p, 0))) ** 2)))


def one_dataset(sigma, seed):
    rng = np.random.default_rng(seed)
    X, y, _ = generate(N_TRAIN, sigma, rng)
    X_new, y_new, f_new = generate(N_NEW, sigma, rng)

    raw = model().fit(X, y).predict(X_new)                     # policy 1: fit y directly, squared error on the raw scale

    log_y = np.log1p(y)                                         # policy 2: fit log1p(y), then undo the transform
    m_new = model().fit(X, log_y).predict(X_new)                # predicted mean of log1p(y)
    naive = np.expm1(m_new)                                     # exponentiate the mean log prediction
    oof = np.zeros(N_TRAIN)                                     # policy 3: smearing. Cross-fitted log residuals e_i give
    for a, b in KFold(5, shuffle=True, random_state=0).split(X):  # mean(expm1(m + e_i)) = exp(m) * mean(exp(e_i)) - 1
        oof[b] = model().fit(X[a], log_y[a]).predict(X[b])
    smear = np.exp(m_new) * np.mean(np.exp(log_y - oof)) - 1.0

    policies = {"raw_fit": raw, "log_naive": naive, "log_smeared": smear, "constant": np.full(N_NEW, y.mean())}
    out = {k: {"rmse": rmse(y_new, p), "rmsle": rmsle(y_new, p), "mean_prediction": float(np.mean(p))}
           for k, p in policies.items()}
    out["known_mean"] = {"rmse": rmse(y_new, np.exp(f_new + sigma ** 2 / 2)), "rmsle": rmsle(y_new, np.exp(f_new + sigma ** 2 / 2))}
    out["known_median"] = {"rmse": rmse(y_new, np.exp(f_new)), "rmsle": rmsle(y_new, np.exp(f_new))}
    out["mean_target"] = float(np.mean(y_new))
    return out


def run(sigma):
    sets = [one_dataset(sigma, SEED * 100 + i) for i in range(REPLICATES)]
    result = {"sigma": sigma, "replicates": REPLICATES, "training_rows": N_TRAIN, "new_rows": N_NEW,
              "mean_target": float(np.mean([s["mean_target"] for s in sets]))}
    for name in ("raw_fit", "log_naive", "log_smeared", "constant", "known_mean", "known_median"):
        result[name] = {m: float(np.mean([s[name][m] for s in sets])) for m in sets[0][name]}
        result[name]["rmse_per_set"] = [s[name]["rmse"] for s in sets]
        result[name]["rmsle_per_set"] = [s[name]["rmsle"] for s in sets]
    return clean(result)

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

VALUES = [0.3, 0.8, 1.2, 1.6]
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.ch08 import draw, NCOLS, HEIGHT

DEFAULT = 0.8
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()