Lecture 6. Monte-Carlo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2d8bea35-37cc-4b98-bdeb-471a07352aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import (\n",
    "    norm, binom, expon, t, chi2, pareto, ttest_ind, sem, beta, laplace\n",
    ")\n",
    "from statsmodels.stats.proportion import proportion_confint\n",
    "import numpy as numpy\n",
    "from seaborn import distplot\n",
    "from matplotlib import pyplot\n",
    "import seaborn\n",
    "\n",
    "import sys\n",
    "sys.path.append('.')\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ac266211-93ed-4151-ab04-01a8598dff79",
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse_plot_colorscheme():\n",
    "    import cycler\n",
    "    def invert(color_to_convert): \n",
    "        table = str.maketrans('0123456789abcdef', 'fedcba9876543210')\n",
    "        return '#' + color_to_convert[1:].lower().translate(table).upper()\n",
    "    update_dict = {}\n",
    "    for key, value in pyplot.rcParams.items():\n",
    "        if value == 'black':\n",
    "            update_dict[key] = 'white'\n",
    "        elif value == 'white':\n",
    "            update_dict[key] = 'black'\n",
    "    \n",
    "    old_cycle = pyplot.rcParams['axes.prop_cycle']\n",
    "    new_cycle = []\n",
    "    for value in old_cycle:\n",
    "        new_cycle.append({\n",
    "            'color': invert(value['color'])\n",
    "        })\n",
    "    pyplot.rcParams.update(update_dict)\n",
    "    pyplot.rcParams['axes.prop_cycle'] = cycler.Cycler(new_cycle)\n",
    "    lec = pyplot.rcParams['legend.edgecolor']\n",
    "    lec = str(1 - float(lec))\n",
    "    pyplot.rcParams['legend.edgecolor'] = lec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d6bcb93b-bbed-4415-aa28-9e5c6052d51c",
   "metadata": {},
   "outputs": [],
   "source": [
    "inverse_plot_colorscheme()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1276fbd9-51d6-4da1-8a4f-93f8d7d93a28",
   "metadata": {},
   "source": [
    "# Лекция 6. Монте-Карло\n",
    "\n",
    "Сегодня мы поговорим о методе Монте-Карло, являющимся очень мощным инструментом в статистике. С его помощью мы ответим с вами на 3 вопроса:\n",
    "- Как проверить ваш критерий? Валиден он на практике или нет?\n",
    "    - Например, работает ли t-test на малых размерах выборок?\n",
    "- У вас есть 2 разных критерия. Как понять, какой критерий лучше подходит для вашей задачи?\n",
    "<!-- - Как жить, если мы не знаем теор.вер? Как придумывать критерии на практике, если мы не хотим ничего считатать? В некоторых случаях у вас появится ответ на этот вопрос) -->"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "068ded1b-f874-43a4-9f24-4f32c9db8353",
   "metadata": {},
   "source": [
    "## 1. Проверка критериев\n",
    "\n",
    "С помощью метода Монте-Карло мы в *общем случае* сможем ответить на вопросы:\n",
    "- **Можно ли использовать данный критерий для нашей задачи?**\n",
    "- **Верно ли вообще реализован критерий?**\n",
    "\n",
    "Вся эта глава в первую очередь будет посвящена AB-тестам и как можно проверять критерии для них. Основным критерием для проверки в этой главе станет t-test. Мы с вами:\n",
    "- Покажем на практике, что t-test работает для выборок не только из нормального распределения\n",
    "- Посмотрим, как определить, с какого размера выборки можно применять t-test. Как мы помним из прошлой лекции, t-test работает теоретически для выборок из любого распределения, если выборка достаточно большая. \n",
    "\n",
    "\n",
    "|                          | маленькая выборка | большая выборка |\n",
    "|--------------------------|-------------------|-----------------|\n",
    "| нормальное распределение | t-test            | t-test |\n",
    "| любое распределение      |                   | t-test |\n",
    "\n",
    "---\n",
    "\n",
    "Что значит, что критерий корректен? Давайте пойдем от определения.\n",
    "- Критерий уровня значимости $\\alpha$ означает, что вероятность неверно отвергнуть нулевую гипотезу $\\le \\alpha$. \n",
    "- А это в свою очередь значит, что если бесконечно много раз повторить один и тот же эксперимент, в котором верна нулевая гипотеза, генерируя заново эксперимент, то число ложноположительных срабатываний будет меньше $\\alpha$ процентов.\n",
    "\n",
    "На самом деле, здесь уже расказана процедура, как проверить ваш критерий :)\n",
    "\n",
    "0. Создаем код критерия, который мы будем проверять.\n",
    "1. Генерируем как можно больше экспериментов, где верна $H_0$. \n",
    "2. Прогоняем на них придуманный критерий.\n",
    "3. Проверяем, правда ли, что только в $\\alpha$ процентов случаев критерий отвергается?\n",
    "\n",
    "А теперь распишем подробнее:\n",
    "1. Первым делом надо выбрать распределение, которое будет описывать наши данные. К примеру, если у нас метрика конверсии, то это бернуллевское распределение, а если метрика — выручка, то лучше использовать экспоненциальное распределение в качестве самого простого приближения.\n",
    "\n",
    "2. Завести счётчик bad_cnt = 0.\n",
    "\n",
    "3. Далее в цикле размера N, где N — натуральное число от 1000 до бесконечности (чем оно больше, тем лучше):\n",
    "   - Симулировать создание выборки из распределения, выбранного на первом шаге. Так, чтобы верна была $H_0$.\n",
    "        - А в случае AB-теста симулировать надо не 1 выборку, а 2: для теста и контроля.\n",
    "\n",
    "    - Запустить на сгенерированных данных проверяемый критерий.\n",
    "\n",
    "    - Далее проверить, `pvalue < alpha`. Если да, то увеличить счётчик bad_cnt на 1. Здесь мы проверяем, ошибся ли критерий на текущей симуляции, или нет. \n",
    "\n",
    "4. Посчитать конверсию bad_cnt / N.\n",
    "    - Если она примерно совпадает с $\\alpha$, то все хорошо.\n",
    "    - Если она меньше $\\alpha$, то в принципе это адекватный критерий на практике, просто он будет менее мощный, чем критерий, который ошибается ровно в $\\alpha$ проценте случаев.\n",
    "        - Но на практике стоит проверить: а теоретически такая ситуация возможна? Или это ошибка в коде критерия?\n",
    "    - Если критерий ошибается больше, чем в $\\alpha$, то значит он некорректен и им нельзя пользоваться. Используя такой критерий, вы будете ошибаться чаще положенного, а значит ваша  компания будет терять больше денег.\n",
    "    \n",
    "    \n",
    "Рассмотрим процедуру на примере: проверим, можно ли использовать t-test для выборок из нормального распределения?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d7e9da0f-3161-49da-97ee-252c2619a077",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FPR: 0.0519\n"
     ]
    }
   ],
   "source": [
    "numpy.random.seed(42)\n",
    "\n",
    "bad_cnt = 0\n",
    "N = 10000\n",
    "alpha=0.05\n",
    "\n",
    "sample_dist = norm(loc=2, scale=3)\n",
    "mu0=sample_dist.expect()\n",
    "for i in range(N):\n",
    "    # Генерирую выборку теста и контроля\n",
    "    test    = sample_dist.rvs(5)\n",
    "    control = sample_dist.rvs(5)\n",
    "\n",
    "    # Запускаю критерий и считаю p-value\n",
    "    pvalue = ttest_ind(test, control, alternative='two-sided').pvalue\n",
    "    \n",
    "    # Проверяю, что pvalue < alpha\n",
    "    bad_cnt += (pvalue < alpha)\n",
    "\n",
    "\n",
    "print(f\"FPR: {round(bad_cnt / N, 4)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05713d6c-0e68-4ee3-8e96-151392c61ef9",
   "metadata": {},
   "source": [
    "Хм, мы получили, что FPR=0.0519, хотя он должен был равняться 5%. Правда ли, что критерий некорректен? Ну конечно нет, мы просто не учли шум  конверсии: мы вряд ли сможем получить на конечном числе экспериментов точное равенство `FPR=alpha`.\n",
    "\n",
    "Поэтому надо чуть улучшить 4 шаг:\n",
    "\n",
    "4. Посчитать полученный FPR и построить доверительный интервал для него. Если $\\alpha$ лежит в нем, значит все хорошо, а иначе разбираемся, что пошло не так.\n",
    "    - Доверительный интервал можно построить разными способами: например, используя идеи построения доверительных интервалов из второй лекции.\n",
    "    - Но можно сделать проще: в питоне есть функция, которая строит [доверительный интервал Уилсона](https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval#Wilson_score_interval): он не такой точный, как мы выводили ранее, зато он более быстрый и работает из \"коробки\". Его не надо реализовывать самому."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8ca48be8-2429-43ee-8ab9-e5828a725997",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.04772180742973847, 0.05642233191006188)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proportion_confint(count = bad_cnt, nobs = N, alpha=0.05, method='wilson')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbb87da4-a54a-4c67-82e0-09a181f1c0ca",
   "metadata": {},
   "source": [
    "Мы видим, что 5% попали в доверительный интервал, а значит мы можем считать, что критерий валиден для нашей задачи.\n",
    "\n",
    "А что, если бы распределение было сложнее?\n",
    "\n",
    "Расссмотрим пример, когда мат. ожидания в тесте и в контроле равны, но выборки из разных распределений. То есть $H_0$ верна, но распределения разные."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "83115ac4-49da-41f6-99c0-ce58a1b0aa8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAFNCAYAAAC5eOMWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABF/klEQVR4nO3deXxU1f3/8ddkspGEkEBYE1ZlERVwAxdULBWFQlFbXCu49It+W5daq7V2r9r228X+6lLXWkGt1LUi4r7vbCKrCARCAkkI2fdt7u+PcwPDMElmwtxZkvfz8biPe+eun8nB5OM5557jsiwLEREREYl+cZEOQEREREQCo8RNREREJEYocRMRERGJEUrcRERERGKEEjcRERGRGKHETURERCRGKHETEZFQGg/MAtzAxcCQyIYj0r0ocRMRXzuBeqDGa2kAPopgTBI7SoBfAPuA+UBZZMMR6V6UuImIP3OANK/l2siGIzGkBDgVyMTUvDVENhyR7kWJm4h0xU7gZ8AmoBz4F5BsH5sGFHideyFgAd+3P19hf77J65xZ9r47vfbNBtYCFcAnwIQAn+/rCuBj4F6gEvgKmO51/EpgM1AN5ALX+Fw/146jCtgOnGvvfw+TlLTVStbbcQUaY0ffD+BJoMnr3t4/0yTgL8AuoBh4EOjldXwE5ufZFlsrB37+ccBt9ncpBZ4B+vpcF+8Tx2/s7Wl0XrbeNbO32se/iYiEhBI3Eemqy4BzgCOAMZjmMV8JwB1Aoc/+bcACr8/fxyRPbY4HHsMkUf2Ah4ClmIQlmOe3mYJJyrKAXwMvcCBZ2YtJotIxSdzf7OcDTAYWA7cAGcAZHJycXceBWsk5fp7bXoyBfD8XcJd975k+9/0/+36TgCOBbOBXXsfbfrf3sa//0OvYDcB5wJmY/mflwP1+Yu9Me2XbJtN+VkUX7i0i7VDiJiJddR+Qj+nDdBdwiZ9zrgE+B7722V+MSYBOAQYAw4EVXsf/B5PMfI6pLVoENAInB/n8NnuB/wc0A/8BtgDfso+9gql9soD3gTeA0+1jV2MSrDcBD7AbU2MXqPZiDOT79cLUuPly2dffZN+3Gvg95kWANol2vK1+rr8G+Dmm5qwRU5v2XQ6uZQtEe2Xb5ueYn11lkPcVkQ4ocRORrsr32s7j0LcHe2Oayn7ZzvWPYmrarsDUankbDtyMqa1pW4b6PKOz53vbjUnM/J0/E/gMkwRVYJpts+xjQzFJXVe1F2Mg328Qpr+Yr/5ACrDa69rX7P1t+mJq0vwZDrzode1mTII30OucfV7HL/Rzj87Kdph93Z/bOS4iXaTETUS6aqjX9jBgj8/xWzD9p/Lauf5V4DRMk+kTPsfyMTVUGV5LCvB0EM/3lo2pqfI9Pwl4HtNfbKD9nOVe5+Zjmjm7qr0YO/t+CcAxwJd+7rkP0+ftaK9r25pE24yh/ZqwfEyy6v3sZExy2ybL69gzfu7RWdneCfwJUxsoIiGkxE1EuuqHQA6mdud2TBNkm96Y/mJ3dXB9K6av1pMcOmTEI5g3WadgkqhUTNNm7wCf72sApr9VAjAPOAqToCVikrcSoAWT0Mzwuu6f9veYjvl9mQ2M6+A5vtqLsbPvdyVQBKzyc0+Pff3f7O+FHdc59vZQ4Ebgv+3E9CCmXIbbn/tjXsAIVGdleyTmez0UxD1FJEBK3ESkq/6N6Q+Way/eb4SmA/fQfnNdm38Bf/CzfxWmH9d99j22YZpUA32+r8+B0ZjaqrswfbpKMTVCN2Bqj8qBSzEvCbRZwYEXFioxfeCGE7j2Yuzo+12GSXpG2vHVYGonh2CSLoCf2td8hnnb9S1grH3sdcwbr39rJ6a/29/xDfv+n2ESrUB1VrYDMS9hNAdxTxEJkMuyrM7PEhE52E5M/7S3YuD5V9jnTnUwHn920rWf0RWYYTl+47M/B5P4XXFYUYlITFONm4hIdKnF1KL5akGzEIj0eMG+/i0iIs56tp39RcCPwxmIiEQfNZWKiIiIxAg1lYqIiIjECCVuIiIiIjGiR/RxKykpsfLy2hsn8vC53W4AWlv9zS4jkaJyiT4qk+ikcok+KpPoFK5yOfHEE/dx8Gwo+/WIxC0vL4+TTjrJsftnZmYCUF7e2ZBVEk4ql+ijMolOKpfoozKJTuEqF8uy2q1tUlOpiIiISIxQ4iYiIiISI5S4iYiIiMQIJW4iIiIiMaJHvJwgIiIih8flcpGVlUVGRsb+tyt7mrg4U981cODAw75XQ0MDBQUFtLS0BHWdEjcRERHpVE5ODpZlsXPnTpqbmyMdTkSEcjiQfv36kZOTw86dO4O6zumm0nOBLcA24DY/xy8D1tnLJ8DEAK7tC7wJbLXXmSGPWkRERA6SmprK7t27e2zSFmqlpaUkJycHfZ2TiZsbuB+YCYwHLrHX3nYAZwITgDuAhwO49jbgbWC0vfaXEIqIiEiIaX7zyHMycZuMqS3LBZqAJcBcn3M+AdpGsfsMyAng2rnAInt7EXBe6EMXERERiT5O9nHLBvK9PhcAUzo4/2rg1QCuHQgU2tuFwIDOAnG73ftHO3ZCRkaGY/f21jLkRNwlm3A114XlebEuXOUigVOZRCeVS/SJxjKJi4vrsS8ltAn194+Liws6P3Gyxs3lZ197daxnYRK3n3bh2vYsBFYBq7KysoK8NPpY7iSqL3iaxmMvi3QoIiIiUWfbtm1Mnz79sO4xf/583n///RBF5Awna9wKgKFen3OAPX7OmwA8iunPVhrAtcXAYExt22BgbzvPf9heKC4utsIx35ujz0hMA3cC9UkDqNfcdUHRXH/RR2USnVQu0SeaymTgwIFRP+l9a2vrYcXo8XiwLKvTe4Tq5+DxeIIuYydr3FZiXiAYCSQCFwNLfc4ZBrwAXA58HeC1S4EF9vYC4CUHYo8+LrsSMn1wZOMQERGJMosXL2bYsGG8/PLLVFdXc8sttzBlyhQ+/vhjysvLWbt2LWeeeeb+8xcsWMD27dupqqoiNzeXSy+9lHHjxvHggw9yyimnUF1dHVVJ80Esy3JymWVZ1teWZW23LOvn9r5r7QXLsh61LKvcsqy19rKqk2uxLKufZVlvW5a11V737SyOlStXWpimVkeWzMxMKzMz09FnkNTb4jeVFgvfc/Y53WgJS7loUZl0g0XlEn1LNJbJuHHjIh5DR8uOHTus6dOnW4A1ZMgQa9++fdbMmTMtl8tlffOb37T27dtnZWVlWSkpKVZlZaU1ZswYC7AGDRpkjR8/3gKsBQsWWB9++GG7z3C73Zbb7Xb8Z+qTDx20OD0A73J78fag1/b37SXQa8E0px5eI3ZMsmvceqvGTUREIu9Xs8czfki6o8/YtKeK3y3bFPR13/ve91i+fDmvvmreeXzrrbdYtWoVs2bN4rnnnsPj8XDMMcewa9cuioqKKCoqCnXojtFcpbEmbQC4EyIdhYiISNQaPnw48+bNo7y8fP8ydepUBg8eTF1dHRdddBHXXnsthYWFLFu2jLFjx0Y65IBpyqtY0dbHzRUHaYOgMr/j80VERBzUlZowJ3kPDpyfn88TTzzBwoUL/Z77xhtv8MYbb5CcnMydd97JI488whlnnBETAwyrxi1meI2QohcUREREDlJcXMyoUaMAePLJJ5kzZw4zZswgLi6OpKQkzjzzTLKzsxkwYABz5swhJSWFxsZGampq9r8lWlxcTE5ODgkJ0duypcQtFqVnRzoCERGRqPKHP/yBX/ziF5SXl3PRRRcxd+5cbr/9dkpKSsjPz+eWW24hLi6OuLg4br75Zvbs2UNZWRlnnnkmP/jBDwB455132LhxI0VFRZSUlET4G/mnptJY4fKucRsSuThERESi0NKlS1m69OBRx6ZNm+b33Pb2Nzc3M3v27BBHFlqqcYtFerNURESkR1LiFovUVCoiItIjKXGLRXo5QUREpEdS4hYr1MdNRESkx1PiFmuaaqD3kIMTOREREekRlLjFDDtRq9xtZk5IyYpsOCIiIhJ2StxiTdUes9YLCiIiIj2OErdY0dY0WrXbrNXPTUREpMdR4hZr9te4KXETERGJZu+++y5XX311SO+pxC1m2DVutXuhtVmJm4iIiIN+/etf88QTT0Q6jEMocYs1Hg/UFGn2BBERkQhzRWCEByVuscL7H0fVHr2cICIi4iUnJ4fnn3+evXv3sm/fPu69915cLhc///nP2blzJ8XFxSxatIj09HQAhg8fjmVZzJ8/n7y8PEpKSrj99tsBOOecc7j99tu56KKLqK6uZu3atQC8/fbb3HHHHXz00UfU1dUxatQoTjnlFFasWEFFRQUrVqzglFNOcfR7KnGLOZaduKmpVEREBCAuLo5ly5aRl5fHiBEjyM7OZsmSJVxxxRVcccUVnHXWWYwaNYq0tDTuu+++g66dOnUqY8eOZfr06fzqV79i3LhxvP766/z+97/nP//5D71792bSpEn7z7/ssstYuHAhvXv3prq6mldeeYV77rmHfv36cffdd/PKK6/Qt29fx75rvGN3lhDzqXEbPSNyoYiIiJz7Bxh0rLPPKFoPr/2s09MmT57MkCFDuOWWW2htbQXg448/5re//S133303O3bsAOBnP/sZGzZs4Morr9x/7W9/+1saGhpYt24dX375JRMnTuSrr75q91mLFy9m06ZNAMyYMYOtW7fy5JNPArBkyRJuuOEG5syZw6JFi7r8tTuiGrdY0dZUallmSJDEVOiVGdmYREREosDQoUPJy8vbn7S1GTJkCHl5efs/5+XlkZCQwMCBA/fvKyoq2r9dV1dHWlpah8/Kz89v9/5tz8jOdq47k2rcYlGF/Y+mTw7Ul0c2FhER6ZkCqAkLl/z8fIYNG4bb7T4oeduzZw/Dhw/f/3nYsGE0NzdTXFxMTk5Oh/e0LKvT/b73b3vGa6+91pWvERDVuMWMtqZSCyoLzGafoRGLRkREJFqsWLGCwsJC/vjHP5KSkkJSUhKnnnoqTz/9NDfddBMjRowgNTV1f78135o5f4qLixkxYkSHb44uX76cMWPGcMkll+B2u7nwwgsZP348y5YtC+XXO4jTidu5wBZgG3Cbn+PjgE+BRuAnXvvHAmu9lirgR/ax3wC7vY7NCmnEsaByl1n36fj/FkRERHoCj8fDnDlzOPLII9m1axcFBQVcdNFFPPbYYzzxxBN88MEH7Nixg4aGBq6//vqA7vnss88CUFpayurVq/2eU1ZWxuzZs7n55pspLS3l1ltvZfbs2ZSWlobsux3CsiynFrdlWdstyxplWVaiZVlfWpY13uecAZZlnWRZ1l2WZf2kg/sUWZY13P78mw7O9busXLnSAhxbMjMzrczMTEefQfoQi99UWhy/wHz+eZHFjDudfWaML2EpFy0qk26wqFyib4nGMhk3blzEY4j04na7Lbfb7fjP1LKsVe3lNE7WuE3G1LTlAk3AEmCuzzl7gZVAcwf3mQ5sB/I6OKfnqSxQjZuIiEgP4+TLCdlAvtfnAmBKF+5zMfC0z77rgPnAKuBmoLyjG7jdbjIznXsDMyMjw7F7t/GkZVAJpKT0Iikzk+q6Iqx+I0l38HvFunCUiwRHZRKdVC7RJxrLJC4uDrfbHekwIirU3z8uLi7o/MTJGjd/vfmsIO+RCHwbeNZr3wPAEcAkoBD4azvXLsQkdquysrKCfGz0i6vajae3BuEVERHpSZyscSsAvF97zAH2BHmPmcAaoNhrn/f2I0B7r248bC8UFxdb5eUdVsqFhKPP8KQCUFdbS115OZRsh2Mupry6DloanXtuNxCOspfgqEyik8ol+kRTmQwcODCgtzF7glD9HDweT9Bl7GSN20pgNDASU3N2MbA0yHtcwqHNpN6zq58PbOhqgDGtbSw3zVkqIiLSYzhZ49aC6Yv2OuAGHgM2Atfaxx8EBmGaM9MBD2bIj/GY4T9SgLOBa3zu+ydMM6kF7PRzvJvyGscNDh7LrSw3IhGJiEjP0dzcTK9evaivr490KN1CQkICLS0tQV/n9MwJy+3F24Ne20WYJlR/6oB+fvZfHoK4YlfbiM2VXrMniIiIOGzv3r1kZ2eze/duJW+HyeVyMXDgQCorK4O+VlNexQrfkZurdoPlgQzNniAiIs6rrq4GzPycCQkJEY4mMuLiTA8zj8dz2Peqra1l3759QV+nxC3m2DVurc1QXaQaNxERCZvq6ur9CVxP1DZ0RyRfGtFcpTHDz+gqlQXQZ1j4QxEREZGIUOIWayyvofAq81XjJiIi0oMocYsVvn3c4MC0V/6OiYiISLejxC3meNW4VeyC+CRI7R+5cERERCRslLjFjHZq3MCM5SYiIiLdnhK3WNHWHOrbxw2UuImIiPQQStxiWdu0VxrLTUREpEdQ4hZzvGrcGqugoQIyNCSIiIhIT6DELdaV50HG8EhHISIiImGgxC1W+OvjBiZxyxwR9nBEREQk/JS4xbrynZA5XGO5iYiI9ABK3GJGW2LmW+O2E+KTIW1guAMSERGRMFPiFuvKd5q1mktFRES6PSVusaLdPm47zFqJm4iISLenxC3WVRaA5VHiJiIi0gMocYsZ7bx80NoEVbs1JIiIiEgPoMQt1vg2lYL9ZumIcEciIiIiYabELVZ0NNyHEjcREZEeQYlbzPFX45YH6UPMsCAiIiLSbSlxixmd1LiBJpsXERHp5pxO3M4FtgDbgNv8HB8HfAo0Aj/xObYTWA+sBVZ57e8LvAlstdeZoQw46rXXxw3UXCoiItLNOZm4uYH7gZnAeOASe+2tDLgB+Es79zgLmASc6LXvNuBtYLS99pcQdj+d9XEDJW4iIiLdnJOJ22RMTVsu0AQsAeb6nLMXWAk0B3HfucAie3sRcN5hRRkz2pnyCqC2BJpqlbiJiIh0c/EO3jsbyPf6XABMCeJ6C3jDXj8EPGzvHwgU2tuFwIDObuR2u8nMdK5FNSMjw7F7t2lNT6cKSE1NJdHPd6msKsA9YAxpDn7PWBOOcpHgqEyik8ol+qhMolM0lIuTiZu/tj0/1UXtOg3Yg0nM3gS+Aj4I4vqF9kJWVlYQl0Uny9VBjRvgrtqFJ10vJ4iIiHRnTiZuBYB3JpGDScQC1XbuXuBFTNPrB0AxMBhT2zbYPu7Pw/ZCcXGxVV5eHsSju8bRZyRWAVBbU0utv+fs3QrZJzsbQ4zSzyT6qEyik8ol+qhMolMky8XJPm4rMS8QjAQSgYuBpQFemwr09tqeAWywPy8FFtjbC4CXQhFs9Ou4xo3ynZDUG1Jjv3ZRRERE/HOyxq0FuA54HfOG6WPARuBa+/iDwCDMUB/pgAf4EebN0yxMLVtbjP8GXrM//xF4Brga2AXMc/A7xI6yXLPuewTU7otsLCIiIuIIJxM3gOX24u1Br+0iTBOqrypgYjv3LAWmH35oMaatj5u/cdwASrebdb8jIP/z8MQkIiIiYaWZE7qLijxobTY1biIiItItKXGLGZ30cfO0muSt35Fhi0hERETCS4lbd1K6HfqNinQUIiIi4hAlbrGisz5uAKXb1FQqIiLSjSlxizkdJG5luZCYCr0Hhy8cERERCRslbjGjg0nm23i/WSoiIiLdjhK3WNNZUymouVRERKSbUuIWK1wB1LhVFUBLg94sFRER6aaUuMWcDmrcLAvKdujNUhERkW5KiVvMCKDGDaBsu5pKRUREuiklbrGmoz5uYF5Q6DsKXCpaERGR7kZ/3WNFgBVulG6H+CTo428KWBEREYllStxiRidTXrUps4cEUXOpiIhIt6PErbvRWG4iIiLdlhK3WBHIlFcA1YXQVKvETUREpBtS4tYdlW6DrDGRjkJERERCTIlbzAiwjxtAyRYlbiIiIt2QErfuaN8WyBhmJpwXERGRbkOJW6zY38ctgHNLvjZr1bqJiIh0K0rcuqOSr8xaiZuIiEi3osQtZgTRx60sF1qbof9YRyMSERGR8FLi1h15Wkzypho3ERGRbsXpxO1cYAuwDbjNz/FxwKdAI/ATr/1DgXeBzcBG4EavY78BdgNr7WVWaEOOUoGO49Zm39eqcRMREelm4h28txu4HzgbKABWAkuBTV7nlAE3AOf5XNsC3AysAXoDq4E3va79G/AXh+KOcgEmbiVbYOxMcCeYZlMRERGJeU7WuE3G1LTlAk3AEmCuzzl7MQmdb2ZRiEnaAKoxNW/ZjkUaEwKdZd5W8hXExWvOUhERkW7EyRq3bCDf63MBMKUL9xkBHAd87rXvOmA+sApTM1fe0Q3cbjeZmZldeHRgMjIyHLt3m+bevakB0tJSSQjgu7Q0FlINpI44nsSWYsfji0bhKBcJjsokOqlcoo/KJDpFQ7k4WePmr4oowHa+/dKA54EfAVX2vgeAI4BJmJq5v7Zz7UJMYrcqKysryMdGo+Bq3NzluQC09j3SiWBEREQkApyscSvAvGTQJgfYE8T1CZik7SngBa/93tVHjwDL2rn+YXuhuLjYKi/vsFIuJBx9Rno1ADXV1RDQc8qhIo+GtGE0hOG7R7NwlL0ER2USnVQu0UdlEp0iWS5O1ritBEYDI4FE4GLMywmBcAH/xPRtu9vn2GCv7fOBDYcXZqwIso8bmBkUsvRmqYiISHfhZI1bC6Yv2uuYN0wfwwztca19/EFgEKY5Mx3wYJpExwMTgMuB9ZghPwBuB5YDf8I0k1rATuAaB79D9Ah2OBAwc5aOmAquOLA8zsQlIiIiYeNk4gYm0Vrus+9Br+0iTBOqr49ov4rp8hDE1TOUbIGEXpAxHMp3RDoaEREROUyaOSFmBDHlVZu99rB3A44KeTQiIiISfsEmbpOBj4EVwDdDH46E1N7NZj3w6MjGISIiIiERbFPpX4FfYWY8eAQ4MeQRiX9d6ePWVAtlO5S4iYiIdBPBJm6pwNv2dl2IYxEnFG9Q4iYiItJNBJq4/dheD7C3XfT4KajCrQt93ACKN8LYWRCfDC0NIY9KREREwifQPm697eURe50GLHYqKAmhvZsgzg39x0U6EhERETlMgda4rab9GQokHLrSxw1MUynAwPFQuDakIYmIiEh4BVrj9jtHoxDnlO2A5jr1cxMREekGAq1xSwGO49BBcdeENhxpXxf7uFke2PsVDFDiJiIiEusCTdyyMUOBeCduFvCNkEckobd3I4w+J9JRiIiIyGEKNHHbhpK0yOpqHzcwb5Yedzmk9ofaktDGJSIiImETaB+3ckejkCB0MXED9XMTERGJcYEmbt/2OdeN6fcmYePbvTAIbYnbgPGhCUVEREQiItDE7S0OTtR62fsk3LrSVFpXCtVFMOiY0McjIiIiYRNo4pYM1Hh9rkE1buF1GBVuABSth0HHhiQUERERiYxAE7da4HivzycA9aEPRzrXhRo3gMIvof9REJ8U2nBEREQkbAJ9q/RHwLPAHvvzYOAiJwKS9hxmlVvhl+BOMOO57dHweyIiIrEo0MRtJTAOGIvJIL4Cmp0KSvw4nOFA4MB0V4MnKnETERGJUYEmbgnA/wJn2J/fAx5CyVvsqNgF9eUwZJKZeVZERERiTqCJ2wOY5O0f9ufL7X3fdyIo8aeLU155K/zS1LiJiIhITAo0cTsJ8P6L/w7wZejDEUftWQsn/6/p69aqylIREZFYE+hbpa3AEV6fR9n7JFz293E7jHsUfmneKu1/VEhCEhERkfAKNHG7BXgX07ftfUyN280BXHcusAUz1+ltfo6PAz4FGoGfBHhtX+BNYKu9zgzwO4j3CwoiIiIScwJN3N4GRgM32MtYTCLXETdwPzATGA9cYq+9ldn3+0sQ197mFc/b+E8Iu6EQ9HEr3wENleYFBREREYk5gfZxm+/zua3KZnEH10zG1Jbl2p+XAHOBTV7n7LWXbwVx7Vxgmr1/EaYW8KedfwVnpSW5KY90EJ2xLChapxo3ERGRGBXMywkAFwLP2NsWHSdu2UC+1+cCYEqAz+vo2oFAob1dCAzo7GZut5vMTOdaVK+aOoorT+zPzH+soqbJma5/TWm9qQV6904j/jC+S13ZVzROuJyMvlm4rO7dTTEjIyPSIYgPlUl0UrlEH5VJdIqGcgk0cbveXk/12u6Mv6H+A23nO5xr2yy0F7KysoK8NDjri+pIS3Iz59gBPL26sPMLIsi9dwPEJ9PabzTx+76KdDgiIiIShEATtzbBJE8FwFCvzzkcmDLrcK4txky5VWiv97Zzj4ftheLiYqu83LmGzJXbYG1BFvMmDeCBtzd1eXKDDg2oBaC6qgoO57ts/RDOhereo2HrpyEKLro5WfbSNSqT6KRyiT4qk+gUyXIJ9OWEe4F7MAnUPV5LR1ZiXiAYCSQCFwNLA3xeR9cuBRbY2wuAlwK8p6P+s6aIkVmpnDm6v7MPOtyssHS7mUEh+8TQxCMiIiJhE2iN2yp7HcxkSS3AdcDrmLdEHwM2Atfaxx8EBtn3Tgc8mMnsxwNV7VwL8EdMP7urgV3AvCBicsxbW0q5qaqBBaeO4L2vS0L/ANdhTjLvrWAV5ChxExERiTWBJm5J2M2OQVpuL94e9NouwtTiBXotQCkwvQuxOKrFY/HU57u46ewxDO+XQl5pnUNPCkE7bMFKmHYbJKZBU83h309ERETCItCm0ms7P0X+vWIXza0e5p8y3IG7h7DGbfcqcMVB9vGhu6eIiIg4LtDELQO4wM8iXkqqG1m+vpB5Jw4lJdHtzENC8ebD7jVmnXNSx+eJiIhIVAk0cesDzAbmeC2znQoqli36JI/05AQuOC47tDcOZR+3+nLY9zVknxC6e4qIiIjjAu3jtgu4yslAuos1u8pZV1DBlVNH8tSKXSEcGiQEU155K1gFR34zNPcSERGRsAi0xm1j56dIm0c+yOWI/ml886iBkQ6lfQWrIG0AZDjRH09EREScEGji9j3MVFOz7aXTaaZ6suUbisgvq2PhGaNCd9O2ptJQVeHttkd40bAgIiIiMSPQxG0esMJeXwh8DnzXqaBiXavH4p8f7eCkEX05flhGpMPxr3gjNNcpcRMREYkhgSZuv8BMNL8AmA9MBn7pVFDdwTOr8qmoa+J/QlbrFuI+bp4W2PMF5EwOzf1ERETEcYEmbnEcPCdoaRDX9kh1Ta08+dkuzhk/iBH9UiIdjn+7PoPBEyEhSuMTERGRgwSafL2GmX7qCnt5BXjVmZC6j0Wf7KTZ4+HqqSGodQt1HzeAvE/AnaDx3ERERGJEoInbLcBDwARgImb6q1udCqq7KKlp5MUvdjPvxBz6pSZGOpxD5a8ATysMPzXSkYiIiEgAAk3c+gLvAXcCdwDv2/ukEw+/n0uiO46rp448zDuFuI8bQGMVFG9Q4iYiIhIjAk3cCoFVXstqey2dyN1XyyvrC5l/6gj69EqIdDiHyvvYNJW6ozA2EREROUigidsmYJTXMtJeSwDue2cbaUnxXHXaiK7fxIk+bgB5n0JCLxhyXGjvKyIiIiEXzFylc4FzMf3cAp0qS4AtxdW8tqGIK08bSe+kKPvR7frErIepuVRERCTaBZq4vQ98B7gaeADYCcx0KKZu6d53tpLeK4H5p444zDuFuMatdp+ZcF793ERERKJeoInblZiBd+cBpwHTgD87FFO3tHFPFW9vLubqqSNJSXR3/UahbioFMyzIsCng0tB8IiIi0ayrf6m3AWeHMpCe4N53ttE3NZHvndyFid3b+rg5Ie8TSM6AgUc79wwRERE5bJ11uLqnk+M3hCqQnmBtfgUffF3CNWeM4qnP8qhtau3CXRyocdv5kVmPOB2K1of+/iIiIhISndW4zcUM/dHeIkH6yxtb6JeWxNWnBzuum4M1blW7Yd9WGDXNuWeIiIjIYeusxq0MWBSOQHqKdQWVvLahiP85fRSLP82joq45uBs40ccNIPc9mHSJGc+tNciYREREJCw6q3FzKEvo2f76xhZSE+O59swjAr/IyT5uYBK3xDTIPtHZ54iIiEiXOf0a4bnAFszLDLf5Oe7C9KPbBqwDjrf3jwXWei1VwI/sY78BdnsdmxXqoJ22dW8NL67dzRWnjmBA76QAr3I4cdv5kZm3VM2lIiIiUauzxG0iJmnyXartdUfcwP2Y8d7GA5fYa28zgdH2shAzRhyYZG+SvZwA1AEvel33N6/jyzuJIyr97c2viXO5uP4boyMditFQAXu+gFFnRjoSERERaUdniZsbSPez9LbXHZmMqUnLBZqAJZiXHbzNBRZjmmQ/AzKAwT7nTAe2A3mdPC+mFJTXs2TlLi6ePJRhfVM6v8CpKa+85b5n5i1N6u3cM0RERKTLnJx/KRvI9/pcAEwJ4JxszKT2bS4Gnva57jrMgMCrgJuB8o4CcbvdZGZmBhx4sDIyMrp03ROrS5h3wlB++e1jufWlrzs8tzEllTogvU86bpz5Ls0lq6mJiyf16Bkk7njHkWeEU1fLRZyjMolOKpfoozKJTtFQLk72cfPXKcu3uqizcxKBbwPPeu17ADgC00xaCPy1necvxCR2q7KysgIIN/z21Tbz+Oe7OXtcFpOyO6nlCkONW3zhGmiup2XoVMeeISIiIl3nZI1bATDU63MOsCfIc2YCa4Bir33e248Ay9p5/sP2QnFxsVVe3mGlXEh05Rn3vFHFeRP6c+OZQzn/Hx+3n5fV1gJQVVUFTn6XvE9ozDmVxjD8vMIlHGUvwVGZRCeVS/RRmUSnSJaLkzVuKzEvHYzE1JxdDCz1OWcppsnTBZwMVHJwM+klHNpM6t0H7nxgQ+hCDr/65lb+/PoWJg3N4NsTh7R/Yjj6uAFsexP6j4OMYc4+R0RERILmZOLWgumL9jqwGXgG2Ahcay9g3gjNxbzE8AjwA6/rUzDzob7gc98/Aesxw4ecBdzkTPjh8+IXu1lfUMlPzx1HckKEJ3r/+nWzHj0jsnGIiIjIIZxsKgWTmPkO1/Gg17YF/LCda+uAfn72Xx6CuKKKZcGdr2ziP9ecwvdPH8V972zzc1Zbd0CHa9zKcqF0G4w5B1Y+6uyzREREJCgRrt6RNp/vKOO1DUX8YNoRDO6THNlgvn4dRp4BCQEMUyIiIiJho8Qtitz5yibiXC5+8S3fcYoJXx83gK1vQHyySd5EREQkaihxiyIF5fXc9+42vjVhMKePjuAQJnkfQ2O1aS4VERGRqKHELco88kEuO/bV8ttvH02i27t4wtTHDaC12cyioBcUREREoooStyjT2OLhN0s3Mqp/Gt8/feShJ4SjqRRMP7c+OTDw6PA8T0RERDqlxC0Kvf91Ca9uKOT6b4wmJ7OX2enyN8mEg7a+AZYHxs4K73NFRESkXUrcotQdL2/CY1nced4xPkfCVONWUwz5K2D8t8PzPBEREemUErcotaeygT+/voVpYwdw3qRs/E/r6rBNL8GgCZDpp8lWREREwk6JWxRb/OlO1uSV86s540lNcpud4erjBrD5ZbNWrZuIiEhUUOIWxTwW3Pr8OlKT3MyekB3+ACrzYfcaOEqJm4iISDRQ4hbltu2t4f53tzFhaIa9J4w1bgCbl0LOiZAegcRRREREDqLELQY88N52iqoaAEhLdnp6WR+bXjLro+aE97kiIiJyCCVuMaC51eL5NbsBuHXGmPA+vCwXijeon5uIiEgUUOIWI/ZUmBq3c48ZzJwJg8P78I3/hWGnqLlUREQkwpS4xQwzHMj6ggruPP9YBqYnhe/R658DVxwc+93wPVNEREQOocQtxvx66QYS3C7+/N2J4ZtMoXwH5H8OEy4M0wNFRETEHyVuscLO0grK6rhz2WbOGNOfK04dEb7nr3sGBh6juUtFREQiSIlbDPr3il28uamYn808imOz+4TnoRtfgNZmmHBReJ4nIiIih1DiFnPMOG63PPclJTWN3HfpcfROCsMQIXVlsO1N088t3BPei4iICKDELWZV1DVzw9NfkJ3Riz9ccGx4HrruGfNm6YjTw/M8EREROYgSt1jRVsvlNVfp6rxy/vrG18yeOIRLJw9zPoYtr0JDBUy6zPlniYiIyCGUuMW4Bz/Yzvtb9vLrOeOd7+/W0gDrnoWjz4Nemc4+S0RERA7hdOJ2LrAF2Abc5ue4C7jHPr4OON7r2E5gPbAWWOW1vy/wJrDVXveQDKKtX9nBc5VaFtz0jOnv9tDlJ9AvNdHZMNYsgvhkDQ0iIiISAU4mbm7gfmAmMB64xF57mwmMtpeFwAM+x88CJgEneu27DXjbvuZt/CeE3Zd16CTzZbVNXPPEavqmJvKPy44nPs7BlweK1sPu1XD8AueeISIiIn45mbhNxtSk5QJNwBJgrs85c4HFmGqkz4AMoLP5nOYCi+ztRcB5IYk22nXyJufGPVX89Pl1TBnVj1/M9s2PQ2z1IjOeW85Jzj5HREREDuLkOBLZQL7X5wJgSgDnZAOFmGTuDXv9EPCwfc5A+zj2ekBngbjdbjIznWtRzcjIcOzebRp6pVAP9OnTh7hk/+d8kFfPEyv2cMWpI8gtb+blDSWOxGIVvE1FUw2Jp15D6pvbHHlGKISjXCQ4KpPopHKJPiqT6BQN5eJkjZu/KiLfdr6OzjkN0+dtJvBD4Iwgn78Q0zduVVZWVpCXRqEAx077+3s7+XxnBT8/5wgmZvd2JpTmWhK3LKVp9Gw8SWEaAFhEREQcrXErAIZ6fc4B9gRxTtt6L/Aipun1A6AY05xaaK/3tvP8h+2F4uJiq7y8vEtfIhiOPqO+HoDKigozJEcHrlm8ghf+91TuPn8sF/zjY3aW1oU+no/uh2MvpfKIufDx30N//xAKR9lLcFQm0UnlEn1UJtEpkuXiZI3bSswLBCOBROBiYKnPOUuB+Ziat5OBSkxClgq0VRelAjOADV7XtPWMXwC85Ez4sauirpkrH1+JZVn868rJZKYkhP4hxRtgxwcweSHEhWHmBhEREXE0cWsBrgNeBzYDzwAbgWvtBWA55uWFbcAjwA/s/QOBj4AvgRXAK8Br9rE/AmdjhgM52/7cA/gfDqQ9eaV1/M/iVQzpk8yjC04iKd6Bov70fuiTA0fNCf29RURE5BBOV5UstxdvD3ptW5j+a75ygYnt3LMUmH74oXV/a3ZV8KP/rOX+S4/n7gsncd3Ta/yNJtJ1W1+Hslw4+X9h44shvLGIiIj4o5kTYoWfKa8C8eqGIn7/6ma+NWEwv/v2MaGNybLgswdg6BTIPiG09xYREZFDKHHrAR79cAcPvredy08Zzq3njA3tzdf+27wscdqNob2viIiIHEK9ymNGcH3cfP3xta9IS47nB2cdSU1jC/94b3towmqqgc8fhjN+Av3HQsmW0NxXREREDqEatx7kly9t4L9f7ObWc8dx+cnDQ3fjzx+A5nqY+uPQ3VNEREQOocQtVnSxj5s3y4KfPPslb24q4o7zjuGik4Z2flEg6spg1WNw7DzIHBmae4qIiMghlLj1MC0ei+v+/QXvbdnL/31nApdNGRaaG39yL3iaYepNobmfiIiIHEKJW8w4vD5u3hpbPCxcvJq3Nhdz1/nHMv+UEDSb1hTDmsUw6RLICGEzrIiIiOynxK2Hamr18L9Prub1jUX8bu4xXD01BE2cH94NnlY46/bDv5eIiIgcQolbrAhBHzdfza0WP3xqDa+sK+SXs8dzw/QjD++G1YXw+UMw4UIYeHRoghQREZH9lLj1cC0eixuWfMFzq/P58dljuWPuMcS5Or+uXR/9DRqqYPqvQxajiIiIGErcYkbo+rj5avVY/OTZdTzw3jYuP2U49116fNfnNm2ogI/uhjHnwPBTQxqniIhIT6fELdaEdLLRg/3fa1u4Y9kmZh07mEVXTSY9uYvjM3/+EFTthnN+Dy79ExMREQkV/VWNFa7Dab8M3D8/2sGNS77ghOGZvPCD0xjeLyX4m7Q0wBu/hCHHwXGXhz5IERGRHkqJW8xxrsatzUtr93D5Pz+nX2oi//3BaZw8qm/wN9nwPOz8CL75a+iVGfogRUREeiAlbjEjPDVubT7LLWPu/R9TWtvEE1dP6dosC6/+FJIzNDyIiIhIiChxixUODAfSmV1ldZx//8d8sm0f//edCfxu7tEkuoP4J1O8AVY+CideDYMnORaniIhIT6HETTpU3djCVYtW8fAHucw/ZQTPXHsK2Rm9Ar/Bu3eZWRXm3gfuBOcCFRER6QGUuMUM54YD6Uyrx+L3yzdzzROrGdU/lWXXT2XamP6BXdxQCa/8GAYdq3lMRUREDpMSNwnY6xuLmHPvRxRW1vP4VZO5ecYY3IGM1rvlVVj/LJxxCww4yvlARUREuiklbrFif4Vb+GvcvOWV1nH+Pz7hPyvzuf4bo3nu2lMYEciQIa/+1MyocP5D4E50PlAREZFuSImbBK2xxcNPn1/Hdf9ew6j+aSy/8XQu7uyt07pSWHo9DJ4I038VnkBFRES6GSVuMSNyfdzas2xdIef87QPW5FXwx+9M4JH5J9AvtYPatC3LzVump14PR0wPX6AiIiLdhNOJ27nAFmAbcJuf4y7gHvv4OuB4e/9Q4F1gM7ARuNHrmt8Au4G19jIr5FFLwIqqGrj8sc+5Y9kmzhjdn7d+fCYXHJ/d/gWv/xz2boLzH4C0AeELVEREpBtwMnFzA/cDM4HxwCX22ttMYLS9LAQesPe3ADcDRwEnAz/0ufZvwCR7We5E8FEnAuO4BcqyzFRZ37r3I7aX1HD3hZNYfNVkcjL9DBvS0gDPXQVJvWHeIg0RIiIiEgQnE7fJmJq0XKAJWALM9TlnLrAY0/73GZABDAYKgTX2OdWYmrcOqnEkGmzbW8O8hz7lF//dwPHDM3njpjO4eurIQ9883bsZXroOhp9qJqIXERGRgMQ7eO9sIN/rcwEwJYBzsjGJW5sRwHHA5177rgPmA6swNXPlHQXidrvJzHRuvsyMjAzH7t2mvlcKDfazXFHUz82fV7ZUsWrPWn42YxS/nD2eS6aM4E9v5bJyV9WBk3a/Q93qh2icfA0plVtJ2vRsyOMIR7lIcFQm0UnlEn1UJtEpGsrFyRo3fwN8+WYcnZ2TBjwP/Aho+6v/AHAEppm0EPhrO89fiEnsVmVlZQUUsIROcXUTP3r+K256fjPJCXE8fMkx/Pm8sQxOT9p/Tq+P/0T8rg+p+8bvaR56agSjFRERiQ1O1rgVYF4yaJMD7AninARM0vYU8ILXOcVe248Ay9p5/sP2QnFxsVVe3mGlXEg4+oyGBgAqysuce4YDXlxZzvIvdvL900fxw7OOYOr3J/HwB7k88kEu1Y0t8NSlcNWr1Mx6AP41E4o3hjyGcJS9BEdlEp1ULtFHZRKdIlkuTta4rcS8dDASSAQuBpb6nLMU0+TpwryEUImpRXMB/8T0bbvb55rBXtvnAxtCHbiEVmOLh/vf3cb0v77PGxuLuGH6aN6/9SyunjqSpNYaeGoeNNXAZc9Cn07GgxMREenBnEzcWjB90V7HJGDPYIb2uNZewLwRmot5ieER4Af2/tOAy4FvcOiwH38C1mOGDzkL6CETYAYwtVSUK6xs4IYla5l974ds2F3JL2eP552fTGPe6DhcT30XElLhimWQrvdQRERE/HGyqRRMYuY7XMeDXtsWZqgPXx/RfqZyeQjikgjasLuK+Y+t4JQj+vHTc8fx53kTuXZvDXeuXMi7kx8xydvj34Iq35Z1ERGRnk0zJ8QKV+zXuPn6dHsp593/Mdc8sZrGllb+NSOeB7mL+PQBuBa8DOlDIh2iiIhIVFHiFjNcYHkiHYQjXt9YxKx7PuKqx1cyoHI9/0n5M736DiH1B++QMsR3zGYREZGeS4mbRI13vtrLBQ98wl/++QS/KLmZXsm9iF/4JlddchEj+qVEOjwREZGIc7qPm4SKyxWV01054dPtpXz6j38zeuw6yi9Ywr/H3stfjz6K5K3LePyTnXy0bV9P+VGIiIgcRImbRK2tWzbAPdOIv2wJPxxyE1cMG8bj415gd1kt/1mVz3Or8ymuaox0mCIiImGjptKY4eLQiSd6gNp9tPxzFqx+nMfd3+HkXdeyucLNLeeM5ZPbpvPoghM5e/xA4n3nQxUREemGVOMm0a+1CV6+EfZ8QcmsP3NN7T0M+PeNzB+0i3kn5vDNowZSWtPIK+sLeWntHtbsKldTqoiIdEtK3GJFD+rj1q7Vj8OetfCdR9h7wXP85ZN7ufvPd/CNI/swd1I2804YyvxTRlBQXsdLa/fwXm4N2/bVRTpqERGRkFHiJrGlcC08dAbMuAtOvR7P6Bm89fKNvPX0p6Qmujl7/CDmThrCNWeM4odnxbGrrJ5X1+/hjU3FfLGrHE8Pz31FRCS2KXGLGT20j5s/zfXwyo9hyyvwrb/BVa/BmsXUvvkr/rt2N/9du5u+qYlccNIozhrTlytPG8k1Zx5BSXUDb27ay5ubivk0dx8Nzd1zXDwREem+lLhJ7Nr2NvzjZJh2G5zyQzjq2/DhX2HFQ5TVNvL8l8U8/2UxLXXVTBvbnxlHD2LOxMFcOmUYjc2trNhZxodb9/HB1yV8VVQd6W8jIiLSKSVusUJ93PxrroM3fwVfPg1n/w5m3AGTvw/v3IlV8DYuy0N1Ywsvryvk5XWFJLrjmDKqL2eM6c8Zo/tz+6yjuH3WUeytauDDrfv4cOs+VuwoZU9lQ6S/mYiIyCGUuEn3sHczPDUPRp4JM+6ECx6hqjyX5FUPwGf/gtZmAJpaPfsTtLvYzKD0ZE4fncUZY/rzjXED+M4JOQAUlNfx+Y4yVtjLjn21kfx2IiIigBK3GKI+bgHZ8T48fAYc9W1c026l7uw/w+Qb4LMH4IunoKHioNOLqhp4dnUBz64uIM4F4walM3lkX04a0ZczRvfnO8ebRK6kuoGVO8v5Mr+CLwsqWF9QSW1Ta/i/n4iI9GhK3KT7sSzY9BK9Cz+gZfg0ao5bCOf8Hr7xS9jwHKx8DPasOeQyjwWbCqvYVFjF45/sBGBUViqTR/Zl8si+nDA8k1nHDjbneiy2ldSYRC6/gnW7K9lSVE1ji154EBER5yhxixXq4xY0F5CQ9x6sfREGHgMnXQ0TLoTjLoeidbDuWdjwPFTtbvceuftqyd1Xy5KV+QBkpiQwISeDiUP7MDEng7PGDWDeiUMBaPVY7NhXy1eFVWwuqmZzYRVfFVapv5yIiISMEreYosSty4o3wLKb4M1fm+Rt4sXmRYYZd0Dex7DxRdjyGlTmd3ib8rpm3v+6hPe/Ltm/LyezF8dk9+GoQb05anA6E3IymD1xyP7jVfXNfFVUzfaSGraX1JBbUktuSQ355fW0amA5EREJghK3mKG5OEOisQpWPmqWvqPgmAvg2Ath1l/MUrwRvn4dvn4Ndq8GT0untywor6egvJ7XNhTt35eWFM+Ygb0ZP7g34wanM3ZQb2aMH0i/tGH7z2lq8ZBXamr0cktq2LGvjvzyOvLL6iisbFBSJyIih1DiFkvUVBpaZbnwwV/M0u8IGHOuWU69Hk7/MTTVwK7PTY3czg9hzxf7307tTE1jC2t2lbNmV/lB+/v0SuCI/qkc0T+NUf1TGdU/jVFZqZw1dgCJ8XH7z2tp9VBY2UBBef3+ZC6/rJ6C8jqKqhrYW9VIU6v604mI9DRK3GKFSzVujirdDp/eb5akdDjiLBgxFYafBtN/Zc5paYCi9SaB273GrPd9DVbgCVRlfTNrdlWwZlfFQfvdcS4G90lmaN8UhmamMLRvL3udwrQx/RmQnnzIvfbVNFJc1UBRpb1U2UtlA8VVDeyraaKirknTfImIdCNK3GKK/gKHRWMVbHrJLAApfU0CN3QKDDkOJl4CkxeaY811UPI1lHx1YNn7FVTkBZXQtXqs/U2un1J6yPGk+DhyMlPIyezFwPRkBvVJYlB6sr2dzKShGfRLS/J737LaJkprGymtaaK0ppHS2ib2eW2X2gleZX0zlfXNtCjTExGJWkrcRDpTVwabXzYLgCsOskabJG7QBOg/ztTOTbz4wDWtzeZFh/I8k8SV77S3d0J1EdTsDaj/XJvGFs/+lxvak+iOY0C6SegG9Ummb2oi/dKSyLLX/dISOSa7D1m9k0hPTmj3PjWNLVTUNVFV30yFncxV1Jl1Zdu6vpnaxhZqfJbaxhaaW5X4iYg4xenE7Vzg74AbeBT4o89xl318FlAHXAGs6eTavsB/gBHATuBC4OCORN2RhgOJHpYHSraY5cslB/YnpUP/MdD/KOg7EjKGQ+YIGDcbUrMOvUftPjuJK4LqYqguhNoSqC+D+nKoKzfr+nJorOy0/JtaPftr7TqTFB9nJ3aJZKUmkd4rgYyUBPp4rfv0SrT75KWZ/b0SSEpwd3rvxpZWahpaqG1spbapheqGlv1JXjNuGppbqageQH1zK/XNrTQ0tdLQ4qG+qfWgfW3b9U2tNNjbSgpFpKdzMnFzA/cDZwMFwEpgKbDJ65yZwGh7mQI8YK87uvY24G1MInebvfzUwe8hEpjGKihYZRZfiWmQORwyhkHaIOg9EHoPhrSB0HuQGWcubSDEtZMYeVrNrA/15dBYbZamWvMCRWONWXtvN1ZDcz20NJq+eS0NXtuNNLY0UNjQQGFNPbRWBfwVk+Lj7KQugbSkeFKT4klLjictKf7A5yQ3aUkJpCa56Z1s9vVLS2RYvxTSkxNJTogjKT7uoJcxAtXS6qGxxUNTi4emVnvtve1vX3vHvT63eiyaPWbd0mrR4vHYa4uWVo9Zex1r9Vg0t7atrYM/ezy0tl3rsfR2sIiElJOJ22RgG5Brf14CzOXgxG0usBjTeeszIAMYjKlNa+/aucA0e/8i4D06SdzcbjeZmZld/yadSB16DE0ZR5Ba49x8lk0DxtIMjn6P7iYjIyPSIRysaQ/s3QN7/R+2XHFYSX2wkjOxkjOwkvvg8dq2kjPxJGdgJaZBYipWaj+sxBSshDSsxFRISOlaXJYHWhpxtTaCpxk8rbg8LaYp19OCy96Hp4UWTwtlnhbK9p934Nj+bctjllYP1Fq4ajz79yUlJoLlobGhnjg8uF0Q7/IQHwcJLuvgdZxFQhz2YpG4f58LdxzEu8xLHfFx9joB3IkuEtwQ73IdfMwF8W5IjHOREgdul2v/ADsur76jbdsul3XI8QNr72sO3ed9Xttny7LwWBaWhb1trmzbd2ANFgf2eSzzCA8dn2fu29555pnYx9rCbuuBaVkW8fHxWBY0t7RgtX03+3oOXLI//v2fve5p1tb+89r2HXikCXL/8bbY2o5aXtcdiMLn3+pBK6/n+j3t4C3LLotOrsXr+/nepaOYDjrPX0xWh3c55NqkJNNntbGxsZ3veOiDA/tfBKvDj05pp0Q7uaYLz3H4+yQlJ1Ox40vI3ezsgzrgZOKWDXiPZlqAqU3r7JzsTq4dCBTa24XAgHaev9BeyMrKaueU0KjLOZOyybc4+gwAV3X7I/xL7HNZHlwN5dDQtZZ/yxUHCSlYCalYiWlY8cngTsKKT4L4JCz3gbUVn3TgmLvtWDJWXDzExUNcAsS5vT7HH3TMiu+1f5/ldT5x8VguFxBn+gK64kwzv73dHOfGwoWFy65dPHCs3drGULKAaJpi1oWGaBSJMdPcv2ZXN03c/P068s2F2zsnkGs787C9UFxcbJWXO9cNrs/aJfTe9h7V1dWOPQPAqtpNeW33784Xak6WffQ59I3UaJJh1xh3WCauuE4Wl9fwOPba+7PfY37O6+hYu/cP4h6+xyKpk+GEeqenA1BdFXizeReCOORjnF3zacI7sO19pm/oLq8drnbOOXC3A8ctl8vnvq5DojpQeq5Dd/qJyV+cLj8nunxP8om57Xne16b1TgOg1qslp8PvG8A/M8vny3TlX2ZXRqZydeFJXRoBKwyxpaamsG1PbkT/rjiZuBUAQ70+5wB7AjwnsYNrizHNqYX2up2Gp/CJqy8lrr4UelSCIOKQtqZWCZv4JrsLRph/h6mU25cZyP/kSNhFQ3el4HsHB24l5qWDkZhE7GLMCwbelgLzMXnyyUAlJiHr6NqlwAJ7ewHwkmPfQERERCSKOFnj1gJcB7yOeUv0MWAjcK19/EFgOWYokG2Y4UCu7ORaMG+TPgNcDewC5jn4HURERESihtPjuC23F28Pem1bwA+DuBZMJ57phx+aiIiISGxxsqlUREREREJIiZuIiIhIjFDiJiIiIhIjlLiJiIiIxAglbiIiIiIxQombiIiISIxQ4iYiIiISI1yWFewUoDGpBMhz+BlZwD6HnyHBU7lEH5VJdFK5RB+VSXQKR7kMB/r7O9BTErdwWAWcGOkg5BAql+ijMolOKpfoozKJThEtFzWVioiIiMQIJW4iIiIiMUKJW+g8HOkAxC+VS/RRmUQnlUv0UZlEp4iWi/q4iYiIiMQI1biJiIiIxAglbqFxLrAF2AbcFuFYeqqhwLvAZmAjcKO9vy/wJrDVXmdGJLqezQ18ASyzP6tMIi8DeA74CvPfzCmoXCLtJszvrg3A00AyKpNIeAzYiymHNh2Vw88wf/u3AOeEI0AlbofPDdwPzATGA5fYawmvFuBm4CjgZOCHmHK4DXgbGG2vlViH342Y5KCNyiTy/g68BowDJmLKR+USOdnADZghJo7B/F25GJVJJDyOqYzx1l45jMeU09H2Nf/AlJ2jlLgdvsmYbDsXaAKWAHMjGlHPVAissberMX+IsjFlscjevwg4L+yR9Ww5wLeAR732qUwiKx04A/in/bkJqEDlEmnxQC97nQLsQWUSCR8AZT772iuHuZi/+Y3ADkwuMNnpAJW4Hb5sIN/rc4G9TyJnBHAc8DkwEJPUYa8HRCimnur/AbcCHq99KpPIGoWZTeZfmCbsR4FUVC6RtBv4C7AL87OvBN5AZRIt2iuHiPz9V+J2+Fx+9ulV3chJA54HfgRURTaUHm82pq/I6kgHIgeJB44HHsD8D04taoKLtExM7c1IYAgmkf5eRCOSQETk778St8NXgOkY3yYHU8Ut4ZeASdqeAl6w9xUDg+3twZhEQsLjNODbwE5Mc8I3gCdRmURagb18bn9+DpPIqVwi55uYprYSoBnz++tUVCbRor1yiMjffyVuh28lpsPiSCAR01FxaUQj6plcmD47m4G7vfYvBRbY2wuAl8IcV0/2M8wvshGY/y7ewdQiqEwiqwjTvDPW/jwd2ITKJZJ2YV6qSsH8LpuO+V2mMokO7ZXDUszvtiRMDjAaWOF0MBqANzRmYfryuDGvEt8V0Wh6pqnAh8B6DvSnuh1Tq/AMMAzzy3Eeh3Y8FedNA36CaT7th8ok0iZh+rYlYl6suhLzP/Iql8j5LXAR5g35L4DvY7p+qEzC62nM76ssTE3br4H/0n45/By4ClNuPwJedTpAJW4iIiIiMUJNpSIiIiIxQombiIiISIxQ4iYiIiISI5S4iYiIiMQIJW4iIiIiMUKJm4j0BK3AWuBLzJy2p0Y0GhGRLtJwICLSE9RgxsQCOAczxt+ZkQtHRKRrVOMmIj1NOlBub08Dlvkcz8JM0wVwBWaU9NeALZjBOAHuAG70uuYu4Ab7/Pu89tf4ec6ZmIGh+2CSybcxtYDrMfNVioi0Kz7SAYiIhEEvTFNpMmauwW8Ece1k4BigDjPF3SuY6dVeAP6O+R/gi+3zvoX/iafbHGtfMwuoxPwOPh+owiSMn2Gm0VFTiIj4pRo3EekJ6jHTPI0DzgUWcyDBOh2T1H2BmbrG15tAqX2PFzDTq+209x0HzLCvLcVMOj0R/79bh2Cmw1nEgYmoXcDvgXXAW0A2MLBrX1FEegLVuIlIT/Mppnarv/35Q8wcqlnAV5gEyptv7Vfb50cxTaODMHMUA7yHaVJdDzT7XDcOUzP3J+BJoAS4zI7jBPv8nZhaQRERv5S4iUhPMw5wY2rIvFVjJop2++w/G+iLqXE7jwO1ci8CvwMSgEvtfR7gaq9ra7y238E0g2ZgmksvxfRz24tJ2s4ChnfpG4lIj6HETUR6grY+bmCaJxdghggBMzTIR0Aq8DdMAuftI+AJ4Ejg38Aqe38T8C5Q4XWvQCzG1LTNAp4CXrbvuRZT4yci0i4NByIi0r4rgBOB6/wci8O8DToP2BrGmESkB9PLCSIiwRsPbMMM5aGkTUTCRjVuIiIiIjFCNW4iIiIiMUKJm4iIiEiMUOImIiIiEiOUuImIiIjECCVuIiIiIjFCiZuIiIhIjPj/w2Totyr53fgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "dark"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_dist = expon(scale = 10)\n",
    "control_dist = expon(loc=5, scale = 5)\n",
    "\n",
    "x = numpy.linspace(0, 100, 1000)\n",
    "\n",
    "pyplot.figure(figsize=(10, 5))\n",
    "pyplot.title('Пример распределений', fontsize=12)\n",
    "pyplot.plot(x, test_dist.pdf(x), label='test')\n",
    "pyplot.plot(x, control_dist.pdf(x), label='control')\n",
    "pyplot.xlabel('Выручка')\n",
    "pyplot.ylabel('Плотность')\n",
    "\n",
    "pyplot.legend(fontsize=12)\n",
    "pyplot.grid(linewidth=0.2)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a6c6226-aa1f-4506-8c3a-c5b541e80006",
   "metadata": {},
   "source": [
    "Например, раньше в среднем выручка от пользователя была примерно 10 руб и она была не меньше 5 рублей, а после введения эффекта часть пользователей стала меньше платить, но средний чек остался таким же: 10 руб."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "56d48d7d-ced4-4b13-9375-7f4ad4590975",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_criterion(test_dist, control_dist, sample_size, N_exps=10000, to_print=True):\n",
    "    \"\"\"\n",
    "        Функция для проверки t-test критерия для AB-теста\n",
    "        Возвращает доверительный интервал для FPR, если флаг to_print = False. Иначе печатает результат.\n",
    "    \n",
    "        Параметры:\n",
    "            - test_dist: Распределение тестовой выборки в эксперименте\n",
    "            - control_dist: Распределение контрольной выборки в эксперименте\n",
    "            - sample_size: размер выборки теста и контроля\n",
    "            - N_exps: число экспериментов, по которым потом считается FPR\n",
    "            - to_print: печатать результат или нет. Если нет, то функция возвращает дов. интервал для FPR.\n",
    "    \"\"\"\n",
    "    \n",
    "    numpy.random.seed(35)\n",
    "    bad_cnt=0\n",
    "    alpha=0.05\n",
    "\n",
    "    for i in range(N_exps):\n",
    "        # Генерирую выборку\n",
    "        test    = test_dist.rvs(sample_size)\n",
    "        control = control_dist.rvs(sample_size)\n",
    "\n",
    "        # Запускаю критерий и считаю p-value\n",
    "        pvalue = ttest_ind(test, control, equal_var=False, alternative='two-sided').pvalue\n",
    "\n",
    "        # Проверяю, что pvalue < alpha\n",
    "        bad_cnt += (pvalue < alpha)\n",
    "\n",
    "    if to_print:\n",
    "        print(f\"FPR: {round(bad_cnt / N_exps, 4)}\")\n",
    "        print(f\"CI={proportion_confint(count = bad_cnt, nobs = N_exps, alpha=0.05, method='wilson')}\")\n",
    "    else:\n",
    "        return proportion_confint(count = bad_cnt, nobs = N_exps, alpha=0.05, method='wilson')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "ecfd1be8-5ea2-492f-8996-bbd029ca1645",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FPR: 0.0443\n",
      "CI=(0.04043912932872393, 0.04851084678631071)\n"
     ]
    }
   ],
   "source": [
    "check_criterion(test_dist=test_dist, control_dist=control_dist, sample_size=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78390350-13e2-4080-9443-8b5f54f27865",
   "metadata": {},
   "source": [
    "Что ж, мы видим, что t-test не сработал. Истинное $\\alpha$ не лежит в доверительном интервале. Но с какого размера выборк t-test начнет работать правильно?\n",
    "\n",
    "### Как проверить минимальный размер выборки при котором t-test работает?\n",
    "\n",
    "Для этого давайте просимулируем эксперимент с разным размером выборок и посмотрим, с какого размеры выборки у нас $\\alpha$% лежит в дов. интервале."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f240f758-29d8-4203-a7ae-fe196b77744a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min sample size: 40\n"
     ]
    }
   ],
   "source": [
    "scale = numpy.arange(20, 110, 20)\n",
    "for N in scale:\n",
    "    left, right = check_criterion(test_dist=test_dist, control_dist=control_dist, sample_size=N, N_exps=10000, to_print=False)\n",
    "    if left < alpha < right:\n",
    "        print(f\"Min sample size: {N}\")\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "aec8451a-9cf6-4964-8832-c0060c32ab52",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FPR: 0.0488\n",
      "CI=(0.04474778133939989, 0.05319873879029866)\n"
     ]
    }
   ],
   "source": [
    "check_criterion(test_dist=test_dist, control_dist=control_dist, sample_size=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d83355b-58ab-4ba8-8d15-5bd883f53126",
   "metadata": {},
   "source": [
    "Так мы видим, что для выборки размера 60 &mdash; 5% уже попали в дов. интервал. Так что мы можем считать, что для таких распределений нам достаточно выборки размера 60, чтобы использовать t-test. \n",
    "\n",
    "Но надо понимать, что реальный FPR может быть не 5%: он лежит в доверительном интервале (0.045, 0.053). Если нужна большая точность &mdash; нужно провести больше экспериментов (`N_exps -> infinity`)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3bde783-98a5-4462-9940-5f304036046c",
   "metadata": {},
   "source": [
    "\n",
    "## Итого\n",
    "\n",
    "Чтобы проверить критерий, надо уметь много раз проводить один и тот же эксперимент.\n",
    "- Правильно ли реализорван критерий?\n",
    "    - Проверьте его! Можно на специально смоделированных данных.\n",
    "- Можно ли использовать данный критерий для нашей задачи?\n",
    "    - Проверьте его! Но только нужно **правильно** сгенерировать эксперимент.\n",
    "- Как найти минимальный размер выборки у t-test?\n",
    "    - Прверьте t-test на разных размерах выборки. С того момента, как $\\alpha$% лежит в доверительном интервале &mdash; можем считать, что t-test будет работать.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af342260-70fa-4258-af45-a6ce3cfbaf42",
   "metadata": {},
   "source": [
    "# Как смоделировать эксперимент?\n",
    "\n",
    "Есть 2 ответа на этот вопрос:\n",
    "1. Генерация теста и контроля через искусственное моделирование. С помощью разных распределений можно попытаться приблизить реальное распределение на данных. Например:\n",
    "    - Для генерации выручки использовать экспоненциальное распределение. Чем больше выручка от пользователя &mdash; тем меньше таких людей.\n",
    "    - Для генерации конверсионных выборок (например, кликнет/не клинкет) использовать бернуллиевскую выборку.\n",
    "    - Иногда можно брать смесь распределений: пусть 90% пользователй нашего сайта приносят нулевую выручку. Тогда можно перемножить бернуллиевское распределение на экспоненциальное для моделирования выручки от пользователя.\n",
    "    - Также для проверки криетрия не обязательно распределения в тесте и в контроле должны совпадать. Для проверки критерия равенства средних не обязательно должны совпадать распределения в тесте и в контроле. Они могут быть разными, но мат. ожидание совпадет, как было в примере выше.\n",
    "    \n",
    "---     \n",
    "    \n",
    "2.  Датасеты на исторических данных компании. У многих компаний есть логирование событий. Тогда мы сможем прямо на реальных данных оценить работоспособность критерия! И не попасться в ловушку того, что на искуственных выборках критерий валиден, а на реальных данных нет. Например, у нас есть данные о транзакциях пользователей за несколько лет. Это уже один готовый датасет: вы делите всех пользователей на тест и контроль и получаете один «эксперимент» для проверки вашего критерия. \n",
    "\n",
    "Осталось понять, как из одного большого датасета сделать N маленьких датасетов. Я расскажу, как мы это делаем в Авито, но описанная механика применима практически к любой компании.\n",
    "\n",
    "Наши пользователи размещают объявления. Каждое объявление относится только к одной категории товаров и размещено только в одном регионе. Отсюда возникает незамысловатый алгоритм:\n",
    "\n",
    " - Разобьём все размещения пользователей на четыре (или N в общем случае) категории: автомобили, спецтехника, услуги и недвижимость. Теперь наш датасет можно разбить на эти подкатегории: к примеру, в одном датасете смотреть выручку пользователя только в этой подкатегории.\n",
    "\n",
    " - Поделим датасеты по месяцам: датасет трат пользователя за ноябрь, за декабрь и так далее.\n",
    "\n",
    " - Ещё все метрики можно поделить по субъектам РФ или по группе субъектов: датасет трат пользователя из Москвы, датасет трат пользователя из Хабаровска и так далее.\n",
    " \n",
    " - Объедним все 3 правила в одно. Например: датасет трат пользователя в Авто за ноябрь в Москве.\n",
    "\n",
    " - Теперь у нас есть большое число датасетов и в каждом из них есть пользователи. Поделим пользователей случайно на тест и контроль и получим финальные датасеты для валидации придуманных статистических критериев.\n",
    "\n",
    "Давайте посмотрим на картинках, как такая схема увеличивает количество датасетов:\n",
    "    <img src=\"https://habrastorage.org/getpro/habr/upload_files/71a/497/8cb/71a4978cbbb1f177c5edb360456f3e90.png\" width=\"1500\" height=\"200\" />\n",
    "\n",
    "\n",
    "Здесь мы смогли разбить 1 датасет на 16 датасетов. А если добавить ещё и разделение по субъектам РФ, которых больше 80, то мы получим уже 16×80 = 1280 датасетов для проверки. И это всего за 5 месяцев! При этом, как показывает практика, 1000 датасетов достаточно, чтобы отделить некорректный критерий от хорошего.\n",
    "\n",
    "**Сравним 2 метода**\n",
    "\n",
    "Главные плюсы искусственных данных в том, что их сколько угодно, они генерируются быстро, и вы полностью контролируете распределение. Можно создать бесконечно много датасетов и очень точно оценить ошибку первого рода вашего критерия. Также, мой опыт говорит, что на начальных этапах дебага нового критерия искусственные данные сильно лучше реальных. Главный минус — вы получили корректность вашего критерия только на искусственных данных! На реальных же данных критерий может работать некорректно.\n",
    "\n",
    "У датасетов, полученных на настоящих данных, всё наоборот: собрать большое количество датасетов сложно, да и не всегда нормально построен процесс сбора логов. Но адекватная оценка корректности критерия для проверки гипотез в вашей компании возможна только таким способом. Всегда можно реализовать такой критерий, который будет правильно работать на искусственных данных. Но, столкнувшись в реальности с более шумными данными, он может начать ошибаться чаще, чем в 5% случаев. Поэтому важно убедиться, что именно на настоящих данных метод будет работать верно. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c004e55f-8ed6-489a-9083-e70e3087516b",
   "metadata": {},
   "source": [
    "\n",
    "----------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de19cddb-efcf-4793-8dab-024bf7c7e2b8",
   "metadata": {},
   "source": [
    "Посмотрим еще раз на табличку:\n",
    "\n",
    "|                          | маленькая выборка | большая выборка |\n",
    "|--------------------------|-------------------|-----------------|\n",
    "| нормальное распределение | t-test            | t-test |\n",
    "| любое распределение      |                   | t-test |\n",
    "\n",
    "\n",
    "Мы уже поняли, как отличить маленькую выборку от большой. Но остался последний вопрос: чем заполнить последнюю пустующую ячейку?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddf67821-3776-4d84-b849-2236b0f6f46e",
   "metadata": {},
   "source": [
    "### Можно ли использовать t-test, если выборка мала и не из нормального распределения?\n",
    "\n",
    "\n",
    "На самом деле ответ простой: проверьте критерий на интересующем вас размере выборки: если FPR $\\leq \\alpha$ значит все хорошо, а если нет &mdash; критерий использовать нельзя.\n",
    "\n",
    "\n",
    "|                          | маленькая выборка | большая выборка |\n",
    "|--------------------------|-------------------|-----------------|\n",
    "| нормальное распределение | t-test               | t-test|\n",
    "| любое распределение      | Монте-Карло проверка | t-test |\n",
    "\n",
    "\n",
    "Например, в случае выше с двумя экспоненциальными распределениями &mdash; нельзя. А если бы ошибка была бы меньше 5% то можно. Например:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "733fcb69-5119-4b65-9516-7bee55662b0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FPR: 0.0396\n",
      "CI=(0.03595116606388123, 0.04360241963583777)\n"
     ]
    }
   ],
   "source": [
    "test_dist    = expon(scale=20)\n",
    "control_dist = expon(scale=20)\n",
    "\n",
    "check_criterion(test_dist=test_dist, control_dist=control_dist, sample_size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b30b8197-1331-488d-a3d0-2478ca8eecdd",
   "metadata": {},
   "source": [
    "Здесь FPR стат. значимо меньше 5%, а значит использовать t-test **можно**. Только надо быть готовым, что он будет не очень мощным. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc0a4522-6eeb-433a-bc6f-93d088c78b88",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b8070d4-94ef-443d-aab5-4f0b0ec3e308",
   "metadata": {},
   "source": [
    "## 2. Какой критерий лучше?\n",
    "\n",
    "Пусть у вас есть 2 критерия, и оба валидны на наших данных. Как понять на практике, какой из них лучше?\n",
    "\n",
    "Правильный ответ &mdash; надо сравнить мощность 2 критериев! Но как ее узнать?\n",
    "\n",
    "Предлагается повторить ту же процедуру, что мы делали выше, только вместо генерации эксперимента, когда верна $H_0$, генерировать эксперимент, когда верна альтернатива. В случае сравнения средних &mdash; надо добавить эффект к тесту. И вместо FPR считать TPR &mdash; сколько раз мы отвергли нулевую гипотезу. Чем больше &mdash; тем лучше.\n",
    "\n",
    "Точно также проверим на t-test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "aa6ae952-8e97-49f1-9037-9b340b43d141",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TPR или мощность: 0.1938\n"
     ]
    }
   ],
   "source": [
    "numpy.random.seed(42)\n",
    "\n",
    "rej_cnt = 0\n",
    "N = 10000\n",
    "alpha=0.05\n",
    "\n",
    "sample_dist = norm(loc=2, scale=3)\n",
    "mu=sample_dist.expect()\n",
    "\n",
    "for i in range(N):\n",
    "    # Генерирую выборку теста и контроля\n",
    "    test    = sample_dist.rvs(15)\n",
    "    control = sample_dist.rvs(15) * 2\n",
    "\n",
    "    # Запускаю критерий и считаю p-value\n",
    "    pvalue = ttest_ind(test, control, equal_var=False, alternative='two-sided').pvalue\n",
    "    \n",
    "    # Проверяю, что pvalue < alpha\n",
    "    rej_cnt += (pvalue < alpha)\n",
    "\n",
    "\n",
    "print(f\"TPR или мощность: {round(rej_cnt / N, 4)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de023f81-c39c-42d3-a1d5-25e39bec8c75",
   "metadata": {},
   "source": [
    "Видим, что мощность критерия в данном случае раняется 19%. Если есть второй критерий &mdash; надо запустить такую проверку для 2го критерия и оценить, какой критерий лучше или хуже, не забыв о стат. значимости. Подробнее вы рассмотрите эту процедуру в домашнем задании.\n",
    "\n",
    "Еще есть вопрос: вы оценили 2 критерия лишь при добавлении одного эффекта, например в случае выше, когда $\\mu_T = \\mu_C * 2$. А если бы было другое изменение, сохранились бы результаты, что этот критерий лучше? Не факт, поэтому поэтому надо идеалогически подбирать такой эффект, который чаще всего встретится на практике. Поэтому ваша задача еще верно сымитировать эффект, похожий на настоящий.\n",
    "Логика здесь точно такая же, как и почему лучше генерировать эксперименты на исторических данных, а не на настоящих.\n",
    "\n",
    "То есть, ваша задача для оценки мощности критерия состоит в:\n",
    "1. Создании 1000 экспериментов, на исторических данных, или на симулированных\n",
    "2. Подборе эффекта, который будет лучше всего имитировать истинный проверяемый эффект в гипотезе.\n",
    "\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2df78501-eedc-4f36-9ba0-869676a41e7d",
   "metadata": {},
   "source": [
    "# Итог\n",
    "\n",
    "В чем заключается метод Монте-Карло в каждой из секции? В генерации большого количества экспериментов и подсчета каких-то статистик на нем.\n",
    "На текущей лекции мы с вами посмотрели разные способы применения Монте-Карло метода:\n",
    "- Как проверять критерий на практике.\n",
    "    - Как генерировать эксперименты для проверки\n",
    "- Как оценивать мощность критерия.\n",
    "<!-- - И как самому приудмать свой критерий, не заниаясь сложной математикой. Метод хороший, но не всегда рабочий. -->\n"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}