bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b34304c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import namedtuple\n",
    "import scipy.stats as sps\n",
    "import statsmodels.stats.api as sms\n",
    "from tqdm.notebook import tqdm as tqdm_notebook # tqdm – библиотека для визуализации прогресса в цикле\n",
    "from collections import defaultdict\n",
    "from statsmodels.stats.proportion import proportion_confint\n",
    "import numpy as np\n",
    "import itertools\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set(font_scale=1.5, palette='Set2')\n",
    "ExperimentComparisonResults = namedtuple('ExperimentComparisonResults', \n",
    "                                        ['pvalue', 'effect', 'ci_length', 'left_bound', 'right_bound'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c62a7e7d",
   "metadata": {},
   "source": [
    "# Bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "428942ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap(control, test, test_type='absolute'):\n",
    "    # Функция от средних, которую надо посчитать на каждой выборке\n",
    "    absolute_func = lambda C, T: T - C\n",
    "    relative_func = lambda C, T: T / C - 1\n",
    "    \n",
    "    boot_func = absolute_func if test_type == 'absolute' else relative_func\n",
    "    stat_sample = []\n",
    "    \n",
    "    batch_sz = 100\n",
    "    \n",
    "    #В теории boot_samples_size стоить брать не меньше размера выборки. Но на практике можно и меньше.\n",
    "    boot_samples_size = len(control)\n",
    "    for i in range(0, boot_samples_size, batch_sz):\n",
    "        N_c = len(control)\n",
    "        N_t = len(test)\n",
    "        # выбираем N_c элементов с повторением из текущей выборки. \n",
    "        # И чтобы ускорить этот процесс, делаем это сразу batch_sz раз\n",
    "        # Вместо одной выборки мы получим batch_sz выборок\n",
    "        control_sample = np.random.choice(control, size=(len(control), batch_sz), replace=True)\n",
    "        test_sample    = np.random.choice(test, size=(len(test), batch_sz), replace=True)\n",
    "\n",
    "        C = np.mean(control_sample, axis=0)\n",
    "        T = np.mean(test_sample, axis=0)\n",
    "        assert len(T) == batch_sz\n",
    "        \n",
    "        # добавляем в массив посчитанных ранее статистик batch_sz новых значений\n",
    "        # X в статье – это boot_func(control_sample_mean, test_sample_mean)\n",
    "        stat_sample += list(boot_func(C, T))\n",
    "\n",
    "    stat_sample = np.array(stat_sample)\n",
    "    # считаем истинный эффект\n",
    "    effect = boot_func(np.mean(control), np.mean(test))\n",
    "    left_bound, right_bound = np.quantile(stat_sample, [0.025, 0.975])\n",
    "    \n",
    "    ci_length = (right_bound - left_bound)\n",
    "    # pvalue - процент статистик, которые лежат левее или правее 0.\n",
    "    pvalue = 2 * min(np.mean(stat_sample > 0), np.mean(stat_sample < 0))\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1053585b",
   "metadata": {},
   "source": [
    "AB-тест"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c00096ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c98465c973b34aa09b3ad28101587889",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Реальный уровень значимости: 0.0525; доверительный интервал: [0.0483, 0.057]\n"
     ]
    }
   ],
   "source": [
    "# 3. Заводим счетчик\n",
    "bad_cnt = 0\n",
    "\n",
    "# 4. Цикл проверки\n",
    "N = 10000\n",
    "for i in tqdm_notebook(range(N)):\n",
    "    # 4.a. Тестирую AB - тест\n",
    "    control = sps.expon(scale=1000).rvs(2000)\n",
    "    test = sps.expon(scale=1000).rvs(2300)\n",
    "    test *= 1.1\n",
    "\n",
    "    # 4.b. Запускаю критерий\n",
    "    _, _, _, left_bound, right_bound = bootstrap(control, test, 'relative')\n",
    "    \n",
    "    # 4.c. Проверяю, лежит ли истинная разница средних в доверительном интервале\n",
    "    if left_bound > 0.1 or right_bound < 0.1:\n",
    "        bad_cnt += 1\n",
    "\n",
    "# 5. Строю доверительный интервал для конверсии ошибок у критерия.\n",
    "left_real_level, right_real_level = proportion_confint(count = bad_cnt, nobs = N, alpha=0.05, method='wilson')\n",
    "# Результат\n",
    "print(f\"Реальный уровень значимости: {round(bad_cnt / N, 4)};\"\n",
    "      f\" доверительный интервал: [{round(left_real_level, 4)}, {round(right_real_level, 4)}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "611b673a",
   "metadata": {},
   "source": [
    "# Пост-нормировка"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6ac62848",
   "metadata": {},
   "outputs": [],
   "source": [
    "def post_normed_bootstrap(control, test, control_before, test_before, test_type='absolute'):\n",
    "    # Функция от средних, которую надо посчитать на каждой выборке\n",
    "    absolute_func = lambda C, T, C_b, T_b: T - (T_b / C_b) * C\n",
    "    relative_func = lambda C, T, C_b, T_b: (T / C) / (T_b / C_b) - 1\n",
    "    \n",
    "    boot_func = absolute_func if test_type == 'absolute' else relative_func\n",
    "    stat_sample = []\n",
    "    \n",
    "    batch_sz = 100\n",
    "    \n",
    "    #В теории boot_samples_size стоить брать не меньше размера выборки. Но на практике можно и меньше.\n",
    "    boot_samples_size = len(control)\n",
    "    for i in range(0, boot_samples_size, batch_sz):\n",
    "        N_c = len(control)\n",
    "        N_t = len(test)\n",
    "        # Надо помнить, что мы семплируем именно юзеров\n",
    "        # Поэтому, если мы взяли n раз i элемент в выборке control\n",
    "        # То надо столько же раз взять i элемент в выборке control_before\n",
    "        # Поэтому будем семплировать индексы\n",
    "        control_indices = np.arange(N_c)\n",
    "        test_indices = np.arange(N_t)\n",
    "        control_indices_sample = np.random.choice(control_indices, size=(len(control), batch_sz), replace=True)\n",
    "        test_indices_sample    = np.random.choice(test_indices, size=(len(test), batch_sz), replace=True)\n",
    "\n",
    "        C   = np.mean(control[control_indices_sample], axis=0)\n",
    "        T   = np.mean(test[test_indices_sample], axis=0)\n",
    "        C_b = np.mean(control_before[control_indices_sample], axis=0)\n",
    "        T_b = np.mean(test_before[test_indices_sample], axis=0)\n",
    "        assert len(T) == batch_sz\n",
    "        stat_sample += list(boot_func(C, T, C_b, T_b))\n",
    "\n",
    "    stat_sample = np.array(stat_sample)\n",
    "    # считаем истинный эффект\n",
    "    effect = boot_func(np.mean(control), np.mean(test), np.mean(control_before), np.mean(test_before))\n",
    "    left_bound, right_bound = np.quantile(stat_sample, [0.025, 0.975])\n",
    "    \n",
    "    ci_length = (right_bound - left_bound)\n",
    "    # pvalue - процент статистик, которые лежат левее или правее 0.\n",
    "    pvalue = 2 * min(np.mean(stat_sample > 0), np.mean(stat_sample < 0))\n",
    "    return ExperimentComparisonResults(pvalue, effect, ci_length, left_bound, right_bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d33fceb2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "41d6fd8298c3460b996b5e064af764ed",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/30000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Реальный уровень значимости: 0.0513; доверительный интервал: [0.0489, 0.0539]\n"
     ]
    }
   ],
   "source": [
    "# 3. Заводим счетчик\n",
    "bad_cnt = 0\n",
    "\n",
    "# 4. Цикл проверки\n",
    "N = 30000\n",
    "cis = []\n",
    "for i in tqdm_notebook(range(N)):\n",
    "    # 4.a. Тестирую AB - тест\n",
    "    control_before = sps.expon(scale=1000).rvs(1000)\n",
    "    control = control_before + sps.norm(loc=0, scale=100).rvs(1000)\n",
    "\n",
    "    test_before = sps.expon(scale=1000).rvs(1000)\n",
    "    test = test_before + sps.norm(loc=0, scale=100).rvs(1000)\n",
    "    test *= 1.1\n",
    "\n",
    "    # 4.b. Запускаю критерий\n",
    "    _, _, ci, left_bound, right_bound = post_normed_bootstrap(control, test, control_before, test_before, 'relative')\n",
    "    cis.append(ci)\n",
    "    # 4.c. Проверяю, лежит ли истинная разница средних в доверительном интервале\n",
    "    if left_bound > 0.1 or right_bound < 0.1:\n",
    "        bad_cnt += 1\n",
    "\n",
    "# 5. Строю доверительный интервал для конверсии ошибок у критерия.\n",
    "left_real_level, right_real_level = proportion_confint(count = bad_cnt, nobs = N, alpha=0.05, method='wilson')\n",
    "# Результат\n",
    "print(f\"Реальный уровень значимости: {round(bad_cnt / N, 4)};\"\n",
    "      f\" доверительный интервал: [{round(left_real_level, 4)}, {round(right_real_level, 4)}]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3af3140",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}