Árboles de decisión y ensambles

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

titanic = sns.load_dataset("titanic")

# Las columnas que vamos a usar, con nombres en español
columnas = {"pclass": "clase", "sex": "sexo", "age": "edad",
            "sibsp": "hermanos_pareja", "parch": "padres_hijos",
            "fare": "tarifa", "embarked": "puerto"}
X = titanic[list(columnas)].rename(columns=columnas)
X["sexo"] = X["sexo"].map({"male": "hombre", "female": "mujer"})
y = titanic["survived"]               # 1 = sobrevivió, 0 = murió

print(X.shape)
print(X.head(3))
print("Faltantes:", X.isna().sum()[X.isna().sum() > 0].to_dict())
print(f"Sobrevivió el {y.mean():.1%} de los pasajeros")

from sklearn.model_selection import train_test_split
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline, make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
from sklearn.base import clone

X_ent, X_prueba, y_ent, y_prueba = train_test_split(
    X, y, test_size=0.25, stratify=y, random_state=0)

numericas = ["clase", "edad", "hermanos_pareja", "padres_hijos", "tarifa"]
categoricas = ["sexo", "puerto"]
categoricas_a_numeros = make_pipeline(
    SimpleImputer(strategy="most_frequent"),
    # drop="if_binary": «sexo» queda en una sola columna, sexo_mujer
    OneHotEncoder(drop="if_binary", handle_unknown="ignore",
                  sparse_output=False))
preparar = ColumnTransformer(
    [("num", SimpleImputer(strategy="median"), numericas),
     ("cat", categoricas_a_numeros, categoricas)],
    verbose_feature_names_out=False)

def con_preparacion(modelo):
    """Un Pipeline: preparar los datos y después el modelo."""
    return Pipeline([("preparar", clone(preparar)), ("modelo", modelo)])

A = preparar.fit_transform(X_ent)     # solo con entrenamiento
nombres = preparar.get_feature_names_out()
y_arr = y_ent.to_numpy()
print(len(X_ent), "pasajeros para entrenar,",
      len(X_prueba), "para la prueba final")
print("Columnas:", ", ".join(nombres[:5]))
print("         ", ", ".join(nombres[5:]))

def gini(etiquetas):
    """Impureza de Gini de un grupo de etiquetas."""
    _, cuentas = np.unique(etiquetas, return_counts=True)
    p = cuentas / len(etiquetas)      # proporción de cada clase
    return 1 - np.sum(p ** 2)

print("Puro:          ", gini([1, 1, 1, 1]))
print("Mitad y mitad: ", gini([0, 1, 0, 1]))
print("Tres a uno:    ", gini([0, 0, 0, 1]))
print(f"Titanic entero: {gini(y_arr):.3f}")

def mejor_corte(x, etiquetas):
    """El umbral u de la pregunta «¿x <= u?» con mayor ganancia.

    Devuelve el mejor umbral, su ganancia, y todos los umbrales
    probados con su impureza ponderada (para graficarlos).
    """
    valores = np.unique(x)            # ordenados y sin repetir
    umbrales = (valores[:-1] + valores[1:]) / 2
    n = len(etiquetas)
    ponderadas = []
    for u in umbrales:
        izq = etiquetas[x <= u]
        der = etiquetas[x > u]
        ponderadas.append((len(izq) * gini(izq)
                           + len(der) * gini(der)) / n)
    ponderadas = np.array(ponderadas)
    k = np.argmin(ponderadas)
    ganancia = gini(etiquetas) - ponderadas[k]
    return umbrales[k], ganancia, umbrales, ponderadas

# El mejor corte de cada columna, para la raíz del árbol
for j, nombre in enumerate(nombres):
    u, ganancia, _, _ = mejor_corte(A[:, j], y_arr)
    print(f"{nombre:16s} ¿<= {u:6.2f}?   ganancia = {ganancia:.3f}")

hombres = A[:, list(nombres).index("sexo_mujer")] == 0
edad_h = A[hombres, list(nombres).index("edad")]
u, ganancia, umbrales, ponderadas = mejor_corte(edad_h, y_arr[hombres])
print(f"Entre los hombres: ¿edad <= {u:g}?   ganancia = {ganancia:.3f}")
menores = y_arr[hombres][edad_h <= u]
print(f"Hombres de {u:g} años o menos: {len(menores)},"
      f" sobrevivió el {menores.mean():.0%}")

