Lectura_16_Intro_Metodos_de_Boosting.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "pHPi-jn1-Xy3"
   },
   "source": [
    "# ***Regresion e introduccion método de Boosting***\n",
    "\n",
    "## ***Universidad Central***\n",
    "\n",
    "## ***Maestría en analítica de datos***\n",
    "\n",
    "## ***Métodos estadísticos para analítica de datos.***\n",
    "\n",
    "## ***Docente: Luis Andrés Campos Maldonado.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 1377,
     "status": "ok",
     "timestamp": 1668721652833,
     "user": {
      "displayName": "Luis Andres Campos Maldonado",
      "userId": "15478348060554261231"
     },
     "user_tz": 300
    },
    "id": "PsggJ7H2MABe"
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import pandas as pd\n",
    "import numpy  as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pprint import pprint\n",
    "from sklearn.model_selection import RandomizedSearchCV, train_test_split\n",
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.tree import DecisionTreeRegressor, plot_tree\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor\n",
    "\n",
    "plt.style.use(\"ggplot\")\n",
    "plt.rcParams[\"figure.figsize\"] = (15,6)\n",
    "url_base = \"https://raw.githubusercontent.com/lacamposm/Metodos-Estadisticos/main/data/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "08GUU8AVM2T6"
   },
   "source": [
    "## _**[DecisionTreeRegressor](https://scikit-learn.org/stable/auto_examples/tree/plot_tree_regression.html)**_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qgUVKlq4NqAo"
   },
   "source": [
    "### ***Predicción de una valor continuo.***\n",
    "\n",
    "Predecir un valor continuo (también denominado regresión) con un árbol sigue la misma lógica y procedimiento que en una clasificación binaria, excepto que la impureza se mide por desviaciones al cuadrado de la media (errores cuadráticos) en cada subpartición, y el rendimiento predictivo se juzga por\n",
    "la raíz cuadrada del error cuadrático medio (RMSE) en cada partición. scikit-learn tiene la clase [`sklearn.tree.DecisionTreeRegressor`](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html) para entrenar un modelo de regresión del árbol de decisión."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nmUtbXYROajn"
   },
   "source": [
    "Cuando se usan los árboles, uno de los grandes obstáculos que se enfrentan los modeladores predictivos es la naturaleza percibida de \"caja negra\" de los métodos que utilizan. En este sentido, el modelo de Decision Tree tiene dos\n",
    "aspectos atractivos:\n",
    "\n",
    "• Los Decision Tree proporcionan una herramienta visual para explorar los datos, para tener una idea de que variables son importantes y cómo se relacionan entre sí. Los árboles pueden capturar relaciones no lineales entre variables predictoras.\n",
    "\n",
    "• Los Decision Tree proporcionan un conjunto de reglas que se pueden comunicar de manera efectiva a\n",
    "especialistas, ya sea para la implementación o para “vender” un proyecto de analitica de datos.\n",
    "\n",
    "Sin embargo, cuando se trata de predicción, aprovechar los resultados de múltiples árboles es típicamente más poderoso que usar un solo árbol. En particular, _Random Forrest_  casi siempre brindan una precisión predictiva superior y rendimiento, pero se pierden las mencionadas ventajas de un solo árbol."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jvS3zcUUOkiN"
   },
   "source": [
    "***Ideas importantes:***\n",
    "\n",
    "\n",
    "1. Los Decision Tree producen un conjunto de reglas para clasificar o predecir un resultado.\n",
    "\n",
    "2. Las reglas corresponden a particiones sucesivas de los datos en subparticiones.\n",
    "\n",
    "3. Cada partición, o división, hace referencia a un valor específico de una variable predictora y divide los datos en registros donde el valor del predictor está por encima o por debajo de ese valor dividido.\n",
    "\n",
    "4. En cada etapa, el algoritmo de Decision Tree seleciona la división que minimiza el resultado de impureza dentro de cada subpartición.\n",
    "\n",
    "5. La asiganación de la predicción se realiza siguiendo el camino que el árbol va indicando.\n",
    "\n",
    "6. Un árbol completamente desarrollado sobreajusta (_overfitting_) los datos y debe ser podado para que capture señal y no ruido.\n",
    "\n",
    "7. Los algoritmos como los _Random Forrest_ y árboles potenciados (_boosting_), ofrecen un mejor rendimiento predictivo, pero pierden el poder comunicativo basado en reglas de los algoritmos de Decision Tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 374
    },
    "executionInfo": {
     "elapsed": 445,
     "status": "ok",
     "timestamp": 1668721657603,
     "user": {
      "displayName": "Luis Andres Campos Maldonado",
      "userId": "15478348060554261231"
     },
     "user_tz": 300
    },
    "id": "G-EZUuCoPNRU",
    "outputId": "955f50d9-6411-4f63-b533-08ce4719f618"
   },
   "outputs": [],
   "source": [
    "# Simulemos unos datos con relación cuadrática.\n",
    "np.random.seed(123)\n",
    "m = 200\n",
    "X = 6 * np.random.rand(m, 1) - 3\n",
    "y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)\n",
    "df = pd.DataFrame(np.hstack((X, y)), columns=['regresor', 'y'])\n",
    "plt.scatter(X,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_reg = DecisionTreeRegressor(max_depth=2, random_state=123)\n",
    "tree_reg.fit(X, y)\n",
    "\n",
    "print(f\"R2 score in LR: {LinearRegression().fit(X, y).score(X, y):4f}\")\n",
    "print(f\"R2 score in tree model: {tree_reg.score(X, y):.4f}\")\n",
    "\n",
    "plot_tree(tree_reg, filled=True, feature_names=df.drop(columns=[\"y\"]).columns)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_DDkYgPhQP8z"
   },
   "source": [
    "Este árbol se parece mucho al árbol de clasificación que ya conocemos. La principal diferencia es que en lugar de predecir una clase en cada nodo, predice un valor. Como ejemplo, supongamos que deseamos hacer una predicción para $x_0=2.15$, luego se debe atravezar el árbol comenzando en la raíz, y finalmente llegando al nodo hoja y el `value` que se predice es 5.671. Esta predicción es simplemente el valor objetivo promedio de la variable target en 24 registros del set de train que quedarón allí. Además esta predicción da como media del error (MSE) (`squared_error`) 1.16 sobre estos 24 registros."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EOhZQptYS2uv"
   },
   "source": [
    "El algoritmo CART funciona casi de la misma manera que antes, excepto que en lugar de tratar de dividir el set de train de una manera que minimice la impureza, ahora trata de dividir el set de train de una manera que minimice el MSE. La función de costo que el algoritmo intenta minimizar es:\n",
    "\n",
    "$$J(k,t_k)=\\frac{m_{left}}{m}MSE_{left}+\\frac{m_{right}}{m}MSE_{rigth}$$\n",
    "\n",
    "Donde:\n",
    "\n",
    "1. $MSE_{node}=\\frac{1}{m_{mode}}\\displaystyle\\sum_{i\\in\\ \\ node}(\\hat{y}_{node}-y^{(i)})^{2}$\n",
    "\n",
    "2. $\\hat{y}_{node}=\\frac{1}{m_{node_i}}\\displaystyle\\sum_{i \\in \\ \\ node}y^{(i)}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_titSuPpU9bK"
   },
   "source": [
    "Al igual que para las tareas de clasificación, los árboles de decisión tienden a tener `overfitting` cuando se trata con tareas de regresión. Sin ninguna regularización (es decir, utilizando el valor predeterminado hiperparámetros), es bastante probable que su algoritmo este con `overfitting`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nodo Raiz\n",
    "print(\"\\nPara el nodo raiz:\")\n",
    "y_mean = df[\"y\"].mean()\n",
    "print(f\"value = {y_mean:.3f}\")\n",
    "print(f\"squared_error = {mean_squared_error(df['y'], [y_mean] * len(df)):.3f}\")\n",
    "print(f\"Numero de muestras= {len(df)}\")\n",
    "\n",
    "# Para el primer split\n",
    "print(\"\\nPara el primer split:\")\n",
    "y_left = df.query(\"regresor <= 1.507\")[[\"y\"]]\n",
    "y_right = df.query(\"regresor > 1.507\")[[\"y\"]]\n",
    "\n",
    "y_left_mean = np.mean(y_left)\n",
    "y_right_mean = np.mean(y_right)\n",
    "m_left = len(y_left)\n",
    "m_right = len(y_right)\n",
    "m_total = m_left + m_right\n",
    "\n",
    "print(f\"Instancias nodo izquiero = {m_left:.3f}\")\n",
    "print(f\"Instancias nodo derecho = {m_right:.3f}\")\n",
    "print(f\"Instancias nodo padre = {m_total:.3f}\")\n",
    "\n",
    "\n",
    "mse_left = mean_squared_error(y_left, [y_left_mean] * len(y_left))\n",
    "mse_right = mean_squared_error(y_right, [y_right_mean] * len(y_right))\n",
    "print(f\"\\nMSE nodo izquierdo = {mse_left:.3f}\")\n",
    "print(f\"MSE nodo derecho = {mse_right:.3f}\")\n",
    "\n",
    "mse_total = (m_left / m_total) * mse_left + (m_right / m_total) * mse_right\n",
    "\n",
    "print(f\"\\nMSE ponderado para los nodos izquierdo y derecho: {mse_total:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EFRnf-xlVfeT"
   },
   "source": [
    "***Nota.***\n",
    "\n",
    "Los árboles priorizan los límites de decisión ortogonales (todas las divisiones son perpendiculares a un eje), esto los hace sensibles a la rotación del set de train. Una forma de limitar problemas que cortes innecesarios (piense en rotar un región en 45 grados)\n",
    "es usar PCA, que a menudo resulta en una mejor orientación en el set de train."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uaV-BfqIfLyI"
   },
   "source": [
    "## ***Ejemplo 1.***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ***Modelos con DecisionTreeRegressor.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_house = pd.read_csv(url_base + \"kc_house_data.csv\")\n",
    "df_to_model = df_house.copy()\n",
    "df_to_model = df_to_model.drop(columns=[\"id\", \"sqft_living15\", \"sqft_lot15\", \"date\"])\n",
    "df_to_model[\"renovated\"] = df_to_model[\"yr_renovated\"] > 0\n",
    "df_to_model.drop(columns=[\"yr_renovated\"], inplace=True)\n",
    "df_to_model.drop_duplicates(inplace=True)\n",
    "df_to_model.info()\n",
    "df_to_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = df_to_model.drop(columns=[\"price\"])\n",
    "y = df_to_model[\"price\"]\n",
    "\n",
    "features_categoric = [\"waterfront\", \"zipcode\", \"view\", \"renovated\", \"grade\"]\n",
    "numerical_columns = [col for col in X.columns if col not in features_categoric]\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)\n",
    "\n",
    "preprocessor = ColumnTransformer(\n",
    "    transformers=[\n",
    "        (\"cat\", OneHotEncoder(handle_unknown=\"ignore\"), features_categoric),\n",
    "        (\"num\", \"passthrough\", numerical_columns)\n",
    "    ]\n",
    ")\n",
    "\n",
    "tree_model = DecisionTreeRegressor()\n",
    "\n",
    "pipeline_tree = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model\", tree_model)\n",
    "])\n",
    "\n",
    "pipeline_tree.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_model = pipeline_tree.named_steps[\"model\"]\n",
    "\n",
    "plt.figure(figsize=(15, 6))\n",
    "plot_tree(\n",
    "    tree_model,\n",
    "    feature_names=pipeline_tree.named_steps[\"preprocessor\"].get_feature_names_out(),\n",
    "    max_depth=3,\n",
    "    filled=True,\n",
    "    rounded=True,\n",
    "    fontsize=6\n",
    ")\n",
    "plt.title(\"Árbol de Decisión\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Overfitting!!!\")\n",
    "print(f\"R2 en Train: {pipeline_tree.score(X_train, y_train):.4f}\")\n",
    "print(f\"R2 en Test:  {pipeline_tree.score(X_test, y_test):.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Buscando los \"mejores\" parámetros para Decision Tree\n",
    "params = {\n",
    "    \"model__max_depth\": [n + 2 for n in range(9)],\n",
    "    \"model__criterion\": [\"squared_error\", \"absolute_error\", \"friedman_mse\"] ,\n",
    "    \"model__min_samples_split\" : [10 + 10 * n for n in range(30)]\n",
    "}\n",
    "\n",
    "random_search_tree = RandomizedSearchCV(\n",
    "    pipeline_tree,\n",
    "    param_distributions=params,\n",
    "    cv=5,\n",
    "    n_jobs=-1,\n",
    "    random_state=123,\n",
    ")\n",
    "\n",
    "search_tree = random_search_tree.fit(X_train, y_train)\n",
    "\n",
    "print(f\"'Best' params:\")\n",
    "pprint(search_tree.best_params_)\n",
    "print(f'\\nR^2 en TRAIN: {search_tree.best_estimator_.score(X_train, y_train):.4f}')\n",
    "print(f'R^2 en TEST:  {search_tree.best_estimator_.score(X_test, y_test):.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Q4vZbUGh-uYj"
   },
   "source": [
    "### ***Modelo con RandomForestRegressor.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 mins corriendo\n",
    "rf_model = RandomForestRegressor()\n",
    "\n",
    "pipeline_rf = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model\", rf_model)\n",
    "])\n",
    "\n",
    "pipeline_rf.fit(X_train, y_train)\n",
    "\n",
    "print(\"Overfitting!!!\")\n",
    "print(f\"R2 en Train: {pipeline_rf.score(X_train, y_train):.4f}\")\n",
    "print(f\"R2 en Test:  {pipeline_rf.score(X_test, y_test):.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 min corriendo\n",
    "rf_model1 = RandomForestRegressor(\n",
    "    max_depth=12,\n",
    "    min_samples_split=50,\n",
    "    n_estimators=150    \n",
    ")\n",
    "\n",
    "pipeline_rf1 = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model\", rf_model1)\n",
    "])\n",
    "\n",
    "pipeline_rf1.fit(X_train, y_train)\n",
    "\n",
    "print(pipeline_rf1.score(X_train, y_train))\n",
    "print(pipeline_rf1.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Buscando los \"mejores\" parámetros para Random Forrest\n",
    "base_estimator_rf = RandomForestRegressor()\n",
    "\n",
    "pipeline_rf_best = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model\", base_estimator_rf)\n",
    "])\n",
    "\n",
    "\n",
    "params = {\n",
    "    \"model__max_depth\": [4, 5, 6, 8, 10], \n",
    "    \"model__criterion\": [\"squared_error\", \"absolute_error\", \"friedman_mse\"],\n",
    "    \"model__min_samples_split\": [50 + 25 * n for n in range(10)],\n",
    "    \"model__max_features\": [\"sqrt\", \"log2\", 10],\n",
    "    \"model__n_estimators\": [100 + 50 * n for n in range(10)]\n",
    "}\n",
    "      \n",
    "search_regressor_rf = RandomizedSearchCV(\n",
    "    pipeline_rf_best,\n",
    "    params,\n",
    "    n_iter=20,\n",
    "    cv=5,\n",
    "    n_jobs=-1,\n",
    "    verbose=3,\n",
    "    return_train_score=True\n",
    ")\n",
    "\n",
    "search_rf = search_regressor_rf.fit(X_train, y_train)\n",
    "\n",
    "print(\"'Best' params:\")\n",
    "pprint(search_rf.best_params_)\n",
    "print(f\"\\nR2 en Train: {search_rf.best_estimator_.score(X_train, y_train):.4f}\")\n",
    "print(f\"R2 en Test:  {search_rf.best_estimator_.score(X_test, y_test):.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Guardamos resultamos para analisis posteriores.\n",
    "pd.DataFrame(search_regressor_rf.cv_results_).to_csv(\"data/rf_cv_grid_results.csv\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_model_best = RandomForestRegressor(\n",
    "    n_estimators=200,\n",
    "    max_depth=10,\n",
    "    max_features=15,\n",
    "    min_samples_split=50\n",
    ")\n",
    "\n",
    "rf_pipeline = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model\", rf_model_best)\n",
    "])\n",
    "\n",
    "rf_pipeline.fit(X_train, y_train)\n",
    "\n",
    "print(f\"R2 en Train: {rf_pipeline.score(X_train, y_train):.4f}\")\n",
    "print(f\"R2 en Test:  {rf_pipeline.score(X_test, y_test):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Ejercicio***\n",
    "\n",
    "Reconstruya el mismo modelo con regresion lineal multiple, compare metricas, discuta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sKzU2x9zkMj7"
   },
   "source": [
    "# ***Boosting. (Impulsar)***\n",
    "\n",
    "Los modelos de ensemble son una herramienta estándar para el modelado predictivo. ***Boosting*** es un técnica general para crear un modelo de un conjunto de modelos. Como ***bagging***, ***boosting*** se usa más comúnmente con árboles de decisión. A pesar de sus similitudes, ***boosting*** toma un enfoque muy diferente. Como resultado, mientras que ***bagging*** se puede hacer con relativamente poca afinación, ***boosting*** requiere mucho mayor cuidado en su aplicación. Si estos dos métodos fueran coches, ***bagging*** podría considerarse un Honda Accord (confiable y estable), mientras que impulsar podría ser considerado un Porsche (potente pero requiere más cuidados).\n",
    "\n",
    "En los modelos de regresión lineal, los residuales se examinan para buscar mejorar el ajuste. ***Boosting*** toma este concepto y va mucho más allá, a diferencia del ***bagging***, que promedia predicciones independientes de modelos, ***boosting*** construye modelos secuenciales que corrigen los errores del anterior. El Boosting es un enfoque iterativo que combina múltiples modelos base (habitualmente débiles) para construir un modelo de ***ensemble*** fuerte.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ***Planteamiento del problema de optimización que se va a resolver con Gradient Boosting***\n",
    "\n",
    "#### ***Contexto***\n",
    "\n",
    "En el contexto de la **regresión**, nuestro objetivo es encontrar una función $F(x)$ que, dada una entrada $x$, prediga un valor $y$ lo más cercano posible al valor real observado. En términos de ***optimización***, queremos ajustar esta función para que minimice el error entre las predicciones y los valores reales.\n",
    "\n",
    "Matemáticamente, tenemos:\n",
    "\n",
    "$$\\hat{y_i} = F(x_i)$$\n",
    "\n",
    "Donde:\n",
    "\n",
    "- $y_i$ es el valor real observado.\n",
    "- $\\hat{y_i} = F(x_i)$ es la predicción de la función $F(x)$ para el punto $x_i$.\n",
    "\n",
    "El objetivo es ***minimizar la diferencia*** entre las predicciones $F(x_i)$ y los valores reales $y_i$, lo que se puede expresar como:\n",
    "\n",
    "$$\\text{Minimizar} \\ \\mathcal{L}(F(x), y) = \\frac{1}{n} \\sum_{i=1}^n (y_i - F(x_i))^2$$\n",
    "\n",
    "Este es el problema de optimización clásico de la ***minimización del Error Cuadrático Medio (MSE)***.\n",
    "\n",
    "#### ***¿Cómo resolver el problema de optimización?***\n",
    "\n",
    "Una forma directa de resolver este problema es utilizar un ***algoritmo de optimización*** que actualice las predicciones de manera iterativa para minimizar la función de pérdida.\n",
    "\n",
    "El ***Gradient Boosting*** resuelve este problema ***iterativamente***. Ajusta una ***serie de modelos*** (típicamente árboles de decisión) que, en cada paso, ***corrigen los errores*** cometidos por el modelo previo. En lugar de ajustar todos los parámetros de una vez (como en la regresión lineal), ***Gradient Boosting*** ajusta pequeñas correcciones a las predicciones de manera secuencial, optimizando en cada iteración para reducir el error en la dirección correcta.\n",
    "\n",
    "#### ***La derivada de la función de pérdida***\n",
    "\n",
    "\n",
    "El gradiente es una medida de la pendiente de la función de pérdida y nos indica la dirección en la que debemos movernos para reducir el error. Asi, para resolver el problema de optimización, necesitamos calcular el gradiente de la función de pérdida con respecto a las predicciones $F(x)$. Esto nos permitirá saber cómo actualizar nuestras predicciones en cada iteración del algoritmo de ***Gradient Boosting***.\n",
    "\n",
    "En este caso, la función de pérdida que estamos utilizando es el **Error Cuadrático Medio (MSE)**:\n",
    "\n",
    "$$\\mathcal{L}(F(x), y) = \\frac{1}{n} \\sum_{i=1}^n (y_i - F(x_i))^2$$\n",
    "\n",
    "Donde:\n",
    "\n",
    "- $y_i$ es el valor real observado para el punto $x_i$.\n",
    "- $F(x_i)$ es la predicción que estamos haciendo para el punto $x_i$.\n",
    "\n",
    "El objetivo es **minimizar** esta función de pérdida con respecto a las predicciones $F(x_i)$. Para eso, necesitamos **derivar** la función de pérdida con respecto a las predicciones $F(x_i)$ y obtener el gradiente, que nos indica la dirección en la que debemos mover nuestras predicciones para reducir el error.\n",
    "\n",
    "#### ***Derivada de MSE***\n",
    "\n",
    "Para calcular la derivada de la MSE con respecto a $F(x_i)$, aplicamos las reglas de derivación:\n",
    "\n",
    "$$\\frac{\\partial}{\\partial F(x_i)} \\mathcal{L}(F(x), y) = \\frac{\\partial}{\\partial F(x_i)} \\left( \\frac{1}{n} \\sum_{i=1}^n (y_i - F(x_i))^2 \\right)$$\n",
    "\n",
    "Al derivar la expresión, obtenemos:\n",
    "\n",
    "$$\\frac{\\partial}{\\partial F(x_i)} \\mathcal{L}(F(x), y) \\varpropto - (y_i - F(x_i))$$\n",
    "\n",
    "#### ***Interpretación del gradiente***\n",
    "\n",
    "El gradiente de la función de pérdida en este caso es el **residuo** entre el valor real $y_i$ y la predicción $F(x_i)$:\n",
    "\n",
    "$$\\text{Gradiente} \\varpropto -(y_i - F(x_i))$$\n",
    "\n",
    "Este valor nos indica cuánto y en qué dirección debemos ajustar nuestras predicciones $F(x_i)$ para minimizar el error. Si el gradiente es positivo, eso significa que la predicción $F(x_i)$ es menor que $y_i$, por lo que debemos aumentar la predicción. Si el gradiente es negativo, eso significa que la predicción $F(x_i)$ es mayor que $y_i$, por lo que debemos disminuirla.\n",
    "\n",
    "En **Gradient Boosting**, utilizamos este gradiente para **ajustar nuestras predicciones** en cada iteración, corrigiendo los errores de los modelos anteriores.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ***Algoritmo de Gradient Boosting***\n",
    "\n",
    "#### ***Pasos del algoritmo***\n",
    "\n",
    "1. **Inicialización**\n",
    "   - Definir el número máximo de iteraciones $M$.\n",
    "   - Inicializar el modelo del _ensemble_ $\\hat{F}_0(x)$ como una constante, que típicamente es la media de las observaciones de $y$, es decir, $\\hat{F}_0(x) = \\text{media}(y)$.\n",
    "   - Inicializar los residuos $r_i^{(0)} = y_i - \\hat{F}_0(x_i)$, que son las diferencias entre los valores reales y las predicciones iniciales.\n",
    "\n",
    "2. **Iteración**  \n",
    "   Para $m = 1, 2, \\dots, M$:\n",
    "\n",
    "   - **Entrenar un modelo base** $\\hat{f}_m(x)$:\n",
    "     - Entrenar un modelo base $\\hat{f}_m(x)$ (por ejemplo, un árbol de decisión) para predecir los residuos o gradientes. Este modelo debe minimizar la función de pérdida en los residuos, de modo que:\n",
    "     \n",
    "       $$\\hat{f}_m(x) = \\text{modelo base ajustado sobre los residuos}$$\n",
    "     - En lugar de ajustarse directamente sobre los valores originales de $y$, el modelo base se ajusta sobre los residuos, que son el gradiente de la función de pérdida con respecto a las predicciones actuales $\\hat{F}_{m-1}(x)$.\n",
    "\n",
    "   - **Actualizar el modelo del _ensemble_**:\n",
    "     - Actualizar el modelo $\\hat{F}_m(x)$ combinando la predicción del modelo base $\\hat{f}_m(x)$ con la predicción anterior $\\hat{F}_{m-1}(x)$ utilizando un coeficiente de aprendizaje $\\nu$:\n",
    "     \n",
    "       $$\\hat{F}_m(x) = \\hat{F}_{m-1}(x) + \\nu \\cdot \\hat{f}_m(x)$$\n",
    "     - Aquí, $\\nu$ es el coeficiente de aprendizaje que controla la magnitud de la actualización en cada paso.\n",
    "\n",
    "   - **Actualizar los residuos**:\n",
    "     - Los nuevos residuos se actualizan como las diferencias entre los valores reales $y_i$ y las nuevas predicciones $\\hat{F}_m(x_i)$, es decir:\n",
    "\n",
    "       $$r_i^{(m)} = y_i - \\hat{F}_m(x_i)$$\n",
    "     - Los residuos son utilizados para la próxima iteración, donde se ajusta un nuevo modelo base $\\hat{f}_{m+1}(x)$ sobre estos residuos.\n",
    "\n",
    "3. **Determinación del criterio de parada**\n",
    "   - Si se alcanza $m = M$ o si no hay mejora significativa en la función de pérdida entre iteraciones, el proceso de iteración se detiene.\n",
    "\n",
    "4. **Predicción final**\n",
    "   - La predicción final del modelo de Boosting es la suma ponderada de los modelos base entrenados en cada iteración:\n",
    "     $$\\hat{F}(x) = \\hat{F}_0(x) + \\sum_{m=1}^M \\nu \\cdot \\hat{f}_m(x)\n",
    "     $$\n",
    "   - Aquí, $\\hat{F}_0(x)$ es el modelo inicial (la media en el caso de regresión), y $\\nu \\cdot \\hat{f}_m(x)$ son las contribuciones de los modelos base entrenados en cada iteración.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iK70gAFeJumH"
   },
   "source": [
    "### ***Modelo con GradientBoosting.***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_house = pd.read_csv(url_base + \"kc_house_data.csv\")\n",
    "df_to_model = df_house.copy()\n",
    "df_to_model = df_to_model.drop(columns=[\"id\", \"sqft_living15\", \"sqft_lot15\", \"date\"])\n",
    "df_to_model[\"renovated\"] = df_to_model[\"yr_renovated\"] > 0\n",
    "df_to_model.drop(columns=[\"yr_renovated\"], inplace=True)\n",
    "df_to_model.drop_duplicates(inplace=True)\n",
    "\n",
    "X = df_to_model.drop(columns=[\"price\"])\n",
    "y = df_to_model[\"price\"]\n",
    "\n",
    "features_categoric = [\"waterfront\", \"zipcode\", \"view\", \"renovated\", \"grade\"]\n",
    "numerical_columns = [col for col in X.columns if col not in features_categoric]\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)\n",
    "\n",
    "preprocessor = ColumnTransformer(\n",
    "    transformers=[\n",
    "        (\"cat\", OneHotEncoder(handle_unknown=\"ignore\"), features_categoric),\n",
    "        (\"num\", \"passthrough\", numerical_columns)\n",
    "    ]\n",
    ")\n",
    "\n",
    "gboost = GradientBoostingRegressor()\n",
    "\n",
    "pipeline_gboost = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model_gdb\", gboost)\n",
    "])\n",
    "\n",
    "pipeline_gboost.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5404,
     "status": "ok",
     "timestamp": 1668721262365,
     "user": {
      "displayName": "Luis Andres Campos Maldonado",
      "userId": "15478348060554261231"
     },
     "user_tz": 300
    },
    "id": "SpkTjLbpJuF1",
    "outputId": "4bf86272-1f63-4dcd-80b6-a0939897ec0e"
   },
   "outputs": [],
   "source": [
    "print(pipeline_gboost.score(X_train, y_train))\n",
    "print(pipeline_gboost.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ***Desglose del Algoritmo***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Primeros tres arboles.\n",
    "first_tree = pipeline_gboost.named_steps[\"model_gdb\"].estimators_[0, 0]\n",
    "second_tree = pipeline_gboost.named_steps[\"model_gdb\"].estimators_[1, 0]\n",
    "third_tree = pipeline_gboost.named_steps[\"model_gdb\"].estimators_[2, 0]\n",
    "\n",
    "# Predicciones acumuladas en Train\n",
    "staged_predictions_train = list(pipeline_gboost.named_steps[\"model_gdb\"].staged_predict(\n",
    "    pipeline_gboost.named_steps[\"preprocessor\"].transform(X_train)\n",
    "))\n",
    "\n",
    "# Predicciones acumuladas en Test\n",
    "staged_predictions_test = pipeline_gboost.named_steps[\"model_gdb\"].staged_predict(\n",
    "    pipeline_gboost.named_steps[\"preprocessor\"].transform(X_test)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Estimador basico inicial: la media de y_train***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Primer estimador la media de y_train (paso 0).\n",
    "y_mean = y_train.mean()\n",
    "residuals_initial = y_train - y_mean\n",
    "print(\"\\nResiduos iniciales:\")\n",
    "print(residuals_initial)\n",
    "mse_inicial = mean_squared_error(y_train, residuals_initial)\n",
    "print(f\"\\nMSE inicial = {mse_inicial:,.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Mejorardo la media como primer estimador con un arbol***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Primer Arbol\n",
    "plot_tree(first_tree, filled=True, feature_names=preprocessor.get_feature_names_out())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predicciones del primer arbol estimando los residuos.\n",
    "first_tree.predict(preprocessor.transform(X_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Actualizamos valores\n",
    "predictions_stage_1 = [y_mean]*len(X_train) + 0.1 * first_tree.predict(preprocessor.transform(X_train))\n",
    "display(predictions_stage_1)\n",
    "mse_paso_1 =  mean_squared_error(y_train, predictions_stage_1)\n",
    "print(f\"\\nMSE paso 1 = {mse_paso_1:,.3f}\")\n",
    "print(mse_inicial > mse_paso_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Segundo arbol***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Segundo Arbol\n",
    "plot_tree(second_tree, filled=True, feature_names=preprocessor.get_feature_names_out())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Elementos del arbol 2.\n",
    "y_step2 = y_train -  predictions_stage_1 # Nuevos residuos a modelar\n",
    "predict_tree_2 = second_tree.predict(preprocessor.transform(X_train))\n",
    "# Segundo arbol prediciendo los nuevos residuos y actualizando prediccion\n",
    "predictions_stage_2  = predictions_stage_1 + 0.1 * second_tree.predict(preprocessor.transform(X_train))\n",
    "display(predictions_stage_2)  # coincide con: staged_predictions_train[2]\n",
    "mse_paso_2 = mean_squared_error(y_train, predictions_stage_2)\n",
    "print(f\"\\nMSE paso 2 = {mse_paso_2:,.3f}\")\n",
    "print(mse_paso_1 > mse_paso_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Tercer arbol***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tercer Arbol\n",
    "plot_tree(third_tree, filled=True, feature_names=preprocessor.get_feature_names_out())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Elementos del arbol 3.\n",
    "y_step3 = y_train -  predictions_stage_2 # Nuevos residuos a modelar\n",
    "predict_tree_3 = third_tree.predict(preprocessor.transform(X_train))\n",
    "# Segundo arbol prediciendo los nuevos residuos y actualizando prediccion\n",
    "predictions_stage_3  = predictions_stage_2 + 0.1 * third_tree.predict(preprocessor.transform(X_train))\n",
    "display(predictions_stage_3)  # coincide con: staged_predictions_train[1]\n",
    "mse_paso_3 =  mean_squared_error(y_train, predictions_stage_3)\n",
    "print(f\"\\nMSE paso 3 = {mse_paso_3:,.3f}\")\n",
    "print(mse_paso_2 > mse_paso_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errors_train = np.zeros(100)\n",
    "errors_test = np.zeros(100)\n",
    "\n",
    "for i, (pred_train, pred_test) in enumerate(zip(staged_predictions_train, staged_predictions_test)):\n",
    "    errors_train[i] = mean_squared_error(y_train, pred_train)\n",
    "    errors_test[i] = mean_squared_error(y_test, pred_test)\n",
    "\n",
    "plt.plot(range(1, 101), errors_train, marker=\"o\", label=\"Train\")\n",
    "plt.plot(range(1, 101), errors_test, marker=\"o\", label=\"Test\")\n",
    "plt.title(\"Error en el conjunto de entrenamiento y prueba a través de las iteraciones\")\n",
    "plt.xlabel(\"Número de Iteraciones (Árboles)\")\n",
    "plt.ylabel(\"Error Cuadrático Medio\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ***Busqueda de Hiperparametros en GradientBoostingRegressor***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preprocessor = ColumnTransformer(\n",
    "    transformers=[\n",
    "        (\"cat\", OneHotEncoder(handle_unknown=\"ignore\"), features_categoric),\n",
    "        (\"num\", \"passthrough\", numerical_columns)\n",
    "    ]\n",
    ")\n",
    "\n",
    "gboost = GradientBoostingRegressor()\n",
    "\n",
    "pipeline_gboost = Pipeline(steps=[\n",
    "    (\"preprocessor\", preprocessor),\n",
    "    (\"model_gdb\", gboost)\n",
    "])\n",
    "\n",
    "pipeline_gboost.fit(X_train, y_train)\n",
    "\n",
    "param_grid = {\n",
    "    \"model_gdb__n_estimators\": [100, 200, 300, 350, 400],\n",
    "    \"model_gdb__learning_rate\": [0.075, 0.01, 0.02, 0.05, 0.1],\n",
    "    \"model_gdb__max_depth\": [3, 4, 5, 6, 8, 10],\n",
    "    \"model_gdb__min_samples_split\": [2, 5, 8, 10],\n",
    "    \"model_gdb__min_samples_leaf\": [50, 100, 150],\n",
    "    \"model_gdb__max_features\": [\"sqrt\", \"log2\", 15]\n",
    "}\n",
    "\n",
    "grid_search_gboost = RandomizedSearchCV(\n",
    "    pipeline_gboost,\n",
    "    param_grid,\n",
    "    n_iter=10,\n",
    "    cv=5,\n",
    "    n_jobs=-1,\n",
    "    verbose=3,\n",
    "    return_train_score=True\n",
    ")\n",
    "\n",
    "grid_search_gboost.fit(X_train, y_train)\n",
    "print(grid_search_gboost.best_params_)\n",
    "print(grid_search_gboost.best_estimator_.score(X_train, y_train))\n",
    "print(grid_search_gboost.best_estimator_.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Guardamos resultamos para analisis posteriores.\n",
    "pd.DataFrame(grid_search_gboost.cv_results_).to_csv(\"data/gboost_grid_cv_search.csv\", index=False)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyOV0WwDJiUoZN5zVRTVcyky",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}