{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "from sklearn.exceptions import ConvergenceWarning\n",
    "import statsmodels.api as sm\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_palette(\"colorblind\")\n",
    "palette = sns.color_palette()\n",
    "figsize = (11,6)\n",
    "legend_fontsize = 16\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('font',**{'family':'sans-serif'})\n",
    "rc('text', usetex=True)\n",
    "rc('text.latex',preamble=r'\\usepackage[utf8]{inputenc}')\n",
    "rc('text.latex',preamble=r'\\usepackage[russian]{babel}')\n",
    "rc('figure', **{'dpi': 300})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Байесовский вывод для испытаний Бернулли"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = np.arange(-0.5, 1.5, 0.0025)\n",
    "\n",
    "prior_params = (0.25, 0.25)\n",
    "experimental_data = (1, 0)\n",
    "\n",
    "## Априорное распределение\n",
    "pri = st.beta(prior_params[0], prior_params[1]).pdf\n",
    "norm_prior = 1 #integrate.quad(lambda x: pri(x), -np.inf, np.inf)[0]\n",
    "ys_prior = [ pri(x) / norm_prior for x in xs ]\n",
    "\n",
    "## Правдоподобие\n",
    "n_heads, n_tails = experimental_data\n",
    "lk = lambda x : x ** n_heads * (1 - x) ** n_tails\n",
    "ys_like = [ lk(x) for x in xs ]\n",
    "\n",
    "## Апостериорное распределение\n",
    "post = lambda x : lk(x) * pri(x)\n",
    "norm_post = integrate.quad(lambda x: post(x), 0, 1)[0]\n",
    "ys_post = [ post(x) / norm_post if x > 0 and x < 1 else 0 for x in xs ]\n",
    "\n",
    "## И нарисуем\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(xs, ys_prior, linewidth=2, label=r\"Априорное распределение\")\n",
    "ax.plot(xs, ys_like, linewidth=2, label=r\"Правдоподобие\")\n",
    "ax.plot(xs, ys_post, linewidth=2, label=r\"Апостериорное распределение\")\n",
    "ax.set_xlim((-0.05, 1.05))\n",
    "ax.set_ylim((-0.5, 8.5))\n",
    "ax.set_xlabel(r\"Вероятность орла $\\theta$\", fontsize=legend_fontsize)\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Линейная и полиномиальная регрессия"
   ]
  },
  {
   "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",
    "## X-координаты точек данных\n",
    "xd_large = np.arange(-1.5, 2, 0.05)\n",
    "num_points_l = len(xd_large)\n",
    "\n",
    "## Данные\n",
    "data_large = orig(xd_large) + np.random.normal(0, .25, num_points_l)\n",
    "\n",
    "\n",
    "## Для рисования\n",
    "xs = np.arange(xd[0]-1.5, xd[-1]+1.5, 0.01)"
   ]
  },
  {
   "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 = [ 1, 3]\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 = linear_model.LinearRegression(fit_intercept=True).fit( 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 left\", fontsize=legend_fontsize)\n",
    "ax.set_xlim((-2, 2.5))\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Локальные признаки в линейной регрессии"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.base import BaseEstimator, TransformerMixin\n",
    "\n",
    "nums_gauss = [2 ]\n",
    "gauss_xd, gauss_yd = xd, data\n",
    "\n",
    "class GaussianFeatures(BaseEstimator, TransformerMixin):\n",
    "    def __init__(self, N, width_factor=.5):\n",
    "        self.N = N\n",
    "        self.width_factor = width_factor\n",
    "    \n",
    "    @staticmethod\n",
    "    def _gauss_basis(x, y, width, axis=None):\n",
    "        arg = (x - y) / width\n",
    "        return np.exp(-0.5 * np.sum(arg ** 2, axis))\n",
    "        \n",
    "    def fit(self, X, y=None):\n",
    "        # create N centers spread along the data range\n",
    "        self.centers_ = np.linspace(X.min(), X.max(), self.N)\n",
    "        self.width_ = self.width_factor * (self.centers_[1] - self.centers_[0])\n",
    "        return self\n",
    "        \n",
    "    def transform(self, X):\n",
    "        return self._gauss_basis(X[:, :, np.newaxis], self.centers_, self.width_, axis=1)\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((-1.5, 1.5))\n",
    "ax.scatter(gauss_xd, gauss_yd, marker='*', s=120)\n",
    "ax.plot(xs, orig(xs), linewidth=1, label=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "for num_gauss in nums_gauss:\n",
    "    gauss_model = make_pipeline(GaussianFeatures(num_gauss),\n",
    "                                linear_model.LinearRegression())\n",
    "    gauss_model.fit(xd[:, np.newaxis], data)\n",
    "    yfit = gauss_model.predict(xs[:, np.newaxis])\n",
    "    print(\"Коэффициенты с %d признаками: %s\" % (num_gauss, \" \".join([\"%.4f\" % x for x in gauss_model.get_params()['linearregression'].coef_])))\n",
    "    ax.plot(xs, yfit, linewidth=2, label=\"%d гауссовских признак%s\" % (num_gauss, \"а\" if num_gauss < 5 else \"ов\") )\n",
    "\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_gauss = 5\n",
    "gauss_model = make_pipeline(GaussianFeatures(num_gauss), linear_model.LinearRegression())\n",
    "gauss_model.fit(gauss_xd[:, np.newaxis], gauss_yd)\n",
    "yfit = gauss_model.predict(xs[:, np.newaxis])\n",
    "mfeat = gauss_model.get_params()['gaussianfeatures']\n",
    "mregr = gauss_model.get_params()['linearregression']\n",
    "print(\"Коэффициенты с %d признаками: %s\" % (num_gauss, \" \".join([\"%.4f\" % x for x in gauss_model.get_params()['linearregression'].coef_])))\n",
    "print(\"Свободный член: %.4f\" % mregr.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 7))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlim((xs[0], xs[-1]))\n",
    "# ax.set_ylim((-1.5, 1.5))\n",
    "ax.scatter(gauss_xd, gauss_yd, marker='*', s=120)\n",
    "ax.plot(xs, orig(xs), linewidth=1, label=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "for i in range(mfeat.N):\n",
    "    cur_yfit = [mregr.coef_[i] * mfeat._gauss_basis(x, mfeat.centers_[i], mfeat.width_) for x in xs]\n",
    "    if i == 0:\n",
    "        ax.plot(xs, cur_yfit, color=\"0.4\", linewidth=1, label=\"Один гауссовский признак\")\n",
    "    else:\n",
    "        ax.plot(xs, cur_yfit, color=\"0.4\", linewidth=1)\n",
    "    ax.plot(xs, [mregr.intercept_ for _ in range(len(xs))], color=\"0.6\", linewidth=1)\n",
    "\n",
    "ax.plot(xs, yfit, linewidth=2, label=\"Регрессия с гауссовскими приизнаками\")\n",
    "ax.legend(loc=\"upper center\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlim((xs[0], xs[-1]))\n",
    "\n",
    "for i in range(mfeat.N):\n",
    "    cur_yfit = [mfeat._gauss_basis(x, mfeat.centers_[i], mfeat.width_) for x in xs]\n",
    "    ax.plot(xs, cur_yfit, color=\"0.6\", linewidth=1)\n",
    "    ax.plot(xs, [mregr.intercept_ for _ in range(len(xs))], color=\"0.6\", linewidth=1)\n",
    "\n",
    "ax.plot(xs, yfit, linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Добавим ещё данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Исходная функция\n",
    "orig = lambda x : np.sin(2*x)\n",
    "\n",
    "## Для \n",
    "xs = np.arange(xd_large[0]-.5, xd_large[-1]+.5, 0.01)\n",
    "\n",
    "## Выделение полиномиальных признаков\n",
    "xs_d = np.vstack([xs ** i for i in range(1, num_points+1)]).transpose()\n",
    "xd_d_large = np.vstack([xd_large ** i for i in range(1, num_points+1)]).transpose()\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_large, data_large, marker='*', s=120)\n",
    "ax.plot(xs, orig(xs), linewidth=2, label=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "set_of_powers = [ 1, 3, 10 ]\n",
    "\n",
    "for d in set_of_powers:\n",
    "    cur_model = linear_model.LinearRegression(fit_intercept=True).fit( xd_d_large[:, :d], data_large )\n",
    "    ax.plot(xs, cur_model.predict( xs_d[:, :d] ), linewidth=2, label=\"$d=%d$\" % d)\n",
    "\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_gauss=20\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((-1.5, 1.5))\n",
    "ax.scatter(xd_large, data_large, marker='*', s=120)\n",
    "# ax.plot(xs, orig(xs), linewidth=1, label=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "gauss_model = make_pipeline(GaussianFeatures(num_gauss), linear_model.LinearRegression())\n",
    "gauss_model.fit(xd_large[:, np.newaxis], data_large)\n",
    "yfit = gauss_model.predict(xs[:, np.newaxis])\n",
    "mfeat = gauss_model.get_params()['gaussianfeatures']\n",
    "mregr = gauss_model.get_params()['linearregression']\n",
    "\n",
    "for i in range(mfeat.N):\n",
    "    cur_yfit = [mregr.coef_[i] * mfeat._gauss_basis(x, mfeat.centers_[i], mfeat.width_) for x in xs]\n",
    "    ax.plot(xs, cur_yfit, color=\"0.4\", linewidth=1)\n",
    "    ax.plot(xs, [mregr.intercept_ for _ in range(len(xs))], color=\"0.6\", linewidth=1)\n",
    "\n",
    "ax.plot(xs, yfit, color=\"C1\", linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Регуляризация"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_values = [ 0, .01, 1]\n",
    "use_lasso = True\n",
    "\n",
    "\n",
    "def train_model(xs, ys, alpha, use_lasso):\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 )\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((-3, 3))\n",
    "ax.scatter(xd, data, marker='*', s=120)\n",
    "ax.plot(xs, orig(xs), linewidth=2, label=\"Original function\", color=\"black\")\n",
    "\n",
    "for alpha in alpha_values:\n",
    "    m = train_model(xd_d[:, :10], data, alpha, use_lasso)\n",
    "    print (\"Коэффициенты для alpha=%f:\\n%s\" % (alpha, m.coef_))\n",
    "    ax.plot(xs, m.predict( xs_d[:, :10] ), linewidth=2, label=\"$\\\\alpha=%f$\" % alpha)\n",
    "\n",
    "ax.legend(loc=\"upper center\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}