plt.plot(umbrales, ponderadas, color="C0")
plt.scatter([u], [ponderadas.min()], color="C1", s=120, zorder=3)
# Un texto con flecha que señala el punto
plt.annotate(f"Mejor corte: {u:g} años", xy=(u, ponderadas.min()),
             xytext=(30, 0.284),
             arrowprops=dict(arrowstyle="->", color="C0"))
plt.xlabel("Umbral de edad (años)")
plt.ylabel("Gini ponderado")
plt.title("Hombres: Gini según el corte")
plt.show()

from sklearn.tree import DecisionTreeClassifier, plot_tree

arbol2 = con_preparacion(DecisionTreeClassifier(max_depth=2,
                                                random_state=0))
arbol2.fit(X_ent, y_ent)

plt.figure(figsize=(4.8, 3.6))
cajas = plot_tree(arbol2["modelo"], feature_names=nombres,
                  class_names=["murió", "vivió"], label="none",
                  impurity=False, proportion=True, precision=2,
                  rounded=True, fontsize=12)
for caja in cajas:
    texto = caja.get_text().strip()
    if texto in ("True", "False"):
        caja.set_text("sí" if texto == "True" else "no")
    elif texto.endswith("vivió"):
        caja.get_bbox_patch().set(facecolor="C3", edgecolor="C3")
    else:
        caja.get_bbox_patch().set(facecolor="white", edgecolor="C2")
plt.show()

print(f"Exactitud en entrenamiento: {arbol2.score(X_ent, y_ent):.3f}")

from sklearn.datasets import make_moons

X2, y2 = make_moons(n_samples=300, noise=0.3, random_state=0)
X2_ent, X2_val, y2_ent, y2_val = train_test_split(
    X2, y2, test_size=0.5, random_state=0)

# Una rejilla fina para pintar la región de cada clase
xx, yy = np.meshgrid(np.linspace(-2, 3, 300), np.linspace(-1.5, 2, 300))
rejilla = np.column_stack([xx.ravel(), yy.ravel()])

def dibujar_puntos(eje):
    """Clase 1: puntos llenos. Clase 0: puntos huecos."""
    # .T deja las dos columnas como filas, y * las pasa por separado
    # a scatter, como x e y
    eje.scatter(*X2_ent[y2_ent == 1].T, s=12, color="C0")
    eje.scatter(*X2_ent[y2_ent == 0].T, s=12, facecolors="white",
                edgecolors="C0", linewidths=0.8)
    eje.set_xticks([])
    eje.set_yticks([])

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
dibujar_puntos(ejes[0, 0])
ejes[0, 0].set_title("Los datos")
for eje, prof in zip([ejes[0, 1], ejes[1, 0], ejes[1, 1]], [1, 3, 10]):
    a = DecisionTreeClassifier(max_depth=prof, random_state=0)
    a.fit(X2_ent, y2_ent)
    z = a.predict(rejilla).reshape(xx.shape)
    # Se pinta solo la región donde el árbol predice la clase 1
    eje.contourf(xx, yy, z, levels=[0.5, 1.5], colors=["C0"], alpha=0.15)
    dibujar_puntos(eje)
    eje.set_title(f"Profundidad {prof}")
    print(f"Profundidad {prof:2d}: {a.get_n_leaves():2d} hojas,"
          f" entrenamiento {a.score(X2_ent, y2_ent):.2f},"
          f" validación {a.score(X2_val, y2_val):.2f}")
plt.show()

from sklearn.model_selection import cross_val_score

print("max_depth   entrenamiento   validación cruzada")
for prof in [1, 2, 3, 4, 6, 8, None]:
    modelo = con_preparacion(DecisionTreeClassifier(max_depth=prof,
                                                    random_state=0))
    vc = cross_val_score(modelo, X_ent, y_ent, cv=5).mean()
    ent = modelo.fit(X_ent, y_ent).score(X_ent, y_ent)
    print(f"{str(prof):>9s}   {ent:13.3f}   {vc:18.3f}")

