Приложение-Сравнение_методов_оценки_АБ_тестов.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Приложение: сравнение методов оценки А/Б-тестов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Проведено сравнение методов выбора лучше группы в А/Б-тестах по количеству правильно угаданных вариантов*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [Введение](#Введение)\n",
    "* [Способы выбора лучшей группы в А/Б-тесте](#Способы-выбора-лучшей-группы-в-А/Б-тесте)    \n",
    "&nbsp;&nbsp; - [Сравнение средних в сэмплах](#Сравнение-средних-в-сэмплах)  \n",
    "&nbsp;&nbsp; - [Оценка разности средних и выбор лучшей группы с помощью неравенства Чебышева](#Оценка-разности-средних-и-выбор-лучшей-группы-с-помощью-неравенства-Чебышева)  \n",
    "&nbsp;&nbsp; - [Оценки разности средних с учетом центральной предельной теоремы](#Оценки-разности-средних-с-учетом-центральной-предельной-теоремы)  \n",
    "&nbsp;&nbsp; - [Метод максимального правдоподобия](#Метод-максимального-правдоподобия)  \n",
    "&nbsp;&nbsp; - [Байесовское моделирование](#Байесовское-моделирование)  \n",
    "&nbsp;&nbsp; - Проверка статистических гипотез  \n",
    "&nbsp;&nbsp; - [Бутстрап](#Бутстрап)  \n",
    "* [Количество правильно угаданных вариантов](#Количество-правильно-угаданных-вариантов)\n",
    "* Выводы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Введение"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть разные методы выбора лучшей группы в А/Б-тесте.  \n",
    "Интересно сравнить их на модельных примерах - какие методы чаще и за меньшее количество данных угадывают лучшие варианты.   \n",
    "\n",
    "Также полезно сравнить, как в рамках различных методов даются ответы на вопросы:\n",
    "- Какой вариант лучше и насколько?\n",
    "- Каковы оценки целевой метрики в каждом варианте?\n",
    "- Насколько уверены в оценке?\n",
    "- Сколько должен продолжаться эксперимент?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Способы выбора лучшей группы в А/Б-тесте"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пусть есть 2 распределения случайных величин одинаковой формы, но с разными параметрами.  \n",
    "Загадываются параметры, генерируются выборки.  \n",
    "По выборкам нужно сделать оценку параметров и выбрать \"лучшую\" группу.  \n",
    "\n",
    "Для характеристики распределений часто проще использовать ожидаемое значение.  \n",
    "Т.е. по выборкам нужно определить распределение с большим средним.\n",
    "\n",
    "Конверсии часто моделируют независимыми испытаниями с двумя исходами (испытания Бернулли)\n",
    "[[Bern](https://en.wikipedia.org/wiki/Bernoulli_process)].   \n",
    "Выручку на платящего пользователя иногда моделируют распределением Парето.   \n",
    "Реальное поведение пользователей и метрик может требовать более сложных распределений."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "np.random.seed(7)\n",
    "\n",
    "import scipy.stats as stats\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "from plotly.subplots import make_subplots\n",
    "\n",
    "#todo: update scipy; make venv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "p_a_exact = 0.3\n",
    "sa = stats.binom.rvs(n, p_a_exact)\n",
    "p_b_exact = 0.33\n",
    "sb = stats.binom.rvs(n, p_b_exact)\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Bar(x=[1, 0], y=[sa, n-sa], name='Sample A', marker_color='red'))\n",
    "fig.add_trace(go.Bar(x=[1, 0], y=[sb, n-sb], name='Sample B', marker_color='blue'))\n",
    "fig.add_trace(go.Scatter(x=[p_a_exact, p_a_exact], y=[0, max(sa, n-sa)], \n",
    "                         line_dash='dash', mode='lines', name='Exact Success Prob A', line_color='red'))\n",
    "fig.add_trace(go.Scatter(x=[p_b_exact, p_b_exact], y=[0, max(sb, n-sb)], \n",
    "                         line_dash='dash', mode='lines', name='Exact Success Prob B', line_color='blue'))\n",
    "fig.update_layout(\n",
    "    title='Bernoulli Sequence Samples and Exact Success Probabilities for Two Groups',\n",
    "    xaxis_range=[-0.1, 1.1],\n",
    "    barmode='group', bargap=0.9, bargroupgap=0.0,\n",
    "    height=500, width=850)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 1000\n",
    "\n",
    "c_a_exact = 2.5\n",
    "exact_dist_a = stats.lomax(c=c_a_exact)\n",
    "samp_a = exact_dist_a.rvs(size=n)\n",
    "\n",
    "c_b_exact = 2.8\n",
    "exact_dist_b = stats.lomax(c=c_b_exact)\n",
    "samp_b = exact_dist_b.rvs(size=n)\n",
    "\n",
    "# mean = 1 / (c - 1), c > 1\n",
    "\n",
    "# x = np.linspace(0, 10, 2000)\n",
    "# ymax = np.max([exact_dist_a.pdf(x), exact_dist_b.pdf(x)])\n",
    "ymax = 1\n",
    "\n",
    "fig = go.Figure()\n",
    "col = 'red'\n",
    "# fig.add_trace(go.Scatter(x=x, y=exact_dist_a.pdf(x), \n",
    "#                          mode='lines', line_color=col,\n",
    "#                          name=f\"Exact A: a={exp_info['a']['A']}, b={exp_info['b']['A']}\"))\n",
    "fig.add_trace(go.Histogram(x=samp_a, histnorm='probability density',\n",
    "                           name='Sample A',\n",
    "                           opacity=0.3, marker_color=col))\n",
    "fig.add_trace(go.Scatter(x=[exact_dist_a.mean(), exact_dist_a.mean()], y=[0, ymax], \n",
    "                         mode='lines', line_color=col, line_dash='dash',\n",
    "                         name='Exact Mean A'))\n",
    "col = 'blue'\n",
    "# fig.add_trace(go.Scatter(x=x, y=exact_dist_b.pdf(x), \n",
    "#                          mode='lines', line_color=col,\n",
    "#                          name=f\"Exact B: a={exp_info['a']['B']}, b={exp_info['b']['B']}\"))\n",
    "fig.add_trace(go.Histogram(x=samp_b, histnorm='probability density',\n",
    "                           name='Sample B',\n",
    "                           opacity=0.3, marker_color=col))\n",
    "fig.add_trace(go.Scatter(x=[exact_dist_b.mean(), exact_dist_b.mean()], y=[0, ymax], \n",
    "                         mode='lines', line_color=col, line_dash='dash',\n",
    "                         name='Exact Mean B'))\n",
    "fig.update_layout(title='Pareto Exact Distributions and Samples',\n",
    "                  xaxis_title='x',\n",
    "                  yaxis_title='Prob Density',\n",
    "                  hovermode=\"x\",\n",
    "                  height=550,\n",
    "                  barmode='overlay')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Сравнение средних в сэмплах"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "По закону больших чисел среднее в выборке сходится к среднему в распределении при росте количества элементов в выборке [[LLN](https://en.wikipedia.org/wiki/Law_of_large_numbers)]\n",
    "\n",
    "$$\n",
    "\\lim_{n \\to \\infty} \\sum_{i=1}^{n} \\frac{X_i}{n} = \\bar{X} .\n",
    "$$\n",
    "\n",
    "Сходимость понимается как сходимость по распределению (слабый закон больших чисел).  \n",
    "*Иначе можно сказать, что среднее в сэмпле - состоятельная несмещенная оценка среднего в распределении* https://en.wikipedia.org/wiki/Estimator ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сравнение средних в выборке с точным средним по мере роста количества точек выборки представлено на графике ниже. *Видно, что при размере выборки в несколько сот точек приближение отличия в пределах 5%, а при выбранных параметрах при $n > 1000$ относительное отличие не превышает 5%.* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_a_exact = 0.7\n",
    "p_b_exact = 0.6\n",
    "\n",
    "n = np.arange(10, 5000, 100)\n",
    "k_a = stats.binom.rvs(n, p_a_exact)\n",
    "k_b = stats.binom.rvs(n, p_b_exact)\n",
    "p_a_samp = k_a/n\n",
    "p_b_samp = k_b/n\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(x=n, y=p_a_samp, name='A: p sample'))\n",
    "fig.add_trace(go.Scatter(x=[0, np.max(n)], y=[p_a_exact, p_a_exact],\n",
    "                         mode='lines',\n",
    "                         line_dash='dash', line_color='black', name='A: p exact'))\n",
    "fig.add_trace(go.Scatter(x=n, y=p_b_samp, name='B: p sample'))\n",
    "fig.add_trace(go.Scatter(x=[0, np.max(n)], y=[p_b_exact, p_b_exact],\n",
    "                         mode='lines',\n",
    "                         line_dash='dash', line_color='black', name='B: p exact'))\n",
    "fig.update_layout(\n",
    "    title='Exact and Sample Means for Binomial Distribution',\n",
    "    xaxis_title='Sample Size',\n",
    "    yaxis_title='p',\n",
    "    yaxis_range=[0, 1],\n",
    "    xaxis_range=[0, np.max(n)],\n",
    "    height=550\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Поскольку выборочные средние приближают среднее в распределении, то для выбора распределения с большим средним можно выбрать группу с большим средним в выборке."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_binomial_sample_means(sa, na, sb, nb):\n",
    "    pa = sa / na\n",
    "    pb = sb / nb\n",
    "    if pa > pb:\n",
    "        res = 'A'\n",
    "    elif pa < pb:\n",
    "        res = 'B'\n",
    "    else:\n",
    "        res = None\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_a_exact = 0.7\n",
    "p_b_exact = 0.75\n",
    "\n",
    "n = 1000\n",
    "sa = stats.binom.rvs(n, p_a_exact)\n",
    "sb = stats.binom.rvs(n, p_b_exact)\n",
    "\n",
    "print(f'pa_exact={p_a_exact}, pb_exact={p_b_exact}')\n",
    "print(f'n={n}, sa={sa}, sb={sb}')\n",
    "print(f'Exact best group: {\"A\" if p_a_exact > p_b_exact else \"B\" if p_b_exact > p_a_exact else None}')\n",
    "print(f'Group selected by sample means comparison: {compare_binomial_sample_means(sa, n, sb, n)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предложенный метод позволяет сделать оценки интересующих метрик и выбрать одну из групп в А/Б-тесте, но не дает оценки вероятности, с какой одна группа лучше другой. Также не обсуждался вопрос о длительности эксперимента.\n",
    "\n",
    "*todo: посмотреть скорость сходимости. Если она считается иначе, чем в центральной предельной теореме, можно попробовать получить оценку длительности из нее*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Метод                 | Оценки целевой метрики                              | Критерий выбора группы   | Оценка длительности |\n",
    "|-----------------------|-----------------------------------------------------|--------------------------|---------------------|\n",
    "| **Сравнение средних** | Точечные оценки средних на основе средних в выборке | Больше среднее в выборке | (?)                 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Оценка разности средних и выбор лучшей группы с помощью неравенства Чебышева"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Посчитать средние и дисперсии в каждой группе\n",
    "* Посчитать стандартные ошибки средних\n",
    "* С помощью неравенства Чебышева оценить интервал возможных значений\n",
    "* Построить разность средних и дисперсию разности средних\n",
    "* Посчитать вероятность разности больше 0. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A = [1, 0, 1, 0, 0, ... , 1]\n",
    "\\\\\n",
    "s_A, n_A\n",
    "\\\\\n",
    "E_A \\equiv p_A = s_A / n_A\n",
    "\\\\\n",
    "D_A = p_A (1 - p_A)\n",
    "\\\\\n",
    "\\tilde{p}_A \\equiv \\frac{1}{n} \\left( \\tilde{A} + \\tilde{A} + \\dots + \\tilde{A} \\right)\n",
    "\\\\\n",
    "E[\\tilde{p}_A] = E_A = p_A\n",
    "\\\\\n",
    "D[\\tilde{p}_A] = \\frac{D_A}{n_A}\n",
    "\\\\\n",
    "\\Delta = \\tilde{p}_A - \\tilde{p}_B\n",
    "\\\\\n",
    "E_{\\Delta} = E[\\tilde{p}_A - \\tilde{p}_B] = E[\\tilde{p}_A] - E[\\tilde{p}_B] = p_A - p_B\n",
    "\\\\\n",
    "D_{\\Delta} = D[\\tilde{p}_A - \\tilde{p}_B] = D[\\tilde{p}_A] + D[\\tilde{p}_B]\n",
    "\\\\\n",
    "\\sigma_{\\Delta} = \\sqrt{D[p_A] + D[p_B]}\n",
    "\\\\\n",
    "P(|0 - E_{\\Delta}| \\ge k \\sigma_{\\Delta}) \\le \\frac{1}{k^2}\n",
    "\\\\\n",
    "P \\left( \\frac{|E_{\\Delta}|}{\\sigma_{\\Delta}} \\ge k \\right) \\le \\frac{1}{k^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_a_exact = 0.3\n",
    "n = 1000\n",
    "sa = stats.binom.rvs(n, p_a_exact)\n",
    "\n",
    "p_samp = sa / n\n",
    "stderr = np.sqrt(p_samp * (1-p_samp) / n)\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Bar(x=[1, 0], y=[sa, n-sa], width=0.1, name='Sample'))\n",
    "fig.add_trace(go.Scatter(x=[p_a_exact, p_a_exact], y=[0, max(sa, n-sa)], \n",
    "                         line_dash='solid', mode='lines', name='Exact Mean'))\n",
    "xtmp = np.arange(0, 1, 0.001)\n",
    "ytmp = stats.norm.pdf(x=xtmp, loc=p_samp, scale=stderr)\n",
    "ytmp = ytmp / max(ytmp) * max(sa, n-sa)\n",
    "fig.add_trace(go.Scatter(x=xtmp, y=ytmp,\n",
    "                        line_dash='dash', name='Mean Estimate from Sample'))\n",
    "fig.update_layout(\n",
    "    xaxis_range=[0,1],\n",
    "    height=500, width=850)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Оценку среднего значения неизвестного распределения по выборке можно рассматривать как построение оценки случайной величины. \n",
    "Предыдущий метод давал точечную оценку среднего, но не рассматривал его неопределенность."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если исходное распределение с конечной диспресией, то \n",
    "дисперсия выборочных средних задается выражением \n",
    "(стандартная ошибка среднего [[StErrMn](https://en.wikipedia.org/wiki/Standard_error#Derivation)])\n",
    "\n",
    "$$\n",
    "D[{\\bar{x}}] = D \\left[ \\frac{1}{n} \\sum_{i=1}^n x_i \\right] = \\frac{1}{n^2} n D[x] = \\frac{\\sigma^2}{n} ,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "где $\\sigma$ - стандартное отклонение исходного распределения, $n$ - количество точек в выборочном среднем. \n",
    "\n",
    "Оценку дисперсии $\\sigma$ можно построить по сэмлу.  \n",
    "Для дисперсии нужны поправки к значению в сэмпле $s$ чтобы получить несмещенную оценку [[StdUnbiased](https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation)]  \n",
    "\n",
    "$$\n",
    "\\sigma \\approx \\sqrt{\\frac{n}{n-1}} s, \\quad s = \\sqrt{ \\frac{1}{n} \\sum_{i=1}^n (x_i - \\bar{x}_i)^2 } .\n",
    "$$\n",
    "\n",
    "При больших $n$ можно для оценки стандартного отклонения распределения использовать стандартное отклонение в выборке $s$ без корректирующих правок\n",
    "\n",
    "$$\n",
    "\\sigma \\approx s .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Оценить вероятность отклонения случайной величины от среднего можно с помощью неравенства Чебышева \n",
    "[[Cheb](https://en.wikipedia.org/wiki/Chebyshev%27s_inequality)]   \n",
    "\n",
    "$$\n",
    "P(|X - \\mu| \\ge k \\sigma) \\le \\frac{1}{k^2} .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для двух групп можно получить оценки средних, дисперсий и интервалов вероятных значений.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для сравнения групп можно посчитать вероятность $P( E[A - B] > 0 )$.  \n",
    "*todo: ввести обозначение для оценки распределения средних.*  \n",
    "*что-то типа $\\tilde{p}$*  \n",
    "*тогда нужно оценить $P(\\tilde{p}_A - \\tilde{p}_B > 0)$*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Среднее разности случайных величин линейно и считается как разность средних [[ExpProp](https://en.wikipedia.org/wiki/Expected_value#Properties)]  \n",
    "\n",
    "$$\n",
    "E[X - Y] = E[X] - E[Y] .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Дисперсия разности случайных величин \n",
    "[[VarProp](https://en.wikipedia.org/wiki/Standard_deviation#Identities_and_mathematical_properties)] \n",
    "\n",
    "$$\n",
    "D[X - Y] = D[X] + D[Y] - 2cov[X,Y] .\n",
    "$$\n",
    "\n",
    "В данных примерах величины X и Y независимы.  \n",
    "Для независимых величин $cov[X,Y] = 0$ [[Cov](https://en.wikipedia.org/wiki/Covariance)].  \n",
    "В реальном А/Б-тесте это слагаемое можно оставить."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Таким образом удается оценить ожидаемое значение и стандартное отклонение разности средних двух групп  \n",
    "\n",
    "$$\n",
    "E_{\\Delta} = E[A] - E[B], \n",
    "\\quad\n",
    "\\sigma_{\\Delta} = \\sqrt{D[p_A] + D[p_B]} .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "С помощью неравенства Чебышева можно оценить вероятность среднее в одной группе больше другой"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "P(|0 - E_{\\Delta}| \\ge k \\sigma_{\\Delta}) \\le \\frac{1}{k^2}\n",
    "\\\\\n",
    "P \\left( \\frac{|E_{\\Delta}|}{\\sigma_{\\Delta}} \\ge k \\right) \\le \\frac{1}{k^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Одностороннее неравенство Чебышева: https://en.wikipedia.org/wiki/Chebyshev%27s_inequality#Cantelli's_inequality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def std_err_mean_binom(s, n):\n",
    "    p = s / n\n",
    "    return np.sqrt(p * (1 - p) / n)\n",
    "\n",
    "def compare_binomial_stderr_chebyshev(sa, na, sb, nb, diff_prob=0.9):\n",
    "    pa = sa / na\n",
    "    pb = sb / nb\n",
    "    std_err_mean_a = std_err_mean_binom(sa, na)\n",
    "    std_err_mean_b = std_err_mean_binom(sb, nb)\n",
    "    mean_diff = pb - pa\n",
    "    std_err_mean_diff = np.sqrt(std_err_mean_a**2 + std_err_mean_b**2)\n",
    "    k = np.abs(mean_diff) / std_err_mean_diff\n",
    "    prob = 1.0 / k**2\n",
    "    if pa > pb and prob <= diff_prob:\n",
    "        res = 'A'\n",
    "    elif pa < pb and prob <= diff_prob:\n",
    "        res = 'B'\n",
    "    else:\n",
    "        res = None\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_a_exact = 0.7\n",
    "p_b_exact = 0.75\n",
    "\n",
    "n = 1000\n",
    "sa = stats.binom.rvs(n, p_a_exact)\n",
    "sb = stats.binom.rvs(n, p_b_exact)\n",
    "\n",
    "print(f'pa_exact={p_a_exact}, pb_exact={p_b_exact}')\n",
    "print(f'n={n}, sa={sa}, sb={sb}')\n",
    "print(f'Exact best group: {\"A\" if p_a_exact > p_b_exact else \"B\" if p_b_exact > p_a_exact else None}')\n",
    "print(f'Group selected by means comparison: {compare_binomial_stderr_chebyshev(sa, n, sb, n)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Метод не делает предположений о форме распределений тестируемых величин. Единственное, что требуется - существование конечных средних и дисперсий. Чаще всего эти предположения можно считать выполненными. Однако при конечных среднем и дисперсии обычно применима центальная предельная теорема, на основании которой можно построить более точную форму распределения и использовать его для оценок вероятностей вместо неравенства Чебышева - см. следующий раздел."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Метод                                        | Оценки целевой метрики                                                                        | Критерий выбора группы                                 | Оценка длительности |\n",
    "|----------------------------------------------|-----------------------------------------------------------------------------------------------|--------------------------------------------------------|---------------------|\n",
    "| **Сравнение средних + неравенство Чебышева** | Точечные оценки средних на основе средних в выборке,  интервалы с помощью неравества Чебышева | Вероятность P(A-B > 0) на основе  неравенства Чебышева | (?)                 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Остается вопрос интерпретации вероятности и интервалов вокруг средних значений.*  \n",
    "*Вводится новая случайная величина - сумма N наблюдений.*  \n",
    "*Делаются утверждения о ее распределении. К тому же приближенные.*  \n",
    "*Как это соотносится с исходными вопросами?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Оценки разности средних с учетом центральной предельной теоремы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Если у распределения конечная дисперсия, то можно применять центральные предельные теоремы [[CLT](https://en.wikipedia.org/wiki/Central_limit_theorem)].  \n",
    "Это дает более строгую оценку на область вероятных значений средних, чем неравенство Чебышева."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для серий Бернулли см. также теорему Муавра-Лапласа [[DMLT](https://en.wikipedia.org/wiki/De_Moivre%E2%80%93Laplace_theorem)] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "S_n = \\xi_1 + ... + \\xi_n\n",
    "\\\\\n",
    "a = E[\\xi],\n",
    "\\qquad\n",
    "\\sigma^2 = D[\\xi]\n",
    "\\\\\n",
    "P\\left(\\alpha < \\frac{S_n - na}{\\sqrt{n}\\sigma} < \\beta \\right) \\to \\frac{1}{2 \\pi} \\int_{\\alpha}^{\\beta} e^{-t^2/2} dt \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Центральная предельная теорема позволяет по выборке построить интервал, который будет содержать среднее значение с заданной вероятностью:\n",
    "\n",
    "$$ \n",
    "P\\left( \\frac{S_n - \\sqrt{n}\\sigma \\beta}{n} < a < \\frac{S_n - \\sqrt{n}\\sigma \\alpha}{n} \\right) \\to \\frac{1}{2 \\pi} \\int_{\\alpha}^{\\beta} e^{-t^2/2} dt \n",
    "$$\n",
    "\n",
    "Такие конструкции называются доверительными интервалами [[ConfInt](https://en.wikipedia.org/wiki/Confidence_interval)]. Частотная интерпретация этой вероятности - провести эксперимент и получить выборку из n элементов, по выборке построить границы интервала, с заданной вероятность (например, в 95 повторах из 100) реальное среднее будет лежать внутри интервала. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для оценки разности можно ввести случайную величину $\\xi_{\\Delta} = \\xi_B - \\xi_A$. Для нее использовать выражение выше напрямую.  \n",
    "Если число точек в группе неодинаковое, это менее удобно."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно приближенно считать $\\tilde{p}_A$ и $\\tilde{p}_B$ распределенными нормально\n",
    "\n",
    "$$\n",
    "\\tilde{p}_A \\sim Norm(p_A, s_A), \\quad \\tilde{p}_B \\sim Norm(p_B, s_B) .\n",
    "$$\n",
    "\n",
    "Далее воспользоваться свойством нормальных распределений по которому разность нормальных распределений также будет распределена нормально [[NormSum](https://en.wikipedia.org/wiki/Sum_of_normally_distributed_random_variables)]\n",
    "$$\n",
    "\\tilde{p}_{\\Delta} \\equiv \\tilde{p}_B - \\tilde{p}_A ,\n",
    "\\\\\n",
    "\\tilde{p}_{\\Delta} \\sim Norm(p_B - p_A, \\sqrt{\\sigma_A^2 + \\sigma_B^2}) .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def std_err_mean_binom(s, n):\n",
    "    p = s / n\n",
    "    return np.sqrt(p * (1 - p) / n)\n",
    "\n",
    "def compare_binomial_clt(sa, na, sb, nb, diff_prob=0.9):\n",
    "    pa = sa / na\n",
    "    pb = sb / nb\n",
    "    std_err_mean_a = std_err_mean_binom(sa, na)\n",
    "    std_err_mean_b = std_err_mean_binom(sb, nb)\n",
    "    mean_diff = pb - pa\n",
    "    std_err_mean_diff = np.sqrt(std_err_mean_a**2 + std_err_mean_b**2)\n",
    "    #\n",
    "    pb_gt_pa = 1 - stats.norm.cdf(x=0, loc=mean_diff, scale=std_err_mean_diff)\n",
    "    #\n",
    "    if pa > pb and diff_prob <= 1 - pb_gt_pa:\n",
    "        res = 'A'\n",
    "    elif pa < pb and diff_prob <= pb_gt_pa:\n",
    "        res = 'B'\n",
    "    else:\n",
    "        res = None\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Метод                                                  | Оценки целевой метрики                                                                                   | Критерий выбора группы                                           | Оценка длительности |\n",
    "|--------------------------------------------------------|----------------------------------------------------------------------------------------------------------|------------------------------------------------------------------|---------------------|\n",
    "| **Сравнение средних + центральная предельная теорема** | Точечные оценки средних на основе средних в выборке,  интервалы с помощью центральной предельной теоремы | Вероятность P(A-B > 0) на основе  центральной предельной теоремы | (?)                 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Метод максимального правдоподобия"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предыдущие методы не делали предположений о форме распределения сравниваемых случайных величин.  \n",
    "При учете таких предположений сравнение может быть эффективнее.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В методе максимального правдоподобия [[MxLk](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation)] делается предположение о распределении случайной величины.  \n",
    "Эту функцию называют функцией правдоподобия. Обычно у этой функции есть параметы.   \n",
    "Значения параметров выбирается так, чтобы максимизировать эту функцию при имеющихся данных.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для некоторых моделей можно проверси аналитические вычисления.  \n",
    "Для биномиального распределения максимум правдоподобия совпадет со средним в сэмпле.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "P(s, n) = C^n_s p^s (1-p)^{n-s}\n",
    "\\\\\n",
    "\\frac{d P}{dp} = C s p^{s-1} (1-p)^{n-s} - C p^s (n-s) (1-p)^{n-s-1} = 0\n",
    "\\\\\n",
    "C p^{s-1} (1-p)^{n-s-1} \\Big[ s (1-p) - p (n-s) \\Big] = 0\n",
    "\\\\\n",
    "s - pn = 0\n",
    "\\\\\n",
    "p = \\frac{s}{n}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Т.е. когда $n$ бросков монетки и $s$ успехов моделируются биномиальным распределением, то получить пару $(n,s)$ с наибольшей вероятностью можно при значении $p = s/n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Для нормального распределения можно показать, что параметры, максимизирующие функцию правдоподобия, совпадут со средним и дисперсией в сэмле. (?)* "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Выбрана модель и оценены параметры, дальше можно сравнивать эти распределения.  \n",
    "Если есть возможность посчитать средние по параметрам аналитически, то можно выбрать группу с большим средним.  \n",
    "\n",
    "Для биномиальных распределений его результаты совпадают с простым сравнением средних."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Метод формально говорит о полной уверенности."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_binomial_maxlikelihood(sa, na, sb, nb):\n",
    "    pa = sa / na\n",
    "    pb = sb / nb\n",
    "    if pa > pb:\n",
    "        res = 'A'\n",
    "    elif pa < pb:\n",
    "        res = 'B'\n",
    "    else:\n",
    "        res = None\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Метод максимального правдоподобия можно рассматривать как частный случай байесовского моделирования.  \n",
    "В байесовском подходе параметры интерпретируются как случайные величины с предполагаемым распределением.  \n",
    "В методе максимального правдоподобия - считаются неизвестной постоянной величиной.\n",
    "Оценка параметров методом максимального правдоподобия совпадет с максимум апостериорного распределения в байесовском моделировании, если априорное распределение равномерное [[MxLk](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation#Relation_to_Bayesian_inference), [MAPE](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation)].  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Метод оценивает точечные оценки параметров, а не распределения.  \n",
    "Поскольку вместо распределения параметров используется точечная оценка, метод будет давать завышенную уверенность по сравнению со сравнением апостериорных распределений в байесовском подходе.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Метод                          | Оценки целевой метрики                                                                                                                | Критерий выбора группы                      | Оценка длительности |\n",
    "|--------------------------------|---------------------------------------------------------------------------------------------------------------------------------------|---------------------------------------------|---------------------|\n",
    "| **Максимальное правдоподобие** | Предположение формы распределения;  оценка параметров, максимизирующих  функцию правдоподобия; сравнение характеристик распределения. | Сравнение интересующих параметров распределения. | (?)                 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Байесовское моделирование"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Использование байесовского моделирования для сравнения групп в А/Б-тестах рассмотрено в предыдущих частях."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Основная идея - предположить аналитический вид распределения (\"модель\"). С помощью соотношения Байеса построить плотность вероятности параметров модели. Сравнивать эти распределения. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def posterior_binom_beta(sa, n, alpha=1, beta=1):\n",
    "    alpha_post = alpha + sa\n",
    "    beta_post = beta + (n - sa)\n",
    "    return stats.beta(alpha_post, beta_post)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pa_exact = 0.7\n",
    "pb_exact = 0.65\n",
    "\n",
    "n = 1000\n",
    "sa = stats.binom.rvs(n, pa_exact)\n",
    "sb = stats.binom.rvs(n, pb_exact)\n",
    "\n",
    "post_dist_a = posterior_binom_beta(sa, n)\n",
    "post_dist_b = posterior_binom_beta(sb, n)\n",
    "\n",
    "x = np.linspace(0, 1, 1000)\n",
    "ymax = np.max([post_dist_a.pdf(x), post_dist_b.pdf(x)])\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(x=x, y=post_dist_a.pdf(x), name='A'))\n",
    "fig.add_trace(go.Scatter(x=[pa_exact, pa_exact], y=[0, ymax], \n",
    "                         mode='lines', line_dash='dash', line_color='black', \n",
    "                         name='A: p exact'))\n",
    "fig.add_trace(go.Scatter(x=x, y=post_dist_b.pdf(x), name='B'))\n",
    "fig.add_trace(go.Scatter(x=[pb_exact, pb_exact], y=[0, ymax],\n",
    "                         mode='lines', line_dash='dash', line_color='black', \n",
    "                         name='B: p exact'))\n",
    "fig.update_layout(title='Posterior',\n",
    "                  xaxis_title='p',\n",
    "                  yaxis_title='Prob',\n",
    "                  xaxis_range=[0.5, 1],\n",
    "                  hovermode=\"x\",\n",
    "                  height=550,\n",
    "                  barmode='overlay')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_groups_bayes_binomial(sa, na, sb, nb, diff_prob=0.9):\n",
    "    n_post_sample = 30000\n",
    "    post_a = posterior_binom_beta(sa, na)\n",
    "    post_b = posterior_binom_beta(sb, nb)\n",
    "    post_samp_a = post_a.rvs(n_post_sample)\n",
    "    post_samp_b = post_b.rvs(n_post_sample)\n",
    "    pa_gt_pb = np.sum(post_samp_a > post_samp_b) / len(post_samp_a)\n",
    "    res = None\n",
    "    if pa_gt_pb >= diff_prob:\n",
    "        res = 'A'\n",
    "    elif pa_gt_pb <= (1 - diff_prob):\n",
    "        res = 'B'\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_a_exact = 0.7\n",
    "p_b_exact = 0.65\n",
    "\n",
    "n = 1000\n",
    "sa = stats.binom.rvs(n, p_a_exact)\n",
    "sb = stats.binom.rvs(n, p_b_exact)\n",
    "\n",
    "print(f'pa_exact={p_a_exact}, pb_exact={p_b_exact}')\n",
    "print(f'n={n}, sa={sa}, sb={sb}')\n",
    "print(f'Exact best group: {\"A\" if p_a_exact > p_b_exact else \"B\" if p_b_exact > p_a_exact else None}')\n",
    "print(f'Group selected by bayesian comparison: {compare_groups_bayes_binomial(sa, n, sb, n)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Количество правильно угаданных вариантов"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно предложить несколько способов выбора лучшей группы в А/Б-тесте."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_methods = pd.DataFrame(columns=[\n",
    "    'Метод', \n",
    "    'Оценки целевой метрики', \n",
    "    'Критерий выбора группы', \n",
    "    'Оценка длительности'\n",
    "])\n",
    "df_methods = pd.concat([df_methods, pd.DataFrame([{\n",
    "    'Метод': 'Сравнение средних', \n",
    "    'Оценки целевой метрики': 'Точечные оценки средних на основе средних в выборке', \n",
    "    'Критерий выбора группы': 'Больше среднее в выборке', \n",
    "    'Оценка длительности': '(?)'\n",
    "}])], ignore_index=True)\n",
    "\n",
    "df_methods = pd.concat([df_methods, pd.DataFrame([{\n",
    "    'Метод': 'Сравнение средних + неравенство Чебышева', \n",
    "    'Оценки целевой метрики': 'Точечные оценки средних на основе средних в выборке, интервалы с помощью неравества Чебышева', \n",
    "    'Критерий выбора группы': 'Вероятность P(A-B > 0) на основе неравенства Чебышева', \n",
    "    'Оценка длительности': '(?)'\n",
    "}])], ignore_index=True)\n",
    "\n",
    "df_methods = pd.concat([df_methods, pd.DataFrame([{\n",
    "    'Метод': 'Центральная предельная теорема', \n",
    "    'Оценки целевой метрики': 'Точечные оценки средних на основе средних в выборке, интервалы по центральной предельной теореме', \n",
    "    'Критерий выбора группы': 'Вероятность P(A-B > 0) на основе центральной предельной теоремы', \n",
    "    'Оценка длительности': '(?)'\n",
    "}])], ignore_index=True)\n",
    "\n",
    "df_methods = pd.concat([df_methods, pd.DataFrame([{\n",
    "    'Метод': 'Максимальное правдоподобие', \n",
    "    'Оценки целевой метрики': 'Предположение формы распределения;  оценка параметров, максимизирующих  функцию правдоподобия; сравнение характеристик распределения.', \n",
    "    'Критерий выбора группы': 'Сравнение интересующих параметров распределения', \n",
    "    'Оценка длительности': '(?)'\n",
    "}])], ignore_index=True)\n",
    "\n",
    "df_methods = pd.concat([df_methods, pd.DataFrame([{\n",
    "    'Метод': 'Байесовское моделирование', \n",
    "    'Оценки целевой метрики': 'Построение апостериорных распределений параметров', \n",
    "    'Критерий выбора группы': 'Сравенние апостериорных распределений', \n",
    "    'Оценка длительности': '(?)'\n",
    "}])], ignore_index=True)\n",
    "\n",
    "with pd.option_context('display.max_colwidth', 0):\n",
    "    display(df_methods.set_index('Метод'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Интересно сравнить эти методы по количеству правильно угаданных вариантов."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Первая задача - выбрать биномиальное распределение с большим значением конверсии из двух распределений по выборке фиксированного размера.  \n",
    "\n",
    "Сравнение будет проводиться следующим образом.  \n",
    "Выбирается количество наблюдений $N$.  \n",
    "Для фиксированного количества наблюдений генерируется $k$ пар параметров $p_A, p_B$.  \n",
    "Для каждой пары генерируется число успехов $s_A, s_B$.  \n",
    "По известным $s_A, s_B, N$ нужно выбрать группу с большим значением $p$.  \n",
    "Сравнить выбранную группу с реально лучшей.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Генерация данных для эксперимента"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_binomial_experiments(n_points, k_experiments, pa_base=0.1):\n",
    "    pa = np.full(k_experiments, pa_base)\n",
    "    pb = stats.uniform.rvs(loc=0.9*pa, scale=0.1*pa) #unif[loc, loc+scale]\n",
    "    sa = stats.binom.rvs(n=n_points, p=pa, size=k_experiments)\n",
    "    sb = stats.binom.rvs(n=n_points, p=pb, size=k_experiments)\n",
    "    e = {\n",
    "        'exp_type': 'binomial',\n",
    "        'k_experiments': k_experiments,\n",
    "        'n_points': n_points,\n",
    "        'pa_exact': pa,\n",
    "        'pb_exact': pb,\n",
    "        'sa': sa,\n",
    "        'sb': sb\n",
    "    }\n",
    "    return e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Сравнение экспериментов и подсчет количества правильно угаданных вариантов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_experiments_binomial_exact(exps):\n",
    "    s = np.empty_like(exps['pa_exact'], dtype=object)\n",
    "    s[exps['pa_exact'] > exps['pb_exact']] = 'A'\n",
    "    s[exps['pb_exact'] > exps['pa_exact']] = 'B'\n",
    "    res = {\n",
    "        'comparison_type': 'exact',\n",
    "        'selected': s\n",
    "    }\n",
    "    return res\n",
    "\n",
    "def compare_experiments(exps, method_name, method_fn):\n",
    "    s = np.empty(k_experiments, dtype=object)\n",
    "    n = exps['n_points']\n",
    "    for i, (sa, sb) in enumerate(zip(exps['sa'], exps['sb'])):\n",
    "        s[i] = method_fn(sa=sa, na=n, sb=sb, nb=n)\n",
    "    res = {\n",
    "        'comparison_type': method_name,\n",
    "        'selected': s\n",
    "    }\n",
    "    return res\n",
    "\n",
    "def method_guesses(method, exact, n_points):\n",
    "    e = exact['selected']\n",
    "    m = method['selected']\n",
    "    s = {\n",
    "        'name': method['comparison_type'],\n",
    "        'n_points': n_points,\n",
    "        'experiments': len(e),\n",
    "        'correct': np.sum(e == m),\n",
    "        'not_sure': len(e) - np.count_nonzero(m),\n",
    "        'incorrect': np.count_nonzero(m) - np.sum(e == m)\n",
    "    }\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_points = np.array([10, 100, 1000, 5000, 10000, 50000, 100000, 500000, 1000000])\n",
    "#n_points = 100000\n",
    "k_experiments = 100\n",
    "pa_base = 0.1\n",
    "\n",
    "comparisons = {\n",
    "    'means': compare_binomial_sample_means,\n",
    "    'means_chebyshev': compare_binomial_stderr_chebyshev,\n",
    "    'means_clt': compare_binomial_clt,\n",
    "    'max_lk': compare_binomial_maxlikelihood,\n",
    "    #'bootstrap': compare_groups_bootstrap_binomial,\n",
    "    'bayes': compare_groups_bayes_binomial,\n",
    "}\n",
    "\n",
    "df = pd.DataFrame(columns=['name', 'n_points', 'experiments', \n",
    "                           'correct', 'not_sure', 'incorrect'])\n",
    "\n",
    "for n in n_points:\n",
    "    exps = gen_binomial_experiments(n, k_experiments, pa_base)\n",
    "    cmp_exact = compare_experiments_binomial_exact(exps)\n",
    "    for k, v in comparisons.items():\n",
    "        cmp = compare_experiments(exps, method_name=k, method_fn=v)\n",
    "        stat = method_guesses(method=cmp, exact=cmp_exact, n_points=n)\n",
    "        df = pd.concat([df, pd.DataFrame([stat])], ignore_index=True)\n",
    "\n",
    "df['correct_and_not_sure'] = df['correct'] + df['not_sure']\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#todo: make means and max_lk distinguishable\n",
    "#df['line_dash'] = 0\n",
    "#df.loc[df['name'] == 'max_lk', 'line_dash'] = 1\n",
    "#https://stackoverflow.com/questions/64371174/how-to-change-variable-label-names-for-the-legend-in-a-plotly-express-line-chart\n",
    "\n",
    "fig = px.line(df, x='n_points', y='correct', color='name', markers='markers', log_x=True)\n",
    "fig.update_layout(\n",
    "    title='Correct',\n",
    "    yaxis_title='N Experiments',\n",
    "    xaxis_title='N Points',\n",
    "    legend_title='Method',\n",
    "    yaxis_range=[0, 130],\n",
    "    height=450, width=800\n",
    ")\n",
    "fig.show()\n",
    "\n",
    "fig = px.line(df, x='n_points', y='correct_and_not_sure', color='name', markers='markers', log_x=True)\n",
    "fig.update_layout(\n",
    "    title='Correct and Not Sure',\n",
    "    yaxis_title='N Experiments',\n",
    "    xaxis_title='N Points',\n",
    "    legend_title='Method',\n",
    "    yaxis_range=[0, 130],\n",
    "    height=450, width=800\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*При большом количестве точек решения бутстрапа и байесовского моделирования совпадают.*  \n",
    "*Бета-распределение и биномиальное распределение переходят в нормальное *\n",
    "$$\n",
    "N(p, \\frac{p(1-p)}{n}) ?\n",
    "$$\n",
    "\n",
    "Бета:  https://en.wikipedia.org/wiki/Beta_distribution#Special_and_limiting_cases  \n",
    "Биномиальное: https://en.wikipedia.org/wiki/Binomial_distribution#Normal_approximation  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ссылки\n",
    "\n",
    "[[Bern](https://en.wikipedia.org/wiki/Bernoulli_process)] - Bernoulli process, *in Wikipedia*      \n",
    "[[LLN](https://en.wikipedia.org/wiki/Law_of_large_numbers)] - Law of large numbers, *in Wikipedia*     \n",
    "[[StErrMn](https://en.wikipedia.org/wiki/Standard_error#Derivation)] - Standard error of the mean, *in Wikipedia*  \n",
    "[[StdUnbiased](https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation)] - Unbiased estimation of standard deviation, *in Wikipedia*  \n",
    "[[Cheb](https://en.wikipedia.org/wiki/Chebyshev%27s_inequality)] - Chebyshev's inequality, *in Wikipedia*     \n",
    "[[CLT](https://en.wikipedia.org/wiki/Central_limit_theorem)] - Central limit theorem, *in Wikipedia*   \n",
    "[[DMLT](https://en.wikipedia.org/wiki/De_Moivre%E2%80%93Laplace_theorem)] - De Moivre-Laplace theorem, *in Wikipedia*   \n",
    "[[ExpProp](https://en.wikipedia.org/wiki/Expected_value#Properties)] - Expected value, Properties *in Wikipedia*   \n",
    "[[VarProp](https://en.wikipedia.org/wiki/Standard_deviation#Identities_and_mathematical_properties)] - Standard deviation, Identities and mathematical properties *in Wikipedia*  \n",
    "[[Cov](https://en.wikipedia.org/wiki/Covariance)] - Covariance, *in Wikipedia*  \n",
    "[[ConfInt](https://en.wikipedia.org/wiki/Confidence_interval)] - Confidence interval, *in Wikipedia*   \n",
    "[[NormSum](https://en.wikipedia.org/wiki/Sum_of_normally_distributed_random_variables)] - Sum of normally distributed random variables, *in Wikipedia*   \n",
    "[[MxLk](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation)] - Maximum likelihood estimation, *in Wikipedia*  \n",
    "[[MAPE](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation)] - Maximum a posteriori estimation, *in Wikipedia*  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Бутстрап"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Еще один вариант закона больших чисел говорит, что по мере набора данных частоты возможных значений в сэмле будут приближаться к их вероятности. Это является основой для частотной интерпретации вероятности."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Также на этом основании имеющийся сэмпл данных можно считать приближением к точному распределению."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Общее название - методы повторной выборки данных или ресэмплинг.  \n",
    "Один из методов - метод складного ножа (Jackknife).   \n",
    "Еще один популярный метод - бутстрап.  \n",
    "Далее обсуждается только последний."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Бутстрап:  \n",
    "https://projecteuclid.org/journals/annals-of-statistics/volume-7/issue-1/Bootstrap-Methods-Another-Look-at-the-Jackknife/10.1214/aos/1176344552.full"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*What Teachers Should Know about the Bootstrap: Resampling in the Undergraduate Statistics Curriculum*  \n",
    "https://arxiv.org/abs/1411.5279  \n",
    "https://arxiv.org/pdf/1411.5279.pdf  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Можно оценить неопределенность средних значений."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_means_rawdata(dt, k_bootstrap):\n",
    "    bs = np.random.choice(dt, size=(k_bootstrap, n), replace=True)\n",
    "    means = bs.mean(axis=1)\n",
    "    return means\n",
    "\n",
    "def bootstrap_means_binomial(s, n, k_bootstrap):\n",
    "    p = s/n\n",
    "    bs = stats.binom.rvs(n, p, size=k_bootstrap)\n",
    "    means = bs / n\n",
    "    return means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_exact_a = 0.7\n",
    "p_exact_b = 0.65\n",
    "\n",
    "n = 5000\n",
    "s_a = stats.binom.rvs(n, p_exact_a)\n",
    "s_b = stats.binom.rvs(n, p_exact_b)\n",
    "p_sample_a = s_a/n\n",
    "p_sample_b = s_b/n\n",
    "\n",
    "k_bootstrap = 30000\n",
    "means_a = bootstrap_means_binomial(s_a, n, k_bootstrap)\n",
    "means_b = bootstrap_means_binomial(s_b, n, k_bootstrap)\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Histogram(x=means_a, histnorm='percent',\n",
    "                           name='A',\n",
    "                           opacity=0.3))\n",
    "#fig.add_vline(x=p_sample_a, line_dash='dash')\n",
    "fig.add_vline(x=p_exact_a)\n",
    "fig.add_trace(go.Histogram(x=means_b, histnorm='percent',\n",
    "                           name='B',\n",
    "                           opacity=0.3))\n",
    "#fig.add_vline(x=p_sample_b, line_dash='dash')\n",
    "fig.add_vline(x=p_exact_b)\n",
    "# fig.add_trace(go.Scatter(x=[p_exact_b, p_exact_b], y=[0, 1], \n",
    "#                          mode='lines', line_color='black', line_dash='dash',\n",
    "#                          name='Exact Mean B'))\n",
    "fig.update_layout(title='Means Bootstrap Dists',\n",
    "                  xaxis_title='p',\n",
    "                  yaxis_title='Prob',\n",
    "                  hovermode=\"x\",\n",
    "                  height=550,\n",
    "                  barmode='overlay')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_groups_bootstrap_binomial(sa, na, sb, nb):\n",
    "    k_bootstrap = 30000\n",
    "    means_a = bootstrap_means_binomial(sa, na, k_bootstrap)\n",
    "    means_b = bootstrap_means_binomial(sb, nb, k_bootstrap)\n",
    "    p_ea_gt_eb = np.sum(means_a > means_b) / len(means_a)\n",
    "    p_level = 0.95\n",
    "    res = None\n",
    "    if p_ea_gt_eb >= p_level:\n",
    "        res = 'A'\n",
    "    elif (1 - p_ea_gt_eb) >= p_level:\n",
    "        res = 'B'\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_a_exact = 0.7\n",
    "p_b_exact = 0.67\n",
    "\n",
    "n = 1000\n",
    "sa = stats.binom.rvs(n, p_a_exact)\n",
    "sb = stats.binom.rvs(n, p_b_exact)\n",
    "\n",
    "print(f'pa_exact={p_a_exact}, pb_exact={p_b_exact}')\n",
    "print(f'n={n}, sa={sa}, sb={sb}')\n",
    "print(f'Exact best group: {\"A\" if p_a_exact > p_b_exact else \"B\" if p_b_exact > p_a_exact else None}')\n",
    "print(f'Group selected by bootstrap comparison: {compare_groups_bootstrap_binomial(sa, n, sb, n)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вычислительно может быть затратно.  \n",
    "Есть пара трюков для ускорения вычислений.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пуассоновский бутстрап:  \n",
    "https://www.unofficialgoogledatascience.com/2015/08/an-introduction-to-poisson-bootstrap26.html  \n",
    "\n",
    "Сэмплировать из выборки [x1, x2, ..., xn] с повторением все равно, что генерировать набор\n",
    "мультиномиальных коэффициетов x ~ Multinomial(n, 1/n). Для больших n предлагается заменить Multinomial(n, 1/n) на n сэмплов из биномиального распределение Binom(n, 1/n). Количество точек в каждом векторе при этом перестает быть одинаковым. Еще один шаг - заменить биномиальное распределение Binom(n, 1/n) на распределение Пуассона Poisson(1). Сэмплировать из него. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Байесовский бутстрап:  \n",
    "https://gdmarmerola.github.io/the-bayesian-bootstrap/ - примеры  \n",
    "https://projecteuclid.org/journalArticle/Download?urlId=10.1214%2Faos%2F1176345338 - оригинальный текст  \n",
    "https://www.sumsar.net/blog/2015/04/the-non-parametric-bootstrap-as-a-bayesian-model/  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть данные x_1, ..., x_n.  \n",
    "Сгенерировать (n-1) точку u_1, ..., u_{n-1} из равномерного распределениея U(0,1);\n",
    "Отсортировать их;  \n",
    "Дополнить точками [0, ... , 1]  \n",
    "Посчитать разности g_i = u_i - u_i-1, i>=1.  \n",
    "Использовать g_i как вероятности для [x1, ..., x_n].   \n",
    "Среднее  \n",
    "m = sum g_i x_i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Есть отличия в интерпретации от обычного бутстрапа."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ограничения?  \n",
    "Бутстрап несколько хуже работает для непрерывных величин.  \n",
    "При малых размерах исходных данных.  \n",
    "Если часть данных мало представлена (исходное распределение скошено, в выборке мало значений из хвоста).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}