07. Гипотеза и вывод.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_0 500.0\n",
      "sigma_0 15.811388300841896\n",
      "normal_two_sided_bounds(0.95, mu_0, sigma_0) (469.01026640487555, 530.9897335951244)\n",
      "\n",
      "Мощность критерия\n",
      "95% границы на основе допущения, что p = 0.5\n",
      "lo 469.01026640487555\n",
      "hi 530.9897335951244\n",
      "фактические mu и sigma на основе p = 0.55\n",
      "mu_1 550.0\n",
      "sigma_1 15.732132722552274\n",
      "вероятность ошибки 2-ого рода 0.11345199870463285\n",
      "мощность 0.8865480012953671\n",
      "Односторонний критерий\n",
      "hi 526.0073585242053\n",
      "вероятность ошибки 2-ого рода 0.06362051966928273\n",
      "мощность 0.9363794803307173\n",
      "\n",
      "two_sided_p_value(529.5, mu_0, sigma_0) 0.06207721579598857\n",
      "two_sided_p_value(531.5, mu_0, sigma_0) 0.046345287837786575\n",
      "upper_p_value(525, mu_0, sigma_0) 0.056923149003329065\n",
      "upper_p_value(527, mu_0, sigma_0) 0.04385251499101195\n",
      "\n",
      "Подгонка p-значения\n",
      "46 отклонений из 1000\n",
      "\n",
      "A/B-тестирование\n",
      "a_b_test_statistic(1000, 200, 1000, 180) -1.1403464899034472\n",
      "p-значение 0.254141976542236\n",
      "a_b_test_statistic(1000, 200, 1000, 150) -2.948839123097944\n",
      "p-значение 0.003189699706216853\n"
     ]
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "# hypothesis_and_inference.py\n",
    "\n",
    "import sys\n",
    "sys.path.append(\"../code-python3-ru\")\n",
    "\n",
    "from lib.probability import normal_cdf, inverse_normal_cdf\n",
    "import math, random\n",
    "\n",
    "def normal_approximation_to_binomial(n, p):\n",
    "    \"\"\"находит mu и sigma, которые соответствуют binomial(n, p)\"\"\"\n",
    "    mu = p * n\n",
    "    sigma = math.sqrt(p * (1 - p) * n)\n",
    "    return mu, sigma\n",
    "\n",
    "#####\n",
    "#\n",
    "# вероятности, что нормальная случайная величина лежит в интервале\n",
    "#\n",
    "######\n",
    "\n",
    "# вероятность, что реализованное значение нормальной величины лежит ниже\n",
    "normal_probability_below = normal_cdf\n",
    "\n",
    "# вероятность, что оно лежит выше, если оно не ниже его\n",
    "def normal_probability_above(lo, mu=0, sigma=1):\n",
    "    return 1 - normal_cdf(lo, mu, sigma)\n",
    "\n",
    "# вероятность, что оно лежит между, если оно меньше hi, но не ниже lo\n",
    "def normal_probability_between(lo, hi, mu=0, sigma=1):\n",
    "    return normal_cdf(hi, mu, sigma) - normal_cdf(lo, mu, sigma)\n",
    "\n",
    "# вероятность, что оно лежит за пределами, если оно не внутри\n",
    "def normal_probability_outside(lo, hi, mu=0, sigma=1):\n",
    "    return 1 - normal_probability_between(lo, hi, mu, sigma)\n",
    "\n",
    "######\n",
    "#\n",
    "#  границы нормальной величины\n",
    "#\n",
    "######\n",
    "\n",
    "def normal_upper_bound(probability, mu=0, sigma=1):\n",
    "    \"\"\"возвращает z, для которого P(Z <= z) = probability\"\"\"\n",
    "    return inverse_normal_cdf(probability, mu, sigma)\n",
    "\n",
    "def normal_lower_bound(probability, mu=0, sigma=1):\n",
    "    \"\"\"возвращает z, для которого P(Z >= z) = probability\"\"\"\n",
    "    return inverse_normal_cdf(1 - probability, mu, sigma)\n",
    "\n",
    "def normal_two_sided_bounds(probability, mu=0, sigma=1):\n",
    "    \"\"\"возвращает симметричные (вокруг среднего значения) границы,\n",
    "    в пределах которых содержится указанная вероятность\"\"\"\n",
    "    tail_probability = (1 - probability) / 2\n",
    "\n",
    "    # верхняя граница должна иметь значение хвостовой вероятности\n",
    "    upper_bound = normal_lower_bound(tail_probability, mu, sigma)\n",
    "\n",
    "    # нижняя граница должна иметь значение хвостовой вероятности\n",
    "    lower_bound = normal_upper_bound(tail_probability, mu, sigma)\n",
    "\n",
    "    return lower_bound, upper_bound\n",
    "\n",
    "def two_sided_p_value(x, mu=0, sigma=1):\n",
    "    if x >= mu:\n",
    "        # если x больше среднего значения, то значения в хвосте больше x\n",
    "        return 2 * normal_probability_above(x, mu, sigma)\n",
    "    else:\n",
    "        # если x меньше среднего значения, то значения в хвосте меньше x\n",
    "        return 2 * normal_probability_below(x, mu, sigma)\n",
    "\n",
    "def count_extreme_values():\n",
    "    extreme_value_count = 0\n",
    "    for _ in range(100000):\n",
    "        num_heads = sum(1 if random.random() < 0.5 else 0    # подсчитать число орлов\n",
    "                        for _ in range(1000))                # при 1000 бросках\n",
    "        if num_heads >= 530 or num_heads <= 470:             # подсчитать, как часто\n",
    "            extreme_value_count += 1                         # число - 'предельно'\n",
    "\n",
    "    return extreme_value_count / 100000\n",
    "\n",
    "upper_p_value = normal_probability_above\n",
    "lower_p_value = normal_probability_below\n",
    "\n",
    "##\n",
    "#\n",
    "# Подгонка p-значения\n",
    "#\n",
    "##\n",
    "\n",
    "def run_experiment():\n",
    "    \"\"\"бросить уравновешенную монету 1000 раз, True = орлы, False = решки\"\"\"\n",
    "    return [random.random() < 0.5 for _ in range(1000)]\n",
    "\n",
    "def reject_fairness(experiment):\n",
    "    \"\"\"используя 5%-ые уровни значимости\"\"\"\n",
    "    num_heads = len([flip for flip in experiment if flip])\n",
    "    return num_heads < 469 or num_heads > 531\n",
    "\n",
    "##\n",
    "#\n",
    "# проведение A/B-тестирования\n",
    "#\n",
    "##\n",
    "\n",
    "def estimated_parameters(N, n):\n",
    "    p = n / N\n",
    "    sigma = math.sqrt(p * (1 - p) / N)\n",
    "    return p, sigma\n",
    "\n",
    "def a_b_test_statistic(N_A, n_A, N_B, n_B):\n",
    "    p_A, sigma_A = estimated_parameters(N_A, n_A)\n",
    "    p_B, sigma_B = estimated_parameters(N_B, n_B)\n",
    "    return (p_B - p_A) / math.sqrt(sigma_A ** 2 + sigma_B ** 2)\n",
    "\n",
    "##\n",
    "#\n",
    "# Байесовский статистический вывод\n",
    "#\n",
    "##\n",
    "\n",
    "def B(alpha, beta):\n",
    "    \"\"\"нормализующая константа, благодаря которой сумма вероятностей равна 1\"\"\"\n",
    "    return math.gamma(alpha) * math.gamma(beta) / math.gamma(alpha + beta)\n",
    "\n",
    "def beta_pdf(x, alpha, beta):\n",
    "    if x < 0 or x > 1:          # за пределами [0, 1] нет веса\n",
    "        return 0\n",
    "    return x ** (alpha - 1) * (1 - x) ** (beta - 1) / B(alpha, beta)\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "\n",
    "    mu_0, sigma_0 = normal_approximation_to_binomial(1000, 0.5)\n",
    "    print(\"mu_0\", mu_0)\n",
    "    print(\"sigma_0\", sigma_0)\n",
    "    print(\"normal_two_sided_bounds(0.95, mu_0, sigma_0)\", normal_two_sided_bounds(0.95, mu_0, sigma_0))\n",
    "    print()\n",
    "    print(\"Мощность критерия\")\n",
    "\n",
    "    print(\"95% границы на основе допущения, что p = 0.5\")\n",
    "\n",
    "    lo, hi = normal_two_sided_bounds(0.95, mu_0, sigma_0)\n",
    "    print(\"lo\", lo)\n",
    "    print(\"hi\", hi)\n",
    "\n",
    "    print(\"фактические mu и sigma на основе p = 0.55\")\n",
    "    mu_1, sigma_1 = normal_approximation_to_binomial(1000, 0.55)\n",
    "    print(\"mu_1\", mu_1)\n",
    "    print(\"sigma_1\", sigma_1)\n",
    "\n",
    "    # ошибка 2-ого рода означает, что нам не удастся отклонить нулевую гипотезу,\n",
    "    # которая случится, когда X все еще лежит в нашем исходном интервале\n",
    "    type_2_probability = normal_probability_between(lo, hi, mu_1, sigma_1)\n",
    "    power = 1 - type_2_probability # 0.887\n",
    "\n",
    "    print(\"вероятность ошибки 2-ого рода\", type_2_probability)\n",
    "    print(\"мощность\", power)\n",
    "    print\n",
    "\n",
    "    print(\"Односторонний критерий\")\n",
    "    hi = normal_upper_bound(0.95, mu_0, sigma_0)\n",
    "    print(\"hi\", hi) # равно 526 (< 531, т.к. нам нужно больше вероятности в верхнем хвосте)\n",
    "    type_2_probability = normal_probability_below(hi, mu_1, sigma_1)\n",
    "    power = 1 - type_2_probability # = 0.936\n",
    "    print(\"вероятность ошибки 2-ого рода\", type_2_probability)\n",
    "    print(\"мощность\", power)\n",
    "    print()\n",
    "\n",
    "    print(\"two_sided_p_value(529.5, mu_0, sigma_0)\", two_sided_p_value(529.5, mu_0, sigma_0))\n",
    "\n",
    "    print(\"two_sided_p_value(531.5, mu_0, sigma_0)\", two_sided_p_value(531.5, mu_0, sigma_0))\n",
    "\n",
    "    print(\"upper_p_value(525, mu_0, sigma_0)\", upper_p_value(525, mu_0, sigma_0))\n",
    "    print(\"upper_p_value(527, mu_0, sigma_0)\", upper_p_value(527, mu_0, sigma_0))\n",
    "    print()\n",
    "\n",
    "    print(\"Подгонка p-значения\")\n",
    "\n",
    "    random.seed(0)\n",
    "    experiments = [run_experiment() for _ in range(1000)]\n",
    "    num_rejections = len([experiment\n",
    "                          for experiment in experiments\n",
    "                          if reject_fairness(experiment)])\n",
    "\n",
    "    print(num_rejections, \"отклонений из 1000\")\n",
    "    print()\n",
    "\n",
    "    print(\"A/B-тестирование\")\n",
    "    z = a_b_test_statistic(1000, 200, 1000, 180)\n",
    "    print(\"a_b_test_statistic(1000, 200, 1000, 180)\", z)\n",
    "    print(\"p-значение\", two_sided_p_value(z))\n",
    "    z = a_b_test_statistic(1000, 200, 1000, 150)\n",
    "    print(\"a_b_test_statistic(1000, 200, 1000, 150)\", z)\n",
    "    print(\"p-значение\", two_sided_p_value(z))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}