{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Смещение, разброс и шум"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import json\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats as st\n",
    "import scipy.integrate as integrate\n",
    "from sklearn import linear_model\n",
    "import warnings\n",
    "from sklearn.exceptions import ConvergenceWarning\n",
    "warnings.filterwarnings('ignore', category=ConvergenceWarning)\n",
    "\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_palette(\"colorblind\")\n",
    "palette = sns.color_palette()\n",
    "figsize = (15,8)\n",
    "legend_fontsize = 16\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('font',**{'family':'sans-serif'})\n",
    "rc('text', usetex=False)\n",
    "np.random.seed(30)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Оверфиттинг\n",
    "## Исходная функция\n",
    "orig = lambda x : np.sin(2*x)\n",
    "\n",
    "## X-координаты точек данных\n",
    "xd = np.array([-3, -2, -1, -0.5, 0, 0.5, 1, 1.5, 2.5, 3, 4]) / 2\n",
    "num_points = len(xd)\n",
    "\n",
    "## Данные\n",
    "data = orig(xd) + np.random.normal(0, .25, num_points)\n",
    "\n",
    "## Для рисования\n",
    "xs = np.arange(xd[0]-1.5, xd[-1]+1.5, 0.01)\n",
    "\n",
    "## Обучаем модель с регуляризацией\n",
    "def train_model(xs, ys, alpha=0, use_lasso=False):\n",
    "    if alpha == 0:\n",
    "        return linear_model.LinearRegression(fit_intercept=True).fit( xs, ys )\n",
    "    else:\n",
    "        if use_lasso:\n",
    "            return linear_model.Lasso(alpha=alpha, fit_intercept=True).fit( xs, ys )\n",
    "        else:\n",
    "            return linear_model.Ridge(alpha=alpha, fit_intercept=True).fit( xs, ys )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Выделение полиномиальных признаков\n",
    "xs_d = np.vstack([xs ** i for i in range(1, num_points+1)]).transpose()\n",
    "xd_d = np.vstack([xd ** i for i in range(1, num_points+1)]).transpose()\n",
    "\n",
    "## Какие степени многочлена будем обучать и рисовать\n",
    "set_of_powers = [ 3, 10 ]\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlim((xs[0], xs[-1]))\n",
    "ax.set_ylim((-2, 2))\n",
    "ax.scatter(xd, data, marker='*', s=120)\n",
    "ax.plot(xs, orig(xs), linewidth=1, label=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "for d in set_of_powers:\n",
    "    if d == 0:\n",
    "        print(np.mean(data))\n",
    "        ax.hlines(np.mean(data), xmin=xs[0], xmax=xs[-1], label=\"$d=0$\", linestyle=\"dashed\")\n",
    "    else:\n",
    "        cur_model = train_model( xd_d[:, :d], data )\n",
    "        print(cur_model.coef_)\n",
    "        ax.plot(xs, cur_model.predict( xs_d[:, :d] ), linewidth=2, label=\"$d=%d$\" % d)\n",
    "\n",
    "ax.legend(loc=\"upper center\", fontsize=legend_fontsize)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 500\n",
    "alpha = 10.1\n",
    "use_lasso = False\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlim((xs[0], xs[-1]))\n",
    "ax.set_ylim((-5, 5))\n",
    "\n",
    "res = []\n",
    "d = 10\n",
    "for _ in range(N):\n",
    "    cur_data = orig(xd) + np.random.normal(0, .25, num_points)\n",
    "    cur_model = train_model(xd_d[:, :d], cur_data, alpha, use_lasso)\n",
    "    res.append(cur_model.predict( xs_d[:, :d] ))\n",
    "    ax.plot(xs, res[-1], linewidth=.15, color=\"0.35\")\n",
    "\n",
    "ax.plot(xs, orig(xs), linewidth=2, label=\"Исходная функция\", color=palette[0])\n",
    "ax.scatter(xd, orig(xd), marker='*', s=150, color=palette[0])\n",
    "\n",
    "ax.plot(xs, np.mean( res, axis=0 ), linewidth=2, label=\"Усреднённые предсказания\", color=\"red\")\n",
    "ax.legend(loc=\"upper center\", fontsize=legend_fontsize)\n",
    "plt.show()\n",
    "\n",
    "cur_model = linear_model.LinearRegression(fit_intercept=True).fit( xd_d[:, :d], data )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_set_size = 30\n",
    "test_set_x = np.random.rand(test_set_size) * (2 + 1.5) - 1.5\n",
    "test_set_xs = np.vstack([test_set_x ** i for i in range(1, 12)]).transpose()\n",
    "print(test_set_xs.shape)\n",
    "test_set_y = orig(test_set_xs[:,0]) + np.random.normal(0, .25, test_set_size)\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(xs, orig(xs), linewidth=2, label=\"Исходная функция\", color=palette[0])\n",
    "ax.scatter(test_set_xs[:,0], test_set_y, marker='*', s=150, color=palette[0])\n",
    "ax.set_xlim((-1.5, 2))\n",
    "\n",
    "def test_set_error(model, d):\n",
    "    return np.mean( (test_set_y - model.predict(test_set_xs[:, :d])) ** 2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Усредняем по новым обучающим выборкам; шум нового ответа независим.\n",
    "N = 1000\n",
    "alphas = np.logspace(-4, 2, num=20)\n",
    "noise_variance = 0.25 ** 2\n",
    "errors, biases, variances = [], [], []\n",
    "truth = orig(test_set_x)\n",
    "for alpha in alphas:\n",
    "    predictions = []\n",
    "    noisy_errors = []\n",
    "    for _ in range(N):\n",
    "        cur_data = orig(xd) + np.random.normal(0, .25, num_points)\n",
    "        cur_model = train_model(xd_d, cur_data, alpha, use_lasso=False)\n",
    "        pred = cur_model.predict(test_set_xs)\n",
    "        predictions.append(pred)\n",
    "        fresh_y = truth + np.random.normal(0, .25, test_set_size)\n",
    "        noisy_errors.append(np.mean((fresh_y-pred)**2))\n",
    "    predictions = np.asarray(predictions)\n",
    "    mean_prediction = predictions.mean(axis=0)\n",
    "    biases.append(np.mean((mean_prediction-truth)**2))\n",
    "    variances.append(np.mean((predictions-mean_prediction)**2))\n",
    "    errors.append(np.mean(noisy_errors))\n",
    "    print(f\"alpha={alpha:.4g}: bias²={biases[-1]:.4f}, \"\n",
    "          f\"variance={variances[-1]:.4f}, MSE={errors[-1]:.4f}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=figsize)\n",
    "ax.set_xscale(\"log\")\n",
    "ax.plot(alphas, biases, label=\"Смещение²\", linewidth=2)\n",
    "ax.plot(alphas, variances, label=\"Разброс\", linewidth=2)\n",
    "ax.axhline(noise_variance, label=\"Шум\", linestyle=\":\")\n",
    "decomposition = np.array(biases) + np.array(variances) + noise_variance\n",
    "ax.plot(alphas, decomposition, label=\"Смещение² + разброс + шум\", linewidth=2)\n",
    "ax.plot(alphas, errors, \"--\", color=\"black\", label=\"MSE на новых ответах\")\n",
    "ax.set(xlabel=\"Сила регуляризации α\", ylabel=\"Средняя квадратичная ошибка\", ylim=(0,1))\n",
    "ax.legend(fontsize=14)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.13.5"
  },
  "lecture_source": "soft/model-selection/bias-variance-noise--samsung-05.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}