sem11_AB/AB_handmade_BMM-192.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AB-тесты своими руками\n",
    "\n",
    "Проделаем историю с Дядей Фёдором руками!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_f = np.array([150, 190])\n",
    "x_s = np.array([180, 190])\n",
    "x_m = np.array([150, 160])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "170.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_f.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_f.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28.284271247461902"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_f.std(ddof=1) # ddof=1 тк делим на n - 1, если на n, тогда ddof=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.6448536269514729"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "\n",
    "# так можно найти квантиль - то есть число левее которого alpha% выборки\n",
    "stats.norm().ppf(0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.959963984540054"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# хотим чтобы с каждой сторону купола было по 2.5% (если хотим 95% интервал)\n",
    "# Уточнили правило 2-х сигм! Чтобы было ровно 95%, надо брать не 2, а 1.96\n",
    "stats.norm().ppf(1 - 0.05/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Доверительный интервал для среднего своими руками:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Доверительный интервал: [130.80072030919894;209.19927969080106])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "78.39855938160213"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def conf_int(x, alpha=0.05):\n",
    "    mn = x.mean()                       # нашли среднее\n",
    "    st = x.std(ddof=1)/np.sqrt(x.size)  # нашли стандартное отклонение \n",
    "    \n",
    "    # нужный нам квантиль, в правиле 3-х сигм, q=3\n",
    "    q = stats.norm().ppf(1 - alpha/2)\n",
    "    \n",
    "    left = mn - q * st\n",
    "    right = mn + q * st\n",
    "    \n",
    "    print(f\"Доверительный интервал: [{left};{right}])\")\n",
    "    return left, right\n",
    "\n",
    "l,r = conf_int(x_f)\n",
    "r - l # длина интервала"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пакетный способ сделать то же самое."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(130.80072030919894, 209.19927969080106)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.norm.interval(0.95, x_f.mean(), x_f.std(ddof=1)/np.sqrt(x_f.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь поговорим про гипотезы! Проверим гипотезу, что $\\mu = 160$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-29.199279690801077, 49.19927969080108)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# alpha = 0.05 \n",
    "\n",
    "# H_0: mu = 160\n",
    "# H_a: mu != 160\n",
    "\n",
    "stats.norm.interval(0.95, x_f.mean() - 160, x_f.std(ddof=1)/np.sqrt(x_f.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что ноль попал в доверительный интервал $\\Rightarrow$ гипотеза не отвергается. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Офирмим это в виде функции, чтобы считалось через $t$-статистику."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_norm_test(x, mu0, alpha=0.05):\n",
    "    mn = x.mean() - mu0\n",
    "    st = x.std(ddof=1)/np.sqrt(x.size)\n",
    "    \n",
    "    t = mn/st # нашли наблюдаемое значение t-статистики\n",
    "    \n",
    "    # критическое значение t-статистики \n",
    "    q = stats.norm().ppf(1 - alpha/2)\n",
    "    \n",
    "    if (t < q) & (t > -q): # np.abs(t) < q\n",
    "        return \"Гипотеза не отвергается\"\n",
    "    else:\n",
    "        return \"Гипотеза отвергается\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Гипотеза не отвергается'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_norm_test(x_f, mu0=160, alpha=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus-track для целевой аудитории семинара в Zoom :) \n",
    "\n",
    "Поговорим про ошибку первого рода!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.46695432,  0.07103873,  0.28820446,  0.75078223,  1.11529013,\n",
       "        1.54429739,  0.68892519, -1.8845052 ,  2.33624069,  2.37726645])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_rv = stats.norm( ) # инициализровали нормальное распределение N(0, 1)\n",
    "                        # нулевое мат ожидание, единичное стандартное отклонение\n",
    "\n",
    "norm_rv.rvs(10) # сгенерировали 10 наблюдений "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6448536269514722"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_rv.ppf(0.95) # можно найти квантиль уровня 95%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(-1.859963984540054, 0.4, '-1.96')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = np.linspace(-5, 5, 300)\n",
    "pdf = norm_rv.pdf(x) # плотность распределения вероятностей\n",
    "\n",
    "alpha = 0.05 \n",
    "q = norm_rv.ppf(1 - alpha/2)\n",
    "\n",
    "plt.plot(x, pdf);\n",
    "plt.vlines(q, 0, 0.4)\n",
    "plt.text(q + 0.1, 0.4, f'{q:.3}')\n",
    "\n",
    "plt.vlines(-q, 0, 0.4)\n",
    "plt.text(-q + 0.1, 0.4, f'{-q:.3}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Уровень значимости (ошибка первого рода) - вероятность отвергнуть гипотезу $H_0$, когда она верна. Попробуем с помощью симуляций вычислить её для проверки гипотезы о среднем."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_norm_test(x, mu0, alpha=0.05):\n",
    "    mn = x.mean() - mu0\n",
    "    st = x.std(ddof=1)/np.sqrt(x.size)\n",
    "    \n",
    "    t = mn/st # нашли наблюдаемое значение t-статистики\n",
    "    \n",
    "    # критическое значение t-статистики \n",
    "    q = stats.norm().ppf(1 - alpha/2)\n",
    "    \n",
    "    if (t < q) & (t > -q): # np.abs(t) < q\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Генерируем данные из нормального распределения с математическим ожиданием (`loc=4`). И проверяем по сгенерированной выборке гипотезу, что $\\mu = 4$. В $5\\%$ случаев мы почём-зря отвергнем эту гипотезу."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_obs = 10**4\n",
    "itog = np.zeros(n_obs)\n",
    "\n",
    "for i in range(n_obs):\n",
    "    # задали генератор\n",
    "    norm_rv = stats.norm(loc=4, scale=1)\n",
    "    \n",
    "    # сгенерировали выборку объёма 100\n",
    "    x = norm_rv.rvs(100)\n",
    "    itog[i] = mean_norm_test(x, mu0=4, alpha=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9486"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "itog.mean() # как часто не отвергали (95%)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}