Chapter 37: Stacking Deep

import numpy as np
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.model_selection import KFold
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor

from kaggle_companion.activities._common import clean

SEED = 37
REPLICATES = 5          # independent worlds per library size; results are means over them
N_DEV, N_FRESH, FEATURES = 250, 1000, 10
RIDGE_ALPHA = 50.0      # the "simple regularized meta model" the chapter prefers


def generate(n, rng):
    """A bent, interacting signal that no single base model captures, plus noise with variance 1."""
    X = rng.normal(size=(n, FEATURES))
    signal = 1.2 * np.sin(1.5 * X[:, 0]) + 0.8 * X[:, 1] * X[:, 2] + 0.6 * X[:, 3] ** 2 + 0.5 * X[:, 4] - 0.4 * X[:, 5]
    return X, signal + rng.normal(0, 1.0, n)


def make_library(size, rng):
    """`size` base models: ridge, kNN and shallow trees, each seeing a random subset of 3 to 6 features."""
    return [(i % 3, np.sort(rng.choice(FEATURES, rng.integers(3, 7), replace=False))) for i in range(size)]


def base_predictions(library, X_fit, y_fit, X_new):
    """Fit every base model on (X_fit, y_fit) and return their predictions for X_new as columns."""
    out = np.empty((len(X_new), len(library)))
    for j, (kind, cols) in enumerate(library):
        model = (Ridge(alpha=1.0) if kind == 0 else KNeighborsRegressor(15) if kind == 1
                 else DecisionTreeRegressor(max_depth=4, min_samples_leaf=5, random_state=j))
        out[:, j] = model.fit(X_fit[:, cols], y_fit).predict(X_new[:, cols])
    return out


def oof_matrix(library, X, y, folds):
    matrix = np.zeros((len(y), len(library)))
    for fit, out in folds:
        matrix[out] = base_predictions(library, X[fit], y[fit], X[out])
    return matrix


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


def one_world(size, seed):
    rng = np.random.default_rng(seed)
    X, y = generate(N_DEV, rng)
    X_fresh, y_fresh = generate(N_FRESH, rng)               # new rows: the reference no estimate may touch
    library = make_library(size, rng)
    outer = list(KFold(5, shuffle=True, random_state=1).split(X))

    # Global scheme: one OOF matrix for the whole development set, then meta-CV on that matrix.
    G = oof_matrix(library, X, y, outer)
    global_cv = np.zeros(N_DEV)
    for fit, out in KFold(5, shuffle=True, random_state=2).split(G):
        global_cv[out] = LinearRegression().fit(G[fit], y[fit]).predict(G[out])

    # Nested scheme: rebuild the base stage inside each outer training partition.
    nested = {k: np.zeros(N_DEV) for k in ("ols", "ridge", "average", "best")}
    own = {k: [] for k in nested}                           # each method scored on the rows it was fitted or selected on
    fresh = {"ols": [], "ridge": []}
    for fit, out in outer:
        inner = list(KFold(3, shuffle=True, random_state=3).split(fit))
        G_in = oof_matrix(library, X[fit], y[fit], inner)    # inner OOF base predictions for outer training
        both = base_predictions(library, X[fit], y[fit], np.vstack([X[out], X_fresh]))
        G_out, G_new = both[:len(out)], both[len(out):]
        best = int(np.argmin(((G_in - y[fit][:, None]) ** 2).mean(axis=0)))
        metas = {"ols": LinearRegression().fit(G_in, y[fit]), "ridge": Ridge(alpha=RIDGE_ALPHA).fit(G_in, y[fit])}
        for k, meta in metas.items():
            nested[k][out] = meta.predict(G_out)
            own[k].append(mse(meta.predict(G_in), y[fit]))
            fresh[k].append(mse(meta.predict(G_new), y_fresh))
        nested["average"][out], own["average"] = G_out.mean(axis=1), own["average"] + [mse(G_in.mean(axis=1), y[fit])]
        nested["best"][out], own["best"] = G_out[:, best], own["best"] + [mse(G_in[:, best], y[fit])]
    return {"global_cv": mse(global_cv, y), "nested": {k: mse(v, y) for k, v in nested.items()},
            "own": {k: float(np.mean(v)) for k, v in own.items()}, "fresh": {k: float(np.mean(v)) for k, v in fresh.items()}}


def run(library_size):
    worlds = [one_world(library_size, SEED * 100 + r) for r in range(REPLICATES)]

    def avg(get):
        return float(np.mean([get(w) for w in worlds]))
    return clean({
        "library_size": library_size,
        "ols": {"own": avg(lambda w: w["own"]["ols"]), "global_cv": avg(lambda w: w["global_cv"]),
                "nested": avg(lambda w: w["nested"]["ols"]), "fresh": avg(lambda w: w["fresh"]["ols"])},
        "methods": {k: {"own": avg(lambda w, k=k: w["own"][k]), "nested": avg(lambda w, k=k: w["nested"][k])}
                    for k in ("best", "average", "ols", "ridge")},
        "ridge_fresh": avg(lambda w: w["fresh"]["ridge"]),
        "ridge_beats_ols_nested": float(np.mean([w["nested"]["ridge"] < w["nested"]["ols"] for w in worlds])),
        "noise_variance": 1.0, "replicates": REPLICATES, "development_rows": N_DEV, "fresh_rows": N_FRESH,
    })

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

VALUES = [2, 10, 30, 50]
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.ch37 import draw, NCOLS, HEIGHT

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