print("\nmin_samples_leaf   entrenamiento   validación cruzada")
for hoja in [1, 5, 10, 20]:
    modelo = con_preparacion(DecisionTreeClassifier(
        min_samples_leaf=hoja, random_state=0))
    vc = cross_val_score(modelo, X_ent, y_ent, cv=5).mean()
    ent = modelo.fit(X_ent, y_ent).score(X_ent, y_ent)
    print(f"{hoja:16d}   {ent:13.3f}   {vc:18.3f}")

rng = np.random.default_rng(0)
# Un pliegue de validación dentro del entrenamiento
A_e, A_v, y_e, y_v = train_test_split(A, y_arr, test_size=0.25,
                                      stratify=y_arr, random_state=1)
predicciones = []
for i in range(5):
    muestra = rng.integers(0, len(A_e), size=len(A_e))  # bootstrap
    a = DecisionTreeClassifier(random_state=0)
    a.fit(A_e[muestra], y_e[muestra])
    predicciones.append(a.predict(A_v))
    acierto = np.mean(predicciones[-1] == y_v)
    print(f"Árbol {i + 1}: ve {len(np.unique(muestra))} de"
          f" {len(A_e)} pasajeros, validación {acierto:.3f}")

predicciones = np.array(predicciones)
# Un pasajero tiene desacuerdo si no todos los árboles dicen lo mismo
distintos = predicciones.min(axis=0) != predicciones.max(axis=0)
desacuerdo = np.mean(distintos)
print(f"Los 5 árboles no coinciden en el {desacuerdo:.0%}"
      " de los pasajeros de validación")
votos = (predicciones.mean(axis=0) > 0.5).astype(int)
print(f"Votación de los 5 árboles: validación {np.mean(votos == y_v):.3f}")

from sklearn.ensemble import BaggingClassifier, RandomForestClassifier

bagging = con_preparacion(BaggingClassifier(
    DecisionTreeClassifier(), n_estimators=200, random_state=0))
bosque = con_preparacion(RandomForestClassifier(
    n_estimators=200, oob_score=True, random_state=0))

for nombre, modelo in [("Bagging", bagging), ("Random forest", bosque)]:
    vc = cross_val_score(modelo, X_ent, y_ent, cv=5).mean()
    print(f"{nombre:14s} validación cruzada = {vc:.3f}")

bosque.fit(X_ent, y_ent)
print(f"Random forest, exactitud OOB: {bosque['modelo'].oob_score_:.3f}")

bosque2 = RandomForestClassifier(n_estimators=200, min_samples_leaf=5,
                                 random_state=0)
bosque2.fit(X2_ent, y2_ent)
prob = bosque2.predict_proba(rejilla)[:, 1].reshape(xx.shape)

fig, eje = plt.subplots()
for desde in [0.5, 0.625, 0.75, 0.875]:
    eje.contourf(xx, yy, prob, levels=[desde, 1.01], colors=["C0"],
                 alpha=0.08)
eje.contour(xx, yy, prob, levels=[0.5], colors=["C1"], linewidths=2.5)
dibujar_puntos(eje)
eje.set_title("Bosque de 200 árboles: frontera y votos")
plt.show()

print(f"Bosque: entrenamiento {bosque2.score(X2_ent, y2_ent):.2f},"
      f" validación {bosque2.score(X2_val, y2_val):.2f}")

from sklearn.tree import DecisionTreeRegressor

rng = np.random.default_rng(0)
x = np.sort(rng.uniform(0, 6, 80))
onda = np.sin(x) + rng.normal(0, 0.25, 80)   # datos sintéticos
Xx = x.reshape(-1, 1)

tasa = 0.1
pred = np.full(len(x), onda.mean())   # paso 1: el promedio
historia = {}                         # predicción tras algunos árboles
errores = []
for m in range(1, 101):
    residuos = onda - pred            # paso 2: lo que falta
    a = DecisionTreeRegressor(max_depth=2, random_state=0)
    a.fit(Xx, residuos)               # paso 3: aprender los residuos
    pred = pred + tasa * a.predict(Xx)  # paso 4: corregir un poco
    errores.append(np.mean((onda - pred) ** 2))
    if m in (1, 5, 50):
        historia[m] = pred.copy()

