Last active
January 12, 2022 12:05
-
-
Save jamm1985/eb98dfe4dd8b332d91238b69004cc417 to your computer and use it in GitHub Desktop.
Lab_10_hypotesis_testing_part_III.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "nbformat": 4, | |
| "nbformat_minor": 0, | |
| "metadata": { | |
| "colab": { | |
| "name": "Lab_10_hypotesis_testing_part_III.ipynb", | |
| "provenance": [], | |
| "authorship_tag": "ABX9TyOJjDOrCIh/snDsToXHLiZy", | |
| "include_colab_link": true | |
| }, | |
| "kernelspec": { | |
| "name": "python3", | |
| "display_name": "Python 3" | |
| }, | |
| "language_info": { | |
| "name": "python" | |
| } | |
| }, | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "id": "view-in-github", | |
| "colab_type": "text" | |
| }, | |
| "source": [ | |
| "<a href=\"https://colab.research.google.com/gist/jamm1985/eb98dfe4dd8b332d91238b69004cc417/lab_10_hypotesis_testing_part_iii.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "Видео лабораторной: https://youtu.be/wLw3GZVXJ6M\n", | |
| "\n", | |
| "TG: https://t.me/data_science_news\n", | |
| "\n", | |
| "\n", | |
| "\n", | |
| "---\n", | |
| "\n" | |
| ], | |
| "metadata": { | |
| "id": "MfhwDmL8CU0R" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "id": "XSjs2khbGsgY" | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import numpy as np\n", | |
| "import matplotlib.pylab as plt\n", | |
| "from scipy.stats import norm, binom\n", | |
| "from scipy.stats import fisher_exact" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "# Тетсирование биномиальных данных" | |
| ], | |
| "metadata": { | |
| "id": "HlEr4p6KXbtm" | |
| } | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "## Задача №1 _(асимтотический двусторонний тест на одной выборке)_\n", | |
| "\n", | |
| "В шахматном онлайн сервисе заявлено, что цвет фигур (черные или белые) между двумя игнроками выбирается случайно и независимо в свободныых партиях с вероятностью 50/50. Вы проанализировали 130 собственных партий и выяснили, что 69 раз вы играли черными фигурами. Так случайно вышло, что вы играли черными фигурами больше или существует какой то другой фактор влияющий на выбор цвета фигур?" | |
| ], | |
| "metadata": { | |
| "id": "f1OXr_bfX03S" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "N=130\n", | |
| "black=69\n", | |
| "p_0=0.5" | |
| ], | |
| "metadata": { | |
| "id": "dwR4bxUd1Kx-" | |
| }, | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "## Формализация задачи\n", | |
| "\n", | |
| "Распределение данных по задаче можно выразить как последовательность независимых случайных величин Бернулли. Предположим, что $X_i=1$ если цвет фигур черный, и $X_i=0$ если белый. Таким образом, количество игр за черных в $N$ независимых итерациях выражено суммой $Y = \\Sigma_{i=1}^N X_i$, которое подчиняется биномиальному распределению:\n", | |
| "\n", | |
| "$p_Y(k;p)=P(Y=k)=\\binom{n}{k}p^k(1-p)^{n-k}$, $k=0,1,...,N$\n", | |
| "\n", | |
| "$\\mathbb{E}[Y]=np$\n", | |
| "\n", | |
| "$Var[Y]=np(1-p)$\n", | |
| "\n", | |
| "По условию задачи мы должны протестировать гипотезу:\n", | |
| "\n", | |
| "$H_0:\\ p=0.5\\ \\mathrm{VS}\\ H_A: p \\neq 0.5$, где $p$ это параметр биномиального распределения." | |
| ], | |
| "metadata": { | |
| "id": "Iw-RBRsFgpW3" | |
| } | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Асимтотичность - практическое правило\n", | |
| "\n", | |
| "\n", | |
| "$0 < np_0 - 3 \\sqrt{np_0(1-p_0)} < np_0 + 3 \\sqrt{np_0(1-p_0)} < n$ \n", | |
| "\n", | |
| "_[Larsen, Marx, An introduction to mathmatical statistics and its application 5-th ED p.361]_\n", | |
| "\n", | |
| "Если неравенство выполняется, то мы можем выполнить асимтотический тест.\n", | |
| "\n", | |
| "При этом у нас есть два пути для выбора тест статистики:\n", | |
| "\n", | |
| "1. Рассмотривать сумму $\\Sigma_{i=1}^N X_i \\sim N(np, np(1-p))$. Таким образом, нормализированная статистика будет иметь вид $\\frac{\\Sigma_{i=1}^N X_i-np}{\\sqrt{np(1-p)}} \\sim N(0,1)$\n", | |
| "\n", | |
| "2. Расматривать выборочное среднее $\\bar{X}_n = \\frac{1}{n}\\Sigma_{i=1}^N X_i \\sim N(p,\\frac{p(1-p)}{n})$, где нормализированный вариант будет иметь вид $\\sqrt{N}\\frac{\\bar{X}_N - p}{p(1-p)} \\sim N(0,1)$" | |
| ], | |
| "metadata": { | |
| "id": "ggpiROiamn9f" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "print(N*p_0-3*np.sqrt(N)*p_0)\n", | |
| "print(N*p_0+3*np.sqrt(N)*p_0)\n", | |
| "print(0<N*p_0-3*np.sqrt(N)*p_0<N*p_0+3*np.sqrt(N)*p_0<N)" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "bT_IP2vlG7zM", | |
| "outputId": "9017a1ac-35c8-4c68-c9dd-639ebd6a10b5" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "47.89736862351293\n", | |
| "82.10263137648707\n", | |
| "True\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Статистический критерий\n", | |
| "\n", | |
| "Зафиксируем уровень статистической значимости $\\alpha=0.05$. То есть вероятность того, что мы отклоним нулевую гипотезу, при том что она истинна равна 5% или $P(\\mathrm{reject\\ H_0}|\\mathrm{H_0\\ true})=0.05$.\n", | |
| "\n", | |
| "В условиях нулевой гипотезы тест статистика имеет вид:\n", | |
| "\n", | |
| "$$\\frac{\\Sigma_{i=1}^N X_i-np_0}{\\sqrt{np_0(1-p_0)}}=\\sqrt{N}\\frac{\\bar{X}_N - p_0}{\\sqrt{p_0(1-p_0)}} \\sim N(0,1)$$\n", | |
| "\n", | |
| "Таким образом, статистический критерий выражен:\n", | |
| "\n", | |
| "$$\\psi_\\alpha=\\mathbb{1}\\left[\\frac{|\\Sigma_{i=1}^N X_i-np_0|}{\\sqrt{np_0(1-p_0)}}=\\sqrt{N}\\frac{|\\bar{X}_N - p_0|}{\\sqrt{p_0(1-p_0)}} > q_{\\frac{\\alpha}{2}}=1.96 \\right] $$" | |
| ], | |
| "metadata": { | |
| "id": "zxfEzCIn13BV" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# alpha/2 quantile\n", | |
| "print(np.round(norm.ppf(0.975), decimals=2))\n", | |
| "\n", | |
| "# test statistics 1 \n", | |
| "test_stat_1 = np.abs(black-N*p_0)/np.sqrt(N*p_0*(1-p_0))\n", | |
| "print(\n", | |
| " np.round(test_stat_1, decimals=3))\n", | |
| "\n", | |
| "# test statistics 2\n", | |
| "test_stat_2 = np.sqrt(N)*np.abs(black/N-p_0)/np.sqrt(p_0*(1-p_0))\n", | |
| "print(\n", | |
| " np.round(test_stat_2, decimals=3))\n", | |
| "\n", | |
| "# decision\n", | |
| "print(\n", | |
| " 'Is the null hypothesis rejected? / {}'.format(\n", | |
| " test_stat_1 > norm.ppf(0.975)\n", | |
| " )\n", | |
| " )\n", | |
| "\n", | |
| "# p-value\n", | |
| "print(np.round(1-norm.cdf(test_stat_1), decimals=2))" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "3N_gF1Zd1hVG", | |
| "outputId": "562572bf-3f8e-4ac9-c14b-9f813a8b2fc8" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "1.96\n", | |
| "0.702\n", | |
| "0.702\n", | |
| "Is the null hypothesis rejected? / False\n", | |
| "0.24\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "## Задача №2 _(точный двусторонний тест на одной выбоке)_\n", | |
| "\n", | |
| "Предположим, что 21 пациент с редким заболеванием участвуют в исследовании лекартсва, которое изменяет болевой синдром. Известно, что действующий препарат эффективен в 80% случаев. Допустим у 19 пациентов уменьшился болевой синдром. Имеет ли новый препарат эффективность перед существующим перпаратом? \n" | |
| ], | |
| "metadata": { | |
| "id": "S6bY7UOTzdRy" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "p_0=0.8\n", | |
| "N=21\n", | |
| "alpha = 0.05" | |
| ], | |
| "metadata": { | |
| "id": "CGlBnUlQ1wAs" | |
| }, | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Асимтотичность - практическое правило\n", | |
| "\n", | |
| "\n", | |
| "$0 < np_0 - 3 \\sqrt{np_0(1-p_0)} < np_0 + 3 \\sqrt{np_0(1-p_0)} < n$ \n", | |
| "\n", | |
| "_[Larsen, Marx, An introduction to mathmatical statistics and its application 5-th ED p.361]_" | |
| ], | |
| "metadata": { | |
| "id": "B95Wzlk-I2W1" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "print(N*p_0-3*np.sqrt(N*p_0*(1-p_0)))\n", | |
| "print(N*p_0+3*np.sqrt(N*p_0*(1-p_0)))\n", | |
| "print(0<N*p_0-3*np.sqrt(N*p_0*(1-p_0))<N*p_0+3*np.sqrt(N*p_0*(1-p_0))<N)" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "FWqb1VRBIuvT", | |
| "outputId": "9f2600eb-2fed-42c7-bdc1-d959af747b40" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "11.300909166052993\n", | |
| "22.29909083394701\n", | |
| "False\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Формализация задачи\n", | |
| "\n", | |
| "Распределение данных по задаче можно выразить как последовательность независимых случайных величин Бернулли. Предположим, что $X_i=1$ если у пациента уменьшился болевой синдром, и $X_i=0$ пациент не ощущает изменений. Таким образом, количество пациентиов, которые ощутили эффект в $N$ независимых итерациях выражено суммой $Y = \\Sigma_{i=1}^N X_i$, которое в условиях нулевой гипотезы подчиняется биномиальному распределению:\n", | |
| "\n", | |
| "$p_Y(k;p)=P(Y=k)=\\binom{n}{k}p^k(1-p)^{n-k}=\\binom{21}{k}0.8^k(1-0.8)^{21-k}$, $k=0,1,...,N$\n", | |
| "\n", | |
| "\n", | |
| "По условию задачи мы должны протестировать гипотезу:\n", | |
| "\n", | |
| "$H_0:\\ p=0.8\\ \\mathrm{VS}\\ H_A: p \\neq 0.8$, где $p$ это параметр биномиального распределения." | |
| ], | |
| "metadata": { | |
| "id": "miZJM-gPIkGH" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "rv_binom = binom(N, p_0)\n", | |
| "x = np.arange(0,22)\n", | |
| "pmf_x = np.round(rv_binom.pmf(x), decimals=3)\n", | |
| "print(\"\\n\".join(\"{} {}\".format(x, y) for x, y in zip(x, pmf_x)))" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "dfei5tZZwsxS", | |
| "outputId": "408b529d-6a90-479c-e160-935b0da9e4de" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "0 0.0\n", | |
| "1 0.0\n", | |
| "2 0.0\n", | |
| "3 0.0\n", | |
| "4 0.0\n", | |
| "5 0.0\n", | |
| "6 0.0\n", | |
| "7 0.0\n", | |
| "8 0.0\n", | |
| "9 0.0\n", | |
| "10 0.001\n", | |
| "11 0.003\n", | |
| "12 0.01\n", | |
| "13 0.029\n", | |
| "14 0.065\n", | |
| "15 0.122\n", | |
| "16 0.183\n", | |
| "17 0.216\n", | |
| "18 0.192\n", | |
| "19 0.121\n", | |
| "20 0.048\n", | |
| "21 0.009\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "right_bound = rv_binom.ppf(1 - alpha/2)\n", | |
| "left_bound = rv_binom.ppf(alpha/2)\n", | |
| "print('left bound is {}, right bound is {}'.format(left_bound, right_bound))" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "xk1CJG0A23Wy", | |
| "outputId": "c58a9ad4-5ad2-43bf-ccae-205ec3d013c1" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "left bound is 13.0, right bound is 20.0\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "print('P(Y<13)={}'.format(np.round(rv_binom.cdf(13), decimals=3)))\n", | |
| "print('P(Y>20)={}'.format(np.round(1-rv_binom.cdf(20), decimals=3)))" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "zchRATVjNFg6", | |
| "outputId": "4ff8c3b3-b6db-408d-aa2f-358584b3e114" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "P(Y<13)=0.043\n", | |
| "P(Y>20)=0.009\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "plt.rcParams['figure.figsize'] = [12, 8]\n", | |
| "plt.plot(x+2, rv_binom.pmf(x+2))\n", | |
| "plt.vlines([left_bound, right_bound],-0.00,+0.23, 'red')\n", | |
| "plt.ylabel(r'PMF $Y$', fontsize = 15)\n", | |
| "plt.xlabel(r'$Y$', fontsize = 15)" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/", | |
| "height": 520 | |
| }, | |
| "id": "httEbGmx4Fia", | |
| "outputId": "80cd8ae9-7f93-401d-b822-6a579cdae854" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "Text(0.5, 0, '$Y$')" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 48 | |
| }, | |
| { | |
| "output_type": "display_data", | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAHmCAYAAACvYjIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXjcZ33u//uZ0WbtM5Jsy5Y8I9uyHTteEi1J7NgJCVnYEsqaQMgCNFAOpWXroYe2nNJfKQXKcoC2CWSFQkig0EADWewkcuwkthzbSRzLki1r8apdsnZp5vn9YTk1juN1Zp5Z3q/r0oU0mpHuP4h0+6vP9/MYa60AAAAARI/HdQAAAAAg2VG6AQAAgCijdAMAAABRRukGAAAAoozSDQAAAEQZpRsAAACIsjTXAaKtuLjYBoNB1zEAAACQ5LZs2dJlrS052eeSvnQHg0HV19e7jgEAAIAkZ4xpfbPPMV4CAAAARBmlGwAAAIgySjcAAAAQZZRuAAAAIMoo3QAAAECUUboBAACAKKN0AwAAAFFG6QYAAACijNINAAAARBmlGwAAAIgySjcAAAAQZZRuAAAAIMoo3QAAAECUUboBAACAKKN0AwAAAFFG6QYAAACijNINAAAARBmlGwAAAIgySjcAIHKuvPLoG4BT47+VlEPpBgAAAKKM0g0AAABEGaUbAAAAiDJKNwAAABBllG4AAAAgyijdAAAAQJRRugEAAIAoo3QDAAAAUUbpBgAAAKKM0g0AAABEGaUbAAAAiDJKNwAAABBllG4AAAAgyijdAAAAQJRRugEAAIAoo3QDAAAAUUbpBgAAAKKM0g0AAABEGaUbAAAAiDJKNwAAABBllG4AAAAgyijdAAAAQJRRugEAwEm9vK9PI+Mh1zGApEDpBgAAb/BSW69u+MEGfeahrbLWuo4DJDxKNwAAeINvPb5LXo/Rk68d1i82t7uOAyQ8SjcAAPgjG3Z3aeOebv312xZp1fwi/f1vX9PeriHXsYCERukGAACvs9bqm4/vUmlBlm65NKBvvX+5MtI8+stfbNNEKOw6HpCwKN0AAOB1T+3s0Lb2Pn3m6kplpXtVWjBN//Sepdre3qfvr9vtOh6QsCjdAABAkhQOW/3LE7sULMrW+6rKXn/87UtL9d6Ly/SDdU3a0trrMCGQuCjdAABAkvTblw+o4dARffaaBUr3/nFF+L83LNaswmn67C+2aXBs0lFCIHFRugEAgCZCYX3nyUYtmpmndy2b9YbP52Wl6zsfXKF9vcP66m93OEgIJDZKNwAA0K+27FNL97A+d80CeTzmpM+pCfr1qSvn6+H6ffrDqwdjnBBIbJRuAABS3OhESP9vbZOWlxfqmsUzTvncv3hrpZaVFehL//mKDg+MxighkPgo3QAApLifvdimA/2j+uK1C2XMya9yH5Pu9ei7H1yhsYmwvvDIdoXDnFYJnAlKNwAAKWxobFL/+sxuXTa3SKvmF53Ra+aW5Opv3nmB1jd16YHnW6KaD0gWlG4AAFLY/Rtb1DU4ri9cd/qr3Mf7UO0cXb1ouv7p9w1qPHwkigmB5EDpBgAgRfWPTOiuZ/fo6kXTVRXwndVrjTH65/ctU35Wmj7z860amwxFKSWQHCjdAACkqB/VNWtgdFKfu3bBOb2+ODdT//zeZWo4dET/8kRjhNMByYXSDQBACuoaHNO9G/bqnctKtWRWwTl/nasvmKFbLp2jH61v1sbdXRFMCCQXSjcAACnoX5/eo9GJkD57zbld5T7el9++WBXFOfr8I9vVPzwRgXRA8qF0AwCQYg70jeinL7TqfVVlmleSe95fb1qGV9/94Ap1HhnTl3/ziqxljSBwIko3AAAp5vvrmmRl9ZmrKyP2NZeVFeqz1yzQ714+qN9s2x+xrwskC0o3AAApZG/XkB6u36cPXxJQmS87ol/7k1fMU03Qp7/7zQ619wxH9GsDiY7SDQBACvnuU41K9xp96i3zIv61vR6jb39ghaykzz+8XSFOqwReR+kGACBFNBwa0KPbD+iOVRWanpcVle9R7s/WV29cok0tPbqrbk9UvgeQiCjdAACkiH95olG5GWn6xJq5Uf0+f3LRbL1jWam+/USjXt3fH9XvBSQKSjcAAClga1uvnnztsO5cM1eF2RlR/V7GGP3juy9UcW6mPvPQVo2Mc1ol4KR0G2OuN8bsMsbsNsZ86SSf/5wx5jVjzMvGmLXGmMBxn7vNGNM09XZbbJMDAJCY/uWJRvlzMnTH5RUx+X6F2Rn6lw8sV3PnkL722M6YfE8gnsW8dBtjvJJ+KOltkhZLutkYs/iEp22VVG2tXSbpl5K+MfVav6SvSLpEUq2krxhjfLHKDgBAItq4p0vP7e7Sp66cp9zMtJh931Xzi/Xxyyv0kxda9XRDR8y+LxCPXFzprpW021rbbK0dl/SQpBuPf4K19mlr7bFdQy9IKpt6/zpJT1pre6y1vZKelHR9jHIDAJBwrLX61uO7NDM/S7dcGjj9CyLsi9cv1KKZefriL7era3As5t8fiBcuSvdsSe3Hfbxv6rE38zFJvz+b1xpj7jTG1Btj6js7O88zLgAAiWtdQ4deauvTZ66uVFa6N+bfPzPNq+/etEIDo5P60q84rRKpK65vpDTG3CKpWtI3z+Z11tq7rbXV1trqkpKS6IQDACDOhcNW33qiUYGibL2/uuz0L4iSRTPz9b+vX6Sndh7Wzze1n/4FQBJyUbr3Syo/7uOyqcf+iDHmrZK+LOkGa+3Y2bwWAABI//3KQe08OKDPvnWB0r1ur7PdsTKoy+cX6x9+95qaOwedZgFccPFf4GZJlcaYCmNMhqSbJD16/BOMMRdJuktHC/fxd148LulaY4xv6gbKa6ceAwAAx5kMhfWdJxu1YEau3rV8lus48niMvvX+5cpM9+izv9imiVDYdSQgpmJeuq21k5I+raNleaekh621O4wxXzXG3DD1tG9KypX0iDFmmzHm0anX9kj6Bx0t7pslfXXqMQAAcJz/fGm/mruG9LlrFsrrMa7jSJJmFmTpa3+yVNv39ev7a5tcxwFiKnZ7g45jrX1M0mMnPPZ3x73/1lO89l5J90YvHQAAiW1sMqTvrW3SsrICXbdkhus4f+TtS0v1vqoy/eDp3bpiYYmqAn7XkYCYiOsbKQEAwNn7+Ytt2t83oi9cu1DGxMdV7uN95V2LNds3TX/5i206MjrhOg4QE5RuAACSyPD4pH7w9B5dUuHX6spi13FOKi8rXd/5wArt7x3R3//2NddxgJigdAMAkEQe2NiqrsExffG6+LzKfUx10K//9Zb5+uWWfXrslYOu4wBRR+kGACBJ9I9M6N+f3aO3LCxRdTD+Z6U/c3WllpcV6P/8+hUd6h91HQeIKko3AABJ4p71zeofmdDnr13oOsoZSfd69J0PrtDYRFhf/OV2hcOcVonkRekGACAJdA+O6Z7n9uodS0t14ewC13HO2NySXP3tOxdrfVOX7t/Y4joOEDWUbgAAksC/PbNHIxMhffaaBa6jnLWba8v11gum6+t/aFDDoQHXcYCooHQDAJDgDvaP6MEXWvWei8s0f3qu6zhnzRijr793mfKz0vSXD23T6ETIdSQg4ijdAAAkuO+v2y1rrf7i6krXUc5ZcW6mvvG+ZWo4dESP1Le7jgNEHKUbAIAE1to9pIc3t+vm2jkq92e7jnNerlo0Q2W+adq4p9t1FCDiKN0AACSw7z7VpDSv0affMt91lIioDfq1uaVH1rLJBMmF0g0AQILadeiIfrNtv25bGdT0/CzXcSKipsKvrsFx7e0ach0FiChKNwAACerbT+5SbkaaPrlmnusoEVMzdajP5pYex0mAyKJ0AwCQgLa39+nxHYf18dVz5cvJcB0nYuaV5KgoJ0Ob9va6jgJEFKUbAIAE9K0ndsmXna6PXh50HSWijDGqDvq40o2kQ+kGACDBvNDcrfVNXfrUlfOVl5XuOk7E1VYUqa1nWIcHRl1HASKG0g0AQAKx1upbj+/SjPxMfeSygOs4UVE7Nde9aS9Xu5E8KN0AACSQZ3Z1qr61V39+VaWy0r2u40TFBaV5ysnwUrqRVCjdAAAkiHDY6ltP7FK5f5o+UF3uOk7UpHk9ujjAXDeSC6UbAIAE8ftXD2nHgQF99q0LlJGW3L/Ca4N+7Tp8RP3DE66jABGR3P/FAgCQJCZDYX37yV2qnJ6rG1fMdh0n6moq/LJWqm/lajeSA6UbAIAE8Out+7Wnc0ifv3aBvB7jOk7UrSgvVLrXaBMjJkgSlG4AAOJcKGz1vbVNWjq7QNctmek6TkxkpXu1vKxQm7mZEkmC0g0AQJzb1t6nfb0j+vjqChmT/Fe5j6mp8OuV/f0aGQ+5jgKcN0o3AABxbn1Tp4yR1lSWuI4SU7VBvyZCVlvbORIeiY/SDQBAnKtr7NSyskL5cjJcR4mpiwM+GSNt3kvpRuKjdAMAEMf6hye0rb1PayqLXUeJuYJp6Vo0M5993UgKlG4AAOLYxj1dCltpzYLUGi05pjbo00ttvZoMhV1HAc4LpRsAgDhW19SpvMw0rSgvdB3FiZoKv4bHQ9pxYMB1FOC8ULoBAIhT1lrVNXbpsnlFSvem5q/s2qBfkhgxQcJLzf+CAQBIAM1dQ9rfN5KyoyWSND0/S8GibG1iXzcSHKUbAIA4VdfYKUm6IoVLtyTVBP3a3NKjcNi6jgKcM0o3AABxqq6xU8GibJX7s11Hcaqmwq/e4Qnt6Rx0HQU4Z5RuAADi0NhkSC8096T0aMkxx+a6NzHXjQRG6QYAIA5taenVyEQo5U6hPJlAUbZK8jK1mbluJDBKNwAAcejZpk6leYwunVfkOopzxhjVBv3a3MLJlEhclG4AAOLQ+sYuVQV8ys1Mcx0lLtQEfdrfN6L9fSOuowDnhNINAECc6TwyptcODjDPfZzaiqNX/BkxQaKidAMAEGee2310VSDz3P9j4cw85WWlcTMlEhalGwCAOFPX2KWinAwtmZXvOkrc8HqMqgM+DslBwqJ0AwAQR8Jhq/VNnbq8slgej3EdJ67UVPi1u2NQPUPjrqMAZ43SDQBAHNl5aEBdg+NazWjJGxzb172ZERMkIEo3AABxpK6xS5K0prLYcZL4s7SsQBlpHm6mREKidAMAEEfqGju1aGaepudnuY4SdzLTvFpRXsiVbiQkSjcAAHFieHxS9a0c/X4ql1T49eqBAQ2NTbqOApwVSjcAAHHiheZuTYQsqwJPoSboVyhstbWtz3UU4KxQugEAiBN1jV3KSveoOuhzHSVuXRzwyWOkTXu7XUcBzgqlGwCAOFHX1KlLKoqUle51HSVu5WamacmsAg7JQcKhdAMAEAf29Q6ruXOIee4zUBP0a2tbn8Ynw66jAGeM0g0AQBw4tirwigWsCjyd2gqfxibDemV/v+sowBmjdAMAEAfWN3WqtCBL80pyXUeJe9UckoMEROkGAMCxyVBYz+3u0prKEhnD0e+nU5ybqXklORySg4RC6QYAwLHt+/p0ZHSSee6zUFvh1+aWHoXD1nUU4IxQugEAcKyusUseI62aX+Q6SsKoCfo1MDqpXYePuI4CnBFKNwAAjtU1dWpZWaEKszNcR0kYNcx1I8FQugEAcKh/eELb2/sYLTlLZb5pKi3I0ibmupEgKN0AADi0YU+XwlZaU8mqwLNhjFFN8Ohct7XMdSP+UboBAHCorrFTeVlpWlFe6DpKwqmp8OvwwJjae0ZcRwFOi9INAIAj1lrVNXZq1bxipXn5lXy2aqfmujkSHomA/8IBAHBkT+eQDvSPajWnUJ6Tyum5KsxOZ183EgKlGwAAR+oaOyVJayq5ifJceDxG1QE/V7qRECjdAAA4UtfUqbnFOSr3Z7uOkrBqK3za2zWkjiOjrqMAp0TpBgDAgbHJkF5o7tZqtpacl2P7uutbeh0nAU6N0g0AgAP1Lb0anQizn/s8XTi7QNPSvezrRtyjdAMA4EBdY6fSvUaXzuXo9/OR7vXoojmFnEyJuEfpBgDAgbqmLlUFfMrJTHMdJeHVBP3aeXBAR0YnXEcB3hSlGwCAGOs4MqqdBwcYLYmQSyr8CltpSytz3YhflG4AAGJsfWOXJFYFRspFc3xK8xhGTBDXKN0AAMTY+qZOFeVkaHFpvusoSWFahlcXzi7gZkrENUo3AAAxFA5brW/q0urKYnk8xnWcpFFb4df29n6NToRcRwFOitINAEAMvXZwQN1D48xzR1hN0K/xUFgv7+t3HQU4KUo3AAAxVNd09Oj3yzkUJ6KqAz5JYq4bcYvSDQBADNU1duqC0nxNz8tyHSWp+HIytGBGLnPdiFuUbgAAYmRobFJbWnu1ZgFXuaOhJujXS629CoWt6yjAG1C6AQCIkReauzURsqwKjJLaCr+OjE1q58EB11GAN6B0AwAQI3WNnZqW7lV10Oc6SlKqrfBLEiMmiEtOSrcx5npjzC5jzG5jzJdO8vk1xpiXjDGTxpj3nfC5kDFm29Tbo7FLDQDA+alr6tKlc/3KTPO6jpKUSgumqcw3jZspEZdiXrqNMV5JP5T0NkmLJd1sjFl8wtPaJN0u6Wcn+RIj1toVU283RDUsAAAR0t4zrL1dQ1rNaElU1Qb92tzSI2uZ60Z8cXGlu1bSbmtts7V2XNJDkm48/gnW2hZr7cuSwg7yAQAQccdWBbKfO7pqKvzqGhzX3q4h11GAP+KidM+W1H7cx/umHjtTWcaYemPMC8aYd0c2GgAA0VHX2KnZhdM0ryTHdZSkVhM8OtfNiAniTSLeSBmw1lZL+pCk7xpj5p34BGPMnVPFvL6zszP2CQEAOM5kKKyNu7u1urJYxnD0ezTNK8lRUU6GNu3tdR0F+CMuSvd+SeXHfVw29dgZsdbun/rfZknPSLroJM+521pbba2tLinhz3gAALe2tffpyNgkoyUxYIxRzdRcNxBPXJTuzZIqjTEVxpgMSTdJOqMtJMYYnzEmc+r9YkmrJL0WtaQAAERAXVOXPEZaNY9DcWKhpsKvtp5hHeofdR0FeF3MS7e1dlLSpyU9LmmnpIettTuMMV81xtwgScaYGmPMPknvl3SXMWbH1MsvkFRvjNku6WlJX7fWUroBAHGtrrFTy8sLVZCd7jpKSqidmuvexNVuxJE0F9/UWvuYpMdOeOzvjnt/s46OnZz4uo2SlkY9IAAAEdI3PK6X9/Xpz6+qdB0lZVxQmqecDK827+3RDctnuY4DSErMGykBAEgYG3Z3K2xZFRhLaV6PLg74mOtGXKF0AwAQRXWNncrLStPysgLXUVJKbdCvXYePqH94wnUUQBKlGwCAqLHWqq6pU5fPL1aal1+5sVRT4Ze1Un0rV7sRH/gJAABAlOzpHNTB/lFGSxxYUV6oDK+HmykRNyjdAABEybONXZKk1ZWsCoy1rHSvlpUVaPNeSjfiA6UbAIAoqWvs1NySHJX5sl1HSUk1FX69vK9fI+Mh11EASjcAANEwOhHSi3u7taaS0RJXaoN+TYattrZzJDzco3QDABAF9S29Gp0Ia80CRktcuTjgkzHS5r2UbrhH6QYAIArqmjqV4fXo0rlFrqOkrIJp6Vo0M5993YgLlG4AAKKgrrFT1UGfsjOcHP6MKbVBn15q69VkKOw6ClIcpRsAgAjrGBhVw6EjWs08t3M1FX4Nj4e048CA6yhIcZRuAAAirK7p6KpA5rndqw36JYkREzhH6QYAIMLWN3WqODdTF8zMdx0l5U3Pz1KwKFsvsq8bjlG6AQCIoHDYan1Tl1ZXFsvjMa7jQFJN0K/6lh6Fw9Z1FKQwSjcAABG048CAeobGGS2JIzUVfvUOT2hP56DrKEhhlG4AACKorqlTknT5fG6ijBfH5ro3MdcNhyjdAABEUF1jpxaX5qskL9N1FEwJFGWrJC9Tm5nrhkOUbgAAImRwbFJbWnu1ZgFXueOJMUa1Qb82t3AyJdyhdAMAECEv7OnWZNgyzx2Haiv82t83on29w66jIEVRugEAiJC6pk5NS/eqKuBzHQUnqGFfNxyjdAMAECF1jZ26bF6RMtO8rqPgBAtn5ikvK02b9jJiAjco3QAAREBb97Bauoe1ppLRknjk9RhVB3xc6YYzlG4AACLg2KrA1dxEGbdqKvza3TGonqFx11GQgijdAABEQF1jp2YXTtPc4hzXUfAmapnrhkOUbgAAztNEKKzn93RrzYISGcPR7/FqaVmBMtI87OuGE5RuAADO07b2Ph0Zm2SeO85lpnl1UXkhV7rhBKUbAIDzVNfYKa/HaOV8Sne8q63w69UDAxoam3QdBSmG0g0AwHmqa+rSivJCFUxLdx0Fp1ET9CsUtnqpjdWBiC1KNwAA56F3aFwv7+vTakZLEsLFAZ88Rsx1I+Yo3QAAnIfndnfJWmkNqwITQm5mmpbMKtAm5roRY5RuAADOw/qmTuVnpWnZ7ALXUXCGaoJ+bW3r0/hk2HUUpBBKNwAA58haq7rGLl1eWaw0L79SE0VthU9jk2G9sr/fdRSkEH5CAABwjpo6BnVoYFRrKhktSSTVHJIDB05Zug0b/gEAeFN1jRz9noiKczM1rySHmykRU6e70v0HY0xhTJIAAJBg6pq6NK8kR7MLp7mOgrNUW+HX5pYehcPWdRSkiNOV7nmSXjLGLItFGAAAEsXoREgvNneztSRB1QT9Ghid1K7DR1xHQYo4XemukvSapOeNMR+OQR4AABLC5pYejU2GmedOUDXMdSPGTlm6rbX91tp3SvqWpAeNMd8zxnhjEw0AgPhV19ipDK9Hl8z1u46Cc1Dmm6bSgixtYq4bMXJG20ustV+RdIOkj0haZ4ypjGoqAADi3PqmLtVU+JSdkeY6Cs6BMUY1waNz3dYy143oO5uVgWsl3SVptaQGY0y/MWajMeYuY8ynjTFXRiUhAABx5vDAqBoOHdFqRksSWk2FX4cHxtTeM+I6ClLAaf95PnVV+1OSbpOUI+khSY9LWijpQknXSvpTSVYSoycAgKR3bFUg89yJ7ZKKo6NBm1p6NKco23EaJLtTlm5jzJOSrpLUI+nfJP3QWnvgJM/L09ECDgBA0lvf1KXi3ExdUJrnOgrOw/ySXBVmp2vT3m69r6rMdRwkudNd6Z4h6ROSfmqtHX2zJ1lrj0h6PpLBAACIR+Gw1XO7u3TlghJxhlxi83iMqgN+bW7pdR0FKeCUpdtay35uAACOs31fn3qGxnXFQkZLkkFthU9P7TysjiOjmp6X5ToOktjZ3EgJAEDKW9fQIa/H6AoOxUkKx/Z113O1G1FG6QYA4Cw8tbNDVQGfCrMzXEdBBFw4u0BZ6R4OyUHUUboBADhDB/pGtPPggK5eNN11FERIutejFeWF2tLKlW5EF6UbAIAztK6hQ5J09QWU7mRSE/Rrx4EBDY1Nuo6CJEbpBgDgDK1r6FCgKFvzSnJdR0EEVQV8CoWttrf3uY6CJEbpBgDgDIyMh7Rhd5euWjSdVYFJ5uKAT8ZI9YyYIIpOWbqNMU8YYxae8NhVxpic6MYCACC+bNzTpbHJsK5eNMN1FERYfla6Fs7I42ZKRNXprnS/VVLBsQ+MMV5JT+roEfAAAKSMtQ0dysnwqnbq6HAkl+qgT1vb+hQKW9dRkKTOZbyEv6kBAFKKtVbrdnZozYISZaQxmZmMaoJ+DY5NquHQgOsoSFL85AAA4DReOzigQwOjuopVgUmrKuCTJFYHImrOpHSf7O8s/O0FAJAy1u3skDHSlQsp3clqduE0lRZkaTMnUyJK0s7gOY8bY05cXLn2JI/JWstPIwBA0lnb0KHlZYUqyct0HQVRYoxRVcCnLdxMiSg5Xen++5ikAAAgTnUeGdP2fX363FsXuI6CKKsO+PS7lw9qf9+IZhdOcx0HSeaUpdtaS+kGAKS0p3d1yFrpKk6hTHrVwaObaepbejR7xWzHaZBsTlm6jTF/dxZfy1pr/+E88wAAEFfW7exQaUGWFpfmu46CKFs0M085GV5tae3VjZRuRNjpxkv+r6QRSUM6/apAK4nSDQBIGmOTIa1v6tS7L5rNKZQpIM3r0cUBHzdTIipOt71kj6R0SVskfUHSXGttyZu88Xc3AEBS2bS3R0PjIV3NaEnKqAr4tOvQgAZGJ1xHQZI5Zem21lZKWilph45exT5sjPlPY8z7jTHcYQAASGprd3YoK92jlfOKXUdBjNQE/QpbaWtbn+soSDKn3dNtra231n7BWjtH0vWSDkn6gaQOY8x/GGPWRDskAACxZq3V2obDWjWvWFnpXtdxECMrygvl9RhWByLizupESmttnbX2U5LKJf27pA9K+stoBAMAwKU9nYNq7xlha0mKyclM0+LSfOa6EXFnVbqNMauMMd+X1CrpzyT9UtL3ohEMAACX1u7skCSOfk9BVQGftrX3aSIUdh0FSeS0pdsYc7Ex5hvGmFZJa3X0KvdnJU231t5krX022iEBAIi1tQ0dWlyar9ICbmFKNdVBn0YmQtp5cMB1FCSRU5ZuY8wuSS9IWibpKzpatN9trX3IWjsci4AAAMRa3/C4trT2srUkRVUHjh6Sw4gJIul0V7orJU1KqpL0DUm7jTEdb/YW9bQAAMTAs42dCoUtoyUpamZBlsp807SllZspETmnOxyHY+ABACln7c4OFedmaHlZoesocKQm6Ndzu7tkreVgJETEKUu3tZbSDQBIKZOhsJ7Z1aHrlsyUx0PZSlVVAZ9+vXW/2ntGNKco23UcJIHTXenW1CE4b5cUlHRQ0lpr7eEo5wIAwIktrb0aGJ1knjvF1QSPzXX3ULoREacs3caYuZKe0tHCfcyAMeYD1tonohkMAAAX1jV0KN1rdHlliesocKhyeq7ys9JU39qr91aVuY6DJHC6Gym/ISksabWkbElLJG2VdFeUcwEA4MTahg5dOrdIuZmn/WMwkpjHY3RxwKd6TqZEhJyudF8m6W+stRustaPW2p2SPiFpjjGmNPrxAACIndbuIe3uGGRrCSQdHTFp6hhU3/C46yhIAqcr3aWSmk94bI8kI2lmVBIBAODIugZOocT/qAr4JEkvtbGvG+fvTI6Bt1FPAQBAHFjX0KH503MVKMpxHQVxYHlZodK9hkNyEBFnMrD2uDFm8iSPrz3xcQoxJjcAACAASURBVGstlwYAAAlpcGxSLzR366OrKlxHQZyYluHVklkF2kLpRgQ4ORzHGHO9pO9J8kr6sbX26yd8fo2k7+ro8fM3WWt/edznbpP0N1Mf/n/W2geikREAkFrWN3ZqIsQplPhjNUGfHni+VWOTIWWmeV3HQQKL+eE4xhivpB9KukbSPkmbjTGPWmtfO+5pbZJul/SFE17rl/QVSdU6OvayZeq1/BMUAHBe1jZ0KD8r7fU5XkCSqgJ+/Wj9Xr26f4D/b+C8nMlMd6TVStptrW221o5LekjSjcc/wVrbYq19WUfXFR7vOklPWmt7por2k5Kuj0VoAEDyCoetnm7o0JULpyvN6+JXI+LVsaLN6kCcLxc/WWZLaj/u431Tj0X7tQAAnNT2fX3qHhrnFEq8QUlepiqKc1Tfyh/VcX6S8p/zxpg7jTH1xpj6zs5O13EAAHFuXUOHvB6jKxZwCiXeqCrg05bWXlnLQjecOxele7+k8uM+Lpt6LGKvtdbeba2tttZWl5TwAxQAcGprd3aoKuBTYXaG6yiIQzVBn3qGxtXcNeQ6ChKYi9K9WVKlMabCGJMh6SZJj57hax+XdK0xxmeM8Um6duoxAADOycH+Eb12cEBXs7UEb6Iq4JckVgfivMS8dFtrJyV9WkfL8k5JD1trdxhjvmqMuUGSjDE1xph9kt4v6S5jzI6p1/ZI+gcdLe6bJX116jEAAM7JsVMomefGm5lXkiNfdro2czMlzsOZHI4TcdbaxyQ9dsJjf3fc+5t1dHTkZK+9V9K9UQ0IAEgZ63Z2aI4/W/NKcl1HQZwyxqgq4NcWbqbEeUjKGykBADgTI+MhPbe7S1ctmi5jjOs4iGPVQZ+au4bUNTjmOgoSFKUbAJCyNu7p0thkmNESnFZN8Oi+bq5241xRugEAKWttQ4dyMryqrfC7joI4d+HsAmWkeTgkB+eM0g0ASEnWWq3b2aHVlSXKTPO6joM4l5nm1fKyAg7JwTmjdAMAUtJrBwd0aGBUVzFagjNUFfDr1f39Gp0IuY6CBETpBgCkpHU7O2SM9JaFlG6cmZqgTxMhq+3tfa6jIAFRugEAKWltQ4eWlxWqJC/TdRQkiKrA0ZspGTHBuaB0AwBSTueRMW3f18cplDgrhdkZmj89l5spcU4o3QCAlPPMrg5ZK+a5cdZqgj5tae1VOGxdR0GCoXQDAFLOuoYOzczP0uLSfNdRkGCqAn4NjE6qqWPQdRQkGEo3ACCljE+GVdfYqasu4BRKnL1jh+TUtzJigrND6QYApJRNe3s0NB5inhvnZI4/W8W5mapv4WZKnB1KNwAgpTy187Ay0zxaOa/YdRQkIGOMaoI+rnTjrFG6AQApw1qrtQ2HtWp+saZlcAolzk1VwKf2nhEdHhh1HQUJhNINAEgZezoH1d4zoqsYLcF5qAn6JYkRE5wVSjcAIGWs3dkhSbqaVYE4D4tn5WtaupcRE5wVSjcAIGWsbejQ4tJ8lRZMcx0FCSzd69Hy8gKudOOsULoBACmhb3hcW1p7ucqNiKgJ+vXawQENjU26joIEQekGAKSEZxs7FQpb5rkREVUBn0Jhq23tfa6jIEFQugEAKWFdQ4eKcjK0vKzQdRQkgYsDPhnDzZQ4c5RuAEDSmwyF9cyuTr1l0XR5PJxCifOXn5WuhTPyuJkSZ4zSDQBIei+19al/ZIJTKBFRNUG/trb1KRS2rqMgAVC6AQBJb+3Ow0r3Gl1eySmUiJzqoE+DY5NqODTgOgoSAKUbAJD01jZ06JKKIuVlpbuOgiRSFfBJYq4bZ4bSDQBIaq3dQ9rdMcjWEkTc7MJpKi3IUn0rpRunR+kGACS1dQ2cQonoMMaoKuBTfQs3U+L0KN0AgKS2rqFD86fnKlCU4zoKklBN0K+D/aPa3zfiOgriHKUbAJC0Bscm9UJzN1tLEDX/M9fN1W6cGqUbAJC0nmvq1ESIUygRPYtm5ik3M42bKXFalG4AQNJau7ND+Vlpr1+NBCItzevRRXMKuZkSp0XpBgAkpXDY6uldHbpy4XSlefl1h+ipCvjUcGhAA6MTrqMgjvFTCACQlF7e36+uwXG2liDqaoJ+WSttbetzHQVxjNINAEhK63YelsdIVywocR0FSW5FeaG8HsPNlDglSjcAICk9tbND1QG/CrMzXEdBksvJTNPi0nxupsQpUboBAEnnYP+IXjs4oKsYLUGMVAV82treq4lQ2HUUxClKNwAg6bx+CiWrAhEjNUG/RifCeu3AgOsoiFOUbgBA0lm3s0Nz/NmaPz3XdRSkiOrg1CE5rA7Em6B0AwCSysh4SM/t7tJVi6bLGOM6DlLEjPwslfmmcTMl3hSlGwCQVJ5v7tLYZJhVgYi5mqBf9a29sta6joI4ROkGACSVtTs7lJPhVW2F33UUpJiqgE+dR8bU1jPsOgriEKUbAJA0rLVa19Ch1ZUlykzzuo6DFFMTPPoPPVYH4mQo3QCApLHz4BEd7B9lVSCcqJyeq/ysNNW3MteNN6J0AwCSxrqGw5KktyykdCP2PB6jqoCPK904KUo3ACBpPLWzQ8vLC1WSl+k6ClJUddCvpo5B9Q2Pu46COEPpBgAkhc4jY9q+r48DceBUVeDovu4t7OvGCSjdAICk8MyuDlkrXUXphkPLywqV7jUckoM3oHQDAJLCuoYOzczP0pJZ+a6jIIVNy/BqyawCDsnBG1C6AQAJb3wyrLrGTl11AadQwr2aoE/b9/VrbDLkOgriCKUbAJDwNu3t0dB4iHluxIWqgF/jk2G9ur/fdRTEEUo3ACDhrW04rMw0j1bOK3YdBVB18OjNlKwOxPEo3QCAhGat1dqdHVo1v1jTMjiFEu4V52aqojhHmyndOA6lGwCQ0PZ0DqmtZ5itJYgr1QGfXmrrlbXWdRTECUo3ACChHTuFktKNeFId9KlnaFzNXUOuoyBOULoBAAlt7c4OXVCar1mF01xHAV5XFfBLEqsD8TpKNwAgYR3oG1F9ay9bSxB35pXkyJedzs2UeB2lGwCQsO7bsFeSdFNtueMkwB8zxqgq4OdkSryO0g0ASEj9IxP62YttesfSUpX5sl3HAd6gJujT3q4hdQ2OuY6COEDpBgAkpJ+92Kah8ZDuXDPXdRTgpNjXjeNRugEACWdsMqT7NuzV5fOLdeHsAtdxgJO6cHaBMtI82tLKzZSgdAMAEtB/bTugjiNjXOVGXMtM82p5WQGH5EASpRsAkGDCYau765p1QWm+Vldy7DviW1XArx0H+jU6EXIdBY5RugEACeXpXR3a3TGoO9dUyBjjOg5wSjVBnyZCVtvb+1xHgWOUbgBAQrmrrlmzCrL0zmWzXEcBTqsqMHUzJasDUx6lGwCQMLa29WrT3h599PIKpXv5FYb4V5idocrpuZxMCUo3ACBx3F3XrLysNN1UO8d1FOCMVQd92tLaq3DYuo4ChyjdAICE0NI1pD/sOKRbLg0oNzPNdRzgjFUH/BoYnVRTx6DrKHCI0g0ASAg/fq5Z6R6P7lgZdB0FOCvHDsnZzIhJSqN0AwDiXvfgmB6p36c/uWi2pudnuY4DnJU5/mwV52ZqCzdTpjRKNwAg7j3wfKvGJsP60zUVrqMAZ80Yo5qgT/WcTJnSKN0AgLg2Mh7ST55v0VsvmK750/NcxwHOSVXAp/aeER0eGHUdBY5QugEAce2RLe3qHZ7QJ66Y5zoKcM5qgn5JUj1HwqcsSjcAIG5NhsL68fq9umhOoaqnDhkBEtHiWfmalu7lZsoURukGAMStP+w4pLaeYX1izVyOfEdCS/d6tKK8kJspUxilGwAQl6y1uruuWRXFObpm8UzXcYDzVh306bWDAxoam3QdBQ5QugEAcemF5h69vK9fH19dIa+Hq9xIfFUBn0Jhq23tfa6jwAFKNwAgLt1Vt0dFORl678VlrqMAEXFxwCdjOCQnVVG6AQBxZ9ehI3pmV6duWxlUVrrXdRwgIvKz0rVwRh5z3SnKSek2xlxvjNlljNltjPnSST6faYz5xdTnXzTGBKceDxpjRowx26be/j3W2QEA0Xd3XbOmpXv1kUsDrqMAEVUT9Oul1l5NipGpVBPz0m2M8Ur6oaS3SVos6WZjzOITnvYxSb3W2vmSviPpn4/73B5r7Yqpt0/GJDQAIGYO9o/o0e379cGacvlyMlzHASKqOujT0HhIDdklrqMgxlxc6a6VtNta22ytHZf0kKQbT3jOjZIemHr/l5KuNuyKAoCUcN+GFoXCVh+7nCPfkXyqpw7J2ZI323ESxJqL0j1bUvtxH++beuykz7HWTkrql1Q09bkKY8xWY8yzxpjV0Q4LAIidgdEJ/ezFNr19aanK/dmu4wARN7twmkoLsrSZ0p1y0lwHOEsHJc2x1nYbY6ok/cYYs8RaO3D8k4wxd0q6U5LmzJnjICYA4Fz8/MU2DY5N6hNrOPIdyasq4FN9R5msxGR3CnFxpXu/pPLjPi6beuykzzHGpEkqkNRtrR2z1nZLkrV2i6Q9khac+A2stXdba6uttdUlJcxMAUAiGJ8M694Ne7VyXpGWlhW4jgNETU3Qr0OZedqfke86CmLIReneLKnSGFNhjMmQdJOkR094zqOSbpt6/32S1llrrTGmZOpGTBlj5kqqlNQco9wAgCj6r237dXhgTHeumes6ChBVVQGfJOmF/PLTPBPJJOale2pG+9OSHpe0U9LD1todxpivGmNumHraPZKKjDG7JX1O0rG1gmskvWyM2aajN1h+0lrLhnkASHDWWv1ofbMWzczTFQv4CyWS2+LSfAVHevXzGctdR0EMOZnpttY+JumxEx77u+PeH5X0/pO87leSfhX1gACAmHpmV6caDw/q2x9YLpZVIdl5PEa3Hn5JXw1erZf39WlZWaHrSIgBTqQEADh3V90elRZk6V3LZ7mOAsTE+zpfVU5oXPdvbHEdBTFC6QYAOLW9vU8vNPfoo6sqlO7l1xJSQ35oXO/veEW/235QnUfGXMdBDPDTDQDg1N11zcrLTNNNtdxUhtRy6+GtGg+F9fNNba6jIAYo3QAAZ9q6h/X7Vw/qw5cGlJeV7joOEFNzR3t15cIS/fSFVo1Phl3HQZRRugEAzvz4uWZ5PUZ3rAq6jgI4cfvKoDqOjOn3rx50HQVRRukGADjRMzSuh+vb9e4VszUjP8t1HMCJNZUlmluco/s2tLiOgiijdAMAnHjw+RaNToQ5DAcpzeMxum1lUNva+7S1rdd1HEQRpRsAEHMj4yE9+Hyrrl40XZUz8lzHAZx6b1WZcjPT9ADrA5MapRsAEHO/fGmfeobGucoNSMrNTNP7q8v0368cVMfAqOs4iBJKNwAgpkJhqx+vb9by8kLVVvhdxwHiwm2XBTUZtvqPF1kfmKwo3QCAmHp8xyG1dg/rE2vmcuQ7MCVYnKOrFk7Xf7zYprHJkOs4iAJKNwAgZqy1uquuWYGibF23ZKbrOEBcuX1VUF2DY/rvl1kfmIwo3QCAmNm0t0fb2/v08dVz5fVwlRs43uXzizV/eq7u29Aia63rOIgwSjcAIGbuqmuWPydD768qcx0FiDvGHF0f+Mr+fr3U1uc6DiKM0g0AiImmw0e0rqFDt14WUFa613UcIC6956LZystK0/2sD0w6lG4AQEzcXdesrHSPbr0s6DoKELdyMtP0wepy/f6VgzrUz/rAZELpBgBE3eGBUf1m2359oLpc/pwM13GAuHbrZUGFrNV/vNjqOgoiiNINAIi6ezfsVShs9fHLOQwHOJ05Rdm6etEM/ezFNo1OsD4wWVC6AQBRdWR0Qj97oU1vu7BUc4qyXccBEsJHVwXVPTSu37E+MGlQugEAUfXQpnYdGZvkyHfgLFw2r0gLZuTqvg17WR+YJCjdAICoGZ8M694Ne3XpXL+Wlxe6jgMkDGOMbl9ZoR0HBlTf2us6DiKA0g0AiJrfbj+gg/2j+sSaea6jAAnn3RfNUsG0dN2/ocV1FEQApRsAEBXWWv1ofbMWzMjVlQtLXMcBEk52RppuqinXH3Yc0oG+EddxcJ4o3QCAqHi2sVMNh47ozjXzZAxHvgPn4pZLA7LW6qcvsD4w0VG6AQBRcXdds2bmZ+mG5bNcRwESVrk/W9csnqGfb2J9YKKjdAMAIu6Vff3auKdbd6wKKiONXzXA+bh9ZYV6hyf06LYDrqPgPPCTEAAQcXfV7VFuZppuvmSO6yhAwrt0rl+LZubpXtYHJjRKNwAgotozC/TYKwf14UvmKD8r3XUcIOEZY3THqqAaDh3Ri3t7XMfBOaJ0AwAi6p6ZVfJ6jO5YVeE6CpA0blwxW4XZrA9MZJRuAEDE9KZl6RfTl+qG5bM1syDLdRwgaWSle3Vz7Rw98doh7esddh0H54DSDQCImJ/MuEgj3gyOfAei4JZLAzLG6CesD0xIlG4AQETs7xvR/TMv1lt692jhzDzXcYCkM7twmq5bMkMPbWrXyDjrAxMNpRsAcN66B8f0kXte1ITx6kttda7jAEnr9pUV6h+Z0G+27XcdBWeJ0g0AOC+DY5O64/7N2t87ont2/UoLR7pcRwKSVk3Qp8Wl+bqP9YEJh9INADhnY5Mh3flgvXYcGNC/fvhi1R7h6hsQTcYY3b4qqMbDg3p+T7frODgLlG4AwDkJha3+4ufbtHFPt775vmW6+oIZriMBKeGG5bPkz8nQfRtbXEfBWaB0AwDOmrVWX/71K/rDjkP623cu1nsuLnMdCUgZWelefah2jp7aeVjtPawPTBSUbgDAWfvG47v00OZ2ffot8/WxyzkEB4i1Wy4NyGOMHny+xXUUnCFKNwDgrNxdt0f/9swefeiSOfr8tQtcxwFS0syCLL3twpl6aHO7hsYmXcfBGaB0AwDO2MP17fraYw16x9JS/cONF8oY4zoSkLLuWBXUkdFJ/XorNzAnAko3AOCMPL7jkL70q5e1urJY3/7gcnk9FG7ApYvn+LR0doHu39jC+sAEQOkGAJzW83u69ec/36qlZYX691uqlJnmdR0JSHnGGN2+MqjdHYN6bjf78eMdpRsAcEqv7u/Xnz5Yrzn+bN1/e41yMtNcRwIw5Z3LS1Wcm6H7N7S4joLToHQDAN5Uc+egbrt3kwqmpesnH6uVLyfDdSQAx8lMO7o+cN2uDrV0DbmOg1OgdAMATupg/4g+cs8mSdJPPlar0oJpjhMBOJlbLg3Ia4wefL7VdRScAqUbAPAGvUPjuvWeTeofmdD9d9Rqbkmu60gA3sT0/Cy9Y1mpHqlv1yDrA+MWpRsA8EeGxiZ1x/2b1dozrB/dWq2lZQWuIwE4jdtXBnVkbFL/+dI+11HwJijdAIDXjU2G9MmfbtHL+/r0/Zsv0mXzilxHAnAGLprj0/LyQt2/sUXhMOsD4xGlGwAgSQqFrT738Hatb+rS19+7TNctmek6EoCzcMfKoJo7h1TX1Ok6Ck6C0g0AkLVWf/tfr+q/Xz6o//P2RfpAdbnrSADO0tuXlqokL1P3b2xxHQUnQekGAOjbTzbqZy+26ZNXzNOda+a5jgPgHGSkefThS+bomV2dau4cdB0HJ6B0A0CKu/e5vfr+ut36YHW5/vf1C13HAXAePnTJHKV7WR8YjyjdAJDCfr11n776u9d03ZIZ+sc/uVDGGNeRAJyH6XlZeteyWXqkvl1HRidcx8FxKN0AkKLWNRzWFx55WZfNLdL3brpIaV5+JQDJ4LaVQQ2Nh/TLLawPjCf8hAWAFLRpb4/+7KcvaXFpvu6+tUpZ6V7XkQBEyPLyQl08p1APsD4wrlC6ASDFvHZgQB97YLNmF07T/XfUKC8r3XUkABF2+6oKtXQP69lG1gfGC0o3AKSQ1u4h3XrvJuVmpuknH79ERbmZriMBiIK3XThTM/Izde+Gva6jYAqlGwBSRMfAqG6550WFwmH95GO1ml04zXUkAFGS7vXolksCWt/Upd0dR1zHgSjdAJAS+ocndOu9m9Q9OK777qjV/Ol5riMBiLKbL5mjDK9HD2xkfWA8oHQDQJIbGQ/pow9s1p7OQd39kWqtKC90HQlADBTnZupdy2fpVy/tU/8I6wNdo3QDQBKbCIX1Z/+xRS+19ep7N12kyyuLXUcCEEN3rApqeDykR+rbXUdJeZRuAEhS4bDVFx7Zrmd2deof371Ub19a6joSgBi7cHaBaoI+Pfh8q0KsD3SK0g0ASah/eEJf/s0r+q9tB/TF6xbqQ5fMcR0JgCO3r6xQW8+wnm7ocB0lpaW5DgAAiJxdh47o/o0t+vXWfRqdCOsTa+bqU1fOcx0LgEPXLpmh0oIs3bdxr966eIbrOCmL0g0ACS4Utlq787Du39iijXu6lZnm0btXzNZtK4NaPCvfdTwAjqV7Pbrl0oC++fgubdzdpZXzubfDBUo3ACSo/uEJPVzfrgeeb9G+3hHNKsjSX12/UDfVzJE/J8N1PABx5ObaOfrpC6360I9f1DuWleqvrluoQFGO61gphdINAAmm6fDREZL/fGm/RiZCqq3w68tvv0DXLJ6hNC+36gB4I39Ohp783BW6u65ZP6pr1hM7DunDlwT0masr+Ud6jFC6ASABhMJW6xo6dP/Gvdqwu1sZaR69e8Us3bYyqCWzClzHA5AAcjPT9LlrFuiWS+boO0816cHnW/SrLfv0ySvn6aOrKjQtw+s6YlKjdANAHOsfmdAjUyMk7T0jmpmfpS9et1A31zJCAuDcTM/P0j+9Z6k+dnlQX//9Ln3z8V36yfOt+ty1C/Tei8vk9RjXEZMSpRsA4tDujqMjJL/acnSEpCbo05euv0DXLpmhdEZIAETA/Ol5+vFt1XqxuVtf+32D/uqXL+ue9Xv1pbcv0pULSmQM5TuSKN0AECdCYaunGzr0wPMtWt/UpYw0j25YPku3rwzqwtmMkACIjkvmFuk3n1qpx145pG883qA77tuslfOK9Ndvu0BLy/jZEymUbgBw7NgIyYPPt6qtZ/j1EZKbaspVlJvpOh6AFGCM0TuWleqaxTP0sxdb9f/W7da7fvCcblwxS1+4dqHK/dmuIyY8SjcAOLK7Y1APbGzRr17ap+HxkKoDPv3V9Qt13ZKZjJAAcCIjzaPbV1XoPVVluuvZPfrx+r36/SuHdOtlAX36qvkqzOZeknNF6QaAGAqHrZ5p7NB9G6ZGSLwevWtqhIQ/4wKIF/lZ6fridYt0y6UBfefJRt2zYa8erm/Xp6+ar1svCyornU0nZ4vSDQAxMDA6oV/W79MDz7eotXtY0/My9flrFujmS+aomBESAHGqtGCavvG+5fro5RX659836GuPNeiBja36/LUL9O4Vs+Vh08kZo3QDQIQMjE6orXtYLd1Dau0efv39tp5hHRoYlbXSxXMK9flrF+r6JTOVkcYICYDEsGhmvu67o1Ybd3fpa7/fqc89vF0/Xr9Xf/32RVpdWeI6XkKgdAPAGbLWqmtwXK1Tpbq1e0itPcOvv987PPFHzy/OzVSgKFuXzStSwJ+jKxeWaHl5oaP0AHD+Vs4v1qP/63L99uUD+ubju/SRezZpdWWx/vptF2jxrHzX8eKak9JtjLle0vckeSX92Fr79RM+nynpQUlVkrolfdBa2zL1ub+W9DFJIUmfsdY+HsPoAJJcKGx1oG9EbT1TV6m7j5bqY1esh8dDrz/XY6RZhdMUKMrW25aWKuDPVqAoW4GiHM3xZysnk+saAJKPx2N044rZuv7CmfrJ8636/rrdesf31+tPLpqtz1+7ULMLp7mOGJdi/hvBGOOV9ENJ10jaJ2mzMeZRa+1rxz3tY5J6rbXzjTE3SfpnSR80xiyWdJOkJZJmSXrKGLPAWhsSAJyBcNhqZCKkg/2jb7hi3dY9rPbeYU2E7OvPz0jzaI4/WwH/0SvWwaIczSk6+nGZL5sREQApKzPNq4+vnqv3V5XrX5/drfs2tOh3Lx/UHauC+tT/3969hdh11XEc//7nkkwmbZKa1KbG5tImLa0Rog6txRBSBbX6EAtS8lLzUEnABhT0QQUxiIIvqSCCktLaUmtrvQdU2hItKcVLU60mMVZDTGlDTGoLtoVcZub8fTg7kz2XXNrJnp195vuBYa+z9to5/4SVxY89+7J2OXNn9dZd4kWljtMwNwL7M/MAQEQ8AqwDyqF7HbClaP8E+E60X4u0DngkM08A/46I/cWf9/spql3SBZaZnBxucXywxYnBYY4Ptjg2OMzx4udY0XdiaJhjJ4v+odbIvhODrdK49tjj5faY404OtcbVcOnMHhbP7+f6K+fwkZULizPWs1kyv5+Fc/q8UUiSzmJufy9fuvV6PnXzUrY+/jzbdh7gR8+8yOZblnPHzUuY2eOTTqCe0L0IeLH0+SXgpjONycyhiPgfML/o/8OYYxdVV+pb94Uf/5Wjr5+ouwxdIJl57kGVfn+pTY7vy9H7xh8ztjHx2KT9d23l6W1r1OfTfTmyL2m1To8fzhx9bKs8doJjJ/FPO6Oni76eLmbN6Kavt5u+nm76ervo6+1mXv8M+nq7mNVb7Bv5ae9fOKePxfP7WTp/Npf19/q6Y0mapEXzZnH37au4c/Uyvvmbf/D1X+3j+08f5Jq3XzLltWxaczUfWL5gyr/3bDrygsOI2AhsBFi8eHEtNbxxfIjXjg2ee6Aao45MVv7KciiMkb7y2Bh3UJS2Ee1GlHpHHV+0uyKICLqi3e4KxnxuH9dV6hu1v6s9vnvUscX+rgmOBWb2djOzp2skGM8qheNyUJ5V/tzT7Rnoi9GTT9ZdgdQMHfx/5V3vmMuDd97Ezn++zD1PHaglD50cHv9bzbrVEboPAVeVPr+z6JtozEsR0QPMpX1D5fkcS2ZuA7YBDAwM1HKK8nt3vK+Or5UkSboorLn2ctZc6+MET6njDqBngBURsSwiZtC+MXL7mDHbgQ1F+5PAb7P9+/3twPqImBkRy4AVwJ+mqG5JkiTpLZnyM93FNdqbgcdoPzLwvszcGxFfA3ZlaoRarwAABHlJREFU5nbgXuDB4kbJV2kHc4pxj9K+6XIIuMsnl0iSJOliF3XfIFa1gYGB3LVrV91lSJIkqcNFxLOZOTDRPh8wK0mSJFXM0C1JkiRVzNAtSZIkVczQLUmSJFXM0C1JkiRVzNAtSZIkVczQLUmSJFXM0C1JkiRVzNAtSZIkVczQLUmSJFXM0C1JkiRVzNAtSZIkVczQLUmSJFXM0C1JkiRVzNAtSZIkVSwys+4aKhURLwMv1F2HLrgFwH/rLkKN5hzSZDmHNFnOoc6zJDMvn2hHx4dudaaI2JWZA3XXoeZyDmmynEOaLOfQ9OLlJZIkSVLFDN2SJElSxQzdaqptdRegxnMOabKcQ5os59A04jXdkiRJUsU80y1JkiRVzNCtRomIgxGxOyKei4hdddejZoiI+yLiaETsKfW9LSKeiIh/FdvL6qxRF68zzJ8tEXGoWIuei4iP1VmjLm4RcVVE/C4i/h4ReyPis0W/69A0YuhWE92Smat8zJLehPuBj47p+yKwIzNXADuKz9JE7mf8/AH4VrEWrcrMX09xTWqWIeDzmXkD8H7groi4AdehacXQLanjZeZO4NUx3euAB4r2A8AnprQoNcYZ5o903jLzcGb+uWi/DuwDFuE6NK0YutU0CTweEc9GxMa6i1GjXZGZh4v2f4Ar6ixGjbQ5Iv5WXH7iZQE6LxGxFHgP8Edch6YVQ7eaZnVmvhe4lfav59bUXZCaL9uPcfJRTnozvgtcA6wCDgNb6y1HTRARlwA/BT6Xma+V97kOdT5DtxolMw8V26PAz4Eb661IDXYkIq4EKLZHa65HDZKZRzJzODNbwD24FukcIqKXduB+KDN/VnS7Dk0jhm41RkTMjohLT7WBDwN7zn6UdEbbgQ1FewPwyxprUcOcCkqF23At0llERAD3Avsy8+7SLtehacSX46gxIuJq2me3AXqAH2bmN2osSQ0REQ8Da4EFwBHgq8AvgEeBxcALwO2Z6c1yGucM82ct7UtLEjgIbCpdmyuNEhGrgaeA3UCr6P4y7eu6XYemCUO3JEmSVDEvL5EkSZIqZuiWJEmSKmboliRJkipm6JYkSZIqZuiWJEmSKmboliRJkipm6JYkERG7I2LHBP3zImJ/RDwREd111CZJncDQLUkC+DbwwYi47lRH8Ra9H9B+GdX6zByuqzhJajpDtyQJ2uH6FeAzpb4twIeA2zLzlTqKkqROYeiWJJGZx4BtwIaI6I+IjwNfof1687/UW50kNZ+vgZckARARi4CDwFZgE/BQZm6utShJ6hCGbknSiIh4GFgPPA3ckpmDNZckSR3By0skSWWPFdtPG7gl6cIxdEuSyq4D3gCer7sQSeokhm5JUtlKYG967aEkXVCGbklS2buB3XUXIUmdxtAtSQIgIuYASzB0S9IFZ+iWJJ2ystjuqbUKSepAPjJQkiRJqphnuiVJkqSKGbolSZKkihm6JUmSpIoZuiVJkqSKGbolSZKkihm6JUmSpIoZuiVJkqSKGbolSZKkihm6JUmSpIr9H15iYk1UZMZSAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 864x576 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| } | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Статистический критерий\n", | |
| "\n", | |
| "Зафиксируем уровень статистической значимости $\\alpha=0.05$. То есть вероятность того, что мы отклоним нулевую гипотезу, при том что она истинна равна 5% или $P(\\mathrm{reject\\ H_0}|\\mathrm{H_0\\ true})=0.05$.\n", | |
| "\n", | |
| "Согласно распределению в условиях нулевой гипотезы $P(Y=k) = \\binom{21}{k}0.8^k(1-0.8)^{21-k}$\n", | |
| "\n", | |
| "$P(Y < 13)=0.043$\n", | |
| "\n", | |
| "$P(Y > 20)=0.009$\n", | |
| "\n", | |
| "Таким образом $P(\\mathrm{reject\\ H_0}|\\mathrm{H_0\\ true}) = P(Y < q_{\\frac{\\alpha}{2}}|\\mathrm{H_0\\ true}) + P(Y> q_{\\frac{\\alpha}{2}}|\\mathrm{H_0\\ true})=P(Y< 13|p=0.8) + P(Y> 20|p=0.8) \\approx 0.05$\n", | |
| "\n", | |
| "Определим критичесткую область: \n", | |
| "\n", | |
| "$R_{\\alpha}=\\{Y:\\ Y< 13\\ \\mathrm{or}\\ Y > 20\\}$\n", | |
| "\n", | |
| "Согласно данным только 19 пациентов ощутили эффект. Значит, $Y=19 \\notin R_{\\alpha}$. Следовательно, **нулевая гипотеза сохраняется**." | |
| ], | |
| "metadata": { | |
| "id": "mjFQ8E8TLIIZ" | |
| } | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "## Задача №3 A/B-testing _(асимтотический, односторонний тест с двумя выборками)_ \n", | |
| "\n", | |
| "При проектировании web-сайта для нового продукта вы хотите принять решение относительно местоположения кнопки для перехода к заказу услуги. В ваших предыдущих проектах кнопка заказа располагется в нижней части страницы, однако, сейчас вы предполагаете, что зафиксировав конпку в другом месте вы сможете повысить количество нажатий. \n", | |
| "\n", | |
| "Для сбора данных вы создаёте две версии сайта с разным местоположением кнопок. Далее, всем новым посетителям сервер случайно выбирает одну из двух версий сайта. Выполняя эксперимент два дня вы собрали следующие данные:\n", | |
| "\n", | |
| "> Общее количество посетителей с традицонным расположением кнопки: 49, количество переходов: 28. Общее количество посетителей с альтернативным расположением элемента дизайна: 51, количество переходов 38. \n", | |
| "\n", | |
| "\n", | |
| "Принимать ли новый дизайн в производство?\n", | |
| "\n", | |
| "\n" | |
| ], | |
| "metadata": { | |
| "id": "YQWAEP1wmCCw" | |
| } | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Формализация задачи\n", | |
| "\n", | |
| "По условию задачи мы имеем две выборки ($X$ - традиционная кнопка, $Y$ - альтернативная), где $\\{X_i\\ \\mathrm{or}\\ Y_i\\}=1$ если пользователь нажал на кнопку при просмотре страницы и $0$ в противном случае.\n", | |
| "\n", | |
| "Пусть $N=49$ и $M=51$ развер выборок для традиционного и альтернативного расположения кнопок. Тогда, $\\bar{X}_N=\\frac{1}{n}\\sum_{i=1}^n X_i$ и $\\bar{Y}_M=\\frac{1}{n}\\sum_{i=1}^M Y_i$ это выборочные средние для каждого наблюдения. \n", | |
| "\n", | |
| "По центральной предельной теореме $\\bar{X}_N \\sim N(p_X, \\frac{p_X(1-p_X)}{N})$ и $\\bar{Y}_M \\sim N(p_Y, \\frac{p_Y(1-p_Y)}{M})$\n", | |
| "\n", | |
| "Для решения задачи необходимо выполить тест:\n", | |
| "\n", | |
| "$H_0:\\ p_X=p_Y\\ \\mathrm{VS}\\ H_A: p_Y > p_X$, где $p$ это параметр биномиального распределения." | |
| ], | |
| "metadata": { | |
| "id": "UvxEjxgqqnEC" | |
| } | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Статистический критерий\n", | |
| "\n", | |
| "Зафиксируем уровень статистической значимости $\\alpha=0.05$. То есть вероятность того, что мы отклоним нулевую гипотезу, при том что она истинна равна 5% или $P(\\mathrm{reject\\ H_0}|\\mathrm{H_0\\ true})=0.05$.\n", | |
| "\n", | |
| "Известно, что $\\mathbb{E}[\\bar{X}_N]=p_X$ и $\\mathbb{E}[\\bar{Y}_M]=p_Y$.\n", | |
| "\n", | |
| "Кроме этого, функция максимального правдоподобия, для двух выборок будет иметь вид:\n", | |
| "\n", | |
| "$L(X_1,...,X_N, Y_1,...,Y_M;\\ p_X,p_Y) =p_X^{\\sum_{i=1}^N X_i}(1-p_X)^{N-\\sum_{i=1}^N X_i} p_Y^{\\sum_{i=1}^M Y_i}(1-p_Y)^{M-\\sum_{i=1}^M Y_i}$\n", | |
| "\n", | |
| "У условиях нулевой гипотезы $p_X=p_Y$, тогда, лог-функция максимального правдоподобия будет выражена:\n", | |
| "\n", | |
| "$(\\sum_{i=1}^N X_i + \\sum_{i=1}^M Y_i)\\log(p)+(N+M-\\sum_{i=1}^N X_i-\\sum_{i=1}^M Y_i)\\log(1-p)$\n", | |
| "\n", | |
| "Минимизация лог-функции максимального правдоподобия относительно $p$ даёт точный результат: $\\hat{p}=\\frac{\\sum_{i=1}^N X_i+\\sum_{i=1}^M Y_i}{N+M}$\n", | |
| "\n", | |
| "Тогда, в условиях нулевой гипотезы, мы имеем следующую тест статистику:\n", | |
| "\n", | |
| "$\\frac{(\\bar{Y}_M - \\bar{Y}_M) - (p_Y - p_X)}{\\sqrt{\\frac{\\hat{p}(1-\\hat{p})}{M} + \\frac{\\hat{p}(1-\\hat{p}}{N}}}=\\sqrt{NM}\\frac{\\bar{Y}_M - \\bar{Y}_M}{\\sqrt{(N+M)(\\hat{p}(1-\\hat{p})}} \\sim N(0,1)$\n", | |
| "\n", | |
| "Таким образом, статистический критерий выражен:\n", | |
| "\n", | |
| "$$\\psi_\\alpha=\\mathbb{1}\\left[\\sqrt{NM}\\frac{\\bar{Y}_M - \\bar{Y}_M}{\\sqrt{(N+M)(\\hat{p}(1-\\hat{p})}} > q_{\\alpha}=1.65 \\right] $$" | |
| ], | |
| "metadata": { | |
| "id": "xqIlKZ8KqydV" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# data\n", | |
| "alpha=0.05\n", | |
| "N = 49\n", | |
| "M = 51\n", | |
| "sum_x = 28\n", | |
| "sum_y = 38\n", | |
| "\n", | |
| "# quantile\n", | |
| "q_alpha = norm.ppf(0.95)\n", | |
| "print('Quantile 5% is {}'.format(np.round(q_alpha, decimals=3)))\n", | |
| "\n", | |
| "# p under null\n", | |
| "mle_p = (sum_x + sum_y)/(N+M)\n", | |
| "print('p under null is {}'.format(np.round(mle_p, decimals=3)))\n", | |
| "\n", | |
| "# test statistic\n", | |
| "test_value = np.sqrt(N*M)*(sum_y/M - sum_x/N)*(1/np.sqrt((N+M)*(mle_p*(1-mle_p))))\n", | |
| "print('test statistic is {}'.format(np.round(test_value, decimals=3)))\n", | |
| "\n", | |
| "# decision\n", | |
| "print(\n", | |
| " 'Is the null hypothesis rejected? / {}'.format(\n", | |
| " test_value > q_alpha\n", | |
| " )\n", | |
| " )\n", | |
| "\n", | |
| "# p-value\n", | |
| "print(np.round(1-norm.cdf(test_value), decimals=2))" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "qGpD5zRr-ro0", | |
| "outputId": "c0f83d1e-fcea-4ecf-bcd8-8c4d8d908b15" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "Quantile 5% is 1.645\n", | |
| "p under null is 0.66\n", | |
| "test statistic is 1.833\n", | |
| "Is the null hypothesis rejected? / True\n", | |
| "0.03\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "source": [ | |
| "### Альтернативное решение - точный тест Фишера\n", | |
| "\n", | |
| "Простое [руководство](https://www.sheffield.ac.uk/polopoly_fs/1.43998!/file/tutorial-9-fishers.pdf) для двух выборок.\n", | |
| "\n", | |
| "Нулевая гипотеза: **никакого эффекта от изменений нет**.\n", | |
| "\n", | |
| "| |Алтернативная кнопка | Традиционная кнопка | Итого_______ \n", | |
| "-------------------|-------------------|------------------|------------------\n", | |
| "| Переходы |38 ($a$) | 28 ($b$) | 66 ($a+b$) \n", | |
| "| Отсутствие перехода |13 ($c$) | 21 ($d$) | 34 ($c+d$)\n", | |
| "| Итого | 51 ($a+c$) | 49 ($b+d$) | 100 ($n$)\n", | |
| "\n", | |
| "Тогда вероятность получить данные как в таблице, в условиях нулевой гипотезы выражены:\n", | |
| "\n", | |
| "$$P=\\frac{(a+b)!(c+d)!(a+c)!(b+d)!}{a!b!c!d!n!}$$\n", | |
| "\n", | |
| "P-value в условиях точного теста Фишера, это сумма вероятностей более экстримальных вариантов, т.е.: $\\mathrm{P-value}=P(a=38,b=27)+P(a=39,b=26)+...+P(a=51,b=15)$\n" | |
| ], | |
| "metadata": { | |
| "id": "MHkC-yeKnuA8" | |
| } | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# a=38, b=28\n", | |
| "(\n", | |
| " np.math.factorial(66)*np.math.factorial(34)*np.math.factorial(51)*np.math.factorial(49)\n", | |
| " )/(\n", | |
| " np.math.factorial(38)*np.math.factorial(28)*np.math.factorial(13)*np.math.factorial(21)*np.math.factorial(100)\n", | |
| " )" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "l8yHY7mxp30w", | |
| "outputId": "5b2b58a9-630e-4010-b863-be0f433b9a22" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "0.03202576668586309" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 50 | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# a=39, b=27\n", | |
| "(\n", | |
| " np.math.factorial(66)*np.math.factorial(34)*np.math.factorial(51)*np.math.factorial(49)\n", | |
| " )/(\n", | |
| " np.math.factorial(39)*np.math.factorial(27)*np.math.factorial(12)*np.math.factorial(22)*np.math.factorial(100)\n", | |
| " )" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "G4dBZq7z392q", | |
| "outputId": "962843fd-7e29-4b93-a5a2-08fcee9fc9de" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "0.013586688897032826" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 51 | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# a=40, b=26\n", | |
| "(\n", | |
| " np.math.factorial(66)*np.math.factorial(34)*np.math.factorial(51)*np.math.factorial(49)\n", | |
| " )/(\n", | |
| " np.math.factorial(40)*np.math.factorial(26)*np.math.factorial(11)*np.math.factorial(23)*np.math.factorial(100)\n", | |
| " )" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "uQn0TtRm4YMv", | |
| "outputId": "a64ce50c-33a1-4049-9c40-38733226ae0e" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "0.00478487739417243" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 52 | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# a=41, b=25\n", | |
| "(\n", | |
| " np.math.factorial(66)*np.math.factorial(34)*np.math.factorial(51)*np.math.factorial(49)\n", | |
| " )/(\n", | |
| " np.math.factorial(41)*np.math.factorial(25)*np.math.factorial(10)*np.math.factorial(24)*np.math.factorial(100)\n", | |
| " )" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "Z_m_7_fW4p0Q", | |
| "outputId": "58f48e9d-b0c3-4c09-9cad-93263495368b" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "0.0013907265596883284" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 53 | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "# a=42, b=24\n", | |
| "(\n", | |
| " np.math.factorial(66)*np.math.factorial(34)*np.math.factorial(51)*np.math.factorial(49)\n", | |
| " )/(\n", | |
| " np.math.factorial(42)*np.math.factorial(24)*np.math.factorial(9)*np.math.factorial(25)*np.math.factorial(100)\n", | |
| " )" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "Z2HCCLkS40tl", | |
| "outputId": "a2a15bb2-ca97-4e5f-f428-6352abf2b6bb" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "0.00033112537135436386" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 54 | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "probs = []\n", | |
| "k=0\n", | |
| "for i in range(38, 51, 1):\n", | |
| " probs.append((\n", | |
| " np.math.factorial(66)*np.math.factorial(34)*np.math.factorial(51)*np.math.factorial(49)\n", | |
| " )/(\n", | |
| " np.math.factorial(i)*np.math.factorial(28-k)*np.math.factorial(51-i)*np.math.factorial(21+k)*np.math.factorial(100)\n", | |
| " ))\n", | |
| " k=k+1" | |
| ], | |
| "metadata": { | |
| "id": "RAw-Qfaw6jTP" | |
| }, | |
| "execution_count": null, | |
| "outputs": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "sum(probs)" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "aXFm6v6LJA2D", | |
| "outputId": "198b8c3e-7821-46ec-dbaa-985c93cb2c03" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "execute_result", | |
| "data": { | |
| "text/plain": [ | |
| "0.05219440145192389" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "execution_count": 56 | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "table = np.array([[38,28],[13,21]])\n", | |
| "oddsr, p = fisher_exact(table, alternative='greater')\n", | |
| "print(p)" | |
| ], | |
| "metadata": { | |
| "colab": { | |
| "base_uri": "https://localhost:8080/" | |
| }, | |
| "id": "Tp7Whp9pta6d", | |
| "outputId": "9d416103-a9fe-456f-e270-aa82cd6b51e3" | |
| }, | |
| "execution_count": null, | |
| "outputs": [ | |
| { | |
| "output_type": "stream", | |
| "name": "stdout", | |
| "text": [ | |
| "0.052194401451926706\n" | |
| ] | |
| } | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "source": [ | |
| "" | |
| ], | |
| "metadata": { | |
| "id": "Y1WaZYNLwrdQ" | |
| }, | |
| "execution_count": null, | |
| "outputs": [] | |
| } | |
| ] | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment