{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-000",
   "metadata": {},
   "source": [
    "# Байесовский вывод: испытания Бернулли и линейная регрессия"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-001",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats as st\n",
    "import scipy.integrate as integrate\n",
    "\n",
    "from sklearn import linear_model\n",
    "from sklearn.base import BaseEstimator, TransformerMixin\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "## Ложные предупреждения matmul от Apple Accelerate (NumPy >= 1.25)\n",
    "warnings.filterwarnings(\"ignore\", message=r\".*encountered in matmul\",\n",
    "                        category=RuntimeWarning)\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}\\usepackage[russian]{babel}')\n",
    "rc('figure', **{'dpi': 200})\n",
    "\n",
    "SEED = 2026\n",
    "np.random.seed(SEED)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-002",
   "metadata": {},
   "source": [
    "# 1. Байесовский вывод для испытаний Бернулли"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_bernoulli_inference(prior_params, experimental_data, ax=None,\n",
    "                             xlim=(-0.05, 1.05), ylim=None, title=None):\n",
    "    a, b = prior_params\n",
    "    n_heads, n_tails = experimental_data\n",
    "    xs = np.arange(-0.5, 1.5, 0.0025)\n",
    "\n",
    "    ## Только внутри (0, 1): у Beta(0.5, 0.5) на концах inf, а правдоподобие 0\n",
    "    inside = (xs > 0) & (xs < 1)\n",
    "    x_in = xs[inside]\n",
    "\n",
    "    ## Априорное распределение\n",
    "    pri = st.beta(a, b).pdf\n",
    "    ys_prior = pri(xs)\n",
    "    ys_prior[~np.isfinite(ys_prior)] = np.nan\n",
    "\n",
    "    ## Правдоподобие\n",
    "    lk = lambda x: x ** n_heads * (1 - x) ** n_tails\n",
    "    ys_like = np.zeros_like(xs)\n",
    "    ys_like[inside] = lk(x_in)\n",
    "\n",
    "    ## Апостериорное распределение\n",
    "    post = lambda x: lk(x) * pri(x)\n",
    "    norm_post = integrate.quad(post, 0, 1)[0]\n",
    "    ys_post = np.zeros_like(xs)\n",
    "    ys_post[inside] = post(x_in) / norm_post\n",
    "\n",
    "    if ax is None:\n",
    "        fig = plt.figure(figsize=figsize)\n",
    "        ax = fig.add_subplot(111)\n",
    "\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",
    "\n",
    "    ax.set_xlim(xlim)\n",
    "    if ylim is not None:\n",
    "        ax.set_ylim(ylim)\n",
    "    ax.set_xlabel(r\"Вероятность орла $\\theta$\", fontsize=legend_fontsize)\n",
    "    if title is not None:\n",
    "        ax.set_title(title, fontsize=legend_fontsize)\n",
    "    ax.legend(loc=\"upper left\", fontsize=legend_fontsize - 2)\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-004",
   "metadata": {},
   "source": [
    "### Сильное априорное распределение против небольшой выборки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-005",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bernoulli_inference((40, 40), (15, 0), ylim=(-0.5, 8.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-006",
   "metadata": {},
   "source": [
    "### Когда данных много, априорное распределение почти не важно"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-007",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_bernoulli_inference((10, 10), (10, 100), ylim=(-0.5, 12))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-008",
   "metadata": {},
   "source": [
    "### Несколько априорных распределений рядом"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-009",
   "metadata": {},
   "outputs": [],
   "source": [
    "cases = [((5, 15), (3, 0), r\"$\\mathrm{Beta}(5, 15)$, 3 орла\"),\n",
    "         ((0.5, 0.5), (3, 0), r\"Распределение Джеффриса $\\mathrm{Beta}(0.5, 0.5)$, 3 орла\"),\n",
    "         ((1, 1), (3, 0), r\"Равномерное $\\mathrm{Beta}(1, 1)$, 3 орла\")]\n",
    "\n",
    "fig, axes = plt.subplots(3, 1, figsize=(11, 18))\n",
    "for ax, (prior, dataset, title) in zip(axes, cases):\n",
    "    plot_bernoulli_inference(prior, dataset, ax=ax, ylim=(-0.3, 5), title=title)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-010",
   "metadata": {},
   "source": [
    "### Последовательное обновление"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-011",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "\n",
    "TRUE_THETA = 0.7\n",
    "N_TOSSES = 200\n",
    "tosses = np.random.rand(N_TOSSES) < TRUE_THETA\n",
    "\n",
    "a0, b0 = 2.0, 2.0\n",
    "xs = np.linspace(0, 1, 500)\n",
    "\n",
    "## По одному наблюдению за раз\n",
    "a, b = a0, b0\n",
    "sequential = []\n",
    "for t in tosses:\n",
    "    a, b = (a + 1, b) if t else (a, b + 1)\n",
    "    sequential.append((a, b))\n",
    "\n",
    "## Одной формулой\n",
    "h, t_ = int(tosses.sum()), int((~tosses).sum())\n",
    "a_batch, b_batch = a0 + h, b0 + t_\n",
    "\n",
    "print(\"Beta(%.1f, %.1f)\" % sequential[-1])\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "for n in [0, 1, 2, 5, 10, 25, 50, 200]:\n",
    "    aa, bb = (a0, b0) if n == 0 else sequential[n - 1]\n",
    "    ax.plot(xs, st.beta(aa, bb).pdf(xs), linewidth=2, label=r\"$n = %d$\" % n)\n",
    "ax.axvline(TRUE_THETA, color=\"black\", linestyle=\"--\", linewidth=1.5,\n",
    "           label=r\"истинное $\\theta = %.1f$\" % TRUE_THETA)\n",
    "ax.set_xlabel(r\"Вероятность орла $\\theta$\", fontsize=legend_fontsize)\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize - 4, ncol=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-012",
   "metadata": {},
   "source": [
    "### Точечные оценки и доверительный интервал"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "a_post, b_post = a_batch, b_batch\n",
    "\n",
    "map_est  = (a_post - 1) / (a_post + b_post - 2)\n",
    "mean_est = a_post / (a_post + b_post)\n",
    "mle      = h / N_TOSSES\n",
    "lo, hi   = st.beta(a_post, b_post).interval(0.95)\n",
    "\n",
    "print(\"истинное theta         : %.4f\" % TRUE_THETA)\n",
    "print(\"MLE (h/n)              : %.4f\" % mle)\n",
    "print(\"MAP (мода)             : %.4f\" % map_est)\n",
    "print(\"Апостериорное среднее  : %.4f\" % mean_est)\n",
    "print(\"95%% доверительный инт. : [%.4f, %.4f]\" % (lo, hi))\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(xs, st.beta(a_post, b_post).pdf(xs), linewidth=2,\n",
    "        label=r\"Апостериорное распределение\")\n",
    "mask = (xs >= lo) & (xs <= hi)\n",
    "ax.fill_between(xs[mask], 0, st.beta(a_post, b_post).pdf(xs[mask]),\n",
    "                alpha=.25, label=r\"95\\% доверительный интервал\")\n",
    "ax.axvline(TRUE_THETA, color=\"black\", linestyle=\"--\", linewidth=1.5,\n",
    "           label=r\"истинное $\\theta$\")\n",
    "ax.axvline(map_est, color=palette[3], linewidth=1.5, label=r\"MAP\")\n",
    "ax.set_xlim((0.5, 0.9))\n",
    "ax.set_xlabel(r\"Вероятность орла $\\theta$\", fontsize=legend_fontsize)\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize - 2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-014",
   "metadata": {},
   "source": [
    "## Hot hand fallacy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-015",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "\n",
    "streak, N, num_experiments = 3, 20, 10000\n",
    "results, results_streak = [], []\n",
    "for _ in range(num_experiments):\n",
    "    x = np.random.randint(2, size=N)\n",
    "    after_streak = np.array([i + streak for i in range(N)\n",
    "                             if np.all(x[i:i + streak] == 1) and (i + streak < N)])\n",
    "    if len(after_streak) > 0:\n",
    "        results_streak.append(np.sum(x[after_streak]) / len(after_streak))\n",
    "    results.append(np.sum(x) / float(N))\n",
    "\n",
    "print(\"Доля орлов вообще            : %.4f\" % np.mean(results))\n",
    "print(\"Доля орлов после %d орлов подряд: %.4f\" % (streak, np.mean(results_streak)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-016",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Серия длины 2...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Серия длины 3...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Серия длины 4...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Серия длины 5...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "готово\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(SEED)\n",
    "\n",
    "num_experiments = 25000\n",
    "res = {}\n",
    "for streak in range(1, 6):\n",
    "    print(\"Серия длины %d...\" % streak)\n",
    "    res[streak] = {}\n",
    "    for N in [i for i in range(streak + 1, 20)] + [i for i in range(20, 101, 5)]:\n",
    "        x = np.random.randint(2, size=(num_experiments, N))\n",
    "        after_streak = [[i + streak for i in range(N)\n",
    "                         if np.all(x[j, i:i + streak] == 1) and (i + streak < N)]\n",
    "                        for j in range(num_experiments)]\n",
    "        after_streak_avg = [np.mean(x[j, a]) for (j, a) in enumerate(after_streak) if len(a) > 0]\n",
    "        res[streak][N] = np.mean(after_streak_avg)\n",
    "print(\"готово\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-017",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "for streak in range(1, 6):\n",
    "    ks = sorted(res[streak].keys())\n",
    "    ax.plot(ks, [res[streak][k] for k in ks], linewidth=2,\n",
    "            label=\"После %d орл\" % streak + (\"а\" if streak == 1 else \"ов\"))\n",
    "ax.axhline(0.5, color=\"black\", linestyle=\"--\", linewidth=1.5,\n",
    "           label=r\"«наивный» ответ $0.5$\")\n",
    "ax.set_xlim((0, 100))\n",
    "ax.set_ylim((.3, 0.55))\n",
    "ax.set_xlabel(r\"Длина последовательности $N$\", fontsize=legend_fontsize)\n",
    "ax.set_ylabel(r\"Доля орлов после серии\", fontsize=legend_fontsize)\n",
    "ax.legend(loc=\"lower right\", ncol=3, fontsize=legend_fontsize - 4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-018",
   "metadata": {},
   "source": [
    "# 2. Линейная и полиномиальная регрессия"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-019",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(6754)\n",
    "\n",
    "## Исходная функция\n",
    "orig = lambda x: np.sin(2 * x)\n",
    "\n",
    "SIGMA_NOISE = .25\n",
    "\n",
    "## Небольшая выборка\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",
    "data = orig(xd) + np.random.normal(0, SIGMA_NOISE, num_points)\n",
    "\n",
    "## Большая выборка\n",
    "xd_large = np.arange(-1.5, 2, 0.05)\n",
    "num_points_l = len(xd_large)\n",
    "data_large = orig(xd_large) + np.random.normal(0, SIGMA_NOISE, num_points_l)\n",
    "\n",
    "## Для рисования\n",
    "xs = np.arange(xd[0] - 1.5, xd[-1] + 1.5, 0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-020",
   "metadata": {},
   "source": [
    "## Оверфиттинг"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-021",
   "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, 2, 3, 5, 8,  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 = linear_model.LinearRegression(fit_intercept=True).fit(xd_d[:, :d], data)\n",
    "        print(\"d = %2d, коэффициенты: %s\" % (d, np.array2string(cur_model.coef_, precision=2)))\n",
    "        ax.plot(xs, cur_model.predict(xs_d[:, :d]), linewidth=2, label=\"$d=%d$\" % d)\n",
    "\n",
    "ax.legend(loc=\"upper right\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-022",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Обусловленность матрицы плана\n",
    "for d in [1, 3, 5, 10]:\n",
    "    print(\"d = %2d:  cond(X) = %.3e\" % (d, np.linalg.cond(xd_d[:, :d])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-023",
   "metadata": {},
   "source": [
    "## Локальные признаки в линейной регрессии"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-024",
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "        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",
    "\n",
    "def plural_features(n):\n",
    "    if 11 <= n % 100 <= 14:\n",
    "        return \"признаков\"\n",
    "    return {1: \"признак\", 2: \"признака\", 3: \"признака\", 4: \"признака\"}.get(n % 10, \"признаков\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-025",
   "metadata": {},
   "outputs": [],
   "source": [
    "nums_gauss = [2 ]\n",
    "gauss_xd, gauss_yd = xd, data\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), linear_model.LinearRegression())\n",
    "    gauss_model.fit(gauss_xd[:, np.newaxis], gauss_yd)\n",
    "    yfit = gauss_model.predict(xs[:, np.newaxis])\n",
    "    coefs = gauss_model.get_params()['linearregression'].coef_\n",
    "    print(\"%d гауссовских %s, коэффициенты: %s\"\n",
    "          % (num_gauss, plural_features(num_gauss), \" \".join(\"%.4f\" % x for x in coefs)))\n",
    "    ax.plot(xs, yfit, linewidth=2,\n",
    "            label=\"%d гауссовских %s\" % (num_gauss, plural_features(num_gauss)))\n",
    "\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-026",
   "metadata": {},
   "source": [
    "### Из чего складывается предсказание"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-027",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_gauss = 10\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",
    "\n",
    "print(\"%d гауссовских %s, коэффициенты: %s\"\n",
    "      % (num_gauss, plural_features(num_gauss), \" \".join(\"%.4f\" % x for x in mregr.coef_)))\n",
    "print(\"Свободный член: %.4f\" % mregr.intercept_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-028",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlim((xs[0], xs[-1]))\n",
    "ax.scatter(gauss_xd, gauss_yd, marker='*', s=120)\n",
    "\n",
    "for i in range(mfeat.N):\n",
    "    cur_yfit = mregr.coef_[i] * np.array(\n",
    "        [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",
    "            label=\"Один взвешенный признак\" if i == 0 else None)\n",
    "ax.axhline(mregr.intercept_, color=\"0.6\", linewidth=1, label=\"Свободный член\")\n",
    "ax.plot(xs, yfit, linewidth=2, label=\"Регрессия с гауссовскими признаками\")\n",
    "ax.legend(loc=\"upper center\", fontsize=legend_fontsize - 2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-029",
   "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",
    "\n",
    "ax.plot(xs, yfit, linewidth=2, label=\"Регрессия\")\n",
    "ax.legend(loc=\"upper center\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-030",
   "metadata": {},
   "source": [
    "## Добавим ещё данных"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-031",
   "metadata": {},
   "outputs": [],
   "source": [
    "xs_big = np.arange(xd_large[0] - .5, xd_large[-1] + .5, 0.01)\n",
    "xs_d_big = np.vstack([xs_big ** 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_big[0], xs_big[-1]))\n",
    "ax.set_ylim((-2, 2))\n",
    "ax.scatter(xd_large, data_large, marker='*', s=120)\n",
    "ax.plot(xs_big, orig(xs_big), linewidth=2, label=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "for d in [1, 3, 10]:\n",
    "    cur_model = linear_model.LinearRegression(fit_intercept=True).fit(xd_d_large[:, :d], data_large)\n",
    "    ax.plot(xs_big, cur_model.predict(xs_d_big[:, :d]), linewidth=2, label=\"$d=%d$\" % d)\n",
    "\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-032",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_gauss = 20\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlim((xs_big[0], xs_big[-1]))\n",
    "ax.scatter(xd_large, data_large, marker='*', s=120)\n",
    "\n",
    "gauss_model = make_pipeline(GaussianFeatures(num_gauss), linear_model.LinearRegression())\n",
    "gauss_model.fit(xd_large[:, np.newaxis], data_large)\n",
    "yfit_big = gauss_model.predict(xs_big[:, np.newaxis])\n",
    "mfeat_b = gauss_model.get_params()['gaussianfeatures']\n",
    "mregr_b = gauss_model.get_params()['linearregression']\n",
    "\n",
    "for i in range(mfeat_b.N):\n",
    "    cur_yfit = mregr_b.coef_[i] * np.array(\n",
    "        [mfeat_b._gauss_basis(x, mfeat_b.centers_[i], mfeat_b.width_) for x in xs_big])\n",
    "    ax.plot(xs_big, cur_yfit, color=\"0.4\", linewidth=1)\n",
    "ax.axhline(mregr_b.intercept_, color=\"0.6\", linewidth=1)\n",
    "\n",
    "ax.plot(xs_big, yfit_big, color=\"C1\", linewidth=2, label=\"Регрессия\")\n",
    "ax.plot(xs_big, orig(xs_big), linewidth=1, color=\"black\", label=\"Исходная функция\")\n",
    "ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-033",
   "metadata": {},
   "source": [
    "## Регуляризация"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-034",
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(features, ys, alpha, use_lasso):\n",
    "    if alpha == 0:\n",
    "        return linear_model.LinearRegression(fit_intercept=True).fit(features, ys)\n",
    "    if use_lasso:\n",
    "        return linear_model.Lasso(alpha=alpha, fit_intercept=True, max_iter=100000).fit(features, ys)\n",
    "    return linear_model.Ridge(alpha=alpha, fit_intercept=True).fit(features, ys)\n",
    "\n",
    "\n",
    "def plot_regularization(alpha_values, use_lasso, d=10):\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=\"Исходная функция\", color=\"black\")\n",
    "\n",
    "    for alpha in alpha_values:\n",
    "        m = train_model(xd_d[:, :d], data, alpha, use_lasso)\n",
    "        print(\"alpha = %-9g -> %s\" % (alpha, np.array2string(m.coef_, precision=3)))\n",
    "        ax.plot(xs, m.predict(xs_d[:, :d]), linewidth=2, label=r\"$\\alpha=%g$\" % alpha)\n",
    "\n",
    "    ax.set_title(\"Lasso ($L_1$)\" if use_lasso else \"Ridge ($L_2$)\", fontsize=legend_fontsize)\n",
    "    ax.legend(loc=\"upper center\", fontsize=legend_fontsize)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-035",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_regularization([0, 1e-6, 1e-3, 1.], use_lasso=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-036",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_regularization([0, .01, 1.], use_lasso=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-037",
   "metadata": {},
   "source": [
    "## Усреднение предсказаний"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-038",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "\n",
    "N = 100\n",
    "alpha = 1e-5\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",
    "for _ in range(N):\n",
    "    cur_data = orig(xd) + np.random.normal(0, SIGMA_NOISE, num_points)\n",
    "    cur_model = train_model(xd_d, cur_data, alpha, use_lasso)\n",
    "    res.append(cur_model.predict(xs_d))\n",
    "    ax.plot(xs, res[-1], linewidth=.1, color=\"0.3\")\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",
    "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-039",
   "metadata": {},
   "source": [
    "## Эквивалентное ядро"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-040",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_IND = 7\n",
    "\n",
    "def get_one_prediction(x_pred, cur_y, d=1, data_ind=DATA_IND):\n",
    "    new_data = np.copy(data)\n",
    "    new_data[data_ind] = cur_y\n",
    "    m = linear_model.LinearRegression(fit_intercept=True).fit(xd_d[:, :d], new_data)\n",
    "    return m.predict(np.array([[x_pred ** i for i in range(1, d + 1)]]))[0]\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=figsize)\n",
    "ax = fig.add_subplot(111)\n",
    "ys_grid = np.arange(-1, 1, 0.05)\n",
    "ax.set_xlim((ys_grid[0], ys_grid[-1]))\n",
    "\n",
    "for x_pred in [-.5, 1, 1.5, 2]:\n",
    "    one_pred = [get_one_prediction(x_pred, y) for y in ys_grid]\n",
    "    slope = np.polyfit(ys_grid, one_pred, 1)[0]\n",
    "    ax.plot(ys_grid, one_pred, linewidth=2,\n",
    "            label=r\"$x_* = %.1f$, наклон $k = %.3f$\" % (x_pred, slope))\n",
    "\n",
    "ax.set_xlabel(r\"Значение $y_{%d}$ в обучающей точке $x_{%d} = %.2f$\"\n",
    "              % (DATA_IND, DATA_IND, xd[DATA_IND]), fontsize=legend_fontsize)\n",
    "ax.set_ylabel(r\"Предсказание $\\hat{y}(x_*)$\", fontsize=legend_fontsize)\n",
    "ax.legend(fontsize=legend_fontsize - 2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-041",
   "metadata": {},
   "source": [
    "# 3. Байесовский вывод в линейной регрессии"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-042",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "\n",
    "## Теперь восстанавливаем прямую\n",
    "TRUE_W = np.array([-0.5, 0.5])\n",
    "orig = lambda x: TRUE_W[0] + TRUE_W[1] * x\n",
    "\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",
    "data = orig(xd) + np.random.normal(0, SIGMA_NOISE, num_points)\n",
    "\n",
    "xs = np.linspace(-3, 3, 250)\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.plot(xs, orig(xs), linewidth=2, label=\"Правильный ответ\")\n",
    "ax.scatter(xd, data, marker='*', s=120, label=\"Данные\")\n",
    "ax.legend(fontsize=legend_fontsize)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-043",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Сетка в пространстве весов (w0, w1)\n",
    "N_GRID = 250\n",
    "W0, W1 = np.meshgrid(np.linspace(-1, 1, N_GRID), np.linspace(-1, 1, N_GRID))\n",
    "pos = np.dstack((W0, W1))\n",
    "\n",
    "\n",
    "def myplot_heatmap(Z, title=None, points=None):\n",
    "    fig = plt.figure(figsize=(6, 6))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.pcolormesh(W0, W1, Z, cmap=plt.cm.jet, shading=\"auto\")\n",
    "    ax.scatter([TRUE_W[0]], [TRUE_W[1]], marker='*', s=200, color=\"white\",\n",
    "               edgecolors=\"black\", zorder=5)\n",
    "    ax.set_xlim((-1, 1))\n",
    "    ax.set_ylim((-1, 1))\n",
    "    ax.set_aspect('equal', adjustable='box')\n",
    "    ax.set_xlabel(r\"$w_0$\", fontsize=legend_fontsize)\n",
    "    ax.set_ylabel(r\"$w_1$\", fontsize=legend_fontsize)\n",
    "    if title is not None:\n",
    "        ax.set_title(title, fontsize=legend_fontsize)\n",
    "    ax.grid(False)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def myplot_sample_lines(mu, sigma, n=20, points=None, title=None):\n",
    "    my_w = np.random.multivariate_normal(mu, sigma, n)\n",
    "    fig = plt.figure(figsize=figsize)\n",
    "    ax = fig.add_subplot(111)\n",
    "    for w in my_w:\n",
    "        ax.plot(xs, w[0] + w[1] * xs, 'k-', lw=.4)\n",
    "    ax.plot(xs, orig(xs), linewidth=2, color=palette[0], label=\"Правильный ответ\")\n",
    "    ax.set_ylim((-3, 3))\n",
    "    ax.set_xlim((-3, 3))\n",
    "    if points is not None:\n",
    "        ax.scatter(points[0], points[1], marker='*', s=200, zorder=5)\n",
    "    if title is not None:\n",
    "        ax.set_title(title, fontsize=legend_fontsize)\n",
    "    ax.legend(loc=\"upper left\", fontsize=legend_fontsize)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-044",
   "metadata": {},
   "source": [
    "## Априорное распределение"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-045",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "\n",
    "PRIOR_MU = np.array([0., 0.])\n",
    "PRIOR_SIGMA = 2 * np.eye(2)\n",
    "\n",
    "cur_mu, cur_sigma = PRIOR_MU.copy(), PRIOR_SIGMA.copy()\n",
    "\n",
    "Z = multivariate_normal.pdf(pos, mean=cur_mu, cov=cur_sigma)\n",
    "myplot_heatmap(Z, title=r\"Априорное распределение $p(w)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-046",
   "metadata": {},
   "outputs": [],
   "source": [
    "myplot_sample_lines(cur_mu, cur_sigma, 200, title=r\"Прямые из априорного распределения\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-047",
   "metadata": {},
   "source": [
    "## Правдоподобие одной точки"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-048",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_likelihood(px, py, sigma=SIGMA_NOISE):\n",
    "    return lambda w: (np.exp(-((w[..., 0] + w[..., 1] * px - py) ** 2) / (2 * sigma ** 2))\n",
    "                      / (sigma * np.sqrt(2. * np.pi)))\n",
    "\n",
    "\n",
    "def likelihood_grid(px, py):\n",
    "    return get_likelihood(px, py)(pos)\n",
    "\n",
    "\n",
    "px, py = xd[5], data[5]\n",
    "print(\"Первое наблюдение: x = %.2f, y = %.4f\" % (px, py))\n",
    "myplot_heatmap(likelihood_grid(px, py), title=r\"Правдоподобие одной точки\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-049",
   "metadata": {},
   "source": [
    "## Байесовское обновление"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-050",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bayesian_update(mu, sigma, x, y, sigma_noise=SIGMA_NOISE):\n",
    "    x_matrix = np.array([[1., x]])\n",
    "    precision = np.linalg.inv(sigma) + (1 / sigma_noise ** 2) * (x_matrix.T @ x_matrix)\n",
    "    sigma_n = np.linalg.inv(precision)\n",
    "    mu_n = sigma_n @ (np.linalg.inv(sigma) @ mu + (1 / sigma_noise ** 2) * x_matrix.T @ np.array([y]))\n",
    "    return mu_n, sigma_n\n",
    "\n",
    "\n",
    "cur_mu, cur_sigma = bayesian_update(cur_mu, cur_sigma, px, py)\n",
    "print(\"mu    =\", np.array2string(cur_mu, precision=4))\n",
    "print(\"sigma =\", np.array2string(cur_sigma, precision=4).replace(\"\\n\", \"\\n        \"))\n",
    "\n",
    "Z = multivariate_normal.pdf(pos, mean=cur_mu, cov=cur_sigma)\n",
    "myplot_heatmap(Z, title=r\"Апостериорное распределение после 1 точки\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-051",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "myplot_sample_lines(cur_mu, cur_sigma, 40, points=[[px], [py]],\n",
    "                    title=r\"Прямые из апостериорного распределения, 1 точка\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-052",
   "metadata": {},
   "source": [
    "## Предсказательное распределение"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-053",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_statistics(mu, sigma, xs, n=2000):\n",
    "    my_w = np.random.multivariate_normal(mu, sigma, n)\n",
    "    return my_w[:, 0][:, None] + my_w[:, 1][:, None] * xs[None, :]\n",
    "\n",
    "\n",
    "def plot_predictions(xs, mu, preds, points, title=None):\n",
    "    mean_pred = mu[0] + mu[1] * xs\n",
    "    std_pred = np.std(preds, axis=0)\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.plot(xs, orig(xs), label=\"Правильный ответ\")\n",
    "    ax.plot(xs, mean_pred, color=\"red\", label=\"MAP гипотеза\")\n",
    "    ax.fill_between(xs, mean_pred - SIGMA_NOISE, mean_pred + SIGMA_NOISE,\n",
    "                    color=palette[1], alpha=.3, label=r\"$\\pm$ дисперсия шума\")\n",
    "    ax.fill_between(xs, mean_pred - std_pred - SIGMA_NOISE, mean_pred + std_pred + SIGMA_NOISE,\n",
    "                    color=palette[5], alpha=.2, label=r\"$\\pm$ дисперсия предсказаний\")\n",
    "    ax.scatter(points[0], points[1], marker='*', s=200, zorder=5)\n",
    "    if title is not None:\n",
    "        ax.set_title(title, fontsize=legend_fontsize)\n",
    "    ax.legend(fontsize=legend_fontsize - 2)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "np.random.seed(SEED)\n",
    "preds = sample_statistics(cur_mu, cur_sigma, xs, n=2000)\n",
    "plot_predictions(xs, cur_mu, preds, [[px], [py]], title=r\"После 1 точки\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-054",
   "metadata": {},
   "source": [
    "## Вторая точка"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-055",
   "metadata": {},
   "outputs": [],
   "source": [
    "px2, py2 = xd[7], data[7]\n",
    "print(\"Второе наблюдение: x = %.2f, y = %.4f\" % (px2, py2))\n",
    "myplot_heatmap(likelihood_grid(px2, py2), title=r\"Правдоподобие второй точки\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-056",
   "metadata": {},
   "outputs": [],
   "source": [
    "cur_mu, cur_sigma = bayesian_update(cur_mu, cur_sigma, px2, py2)\n",
    "Z = multivariate_normal.pdf(pos, mean=cur_mu, cov=cur_sigma)\n",
    "myplot_heatmap(Z, title=r\"Апостериорное распределение после 2 точек\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-057",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "myplot_sample_lines(cur_mu, cur_sigma, n=40, points=[[px, px2], [py, py2]],\n",
    "                    title=r\"Прямые из апостериорного распределения, 2 точки\")\n",
    "\n",
    "preds = sample_statistics(cur_mu, cur_sigma, xs, n=2000)\n",
    "plot_predictions(xs, cur_mu, preds, [[px, px2], [py, py2]], title=r\"После 2 точек\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-058",
   "metadata": {},
   "source": [
    "## Третья точка"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-059",
   "metadata": {},
   "outputs": [],
   "source": [
    "px3, py3 = xd[1], data[1]\n",
    "print(\"Третье наблюдение: x = %.2f, y = %.4f\" % (px3, py3))\n",
    "myplot_heatmap(likelihood_grid(px3, py3), title=r\"Правдоподобие третьей точки\")\n",
    "\n",
    "cur_mu, cur_sigma = bayesian_update(cur_mu, cur_sigma, px3, py3)\n",
    "Z = multivariate_normal.pdf(pos, mean=cur_mu, cov=cur_sigma)\n",
    "myplot_heatmap(Z, title=r\"Апостериорное распределение после 3 точек\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-060",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "myplot_sample_lines(cur_mu, cur_sigma, n=200, points=[[px, px2, px3], [py, py2, py3]],\n",
    "                    title=r\"Прямые из апостериорного распределения, 3 точки\")\n",
    "\n",
    "preds = sample_statistics(cur_mu, cur_sigma, xs, n=2000)\n",
    "plot_predictions(xs, cur_mu, preds, [[px, px2, px3], [py, py2, py3]],\n",
    "                 title=r\"После 3 точек\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-061",
   "metadata": {},
   "source": [
    "## Все точки сразу"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-062",
   "metadata": {},
   "outputs": [],
   "source": [
    "cur_mu, cur_sigma = PRIOR_MU.copy(), PRIOR_SIGMA.copy()\n",
    "for x_i, y_i in zip(xd, data):\n",
    "    cur_mu, cur_sigma = bayesian_update(cur_mu, cur_sigma, x_i, y_i)\n",
    "\n",
    "print(\"истинные веса         :\", np.array2string(TRUE_W, precision=4))\n",
    "print(\"апостериорное среднее :\", np.array2string(cur_mu, precision=4))\n",
    "print(\"апостериорные СКО     :\", np.array2string(np.sqrt(np.diag(cur_sigma)), precision=4))\n",
    "\n",
    "ols = linear_model.LinearRegression().fit(xd[:, None], data)\n",
    "print(\"\\nМНК (для сравнения)   : [%.4f %.4f]\" % (ols.intercept_, ols.coef_[0]))\n",
    "\n",
    "Z = multivariate_normal.pdf(pos, mean=cur_mu, cov=cur_sigma)\n",
    "myplot_heatmap(Z, title=r\"Апостериорное распределение по всем %d точкам\" % num_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-063",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "myplot_sample_lines(cur_mu, cur_sigma, n=200, points=[xd, data],\n",
    "                    title=r\"Прямые из апостериорного распределения, все точки\")\n",
    "\n",
    "preds = sample_statistics(cur_mu, cur_sigma, xs, n=2000)\n",
    "plot_predictions(xs, cur_mu, preds, [xd, data], title=r\"После всех %d точек\" % num_points)"
   ]
  }
 ],
 "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}