for m in (1, 5, 50, 100):
    print(f"Después de {m:3d} árboles: MSE = {errores[m - 1]:.3f}")

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
for eje, m in zip([ejes[0, 0], ejes[0, 1], ejes[1, 0]], (1, 5, 50)):
    eje.vlines(x, historia[m], onda, color="C2", linewidth=0.8)
    eje.scatter(x, onda, s=8, color="C0")
    eje.plot(x, historia[m], color="C1" if m == 50 else "C0")
    eje.set_title(f"{m} {'árbol' if m == 1 else 'árboles'}")
    eje.set_xticks([])
    eje.set_yticks([])
ejes[1, 1].plot(range(1, 101), errores, color="C0")
ejes[1, 1].set_title("MSE")
ejes[1, 1].set_xlabel("Árboles")
plt.show()

from sklearn.ensemble import GradientBoostingRegressor

gb = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1,
                               max_depth=2, random_state=0)
gb.fit(Xx, onda)
diferencia = np.abs(gb.predict(Xx) - pred).max()
print(f"MSE de scikit-learn: {np.mean((onda - gb.predict(Xx)) ** 2):.3f}")
print(f"Mayor diferencia con nuestro bucle: {diferencia:.1e}")

from sklearn.model_selection import StratifiedKFold
from sklearn.ensemble import GradientBoostingClassifier

cantidades = [1, 2, 5, 10, 20, 50, 100, 200, 400, 1000]
pliegues = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
acierto_bosque = np.zeros(len(cantidades))
acierto_boost = np.zeros(len(cantidades))
for ent, val in pliegues.split(A, y_arr):
    for i, n in enumerate(cantidades):
        rf = RandomForestClassifier(n_estimators=n, random_state=0)
        rf.fit(A[ent], y_arr[ent])
        acierto_bosque[i] += rf.score(A[val], y_arr[val]) / 5
    gbc = GradientBoostingClassifier(n_estimators=1000, learning_rate=0.1,
                                     max_depth=3, random_state=0)
    gbc.fit(A[ent], y_arr[ent])
    etapas = list(gbc.staged_predict(A[val]))  # tras cada árbol
    for i, n in enumerate(cantidades):
        acierto_boost[i] += np.mean(etapas[n - 1] == y_arr[val]) / 5

for n, b, g in zip(cantidades, acierto_bosque, acierto_boost):
    print(f"{n:4d} árboles:  bosque {b:.3f}   boosting {g:.3f}")

plt.plot(cantidades, acierto_bosque, marker="o", color="C1",
         label="Random forest")
plt.plot(cantidades, acierto_boost, marker="o", color="C2",
         label="Gradient boosting")
plt.xscale("log")
plt.xlabel("Número de árboles (escala logarítmica)")
plt.ylabel("Exactitud (validación)")
plt.title("Más árboles: el bosque se estabiliza")
plt.legend()
plt.show()

from sklearn.ensemble import HistGradientBoostingClassifier

modelos = {
    "Árbol (profundidad 4)": DecisionTreeClassifier(max_depth=4,
                                                    random_state=0),
    "Bagging": BaggingClassifier(DecisionTreeClassifier(),
                                 n_estimators=200, random_state=0),
    "Random forest": RandomForestClassifier(n_estimators=200,
                                            min_samples_leaf=3,
                                            random_state=0),
    "HistGradientBoosting": HistGradientBoostingClassifier(
        max_depth=3, learning_rate=0.05, random_state=0),
}
for nombre, modelo in modelos.items():
    vc = cross_val_score(con_preparacion(modelo), X_ent, y_ent, cv=5)
    print(f"{nombre:22s} {vc.mean():.3f} ± {vc.std():.3f}")

final = con_preparacion(modelos["Random forest"]).fit(X_ent, y_ent)
print(f"Exactitud en la prueba: {final.score(X_prueba, y_prueba):.3f}")

from sklearn.inspection import permutation_importance

imp = permutation_importance(final, X_prueba, y_prueba, n_repeats=20,
                             random_state=0)
orden = np.argsort(imp.importances_mean)
for j in orden[::-1]:
    print(f"{X.columns[j]:16s} {imp.importances_mean[j]:+.3f}")

colores = ["C1" if j == orden[-1] else "C0" for j in orden]
plt.barh(X.columns[orden], imp.importances_mean[orden], color=colores)
plt.xlabel("Exactitud que se pierde")
plt.title("Importancia")
plt.show()