Example2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From: https://github.com/ksatola/Computational-Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Is a drug effective or not?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "drug = [54, 73, 53, 70, 73, 68, 52, 65, 65] # treatment group\n",
    "placebo = [54, 51, 58, 44, 55, 52, 42, 47, 58, 46] # control group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "H_0 : \\mu_{drug} - \\mu_{placebo} = 0 \\\\\n",
    "H_1 : \\mu_{drug} - \\mu_{placebo} <> 0\n",
    "$$\n",
    "Where:\n",
    "\n",
    "$ H_0 $ -> **The null hypothesis:** there is no difference in the drug and placebo groups on average (the drug does not have any effect). The null hypothesis should always contain `=` sign.\n",
    "\n",
    "$ H_1 $ -> **The alternative hypothesis:** there is a significant difference in the drug and placebo groups on average (the drug works).\n",
    "\n",
    "I will assume **Type I Error / a significance level (alpha)** at 0.01.\n",
    "\n",
    "Remarks:\n",
    "- **Type I Error.** Reject the null hypothesis when there is in fact no significant effect (false positive). The p-value is optimistically small.\n",
    "- **Type II Error.** Not reject the null hypothesis when there is a significant effect (false negative). The p-value is pessimistically large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Significance level\n",
    "alpha = 0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.966666666666661"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Means difference of the samples\n",
    "np.mean(drug) - np.mean(placebo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the means difference statistically significant?\n",
    "\n",
    "We use **T-Student sampling distribution** and calculate **t-test** statistic to find out.\n",
    "\n",
    "**T-test** assumptions:\n",
    "1. Both groups (drug and placebo) should be normally distributed\n",
    "- Samples should come from populations with equal variances\n",
    "- Samples should be of the same size\n",
    "- At best there should be more than 20 observations in a sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assumption 1\n",
    "Both groups (drug and placebo) should be normally distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEetJREFUeJzt3X2Q1dV9x/H3110sKNpEWIwJItBQMQqCs9BkNBilQRIV60QzIYljO3bQGUwUa9LEjDP6h52O04fEMQ+zmgTSmERr65hYySMyRosPkKAxYBpDUdcYWbBBsaIg3/6xF4J0l/2t3oc97Ps1s7P34fD7fc+e8eO55577u5GZSJLKcVCrC5AkDY7BLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSpMeyMOOnbs2Jw4cWIjDi1JB6Q1a9ZszsyOKm0bEtwTJ05k9erVjTi0JB2QIuLJqm1dKpGkwhjcklQYg1uSCtOQNe6+7Nixg+7ubrZv396sUxZn5MiRjB8/nhEjRrS6FElDWNOCu7u7m8MOO4yJEycSEc06bTEyky1bttDd3c2kSZNaXY6kIazSUklEvCUibo+IxyNifUS8Z7An2r59O2PGjDG0+xERjBkzxlckkgZUdcb9BeD7mXleRBwMHPJGTmZo759/H0lVDBjcEXE4MAf4S4DMfBV4tbFlSZL6U2XGPRnoAb4eEScCa4DLMvOlN3PitZvruyQwY+zIAdu0tbUxbdo0duzYQXt7OxdeeCGXX345Bx3k5hpJ5agS3O3AScAnMvPBiPgC8Bng6r0bRcQiYBHAhAkT6l1nXYwaNYq1a9cCsGnTJj760Y+ydetWrr322te127lzJ+3tTXvfVhqUek96VD9VJpD1UGWq2Q10Z+aDtfu30xvkr5OZXZnZmZmdHR2VPm7fUuPGjaOrq4sbb7yRzGTp0qWcf/75nH322cybN4+VK1dy1lln7Wl/6aWXsnTpUgDuvvtupk6dyimnnMInP/nJ17WTpEYbMLgz83fA0xFxbO2hucC6hlbVJJMnT2bXrl1s2rQJgFWrVrFs2TJWrFjR77/Zvn07F198McuXL+e+++6jp6enWeVKElD9k5OfAG6JiEeBGcDfNa6k5srMPbff//73c8QRR+y3/eOPP87kyZP37LVeuHBhQ+uTpH1VWsjNzLVAZ4NraboNGzbQ1tbGuHHjADj00EP3PNfe3s6uXbv23N+9v3rvoJekVhi22yl6enq45JJLuPTSS/vcP33MMcewbt06XnnlFbZu3cpPfvITAKZOncqGDRvYuHEjALfeemszy5ak5n3kfV/Nevd1by+//DIzZszYsx3wggsu4Iorruiz7dFHH82HP/xhpk+fzpQpU5g5cybQuzPlS1/6EvPnz2fs2LHMnj27mV2QJKIRL/07Oztz3y9SWL9+Pccdd1zdz9UK27ZtY/To0WQmixcvZsqUKSxZsqQuxz6Q/k5qDLcDDl1vZkIaEWsys9KS9LBdKnkzbrrpJmbMmMHxxx/P1q1bufjii1tdkqRhxE+ZvAFLliyp2wxbkgbLGbckFcbglqTCGNySVBiDW5IK07o3J7u66nu8RYsGbLL7sq47d+7kuOOOY9myZRxyyCGMHj2abdu21a2Ua665htGjR3PllVfW7ZiStNuwmnHvvqzrY489xsEHH8xXvvKVVpckSYM2rIJ7b+9973t54oknXvfYtm3bmDt3LieddBLTpk3jzjvv3PPcN77xDaZPn86JJ57IBRdcAPR+bP5DH/oQs2bNYtasWdx///172j/yyCOcfvrpTJkyhZtuugnovc7Jpz71KU444QSmTZvmx+UlvSHDch/3zp07Wb58OfPnz3/d4yNHjuSOO+7g8MMPZ/Pmzbz73e9mwYIFrFu3juuuu47777+fsWPH8vzzzwNw2WWXsWTJEk455RSeeuopzjjjDNavXw/Ao48+ygMPPMBLL73EzJkzOfPMM1m1ahVr167lkUceYfPmzcyaNYs5c+Zw1FFHNf1vIKlcwyq4d1+rBHpn3BdddNHrns9MrrrqKu69914OOuggnnnmGZ577jlWrFjBeeedx9ixYwH2XPr1xz/+MevW/eHS5C+88AIvvvgiAOeccw6jRo1i1KhRnHbaaTz00EPcd999LFy4kLa2No488khOPfVUHn74YRYsWNCM7ks6QAyr4N77q8v6csstt9DT08OaNWsYMWIEEydOZPv27WRmn1cQ3LVrF6tWrWLUqFH/77l920eEl4SVVBfDdo27L1u3bmXcuHGMGDGCe+65hyeffBKAuXPnctttt7FlyxaAPUsl8+bN48Ybb9zz7/f+n8Kdd97J9u3b2bJlCytXrtyzLHLrrbfy2muv0dPTw7333uvVBSUNWutm3BW27zXbxz72Mc4++2w6OzuZMWMGU6dOBeD444/nc5/7HKeeeiptbW3MnDmTpUuXcsMNN7B48WKmT5/Ozp07mTNnzp6dKrNnz+bMM8/kqaee4uqrr+btb3875557LqtWreLEE08kIrj++ut529ve1souSyqQl3UdYvw7aSBe1nXo8rKukqQ+GdySVJimBre7KvbPv4+kKpoW3CNHjmTLli2GUz8yky1btjByZPO/i1NSWZq2q2T8+PF0d3fT09PTrFMWZ+TIkYwfP77VZUga4poW3CNGjGDSpEnNOp0kHbB8c1KSCmNwS1JhKi2VRMRG4EXgNWBn1U3ikqT6G8wa92mZublhlUiSKnGpRJIKUzW4E/hhRKyJiD6vDhURiyJidUSsdsufJDVO1eA+OTNPAj4ALI6IOfs2yMyuzOzMzM6Ojo66FilJ+oNKwZ2Zv6393gTcAXgRaUlqkQGDOyIOjYjDdt8G5gGPNbowSVLfquwqORK4o/ZVXO3AtzLz+w2tSpLUrwGDOzM3ACc2oRZJUgVuB5SkwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwlQO7ohoi4ifR8RdjSxIkrR/g5lxXwasb1QhkqRqKgV3RIwHzgRubmw5kqSBVJ1xfx74NLCrgbVIkipoH6hBRJwFbMrMNRHxvv20WwQsApgwYcIbr6ir643/26Fo0aJWV1B/jpHUUlVm3CcDCyJiI/Ad4PSI+Oa+jTKzKzM7M7Ozo6OjzmVKknYbMLgz87OZOT4zJwIfAVZk5scbXpkkqU/u45akwgy4xr23zFwJrGxIJZKkSpxxS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKsyAwR0RIyPioYh4JCJ+GRHXNqMwSVLf2iu0eQU4PTO3RcQI4L6IWJ6ZDzS4NklSHwYM7sxMYFvt7ojaTzayKElS/6rMuImINmAN8E7gi5n5YB9tFgGLACZMmFDPGtUCazdv7/e5Mdt2NLGSxtuyn75KQ1GlNycz87XMnAGMB2ZHxAl9tOnKzM7M7Ozo6Kh3nZKkmkHtKsnM3wMrgfkNqUaSNKAqu0o6IuIttdujgD8HHm90YZKkvlVZ4z4KWFZb5z4IuC0z72psWZKk/lTZVfIoMLMJtUiSKvCTk5JUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUmAGDOyKOjoh7ImJ9RPwyIi5rRmGSpL61V2izE/ibzPxZRBwGrImIH2XmugbXJknqw4Az7sx8NjN/Vrv9IrAeeEejC5Mk9W1Qa9wRMRGYCTzYiGIkSQOrslQCQESMBv4NuDwzX+jj+UXAIoAJEybUrcDidXW1uoI3ZMy2Ha0uQVI/Ks24I2IEvaF9S2b+e19tMrMrMzszs7Ojo6OeNUqS9lJlV0kAXwXWZ+Y/Nb4kSdL+VJlxnwxcAJweEWtrPx9scF2SpH4MuMadmfcB0YRaJEkV+MlJSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhRkwuCPiaxGxKSIea0ZBkqT9qzLjXgrMb3AdkqSKBgzuzLwXeL4JtUiSKqjbGndELIqI1RGxuqenp16HlSTto27BnZldmdmZmZ0dHR31OqwkaR/uKpGkwhjcklSYKtsBvw2sAo6NiO6IuKjxZUmS+tM+UIPMXNiMQiRJ1bhUIkmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMEMvuL/3vaF5rMEcr97nraeKtb384MOV2j068fi6He/wH/xHpWPV2zvPmdeS8w5lDz73v00/56+3vjKkj1fFLf/1+6acZ+gF9113Dc1jDeZ49T5vPVWs7U9vW1ap3S8mnVC34/3xD5dXOla9jf7Pn7bkvEPZ/7y6q+nn/M0LO4b08ap4+qWdTTnP0AtuSdJ+GdySVBiDW5IKY3BLUmEMbkkqTKXgjoj5EfGriHgiIj7T6KIkSf0bMLgjog34IvAB4F3Awoh4V6MLkyT1rcqMezbwRGZuyMxXge8A5zS2LElSf6oE9zuAp/e63117TJLUApGZ+28QcT5wRmb+de3+BcDszPzEPu0WAYtqd48FflX/chtmLLC51UW8SfZh6DgQ+mEfmu+YzOyo0rC9Qptu4Oi97o8Hfrtvo8zsAroqlTfERMTqzOxsdR1vhn0YOg6EftiHoa3KUsnDwJSImBQRBwMfAb7b2LIkSf0ZcMadmTsj4lLgB0Ab8LXM/GXDK5Mk9anKUgmZeTdwd4NraaUil3j2YR+GjgOhH/ZhCBvwzUlJ0tDiR94lqTDDMrgjoi0ifh4Rd9XuL42I/46ItbWfGa2ucSARsTEiflGrd3XtsSMi4kcR8eva77e2us796acP10TEM3uNxQdbXef+RMRbIuL2iHg8ItZHxHsKHIe++lDMOETEsXvVuTYiXoiIy0sbh8EYlkslEXEF0AkcnplnRcRS4K7MvL21lVUXERuBzszcvNdj1wPPZ+bf164p89bM/NtW1TiQfvpwDbAtM/+hVXUNRkQsA36amTfXdl0dAlxFWePQVx8up6Bx2K12iY5ngD8DFlPQOAzGsJtxR8R44Ezg5lbX0gDnALu/I2wZ8BctrOWAFxGHA3OArwJk5quZ+XsKGof99KFUc4HfZOaTFDQOgzXsghv4PPBpYN8v1bsuIh6NiH+OiD9qQV2DlcAPI2JN7VOrAEdm5rMAtd/jWlZdNX31AeDS2lh8bYi/vJ0M9ABfry293RwRh1LWOPTXByhnHPb2EeDbtdsljcOgDKvgjoizgE2ZuWafpz4LTAVmAUcAJbycOjkzT6L3qo2LI2JOqwt6A/rqw5eBPwFmAM8C/9jC+gbSDpwEfDkzZwIvAaVd9ri/PpQ0DgDUlnkWAP/a6loabVgFN3AysKC2tvod4PSI+GZmPpu9XgG+Tu8VEYe0zPxt7fcm4A56a34uIo4CqP3e1LoKB9ZXHzLzucx8LTN3ATcxtMeiG+jOzAdr92+nNwRLGoc++1DYOOz2AeBnmflc7X5J4zAowyq4M/OzmTk+MyfS+5JqRWZ+fK/BDXrXwR5rYZkDiohDI+Kw3beBefTW/F3gwlqzC4E7W1PhwPrrw+6xqDmXITwWmfk74OmIOLb20FxgHQWNQ399KGkc9rKQPyyTQEHjMFjDclcJQES8D7iytqtkBdABBLAWuCQzt7Wyvv2JiMn0zlCh96XutzLzuogYA9wGTACeAs7PzOdbVOZ+7acP/0Lvy/MENgIX716nHIpqW0dvBg4GNgB/Re+EqIhxgH77cANljcMh9F5+enJmbq09Vsx/D4M1bINbkko1rJZKJOlAYHBLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklSY/wNpVFf+DPbmEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "\n",
    "sns.distplot(drug, color=\"skyblue\", label=\"Drug\", kde=False, hist=True, rug=\"True\")\n",
    "sns.distplot(placebo, color=\"red\", label=\"Placebo\", kde=False, hist=True, rug=\"True\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "def normaltest(sample, name, alpha=0.05):\n",
    "    '''\n",
    "    Shapiro-Wilk test for samples with n < 50 observations.\n",
    "    Reference: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html\n",
    "    '''\n",
    "    W, p = stats.shapiro(sample)\n",
    "    #print(\"W = {:g}\".format(W))\n",
    "    #print(\"p = {:g}\".format(p))\n",
    "    if p < alpha:  # null hypothesis: sample comes from a normal distribution\n",
    "        print(\"{} does not come from a normal distribution\".format(name))\n",
    "    else:\n",
    "        print(\"{} comes from a normal distribution\".format(name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drug comes from a normal distribution\n"
     ]
    }
   ],
   "source": [
    "# Is drug normally distributed?\n",
    "normaltest(drug, \"Drug\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Placebo comes from a normal distribution\n"
     ]
    }
   ],
   "source": [
    "# Is placebo normally distributed?\n",
    "normaltest(placebo, \"Placebo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To solve the issue of having too small sample size...\n",
    "# I could use the following:  (just kidding!)\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assumption 2\n",
    "Samples should come from populations with equal variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def varsequals(sample1, sample2, alpha=0.01):\n",
    "    '''\n",
    "    Perform Bartlett’s test for equal variances.\n",
    "    Bartlett’s test tests the null hypothesis that all input samples are from populations with equal variances.\n",
    "    Reference: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.bartlett.html\n",
    "    ''' \n",
    "    T, p = stats.bartlett(sample1, sample2)\n",
    "    #print(\"T = {:g}\".format(T))\n",
    "    #print(\"p = {:g}\".format(p))\n",
    "    if p < alpha:  # null hypothesis: all input samples are from populations with equal variances\n",
    "        print(\"Not all input samples are from populations with equal variances.\")\n",
    "    else:\n",
    "        print(\"All input samples are from populations with equal variances.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All input samples are from populations with equal variances.\n"
     ]
    }
   ],
   "source": [
    "varsequals(drug, placebo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate t-statistic\n",
    "Are our samples expected to have been drawn from the same population?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t=3.9357, p=0.0011\n"
     ]
    }
   ],
   "source": [
    "# p-value: Probability of obtaining a result equal to or more extreme than was observed in the data.\n",
    "t_stat, p = stats.ttest_ind(drug, placebo) # all observations are independent\n",
    "print('t={:.4f}, p={:.4f}'.format(t_stat, p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The p-value returned by the t-test is two-tailed by design (if we were testing one-tailed hypothesis, is such case we would have to divide the p-value by 2). The p-value is far lower than our alpha, so we can conclude we can reject the null hypothesis.\n",
    "\n",
    "If the p-value were not as small, we should calculate the critical value (below):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate degrees of freedom\n",
    "df = len(drug) + len(placebo) - 2\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.01"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.8982305196347173\n"
     ]
    }
   ],
   "source": [
    "# Calculate the critical value (two-tailed test)\n",
    "# PPF (percent point function)\n",
    "cv = stats.t.ppf(1.0 - alpha/2, df)\n",
    "print(cv)\n",
    "\n",
    "# Confirm cv with CDF (cumulative distribution function)\n",
    "#p = stats.t.cdf(cv, df)\n",
    "#print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"images/t_table.jpeg\" width=\"70%\"/></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reject the null hypothesis that the means are equal.\n"
     ]
    }
   ],
   "source": [
    "# Interpret via critical value (abs for symmetric distribution)\n",
    "if abs(t_stat) <= cv:\n",
    "    print('Accept null hypothesis that the means are equal.')\n",
    "else:\n",
    "    print('Reject the null hypothesis that the means are equal.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reject the null hypothesis that the means are equal.\n"
     ]
    }
   ],
   "source": [
    "# Interpret via p-value\n",
    "if p < alpha:\n",
    "    print('Reject the null hypothesis that the means are equal.')\n",
    "else:\n",
    "    print('Accept null hypothesis that the means are equal.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hypothesis test result\n",
    "We have sufficient evidence to reject **the null hypothesis** which means the drug works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Permutation Testing (Computational Method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Python's statitics module provides functions \n",
    "# for calculating mathematical statistics of numeric (real-valued) data\n",
    "from statistics import mean\n",
    "from random import shuffle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m\n",
       "Return the sample arithmetic mean of data.\n",
       "\n",
       ">>> mean([1, 2, 3, 4, 4])\n",
       "2.8\n",
       "\n",
       ">>> from fractions import Fraction as F\n",
       ">>> mean([F(3, 7), F(1, 21), F(5, 3), F(1, 3)])\n",
       "Fraction(13, 21)\n",
       "\n",
       ">>> from decimal import Decimal as D\n",
       ">>> mean([D(\"0.5\"), D(\"0.75\"), D(\"0.625\"), D(\"0.375\")])\n",
       "Decimal('0.5625')\n",
       "\n",
       "If ``data`` is empty, StatisticsError will be raised.\n",
       "\u001b[0;31mFile:\u001b[0m      ~/anaconda3/lib/python3.7/statistics.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "63.666666666666664\n",
      "50.7\n",
      "12.97\n"
     ]
    }
   ],
   "source": [
    "observed_diff = mean(drug) - mean(placebo)\n",
    "\n",
    "print(mean(drug))\n",
    "print(mean(placebo))\n",
    "\n",
    "print(f\"{observed_diff:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEetJREFUeJzt3X2Q1dV9x/H3110sKNpEWIwJItBQMQqCs9BkNBilQRIV60QzIYljO3bQGUwUa9LEjDP6h52O04fEMQ+zmgTSmERr65hYySMyRosPkKAxYBpDUdcYWbBBsaIg3/6xF4J0l/2t3oc97Ps1s7P34fD7fc+e8eO55577u5GZSJLKcVCrC5AkDY7BLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSpMeyMOOnbs2Jw4cWIjDi1JB6Q1a9ZszsyOKm0bEtwTJ05k9erVjTi0JB2QIuLJqm1dKpGkwhjcklQYg1uSCtOQNe6+7Nixg+7ubrZv396sUxZn5MiRjB8/nhEjRrS6FElDWNOCu7u7m8MOO4yJEycSEc06bTEyky1bttDd3c2kSZNaXY6kIazSUklEvCUibo+IxyNifUS8Z7An2r59O2PGjDG0+xERjBkzxlckkgZUdcb9BeD7mXleRBwMHPJGTmZo759/H0lVDBjcEXE4MAf4S4DMfBV4tbFlSZL6U2XGPRnoAb4eEScCa4DLMvOlN3PitZvruyQwY+zIAdu0tbUxbdo0duzYQXt7OxdeeCGXX345Bx3k5hpJ5agS3O3AScAnMvPBiPgC8Bng6r0bRcQiYBHAhAkT6l1nXYwaNYq1a9cCsGnTJj760Y+ydetWrr322te127lzJ+3tTXvfVhqUek96VD9VJpD1UGWq2Q10Z+aDtfu30xvkr5OZXZnZmZmdHR2VPm7fUuPGjaOrq4sbb7yRzGTp0qWcf/75nH322cybN4+VK1dy1lln7Wl/6aWXsnTpUgDuvvtupk6dyimnnMInP/nJ17WTpEYbMLgz83fA0xFxbO2hucC6hlbVJJMnT2bXrl1s2rQJgFWrVrFs2TJWrFjR77/Zvn07F198McuXL+e+++6jp6enWeVKElD9k5OfAG6JiEeBGcDfNa6k5srMPbff//73c8QRR+y3/eOPP87kyZP37LVeuHBhQ+uTpH1VWsjNzLVAZ4NraboNGzbQ1tbGuHHjADj00EP3PNfe3s6uXbv23N+9v3rvoJekVhi22yl6enq45JJLuPTSS/vcP33MMcewbt06XnnlFbZu3cpPfvITAKZOncqGDRvYuHEjALfeemszy5ak5n3kfV/Nevd1by+//DIzZszYsx3wggsu4Iorruiz7dFHH82HP/xhpk+fzpQpU5g5cybQuzPlS1/6EvPnz2fs2LHMnj27mV2QJKIRL/07Oztz3y9SWL9+Pccdd1zdz9UK27ZtY/To0WQmixcvZsqUKSxZsqQuxz6Q/k5qDLcDDl1vZkIaEWsys9KS9LBdKnkzbrrpJmbMmMHxxx/P1q1bufjii1tdkqRhxE+ZvAFLliyp2wxbkgbLGbckFcbglqTCGNySVBiDW5IK07o3J7u66nu8RYsGbLL7sq47d+7kuOOOY9myZRxyyCGMHj2abdu21a2Ua665htGjR3PllVfW7ZiStNuwmnHvvqzrY489xsEHH8xXvvKVVpckSYM2rIJ7b+9973t54oknXvfYtm3bmDt3LieddBLTpk3jzjvv3PPcN77xDaZPn86JJ57IBRdcAPR+bP5DH/oQs2bNYtasWdx///172j/yyCOcfvrpTJkyhZtuugnovc7Jpz71KU444QSmTZvmx+UlvSHDch/3zp07Wb58OfPnz3/d4yNHjuSOO+7g8MMPZ/Pmzbz73e9mwYIFrFu3juuuu47777+fsWPH8vzzzwNw2WWXsWTJEk455RSeeuopzjjjDNavXw/Ao48+ygMPPMBLL73EzJkzOfPMM1m1ahVr167lkUceYfPmzcyaNYs5c+Zw1FFHNf1vIKlcwyq4d1+rBHpn3BdddNHrns9MrrrqKu69914OOuggnnnmGZ577jlWrFjBeeedx9ixYwH2XPr1xz/+MevW/eHS5C+88AIvvvgiAOeccw6jRo1i1KhRnHbaaTz00EPcd999LFy4kLa2No488khOPfVUHn74YRYsWNCM7ks6QAyr4N77q8v6csstt9DT08OaNWsYMWIEEydOZPv27WRmn1cQ3LVrF6tWrWLUqFH/77l920eEl4SVVBfDdo27L1u3bmXcuHGMGDGCe+65hyeffBKAuXPnctttt7FlyxaAPUsl8+bN48Ybb9zz7/f+n8Kdd97J9u3b2bJlCytXrtyzLHLrrbfy2muv0dPTw7333uvVBSUNWutm3BW27zXbxz72Mc4++2w6OzuZMWMGU6dOBeD444/nc5/7HKeeeiptbW3MnDmTpUuXcsMNN7B48WKmT5/Ozp07mTNnzp6dKrNnz+bMM8/kqaee4uqrr+btb3875557LqtWreLEE08kIrj++ut529ve1souSyqQl3UdYvw7aSBe1nXo8rKukqQ+GdySVJimBre7KvbPv4+kKpoW3CNHjmTLli2GUz8yky1btjByZPO/i1NSWZq2q2T8+PF0d3fT09PTrFMWZ+TIkYwfP77VZUga4poW3CNGjGDSpEnNOp0kHbB8c1KSCmNwS1JhKi2VRMRG4EXgNWBn1U3ikqT6G8wa92mZublhlUiSKnGpRJIKUzW4E/hhRKyJiD6vDhURiyJidUSsdsufJDVO1eA+OTNPAj4ALI6IOfs2yMyuzOzMzM6Ojo66FilJ+oNKwZ2Zv6393gTcAXgRaUlqkQGDOyIOjYjDdt8G5gGPNbowSVLfquwqORK4o/ZVXO3AtzLz+w2tSpLUrwGDOzM3ACc2oRZJUgVuB5SkwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwlQO7ohoi4ifR8RdjSxIkrR/g5lxXwasb1QhkqRqKgV3RIwHzgRubmw5kqSBVJ1xfx74NLCrgbVIkipoH6hBRJwFbMrMNRHxvv20WwQsApgwYcIbr6ir643/26Fo0aJWV1B/jpHUUlVm3CcDCyJiI/Ad4PSI+Oa+jTKzKzM7M7Ozo6OjzmVKknYbMLgz87OZOT4zJwIfAVZk5scbXpkkqU/u45akwgy4xr23zFwJrGxIJZKkSpxxS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKsyAwR0RIyPioYh4JCJ+GRHXNqMwSVLf2iu0eQU4PTO3RcQI4L6IWJ6ZDzS4NklSHwYM7sxMYFvt7ojaTzayKElS/6rMuImINmAN8E7gi5n5YB9tFgGLACZMmFDPGtUCazdv7/e5Mdt2NLGSxtuyn75KQ1GlNycz87XMnAGMB2ZHxAl9tOnKzM7M7Ozo6Kh3nZKkmkHtKsnM3wMrgfkNqUaSNKAqu0o6IuIttdujgD8HHm90YZKkvlVZ4z4KWFZb5z4IuC0z72psWZKk/lTZVfIoMLMJtUiSKvCTk5JUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUmAGDOyKOjoh7ImJ9RPwyIi5rRmGSpL61V2izE/ibzPxZRBwGrImIH2XmugbXJknqw4Az7sx8NjN/Vrv9IrAeeEejC5Mk9W1Qa9wRMRGYCTzYiGIkSQOrslQCQESMBv4NuDwzX+jj+UXAIoAJEybUrcDidXW1uoI3ZMy2Ha0uQVI/Ks24I2IEvaF9S2b+e19tMrMrMzszs7Ojo6OeNUqS9lJlV0kAXwXWZ+Y/Nb4kSdL+VJlxnwxcAJweEWtrPx9scF2SpH4MuMadmfcB0YRaJEkV+MlJSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhRkwuCPiaxGxKSIea0ZBkqT9qzLjXgrMb3AdkqSKBgzuzLwXeL4JtUiSKqjbGndELIqI1RGxuqenp16HlSTto27BnZldmdmZmZ0dHR31OqwkaR/uKpGkwhjcklSYKtsBvw2sAo6NiO6IuKjxZUmS+tM+UIPMXNiMQiRJ1bhUIkmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMEMvuL/3vaF5rMEcr97nraeKtb384MOV2j068fi6He/wH/xHpWPV2zvPmdeS8w5lDz73v00/56+3vjKkj1fFLf/1+6acZ+gF9113Dc1jDeZ49T5vPVWs7U9vW1ap3S8mnVC34/3xD5dXOla9jf7Pn7bkvEPZ/7y6q+nn/M0LO4b08ap4+qWdTTnP0AtuSdJ+GdySVBiDW5IKY3BLUmEMbkkqTKXgjoj5EfGriHgiIj7T6KIkSf0bMLgjog34IvAB4F3Awoh4V6MLkyT1rcqMezbwRGZuyMxXge8A5zS2LElSf6oE9zuAp/e63117TJLUApGZ+28QcT5wRmb+de3+BcDszPzEPu0WAYtqd48FflX/chtmLLC51UW8SfZh6DgQ+mEfmu+YzOyo0rC9Qptu4Oi97o8Hfrtvo8zsAroqlTfERMTqzOxsdR1vhn0YOg6EftiHoa3KUsnDwJSImBQRBwMfAb7b2LIkSf0ZcMadmTsj4lLgB0Ab8LXM/GXDK5Mk9anKUgmZeTdwd4NraaUil3j2YR+GjgOhH/ZhCBvwzUlJ0tDiR94lqTDDMrgjoi0ifh4Rd9XuL42I/46ItbWfGa2ucSARsTEiflGrd3XtsSMi4kcR8eva77e2us796acP10TEM3uNxQdbXef+RMRbIuL2iHg8ItZHxHsKHIe++lDMOETEsXvVuTYiXoiIy0sbh8EYlkslEXEF0AkcnplnRcRS4K7MvL21lVUXERuBzszcvNdj1wPPZ+bf164p89bM/NtW1TiQfvpwDbAtM/+hVXUNRkQsA36amTfXdl0dAlxFWePQVx8up6Bx2K12iY5ngD8DFlPQOAzGsJtxR8R44Ezg5lbX0gDnALu/I2wZ8BctrOWAFxGHA3OArwJk5quZ+XsKGof99KFUc4HfZOaTFDQOgzXsghv4PPBpYN8v1bsuIh6NiH+OiD9qQV2DlcAPI2JN7VOrAEdm5rMAtd/jWlZdNX31AeDS2lh8bYi/vJ0M9ABfry293RwRh1LWOPTXByhnHPb2EeDbtdsljcOgDKvgjoizgE2ZuWafpz4LTAVmAUcAJbycOjkzT6L3qo2LI2JOqwt6A/rqw5eBPwFmAM8C/9jC+gbSDpwEfDkzZwIvAaVd9ri/PpQ0DgDUlnkWAP/a6loabVgFN3AysKC2tvod4PSI+GZmPpu9XgG+Tu8VEYe0zPxt7fcm4A56a34uIo4CqP3e1LoKB9ZXHzLzucx8LTN3ATcxtMeiG+jOzAdr92+nNwRLGoc++1DYOOz2AeBnmflc7X5J4zAowyq4M/OzmTk+MyfS+5JqRWZ+fK/BDXrXwR5rYZkDiohDI+Kw3beBefTW/F3gwlqzC4E7W1PhwPrrw+6xqDmXITwWmfk74OmIOLb20FxgHQWNQ399KGkc9rKQPyyTQEHjMFjDclcJQES8D7iytqtkBdABBLAWuCQzt7Wyvv2JiMn0zlCh96XutzLzuogYA9wGTACeAs7PzOdbVOZ+7acP/0Lvy/MENgIX716nHIpqW0dvBg4GNgB/Re+EqIhxgH77cANljcMh9F5+enJmbq09Vsx/D4M1bINbkko1rJZKJOlAYHBLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklSY/wNpVFf+DPbmEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Always visually inspect your data\n",
    "import matplotlib.pylab as plt\n",
    "import seaborn as sns # A library for statical plotting\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "sns.distplot(drug, color=\"skyblue\", label=\"Drug\", kde=False, hist=True, rug=\"True\");\n",
    "sns.distplot(placebo, color=\"red\", label=\"Placebo\", kde=False, hist=True, rug=\"True\");\n",
    "plt.legend(); # Seaborn plots are interoperatable with matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What would it look if there were no difference in the groups?\n",
    "#### That is the same if take the data and shuffled labels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assume there is no difference in means of the treatment and control group. If so, we can merge the observations from both groups and we can check how likely is that the observed_diff value or greater appears if there is no difference in the groups.\n",
    "\n",
    "Let's check how likely is for our `observed_diff` value or greater to exist when doing random groups assignment from the common pool many times. If the null hypothesis is true (the drug does not work, on average group means are close to 0) we should observe many occurences of mean differences equal or greater to `observed_diff`. If not, then we can safely reject the null hypothesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([54, 73, 53, 70, 73, 68, 52, 65, 65],\n",
       " [54, 51, 58, 44, 55, 52, 42, 47, 58, 46])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "drug, placebo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[54, 73, 53, 70, 73, 68, 52, 65, 65, 54, 51, 58, 44, 55, 52, 42, 47, 58, 46]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Unite the data\n",
    "combined = drug + placebo\n",
    "combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m \u001b[0mshuffle\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m\n",
       "Shuffle list x in place, and return None.\n",
       "\n",
       "Optional argument random is a 0-argument function returning a\n",
       "random float in [0.0, 1.0); if it is the default None, the\n",
       "standard random.random will be used.\n",
       "\u001b[0;31mFile:\u001b[0m      ~/anaconda3/lib/python3.7/random.py\n",
       "\u001b[0;31mType:\u001b[0m      method\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "shuffle?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[47, 65, 52, 54, 52, 55, 73, 65, 53, 68, 44, 70, 54, 42, 73, 51, 58, 46, 58]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shuffle(combined) # Rearrange in-place\n",
    "combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[47, 65, 52, 54, 52, 55, 73, 65, 53]\n",
      "[68, 44, 70, 54, 42, 73, 51, 58, 46, 58]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEdtJREFUeJzt3X2Q1dV9x/H3110sKNpEWIwJItBSMcqTs9BkNFilQRIU60QzIQljO3bQGUgUY9LEjDP6h20nk6aJYx6KJoE0JtHaOiZW8oiMxcEHSNAYMBNLUdcYWXCC4oiCfPvHXgiQXfa3cu/ePcv7NbOz9+Hwu9+zZ/fDueee+7uRmUiSynFUswuQJPWNwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqTGsjDjpy5MgcO3ZsIw4tSYPSunXrtmZmW5W2DQnusWPHsnbt2kYcWpIGpYh4umpbl0okqTAGtyQVxuCWpMI0ZI1b0pFp165ddHR0sHPnzmaXMmANHTqU0aNHM2TIkDd9DINbUt10dHRw3HHHMXbsWCKi2eUMOJnJtm3b6OjoYNy4cW/6OJWWSiLiLRFxV0Q8GREbI+Ldb/oRJQ1aO3fuZMSIEYZ2DyKCESNGHPYzkqoz7i8BP8zMSyLiaOCYw3pUSYOWoX1o9fj59BrcEXE8MBP4W4DMfB14/bAfWZL0plSZcY8HOoFvRsQUYB1wVWa+0tDKJBVv/db6vkg5deTQQ97f0tLCpEmT2LVrF62trVx22WVcffXVHHXU4NpAVyW4W4EzgY9l5sMR8SXg08D1+zeKiIXAQoAxY8bUu86i1PuXVb3/wUoAw4YNY/369QBs2bKFD3/4w2zfvp0bb7zxgHa7d++mtbXcvRlV/hvqADoy8+Ha9bvoCvIDZObSzGzPzPa2tkpvt5ekhhk1ahRLly7llltuITNZtmwZl156KRdeeCGzZ89m1apVXHDBBfvaL168mGXLlgFw3333MXHiRM4++2w+/vGPH9BuIOg1uDPzd8CzEXFq7aZZwIaGViVJdTB+/Hj27NnDli1bAFizZg3Lly9n5cqVPf6bnTt3csUVV7BixQpWr15NZ2dnf5VbWdWFn48Bt0fE48BU4B8bV5Ik1U9m7rv83ve+lxNOOOGQ7Z988knGjx+/b5/1/PnzG1rfm1FpkScz1wPtDa5Fkupq06ZNtLS0MGrUKACOPfbYffe1trayZ8+efdf37q3eP+gHqsH1Uqsk1XR2dnLllVeyePHibvdOn3LKKWzYsIHXXnuN7du387Of/QyAiRMnsmnTJjZv3gzAHXfc0Z9lV1Luy6qSBrz+3g306quvMnXq1H3bARcsWMA111zTbduTTz6ZD37wg0yePJkJEyYwbdo0oGtnyle+8hXmzJnDyJEjmTFjRn92oZJoxNOC9vb2PJI/SMHtgPXndsAybNy4kdNOO63ZZRy2HTt2MHz4cDKTRYsWMWHCBJYsWVK343f3c4qIdZlZaUnapRJJOsitt97K1KlTOf3009m+fTtXXHFFs0s6gEslknSQJUuW1HWGXW/OuCWpMAa3JBXG4JakwhjcklQYX5yU1DhLl9b3eAsXHvLuvad13b17N6eddhrLly/nmGOOYfjw4ezYsaNuZdxwww0MHz6ca6+9tm7H7Atn3JIGjb2ndX3iiSc4+uij+drXvtbskhrC4JY0KL3nPe/hqaeeOuC2HTt2MGvWLM4880wmTZrEPffcs+++b33rW0yePJkpU6awYMECoOtt8x/4wAeYPn0606dP58EHH9zX/rHHHuO8885jwoQJ3HrrrUDXeU4++clPcsYZZzBp0qSGvV3epRJJg87u3btZsWIFc+bMOeD2oUOHcvfdd3P88cezdetW3vWudzFv3jw2bNjATTfdxIMPPsjIkSN58cUXAbjqqqtYsmQJZ599Ns888wznn38+GzduBODxxx/noYce4pVXXmHatGnMnTuXNWvWsH79eh577DG2bt3K9OnTmTlzJieddFJd+2dwSxo09p6rBLpm3JdffvkB92cm1113HQ888ABHHXUUzz33HC+88AIrV67kkksuYeTIkQD7Tv3605/+lA0b/vDxAy+99BIvv/wyABdddBHDhg1j2LBhnHvuuTzyyCOsXr2a+fPn09LSwoknnsg555zDo48+yrx58+raT4Nb0qCx/0eXdef222+ns7OTdevWMWTIEMaOHcvOnTvJzG7PILhnzx7WrFnDsGHD/ui+g9tHRL+dEtY1bklHjO3btzNq1CiGDBnC/fffz9NPPw3ArFmzuPPOO9m2bRvAvqWS2bNnc8stt+z79/v/p3DPPfewc+dOtm3bxqpVq/Yti9xxxx288cYbdHZ28sADDzTk7ILOuCU1Ti/b9/rbRz7yES688ELa29uZOnUqEydOBOD000/ns5/9LOeccw4tLS1MmzaNZcuWcfPNN7No0SImT57M7t27mTlz5r6dKjNmzGDu3Lk888wzXH/99bz97W/n4osvZs2aNUyZMoWI4HOf+xxve9vb6t4PT+vaAJ7Wtf48rWsZBstpXRvN07pK0hHG4JakwhjckuqqhA/bbaZ6/HwMbkl1M3ToULZt22Z49yAz2bZtG0OHHt5rNu4qkVQ3o0ePpqOjg87OzmaXMmANHTqU0aNHH9YxDG5JdTNkyBDGjRvX7DIGPZdKJKkwBrckFabSUklEbAZeBt4AdlfdJC5Jqr++rHGfm5lbG1aJJKkSl0okqTBVZ9wJ/DgiEvi3zPyjD5KLiIXAQoAxY8a8+Yrq/Rl1TTBix65mlzD4DB/S7Aqk3vXTSbWqzrjPyswzgfcBiyJi5sENMnNpZrZnZntbW1tdi5Qk/UGl4M7M39a+bwHuBup/gllJUiW9BndEHBsRx+29DMwGnmh0YZKk7lVZ4z4RuLv2MT2twHcy84cNrUqS1KNegzszNwFT+qEWSVIFbgeUpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMJUDu6IaImIX0TEvY0sSJJ0aH2ZcV8FbGxUIZKkaioFd0SMBuYCtzW2HElSb6rOuL8IfArY08BaJEkV9BrcEXEBsCUz1/XSbmFErI2ItZ2dnXUrUJJ0oCoz7rOAeRGxGfgecF5EfPvgRpm5NDPbM7O9ra2tzmVKkvbqNbgz8zOZOTozxwIfAlZm5kcbXpkkqVvu45akwrT2pXFmrgJWNaQSSVIlzrglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVptfgjoihEfFIRDwWEb+KiBv7ozBJUvdaK7R5DTgvM3dExBBgdUSsyMyHGlybJKkbvQZ3Ziawo3Z1SO0rG1mUJKlnVWbcREQLsA74c+DLmflwN20WAgsBxowZU88aJTXAszt2NbuEQefkfnqcSi9OZuYbmTkVGA3MiIgzummzNDPbM7O9ra2t3nVKkmr6tKskM38PrALmNKQaSVKvquwqaYuIt9QuDwP+Gniy0YVJkrpXZY37JGB5bZ37KODOzLy3sWVJknpSZVfJ48C0fqhFklSB75yUpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpML0GtwRcXJE3B8RGyPiVxFxVX8UJknqXmuFNruBT2TmzyPiOGBdRPwkMzc0uDZJUjd6nXFn5vOZ+fPa5ZeBjcA7Gl2YJKl7fVrjjoixwDTg4UYUI0nqXeXgjojhwH8CV2fmS93cvzAi1kbE2s7OznrWKEnaT6XgjoghdIX27Zn5X921ycylmdmeme1tbW31rFGStJ8qu0oC+DqwMTO/0PiSJEmHUmXGfRawADgvItbXvt7f4LokST3odTtgZq4Goh9qkSRV4DsnJakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKozBLUmFMbglqTAGtyQVxuCWpMIY3JJUGINbkgpjcEtSYQxuSSqMwS1JhTG4JakwBrckFabX4I6Ib0TEloh4oj8KkiQdWpUZ9zJgToPrkCRV1GtwZ+YDwIv9UIskqYK6rXFHxMKIWBsRazs7O+t1WEnSQeoW3Jm5NDPbM7O9ra2tXoeVJB3EXSWSVBiDW5IKU2U74HeBNcCpEdEREZc3vixJUk9ae2uQmfP7oxBJUjUulUhSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVBiDW5IKY3BLUmEMbkkqjMEtSYUxuCWpMAa3JBXG4JakwhjcklQYg1uSCmNwS1JhDG5JKszAC+4f/GBgHqsPx3t87On1fdw6Ov5H/12p3U+mnlupXdW+vvrwo3U7Vt19/vPNedwBrOr411PV381mHa+Kpz9xXb88zsAL7nvvHZjH6sPxfjnujPo+bh396Y9XVGq35a2jKrWr2te/uHN53Y5Vd7/5TXMedwCrOv71VPV3s1nHq+KUL/xTvzzOwAtuSdIhGdySVBiDW5IKY3BLUmEMbkkqTKXgjog5EfHriHgqIj7d6KIkST3rNbgjogX4MvA+4J3A/Ih4Z6MLkyR1r8qMewbwVGZuyszXge8BFzW2LElST6oE9zuAZ/e73lG7TZLUBJGZh24QcSlwfmb+fe36AmBGZn7soHYLgYW1q6cCv65/uQ0zEtja7CIOk30YOAZDP+xD/zslM9uqNGyt0KYDOHm/66OB3x7cKDOXAksrlTfARMTazGxvdh2Hwz4MHIOhH/ZhYKuyVPIoMCEixkXE0cCHgO83tixJUk96nXFn5u6IWAz8CGgBvpGZv2p4ZZKkblVZKiEz7wPua3AtzVTkEs9B7MPAMRj6YR8GsF5fnJQkDSy+5V2SCnNEBndEtETELyLi3tr1ZRHxfxGxvvY1tdk19iYiNkfEL2v1rq3ddkJE/CQiflP7/tZm13koPfThhoh4br+xeH+z6zyUiHhLRNwVEU9GxMaIeHeB49BdH4oZh4g4db8610fESxFxdWnj0BdH5FJJRFwDtAPHZ+YFEbEMuDcz72puZdVFxGagPTO37nfb54AXM/Ofa+eUeWtm/kOzauxND324AdiRmUV8nlhELAf+JzNvq+26Oga4jrLGobs+XE1B47BX7RQdzwF/CSyioHHoiyNuxh0Ro4G5wG3NrqUBLgL2fkbYcuBvmljLoBcRxwMzga8DZObrmfl7ChqHQ/ShVLOA/83MpyloHPrqiAtu4IvAp4A9B91+U0Q8HhH/GhF/0oS6+iqBH0fEutq7VgFOzMznAWrf+/+DA/umuz4ALK6NxTcG+NPb8UAn8M3a0tttEXEsZY1DT32AcsZhfx8Cvlu7XNI49MkRFdwRcQGwJTPXHXTXZ4CJwHTgBKCEp1NnZeaZdJ21cVFEzGx2QW9Cd334KvBnwFTgeeBfmlhfb1qBM4GvZuY04BWgtNMe99SHksYBgNoyzzzgP5pdS6MdUcENnAXMq62tfg84LyK+nZnPZ5fXgG/SdUbEAS0zf1v7vgW4m66aX4iIkwBq37c0r8LeddeHzHwhM9/IzD3ArQzssegAOjLz4dr1u+gKwZLGods+FDYOe70P+HlmvlC7XtI49MkRFdyZ+ZnMHJ2ZY+l6SrUyMz+63+AGXetgTzSxzF5FxLERcdzey8Bsumr+PnBZrdllwD3NqbB3PfVh71jUXMwAHovM/B3wbEScWrtpFrCBgsahpz6UNA77mc8flkmgoHHoqyNyVwlARPwVcG1tV8lKoA0IYD1wZWbuaGZ9hxIR4+maoULXU93vZOZNETECuBMYAzwDXJqZLzapzEM6RB/+na6n5wlsBq7Yu045ENW2jt4GHA1sAv6OrglREeMAPfbhZsoah2PoOv30+MzcXrutmL+Hvjpig1uSSnVELZVI0mBgcEtSYQxuSSqMwS1JhTG4JakwBrckFcbglqTCGNySVJj/BwH9SnxAFLgtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Randomly assign to drug and placebo group\n",
    "drug_random = combined[:len(drug)]\n",
    "placebo_random = combined[len(drug):]\n",
    "\n",
    "print(drug_random)\n",
    "print(placebo_random)\n",
    "\n",
    "sns.distplot(drug_random, color=\"skyblue\", label=\"Drug\", kde=False, hist=True, rug=\"True\");\n",
    "sns.distplot(placebo_random, color=\"red\", label=\"Placebo\", kde=False, hist=True, rug=\"True\");\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate it a bunch of times\n",
    "n = 10_000\n",
    "count = 0\n",
    "simulated_means = []\n",
    "\n",
    "for _ in range(n):\n",
    "    shuffle(combined)\n",
    "    shuffled_diff = mean(combined[:len(drug)]) - mean(combined[len(drug):])\n",
    "    simulated_means.append(shuffled_diff)\n",
    "    count += (abs(shuffled_diff) >= observed_diff) # abs() as we test two-tailed alternative hypothesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10,000 label reshufflings produced only 16 instances \n",
      "with a difference at least as extreme as the observed difference of 12.97.\n"
     ]
    }
   ],
   "source": [
    "print(f\"\"\"{n:,} label reshufflings produced only {count} instances \n",
    "with a difference at least as extreme as the observed difference of {observed_diff:.2f}.\"\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0016"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The p-value is a chance of observing the current difference when there is truly no difference.\n",
    "# Probability of obtaining a result equal to or more extreme than was observed in the data.\n",
    "p = count / n\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha=0.005, p-value=0.0016\n"
     ]
    }
   ],
   "source": [
    "print(\"alpha={}, p-value={}\".format(alpha/2,p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to notice, that this p-value value has nothing to do with the p-value calculated with the analytical method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reject the null hypothesis that the means are equal.\n"
     ]
    }
   ],
   "source": [
    "# Interpret via p-value\n",
    "if p < alpha:\n",
    "    print('Reject the null hypothesis that the means are equal.')\n",
    "else:\n",
    "    print('Accept null hypothesis that the means are equal.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center><img src=\"images/p_values.png\" width=\"30%\"/></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7wAAAHzCAYAAAAD0hGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmYZHdV//HPqb2ql5me7p59nxDIhEz2BQQMq6AILqC4gqKCP1FRcENFQFDxUUB/4E9FBAQlKIoCBgGNRCFkJZM9k8yWTM/a+1p7nd8f93ZS6XTPdPd0LV31fj1PP11ddavuqVv33r6nzncxdxcAAAAAAK0m0ugAAAAAAACoBRJeAAAAAEBLIuEFAAAAALQkEl4AAAAAQEsi4QUAAAAAtCQSXgAAAABASyLhBbAoZvZjZvbVGr32J8zsvbV47XnW9QYz+8Z5vsYDZnZ9eNvM7ONmNmpmt4f3/byZnTazKTPrXYGwV42lvnczO2pmL6lHbCvNzNzMLljG83aGz43VIq5wHee9n6M2zOz7zexYeIxcXqN1LGvfrEEcTxzfZvYuM/v0Astdb2YD9Y1ufrU4J5nZO8zsb1byNQEsHgkvgCeY2fPM7BYzGzezETP7ppldLUnu/vfu/rImiPHrZvYzNXrt2URkKvw5bWZfMrOXVi/n7he7+9fDP58n6aWStrr7NWYWl/QBSS9z9053H65FrM2ond87sAR/Iukt4TFy9/m+WC3Pia2uXl+2uvsfuDufEdAgJLwAJElm1i3pS5L+r6R1krZIerekfCPjapC17t4p6VJJX5P0eTN7wwLL7pB01N2nw783SEpJemA5Kzaz6HKe1yTO673PVcsKKGqDz2xRdqg9zw8A0BAkvABmXShJ7v4Zdy+7e9bdv+ru90pPbyIZVkL/j5k9amaTZvb7ZrbHzL5lZhNm9o9mlpjvuVXPf1qTOzPrCauqg2Ez4S+Z2dbwsfdJer6kD4cV2A+H9z/LzL4WVqUPmNkPVb1er5l9IYzpdkl7FrtB3P2Uu/+ZpHdJer+ZRcLXPGpmLzGzN0r6G0nPCeP5jKQD4dPHzOymRcT3CTP7f2Z2o5lNS3qhmSXN7E/M7PGwyvyXZpYOl7/ezAbM7G1mdsbMTprZT1W9XtrM/tTMHgsr9d+oeu51YQV/zMzusbBZdtVndDj8LI+Y2Y/Nt03C2D5kZifCnw+F910433uf5/k/EcY2bGa/Peexd5nZ58zs02Y2IekNcyswNqfpo5ldYWZ3h3H/k5l9dqGKTbh/3hSue8jM/t7M1lY9ftTM3m5m94bb7rNmlqp6/NfC7X3CzH56vnVULft1M/tDM7s9fK1/M7N1Cyz7U2b2UPgeDpvZm+Y8/moz2x/uw4fM7OXh/WvM7GNhTMfN7L321ITIzOz/hut/2MxeXPXA5vC4GDGzg2b2s2d5L71m9sVw/XeE65l7LvgFM3tU0qPhfc8Nlx0Pfz93znZ+SdXfTzR1tSdbWfxcuJ1PmtnbzhLbJ8zsL8zsyxYcg980s43hfjkavu/Lq5bfbGb/bMH55YiZ/VLVY9dYcP4aC9f7YQvPYVXv880WnPNGzewjZmbhYxeY2c3h+x0ys8/OE2vSzKYkRSXdY2aHwvsvCveXMQu6S7xqzvt7yvlhzmvOe04MvWS+WMPn/XS4z42a2VfMbMcC23f283i9BeejIas6bu0cx+dS2TznNTO72oLzYKxquR80s/3h7dnzxmctOIa+bWaXVi077/Y1s5+T9GOSfj3cdl+sCuUyW/g88EoLjscxC86n+6oe+w0LjsVJC871L66KcXYfT1lwjhsOX+MOM9uw3G0GYBHcnR9++OFHkrolDUv6pKRXSOqZ8/gbJH2j6m+X9IXweRcrqAT/l6TdktZIelDS6+d7btXzLwhvf0LSe8PbvZJ+UFJGUpekf5L0r1XP+7qkn6n6u0PSMUk/JSkm6QpJQ5IuDh+/QdI/hss9W9LxubFUvdbOMK7YnPt3h/dfFP59VNJLFtguT3mNRcT3CUnjkr5DwZeQKUkfCrftunAbfFHSH4bLXy+pJOk9kuKSvlvSzOznJekj4TbaouDC+rmSkuHfw+HyEQXNsIcl9YcxTkh6Zvgam2bjm2cbvUfSrZLWh8+9RdLvn237VT13r6QpSS8IY/pA+F5mt+W7JBUlfV8YY1pV+0bV+x8IbyckPSbpl8Nt8QOSCtXLz1n/BeH7Toax/4+kD1U9flTS7ZI2h9v+IUlvDh97uaTTCvahDkn/oKp9eJ51fV3Bvja7/D9L+vQC+8j3KPgixiR9Z/h5XhE+dk24f7w03CZbJD0rfOxfJf1V+Prrw9jfVLVfliT9Srhtfjh8nXXh4zdL+gsF+9tlkgYlvXiB93JD+JMJP8Njevq54GvhNkuHv0cl/YSCff5Hwr975x4/VZ/73G3zmfB9XRLG9pIFYvuEguPpyvC93CTpiKSfVLD/v1fSf4fLRiTdJemdCvad3ZIOS/qu8PErJV0Xxrwz/PzfOud9fknSWknbw7heHj72GUm/rSeP4eed5Vxbfe6LSzoo6R1hTC+SNKknj8VPaM75YYF97WfmWcdCsX5fuM6Lwvf6O5JuOcc58aPhZ3upgnP9RVXxzXt8znOufOJznmc91+vs57UHJb2iavnPS3rbnPPGa8Lnvj3cB+KL3L7vnRPLUS18HrhC0hlJ1yrYv14fLp+U9EwFx8bmqm23Z559/E0KzumZ8DWulNS90P7CDz/8nP8PFV4AkiR3n1DQH3X24mbQggrQ2b55fr+7T7j7A5Lul/RVdz/s7uOSvixpyQOyuPuwu/+zu8+4+6Sk9ylIAhbySgVNij/u7iV3/7aC5OI1FlS7flDSO9192t3vV5DQL9WJ8Pe8FbpzWDC+qmX+zd2/6e4VBReTPyvpV9x9JNwGfyDpdVXLFyW9x92L7n6jgiTymRZUoH9a0i+7+3EPKvW3uHte0o9LutHdb3T3irt/TdKdCi4sJaki6dlmlnb3k+FnOp8fC9d9xt0HFTR7/4lFbovXSPqSu/9PGNPvhuut9i13/9cwxuw5Xm82OfnzcFv8i4IL1Xm5+0F3/5q758PYP6Cn71t/7u4n3H1EwUXpZeH9PyTp4+5+vwfN19+1iPf7qarlf1fSD9k8TVLd/d/d/ZAHbpb0VQVVO0l6o6S/DeOuhJ/rw+Fx+QoFCdm0u5+R9EE9dT85oyChL7r7ZxVU4L/HzLYpONZ/w91z7r5fQUuFp32OVcfQ74XH5IOa/xj6w3B/zSpI4B9190+F+/xnJD0s6XsXsc1mvTt8X/dJ+riCpHkhn3f3u9w9pyARyrn737l7WdJn9eR56GpJ/e7+HncvuPthBee610lS+Bq3hjEfVfBlwtz944/cfczdH5f033py/ygqaKq8Odymix0w7DpJneHrFtz9JgWJavX7feL8EL7HxVoo1jcp+LwecveSgvPLZQtVeUPv9qDVzz2S7lGQ+K60ec9r4WOfVHAOkwUtJb5LwZdOs+5y98+5e1HBcZ1SsG0Xs33ns9B54Gcl/ZW73xaeXz+p4Jx9naSygsR3r5nF3f2oux9a4H32KvjSoxzudxOL3EYAloGEF8ATwgugN7j7VgWVqc0Kqo0LOV11OzvP351LjcHMMmb2VxY0e51QUIVbO1+iENoh6dqwadiYmY0pSMo2KqjixRR86z7rsaXGpKCqJkkjy3ju2eKbVR1fv4Jv/u+qWv4/wvtnDYcXqrNmFGzrPgUXevNdZO2Q9No5cTxP0qYwIfthSW+WdNLM/t3MnrXA+9msp27Dx8L7FmNz9XsN1zt3YKtjWrzNko67uy/m+Wa23sxuCJscTkj6tIJtVu1U1e3Z7fq02LW4/Wju8vF51icze4WZ3WpB8+IxBV9CzC63TQt/nnEFn9fs5/lXCiq9s+Zum9nParOk2S9Tqh/boqeb7xiabxtX3zd3Hznb6y9k7rY72z622PPQDkmb5xwD71DQ91xmdqEFXShOhfvHH2jx+8evK6jQ3x42mz1rk/cqmyUdC7/smjV3Wy3lmFhMrDsk/VnVNhhREPvZPp+FXmslLXRek4Jj9XvNrFPBl0//6+4nq5atPq9UJA3oyX39XNt3Pmfbdm+bsw9tU/BFx0FJb1XwZdiZ8Fwz3377KUlfkXSDBc32/9iCAf8A1AgJL4B5ufvDCpp7PXsFXm5aQRInSTKzjWdZ9m0KvtW/1t27FTR/lYILMimoQFc7Julmd19b9dPp7j+voBlfScEFyazty4j/+xVUyw6ca8F5nC2+WdXvaUjBRfrFVcuv8WAQrXMZkpTT/P2UjymoOFbH0eHufyRJ7v4Vd3+pgubMDyuofM3nhIKLvlnb9WQF/FxOquqzMLOMgkpHtbmf71P2HT31i4KTkrZU903UUz/ruf4wfP194b7143pyvzqXp8Suxe1Hc5cvKviMnmBmSQUV/z+RtMHd10q6sSquY1r488xL6qv6PLvd/eKqZeZum9nP6oSkdWbWNeex4/OsZ/YY2rrA+5pV/bnN3Ufmvv7ZPtP51rGUfexsjkk6MucY6HL32VYO/0/Bvv+McP94hxa5f3jQ3/9n3X2zggrqX9jipgU6IWlb2Dpj1tzPYu4x8bTVLybGKscUNH2v3g5pd79lia8jLe6zPG/uflzStxSci39CQdJYrfq8ElGwv87u62fbvsvZdu+bs+0yYSsGufs/uPvzFOz/Lun987yXoru/2933Kuhy8koFTfAB1AgJLwBJTwys9DZ7coCobQqafd26Ai9/j6SLzeyycPCPd51l2S4FCd9Y2HTt9+Y8flpB37tZX5J0oQWDIcXDn6vN7KKwSeO/SHpXWDneq6DP1aKY2QYze0sYw2/NqRIs1oLxzbdwuI6PSvqgma0P49hiZt91rhWFz/1bSR+wYHCeqJk9J0yqZisk3xXen7JggJmt4ft8lZl1KEiiphQ0z5vPZyT9jpn1m1mfgv6Q886tOY/PSXqlBdNfJRT01zvX/6H9kr7bzNaFX5S8teqxb4VxvsXMYmb2agV9XhfSpeC9jZnZFkm/tsi4paAf+BvMbG+YqM/dL+fz41XLv0fS58J9slpCQTPIQUklM3uFpOrpvz4m6afM7MVmFgn3hWeF1a2vSvpTM+sOH9tjZtVNcNdL+qVwn3utgj6bN7r7MQV9r/8w3A/2KWg6/fdz38A8x9CzdO6L8xsV7PM/Gn4uP6yg7++Xwsf3S3pdGNdVemrz/lm/G67vYgX93582CNQy3C5pwoKBhdLhcfBsC6deU7B/TEiaCt/nzy/4SnOY2Wtnz50K+iu7Fj6Gqt2mIGn89XB7XK+g6fcNi123nn5OPJe/lPRb4badHfzstUt4frWzHZ8r7e8UVNIvUdB0vdqVZvYDFgxs9VYF57Fbde7tu9Rt91FJbzazay3QYWbfY2ZdZvZMM3tReL7NKfg/9rR9wMxeaGaXWNBqaULBF2GL2VcALBMJL4BZkwoG4rjNgtFAb1XQL3fBEVIXy90fUXDB/58KRnE9W/+2DykYHGUojOE/5jz+Zwr6546a2Z+HzTJfpqAf3gkFTdHeryCJkKS3KGiOdkpBxfrjiwh5LNwG9yloXvpad//bRTzvaRYR33x+Q8FAK7da0LTyP/VkX7ZzeXsY9x0Kmiq+X1IkTHJeraBqNaigUvFrCv4PRBR8zifC53ynpP+zwOu/V0Hf33vD9Xw7vO+cPOgX/AsK+t6dVJAYnGtE108p+MLkqIIE74nEx90LCgaqeqOkMQUV2y9p4am03q1g0JlxSf+uIJFbFHf/soJ98yYFn828o1DPE/snFHzmKUm/NHeBcP/4JQUJ9aikH1UwYNns47crSPg+GMZ9s56snv6kgoT5wfC5n1NQoZ91m6RnKDiW3ifpNf7k3Mg/omBQnRMKkoff86Bf93zeomAgulPhe/qMzjJdWbiOVyrYp4YVJCmvdPfZ6vbvKqhajyr4TP5hnpe5WcF2/i9Jf+LuX11ofYsVJu/fq6A/5hEF2+VvwvcmBcfOjyo4F35US0uyr1Zw7pxS8Pn9srsfWURMBUmvUtAfe0jBQGI/GbawWaynnBMXsc7PKzgv3BCeX+4P178cCx6fNfB5Bfv+5/3JaeBm/ZuCbhmzg6X9QFhJPdf2/ZiCPrdjZvav5wrA3e9U0I/3w+G6DioYIE4Kzul/FK7nlIIvnN4xz8tsVHCsTigYEOtmLf5LQwDLYE/t3gMAwOplZrdJ+kt3X8wXG7WM4+sKRmX9m0bGUQtm9n5JG9190a0llvDaOxWOsDunPycgC6ZyepO7/2fVfe9SMADUjzcsMABNjQovAGDVMrPvtGDe1ZiZvV7SPj29VQDOQ9jdYV/YhPMaBRX1uU1KgZoysx9U0FR8Ma0rAOAJsXMvAgBA03qmgubAnQpGM37NnNFbcf66FDRj3qxg8LY/VdCEFKiLsMXEXkk/scyxFAC0MZo0AwAAAABaEk2aAQAAAAAtiYQXAAAAANCSWrIPb19fn+/cubPRYQAAAAA4HwcOBL+fudjZ+dAu7rrrriF37z/Xci2Z8O7cuVN33nlno8MAAAAAcD6uvz74/fWvNzIKNCEze2wxy9GkGQAAAADQkkh4AQAAAAAtiYQXAAAAANCSSHgBAAAAAC2JhBcAAAAA0JJIeAEAAAAALYmEFwAAAADQkkh4AQAAAAAtiYQXAAAAANCSSHgBAAAAAC2JhBcAAAAA0JJIeAEAAAAALYmEFwAAAADQkkh4AQAAAAAtiYQXAAAAANCSSHgBAAAAAC2JhBcAAAAA0JJIeAEAAAAALYmEFwAAAADQkmKNDgAAAJy//UO5uq7vsr5UXdcHAMByUOEFAAAAALQkEl4AAAAAQEsi4QUAAAAAtCQSXgAAAABASyLhBQAAAAC0JBJeAAAAAEBLIuEFAAAAALQk5uEFAABLttx5f5m/FwBQT1R4AQAAAAAtiYQXAAAAANCSaNIMAMAClttsV6LpLgAAzYAKLwAAAACgJZHwAgAAAABaEgkvAAAAAKAlkfACAAAAAFoSCS8AAAAAoCWR8AIAAAAAWhIJLwAAAACgJZHwAgAAAABaEgkvAAAAAKAlxRodAAAArWj/UG5Zz7usL7XCkQAA0L6o8AIAAAAAWhIJLwAAAACgJZHwAgDQQkoVV6nijQ4DAICmQB9eAABWoVLFNZwrazBX0lC2rEMTBU0WK8qVg2Q3ZlI6FlEqakpFTelYRB2xiPrSUUXNGhw9AAD1QcILAMAq4e46Nl3S3YNZPTJeUJjbKmJSJhZRTzKqznhEJilXdmVLQQI8ViirWAmWjUekLZm4tnbG1RmnoRcAoLWR8AIA0ORypYruH81r/1BOQ7myklHTZX0pbeuIqy8dVU8yqvuG82d9jVIlSHwHpkp6bKqoo1NF9SQj2tYR14ZMjKovAKAlkfACANCkhrIl3T6Y1UOjeRUr0qZMTN+9vVMX9SQVjywtQY1FTH2pmPpSMeXLFR2fLmlguqh7R/KKj+W1uyuhHV1xRUh8AQAthIQXAIAmM1Os6H9PzWj/UE6xiLS3J6nL+9LamFmZf9vJaES7uxPa1RXXSL6so5NFHRgv6PhMSRf3JNWTjK7IegAAaDQSXgAAmkTFXbedntEtp7MqlF2X96X0vE0ZZWK16WtrZupNxdSbiun0TEkPjeV125mstnXEdOHapVeRAQBoNiS8AAA0mLvrTLasA+N5zZRcu7vjetHmDvWl6/dvekMmpt5UVI+OF/TYVFGns2Vd1JPQxnRMRjNnAMAqRcILAEADTRcremA0r5F8WZ2xiH5oT5d2dycaEkssYrqoJ6nNHTE9MJLXPcN5nUiVdMm6lBJRkl4AwOpDwgsAQAO4ux6bKuqR8YIiki5am9C2znjDkt1qaxJRPWdDOohvrKBvnp7RZb0p+vYCAFYdEl4AAOpsqljR/SM5jRUq6k9FdXFPUqka9dNdLjPTzq6E1iWj2j+c0+1nsrpwTUI7u+I0cQYArBokvAAA1EnFXUcnizo4XlA0Iu1bl9SmTHP3ke1ORPWcDRndP5LTgfGCRvNlXdKbYkArAMCqQMILAEAdTBUrum8kp/FCRRvSUe3tSSoZba6q7kLiEdNlvSk9NlXUgbGCbjkVNHFeQxNnAECTI+EFAKCG3F0D08GUP1GTLu1NrsqRj2ebOK9NBE2cbz2T1d6epLZ1xhsdGgAACyLhBQCgRgpl1wOjOZ3OltWbjOqS3qRSq6Squ5C1yaieuyGje0ZyemA0r8liRc9am1BklSXwAID2QMILAEANjOTKunckp3zZ9cwWG+wpETVd2ZfSI2MFHZ0qaqpY0WW9i5u6aP9QblnrvKwvtaznAQDa2+r+mhkAgCZTcdej43ndPphVxKTrNqS1qzvRMsnurIiZntWT1CXrkhrNl/Wt0zOaLJQbHRYAAE9BhRcAgBWSL1e0fzin0XxFWzIxXdSTVKzFRzPe0hFXRyyiu4eCfr37elPakObyAgDQHKjwAgCwAsbzZd1yOqvxQkX71iV1SW+q5ZPdWWuTUT1nQ/qJxPfQREHu3uiwAAAg4QUA4HwNTBd125msIpKuW5/W5o72G7k4FYvo2vVpbcrE9Oh4QfeP5FUh6QUANBhtjgAALW+5AyWdS8VdB8YKemyqqHXJ6KIHbmpV0Yhp37qkOmIRHZwoaKZc0eW96bbeJgCAxqLCCwDAMhTKrjsGs3psqqidnXFd1d/eye4sM9MFaxLaty6psXxFt56Z0XSx0uiwAABtioQXAIAlmipWdMvpmSf66z6rJ8k8tHNs7ojrmvVpFSuuW8/MaCTPCM4AgPqjSTMAYNWoVdPkpRjLl3XXUFYRma5dn9aaRLTRITWtnmRUz9mQ0V2DWd1xJqtnr0tqSxv2bwYANA4VXgAAFmkoV9Idg1nFzHTtBpLdxcjEIrpuQ0Y9yajuG8nr8ESh0SEBANoICS8AAItwaqakuwZzSsciunZDWpkY/0IXKx4xXdWf0sZ0TI+MF3RgLM+0RQCAuqjbf2sze7mZHTCzg2b2m/M8njSzz4aP32ZmO8P742b2STO7z8weMrPfqlfMAABI0rGpovYP57QmEUy9k4qS7C5VxEyX9ia1vTOuI5NF3T/KtEUAgNqry39sM4tK+oikV0jaK+lHzGzvnMXeKGnU3S+Q9EFJ7w/vf62kpLtfIulKSW+aTYYBAKi1wxMFPTCaV18qqqv704pHGJxqucxMF61NaE93XMenS9o/nFOZpBcAUEP1+or6GkkH3f2wuxck3SDp1XOWebWkT4a3PyfpxWZmklxSh5nFJKUlFSRN1CdsAEA7e3Q8r0fGC9qUiemKvpSiJLvnzcz0jDVJXbQ2oTPZsu4azKlUIekFANRGvRLeLZKOVf09EN437zLuXpI0LqlXQfI7LemkpMcl/Ym7j9Q6YABAezsyWdChiaK2dsS0bx3TDq20HV3BXL2j+bJuP5NVoUzSCwBYefVKeOe7Spj7n22hZa6RVJa0WdIuSW8zs91PW4HZz5nZnWZ25+Dg4PnGCwBoYwNTRR0YK2hjOqaLe5Iykt2a2NwR1xV9KU2VKrpjMKsilV4AwAqrV8I7IGlb1d9bJZ1YaJmw+fIaSSOSflTSf7h70d3PSPqmpKvmrsDd/9rdr3L3q/r7+2vwFgAA7eD0TEn3j+bVm4pqXy/Jbq31p2O6vDelqWJFdw1mad4MAFhR9Up475D0DDPbZWYJSa+T9IU5y3xB0uvD26+RdJMHcxY8LulFFuiQdJ2kh+sUNwCgjQznSk+Mxnx5b4pmzHXSn47p0t6UxgoV3T2UY/RmAMCKqUvCG/bJfYukr0h6SNI/uvsDZvYeM3tVuNjHJPWa2UFJvyppduqij0jqlHS/gsT54+5+bz3iBgC0j7F8Wd8eyqkjHtGVfWnFGKCqrjZmYnp2T1LD+bLuGSbpBQCsjFi9VuTuN0q6cc5976y6nVMwBdHc503Ndz8AACtlqljRXUNZJSKmq/pTSkRJdhtha2dcJXc9PFbQ/SN5XbKOJuUAgPNTt4QXAIBmlC8HAyZFZLp6fVqpaL16+2A+O7sSKlWkgxMFxSLBvL0kvQCA5SLhBQC0rYq77h4K5oG9dn1amRjJbjPY0x1XqeI6OlVUPCI9Y02y0SEBAFYpEl4AQNt6aDSvsUJFl/Ym1Z2INjochMxMz1ybUNFdhyaK6o7z2QAAloevsgEAbenxqaKOTZe0qyuuTZl4o8PBHGami3uSWpOI6N6RnEZy5UaHBABYhUh4AQBtZzRf1kOjefWlorpwTaLR4WABETNd1ptSxKR/OTKhQpmRmwEAS0PCCwBoK7lSMNdrOma6tDfFgEhNLh2L6NLelIZzZf3HsSk50xUBAJaAhBcA0DbK7rp7OKeyuy7vSynOXLurQl8qphdsyujB0bzuGsw1OhwAwCrCoFUAgLbg7npwJK/xQkWX96bUxUBIq8p1G9I6PlPSTcentSET07bOxfW73j+0vAT5sr7Usp4HAGguVHgBAG1hYLqk4zMl7emOa0OG73tXGzPTK7d3ak0yon87MqmpYqXRIQEAVgESXgBAy5suVvTwWF69yagu6GaQqtUqFYvo+3d1K1eu6F+PTKhCf14AwDmQ8AIAWlrFXfeO5BQx6ZJ1SQapWuXWp2N6+fZODUyXdPuZbKPDAQA0ORJeAEBLOzRR0Hihoot7UkrF+LfXCp69LqUL1yT0jZMzzM8LADgr/vMDAFrWaL6sQxNFbc7EtJF+uy3lpds6FI2YvnxskqmKAAALIuEFALSkUsV173BO6ahpb0+y0eFghXXFo3rRlg4dmyrpnuF8o8MBADQpEl4AQEt6aCyvbNm1rzelGPPttqR965La0RnXfx+f1mSBps0AgKcj4QUAtJxTMyUdnw6mIOpJMt9uqzIzvXx7p8ru+srANE2bAQBPQ8ILAGgpuXJFD4zmtCYR0R6mIGp5Pcmonr8po4PjBT08Vmh0OACAJkPCCwBoGe6u+0fyqri0b11KEaYgagtXr09rYyamrw1MKVuqNDocAEATIeEFALSM09myhnJlPWNNQh1x/sW1i4iZvnt7p3Il138dn250OACAJsLVAACgJZQqrofG8uqOR7S9M97ocFBn69MYs2AjAAAgAElEQVQxXbchrftH8jo8QdNmAECAhBcA0BIeHS8oX3Zd3JOkKXObeu7GjNYlo/rPgWmVGcAKACASXgBAC5golPXYVFHbOmJaw6jMbSsWMb1wS0Yj+bLuGco1OhwAQBMg4QUArGrurgdG80pETBeuTTY6HDTYBd0JbeuM6RunZpQvM4AVALQ7El4AwKp2bLqk8UJFz1qbUDxCU+Z2Z2Z60eYOzZRct57ONjocAECDkfACAFatfLmiR8byWpeMalMm1uhw0CQ2dcS1tyepO85klWOaIgBoayS8AIBV68BYQWWX9vYkZQxUhSov2JSRS3pknBGbAaCdkfACAFal4VxJJ2ZK2t0dVydz7mKOtcmorupP68RMSROFcqPDAQA0CFcIAIBVp+KuB0fzSkdNu7sSjQ4HTeo5G9KKR4KWAM40RQDQlkh4AQCrzrGpoqZLrr09SUUZqAoLSMUi2tOd0HC+rKEcVV4AaEckvACAVaVUcR2aKGpdMqq+FHPu4uy2d8aViZkOjBVUocoLAG2HhBcAsKocmSyoUHE9c02CgapwThEzXbgmqalSRcenS40OBwBQZyS8AIBVI1+u6OhkURvSUa1JUt3F4mxIR7U2EdHB8YLKFaq8ANBOSHgBAKvGoYmiKi5duCbZ6FCwipiZLlybVL7iOjZdbHQ4AIA6IuEFAKwKo/myjk0VtbUjpg6mIcISrUtG1ZuM6vBEkSovALQRrhgAAKvC/56cUcSkPWuYhgjLs2dNQgWqvADQVkh4AQBN79RMSQ+O5rWjK65UlH9dWB6qvADQfrhqAAA0vZtPTCsdNe3qorqL80OVFwDaCwkvAKCpHZ0s6MhkUc/ZmFE8wjREOD9UeQGgvZDwAgCalrvr6ydm1B2P6Iq+VKPDQYugygsA7YOEFwDQtA6MFXRqpqTnb8ooRnUXK4QqLwC0DxJeAEBTcnd989SMelNRXbyOeXexsqjyAkB7IOEFADSlQxNFDebKum59WhGjuouVRZUXANoDCS8AoCndenpG3YmI9lLdRY1Q5QWA1hdrdAAAAMx1bKqogemSXrq1Q1Gqu6iR6irvto64olX9xPcP5Zb1mpcxuBoANBUqvACApvOt0zPKxEz7ekkeUFsXUOUFgJZGwgsAaCqnZko6PFHUVf1p5t1FzfUko+pJRnRksqiK05cXAFoNCS8AoKncenpGyYgx7y7qZk93Qvmy68R0qdGhAABWGAkvAKBpjOTKOjBW0OX9KaVi/ItCffQmo+qOR3R4siCnygsALYWrCQBA07jtzIyiJl3dn250KGgjZqbd3QnNlFynslR5AaCVkPACAJrCZKGs+0by2tebUkecf0+orw3pqDpipsMTRaq8ANBCuKIAADSF289k5S5ds57qLupvtso7WaxoMFdudDgAgBVCwgsAaLhsqaL9wznt7UlqbTLa6HDQpjZlYkpFTYcn6MsLAK2ChBcA0HB3DmZVrEjXbaC6i8aJmGl3V0JjhYpG85VGhwMAWAEkvACAhipWXN8ezOmCNQn1p2ONDgdtbktHTImI6dBkodGhAABWAAkvAKChHhzNK1t2XcPIzGgC0YhpZ1dcw7myxgv05QWA1Y6EFwDQMO6uuwaz6k9Fta2T6i6aw/bOuGImHZ6gygsAqx1XFwCAhhmYLulMtqyXb+uUmTU6nKawfyjX6BDaXixi2tEV16GJoqaKFXUyTRYArFqcwQEADXPnYFapqOnidclGhwI8xY7OhKJUeQFg1aPCCwBYtuVWIy/rS2miUNYjYwVdsz6teITqLppLImra2hHX41NFXbCmokyMGgEArEacvQEADXF3mCxf3pdqcCTA/HZ1xSVJRyeLDY4EALBcJLwAgLorVVz7h4OpiNYmo40OB5hXKhbRlo6YBqaKypeZlxcAViMSXgBA3T04mle25Lqyn+oumtuuroQqosoLAKsVCS8AoK5mpyLqS0W1ozPe6HCAs+qIR7QxE9PjU0UVK97ocAAAS0TCCwCoq7FCRaezZV3Zn2IqIqwKe7riKrv0OFVeAFh1SHgBAHX12GRRyajp4h6aM2N16EpE1Z+K6uhUQSWqvACwqpDwAgDqJleu6HS2pH3rkkpEqe5i9djdnVCxIg1MU+UFgNWEhBcAUDfHpopySVf2pxsdCrAkPcmo1iWjOjJZVMWp8gLAakHCCwCoi4q7jk2V1J+KMhURVqXd3XHly67j06VGhwIAWCQSXgBAXZyaKalQce3oYmRmrE69yai64xEdmSxQ5QWAVYKEFwBQF8emi8rETL1Ud7FKmZl2dyc0U3KdzlLlBYDVgIQXAFBzU8WKRvMVbe2IMxURVrUN6ag6YqbDE0U5VV4AaHokvACAmhuYLsokbemINToU4LzMVnknixUN5sqNDgcAcA4kvACAmqq46/h0UevTUSWj/NvB6rcpE1Mqajo8UaDKCwBNjisPAEBNnc6WVKxI2zoZrAqtIWKm3V0JjRWCpvoAgOZFwgsAqKljUyWlowxWhdaypSOmRMR0aLLQ6FAAAGdBwgsAqJnpYkUj+bK2djJYFVpLNGLa2RXXcK6s8QJ9eQGgWZHwAgBqhsGq0Mq2d8YVM+nwBFVeAGhWJLwAgJqouGsgHKwqxWBVaEGxiGlHV1yns2VNFenLCwDNiCsQAEBNzA5WtbWDwarQunZ0JhSlygsATYuEFwBQEwNTJaWipr4Ug1WhdSWipq0dcZ2cKWmmRJUXAJoNCS8AYMXNlCoazpe1rYPBqtD6dnUFrRiOThYbHAkAYK66Jbxm9nIzO2BmB83sN+d5PGlmnw0fv83MdlY9ts/MvmVmD5jZfWaWqlfcAIClOzYVXPgzWBXaQSoW0ZaOmAamivTlBYAmU5eE18yikj4i6RWS9kr6ETPbO2exN0oadfcLJH1Q0vvD58YkfVrSm939YknXS+IrVABoUhV3HZ8uaX0qqlSMhkRoD7u6EqpIuuNMttGhAACq1OtK5BpJB939sLsXJN0g6dVzlnm1pE+Gtz8n6cUWtIN7maR73f0eSXL3YXdnwjsAaFJnsmUVKq6tnQxWhfbREY9oYyamu4dyytGXFwCaRr3amm2RdKzq7wFJ1y60jLuXzGxcUq+kCyW5mX1FUr+kG9z9j2sfMgBgOQami0pFTf0MVoU2s6crrlMzJf37Y1Pasyax5Odf1kePLQBYafWq8M43YokvcpmYpOdJ+rHw9/eb2YuftgKznzOzO83szsHBwfONFwCwDNlSRUO5srZ0xBisCm2nKxFVfyqqo1MFlSpzL3MAAI1Qr4R3QNK2qr+3Sjqx0DJhv901kkbC+2929yF3n5F0o6Qr5q7A3f/a3a9y96v6+/tr8BYAAOdyfLokSdrC3LtoU7u7EypWgpYOAIDGq1fCe4ekZ5jZLjNLSHqdpC/MWeYLkl4f3n6NpJvc3SV9RdI+M8uEifB3SnqwTnEDABbJ3XV8uqjeZFQZBqtCm+pJRrUuGdWRyaIqTpUXABqtLlck7l6S9BYFyetDkv7R3R8ws/eY2avCxT4mqdfMDkr6VUm/GT53VNIHFCTN+yV9293/vR5xAwAWbyRfVrbsTEWEtre7O6582Z9o8QAAaJy6XZW4+40KmiNX3/fOqts5Sa9d4LmfVjA1EQCgSQ1MlxQzaUOahBftrTcZVXc8oiOTBW3piClCf3YAaBjanAEAzlux4jo9U9LmjriiES7u0d7MTLu7E5opuU5nqfICQCOR8AIAztvJ6ZIqEs2ZgdCGdFQdMdOhiaKcvrwA0DAkvACA8zYwXVRXPKI1CebeBaQnq7xTxYoGc+VGhwMAbYuEFwBwXiYKZU0UK9pKdRd4ik2ZmFJR06GJAlVeAGgQEl4AwHkZmC7JJG3KMPcuUC1ipt1dCY0XKhrJU+UFgEYg4QUALFvZXSdnitqQjikRZbAqYK4tHTElIqbDE8VGhwIAbYmEFwCwbGeyJRUr0tZOmjMD84lGTDu74hrOlzVeoMoLAPVGwgsAWLaBqZJSUVNvksGqgIVs74wrZtLhiUKjQwGAtkPCCwBYlvFCWcP5srZ0xGRGc2ZgIbGIaUdXXKezZU0VK40OBwDaCgkvAGBZ7hvOS5K2djBYFXAuOzoTilLlBYC6I+EFACyZu+vekZx6k1GlY/wrAc4lETVt7Yjr5ExJMyWqvABQL4wyAgBYsscmi5ooVHRpb7LRoaBN7B/KNTqE87arK67Hp4o6OlnU3h6OHQCoB76WBwAs2T3DOaWipvVpvjcFFisVi2hLR0wDU0Xly1R5AaAeSHgBAEuSLVX0yHhBe3uSijJYFbAku7oSqkg6Osm8vABQDyS8AIAleXA0r7JLl/amGh0KsOp0xCPamInp8amiihVvdDgA0PJIeAEAS3LvcE4b0lFtyNCcGViOPV1xlT3oCw8AqC0SXgDAop2aKel0tqx9VHeBZetKRNWfiuqxqYJKVHkBoKb4eh4AsOgRcB8czSsiqeKtMWou0Ci7uxO67UxWA9NF7exKNDocAGhZVHgBAItSdtfJmaI2ZGKKRxisCjgfPcmo1iWjOjJZVMWp8gJArZDwAgAW5fRMScWKtLWDxkHAStjdHVe+7Do+XWp0KADQsrhqAQAsyvHpktJR07pktNGhAC2hNxlVdzyiI5MFbemILbubwGV99KkHgIVQ4QUAnNNMqaLhfFlbOuIy5t4FVoSZaXd3QjMl1+ksVV4AqAUSXgDAOR2fDqZP2UJzZmBFbUhH1REzHZooyunLCwArjoQXAHBW7kEfw95UVOkY/zaAlTRb5Z0qVjSYKzc6HABoOVy5AADOajhfVq7sDFYF1MimTEypqOnQRIEqLwCsMBJeAMBZDUyVFI9IG9IkvEAtRMy0uyuh8UJFI3mqvACwkkh4AQALKpSDwXQ2Z+KKMFgVUDNbOmJKREyHJ4qNDgUAWgoJLwBgQSdminIxWBVQa9GIaWdXXMP5ssap8gLAill0wmtmb1/g/l9duXAAAM3C3TUwXVJ3PKLuBHPvArW2vTOumEmHJwuNDgUAWsZSKrzvXOD+31mJQAAAzWWsUNFUsaJtnfFGhwK0hVjEtKMrrtPZsqaKlUaHAwAt4Zxt1MzsReHNqJm9UFJ1J67dkiZrERgAoLEGpoqKWjCCLID62NGZ0NHJog5PFLSvN9XocABg1VvMVczHwt8pSX9bdb9LOiXpF1c6KABAYxUrrpPZkrZkYopFGKwKqJdE1LS1I67Hp4q6YE1FGea+BoDzcs6E1913SZKZ/Z27/2TtQwIANNqJ6aIqLm2lOTNQd7u6goT36GRRe3uSjQ4HAFa1RX9tWJ3smlmk+qc2oQEAGsHddSwcrGoNg1UBdZeKRbSlI6aBqaLyZfryAsD5WMoozVeY2bfMbFpSMfwphb8BAC1inMGqgIbb1ZVQRdLRSS6zAOB8LGUkkk9K+qKkn5Y0U5twAACNdmyawaqARuuIR7QxE9PjU0Xt7k4oTl96AFiWpVzN7JD02+7utQoGANBYxYrr5ExJmxmsCmi4PV1xnZop6bHJoi5Yk2h0OACwKi2l/+3nJb2sVoEAABrv5HRJFRfNmYEm0JWIqj8V1WNTBZUq1BsAYDmWUuFNSfq8mX1DwXRET2D0ZgBY/YLBqooMVgU0kd3dCd12JquB6aJ2dlHlBYClWkrC+2D4AwBoQeOFiiaLFaZBAZpITzKqdcmojkwWtb0zrojR1QAAlmLRCa+7v7uWgQAAGmsgHKxqM4NVAU1ld3dcdw7mdHy6RHcDAFiiJc2ha2YvNbOPmdkXw7+vMrMX1SY0AEC9lMLBqjYxWBXQdHqTUXXHIzoyWVCFsUMBYEkW/TW+mf2ipF+W9DeSXhPenZX055Keu/KhAQDq5cRMSWWXtnXUp3q0fyhXl/UArcDMtKc7obuHczo1U9LmOh2nANAKllLhfaukl7j7H0mqhPc9LOmZKx4VAKBu3F3HporqikfUnVhSwx8AdbI+HVVHzHR4sihmiASAxVvKlU2XpGPh7dkzbVxSYUUjAgDU1fHpkiaLFW3rjMsYEAdoSmam3d0JTRUrGsyVGx0OAKwaS0l4/0fSb86575ck/ffKhQMAqLe7BrOKMVgV0PQ2ZWJKRU2HJgpUeQFgkZaS8P6ipO83s6OSuszsgKTXSvrVWgQGAKi9yWJZB8YK2toRZ7AqoMlFzLS7K6HxQkUjeaq8ALAYS5mW6KSZXS3pakk7FDRvvt3dK2d/JgCgWe0fyqkiaXsXg+AAq8GWjpgOThR0eKKo3hStMgDgXJYySvNlkobd/XZJt4f3bTOzde5+T60CBADURrni2j+U057uuDIxBqsCVoNoxLSzK65Hxgsay5e1NhltdEgA0NSWcoXzaQWDVFVLSPrUyoUDAKiXh8fymi65ruxPNzoUAEuwvTOumEmHJxk3FADOZSkJ73Z3P1x9h7sfkrRzRSMCANTFt4dy6klGtIvmzMCqEouYdnTFdSZb1mSRvrwAcDZLSXgHzOyK6jvCv0+sbEgAgFo7NVPS8emSruhLMxURsArt6EwoatLhiWKjQwGApraU0Q4+KOnfzOyPJR2StEfS2yW9rxaBAQBq567BrOIR6ZLeZKNDAbAMiahpW2dcRyeL9OUFgLNYyijNHzWzMUlvlLRNwSjNb3P3z9UqOADAypspVvTgaF6X9qaUijJYFbBa7eyK67HJom49ndXLt3c2OhwAaEqLSnjNLCrp9yS9z93/qbYhAQBq6Z7hnMouXdGXanQoAM5DKhrR1o6Y7hvJ6Ts2ptWVoMoLAHMt6qt9dy9L+gVJdBQBgFWs4q67h3La0RlXX5o5PIHVbld3QhWXbj+TbXQoANCUltKW7ZOS3lyrQAAAtffoeEETxYqu7Ke6C7SCTCyivT1J7R/OaaZUaXQ4ANB0lpLwXiPpz8zsqJn9r5n9z+xPrYIDAKysuwZz6o5HdMGaRKNDAbBCnrMhrWJFunOQKi8AzLWU9mwfDX8AAKvQYLakx6eKun5zRhGmIgJaRl86pgvXJHTXYE7Xrk8ryWB0APCEpYzS/MlaBgIAOH/7h3ILPnbfcE5RkyKysy4HYPV57saMHjkwpruHcrpuQ6bR4QBA01j0V4AW+Fkzu8nM7g3ve4GZ/VDtwgMArIRcqaITMyVt6YgrEaW6C7SajZmYdnXFdfuZrIoVb3Q4ANA0ltLm5T0K5uD9a0nbw/sGJP3GSgcFAFhZj00V5Qrm7QTQmp6zMaOZkuueYVpwAMCspSS8b5D0Sne/QdLsV4dHJO1e6aAAACunVHEdmypqYzqmTIy+fUCr2t4Z19aOmG47nVWJKi8ASFpawhuVNBXenj2LdlbdBwBoQsemiiq5tKub6i7Q6p63MaPJYkX3jVDlBQBpaQnvlyV9wMySUtCnV9LvS/piLQIDAJy/irsemypqXTKqNYloo8MBUGM7uuLanInpW6eyKlPlBYAlJby/ImmTpHFJaxRUdneIPrwA0LROzpSUK7t20XcXaAtmpu/YmNFEsaL7R/ONDgcAGu6c0xKZWUbS70h6tqS7Jf2ypA2Sjrn7qdqGBwBYLnfXkcmiOmMR9aWo7gLtYnd3XBszMd1yakbPXpdUlHm3AbSxxVR4PyzpeyU9LOkHJf2au99BsgsAzW0oV9ZUsaKd3XEZF7xA2wiqvGmNFyp6cIQqL4D2tpiE9xWSXubuvx7efmVtQwIArISjk0Ulo6bNmXM25gHQYi7oTmh9OqpbTs+o4vTlBdC+FpPwdrj7SUly92MK+u8CAJrYRKGs4XxZOzrjilDdBdrObF/e0XxFD9GXF0AbW8zX/jEze6EkW+BvuftNtQgOALA8RyaLipq0rZPBqoB2deGahPpTUd1yOquLepJ8+QWgLS0m4T0j6W+r/h6e87dL2r2SQQEAli9bqujUTEk7OuOKR7jABVrd/qGF59zd3BHTPcN5feXYlDZlnvoF2GV9qVqHBgANd86E19131iEOAMAKOTpZlBTMxwmgvW1Mx3QwVtChiaI2pmMMYAeg7SxlHl4AQJPLlys6Nl3UpkxM6RineKDdmZn2dCc0VazoTLbc6HAAoO64GgKAFnJ0sqiKS3u6E40OBUCT2JiJKRMzHZwoyBmxGUCbIeEFgBYxU6ro8amgutsR5/QOIBAJq7yTVHkBtCGuiACgRdxxJqsy1V0A89gUVnkfpcoLoM2Q8AJAC8iWKrprMKeN6Zg6qe4CmCNipgvCvrynqfICaCNcFQFAC7hzMKtCxbWnm5GZAcxvUyamjpjp4DhVXgDtg4QXAFa5XLmiOwdzunBNQl2JaKPDAdCkzEwXrEloqlTRyZlSo8MBgLog4QWAVe6uwZzyZddzN2YaHQqAJjfb7eHgREEVqrwA2kDdEl4ze7mZHTCzg2b2m/M8njSzz4aP32ZmO+c8vt3Mpszs7fWKGQCaXb5c0R1nstrTHdfGTKzR4QBocmamZ3QnNFNyPTCSb3Q4AFBzdUl4zSwq6SOSXiFpr6QfMbO9cxZ7o6RRd79A0gclvX/O4x+U9OVaxwoAq8ndQznlyq7voLoLYJHWp6Pqjkf0zVMzKlPlBdDi6lXhvUbSQXc/7O4FSTdIevWcZV4t6ZPh7c9JerGZmSSZ2fdJOizpgTrFCwBNr1B23XYmq11dcW3uYLAqAIsz25d3rFDR/VR5AbS4eiW8WyQdq/p7ILxv3mXcvSRpXFKvmXVI+g1J765DnACwauwfzilboroLYOn6U1FtzMR0y6kZlStUeQG0rnolvDbPfXPPrgst825JH3T3qbOuwOznzOxOM7tzcHBwmWECwOpQrLhuOz2jHZ1xbe2kugtgacxMz9+Y0Xihovuo8gJoYfVKeAckbav6e6ukEwstY2YxSWskjUi6VtIfm9lRSW+V9A4ze8vcFbj7X7v7Ve5+VX9//8q/AwBoIt8ezGq65HreJqq7AJZnd3dcm8Mqb4kqL4AWVa+E9w5JzzCzXWaWkPQ6SV+Ys8wXJL0+vP0aSTd54PnuvtPdd0r6kKQ/cPcP1yluAGg6+XJFt54O+u5uo7oLYJnMTM/flNFEsaL9w7lGhwMANVGXhDfsk/sWSV+R9JCkf3T3B8zsPWb2qnCxjynos3tQ0q9KetrURQAA6Y4zOWXLrhdsproL4Pzs7IprW2dQ5S2UqfICaD11m7TR3W+UdOOc+95ZdTsn6bXneI131SQ4AFglsqVg3t0L1yS0KUN1F8D5MTNdv7lDn3pkXHcOZvVcBsED0GLq1aQZALACbjudVb7iej59dwGskC0dcV3QndBtZ7LKliqNDgcAVhQJLwCsElPFiu4czGpvT1L96bo10AHQBl6wOaN82XXb6WyjQwGAFUXCCwCrxLdOz6jsoroLYMWtT8d0cU9Sdw5mNVWkygugdZDwAsAqMF4oa/9QTvt6k+pJRhsdDoAW9LxNGVVc+uapmUaHAgArhoQXAFaBW8ILUAaUAVArPcmoLu1L6Z6hnEbz5UaHAwArgoQXAJrcSK6se4fzuqwvpTUJqrsAauc7NmYUMekbJ6nyAmgNJLwA0OS+cWpGsYj0nA1UdwHUVmc8oiv703pgNK8z2VKjwwGA80bCCwBNbDBb0oOjeV3Zl1ZnnFM2gNq7bkNayajpf6jyAmgBXD0BQBO7+cSMklHTtRvSjQ4FQJtIxyK6dn1aB8cLOj5dbHQ4AHBeSHgBoEkdmyrq4ERB161PKx3jdA2gfq7qTysTM918Ykbu3uhwAGDZuILC/2/vzuMju8s733+fU3tpKbXUanfLva9ecRtvSQzGYYfAGBgbTGaCk4EQ8gqTm8vM3CR37iQMEzIkMxkyL7LwSsLuJGYLHhMMZosxeMFbt5e2e9/Um1r7VnvV7/5xjhq1rG6r1VKdqtLn/XqpVXXOqaqn6nSVzlPP7zw/AHXIOacHT0yqNebp+hVUdwHUVjxi+oWVaR2dKOnwOFVeAI2LhBcA6tC+0aKOT5b16pVpxTwLOxwAS9D2rqTa4x5VXgANjYQXAOpM1Tn96ERWXYmIru5KhB0OgCUq6plevTKtU7my9owUww4HAOaFhBcA6syzgwUNFip6TU9anlHdBRCeKzsT6kpG9NDJrKpUeQE0IBJeAKgjparTT05mdWlLVFsy8bDDAbDEeWa6ZVVaQ4WKnhsqhB0OAFwwEl4AqCNPns5polzVrT0tMqq7AOrA1kxcq9JRPXwyq3KVKi+AxkLCCwB1Ileu6rG+nDZn4lrTGgs7HACQJJmZXrMqrbFSVTsG8mGHAwAXhIQXAOrEI6eyKladXrMqHXYoAHCW9e1xrWuN6dG+rAqVatjhAMCcRcMOAACa2c45VkOy5aqe7M+rpyWq45Nldaf4eAZQX17Tk9YX947qyf68bl7JF3MAGgMVXgCoA/tHizKTNrfTqApAfeppiWlLJq7H+3LKlanyAmgMJLwAELLRYkUnsmWta40pFeVjGUD9umVVWoWq06N9ubBDAYA54cgKAELknNOekaJinrSR6i6AOtediuqqzoSe7s9pvFgJOxwAeFkkvAAQov58RUOFija3xxXzmIYIQP171cq0qpIePkWVF0D9oysKAISk6pz2jBTUEjWmIQLQMA6Pl7S6JaZnBvNqi3lqic2tfrJ9eXKRIwOAl6LCCwAh6Z0oabLstK0jIc+o7gJoHJvaY/JM2jdWDDsUADgvEl4ACEGp6rR/rKjORETdyUjY4QDABUlEPK1rjelUtqwxzuUFUMdIeAEgBAfHiipVpW0dcRnVXQANaEN7XDFP2jtKlRdA/SLhBYAay5arOjJeUk86qkyc6i6AxhTzTBva4hoImu8BQD0i4QWAGts3WpRM2pJhGiIAjW1da0wJz7R3pCDnXNjhAMBLkPACQA2NFCo6mS1rfVtMqSgfwQAaW8QzbcrENVKsqj9PlRdA/eFoCwBqxAXTEMU908Y2qrsAmsPqlqjSUdO+0SJVXgB1h4QXAGqkL1fRcLGqzZm4oh6NqmEVmWYAACAASURBVAA0B89Mm9vjGi9VdTJbDjscADgLCS8A1ECl6ld3W2OeVrdEww4HABbUqnRUbTFP+8eKqlLlBVBHSHgBoAYOT5SUqzhd3hGXxzREAJqMmWlLJq5s2enYJFVeAPWDhBcAFlm+XNXBsaIuSUXUlaS6C6A5dScj6oh7OjBaVKVKlRdAfSDhBYBFtne0KOekbR2JsEMBgEVjZtrakVCh6nRkohR2OAAgiYQXABbVcKGiE8E0RGmmIQLQ5DoTES1PRnRovKgSVV4AdYCjLwBYJM457R4pKOGZNrYzDRGApWFrJq5SVTo0Xgw7FAAg4QWAxfL8UEGjxaq2djANEYCloz0e0cp0VEfGSypUqmGHA2CJo3sKAMzBzoH8BW1frjr9+GRWmbinnjQftQCWli3tcfVlyzowVtIVy+hfACA8VHgBYBEcGCuqUHW6vCMhYxoiAEtMS8zTpS1R9U6UlC1T5QUQHhJeAFhg2XJVh8dL6klH1ZGIhB0OAIRic3tcZtL+Uc7lBRAeEl4AWGC7hwvyTNraQaMqAEtXMuppbWtMJ7JljZcqYYcDYIki4QWABXQ6V9bpfEWb2+NKRviIBbC0bWyLK2rSPqq8AEJCJxUAWCCVqtOLwwW1Rj2ta4uFHQ4AhC4eMW1oi2vfWFEPHp+c12ke25cnFyEyAEsFCS8ALJCD40XlKk43dCfkXWSjqgvtCg0A9WpdW0xHJkraO1rUDd1JGvkBqCnG2wHAApgsVXVwrKRV6ai6knyXCABTop5pY3tMQ4WKBgucywugtkh4AeAiOef0wkhBEZMuo1EVALzE2taYkhHT3pGinHNhhwNgCaEMAQAXqS9X0WC+oss74krQqApAg6jlqROembZk4npuqKC+XEUr0xyCAqgNjswA4CKUq04vjhTUFvO0ppVGVQBwLj3pqFqjnvaNFlSlygugRkh4AeAi7B8rqlBxumLZxTeqAoBmZkGVd7LsdGKyHHY4AJYIEl4AmKfxYkVHxkta3RLVsnlMtQEAS82KVESZuKf9Y0VVqPICqAESXgCYh6lGVVFP2ppJhB0OADQEM9PWTFz5ilPvRCnscAAsASS8ADAPxybLGi5UtTWTUDzCUGYAmKuuZFRdiYgOjBVVrlLlBbC4SHgB4ALlK1XtGSmoMxHR6hY6jQLAhdqaiatUlQ6PU+UFsLhIeAHgAjjn9MJwQVVJVy5LyGhUBQAXLJOI6JJURIfGiypWqPICWDwkvABwAfpyFZ3OVbS5Pa6WGB+hADBfWzIJVZx0cKwYdigAmhhHawAwR8WKX91tj3la38acuwBwMVpjni5NR3V0oqRcuRp2OACaFAkvAMzRnpGCSlWnqzqZcxcAFsLmTFyStH+UKi+AxUHCCwBzMJAv63i2rA1tMbXHmXMXABZCKuppXVtMx7NljRcrYYcDoAmR8ALAyyhWnHYNFZSOmjYF1QgAwMLY2B5XzJP2UOUFsAhIeAHgZTx0clK5itNVnUlFGMoMAAsq5pk2tsc1kK9oIF8OOxwATYaEFwDO4/hkSU/257WmJarOBEOZAWAxrGuNKRUx7R0pyjmmKQKwcEh4AeAcihWnfz4yrvaYp20dibDDAYCm5ZlpSyausVJVJ7NUeQEsHBJeADiHfzkxqeFCVW9b16aox1BmAFhMq9JRtcc87R0tqkKVF8ACIeEFgFkcGC1qx0BeN65IaS1z7gLAojMzbeuIK19xOjpeCjscAE2ChBcAZsiWq7r/6Li6kxHdsioddjgAsGR0JaNanozowFhRxQpVXgAXj4QXAKZxzuk7RyeUrzi9fT1DmQGg1rZl4io76eA40xQBuHgkvAAwzXNDBe0dLeqWVWmtSEXDDgcAlpy2eESXpqM6Ml5StlwNOxwADY6EFwACI4WKvn9sUmtao7phRSrscABgydqcictM2jtKlRfAxSHhBQBJVedPQSRJv7S2TZ4xlBkAwpKKetrQFtOpbFlHJ2hgBWD+SHgBQNLjp3M6NlnWG1a3qCMRCTscAFjyNrbFlYyYvtc7oSrTFAGYJxJeAEve8cmSHjqR1baOuK7qTIQdDgBAUsQzXdaRUH++oh0D+bDDAdCg6MgCYMnYOcsBU7Hi9EhfVomIaXVLTM8MFkKIDAAwm0tSEa1rjenHJ7O6fFlC6Si1GgAXhk8NAEuWc07PDeVVqDht70oqxhREAFBXzEyvX92iQsXpoRPZsMMB0IBIeAEsWYcnSurPV7StI64M5+0CQF3qTkV1XXdSOwfzOpUthx0OgAZDwgtgSRopVLR3pHhmuBwAoH69alVa6ajpu70TcjSwAnABSHgBLDnFitPOwbySEdNVnUkZUxABQF1LRjzd2tOiE9mynh+i1wKAuSPhBbCkcN4uADSmqzsT6klH9eCJSeUr1bDDAdAgSHgBLCmctwsAjcnM9IbVLZosOz18kgZWAOaGhBfAkjHMebsA0NBWtcR0TVdCT/bTwArA3NQs4TWzN5vZHjPbb2a/N8v6hJl9OVj/UzNbHyx/g5k9ZWbPBb9fW6uYATSP8VJFOwfySkU5bxcAGtmtPS1KR033Hx1XpUoDKwDnV5OE18wikv5S0lskXSHpvWZ2xYzN3i9p2Dm3WdInJf1JsHxA0tudc1dLukvSl2oRM4DmUa463XtoXGXndC3n7QJAQ0tFPb15batO5yp6tC8XdjgA6ly0Ro9zo6T9zrmDkmRm90i6TdIL07a5TdJHg8tfk/QXZmbOuR3TttklKWlmCeccLfqAJWrnQP6Ctt81lNfxybK2dyXVFue8XQBodFsyCV25rKhHTmW1JRPXJelaHdICaDS1GtJ8qaTeadePBctm3cY5V5Y0Kqlrxjb/WtIOkl0Ac9U7UVLvZFkb2mJayQERADSN169uUSpq+tbRcVWYmxfAOdQq4Z1t/ODMT6bzbmNmV8of5vwbsz6A2QfN7Ekze7K/v3/egQJoHiOFil4YLqgrGdHWTDzscAAACygV9fSmNf7Q5scY2gzgHGqV8B6TtGba9dWSTpxrGzOLSspIGgqur5b0DUnvc84dmO0BnHN/45y73jl3fXd39wKHD6DRFCpV7RjMKxkxXUOTKgBoSls7ErpiWUIPn8rqdI6uzQBeqlYJ7xOStpjZBjOLS7pT0n0ztrlPflMqSbpd0g+dc87MOiR9S9LvO+cerlG8ABpY1TntHMirVHW6dnlS8QjJLgA0q9evblEyYvrWEYY2A3ipmiS8wTm5H5b0gKQXJX3FObfLzD5mZv8q2OwzkrrMbL+kj0iamrrow5I2S/ovZrYz+FlRi7gBNKbdI0UNF6u6allC7TSpAoCmlg6GNvflKvopQ5sBzFCzDi7Oufsl3T9j2R9Mu5yXdMcst/sjSX+06AECaApHx0s6OlHS+raYelpiYYcDAKiBbR0JXd5R0E9OZbU5E9eKFE0KAfhqNaQZABbdYL6sF0cK6k5GtI0mVQCwpLxhdatSEdO9h8ZVqFTDDgdAneDrLwBNYbJU1c7BvFqinq7pokkVADSLC5l7/cplST3Rn9M/7hvTXdsy/C0AQIUXQOMrVZ2eHvDP23pld1JRjwMcAFiKOpMRbcnEdSpX1lP9c0+UATQvEl4ADa3qnHYO5pUtO13blVI6yscaACxlG9piWpGM6IfHJ3V8shR2OABCxpEhgIa2Z6SowXxFVy5LqDNJR2YAWOrMTFd3JdUW93TvoXFlS5zPCyxlJLwAGlbvRElHJkpa1xrT6lY6MgMAfDHP9M4N7cqWq7rvyLiqzM8LLFk0rQLQkAbzZb0wXNDyZESXddCRGQBwtlPZsi7rSGjXcEFfPzimLZnEnG63fXlykSMDUEtUeAE0nKmOzOmo0ZEZAHBOq1ui6klHdWCspP5cOexwAISAhBdAQ8lXqmc6Ml/XnVKMjswAgHMwM125LKHWmKdnBvMaL1XCDglAjZHwAmgYVed03+FxZctO27uSdGQGALysiGe6bnlSETM91Z9XvkITK2Ap4WgRQMP4l+OTOjhW0uXLEupK0oIAADA3qain67qTKlWdnurPq1yliRWwVJDwAmgIzwzm9UR/Xtd1J7WWjswAgAvUHo/o2uVJTZSq2jGQp3MzsESQ8AKoe0cnSnqgd0Lr22J63aUtYYcDAGhQy5NRXdmZ0GChoueHCnIkvUDTY0wggLo2UqjoG4fG1BGP6B3r2+TRkRkAcBFWt8SULzvtHysqFbU5T1cEoDGR8AKoW/lyVV89MCbnpNs3titJkyoAwALY1B5TvlLVgbGSkhFPazhVBmhaJLwA6lKl6vRPh8Y1XKzozk0ZdSYjYYcEAGgSZqYrliVUqDjtGi4o5plWpjksBpoR5RIAdcc5p+/0TujoRElvXduqtW188w4AWFiema7pSqoj7s/R25cthx0SgEVAwgug7jzSl9NzQwXdvDKlqzqTYYcDAGhSUc90fXdK7XFPOwfzOp0j6QWaDQkvgLqyayivH5/M6splCb1qZTrscAAATW4q6W2LedoxkNfBsWLYIQFYQCS8AOpG70RJ9x+d0JrWqN6ytlVGR2YAQA3EgqS3Nebp6wfHdJikF2gaJLwA6sJQvqJ/OjimTDyid21oV9Qj2QUA1E48YrqhO6XORERfOzimI+MkvUAzIOEFELrJUlVfOTAqSbpjU7tSTD8EAAhBPGK6c3NGHUHS2ztRCjskABeJo0oAoSpU/GR3slzV7ZvatSzB9EMAgPC0xDzduTmj9lhEXzkwqsNUeoGGxoRjAELzVH9OT/XnNVyo6JXLk+rPVdSfq4QdFgBgiWuNefrlLRnds39UXz0wpndtaNemTDzssADMAxVeAKGoOqdnB/MaKlR0dWdC3Sm+fwMA1I+WmKf3bsloeTKirx8a056RQtghAZgHEl4ANeec0/eOTaovV9G2jrh6WmJhhwQAwEuko57euzmjlamo7j00rheGSXqBRkPCC6DmfnIqqx0DeW1oi2lDG0PEAAD1Kxn19J7N7VrdGtU3D4/r2cF82CEBuAAkvABq6un+nB4+ldPVnQlt5XwoAEADSEQ8vXtTRuvaYrr/6ISe7s+FHRKAOSLhBVAzzw3m9d1jk9rUHtNb1rbKjLl2AQCNIeaZbt/Yrk3tMX332KQe68uGHRKAOSDhBVATzw3m9a2jE1rfFtM7NrTLI9kFADSYqGd614Z2Xd4R14MnsnrwxKScc2GHBeA8aIsKYNE9P+Qnu+taY/rXG9sV80h2AQCNKeKZ3r6+TYneST3Wl1Oh4vTG1S16ZnB+Da22L08ucIQApiPhBbConh/K65+P+Mnu7ZtIdgEAjc8z05vWtCgZMT122k96V7dEGb0E1CESXgCLhmQXANCszEy3XtqiRMT0o5NZ9efK2t6VVIS/dUBd4RxeAIvi+aG8vkWyCwBocj+/Mq03rm5Rf76iJwdyKlc5pxeoJyS8ABbcs4N+sruGZBcAsAS8sjulV3QmNFKo6vHTORUq1bBDAhAg4QWwYJxzeuRUVvcfndC6tpjuINkFACwRPS0xvXJ5UpPlqn56OqdsmaQXqAckvAAWRNU5fe/YpB46mdUVyxK6g27MAIAlpjsV1Q3dKZWqTo/15TRWrIQdErDkkfACuGjlqtO9h8b19EBeN65I6e3rWmnaAQBYkjoSEd20Ii3PpJ+ezmkwXw47JGBJI+EFcFHy5aru2T+qvaNFve7SFr320hYZ0zIAAJaw1pinn1uRUiri6cn+vE5lSXqBsJDwApi3sWJFd+8b1clsWbetb9MNK1JhhwQAQF1IRj3ddElKmbinnYN5HRkvhh0SsCSR8AKYlxOTJX1x76jGi1W9e1O7Ll+WCDskAADqSswz3dCd0opkRC+OFLV7uCDnmLYIqKVo2AEAaDzPDOb13d4JtcY8vXtrRitSfJQAADCbiGe6dnlSL44UdXiipFylqld0Jul1AdQIR6kA5qzinH5wbFJPD+S1vi2m29a3KRVloAgAAOdjZrpiWULpqGn3SFGP9+f0yuVJJSL8DQUWGwkvgDmZLFV17+Ex9U6UdeOKlG7tScujORUAAHO2vi2uVMTTM0N5PdaX03Xd9L4AFhsJL4Azdg7kZ10+Wqxox0BexarTKzoT6kxE9Oxg4cz67cuTtQoRAICGdkk6qpsiKT01kNdjfVmtaY1qXVs87LCApsU4CgDn5JxT70RJP+3LSZJuWpFST0ss5KgAAGhsmUREP3dJSsmIpy8fGNOOgVzYIQFNiwovgFmVq067hgs6mS2rKxHRNV1JxSMMYQYAYCGkg2mLDo0X9UDvpPqyFb1+dYuiNLMCFhQJL4CXGCtWtHMwr2zZaUsmro1tMRnn6wIAsKBinun2je368cmsHu3LqT9f1js3tKs1xiBMYKGQ8AI4wx/CXNaLIwXFI6Ybu1PqTEZe9nbnOvcXAACcn2em1/S0aEUqqvuPjuvze0b0rg1tnEIELBC+PgIgScpXqto5mNcLIwV1JSO6+ZL0nJJdAABw8S5fltCvbO1QxKS/3zeqZwb5MhlYCCS8ANQ7UdJnd4/odK6irZm4rlvO+boAANTailRUv7qtQ2taY/r20Ql9++i4SlUXdlhAQ2NIM7CEVapOPzmV1WN9OWXinm5ckdKyBFVdAADCkop6evemdj10IqvHTud0fLKs29a3qTvFYTswH7xzgCVqMF/WNw9P6FSurFd0JvS61S16cbgYdlgAACwp5+qD0ZGI6LrlST03VNDn9ozo8o6EVrdEzzSR3L48WcswgYZFwgs0ofM1kXLOqXeyrN0jBUVM2t6V1Mp0lGQXAADVVyPG7lRUN6/09OxgQbuGCxrMV3RlZ0Ixpi4C5oyEF1hCCpWqnh8qqD9fUVcyoqs7E0pGOJUfAIB6lYh4ur47qYPjJe0fLWr0VEXXdFHdBeaKhBdYIk7nynp+qKBy1emyjrjWtTK3LgAAjcDMtKk9rs5ERM8M5vXT0zmZSTevTCtKtRc4LxJeoMmVq057RgrqnSyrLebphhVJtcVoTAUAQKNZlojo5pVp7R4u6NG+nPaNFvVLa1u1ijl7gXNiLCPQxEaLFT3Sl1XvZFnr22L6+UtSJLsAADSwmGe6uiupOza2q1Bx+uLeUT14fFJlpi8CZkWFF2hCzrkz5/okIqYbupPqSvJ2BwCgWWzKxPX+1g798PikHjvtV3vfuq5Vl1LtBc5ChRdoMgP58pk/fJekorp5ZZpkFwCAJpSMeHrr2ja9e1O7ilWnu/eO6ru9E8qVq2GHBtQNjoKBJlF1To+fzunHJ7PyTHpFZ0Kr0lEaUwEA0OQ2tsf1gcs79KMTWe0YyOvF4YJu6Unrmq6kPI4DsMSR8AJNoD9X1v1HJ3QyW9bWTFyXtkSVYLohAACWjETE0xvXtGr78qS+d2xCD/ROaudAXm9cwzBnLG0kvEADqzqnx/pyevhUVvGI6bb1bbqsI65nBgthhwYAAEKwIhXVL2/O6MWRov7l+KS+tHdUV3UmdGtPi1pjfBmOpYf/9VgwPz45OadlF3N/5/JXzw++ZNnf7x2Z97L/9czAnOKZbdmf7njpbWd7jE/Mst1s9zdbLJ/YMaATkyV9Yc+IHjqZ1ZZMXLmy0+XLEmcNYf5O78RLbsuy2ZfVWzwsa8xl9RYPyxpzWb3Fw7L6XLZzIH/m5xM7Bs66/sxgQcWK03ipqg1tMe0aKuivdw3py/tHZz3+mOuy2Y63ZttutuOeuR4fzWakWJnTdnO9v4V2Mce7WHwkvFgwD5/KzWnZxdzfuYyVXtqKv3eyPO9lxVl6Pcz1+c3WJmK2x5jNbPc3M5apRhRf3DuqybLTOza06R0b2ud0/wAAYGnZ1pHQq1eltSIV1aHxkiTp4VNZFSoX3thqtuOt2cx23HMxx0djsx2YXcT9LbSLOd7F4mNIM1DHdg7kz/x2zun4ZFl7Rv3hyutbY9qciStfdmdtBwAAMF066umarqQ2tlX0cJ/f4PKp/px+/pK0rl2eDDs8YFGR8AINYLxY0a7hgkaKVXXEPY0Uq7psWSLssAAAQANpi0ckSe/bmtGPTmT1g+OTevy0X50sVKo0vERT4n81UMemhi8/0pfTZLmqq5YldNOKVMhRAQCARtbTEtN7t2R05+Z2dST8dOCvdg3rRycmNVFiDl80Fyq8QB3Klqp6tC+rp/r9IcprW2Pa1B5XPMJcegAAYGGsb4trfVtcn9gxoPVtMT3al9Pjp3O6qpNRZGgeJLxAnfnxyUk9cTqvUtWpJx3V8WxZlzN8GQAALKJ3bmjXcKGix0/n9Nyg/4X7Vw+M6trlKW1sj8kzvnRHYyLhBepAtlw903Dq4VM5bc3EdcuqtI5NlnU8O/uUFQAAAAtpWSKiN61p1atWpvWp54fUl63oawfH1B7ztJ3mVmhQJLxAiPqyZT01kNMLQwWVg07/d23NaFVLTJJ0LKT2+gAAYOlqifnn9f7mVcu0b7SoHf15PXQyK0m699CYti9Pal1rTEbVFw2AhBeosapz2jdalCR9bs+IoiZd1ZnUdd1JfWb3yJlkFwAAIEwRM13WkdBlHQkN5sv62xdHdHi8pN0jRbXHPM71RUMg4QVq5PhkSbuHC9o9UtR40AHxF3vSekVXUqkoDdMBAED96kr6acOrV6bVlyvrRLasR/r8KY0+vWtIPS0xrUpHFfP8qu/UqVqSGA6NUJHwAovEOafjk2XtHilIkr60d1QRk9a3xfT61S36xqFx3XRJOuQoAQDAUjM9Gb3QZRHP1NMSU09LTPlKVQ+eyKpclV4YLmj3cEHLU/5cv+WqU9RjyDPCR8ILLKCRQkWS9K0j4zo8XtJ4qaqpmYTetq5VmzNxJc9M6j4eTpAAAAALYOqY5uaVKY2VqjoxWdapnN9/5IcnJtWdjGhlOqorKgx9RnhIeIF5qjq/y9Qzg3kdHS+pd6KksWCo8r7Rota2xnRrR1ybM3F98tkhlavS7uHiWfcx129TAQAA6pWZKROPKBOP6LKOuB44NqnVLTH1ZcvqyxW0a8gf7fbsYF6b2uNnmmIBtUDCC8xBtuwnsk+czul0rqz+XEUDef8bzG8fnVA6alrTGtNNrTF979ik/q+rO+lcCAAAlpyp458rliV0eUdcw4Wqys7p6YG87j/qT7XYk45qcyYuyT8FjGMmLCYSXkBSxTlNlKoaLVQ1XPSHJd97aEwjwfVCxa/m/uD4pNJR04pUVNcuT+qJ/rw+cFmHupKRMx/W3zs2yQc3AABY8sxMncmIti9P6umBvH51W4cOjBW1f7R4Zpqjv941rLVtMa1rjWldW0zt8UjIUaPZkPCiYc0c+nuu4cGVqlOh6pQPktZ7D40pV3bKV6pnlv2PnYMvue3RiZJSUU8rkhGlo572jBb14as61TptGM4T/XktT/E2AgAAOJepY7RT2bJaop6u6UpqW4ff8CoVNe0ZKej5YNhzOmrqTETUmYjoF1amlQkzcDQFjtTRsKrOT2Knqq9Hxotnrk8t+/6xCZXd2bfbPVKUZ1IyYkoFzRY2tceUjHhKRU3pqKeHTmZ1y6qWs263Z9T/RnImzsMFAAC4MFMNr65dnpJzTuOlqoYKFQ3lKzqVLevYZFnPDhX0bydLkqSHTkxqZTqqVemoWmMeo+kwZzVLeM3szZL+t6SIpL9zzn1ixvqEpC9Kuk7SoKT3OOcOB+t+X9L7JVUk/bZz7oFaxY3aqlSdsuWqskGWumsor2zZXzZROvsnVzk7k31xpCiTlIiYEkFr5J6WmBIRUzJY9mR/Xq/taVHM+9k5Jt/pndCWDN0DAQAAwmBmao9H1B6PaH2bf17vWKmq9rinVNTTZKmqR/tymjryS0X8odKdiYi6gt+dSYZCY3Y1SXjNLCLpLyW9QdIxSU+Y2X3OuRembfZ+ScPOuc1mdqekP5H0HjO7QtKdkq6U1CPp+2a21TlXqUXsuHBVd3aVtXeipGLFH1ZcrDgVKn6ymi1VzySz2XJVubK/zXTfPOI3NzBJLTFPrTFPmUREq1tjGi9WzySyTw3MnsheseyliWw8wjeCAAAA9Wqq6/P25UkpEdFkqaqPXNOl07myTmbL6s+VNVSo6OBYUc8NnX3s+KnnBtUWj6gt5qktOHZsi3lKRz0lo/7ovmTUL4Z4VImXhFpVeG+UtN85d1CSzOweSbdJmp7w3ibpo8Hlr0n6C/Mzl9sk3eOcK0g6ZGb7g/t7tEaxh8Y5d+abLBf846YuSwpmxZGTO7N+6tfU9Yqcqk6qOr8x08zLFefknFRxfqJacf5E4WXnVKoGl6tOZXf25dKZy06lilSo+l2M/3TngGbkrPr7faMveW4mKe6Z4hH/PI1V6ajSwYdROmpKRT3de2hcv355h/8BFbGXDF2ZOWyYRBYAAKB57BzIa3Mw5ePU1EZRM61Kx7QqHZPkH5NmS1VNlqt6dqigjkRECc80Uqiod6J0pl/LbBKeXziJeqaoJ8U8U9SmXQ6uxzwF25ii5l/2TPJkmjo8fWGoIM8kM51Z97PrJk/BcjPZmcs/uy7525592V9n09cFv23aY09f/7Pbc1w8pVYJ76WSeqddPybppnNt45wrm9mopK5g+WMzbnvp4oW6+D67e1hD+cp5k9l6E/OkiFnw5lfwYeBfTsdMy7yohgv+3LNTHxRRMz0/XND13cmfbR98SETs/G/EfHmqOlyu1VMEAABAg4l5pkwiokwiomeHCrq6M+lXhgOlqj8TRy4YTZgPRhrmpzUw9Ys60nChooqrnikKVWYUis53nH7fkfHFf7IXaHqyPNu6WW8zbcU7N7RrY3t8ESKrLXNu8VMsM7tD0puccx8Irv+KpBudc/9+2ja7gm2OBdcPyK/kfkzSo865u4Pln5F0v3Pu6zMe44OSPhhc3SZpz+I+K0yzXNJA2EFgUbBvmxf7trmxf5sX+7Z5sW+bF/t2caxzznW/3Ea1qvAek7Rm2vXVkk6cY5tjZhaVlJE0NMfbyjn3N5L+ZgFjFMVFEQAADKFJREFUxhyZ2ZPOuevDjgMLj33bvNi3zY3927zYt82Lfdu82Lfh8l5+kwXxhKQtZrbBzOLym1DdN2Ob+yTdFVy+XdIPnV9+vk/SnWaWMLMNkrZIerxGcQMAAAAAGlRNKrzBObkflvSA/GmJPuuc22VmH5P0pHPuPkmfkfSloCnVkPykWMF2X5Hf4Kos6bfo0AwAAAAAeDk1m4fXOXe/pPtnLPuDaZfzku44x20/LunjixogLgZDyZsX+7Z5sW+bG/u3ebFvmxf7tnmxb0NUk6ZVAAAAAADUWq3O4QUAAAAAoKZIeDEvZnaHme0ys6qZXT9t+Xozy5nZzuDn02HGifk51/4N1v2+me03sz1m9qawYsTFM7OPmtnxae/Xt4YdEy6Omb05eG/uN7PfCzseLCwzO2xmzwXv1yfDjgfzZ2afNbPTZvb8tGWdZvY9M9sX/F4WZoyYn3PsW/7ehoiEF/P1vKR3SXpolnUHnHPbg58P1TguLIxZ96+ZXSG/odyVkt4s6a/MLFL78LCAPjnt/Xr/y2+OehW8F/9S0lskXSHpvcF7Fs3lF4P3K1OcNLbPy/87Ot3vSfqBc26LpB8E19F4Pq+X7luJv7ehIeHFvDjnXnTO7Qk7DiyO8+zf2yTd45wrOOcOSdov6cbaRgfgHG6UtN85d9A5V5R0j/z3LIA645x7SP6sJNPdJukLweUvSHpHTYPCgjjHvkWISHixGDaY2Q4z+5GZvTrsYLCgLpXUO+36sWAZGteHzezZYAgWw+caG+/P5uckfdfMnjKzD4YdDBbcJc65k5IU/F4RcjxYWPy9DQkJL87JzL5vZs/P8nO+isFJSWudc9dK+oikfzCz9tpEjAsxz/1rsyyj1Xsde5n9/NeSNknaLv+9+2ehBouLxfuz+d3snHul/GHrv2Vmt4QdEIA54e9tiGo2Dy8aj3Pu9fO4TUFSIbj8lJkdkLRVEs016sx89q/8itGaaddXSzqxMBFhMcx1P5vZ30r650UOB4uL92eTc86dCH6fNrNvyB/GPlsvDTSmPjNb5Zw7aWarJJ0OOyAsDOdc39Rl/t7WHhVeLCgz655qYmRmGyVtkXQw3KiwgO6TdKeZJcxsg/z9+3jIMWGeggOqKe+U36wMjesJSVvMbIOZxeU3mLsv5JiwQMysxczapi5LeqN4zzab+yTdFVy+S9L/CTEWLCD+3oaLCi/mxczeKelTkrolfcvMdjrn3iTpFkkfM7OypIqkDznnOHG/wZxr/zrndpnZVyS9IKks6becc5UwY8VF+VMz2y5/2OthSb8Rbji4GM65spl9WNIDkiKSPuuc2xVyWFg4l0j6hplJ/vHbPzjnvhNuSJgvM/tHSbdKWm5mxyT9oaRPSPqKmb1f0lFJd4QXIebrHPv2Vv7ehsec4/QeAAAAAEDzYUgzAAAAAKApkfACAAAAAJoSCS8AAAAAoCmR8AIAAAAAmhIJLwAAAACgKZHwAgAalpl91Mzunudtf9XMfnKe9d82s7tm29bMJoK5xjGDmX3azP7LtOu/aWZ9wWvWZWY3m9m+4Po7wowVAND8SHgBADVlZofNLBckPH1m9jkzaw07rpmcc29xzn3hHOtanXMHJcnMPm9mfzTfxwlej6KZLZ+xfKeZOTNbP9/7XmjT9t24mY2Y2SNm9iEzO3M84Zz7kHPuvwXbxyT9L0lvDF6zQUkfk/QXwfV7w3kmAIClgoQXABCGtzvnWiW9UtINkv6/mRuYb6n8nTok6b1TV8zsakmp8MI5r7c759okrZP0CUm/K+kz59j2EklJSbumLVs34/qcmVl0PrcDACxdS+VAAgBQh5xzxyV9W9JVkmRmD5rZx83sYUlZSRvNrMfM7jOzITPbb2a/PuNukmb25aDq+LSZXTO1wsx+z8wOBOteMLN3zritmdmnzGzUzHab2eumrXjQzD4wW9xB5XWzmX1Q0r+R9P8EFetvmtl/MrOvz9j+U2b25+d5Kb4k6X3Trt8l6Ysz7iNhZv/TzI4GlfFPm1kqWLfMzP7ZzPrNbDi4vHrGc/lvZvZw8Fp8d6qibGZJM7vbzAaDqu0TZnbJeWKVJDnnRp1z90l6j6S7zGxqH37ezP7IzLZK2hNsPmJmPzSzA5I2Svpm8HolzCxjZp8xs5Nmdjy4bSS4r18NYv6kmQ1J+miw/N+Z2YvBc33AzNbN2DcfCoZND5vZX5qZTVv/68Ftp/5PvDJY3mNmXw9ew0Nm9tsv9xoAAOofCS8AIDRmtkbSWyXtmLb4VyR9UFKbpCOS/lHSMUk9km6X9MfTE1NJt0n6qqROSf8g6d5gKK0kHZD0akkZSf9V0t1mtmrabW+SdFDSckl/KOmfzKxzrvE75/5G0t9L+tNgiO7bJd0t6c1m1hE8x6j8pPBL57mrxyS1m9nlQbL3nuB+pvsTSVslbZe0WdKlkv4gWOdJ+pz86ulaSTlJfzHj9r8s6dckrZAUl/Qfg+V3yX991kjqkvSh4PZz4px7XP7+efWM5XslXRlc7XDOvdY5t0nSUQUVfudcQdIXJJWD53StpDdKmv5Fw9Q+WiHp4+af9/v/SnqXpG5JP5b/f2S6t8kfOXCNpHdLepMkmdkd8pPm90lql/SvJA0GIwm+KekZ+a/r6yT9jpm9aa6vAwCgPpHwAgDCcK+ZjUj6iaQfSfrjaes+75zb5ZwrS1op6VWSftc5l3fO7ZT0d/KT4ilPOee+5pwryT9fNCnp5yTJOfdV59wJ51zVOfdlSfsk3Tjttqcl/blzrhSs3yPply7miTnnTkp6SNIdwaI3Sxpwzj31MjedqvK+QdJuScenVgQVyl+X9H8754acc+PyX7M7g8ccdM593TmXDdZ9XNJrZtz/55xze51zOUlfkZ84S1JJfqK72TlXcc495Zwbu8CnfUL+Fw4XJKgkv0XS7zjnJp1zpyV9cup5Td23c+5TzrlyEPtvSPrvzrkXg/8jfyxp+/Qqr6RPOOdGnHNHJf3LtOf6AflfTjzhfPudc0fkJ8fdzrmPOeeKwfnZfzsjDgBAA+JcGABAGN7hnPv+Odb1TrvcI2kqwZtyRNL1s23vnKua2VQ1WGb2PkkfkbQ+2KRVfjV3ynHnnJtx3z0X8DzO5QuSflN+0vRvdf7q7pQvyU+UN2jGcGb5lcy0pKemj86VNDX0Ny0/UXyzpGXB+jYzizjnKsH1U9PuLyv/tZh63DWS7gmq0ndL+s/BFwhzdamkoQvYfso6STFJJ6c9L09n/x/oneU2/9vM/mzaMgtiOBJcP9dzXSO/6j9bHD3BlzBTIvKrxwCABkaFFwBQb6YnoCckdZpZ27RlazWt+ik/iZEkBUNTV0s6EVT8/lbShyV1Oec6JD0vPzmacun08zuD+z5xEfFOuVfSK4LzWt8mf9jz+e/ErzQekj/E+59mrB6QP8z4SudcR/CTCRp/SdJ/kLRN0k3OuXZJtwTLTS8jqG7/V+fcFZJ+IYj3fS9zszPM7Ab5yeY5p3g6j15JBUnLpz2vdufcldO2mfn69kr6jWnbdzjnUs65R+b4eJvOsfzQjPtsc869dR7PCQBQR0h4AQB1yznXK+kRSf89aK70Cknv19kJ5HVm9q7gXNnfkZ9APSapRX6y1C9JZvZrCppjTbNC0m+bWSw4v/NySfdfYJh98hsxTY87L+lr8s8pfjwYWjsX75f0Wufc5Iz7q8pP3j9pZiuC53PptHNM2+QnxCPBOch/ONfgzewXzezq4NzhMflDnCsvczOZWbuZvU3SPZLuds49N9fHnBIM//6upD8L7s8zs01mNnM49nSflvT7ZnZlEEcm2Hdz8XeS/qOZXWe+zcEXI49LGjOz3zWzlJlFzOyqIJkHADQwEl4AQL17r/whySckfUPSHzrnvjdt/f+R3+RpWP65ve8KqpYvSPozSY/KT0qvlvTwjPv+qaQt8iuoH5d0ezBX7IX4jKQrgg7H0+eV/ULwmHMZzixJcs4dcM49eY7Vvytpv6THzGxM0vflV3Ul6c/lT2M0ID/Z/84FxL9SfnI+JulF+edUz2yYNd03zWxcflX0P8s/b/rXLuDxZnqf/CZaL8jfh1+TtOpcGzvnviG/gdc9wevwvPzzgF+Wc+6r8vfzP0gal1+J7wyGfb9d/rm+h+S/jn8nv5kXAKCB2dmnLgEAgIVgZmvlN59aOY8mUAAAYAFQ4QUAYIEF5xJ/RNI9JLsAAISHLs0AACwgM2uRP4T6iPyuyQAAICQMaQYAAAAANCWGNAMAAAAAmhIJLwAAAACgKZHwAgAAAACaEgkvAAAAAKApkfACAAAAAJoSCS8AAAAAoCn9/5Z6DX85WPU5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Interpret visually\n",
    "\n",
    "simulated_means = np.asarray(simulated_means)\n",
    "plt.figure(figsize=(16, 8))\n",
    "sns.distplot(simulated_means, color=\"skyblue\", kde=True, hist=True, rug=\"True\")\n",
    "plt.title(\"Simulated Differences of drug and placebo group means for the null hypothesis\", fontsize=12)\n",
    "plt.xlabel(\"Probability Means Difference\", fontsize=12)\n",
    "plt.ylabel(\"Percent\", fontsize=12)\n",
    "plt.axvline(observed_diff, color='r'); # observed mean difference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Permutation test result\n",
    "It is very unlikely that the observed_diff value or greater appeared if there is no difference in the groups, so **there is a difference between the drug and the placebo groups**. We have sufficient evidence to reject the null hypothesis, and conclude that the drug is effective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Required number of permutations for a permutation-based p-value](https://stats.stackexchange.com/questions/80025/required-number-of-permutations-for-a-permutation-based-p-value)\n",
    "- [What is a Permutation Test](https://www.r-bloggers.com/what-is-a-permutation-test/amp/)\n",
    "- [Permutation Tests, R example](https://thomasleeper.com/Rcourse/Tutorials/permutationtests.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "def tstatistic(sample1, sample2, alpha, one_tail):\n",
    "    '''\n",
    "    Calculate t-statistic for 2 samples\n",
    "    Are the two samples expected to have been drawn from the same population?\n",
    "    '''\n",
    "    # Calculate means\n",
    "    mean1, mean2 = np.mean(sample1), np.mean(sample2)\n",
    "    \n",
    "    # Calculate sample standard deviations\n",
    "    std1, std2 = np.std(sample1, ddof=1), np.std(sample2, ddof=1)\n",
    "    \n",
    "    # Calculate standard errors\n",
    "    n1, n2 = len(sample1), len(sample2)\n",
    "    se1, se2 = std1/np.sqrt(n1), std2/np.sqrt(n2)\n",
    "\n",
    "    # Standard error on the difference between the samples\n",
    "    sed = np.sqrt(se1**2.0 + se2**2.0)\n",
    "    \n",
    "    # Calculate degrees of freedom\n",
    "    df = len(sample1) + len(sample2) - 2\n",
    "    \n",
    "    # Calculate the t statistic\n",
    "    t_stat = (mean1 - mean2) / sed\n",
    "    \n",
    "    # Calculate the critical value\n",
    "    if one_tail:\n",
    "        cv = stats.t.ppf(1.0 - alpha, df)\n",
    "    else:\n",
    "        cv = stats.t.ppf(1.0 - alpha/2, df) # 2-tailed\n",
    "    \n",
    "    # Calculate the p-value\n",
    "    p = (1.0 - stats.t.cdf(abs(t_stat), df)) * 2.0\n",
    "    \n",
    "    return t_stat, df, cv, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t, df, cv, p = tstatistic(drug, placebo, alpha, one_tail=False)\n",
    "t, df, cv, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}