2023/seminars/seminar07/sem07_boosting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fitted-subscriber",
   "metadata": {
    "id": "fitted-subscriber"
   },
   "source": [
    "# Семинар 7. Бустинги"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "canadian-design",
   "metadata": {
    "id": "canadian-design"
   },
   "source": [
    "### 1. Бэггинг и бустинг: идеи подходов"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "patient-watson",
   "metadata": {
    "id": "patient-watson"
   },
   "source": [
    "Основная идея **бэггинга** состоит в том, чтобы агрегировать ответы большого числа независимых переобученных моделей, имеющих низкое смещение и большую дисперсию. Агрегация позволяет добиться значительного снижения дисперсии при сохранении низкого смещения. Бэггинг над деревьями, в котором при построении каждого дерева при разбиениях используются случайные подмножества признаков, называется **случайным лесом**.\n",
    "\n",
    "Идея бустинга заключается в последовательном построении алгоритмов, каждый из которых учитывает ошибки построенной до сих пор композиции."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff9d0e7d",
   "metadata": {
    "id": "ff9d0e7d"
   },
   "source": [
    "$$\n",
    "    \\sum_{i = 1}^{\\ell}\n",
    "        L(y_i, a_{N - 1}(x_i) + \\gamma_N b_N(x_i))\n",
    "    \\to\n",
    "    \\min_{b_N, \\gamma_N}\n",
    "$$\n",
    "\n",
    "\n",
    "**Вопросики:**\n",
    "- Что происходит на шаге $0$?\n",
    "- Что происходит на шаге $N$?\n",
    "\n",
    "$$ s_i\n",
    "    =\n",
    "    -\n",
    "    \\left.\n",
    "    \\frac{\\partial L(y_i, z)}{\\partial z}\n",
    "    \\right|_{z = a_{N - 1}(x_i)}\n",
    "$$\n",
    "\n",
    "$$\n",
    "b_N(x)\n",
    "    =\n",
    "    argmin\n",
    "        \\sum_{i = 1}^{l}\n",
    "            \\left(\n",
    "                b(x_i) - s_i\n",
    "            \\right)^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "working-illustration",
   "metadata": {
    "id": "working-illustration"
   },
   "source": [
    "#### 1.1 Простой пример.\n",
    "\n",
    "Пусть мы решаем задачу регрессии на тренировочной выборке $(X, y)$. Мы обучаем модель $a(X)$, которая является композицией базовых моделей $b_n(X)$. Будем считать, что композиция строится простым суммированием ответов базовых моделей:\n",
    "\n",
    "$$\n",
    "a(X) = \\sum_{n = 1}^{k} b_n(X).\n",
    "$$\n",
    "\n",
    "В качестве базовой модели будем использовать решающее дерево глубины один (decision stump).\n",
    "\n",
    "Будем решать задачу при помощи минимизации среднеквадратичной ошибки:\n",
    "\n",
    "$$\n",
    "\\dfrac{1}{\\ell}\\sum_{i = 1}^{\\ell} (a(x_i) - y_i)^2 \\to \\min_a.\n",
    "$$\n",
    "\n",
    "В этом примере будем работать только с обучающей выборкой."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "developed-salad",
   "metadata": {
    "id": "developed-salad"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import mean_squared_error, mean_absolute_error\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adjustable-absolute",
   "metadata": {
    "id": "adjustable-absolute"
   },
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "N = 100\n",
    "X = np.linspace(0, 1*np.pi, N).reshape(-1, 1)\n",
    "y = np.sin(X)[:, 0] + np.random.normal(0, 0.1, size = N)\n",
    "\n",
    "# Функция для визуализации выборки и предсказаний\n",
    "def plot_sample_model(X, y, plot_predictions=False, y_pred=None, y_pred_label=None, loss='mse'):\n",
    "    plt.figure(figsize = (10, 7))\n",
    "    plt.scatter(X, y, label = 'Обучающая выборка', alpha = 0.7)\n",
    "    if plot_predictions:\n",
    "        plt.plot(X, y_pred, label = y_pred_label, c = 'r')\n",
    "        if loss == 'mse':\n",
    "            plt.title('MSE: ' + str(mean_squared_error(y, y_pred)))\n",
    "        elif loss == 'mae':\n",
    "            plt.title('MAE: ' + str(mean_absolute_error(y, y_pred)))\n",
    "    plt.xlabel('X')\n",
    "    plt.ylabel('y')\n",
    "    plt.legend();\n",
    "\n",
    "plot_sample_model(X, y, plot_predictions = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sustained-logistics",
   "metadata": {
    "id": "sustained-logistics"
   },
   "source": [
    "**Шаг 0.** Имеем пустую композицию решающих деревьев $a(X) = \\{\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "jewish-sharing",
   "metadata": {
    "id": "jewish-sharing"
   },
   "outputs": [],
   "source": [
    "a = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "naughty-interim",
   "metadata": {
    "id": "naughty-interim"
   },
   "source": [
    "**Шаг 1.** Обучим первое решающее дерево $DT_1$ и включим его в композицию: $a^{(1)}(X) = DT_1$. Получим предсказания композиции на обучающей выборке."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d22a747",
   "metadata": {
    "id": "5d22a747"
   },
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeRegressor\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "coastal-thesis",
   "metadata": {
    "id": "coastal-thesis"
   },
   "outputs": [],
   "source": [
    "# TODO: обучите решающее дерево глубины 1 и сделайте предсказания на обучающей выборке\n",
    "...\n",
    "\n",
    "# TODO: включите предсказания первого дерева в композицию (простым суммированием)\n",
    "a = ...\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a, y_pred_label = 'Предсказания a = DT1')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12c7c794",
   "metadata": {
    "id": "12c7c794"
   },
   "source": [
    "**Шаг 2.** В качестве ошибок композиции будем использовать **сдвиги**, в случае MSE — это то же самое, что и остатки или ошибки\n",
    "\n",
    "\n",
    "**Задание:** вычислите градиент $MSE$\n",
    "\n",
    "$$\n",
    "L(y, p) = \\dfrac{1}{2} \\sum_{i = 1}^{\\ell} (y_i - p_i)^2\n",
    "$$\n",
    "\n",
    "по $p_i$ в точке $p_i = a(x_i)$.\n",
    "\n",
    "**Решение:**\n",
    "\n",
    "$$\n",
    "\\dfrac{\\partial L(y, p)}{\\partial p_i} = -(y_i - p_i) = \\{p_i = a(x_i)\\} = -(y_i - a(x_i)).\n",
    "$$\n",
    "\n",
    "Таким образом, антиградиент функции потерь равен:\n",
    "\n",
    "$$\n",
    "s_i = - (-(y_i - a(x_i))) = y_i - a(x_i),\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "familiar-sessions",
   "metadata": {
    "id": "familiar-sessions"
   },
   "outputs": [],
   "source": [
    "# TODO: рассчитайте сдвиги\n",
    "s1 = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "improved-brake",
   "metadata": {
    "id": "improved-brake"
   },
   "source": [
    "**Шаг 3.** Обучим второе решающее дерево $DT_2$, предсказывающее сдвиги. Добавим предсказания второго дерева в композицию: $a^{(2)}(X) = DT_1 + DT_2$. Получим предсказания композиции на обучающей выборке."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "knowing-albert",
   "metadata": {
    "id": "knowing-albert"
   },
   "outputs": [],
   "source": [
    "# TODO: обучите второе решающее дерево глубины 1, предсказывающее остатки\n",
    "...\n",
    "\n",
    "# TODO: сделайте предсказания на обучающей выборке\n",
    "...\n",
    "\n",
    "# TODO: включите предсказания второго дерева в композицию (простым суммированием)\n",
    "a += ...\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a, y_pred_label = 'Предсказания a = DT1 + DT2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "funded-extreme",
   "metadata": {
    "id": "funded-extreme"
   },
   "source": [
    "Как мы видим, решающая поверхность стала более сложной, и теперь она более точно приближает обучающую выборку."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "appropriate-passion",
   "metadata": {
    "id": "appropriate-passion"
   },
   "source": [
    "**Шаги 4 – ...** Повторяем шаги 2-3, пока не надоест."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "saved-longer",
   "metadata": {
    "id": "saved-longer"
   },
   "outputs": [],
   "source": [
    "# TODO: рассчитайте остатки\n",
    "s2 = ...\n",
    "\n",
    "# TODO: обучите третье решающее дерево глубины 1, предсказывающее остатки\n",
    "...\n",
    "\n",
    "# TODO: сделайте предсказания на обучающей выборке\n",
    "...\n",
    "\n",
    "# TODO: включите предсказания третьего дерева в композицию (простым суммированием)\n",
    "a += ...\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a, y_pred_label = 'Предсказания DT1 + DT2 + DT3')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "certified-applicant",
   "metadata": {
    "id": "certified-applicant"
   },
   "source": [
    "Как мы видим, при добавлении базовых моделей решающая поверхность становится более сложной и всё точнее приближает обучающую выборку."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "integrated-northeast",
   "metadata": {
    "id": "integrated-northeast"
   },
   "source": [
    "**Важный момент:** заметим, что при обучении очередной базовой модели композиция предыдущего шага считается фиксированной. Это означает, что обучение новой базовой модели не влияет на уже обученные модели, содержащиеся в композиции."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "average-tumor",
   "metadata": {
    "id": "average-tumor"
   },
   "source": [
    "### 2. Функции потерь градиентного бустинга для регрессии и классификации"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "becoming-causing",
   "metadata": {
    "id": "becoming-causing"
   },
   "source": [
    "#### 2.1 Регрессия."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "referenced-feeding",
   "metadata": {
    "id": "referenced-feeding"
   },
   "source": [
    "В задачах регрессии в качестве целевой функции для композиции обычно берутся:\n",
    "\n",
    "1. $MSE$, которую мы рассматривали ранее.\n",
    "\n",
    "2. $MAE$:\n",
    "\n",
    "$$\n",
    "MAE(y, p) = \\sum_{i = 1}^{\\ell} |y_i - p_i|\n",
    "$$\n",
    "\n",
    "Попробуем сделать шаг градиентного бустинга, используя в качестве функции потерь $MAE$.\n",
    "\n",
    "**Задание:** рассчитайте $s_i$ для $MAE$.\n",
    "\n",
    "**Решение:** $s_i^{(N)} = -\\mathrm{sign}(a^{(N-1)}(x_i) - y_i)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "comic-above",
   "metadata": {
    "id": "comic-above"
   },
   "outputs": [],
   "source": [
    "# Инициализация\n",
    "a = 0\n",
    "\n",
    "# Обучение решающего дерева глубины 1 и получение предсказаний на обучающей выборке\n",
    "dt1 = DecisionTreeRegressor(max_depth = 1)\n",
    "dt1.fit(X, y)\n",
    "dt1_pred = dt1.predict(X)\n",
    "\n",
    "# Включение предсказаний первого дерева в композицию (простым суммированием)\n",
    "a = dt1_pred\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a, y_pred_label = 'Предсказания a = DT1', loss = 'mae')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "breeding-warehouse",
   "metadata": {
    "id": "breeding-warehouse"
   },
   "outputs": [],
   "source": [
    "# TODO: рассчитайте антиградиент\n",
    "s1 = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "reliable-actor",
   "metadata": {
    "id": "reliable-actor"
   },
   "outputs": [],
   "source": [
    "# Обучение второго решающего дерева глубины 1, предсказывающего антиградиент\n",
    "dt2 = DecisionTreeRegressor(max_depth = 1)\n",
    "dt2.fit(X, s1)\n",
    "\n",
    "# Получение предсказаний на обучающей выборке\n",
    "dt2_pred = dt2.predict(X)\n",
    "\n",
    "# Включение предсказаний второго дерева в композицию (простым суммированием)\n",
    "a += dt2_pred\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a,\n",
    "                  y_pred_label = 'Предсказания a = DT1 + DT2', loss = 'mae')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "needed-farmer",
   "metadata": {
    "id": "needed-farmer"
   },
   "source": [
    "Заметим, что $MAE$ увеличилась, а решающая поверхность имеет \"виток\", уходящий за пределы обучающей выборки. Эту проблему мы решим позже."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acting-prerequisite",
   "metadata": {
    "id": "acting-prerequisite"
   },
   "source": [
    "#### 2.2 Классификация."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "explicit-breakdown",
   "metadata": {
    "id": "explicit-breakdown"
   },
   "source": [
    "В задачах бинарной классификации обычно используется логистическая функция потерь, с которой мы уже сталкивались в логистической регрессии."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modular-southwest",
   "metadata": {
    "id": "modular-southwest"
   },
   "source": [
    "$$\n",
    "L(y, p) = \\log(1 + e^{-yp}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "allied-consistency",
   "metadata": {
    "id": "allied-consistency"
   },
   "source": [
    "Плиз, умейте показывать, что в этом случае антиградиенты можно вычислить по формуле\n",
    "\n",
    "$$\n",
    "s^{(N)}_i = \\dfrac{y_i}{1 + e^{y_ia^{(N-1)}(x_i)}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adult-operation",
   "metadata": {
    "id": "adult-operation"
   },
   "source": [
    "Попробуем сделать шаг градиентного бустинга в задаче бинарной классификации с логистической функцией потерь."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12227b2e-3a30-4c52-b882-3e6581e311a4",
   "metadata": {
    "id": "12227b2e-3a30-4c52-b882-3e6581e311a4"
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_moons, load_iris # import function from the library\n",
    "\n",
    "N = 1000\n",
    "\n",
    "X, y = make_moons(n_samples=N, noise=0.1, random_state=11) # generate data sample\n",
    "y = y * 2 - 1 # {-1, +1}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5b22299-b22a-4f7b-b38d-f037736c8ab4",
   "metadata": {
    "id": "b5b22299-b22a-4f7b-b38d-f037736c8ab4"
   },
   "outputs": [],
   "source": [
    "# Create an figure with a custom size\n",
    "plt.figure(figsize=(9, 6))\n",
    "\n",
    "# Plot all objects with y == 0 (class 0)\n",
    "plt.scatter(X[y == -1, 0],     # selects all objects with y == 0 and the 1st column of X\n",
    "            X[y == -1, 1],     # selects all objects with y == 0 and the 2nd column of X\n",
    "            color='C3',        # points color\n",
    "            label='-1')        # label for the plot legend\n",
    "\n",
    "\n",
    "# Plot all objects with y == 1 (class 1)\n",
    "plt.scatter(X[y == 1, 0],     # selects all objects with y == 1 and the 1st column of X\n",
    "            X[y == 1, 1],     # selects all objects with y == 1 and the 2nd column of X\n",
    "            color='C0',        # points color\n",
    "            label='1')        # label for the plot legend\n",
    "\n",
    "plt.xlabel('X1') # set up X-axis label\n",
    "plt.ylabel('X2') # set up Y-axis label\n",
    "\n",
    "plt.legend(loc='best') # create the plot legend and set up it position\n",
    "plt.grid() # create grid on the plot\n",
    "\n",
    "plt.show() # display the plot\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "standard-compensation",
   "metadata": {
    "id": "standard-compensation"
   },
   "outputs": [],
   "source": [
    "# Инициализация\n",
    "a = 0\n",
    "\n",
    "# Обучение решающего дерева глубины 1 и получение предсказаний на обучающей выборке\n",
    "dt1 = DecisionTreeRegressor(max_depth = 1)\n",
    "dt1.fit(X, y)\n",
    "dt1_pred = dt1.predict(X)\n",
    "\n",
    "# Включение предсказаний первого дерева в композицию (простым суммированием)\n",
    "a = dt1_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unusual-turtle",
   "metadata": {
    "id": "unusual-turtle"
   },
   "outputs": [],
   "source": [
    "def log_loss(y, y_pred):\n",
    "    return np.log(1 + np.exp(- y * y_pred)).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "interesting-picture",
   "metadata": {
    "id": "interesting-picture"
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "# Вычисление ошибки для одного дерева\n",
    "ll = log_loss(y, a)\n",
    "print(\"LogLoss = %f\" % (ll))\n",
    "\n",
    "# Вычисление ROC AUC\n",
    "auc = roc_auc_score(y, a)\n",
    "print(\"ROC AUC = %f\" % (auc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acting-concentration",
   "metadata": {
    "id": "acting-concentration"
   },
   "outputs": [],
   "source": [
    "# TODO: рассчитайте антиградиент\n",
    "s1 = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "announced-allocation",
   "metadata": {
    "id": "announced-allocation"
   },
   "outputs": [],
   "source": [
    "# Обучение второго решающего дерева глубины 1, предсказывающего антиградиент\n",
    "dt2 = DecisionTreeRegressor(max_depth = 1)\n",
    "dt2.fit(X, s1)\n",
    "\n",
    "# Получение предсказаний на обучающей выборке\n",
    "dt2_pred = dt2.predict(X)\n",
    "\n",
    "# Включение предсказаний второго дерева в композицию (простым суммированием)\n",
    "a += dt2_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "isolated-staff",
   "metadata": {
    "id": "isolated-staff"
   },
   "outputs": [],
   "source": [
    "# Вычисление ошибки для одного дерева\n",
    "ll = log_loss(y, a)\n",
    "print(\"LogLoss = %f\" % (ll))\n",
    "\n",
    "# Вычисление ROC AUC\n",
    "auc = roc_auc_score(y, a)\n",
    "print(\"ROC AUC = %f\" % (auc))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "undefined-cathedral",
   "metadata": {
    "id": "undefined-cathedral"
   },
   "source": [
    "### 3. Особенности градиентного бустинга"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "innocent-bride",
   "metadata": {
    "id": "innocent-bride"
   },
   "source": [
    "#### 3.1 Построение композиции."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "united-commission",
   "metadata": {
    "id": "united-commission"
   },
   "source": [
    "В примерах выше мы строили композицию моделей простым суммированием. Понятно, что в общем виде можно строить взвешенную сумму базовых моделей:\n",
    "\n",
    "$$\n",
    "a(X) = \\sum_{n = 1}^{k} w_nb_n(X),\n",
    "$$\n",
    "\n",
    "где коэффициенты $w_n$ можно подобрать, например, при помощи градиентного спуска в задаче\n",
    "\n",
    "$$\n",
    "w_n = \\arg\\min_{w} \\sum_{i = 1}^{\\ell} L(y_i, a^{(N-1)}(x_i) + w \\times b_N(x_i)),\n",
    "$$\n",
    "\n",
    "что имеет смысл, так как $b_N(x_i)$ уже обучена и фиксирована.\n",
    "\n",
    "**Замечание:** длина шага является одним из ключевых параметров градиентного бустинга, и может очень сильно повлиять на результат.\n",
    "\n",
    "В нашем случае подбор правильного коэффициента поможет решить увеличение MAE в примере 3.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "coordinated-freeze",
   "metadata": {
    "id": "coordinated-freeze"
   },
   "outputs": [],
   "source": [
    "# Инициализация\n",
    "np.random.seed(123)\n",
    "N = 100\n",
    "X = np.linspace(0, 1*np.pi, N).reshape(-1, 1)\n",
    "y = np.sin(X)[:, 0] + np.random.normal(0, 0.1, size = N)\n",
    "a = 0\n",
    "\n",
    "# Обучение решающего дерева глубины 1 и получение предсказаний на обучающей выборке\n",
    "dt1 = DecisionTreeRegressor(max_depth = 1)\n",
    "dt1.fit(X, y)\n",
    "dt1_pred = dt1.predict(X)\n",
    "\n",
    "# Включение предсказаний первого дерева в композицию (простым суммированием)\n",
    "a = dt1_pred\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a, y_pred_label = 'Предсказания a = DT1', loss = 'mae')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "incorporate-stockholm",
   "metadata": {
    "id": "incorporate-stockholm"
   },
   "outputs": [],
   "source": [
    "# Вычисление антиградиента\n",
    "s1 = -np.sign(a - y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "controlled-cameroon",
   "metadata": {
    "id": "controlled-cameroon"
   },
   "outputs": [],
   "source": [
    "# TODO: подберите (угадайте) вес так, чтобы уменьшить MAE\n",
    "w = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "premier-tunisia",
   "metadata": {
    "id": "premier-tunisia"
   },
   "outputs": [],
   "source": [
    "# Обучение второго решающего дерева глубины 1, предсказывающего антиградиент\n",
    "dt2 = DecisionTreeRegressor(max_depth = 1)\n",
    "dt2.fit(X, s1)\n",
    "\n",
    "# Получение предсказаний на обучающей выборке\n",
    "dt2_pred = dt2.predict(X)\n",
    "\n",
    "# Включение предсказаний второго дерева в композицию (простым суммированием)\n",
    "a += w * dt2_pred\n",
    "\n",
    "# Визуализация выборки и предсказаний\n",
    "plot_sample_model(X, y, plot_predictions = True, y_pred = a,\n",
    "                  y_pred_label = 'Предсказания a = DT1 + DT2', loss = 'mae')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "personal-winner",
   "metadata": {
    "id": "personal-winner"
   },
   "source": [
    "#### 3.2 Переобучение."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incorrect-clock",
   "metadata": {
    "id": "incorrect-clock"
   },
   "source": [
    "Вспомним, что в случае бэггинга увеличение числа базовых моделей не приводит к переобучению. Мы могли бы наблюдать, что вслучайном лесе ошибки на тренировочной и тестовой выборках выходили на плато, начиная с какого-то числа деревьев. Понятно, что бустинг не будет обладать этим свойством *по построению*: при добавлении новых базовых моделей композиция всё точнее будет приближать обучающую выборку, что в конечном итоге (зачастую быстро) может привести к переобучению. Из этого можно сделать вывод, что в случае бустинга ошибка на обучающей выборке в зависимости от числа деревьев является убывающей функцией, а ошибка на тестовой выборке, скорее всего, имеет U-образный вид. Убедимся в этом и сравним поведение ошибок и решающих поверхностей бустинга и случайного леса.\n",
    "\n",
    "Здесь и далее будем использовать [реализацию бустинга](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html) из `sklearn`, которая представляет собой градиентный бустинг над деревьями."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "special-entity",
   "metadata": {
    "id": "special-entity"
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tribal-grammar",
   "metadata": {
    "id": "tribal-grammar"
   },
   "outputs": [],
   "source": [
    "# Генерация выборки\n",
    "np.random.seed(123)\n",
    "X = np.linspace(0, 1, 300).reshape(-1, 1)\n",
    "\n",
    "def target(a):\n",
    "    return a > 0.5\n",
    "\n",
    "y = target(X) + np.random.normal(size = X.shape) * 0.1\n",
    "y = y[:, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cordless-israeli",
   "metadata": {
    "id": "cordless-israeli"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "linear-necklace",
   "metadata": {
    "id": "linear-necklace"
   },
   "outputs": [],
   "source": [
    "# Визуализация выборки\n",
    "plt.figure(figsize = (10, 7))\n",
    "plt.scatter(X_train, y_train, label = 'Обучающая выборка')\n",
    "plt.scatter(X_test, y_test, label = 'Тестовая выборка')\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "duplicate-patient",
   "metadata": {
    "id": "duplicate-patient"
   },
   "outputs": [],
   "source": [
    "# Визуализация решающих поверхностей\n",
    "trees = [1, 2, 5, 20, 100, 500, 1000]\n",
    "\n",
    "fig, ax = plt.subplots(len(trees), 2, figsize = (30, 40))\n",
    "\n",
    "loss_rf_train = []\n",
    "loss_gb_train = []\n",
    "loss_rf_test = []\n",
    "loss_gb_test = []\n",
    "\n",
    "for i, ts in enumerate(trees):\n",
    "    rf = RandomForestRegressor(n_estimators = ts, max_depth = 1)\n",
    "    gb = GradientBoostingRegressor(n_estimators = ts, max_depth = 1, learning_rate = 1)\n",
    "\n",
    "    rf.fit(X_train, y_train)\n",
    "    ax[i, 0].scatter(X_train, y_train, label = 'Обучающая выборка')\n",
    "    ax[i, 0].scatter(X_test, y_test, label = 'Тестовая выборка')\n",
    "    ax[i, 0].plot(sorted(X_test), rf.predict(sorted(X_test)), lw = 3, c = 'red', label = 'Предсказания на тестовой выборке')\n",
    "    ax[i, 0].set_xlabel('X')\n",
    "    ax[i, 0].set_ylabel('Y')\n",
    "    ax[i, 0].set_title('Случайный лес, число деревьев = ' + str(ts) + ', MSE = ' + str(mean_squared_error(y_test, rf.predict(X_test))))\n",
    "    ax[i, 0].legend();\n",
    "\n",
    "    loss_rf_train.append(mean_squared_error(y_train, rf.predict(X_train)))\n",
    "    loss_rf_test.append(mean_squared_error(y_test, rf.predict(X_test)))\n",
    "\n",
    "    gb.fit(X_train, y_train)\n",
    "    ax[i, 1].scatter(X_train, y_train, label = 'Обучающая выборка')\n",
    "    ax[i, 1].scatter(X_test, y_test, label = 'Тестовая выборка')\n",
    "    ax[i, 1].plot(sorted(X_test), gb.predict(sorted(X_test)), lw = 3, c = 'red', label = 'Предсказания на тестовой выборке')\n",
    "    ax[i, 1].set_xlabel('X')\n",
    "    ax[i, 1].set_ylabel('Y')\n",
    "    ax[i, 1].set_title('Градиентный бустинг, число деревьев = ' + str(ts) + ', MSE = ' + str(mean_squared_error(y_test, gb.predict(X_test))))\n",
    "    ax[i, 1].legend();\n",
    "\n",
    "    loss_gb_train.append(mean_squared_error(y_train, gb.predict(X_train)))\n",
    "    loss_gb_test.append(mean_squared_error(y_test, gb.predict(X_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "isolated-czech",
   "metadata": {
    "id": "isolated-czech"
   },
   "source": [
    "Как мы видим, решающая поверхность случайного леса не изменяется (уже со второго дерева), в то время как для бустинга она становится всё более сложной."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cleared-accused",
   "metadata": {
    "id": "cleared-accused"
   },
   "outputs": [],
   "source": [
    "# Визуализация ошибок\n",
    "fig, ax = plt.subplots(1, 2, figsize = (20, 7))\n",
    "\n",
    "ax[0].plot(trees, loss_rf_train, label = 'MSE_Train')\n",
    "ax[0].plot(trees, loss_rf_test, label = 'MSE_Test')\n",
    "ax[0].set_xlabel('Число деревьев')\n",
    "ax[0].set_ylabel('MSE')\n",
    "ax[0].set_title('Случайный лес')\n",
    "ax[0].legend()\n",
    "\n",
    "ax[1].plot(trees, loss_gb_train, label = 'MSE_Train')\n",
    "ax[1].plot(trees, loss_gb_test, label = 'MSE_Test')\n",
    "ax[1].set_xlabel('Число деревьев')\n",
    "ax[1].set_ylabel('MSE')\n",
    "ax[1].set_title('Градиентный бустинг')\n",
    "ax[1].legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "handy-standard",
   "metadata": {
    "id": "handy-standard"
   },
   "source": [
    "В случае бустинга ошибка на тренировочной выборке стабильно убывает, а на тестовой – возрастает, что вероятно говорит об (очень быстром) переобучении. Если же мы аккуратно подберём гиперпараметры для деревьев и длину шага, то сможем добиться типичной U-образной формы функции потерь на тестовой выборке.\n",
    "\n",
    "В качестве способа регуляризации также может использоваться сокращение шага:\n",
    "\n",
    "$$\n",
    "a^{(N)}(X) = a^{(N-1)}(X) + \\alpha w_N b^{(N)}(X),\n",
    "$$\n",
    "\n",
    "где $\\alpha \\in (0, 1]$ – темп обучения, или стохастический градиентный бустинг.\n",
    "\n",
    "**Важный момент:** примеры выше демонстрируют, насколько быстро бустинг может переобучаться, и как в его случае важен подбор гиперпараметров. В связи с этим может сложиться ситуация, что градиентный бустинг с дефолтными параметрами показывает более плохое качество, чем, например, случайный лес с дефолтными параметрами, однако хорошо настроенный градиентный бустинг [обычно превосходит](https://www.quora.com/How-can-the-performance-of-a-Gradient-Boosting-Machine-be-worse-than-Random-Forests) случайный лес по качеству. Важна и конкретная реализация бустинга: например, [здесь](https://towardsdatascience.com/boosting-showdown-scikit-learn-vs-xgboost-vs-lightgbm-vs-catboost-in-sentiment-classification-f7c7f46fd956), представлено очень хорошее сравнение различных имплементаций."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "proper-theory",
   "metadata": {
    "id": "proper-theory"
   },
   "source": [
    "#### 3.3 Пара слов о смещении и дисперсии."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "referenced-insulin",
   "metadata": {
    "id": "referenced-insulin"
   },
   "source": [
    "Вспомним, что в случайном лесе в качестве базовых моделей мы брали переобученные деревья, мотивируя это тем, что бэггинг позволит сохранить низкое смещение и при этом уменьшить дисперсию. Бустинг по построению работает проивоположно: композиция обладает более низким, чем базовые модели, смещением, но такой же или большей дисперсией (ещё раз проговорите последовательность построения бустинга, и вы увидите, почему это так). Поэтому в качестве базовых моделей для бустинга обычно используются модели с высоким смещением и низкой дисперсией – как мы знаем, такие модели являются недообученными – например, неглубокие решающие деревья."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "nCCSHAlVJs0M",
   "metadata": {
    "id": "nCCSHAlVJs0M"
   },
   "source": [
    "### 4. Catboost vs. XGBoost vs. LightGBM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbvyq4YYKMPK",
   "metadata": {
    "id": "cbvyq4YYKMPK"
   },
   "source": [
    "[Source](https://github.com/Murcha1990/ML_econom_2022-2023/blob/main/%D0%A1%D0%B5%D0%BC%D0%B8%D0%BD%D0%B0%D1%80%D1%8B/%D0%A1%D0%B5%D0%BC%D0%B8%D0%BD%D0%B0%D1%80%2010*/Seminar10_%D0%B1%D0%BE%D0%BD%D1%83%D1%81.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "N2igDAelKEDT",
   "metadata": {
    "id": "N2igDAelKEDT"
   },
   "source": [
    "Ранее мы смотрели простую версию градиентного бустинга из scikit-learn, [придуманную в 1999 году Фридманом](https://projecteuclid.org/download/pdf_1/euclid.aos/1013203451). Прогресс не стоит на месте, и на сегодняшний день есть три популярные библиотеки с разными имплементациями градиентного бустинга, которые на практике показывают лушие результаты:\n",
    "*  **XGBoost**. Появилась в 2014 году, [статья автора](https://arxiv.org/pdf/1603.02754.pdf) вышла в 2016. После выхода быстро набрала популярность и оставалась стандартом до конца 2016 года.\n",
    "* **CatBoost** от компании Яндекс с релизом в 2017 году. Алгоритм можно запускать с дефолтными гиперпараметрами, потому тчо он является менее чувствительным к выбору их конкретных значений. Отлично умеет работать с категориальным признаками, при этом автоматически обрабатывая полученные на вход непредобработанные фичи.\n",
    "* **LightGBM**. Релиз в один год с Catboost, библиотека от Microsoft. Отличается очень быстрым построением композиции. Например, при построении узла дерева, вместо перебора по всем значениям признака, производится перебор значений гистограммы этого признака. Таким образом, вместо $O(N)$ требуется $O$(m), где $m$ - число бинов гистограммы. В отличие от других библиотек, строит деревья в глубину, при этом на каждом шаге строит вершину, дающую наибольшее уменьшение функционала."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "TnNK1A5wMIVX",
   "metadata": {
    "id": "TnNK1A5wMIVX"
   },
   "source": [
    "**Техническое отступление**\n",
    "\n",
    "Данные библиотеки необходимо сначала устанавливать (можно через pip / conda или brew, если Вы работаете на MAC OS).\n",
    "Чтобы у Вас точно вопроизводился ноутбук и не было проблем из-за несовпадающих версий библиотек, рекомендуется через python создавать виртуальную среду."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "qo3pShv-MhOD",
   "metadata": {
    "id": "qo3pShv-MhOD"
   },
   "outputs": [],
   "source": [
    "#!pip install catboost lightgbm cmake xgboost # без cmake xgboost установится, но не будет импортироваться\n",
    "#!pip install catboost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tx_iiQ7wNGOo",
   "metadata": {
    "id": "tx_iiQ7wNGOo"
   },
   "outputs": [],
   "source": [
    "import catboost\n",
    "import lightgbm\n",
    "import xgboost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "an-hatr-OP9q",
   "metadata": {
    "id": "an-hatr-OP9q"
   },
   "outputs": [],
   "source": [
    "#import warnings\n",
    "#warnings.filterwarnings('ignore')\n",
    "\n",
    "from sklearn.datasets import make_classification\n",
    "from matplotlib.colors import ListedColormap\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5LpSQuaNJ2aM",
   "metadata": {
    "id": "5LpSQuaNJ2aM"
   },
   "source": [
    "#### 4.1 Catboost"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "zhhLcTaYNezJ",
   "metadata": {
    "id": "zhhLcTaYNezJ"
   },
   "source": [
    "В алгоритме сделаны улучшения и выбор разных опций для борьбы с переобучением, подсчету среднего таргета на отложенной выборке, подсчету статистик по категориальным фичам, бинаризацией фичей, рандомизации скора сплита, разные типы бутсрапирования.\n",
    "\n",
    "Давайте сначала зафиксируем все гиперпараметры со значениями по умолчанию, кроме количества деревьев в композиции - `n_estimators`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "certified-scottish",
   "metadata": {
    "id": "certified-scottish"
   },
   "outputs": [],
   "source": [
    "def plot_surface(X, y, clf):\n",
    "    h = 0.2\n",
    "    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1\n",
    "    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1\n",
    "    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),\n",
    "                         np.arange(y_min, y_max, h))\n",
    "    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "\n",
    "    Z = Z.reshape(xx.shape)\n",
    "    cmap_light = ListedColormap(['#FFAAAA', '#AAFFAA', '#AAAAFF'])\n",
    "    plt.pcolormesh(xx, yy, Z, cmap=cmap_light)\n",
    "\n",
    "    # Добавим на график сами наблюдения\n",
    "    cmap_bold = ListedColormap(['#FF0000', '#00FF00', '#0000FF'])\n",
    "    plt.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold)\n",
    "    plt.xlim(xx.min(), xx.max())\n",
    "    plt.ylim(yy.min(), yy.max())\n",
    "\n",
    "X, y = make_classification(n_samples=500, n_features=2, n_informative=2,\n",
    "                           n_redundant=0, n_repeated=0,\n",
    "                           n_classes=2, n_clusters_per_class=2,\n",
    "                           flip_y=0.05, class_sep=0.8, random_state=241)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=241)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1VvnWMoIOnSJ",
   "metadata": {
    "id": "1VvnWMoIOnSJ"
   },
   "outputs": [],
   "source": [
    "from catboost import CatBoostClassifier\n",
    "\n",
    "catboost = CatBoostClassifier(n_estimators=300, logging_level='Silent')\n",
    "catboost.fit(X_train, y_train)\n",
    "plot_surface(X_test, y_test, catboost)\n",
    "\n",
    "print(roc_auc_score(y_test, catboost.predict_proba(X_test)[:, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "XSZ9dBadPFDG",
   "metadata": {
    "id": "XSZ9dBadPFDG"
   },
   "source": [
    "#### 4.2 XGBoost"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "__1bbzmcPs7-",
   "metadata": {
    "id": "__1bbzmcPs7-"
   },
   "source": [
    "1. Базовый алгоритм приближает направление, посчитанное с учетом второй производной функции потерь\n",
    "\n",
    "2. Функционал регуляризуется – добавляются штрафы за количество листьев и за норму коэффициентов\n",
    "\n",
    "3. При построении дерева используется критерий информативности, зависящий от оптимального вектора сдвига\n",
    "\n",
    "4. Критерий останова при обучении дерева также зависит от оптимального сдвига\n",
    "\n",
    "Ссылка на [источник](https://github.com/esokolov/ml-course-hse/blob/master/2021-fall/lecture-notes/lecture11-ensembles.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "UkFeG5bqO1t0",
   "metadata": {
    "id": "UkFeG5bqO1t0"
   },
   "outputs": [],
   "source": [
    "from xgboost import XGBClassifier\n",
    "\n",
    "xgb = XGBClassifier(n_estimators=300, verbosity=0)\n",
    "xgb.fit(X_train, y_train)\n",
    "plot_surface(X_test, y_test, xgb)\n",
    "\n",
    "print(roc_auc_score(y_test, xgb.predict_proba(X_test)[:, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "TgNkZeI3Qe7s",
   "metadata": {
    "id": "TgNkZeI3Qe7s"
   },
   "source": [
    "#### 4.3 LightGBM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "LvygU_3uP7es",
   "metadata": {
    "id": "LvygU_3uP7es"
   },
   "outputs": [],
   "source": [
    "from lightgbm import LGBMClassifier\n",
    "\n",
    "lightgbm = LGBMClassifier(n_estimators=300)\n",
    "lightgbm.fit(X_train, y_train)\n",
    "plot_surface(X_test, y_test, lightgbm)\n",
    "\n",
    "print(roc_auc_score(y_test, lightgbm.predict_proba(X_test)[:, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "FBaHmQivRupR",
   "metadata": {
    "id": "FBaHmQivRupR"
   },
   "source": [
    "Когда мы обучили лучшие версии моделей, можно их сохранить и использовать для получения предсказаний, например, на новом батче данных."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "v-rD2hijRWGH",
   "metadata": {
    "id": "v-rD2hijRWGH"
   },
   "outputs": [],
   "source": [
    "'''\n",
    "# Сохранить\n",
    "lightgbm.booster_.save_model('lightgbm.txt')\n",
    "catboost.save_model('catboost.cbm', format='cbm')\n",
    "xgboost.save_model('xgboost.json')\n",
    "\n",
    "# Загрузить\n",
    "lightgbm = LGBMClassifier(model_file='mode.txt')\n",
    "catboost = catboost.load_model('catboost.cbm')\n",
    "xgboost = xgboost.load_model('xgboost.json')\n",
    "''';"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "_CNQDrPbSC0S",
   "metadata": {
    "id": "_CNQDrPbSC0S"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}