{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Lab 6: Comparing Three or More Groups with ANOVA\n",
    "## Created by Jane and Kristin, modified by Xinyue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Often, we have to compare data from three or more groups to see if one or more of them is different from the others. To do this, scientists use a statistic called the _F_- statistic, defined as $F=\\frac{Var_{between}}{Var_{within}}$. The variability between groups is the sum of the squared differences between the group means and the grand mean. The variability within groups is the sum of the group variances. In LS 40, we will use a variation of the F-statistic that does not require squaring, which is called the _F_-like statistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "1. Import pandas, Numpy, Seaborn and Pyplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "In this lab, we will examine the results of a pharmaceutical company's study comparing the effectiveness of different pain relief medications on migraine headaches. For the experiment, 27 volunteers were selected and 9 were randomly assigned to one of three\n",
    "drug formulations. The subjects were instructed to take the drug during their next migraine headache episode and to report their pain on a scale of 1 = no pain to 10 = extreme pain 30 minutes after taking the drug."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "2. Using the pandas `read_csv` function, import the file `migraines.csv` and show the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Drug A</th>\n",
       "      <th>Drug B</th>\n",
       "      <th>Drug C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4</td>\n",
       "      <td>7</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>8</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4</td>\n",
       "      <td>6</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>3</td>\n",
       "      <td>8</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>3</td>\n",
       "      <td>7</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Drug A  Drug B  Drug C\n",
       "0       4       7       6\n",
       "1       5       8       7\n",
       "2       4       4       6\n",
       "3       2       5       6\n",
       "4       2       4       7\n",
       "5       4       6       5\n",
       "6       3       5       6\n",
       "7       3       8       5\n",
       "8       3       7       5"
      ]
     },
     "execution_count": 8,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "migraines = pd.read_csv(\"migraines.csv\")\n",
    "migraines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "3. Visualize the data. Use as many different plots as you need to get a sense of the distributions. What are your initial impressions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7fe75c14aca0>"
      ]
     },
     "execution_count": 11,
     "metadata": {
     },
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAALACAYAAABhH7gMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAABYlAAAWJQFJUiTwAAA7s0lEQVR4nO3deZhuV1kn7N8TIlOABFGIzWCYQiL4ASEyGIRANB1BZFBEBYS0oHyAIA2NNghJaBFQlHmylUSitgN8iANDBBIghKElICKQMCUEIQyBEAhz8nx/7F1yqFSdU3XOOlXnrbrv66pr1bvX3ms/75u6cupXa++1q7sDAADAGPttdgEAAABbiZAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJC1l1TVX1TVX2x2HQAAwMbaf7ML2MIOO+KII45I8subXQgAwIKqzS4AdoeZLAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIG2RMiqqgdVVc9fD1vnsT9SVX9TVZ+rqm9U1TlVdVJVXW1v1QsAAGxdCx+yquqGSV6Y5Ku7cewdkvzfJPdJ8sYkz0tySZKnJvnnqrrKuEoBAIDtYKFDVlVVkpOTXJTkpes89krzsVdP8vPd/cvd/VtJ7pDkVUmOSvK4sRUDAABb3UKHrCSPSXL3JMcnuXSdx941yeFJ3trdf7+0sbsvT/LE+eUj5iAHAACwJvtvdgG7q6oOT/LMJM/r7rdW1d3XOcTS/q9f3tHdH6+qc5McmuQmST62kzres0rXYeusB4AtyN/q2Iq6e7NLgH3aQs5kVdX+SU5N8skkT9rNYW4xt+eu0v+RuT10N8cHAAC2oUWdyXpqktsmuXN3f303xzhwbr+8Sv/S9oN2Nkh3326l7fMM1xG7VRkAW87xrz9+s0vY0k4+7uQkSZ98z02uZGur4/9ps0uAhbBwM1nzioBPSvKH3f2Oza4HAABgRwsVsubLBF+R6RK/p+zhcEszVQeu0r+0/eI9PA8AALCNLFTISnKNTPdIHZ7kGzs8gLiTnDDv87/nbc/dxVjnzO1q91zdfG5Xu2cLAADgChbtnqxvJvnTVfqOyHSf1pmZAtSuLiV8c5InJzkuyTN27Kiqm2QKX+cn+fge1AsAAGwzCxWy5kUuHrZSX1WdmClk/Vl3/8kO26+e5EZJvtbdn9zhkLck+VCSu1TVzy49K6uq9kvyrHmfl7Y1SgEAgHVYqJC1m26f5PRMoeropY3dfVlVHZ9pRuuVVfXKTEvCH5PkyCRvT/KcDa8WAABYaIt2T9ZQ3f2uJD+W5DVJjk3yuEwLXjwtyU919zc3sTwAAGABbZmZrO4+McmJK2w/I0nt5LgPJrn/3qoLAADYXrb1TBYAAMBoQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBACxeyqupZVfWmqrqgqr5eVV+sqvdW1QlVdZ11jHNeVfUqXxfuzfcAAABsXftvdgG74XFJzk7yz0k+l+SAJHdMcmKSX6uqO3b3BWsc68tJnrvC9q/ueZkAAMB2tIgh61rd/Y3lG6vq6UmelOR/JnnkGse6uLtPHFgbAACwzS3c5YIrBazZ38ztzTeqFgAAgOUWcSZrNfea2/ev45irVNWDktwoyaXzsW/t7stGFwcAAGwPCxuyquoJSa6R5MAkRya5c6aQ9Mx1DHNwklOXbftEVR3f3W9ZYx3vWaXrsHXUAQAAbBELG7KSPCHJ9XZ4/fokD+3uz6/x+JOTvC3Jvyf5SpKbJHl0kl9L8rqqulN3/+vAegEAgG1gYUNWdx+cJFV1vSQ/nmkG671V9TPdffYajj9p2aYPJHlEVX01yeMzrVZ43zWMc7uVts8zXEfs6ngAAGBrWbiFL5br7s9296uTHJvkOklesYdDvnRu77KH4wAAANvQwoesJd19fpIPJrllVf3AHgy1dLnhAXteFQAAsN1smZA1+y9zuyerA95xbj++h7UAAADb0EKFrKo6tKoOXGH7fvPDiK+b5Kzu/tK8/fuq6rCquumy/Q+vqivMVFXVIUleOL/88+FvAAAA2PIWbeGLeyR5RlWdmeQTSS7KtMLgXTOtDnhhkofvsP/1k3woyflJDtlh+wOSPL6q3jr3fSXJTZPcM8lVk7w2ybP35hsBAAC2pkULWW9McrNMz8S6bZKDMj1E+NxMz7t6fnd/cQ3jnJ7kFvMYR2W6/+riJGfO45za3T24dgAAYBtYqJDV3R/I9Cyrte5/XpJaYftbkqzpYcMAAADrsVD3ZAEAAOzrhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBhCwAAICBFi5kVdWzqupNVXVBVX29qr5YVe+tqhOq6jrrHOsGVfXyqvp0VX2zqs6rqudW1bX3Vv0AAMDWtnAhK8njkhyQ5J+TPC/JXyT5TpITk7y/qm64lkGq6qZJ3pPk+CTvTvKcJB9P8tgk71hvYAMAAEiS/Te7gN1wre7+xvKNVfX0JE9K8j+TPHIN47w4yXWTPKa7X7DDOH+UKcg9PckjhlQMAABsGws3k7VSwJr9zdzefFdjzLNYxyY5L8mLlnWfkOTSJA+uqgN2s0wAAGCbWsSZrNXca27fv4Z97za3p3X35Tt2dPdXqurtmULYHZO8aWcDVdV7Vuk6bA11sECqarNLAABgASxsyKqqJyS5RpIDkxyZ5M6ZAtYz13D4Leb23FX6P5IpZB2aXYQsAACAHS1syEryhCTX2+H165M8tLs/v4ZjD5zbL6/Sv7T9oF0N1N23W2n7PMN1xBpqYcE84GVnbXYJW9pf//qPJ/E5721Ln/Pxrz9+kyvZ+k4+7uTNLgGADbZw92Qt6e6Du7uSHJzkfklukuS9VSXYAAAAm2ZhQ9aS7v5sd7860+V910nyijUctjRTdeAq/UvbL96z6gAAgO1m4UPWku4+P8kHk9yyqn5gF7ufM7eHrtK/tELhavdsAQAArGjLhKzZf5nby3ax3+lze2xVfc9nUFXXTHJUkq8leefY8gAAgK1uoUJWVR1aVVe4xK+q9psfRnzdJGd195fm7d9XVYfNz8X6T939sSSnJTkkyaOWDXdSkgOSnNrdl+6FtwEAAGxhi7a64D2SPKOqzkzyiSQXZVph8K6ZFr64MMnDd9j/+kk+lOT8TIFqR49MclaS51fVMfN+d8j0DK1zkzx5r70LAABgy1q0kPXGJDfL9Eys22ZaYv3STKHo1CTP7+4vrmWg7v5YVR2Z5GlJjssU4D6T5HlJTlqaDQMAAFiPhQpZ3f2BJI9ex/7nJamd9F+QxENiAACAYRbqniwAAIB9nZAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAw0EKFrKq6TlU9rKpeXVUfraqvV9WXq+rMqvrVqlrz+6mq86qqV/m6cG++DwAAYOvaf7MLWKf7J3lJks8kOT3JJ5NcL8n9kvxJkp+uqvt3d69xvC8nee4K27+656UCAADb0aKFrHOT/GySf+ruy5c2VtWTkrw7yc9lClyvWuN4F3f3iaOLBAAAtq+Fulywu9/c3f+wY8Cat1+Y5KXzy6M3vDAAAIDZos1k7cy35/Y76zjmKlX1oCQ3SnJpkvcneWt3Xza6OAAAYHvYEiGrqvZP8ivzy9ev49CDk5y6bNsnqur47n7LGs/9nlW6DltHHQAAwBaxUJcL7sQzk9wqyWu7+w1rPObkJMdkCloHJPnRJC9LckiS11XVrfdCnQAAwBa38DNZVfWYJI9P8uEkD17rcd190rJNH0jyiKr66jzeiUnuu4ZxbrdKXe9JcsRa6wEAALaGhZ7JqqpHJ3lekg8muVt3f3HAsEsLaNxlwFgAAMA2s7Ahq6p+M8kLMs1A3W1eYXCEz8/tAYPGAwAAtpGFDFlV9VtJnpPkfZkC1ucGDn/Huf34wDEBAIBtYuFCVlU9JdNCF+9Jckx3f2En+35fVR1WVTddtv3wqrrCTFVVHZLkhfPLPx9XNQAAsF0s1MIXVfWQJE9LclmStyV5TFUt3+287j5l/v76ST6U5PxMqwYueUCSx1fVW+e+ryS5aZJ7JrlqktcmefZeeRMAAMCWtlAhK8mN5/ZKSX5zlX3ekuSUXYxzepJbJLltkqMy3X91cZIzMz0369Tu7j0rFQAA2I4WKmR194mZllZf6/7nJbnCVNf8oOE1PWwYAABgPRbuniwAAIB9mZAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAwkJAFAAAw0NCQVVU3qqpr7WKfa1bVjUaeFwAAYF8xeibrE0keu4t9HjPvBwAAsOWMDlk1fwEAAGxLm3FP1sFJLt2E8wIAAOx1++/pAFX1K8s23WaFbUlypSQ3SvKgJP+2p+cFAIC9pap6s2tYi+52Fdk+aI9DVpJTkiz9EHaSe89fyy39AHwtyUkDzgsAALDPGRGyjp/bSvLyJH+X5DUr7HdZkouSvKO7Lx5wXgAA2Kt++Lf+8YzNrmEl5z/rZ44eNdYKs3bfSnJJkguSnJ3kVUlO6+7LRp1zs1TVuUlunimT/PjeOs8eh6zu/rOl76vqIUn+rrtfsafjAgAAG2rparMrJTkoyS2TPDjJryb5l6p6YHefu0m17bGqulumgNVJ7lRVt+ruD+yNc42YyfpP3X23keMBAAAbo7tPXL6tqq6X5AVJ7p/kjVV1ZHd/bqNrG+TX5vZZSX57fv2YvXGizVhdEAAAWADd/dkkv5jkjCQ3TPKkHfur6oyq6qq6clU9tarOqapvVtUpc/+Jc//Ry8euqkPmvlNW6Du0ql5VVV+qqkur6qyqumdVPXQ+5qHreR9VdZ0k903ykSRPSXJhkgdV1VXXM85aDQ9ZVXXXqvrHqvpcVX27qi5b4es7o88LAACM192XJ/nd+eUvVdVKKxq+Kskjk5yV5LnZg9XEq+qwJO9Mcr8kb0/yvCSfTPLqJPfZzWEfkuQqSU7p7u8k+Ysk1840Qzfc0MsFq+qemRa+uFKmD+KcJAIVAAAstjMz/V5/3SSHJPnEsv4fTnKr7v7CgHO9KFMAemR3v2RpY1X9dJLX7uaYD09yeZKltSNOSfL4TJcMnrrbla5iaMhKcmKSbye5Z3efNnhsAABgE3T3N6vqoiTXS/KDuWLIesqIgFVVN0xy9yQfTfKyZTW8rqremOQn1znmTyQ5LNMKiZ+ax/pAVb0nyZ2r6vDu/tCe1r6j0ZcL3irJXwtYAACw5SxdJrjSg5rfPegct5nbd8yXKS535m6MubTgxcnLtp8ytw/fjTF3anTI+mqSLw4eEwAA2ETzAhHfP7/8/Aq7XDjoVAfO7WdX6V9t+4qq6tpJfj7JxZlua9rRX2Z6JtivVNVV1jPurowOWW9KcqfBYwIAAJvrzpluNfpsd5+3vLO7V5rdSqb7oJKVb1M6aIVtl8zt9VYZb7Xtq/mVJFedz/X1eWXCnh/AfFGSKye5TpKfW+e4OzX6nqzfSvLuqvqdJE/fyYcNAAAsgKraL8mT55d/uc7DvzS3N1yh78gVtr1vbu9UVfutcMngndd5/qVLAf9Pkq+t0H9gppmuh2f9721Vo0PWCUn+PdPTov9bVb0v09Tcct3dvzr43AAAwEBVdd0kL0xydKbVw39vnUMs3at1fFWdOi+fvrTAxVOX79zdn6yqM+bz/XqSHVcXPC7rWPSiqn48yS2TfLC7f3mVffZL8vEkR1fVzbv7I2sdf2dGh6yH7vD9IfPXSjqJkAUAwD7t/Gf9zNGbXcNGqaoT52/3y3R53S0zzRxdOVNYeuB6VxDs7ndV1VuT3CXTFW9vznTJ372SvCErz3A9KtPzsV5cVfdI8v4kN8l0Sd9rktw7370McWeWFrz4053Ud3lVnZxplfRfS/I/1jDuLo0OWTcePB4AALAxTpjbbyX5SpLzMz1X6lWZlj9fS7BZyb2T/MHc/kaSjyR5YpLTkvzC8p27+4NVdadMs2Z3n7/en+S+SQ6fx7lk+XE7qqoDMz1o+Fv57rOxVvPyTLNqD6mqJ3f3t9b8zlYxNGR19/kjxwMAgM3Q3bXrvbaGPXmv3X30Gva5ONM9Tystlb7iubv7w0nud4Wdq5Yu+9vpc626+8tJDthVbfO+FyS50lr2XavRqwsCAADstqrar6oOXmH7MUkekOkeq3M2vrK1GzqTVVU3Wuu+3f3JkecGAAC2hCsnuaCqTk/y4STfyXR/2E9luvzvUZtY25qMvifrvKz8BOjlei+cGwAAWHzfTvLSTPdi3SHJ1ZN8IcnfJnlmd793E2tbk9FB5xVZOWQdlOQ2SX44yRmZbqIDAAD4Ht19WaYFMhbW6IUvHrpa37wG/VOSPCLJQ0aeFwAAYF+xYQtfdPfl3X1SpksKn7lR5wUAANhIm7G64FlJjt2E8wIAAOx1mxGyvj9rXLMeAABg0WxoyKqqn8y0tv0HNvK8AAAAG2X0c7LevJPz3DDJ0nO0njbyvAAAAPuK0Uu4H73K9k7ypSRvSPLs7l4tjAEAACy00Uu4b8Y9XgAAAPuM0TNZAACw8KqqN7uGteju2uwauKK9GrKq6ppJDkry5e6+ZG+eCwAAYF8wPGRV1f5JnpDkYUluvMP2TyT5k0z3ZH1n9HkBAGC0PuFaZ2x2DSupky45ethYV5y1+1aSS5JckOTsJK9Kclp3XzbqnBulqh6a5OQVur6a5CNJXpnkud39tZHnHb264JWTvD7JXTMtdnFBks8k+aEkhyR5epLjqurY7v7Wbox/nST3TXLPJD+a5PqZfgj+LdOHd3J3X76O8W6QaaXD45JcZ67175Kc1N1fWm99AACwwE6a2ytluhrtlkkenORXk/xLVT2wu8/dpNr21L9m+j0/mR5jdXCSe+W7+eRuI0Pk6Jms/55phcF/TPL47v7IUkdV3TTJH2Z6M/89yTN3Y/z7J3lJpjB0epJPJrlekvtlmiX76aq6f3fv8hrauZ6zklw3yWuSfDjJ7ZM8NtMHfVR3X7QbNQIAwMLp7hOXb6uq6yV5Qabfw99YVUd29+c2urYB3rf8/VXVQUnen+Qn5q8zRp1s9GqAv5zpQcP32TFgJUl3fyxTGPr3JA/czfHPTfKzSW7Q3Q/s7v/Z3f8tyWGZZs1+bj7HWrw4U8B6THffp7t/u7vvnuQ5SW6RKdUCAMC21d2fTfKLmQLIDZM8acf+qjqjqrqqrlxVT62qc6rqm1V1ytx/4tx/9PKxq+qQue+UFfoOrapXVdWXqurSqjqrqu5ZVQ+dj3nogPd2cZL/O7/8wT0db0ejQ9bNkrxutUv25u2vS3LT3Rm8u9/c3f+wfPzuvjDJS+eXR+9qnHkW69gk5yV50bLuE5JcmuTBVXXA7tQJAABbxfy79+/OL3+pqlZa0fBVSR6Z6Uqx52a6nWe3VNVhSd6ZafLk7Umel+kKtlcnuc/ujrvCeQ5M8mNJLk/y3lHjJuMvF/xWkmvsYp8Dknx78Hmzw5hrWVTjbnN72gqB7StV9fZMIeyOSd40rkQAAFhIZ2b6Pfu6mdZa+MSy/h9Ocqvu/sKAc70oybWTPLK7X7K0sap+Oslrd3PM21TVifP3+2W65ehnkhyY6cq2j+5+uVc0OmS9P8nPV9WJ3f355Z1V9QNJfj7TjWfDzCsa/sr88vVrOOQWc7vajXsfyRSyDs0uQlZVvWeVrsPWUMcQK/8xARbbGaefsdklbAsnH7fSgkuwuE4//fTNLgG2pO7+ZlVdlCmc/GCuGLKeMiJgVdUNk9w9yUeTvGxZDa+rqjcm+cndGPrW89dy/yfJm3djvJ0afbngCzN96O+uql+tqptU1dWq6sZVdXySd839Lxx83mcmuVWS13b3G9aw/4Fz++VV+pe2H7SHdQEAwFax9Jf9lRaZe/egc9xmbt+xyi1IZ+7muH/W3bX0lWl1wQdlmlh5V1UdsZvjrmjoTFZ3/01V3SbJbyf54xV2qSS/391/M+qcVfWYJI/PtDrgg0eNu1bdfbuVts8zXEP/Y+3KA1521kaebtv561//8c0uYVs5+vDrbHYJW9pfz+0H/3C3bpFlHX7k8R/b7BK2lYNv/UObXcIW5+d5u6qqqyb5/vnlFa5YS3LhoFMtTYZ8dpX+1bavy7ygx19U1dWS/O8kz0jyX0eMneyFhxF395Oq6u8zrad/20wf1Jcz3Uz28u5+x6hzVdWjM90I98Ekx3T3F9d46NJM1YGr9C9tv3j3qwMAgC3jzpmyw2e7+7zlnTt5hNLSbNRKueOgFbZdMrfXW2W81bbvrnfN7e1HDjo8ZCVJd78z04oge01V/Wam5dY/kClgrWe9/nPm9tBV+m8+t4v6sDUAABiiqvZL8uT55V+u8/Avze0NV+g7coVt75vbO1XVfitcMnjndZ5/V649t0Nvo9rjweY18d9dVW+qqu/bxX5vqqp37my/NZ7ztzIFrPcludtuPBBt6a7YY+cfmh3HvmaSo5J8LXs5KAIAwL6sqq6b5K8yPSbpk0l+b51DLN2rdfy8WN3SuDdM8tTlO3f3JzM9k+tmSX59WS3HZfcWvVhRVV0pyWPnl2eMGjcZM5P1oCS3S3Kv7l51afbu/lZV/UGmZRcfmOSU3TlZVT0lydOSvCfJsTu7RHAOczdN8u35YchLtXysqk7LdKPbozI9xXrJSZmWmX9Zd1+6OzUCALA11EmXHL3ZNWyUZUucH5Tklplmjq6cKSw9cL0rCHb3u6rqrUnukmlxvDdnuuTvXknekJVnuB6V6flYL66qe2RawfwmSX4uyWuS3DvfvQxxrXZcwj2ZlqK/e6ZVx7+Q5InrHG+nRoSs+yX5eHfvcs367n59VX0kyf2zGyGrqh6SKWBdluRtSR6zwvLl53X30tjXT/KhJOdnWs9/R0sPS3t+VR0z73eHTM/QOjffnRIFAIDt4IS5/VaSr2T6HfoVmR40fIXny67DvZP8wdz+RqbHJT0xyWlJfmH5zt39waq6U6ZZs7vPX+9Pct8kh8/jXLL8uF1YvoT7N5Kcl2l9h9/v7k+vc7ydGhGybpv1PRTsrUnusZvnuvHcXinJb66yz1uyhgA3z2YdmSm0HTfX9JlMH/RJ3f2lnR0PAMDWNS/zvS3syXvt7qPXsM/FSR4+fy234rm7+8OZJnO+d+eqX56//dAa6zslu3kF3Z4YEbJ+IOtbSvGzSXZrbebuPjHJievY/7ys8h9u7r8gyfG7UwsAADDevGbCdbv7wmXbj0nygCQf7O5zVjx4HzEiZH09yTXWsf81Mk3PAQAALHflJBdU1emZnoX7nUz3h/1UpksZH7WJta3JiJB1QVZefnE1R2ZamQQAAGC5byd5aaZ7se6Q5OqZFqf42yTP7O73bmJtazIiZJ2R5JFVdWR3/8vOdqyq2yX58Xzvan4AAABJku6+LNMCGQtrxEO3Xpikk/xtVR2+2k5VdVim9HlZkhcPOC8AAMA+Z49nsrr7nKp6WqYFKd5bVa9M8uYkn5p3uX6SYzKta3+VJE/d129UAwAA2F0jLhdMdz+tqr6TaW39X07yS8t2qUzXVj65u58x4pwAAAD7oiEhK0m6+/eq6i+S/LckRyX5obnrM0nOTHJyd58/6nwAAAD7omEhK0nmEHXCLncEAADYokYsfAEAAMBMyAIAABhIyAIAABho6D1ZAACwFVRVb3YNa9Hdtdk1cEVmsgAAAAYykwUAAKu41Sm3OmOza1jJBx76gaNHjbXCrN23klyS5IIkZyd5VZLTuvuyUefcaFW1X5L7ZXqe7+2T/GCSy5J8MsnbkvxZd7991PmELAAAIElOmtsrJTkoyS2TPDjJryb5l6p6YHefu0m17baqOjjJKzM9y/crSf45yceSVJKbZwpeD6+q3+juF444p5AFAACku09cvq2qrpfkBUnun+SNVXVkd39uo2vbXVV19SSvT3LrJH+V5JHd/aVl+1wryROSHDjqvO7JAgAAVtTdn03yi0nOSHLDJE/asb+qzqiqrqorV9VTq+qcqvpmVZ0y95849x+9fOyqOmTuO2WFvkOr6lVV9aWqurSqzqqqe1bVQ+djHrrGt/C4TAHr7UkeuDxgze/xku5+apJnr3HMXTKTBQAArKq7L6+q301ydJJfqqrHdffy+7heleTHkrwuyd8l2e3Zrqo6LMlZSa6d5J+SvD/JTZK8Oslr1zncr83t/+ruy3e2Y3d/c51jr0rIAgAAduXMJN9Jct0khyT5xLL+H05yq+7+woBzvShTwHpkd79kaWNV/XTWEbKq6oZJbpSp7rcMqGvNXC4IAADs1DzLc9H88gdX2OUpIwLWHIzunuSjSV62rIbXJXnjOob7obm9qLu/sae1rYeQBQAArMXSg49XelDzuwed4zZz+45VLu87c9B59iohCwAA2KmqumqS759ffn6FXS4cdKqlFf4+u0r/attX8pm5vc5c/4YRsgAAgF25c6b1HD7b3ect71xhIYwlS7NRK60FcdAK2y6Z2+utMt5q26+guy/I9LDh/ZPcZa3HjSBkAQAAq6qq/ZI8eX75l+s8fGnJ9Buu0HfkCtveN7d3ms+73J3Xef4/ntvfWWW8/1RVV1nn2KsSsgAAgBVV1XUzPcT36EyzQr+3ziGW7tU6vqr+czZrXuDiqct37u5PZnom182S/PqyWo5L8pPrPP9zkvxrkp9I8oqqOmj5DlV1jao6IdMDiYewhDsAAKziAw/9wNGbXcNGqaoT52/3y3Qp3y0zzRxdOVNYeuB6VxDs7ndV1VszXa737qp6c6ZL/u6V5A1ZeYbrUZkeHvziqrpHvvucrJ9L8pok9853L0Pc1fm/NoezVyZ5YJJ7VdU/J/lYpoU8bpbkmCTXSvLo9by3nRGyAACAJDlhbr+V5CtJzk/yikwPGj5tVw/z3Yl7J/mDuf2NJB9J8sQkpyX5heU7d/cHq+pOmWbN7j5/vT/JfZMcPo9zyfLjVtPdF1bVXTKFtF9KcsckP5MpqH0yyd8meXl3n7Wb7+8KhCwAAFimu2vXe20Ne/Jeu/voNexzcZKHz1/LrXju7v5wkvtdYeeqX56//dCai5zGuzxTmPrb9Ry3u9yTBQAA7DOqar+qOniF7cckeUCSD3b3ORtf2dqZyQIAAPYlV05yQVWdnuTDSb6T6f6wn8p0KeOjNrG2NRGyAACAfcm3k7w0071Yd0hy9SRfyHSp3zO7+72bWNuaCFkAAMA+o7svy7RAxsJyTxYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBACxeyqurnq+oFVfW2qrqkqrqq/nw3xjlvPnalrwv3Ru0AAMDWt/9mF7AbfifJrZN8Ncmnkhy2B2N9OclzV9j+1T0YEwAA2MYWMWQ9LlO4+miSuyY5fQ/Guri7TxxRFAAAQLKAIau7/zNUVdVmlgIAAHAFCxeyBrtKVT0oyY2SXJrk/Une2t2XbW5ZAADAotruIevgJKcu2/aJqjq+u9+ylgGq6j2rdO3JvWIAAMCCWrjVBQc6OckxmYLWAUl+NMnLkhyS5HVVdevNKw0AAFhU23Ymq7tPWrbpA0keUVVfTfL4JCcmue8axrndStvnGa4j9rBMAABgwWznmazVvHRu77KpVQAAAAtJyLqiz8/tAZtaBQAAsJCErCu649x+fFOrAAAAFtKWDllV9X1VdVhV3XTZ9sOr6gozVVV1SJIXzi//fANKBAAAtpiFW/iiqu6T5D7zy4Pn9k5Vdcr8/Re6+wnz99dP8qEk52daNXDJA5I8vqreOvd9JclNk9wzyVWTvDbJs/fKGwAAALa0hQtZSW6T5CHLtt1k/kqm0PSE7NzpSW6R5LZJjsp0/9XFSc7M9NysU7u7x5QLAABsJwsXsrr7xEzLq69l3/OS1Arb35JkTQ8bBgAAWI8tfU8WAADARhOyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABhKyAAAABlq4kFVVP19VL6iqt1XVJVXVVfXnuznWDarq5VX16ar6ZlWdV1XPraprj64bAADYHvbf7AJ2w+8kuXWSryb5VJLDdmeQqrppkrOSXDfJa5J8OMntkzw2yXFVdVR3XzSkYgAAYNtYuJmsJI9LcmiSayX5f/dgnBdnCliP6e77dPdvd/fdkzwnyS2SPH2PKwUAALadhQtZ3X16d3+ku3t3x5hnsY5Ncl6SFy3rPiHJpUkeXFUH7HahAADAtrRwIWuQu83tad19+Y4d3f2VJG9PcvUkd9zowgAAgMW2iPdkjXCLuT13lf6PZJrpOjTJm3Y2UFW9Z5Wu3bpXDGAj/cjjP7bZJWwbJx938maXAMAG2a4zWQfO7ZdX6V/aftDeLwUAANhKtutM1jDdfbuVts8zXEdscDkA63Kv33nOZpew5f3D7z4uSfLBP7zpJleytZmVBfYl23Uma2mm6sBV+pe2X7z3SwEAALaS7RqyzpnbQ1fpv/ncrnbPFgAAwIq2a8g6fW6Prarv+Qyq6ppJjkrytSTv3OjCAACAxbalQ1ZVfV9VHTY/F+s/dffHkpyW5JAkj1p22ElJDkhyandfuiGFAgAAW8bCLXxRVfdJcp/55cFze6eqOmX+/gvd/YT5++sn+VCS8zMFqh09MslZSZ5fVcfM+90h0zO0zk3y5PHVAwAAW93Chawkt0nykGXbbjJ/JVOgekJ2obs/VlVHJnlakuOS3CPJZ5I8L8lJ3f2lUQUDAADbx8KFrO4+McmJa9z3vCS1k/4Lkhw/oi4AAIBki9+TBQAAsNGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIGELAAAgIEWMmRV1Q2q6uVV9emq+mZVnVdVz62qa69jjDOqqnfyddW9+R4AAICtaf/NLmC9quqmSc5Kct0kr0ny4SS3T/LYJMdV1VHdfdE6hjxple3f2aNCAQCAbWnhQlaSF2cKWI/p7hcsbayqP0ryuCRPT/KItQ7W3SeOLhAAANi+FupywXkW69gk5yV50bLuE5JcmuTBVXXABpcGAACQZPFmsu42t6d19+U7dnT3V6rq7ZlC2B2TvGktA1bVA5LcOMm3knwoyZu7+5vjSgYAALaTRQtZt5jbc1fp/0imkHVo1hiykvzVstefq6pHdfcr13JwVb1nla7D1nh+AABgC1moywWTHDi3X16lf2n7QWsY6zVJ7pXkBkmulikUPWM+9q+r6rjdrhIAANi2Fm0ma5jufs6yTeckeVJVfTrJCzIFrtevYZzbrbR9nuE6Yk/rBAAAFsuizWQtzVQduEr/0vaL9+Acf5Jp+fbbVNU192AcAABgG1q0kHXO3B66Sv/N53a1e7Z2qbu/keQr80urFAIAAOuyaCHr9Lk9tqq+p/Z51umoJF9L8s7dPUFV3SLJtTMFrS/s7jgAAMD2tFAhq7s/luS0JIckedSy7pMyzTyd2t2XLm2sqsOq6ntW+quqG1fV9y8fv6p+MMnJ88u/6u7vDCwfAADYBhZx4YtHJjkryfOr6phMz7a6Q6ZnaJ2b5MnL9v/Q3NYO2+6a5KVVdWaSjyf5YpIbJblHpvu6/iXJE/fWGwAAALauhQtZ3f2xqjoyydOSHJcpGH0myfOSnNTdX1rDMO/J9Hys2yW5bZJrZbo88N+S/E2Sl3X3t/ZC+QAAwBa3cCErSbr7giTHr3HfWmHbvyV56OCyAAAAFuueLAAAgH2dkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADCQkAUAADDQQoasqrpBVb28qj5dVd+sqvOq6rlVde11jvP983HnzeN8eh73BnurdgAAYGvbf7MLWK+qummSs5JcN8lrknw4ye2TPDbJcVV1VHdftIZxrjOPc2iSNyf5qySHJTk+yT2r6k7d/fG98y4AAICtahFnsl6cKWA9prvv092/3d13T/KcJLdI8vQ1jvN7mQLWH3X3MfM498kU1q47nwcAAGBdFipkzbNYxyY5L8mLlnWfkOTSJA+uqgN2Mc41kjx43v/EZd0vTHJ+kv9aVTfZ86oBAIDtZKFCVpK7ze1p3X35jh3d/ZUkb09y9SR33MU4d0xytSRvn4/bcZzLk7xh2fkAAADWpLp7s2tYs6r6gyRPSPKE7v7DFfpfmORRSR7Z3S/ZyTiPyjRj9cLu/o0V+p+Q5A+S/H53/9YuanrPKl23vtrVrnalww8/fGeHD3H22Wfv9XMAACw54ogjNuQ8Z5999l929wM35GQw0KItfHHg3H55lf6l7Qdt0Dg7c9nXv/71L5999tnn7cEYa3XY3H54A861nfmcN4bPeWP4nDeOz3pj+Jw3xmFJcvbZZ/ucYScWLWTtc7r7dptdw9Js2r5Qy1bmc94YPueN4XPeOD7rjeFz3hg+Z1ibRbsna2mG6cBV+pe2X7xB4wAAAHyPRQtZ58ztoav033xuz92gcQAAAL7HooWs0+f22Kr6ntqr6ppJjkrytSTv3MU470zy9SRHzcftOM5+mZaJ3/F8AAAAa7JQIau7P5bktCSHZFpFcEcnJTkgyandfenSxqo6rKoO23HH7v5qklPn/U9cNs6j5/Hf0N0fH1g+AACwDSziwhePTHJWkudX1TFJPpTkDpmeaXVukicv2/9Dc1vLtj8pydFJ/ntV3SbJu5McnuTeST6XK4Y4AACAXVqo52QtqaobJnlakuOSXCfJZ5K8OslJ3f2lZft2knT38pCVqvr+JCckuU+SH0pyUZLXJXlqd39qL74FAABgi1rIkAUAALCvWqh7sgAAAPZ1QhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQhYAAMBAQtaCqqrrVNXDqurVVfXRqvp6VX25qs6sql+tKv9tB6mqZ1XVm6rqgvlz/mJVvbeqTqiq62x2fVtVVT2oqnr+ethm17NVVNV5O3yuy78u3Oz6tpqqOmb+//SFVfXNqvp0Vb2hqu6x2bUtuqp66E5+lpe+LtvsOreKqrpnVZ1WVZ+a/y38eFX9bVXdabNrg32R52QtqKp6RJKXZHoQ8+lJPpnkeknul+TAJK9Kcv/2H3iPVdW3kpyd5INJPpfkgCR3THJkkk8nuWN3X7B5FW498wPH/y3JlZJcI8nDu/tPNreqraGqzktyUJLnrtD91e5+9kbWs5VV1e8n+R9JPpXpQfdfSPKDSW6X5I3d/cRNLG/hVdVtktxnle6fSHL3JP/U3T+zUTVtVVX1rCRPTHJRkr/L9LN8syQ/m2T/JL/S3X++aQXCPkjIWlBVdfdMv+z/U3dfvsP2g5O8O8kNk/x8d79qk0rcMqrqqt39jRW2Pz3Jk5K8pLsfufGVbU1VVUn+OcmNk/x/SZ4QIWuYOWSluw/Z3Eq2tqp6eJI/TvJnSX6tu7+1rP/7uvvbm1LcNlBV78j0x7B7d/ffb3Y9i2z+veI/knw+yf/T3Z/boe9uSd6c5BPdfZNNKhH2SS4pW1Dd/ebu/ocdA9a8/cIkL51fHr3hhW1BKwWs2d/M7c03qpZt4jGZ/gJ9fJJLN7kWWLequkqSp2e6wuAKAStJBKy9p6p+NFPA+o8k/7TJ5WwFP5zp98V37RiwkqS7T0/ylUwztMAO9t/sAtgrlv7x/s6mVrH13Wtu37+pVWwhVXV4kmcmeV53v3WesWW8q1TVg5LcKFOQfX+St3a3+1fG+KlMv3Q+N8nlVXXPJLdK8o0k7+7ud2xibdvBr83tn/qZHuIjSb6V5PZV9QPd/YWljqq6S5JrZrqEENiBkLXFVNX+SX5lfvn6zaxlq6mqJ2S6P+jATPdj3TnTL6fP3My6tor5Z/fUTH/9f9Iml7PVHZzps97RJ6rq+O5+y2YUtMX82Nx+I8l7MwWs/1RVb810OffnN7qwra6qrpbkQUkuS+IS4wG6+4tV9VtJ/ijJB6vq7zLdm3XTTPdk/XOSX9+8CmHfJGRtPc/M9A/6a7v7DZtdzBbzhEyLiyx5fZKH+kVpmKcmuW2SO3f31ze7mC3s5CRvS/LvmS7zuUmSR2f66//rqupO3f2vm1jfVnDduf0fmRbM+Ykk78t0n+Gzkxyb5G/jku694RcyLezyTxYkGqe7nzvfz/nyJA/foeujSU5Zfhkh4J6sLaWqHpPk8Uk+nOTBm1zOltPdB3d3ZZoFuF+mX07fW1VHbG5li6+q7pBp9uoPXUq1d3X3SfM9nZ/t7q919we6+xGZ/kp9tSQnbm6FW8LSv63fSfKz3X1md3+1u/8tyX0zrTZ4V0tf7xVLlwq+bFOr2GKq6olJXpnklEwzWAdkWiXz40n+Yl5JE9iBkLVFVNWjkzwv019N79bdX9zkkras+ZfTV2f6a/R1krxik0taaPNlgq9Icm6Sp2xyOdvZ0oI5d9nUKraGi+f2vd193o4d3f21JEtXGdx+A2va8qrqlkl+PFOIfe0ml7NlVNXRSZ6V5O+7+79398fnP9CcnemPBv+R5PFVZXVB2IGQtQVU1W8meUGSD2QKWB4ougG6+/xMofaWVfUDm13PArtGkkOTHJ7kGzs+SDTJCfM+/3ve9tzNKnIbWLrs9YBNrWJrOGduL16l/0tze7W9X8q2YsGLvWPpOWOnL++Y/2jw7ky/T952I4uCfZ17shbcfDPqMzNd7/9TO676w4b4L3PrH/Td980kf7pK3xGZ/uE+M9Mvri4l3HvuOLcf39QqtoY3JekkP1JV+y1/1Ea+uxDGJza2rK2rqq6a6TL5y7L6/0/YPVeZ29WWaV/afoVHFcB2JmQtsKp6SpKnJXlPkmNdIjheVR2a5LPd/eVl2/dL8r8y3eB+Vnd/aaXj2bV5kYuHrdRXVSdmCll/5mHEe25eIv+T3X3psu2HJHnh/PLPN7quraa7z6+qf8i08tpjkzxnqa+qjk3yXzPNclkBdpz7J7l2kn+04MVwb8u8OE5Vvay7/2Opo6p+OslRmVbSPGuT6oN9kpC1oKrqIZkC1mWZ/gf4mKpavtt53X3KBpe21dwjyTOq6sxMf3W+KNMKg3fNtPDFhfnelZZgX/aATPdOvDXJ+ZlWF7xpknsmuWqm+1ievXnlbSmPyvQHgj+an5P13kyrC94n0/+3H7b8jzfskaVLBf94U6vYml6Z5I1JfjLJh6rq1Zn+7Ts806WEleS3u/uizSsR9j1C1uK68dxeKclvrrLPWzKtBMTue2OSm2V6JtZtMy0NfGmmRRpOTfJ8M4gskNOT3CLTz/JRme6/ujjT5ZinJjm1u3vTqttCuvtTVXW7TI8m+NlMC4pckuQfkjyju9+9mfVtJfMM7Z1jwYu9orsvr6p7ZPrDwS9mWuzi6km+mOnzfn53n7aJJcI+qfx7CgAAMI7VBQEAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAYSsgAAAAb6/wE7dM1cErB//gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 430.625x360 with 1 Axes>"
      ]
     },
     "execution_count": 11,
     "metadata": {
      "image/png": {
       "height": 352,
       "width": 428
      },
      "needs_background": "light"
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sns.displot(migraines)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Computing the _F_-like statistic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "In the next several exercises, you will compute the _F_-like statistic for your data. \n",
    "\n",
    "<h2>\n",
    "<p align=\"center\">\n",
    "    <i>F</i>-like=\\[\\frac{n_{a}|\\widetilde{a}-\\widetilde{G}|+n_{b}|\\widetilde{b}-\\widetilde{G}|+n_{c}|\\widetilde{c}-\\widetilde{G}|}{\\Sigma|a_{i}-\\widetilde{a}|+\\Sigma|b_{i}-\\widetilde{b}|+\\Sigma|c_{i}-\\widetilde{c}|}\\]\n",
    "</p>\n",
    "</h2>\n",
    "    \n",
    "Each quantity you compute in this section should be assigned to a variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "4. Find the median of each column. HINT: pandas data frames have a built-in function for columnwise medians, so you can just use `df.median()`, where `df` is the name of your data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Drug A    3.0\n",
       "Drug B    6.0\n",
       "Drug C    6.0\n",
       "dtype: float64"
      ]
     },
     "execution_count": 92,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meds = migraines.median()\n",
    "meds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "5. Find the grand median (the median of the whole sample). HINT: Use `np.median`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "grand = np.median(migraines)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "6. Find the numerator of the _F_-like statistic (variation among groups).  HINT: When working with data frames and NumPy arrays, you can do computations like addition and multiplication directly, without for loops (unlike in regular Python). Also, Numpy has `abs` and `sum` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36.0"
      ]
     },
     "execution_count": 27,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num = np.sum(len(migraines) * abs(meds - grand)) \n",
    "num"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "7. Compute the denominator of the _F_-like statistic. This represents variation within groups. HINT: Numpy and pandas will handle columns automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24.0"
      ]
     },
     "execution_count": 29,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "den = np.sum(np.sum(np.abs(migraines - meds)))\n",
    "den"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "8. Compute _F_-like, which is the ratio $\\frac{\\text{variation among groups}}{\\text{variation within groups}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5"
      ]
     },
     "execution_count": 30,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num/den"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Bootstrapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We now want to find a p-value for our data by simulating the null hypothesis. This, of course, means computing the _F_-like statistic each time, which takes a lot of code and would make a mess in the bootstrap loop. Instead, we will package our code into a function and call this function whereever necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "9. Write a function that will compute the _F_-like statistic for this dataset or one of the same size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "def Flike(df):\n",
    "    meds = df.median()\n",
    "    grand = np.median(df)\n",
    "    num = np.sum(len(df) * abs(meds - grand))\n",
    "    den = np.sum(np.sum(np.abs(df - meds)))\n",
    "    result = num/den\n",
    "    return(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "We now want to simulate the null hypothesis that there is no difference between the groups. To do this, we have to make all the data into one dataset, sample pseudo-groups from it, and compute the _F_-like statistic for the resampled data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "10. Use the code `alldata = np.concatenate([migraine[\"Drug A\"], migraine[\"Drug B\"], migraine[\"Drug C\"]])` (this assumes your data frame is called \"migraine\") to put all the data into one 1-D array. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "alldata = np.concatenate([migraines[\"Drug A\"], migraines[\"Drug B\"], migraines[\"Drug C\"]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "11. Sample three groups of the appropriate size from `alldata`. Assign each to a variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "A = np.random.choice(alldata, len(\"Drug A\"))\n",
    "B = np.random.choice(alldata, len(\"Drug B\"))\n",
    "C = np.random.choice(alldata, len(\"Drug C\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "12. Make the three samples into a data frame. To do this, use the NumPy function `column_stack` to put the 1-D arrays side by side and then use the pandas function `DataFrame` to convert the result into a data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8</td>\n",
       "      <td>5</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3</td>\n",
       "      <td>4</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   0  1  2\n",
       "0  3  5  5\n",
       "1  8  5  6\n",
       "2  5  3  4\n",
       "3  7  5  2\n",
       "4  3  4  6\n",
       "5  4  3  5"
      ]
     },
     "execution_count": 69,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resamps = np.column_stack([A, B, C])\n",
    "pd.DataFrame(resamps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "13. Compute the _F_-like statistic for your resampled data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "14. Do the above steps 10,000 times to simulate the null hypothesis, storing the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "sims = np.zeros(10000)\n",
    "\n",
    "for i in range(10000):\n",
    "    A = np.random.choice(alldata, len(migraines))\n",
    "    B = np.random.choice(alldata, len(migraines))\n",
    "    C = np.random.choice(alldata, len(migraines))\n",
    "    resampling = np.column_stack([A, B, C])\n",
    "    result = Flike(pd.DataFrame(resampling))\n",
    "    sims[i] = result\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7fe75c0b7610>"
      ]
     },
     "execution_count": 82,
     "metadata": {
     },
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAALACAYAAABy7a/RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAABYlAAAWJQFJUiTwAAA+E0lEQVR4nO3deZwlVX3//9dHkGFRBxASJxJFVBxElAAuDF8BITGgIKjwhbghKIpfFmXR8ENUNC4kIi6QSILKkJAEEL/idwygicMEFJSACwiyMyDaimyDM+OwzHx+f1S1XC/3dt/urr63b5/X8/G4j+JWnVN1qqjpeU/1qXMiM5EkSZJK8qRBN0CSJEnqN0OwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4QxeCI2K/iDgtIi6PiIciIiPinHHqrBUR74yIyyLigYj4XUTcHhHnRcSWXeocFBFXRcTyiFgWEUsiYq9xjnF0RFxb7//+iLgoIhZM9Zzr/f9rRPxrE/uSJEkq3dqDbsAknAi8BFgO3A3MH6twRDwF+AawG/Bj4GxgFfBM4JXAlsDNbXVOAY6t938msA5wILAoIo7MzNPbygdwLrAfcBNwOrAxcABwWUS8MTO/MekzrszfbrvttgPeNMX9SJIkzSYxmUrDGIKPpgqntwK7AJeOU/4fqQLwYZn5j+0bI+LJbd8XUAXg24CXZuYD9fpPA9cAp0TENzNzaUu1A6kC8BXA7pm5qq5zBvBd4MyIWJyZv53guUqSJGkaDF13iMy8NDNvycwcr2xEjD45Pa9TAK7392jbqsPq5SdGA3Bdbinw98Ac4OC2Ou+plyeOBuC6zv8A5wGbUoVkSZIkzQBDF4InaLTrwL9HxNyIeEtE/H8R8a6IeF6XOrvVy0s6bLu4rQwRsS6wAFgJXN5LHUmSJA3WMHaHmIiX1stnU3VveHrLtoyILwJHZeZqgIjYgKqv8PLMHOmwv1vqZevLdM8F1gJuz8zHeqwjSZKkAZrtIfiP6uWpwIVUL9XdDbwcOAP4P8BvgJPqcnPr5bIu+xtdv2HLusnU6SoirumyacwXACVJktS72d4dYvT8bgQOyMwbM3N5Zn6Hqo/uGuCYiFhnYC2UJElS3832J8EP1stFo10eRmXmTyLiDqruDFsBP+Hxp7Zz6Wx0/YMt6yZTp6vM3L7T+voJ8Xa97EOSJEljm+1Pgm+qlw922T46+sN6AJm5AvgF8JSImNeh/PPrZeu4wrcBq4EtIqLTPyo61ZEkSdIAzfYQ/F/18kXtGyJiDo8H1KUtmxbXyz067G/PtjLUQ6JdAaxPNfnGuHUkSZI0WLM9BH8N+CVwQES8rG3bh6i6Klyamb9qWX9GvfxgRGw0ujIiNgcOBx4Gzmrb1xfr5cfrIdNG67yUata439RtkSRJ0gwwdH2CI2JfYN/66zPq5Y4RsbD+73sz8zioujdExNuBbwKXR8T/peru8HLgfwH3AO9u3X9mXhERpwLHANdGxAVU0yYfQDUV8pFts8VBNWXyG6hetvtRRCyiGo7tAKrh0w7NzIemeu6SJElqxtCFYGBb4KC2dVvUH4A7geNGN2Tmf9ZPgT8E/DnV099fUT3x/ZvM/GX7ATLz2Ii4jurJ77uoRpH4IfDpzPxmh/IZEX9F1S3iEOBIYBVwGfDxzLxi0mcrSZKkxkUPsw9rBoiIa7bbbrvtrrmm2zDCkiRJRYrJVJrtfYIlSZKkJzAES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKs7ag26AZp75W2/DyMjImGXmzZvHjddf16cWSZIkNcsQrCcYGRlhz5MXjVnm4uP37lNrJEmSmmd3CEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOGsPugHSsJi/9TaMjIyMWWbevHnceP11fWqRJEmaLEOw1KORkRH2PHnRmGUuPn7vPrVGkiRNhd0hJEmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFGboQHBH7RcRpEXF5RDwUERkR50yg/pfqOhkRz+tSZq2IODoiro2I30XE/RFxUUQsGGO/60XERyPipohYFRH3RMT5EbHVZM5TkiRJ02foQjBwInAEsC3wi4lUjIi9gXcAy8coE8C5wKnAOsDpwNeBnYHLImKfDnXmAP8JfBh4CPg88F/A64GrI+LlE2mnJEmSptfag27AJBwN3A3cCuwCXNpLpYjYFDgTOA94Rl23kwOB/YArgN0zc1Vd/wzgu8CZEbE4M3/bUucYYCfgAuCAzFxT1zkPuBD4SkRsM7pekiRJgzV0T4Iz89LMvCUzc4JV/6leHj5OuffUyxNHA3B93P+hCtCbUoVk4PdPjg+rv36gNehm5jeAy4EX0j10S5Ikqc+GLgRPRkS8HdgXeHdm3jdGuXWBBcBKqvDa7uJ6uVvLuucCzwJuzsw7eqwjSZKkARrG7hATEhHPpuqje079ZHYszwXWAm7PzMc6bL+lXm7Zsu4F9fLmLvvsVKeriLimy6b5vdSXJEnS+Gb1k+CIeBJwNtWLcEf1UGVuvVzWZfvo+g2nWEeSJEkDNNufBB9N1Rf3tZn5wKAb04vM3L7T+voJ8XZ9bo4kSdKsNGufBEfElsAngLMy86Ieq40+tZ3bZfvo+genWEeSJEkDNGtDMNWIDHOAg1smx8iISB4fqeGWet2+9ffbgNXAFhHR6Sn58+tla//fm+pltz6/nepIkiRpgGZzd4ilwJe7bHst1VjBX6Wa3GIpQGauiogrgFfWn/YxiPesl4tb1t0G3AVsGRHP6TBCRKc6kiRJGqBZG4Iz88fAOztti4glVCH4hMy8tW3zF6kC8McjonWyjJcCBwC/Ab7WcpysJ9L4JPB3EdE6WcY+9b5uAP67ubOTJEnSVAxdCK67Luxbf31GvdwxIhbW/31vZh43hUOcC7yBakKMH0XEIuDpVAF4LeDQzHyorc6pwF51nR9ExHeoxg7en2rM4UOcLU6SJGnmGLoQDGwLHNS2bov6A3AnMOkQXD/Z/SuqaZMPAY4EVgGXAR/PzCs61Hk4Iv4COB74K6pRKR6imjL5I5l5w2TbI0mSpOYNXQjOzJOAk6a4j13H2f4Y8Nn60+s+VwIfrj+SJEmawWbz6BCSJElSR4ZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxRm6EBwR+0XEaRFxeUQ8FBEZEed0Kfv8iPjriFgcET+PiEci4tcR8Y2IeNU4xzkoIq6KiOURsSwilkTEXmOUXysijo6IayPidxFxf0RcFBELpnrOkiRJatbQhWDgROAIYFvgF+OU/RvgZOCPgYuAzwDfA14LLI6IozpViohTgIXAPOBM4BxgG2BRRBzRoXwA5wKnAusApwNfB3YGLouIfSZygpIkSZpeaw+6AZNwNHA3cCuwC3DpGGUvAf42M3/UujIidgH+E/h0RHw1M0dati0AjgVuA16amQ/U6z8NXAOcEhHfzMylLbs8ENgPuALYPTNX1XXOAL4LnBkRizPzt5M/bUmSJDVl6J4EZ+almXlLZmYPZRe2B+B6/X8DS6ie2rZ3VzisXn5iNADXdZYCfw/MAQ5uq/OeenniaACu6/wPcB6wKVVIliRJ0gwwdCG4QY/Wy8fa1u9WLy/pUOfitjJExLpUQXolcHkvdSRJkjRYRYbgiHg2sDtVcL2sZf0GwDOB5a1dJFrcUi+3bFn3XGAt4PbMbA/U3epIkiRpgIaxT/CURMQc4F+pujV8oLXLAzC3Xi7rUn10/YZTrDNW+67psml+L/UlSZI0vqKeBEfEWsC/ADtR9dU9ZbAtkiRJ0iAU8yS4DsDnAPsD5wNv6fBy3ehT27l0Nrr+wSnW6Sozt++0vn5CvF0v+5AkSdLYingSHBFPBv6daiizfwPe1Kn/bmauoBp7+CkRMa/Drp5fL29uWXcbsBrYIiI6/aOiUx1JkiQN0KwPwRGxDvBVqifA/wy8NTNXj1Flcb3co8O2PdvKUA+JdgWwPvDKXupIkiRpsGZ1CK5fgvs6sA/wZeDgzFwzTrUz6uUHI2Kjln1tDhwOPAyc1Vbni/Xy4/WQaaN1XgocAPwG+NokT0OSJEkNG7o+wRGxL7Bv/fUZ9XLHiFhY//e9mXlc/d9nAK8B7qXq5vDhaobjP7AkM5eMfsnMKyLiVOAY4NqIuIBqUo0DgI2BI9tmi4NqyuQ3UE2I8aOIWAQ8va6zFnBoZj40uTOWJElS04YuBAPbAge1rdui/gDcCYyG4OfUy02AD4+xzyWtXzLz2Ii4jurJ77uANcAPgU9n5jfbK2dmRsRfUXWLOAQ4ElhFNQbxxzPzil5OTJIkSf0xdCE4M08CTuqx7K5TOM5CYOEEyj8GfLb+SJIkaQab1X2CJUmSpE4MwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUnLUH3QBpNlmxYiVzN96k6/Z58+Zx4/XX9bFFkiSpE0Ow1KA1a1az58mLum6/+Pi9+9gaSZLUjd0hJEmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFWfoQnBE7BcRp0XE5RHxUERkRJwzTp0FEXFRRNwfEb+LiGsj4n0RsdYYdfaKiCURsSwilkfEDyLioHGOc1BEXFWXX1bX32uy5ypJkqTpMXQhGDgROALYFvjFeIUjYh/gMmBn4OvA6cA6wGeBc7vUOQJYBLwIOAc4E/gTYGFEnNKlzinAQmBeXf4cYBtgUb0/SZIkzRDDGIKPBrYEnga8Z6yCEfE0qkC6Gtg1M9+Rme+nCtBXAvtFxIFtdTYHTgHuB3bIzMMz82jgxcBtwLERsWNbnQXAsfX2F2fm0Zl5OLB9vZ9T6v1KkiRpBhi6EJyZl2bmLZmZPRTfD9gUODczr27ZxyqqJ8rwxCB9CDAHOD0zl7bUeQD4ZP31sLY6o98/UZcbrbMU+Pt6fwf30F5JkiT1wdCF4AnarV5e0mHbZcBKYEFEzOmxzsVtZaZSR5IkSQMy20PwC+rlze0bMvMx4A5gbWCLHuuMACuAzSJifYCI2AB4JrC83t7ulnq55WROQJIkSc1be9ANmGZz6+WyLttH1284wTob1OVWTvIYXUXENV02ze+lviRJksY3258ES5IkSU8w258Ejz6Fndtl++j6B9vqbFJvu2+MOsvalhM5RleZuX2n9fUT4u162YckSZLGNtufBN9UL5/QHzci1gaeAzwG3N5jnXlUXSHuzsyVAJm5gmq84qfU29s9v14+oY+xJEmSBmO2h+DF9XKPDtt2BtYHrsjMh3uss2dbmanUkSRJ0oDM9hB8AXAvcGBE7DC6MiLWBT5ef/1iW52zgIeBI1onuIiIjYAT6q9ntNUZ/f7Butxonc2Bw+v9nTWVE5EkSVJzhq5PcETsC+xbf31GvdwxIhbW/31vZh4HkJkPRcShVGF4SUScSzWD2+uohkK7ADivdf+ZeUdEvB/4AnB1RJwHPEI18cZmwGcy88q2OldExKnAMcC1EXEB1dTMBwAbA0e2TrwhSZKkwRq6EEw15fFBbeu24PGxfu8EjhvdkJkXRsQuwAeBNwLrArdSBdYvdJp5LjNPi4il9X7eRvXE/AbgxMw8u1OjMvPYiLiO6snvu4A1wA+BT2fmNyd1ppIkSZoWQxeCM/Mk4KQJ1vke8JoJ1lkELJpgnYXAwonUkSRJUv/N9j7BkiRJ0hMYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiNBqCI+JZEfG0cco8NSKe1eRxJUmSpIlo+knwHcB7xylzVF1OkiRJGoimQ3DUH0mSJGnGGkSf4GcAKwZwXEmSJAmAtae6g4h4W9uqbTusA1gLeBbwFuC6qR5XkiRJmqwph2BgIZD1fyewT/1pN9pNYiXw0QaOK0mSJE1KEyH44HoZwFeAC4FvdCi3GrgPuDIzH2zguJIkSdKkTDkEZ+bZo/8dEQcBF2bmP091v5IkSdJ0aeJJ8O9l5qua3J8kSZI0HZwxTpIkScVpPARHxC4R8c2IuCciHo2I1R0+jzV9XEmSJKlXjXaHiIjXUr0YtxZwF3ATYOCVJEnSjNJoCAZOAh4FXpuZ325435IkSVIjmu4O8SLgPAOwJEmSZrKmQ/By4P6G9ylJkiQ1qukQ/B1gx4b3KUmSJDWq6RD818BzI+LEiIhxS0uSJEkD0PSLcR8Brgc+ChwSET8GHuxQLjPzHQ0fW5IkSepJ0yH47S3/vXn96SQBQ7AkSZIGoukQ/JyG9ydJkiQ1rtEQnJl3Nrk/SZIkaTo0Pm2yJEmSNNM1PW3ys3otm5l3NXlsSZIkqVdN9wleSvXS23hyGo4tSZIk9aTpIPrPdA7BGwLbAs8GlgB97zscEa8F3gu8EHg6MAJcA5yamVd2KL8AOBF4BbAecAvwFeC0zFzd5Rh7AccBfwasRTVc3D9k5tmNn5AkSZImrekX497ebVtEPAn4EHAYcFCTxx1PRPwt8AHgPuBC4F7gecA+wBsj4m2ZeU5L+X2ArwGrgPOopoLeG/gssBOwf4djHAGcVh/jHOARYD9gYURsk5nHTdf5SZIkaWL69mJcZq7JzI9SdZk4uV/HjYhnUD2d/TXwwsx8Z2Yen5n7AX8JBPCxlvJPA84EVgO7ZuY7MvP9VE+yrwT2i4gD246xOXAKVVjeITMPz8yjgRcDtwHHRoTTSUuSJM0Qgxgd4grg1X083rOpzvMHmXlP64bMvBT4LbBpy+r96u/nZubVLWVXUXWPAHhP2zEOAeYAp2fm0pY6DwCfrL8eNuUzkSRJUiMGEYI3Bjbo4/Fuoeqa8LKI2KR1Q0TsDDwV+K+W1bvVy0s67OsyYCWwICLm9Fjn4rYykiRJGrC+jtAQEX8OHAD8tF/HzMz7I+KvgVOBGyLiQqp+u88FXgf8J/DuliovqJc3d9jXYxFxB7A1sAXwsx7qjETECmCziFg/M1eO1d6IuKbLpvlj1ZMkSVLvmh4nePEYx/lTYHQc4Y91KTctMvNzEbGUanSHQ1s23QosbOsmMbdeLuuyu9H1G06wzgZ1uTFDsCRJkqZf00+Cd+2yPoEHgG8Bp2Rmt7A8LSLiA1R9c78AnA78iurJ6qeAf42IbTPzA/1sUzeZuX2n9fUT4u363BxJkqRZqekh0mbcNMwRsSvwt8DXM/OYlk0/jIjXU3VhODYizsjM23n8ae5cOhtd/2DLumXAJvW2+8ao0+1JsSRJkvpoxoXWabBXvby0fUPdP/cqquvwZ/Xqm+rllu3lI2Jt4DnAY8DtLZvGqjOPqivE3eP1B5YkSVJ/TGsIjoinRsSf1mPvDsroKA6bdtk+uv6RejnaVWOPDmV3BtYHrsjMh1vWj1Vnz7YykiRJGrDGQ3BErB0Rx0fErVRdBpYCD0TErfX6vo5IAVxeL98VEc9s3RARe1LNALeKavxigAuoZpQ7MCJ2aCm7LvDx+usX245xFvAwcEQ9ccZonY2AE+qvZ0z5TCRJktSIpkeHWIdqrNxdqF6G+zkwAswDNgc+AewREa/OzEe67adhF1CNA/znwM8i4utUL8ZtRdVVIoDjM/M+gMx8KCIOrestiYhzqWaCex3VUGgXUE2l/HuZeUdEvJ/qxburI+I8Hp82eTPgM5l55bSfqSRJknrS9JPgY6hGiPgPYKvM3Dwzd8zMzakC5CLglXW5vsjMNcBrgKOBG4DXA8cCrwAuAv4yMz/fVudCqiB/GfBG4Ejg0brdB2ZmdjjOaVRB+XrgbcC7qML22zPzuOk4N0mSJE1O010T3kQ1Eca+dfj8vcy8LSLeAPwYeDNwcsPH7iozHwU+V396rfM9qvA8keMsogr6kiRJmsGafhL8PODi9gA8ql5/MdVsbZIkSdJANB2CHwGeMk6ZDai6FkiSJEkD0XQIvhbYLyI6DkcWEZtQvSz2k4aPK0mSJPWs6RB8OtW4u1dFxDsiYouIWC8inhMRBwM/qLef3vBxJUmSpJ41PW3y+RGxLXA88E8digTwd5l5fpPHlSRJkiai8YkrMvOEiPh/wDuopiKeCywDfgR8xfFyJUmSNGjTMntbZn4f+P507FuSJEmaqin3CY6IdSLiqoj4TkQ8eZxy34mI749VTpIkSZpuTbwY9xZge6qpgbsOfVZPk/xp4GVUk2VIkiRJA9FEd4g3ALdn5kXjFczMSyLiFmB/YGEDx5aGyooVK5m78SZjlpk3bx43Xn9dn1okSVKZmgjBfwaMG4BbXMYEpyOWZos1a1az58ljz6x98fF796k1kiSVq4nuEJsAv55A+V8DT2/guJIkSdKkNBGCf8f4UyW3egqwqoHjSpIkSZPSRAj+ObDDBMrvANzVwHElSZKkSWkiBC8BdoyIcYNwRGwPLAAubeC4kiRJ0qQ0EYJPBxL4akRs1a1QRMwHvgqsBv6hgeNKkiRJkzLl0SEy86aI+BhwEvCjiLgAWAzcXRd5JrA78EZgDvDhzLxpqseVJEmSJquRaZMz82MR8RjwEeBNwF+1FQngUeCDmfmpJo4pSZIkTVYjIRggMz8ZEf8KHALsBMyrN40A3wXOysw7mzqeJEmSNFmNhWCAOuR+pMl9SpIkSU1r4sU4SZIkaagYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnGKCsERsXtEfD0ifhURD0fELyPiWxHxmg5lF0TERRFxf0T8LiKujYj3RcRaY+x/r4hYEhHLImJ5RPwgIg6a3rOSJEnSRBUTgiPi74D/AnYA/h/wGeA/gE2BXdvK7gNcBuwMfB04HVgH+Cxwbpf9HwEsAl4EnAOcCfwJsDAiTmn8hCRJkjRpaw+6Af0QEYcC7wfOBt6VmY+0bX9yy38/jSrArgZ2zcyr6/UfAhYD+0XEgZl5bkudzYFTgPuBHTJzab3+Y8D/AMdGxNcy88ppO0lJkiT1bNY/CY6IOcAngLvoEIABMvPRlq/7UT0dPnc0ANdlVgEn1l/f07aLQ4A5wOmjAbiu8wDwyfrrYVM7E0mSJDWlhCfBf0EVaj8HrImI11J1WVgFXNXh6exu9fKSDvu6DFgJLIiIOZn5cA91Lm4rI0mSpAErIQS/tF6uAn5EFYB/LyIuA/bLzN/Uq15QL29u31FmPhYRdwBbA1sAP+uhzkhErAA2i4j1M3PlVE5GkiRJU1dCCP6jevl+4AbglcCPgedQ9eN9NfBVHn85bm69XNZlf6PrN2xZ10udDepyY4bgiLimy6b5Y9WTJElS72Z9n2AeP8fHgNdl5nczc3lmXge8Hrgb2CUidhxYCyVJktRXJTwJfrBe/qj1pTWAzFwZEd8C3gG8DLiSx5/mzqWz0fUPtqxbBmxSb7tvjDrdnhS3tmn7TuvrJ8TbjVdfkiRJ4yvhSfBN9fLBLtsfqJfrtZXfsr1gRKxN1Y3iMeD2DsfoVGceVVeIu+0PLEmSNDOUEIK/AyTwwojodL6jL8rdUS8X18s9OpTdGVgfuKJlZIjx6uzZVkaSJEkDNutDcGbeSTWT27OA97Zui4hXA39J9ZR4dHizC4B7gQMjYoeWsusCH6+/frHtMGcBDwNH1BNnjNbZCDih/nrG1M9GkiRJTSihTzDA4cCfAafW4wT/iKpbw75UM8O9MzOXAWTmQ/UMcxcASyLiXKqZ4F5HNRTaBcB5rTvPzDsi4v3AF4CrI+I84BGqiTc2Az7jbHGSJEkzRxEhODPvjojtgQ9ThdmdgYeonhB/KjOvait/YUTsAnwQeCOwLnArcAzwhczMDsc4LSKWAscBb6N6yn4DcGJmnj1d5yZJkqSJKyIEA9STYRxZf3op/z3gNRM8xiKqYC0VY/7W2zAyMtJ1+7x587jx+uv62CJJksZXTAiWND1GRkbY8+Tu//a7+Pi9+9gaSZJ6M+tfjJMkSZLaGYIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4qw96AZIUi/mb70NIyMjXbfPmzePG6+/ro8tkiQNM0OwpKEwMjLCnicv6rr94uP37mNrJEnDzu4QkiRJKo4hWJIkScWxO4SkabVixUrmbrzJmGXszytJ6jdDsKRptWbN6jH78oL9eSVJ/Wd3CEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEcHUJi/Cl5AZYvX9Gn1kiSpOlmCJYYf0pegPMP37U/jZEkSdPO7hCSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSpOkSE4It4SEVl/3tmlzF4RsSQilkXE8oj4QUQcNM5+D4qIq+ryy+r6e03PWUiSJGmyigvBEfGnwOnA8jHKHAEsAl4EnAOcCfwJsDAiTulS5xRgITCvLn8OsA2wqN6fJEmSZoiiQnBEBHAWcB9wRpcymwOnAPcDO2Tm4Zl5NPBi4Dbg2IjYsa3OAuDYevuLM/PozDwc2L7ezyn1fiVJkjQDFBWCgaOA3YCDgRVdyhwCzAFOz8yloysz8wHgk/XXw9rqjH7/RF1utM5S4O/r/R08xbZLkiSpIcWE4IjYCjgZ+HxmXjZG0d3q5SUdtl3cVmYqdSRJkjQgRYTgiFgb+BfgLuCEcYq/oF7e3L4hM0eoniBvFhHr1/veAHgmsLze3u6WernlJJouSZKkabD2oBvQJx8G/gz4X5n5u3HKzq2Xy7psXwZsUJdb2WN5gA17aWhEXNNl0/xe6kuSJGl8s/5JcES8nOrp72cy88pBt0eSJEmDN6ufBNfdIP6ZqmvDh3qstgzYhOoJ730dtrc/+V3Wtr5b+Qd7OXhmbt9pff2EeLte9iFJkqSxzfYnwU+h6ou7FbCqZYKMBD5SlzmzXve5+vtN9fIJfXgjYh5VV4i7M3MlQGauAH4BPKXe3u759fIJfYwlSZI0GLP6STDwMPDlLtu2o+on/F2q4DvaVWIxsBOwR8u6UXu2lGm1GHhrXeesHutIkiRpQGZ1CK5fgus2LfJJVCH47Mz8Usums4APAEdExFmjYwVHxEY8PrJE+0QbZ1CF4A9GxIWjYwXXE2QcThXG28OxJEmSBmRWh+DJyMw7IuL9wBeAqyPiPOARYD9gMzq8YJeZV0TEqcAxwLURcQGwDnAAsDFwZOvEG5IkSRosQ3AHmXlaRCwFjgPeRtV3+gbgxMw8u0udYyPiOqonv+8C1gA/BD6dmd/sS8MlSZLUk2JDcGaeBJw0xvZFwKIJ7nMhsHAKzZIkSVIfzPbRISRJkqQnMARLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqztqDboAk9cv8rbdhZGRkzDLz5s3jxuuv61OLJEmDYgiWVIyRkRH2PHnRmGUuPn7vPrVGkjRIdoeQJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBVn7UE3QNIfWrFiJXM33mTMMvPmzePG66/rU4skSZp9DMHSDLNmzWr2PHnRmGUuPn7vPrVGkqTZye4QkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOIViSJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKs/agGyBJK1asZO7Gm4xZZvnyFX1qjSSpBIZgSQO3Zs1q9jx50Zhlzj981/40RpJUBLtDSJIkqTiGYEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSrOrA/BEfH0iHhnRHw9Im6NiN9FxLKI+G5EvCMiOl6DiFgQERdFxP11nWsj4n0RsdYYx9orIpbU+18eET+IiIOm7+wkSZI0GSVMlrE/8EVgBLgUuAv4Y+ANwJeAPSNi/8zM0QoRsQ/wNWAVcB5wP7A38Flgp3qffyAijgBOA+4DzgEeAfYDFkbENpl53HSdoCRJkiamhBB8M/A64D8yc83oyog4AbgKeCNVIP5avf5pwJnAamDXzLy6Xv8hYDGwX0QcmJnntuxrc+AUqrC8Q2Yurdd/DPgf4NiI+FpmXjm9pypJkqRezPruEJm5ODMXtQbgev2vgDPqr7u2bNoP2BQ4dzQA1+VXASfWX9/TdphDgDnA6aMBuK7zAPDJ+uthUzsTSWNZsWIlczfeZMzP8uUrBt1MSdIMUcKT4LE8Wi8fa1m3W728pEP5y4CVwIKImJOZD/dQ5+K2MpKmwZo1q9nz5EVjljn/8F370xhJ0oxXbAiOiLWBt9VfW8PrC+rlze11MvOxiLgD2BrYAvhZD3VGImIFsFlErJ+ZK8dp1zVdNs0fq54kSZJ6N+u7Q4zhZOBFwEWZ+a2W9XPr5bIu9UbXbziJOnO7bJckSVIfFfkkOCKOAo4FbgTeOuDm/IHM3L7T+voJ8XZ9bo4kSdKsVNyT4Hoos88DNwCvysz724qM99R2dP2Dk6jT7UmxJEmS+qioJ8ER8T6qsX5/Cuyemfd0KHYTsAOwJfAH/XPrfsTPoXqR7va2OpvUda5sqzMP2AC4e7z+wJLUav7W2zAyMtJ1+7x587jx+uv62CJJmj2KCcER8ddU/YB/DPxFZt7bpehi4M3AHsC/t23bGVgfuKxlZIjROjvVddrHAt6zpYw0VMYLYYDDjk2jkZGRMUe8uPj4vfvYGkmaXYoIwfVEFx+jerL76g5dIFpdAPwtcGBEnNYyWca6wMfrMl9sq3MW8AHgiIg4q2WyjI2AE+oyZyANmfFCGDjsmCRpOM36EBwRB1EF4NXA5cBREdFebGlmLgTIzIci4lCqMLwkIs6lmgnudVRDoV1ANZXy72XmHRHxfuALwNURcR6PT5u8GfAZZ4uTJEmaOWZ9CKbqwwuwFvC+LmX+G1g4+iUzL4yIXYAPUk2rvC5wK3AM8IXMzPYdZOZpEbEUOI5q/OEnUb18d2Jmnt3EiUiSJKkZsz4EZ+ZJwEmTqPc94DUTrLMIGPt3x5IkSRq44oZIkyRJkgzBkiRJKo4hWJIkScUxBEuSJKk4s/7FOEmaiBUrVjJ3403GLONMbZI0/AzBktRizZrV404Q4kxtkjT87A4hSZKk4hiCJUmSVBxDsCRJkopjn2BJatj8rbdhZGRkzDK+XCdJg2UIloaQIxjMbCMjI75cJ0kznCFYGkKOYCBJ0tTYJ1iSJEnFMQRLkiSpOIZgSZIkFcc+wVKhehnBYPnyFX1qjSRJ/WUIlgrVywgG5x++a38aI0lSn9kdQpIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxHh5CkCVqxYiVzN96k63aHlpOkmc8QLEkTtGbN6jGHl3NoOUma+ewOIUmSpOL4JFiapfyV/ew33v9jgHnz5nHj9df1qUWSNDwMwdIs5a/sZ7/x/h8DXHz83n1qjSQNF0OwJA1AL09xfVovSdPHECxJA9DLU1yf1kvS9PHFOEmSJBXHECxJkqTiGIIlSZJUHEOwJEmSimMIliRJUnEMwZIkSSqOQ6RJ0izmrHKS1JkhWJJmMWeVk6TO7A4hSZKk4vgkWJIKZ5cJSSUyBEtS4ewyIalEdoeQJElScQzBkiRJKo4hWJIkScUxBEuSJKk4vhgnSRrXeCNIrHr4YdadM2fMfTjChKSZxBAsSRrXeCNInH/4rrz+s98ecx+OMCFpJrE7hCRJkorjk2AVYf7W2zAyMtJ1+/LlK/rYGqlMvUzKMV63CrtUSGqKIVhFGBkZGfdXuZKmVy+TcozXrcIuFZKaYncISZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScVxnGBJ0qwy3uQ4MP6kHODEHNJsZwiWJM0q402OA+NPygFOzCHNdnaHkCRJUnEMwZIkSSqOIViSJEnFsU+wJEmT0MsLeL5cJ81chmBJkiahlxfwLjhyd+ZuvMmYZQzK0mAYgiVJmiZr1qweNyg7CoU0GPYJblBEbBYRX4mIX0bEwxGxNCI+FxEbDbptkiRJepxPghsSEc8FrgD+CPgGcCPwMuC9wB4RsVNm3jfAJkrS0FuxYuW43QuWL1/Rp9Y0Y7xzsruEND0Mwc35B6oAfFRmnja6MiJOBY4GPgEcNqC2SdKs0Ev3gvMP37U/jWnIeOdkdwlpehiCG1A/BX41sBT4+7bNHwHeBbw1Io7NzOF6RCFJmhXGG83CJ84qjSG4Ga+ql9/OzDWtGzLztxHxPaqQ/ArgO/1unCRpdutluLbly1ew/+mLu27v5YlzL8dZ9fDDrDtnzphlDNyaCQzBzXhBvby5y/ZbqELwlhiCJUkT0Gs/6LECLozfTaTJ47z+s98es8x4gbupsD1eGcN42SIzB92GoRcR/wQcChyamV/qsP0TwAnACZn5qXH2dU2XTS9Zb7311tpqq62m3N7x/PgnP2HuM583Zpllv7iVbV/ykmlvS1PGO6cH7rqZjZ615Zj7aKJMv44zk9oy244zk9riOU9zmZ/fzJOe1H0QpTWr18yqcy7x2jb1d9n1N9zAo48+OmaZJz/5yWz9whdO+VgzxUw65x/+8If/lplvnmg9Q3AD+hSCXwQsp+p3PJ3m18sbp/k4w87r1BuvU2+8TuPzGvXG69Qbr1NvhuU63TiZEGx3iGYsq5dzu2wfXf/geDvKzO2baNBkjYbwQbdjpvM69cbr1Buv0/i8Rr3xOvXG69Sb2X6dnCyjGTfVy26/c3l+vezWZ1iSJEl9ZAhuxqX18tUR8QfXNCKeCuwErAS+3++GSZIk6YkMwQ3IzNuAbwObA4e3bf4osAHwL44RLEmSNDPYJ7g5/4dq2uQvRMTuwM+Al1ONIXwz8MEBtk2SJEktfBLckPpp8A7AQqrweyzwXODzwCsy877BtU6SJEmtHCJNkiRJxfFJsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxTEES5IkqTiGYEmSJBXHECxJkqTiGIJniYjYLCK+EhG/jIiHI2JpRHwuIjaa4H42rustrffzy3q/m033sfthqm2NiA0i4s0R8W8RcWNErIiI30bE1RFxbESs06VejvH5frNnOXVN/D+NiCXjnPe6Xeq9MCLOj4h7ImJVRNwUER+NiPWaO8NmNHA/7TrONRr9/GlbvaG4nyJiv4g4LSIuj4iH6vadM8l9TfhaD8u91MR1ioinR8Q7I+LrEXFrRPwuIpZFxHcj4h0R8YS/7yNi83HupXObO8upa+p+qu+dbuf8qzHqLYiIiyLi/vr6XhsR74uItaZ2Zs1q6H56ew8/l1a31Rmq+wmcNnlWiIjnUk3Z/EfAN4AbgZcB7wX2iIidepmxLiKeXu9nS2AxcC4wHzgYeG1E7JiZt0/Hsfuhoba+EjgHuB+4FLgQ2Ah4HXAK8IaI2D0zV3WoeyfVjILt7p7wyUyjafh/+tEu6x/rcOyXU917TwYuAH4O7AZ8GNi9vrYPT+DY06ah67SU7tdnG+ANwE8z8+cdtg/D/XQi8BJgOVW75k9mJ5O51sN0L9HMddof+CIwQvWz6S7gj6nuoS8Be0bE/tl5hqyfUP0sa/fTSbRjOjVyP9WWAZ/rsH55p8IRsQ/wNWAVcB7V3wF7A58FdqK6/jNFE9fpx3T/2fRKqj9LF3fZPiz3E2SmnyH/AN8CEjiybf2p9fozetzPP9blP9O2/qh6/SXTdexhuU7AtsCbgXXa1j8VuKbez7Ed6iWwZNDXoM/305LqR0zPx10LuKE+xuta1j+JKsQkcPygr0/T12mM/f97vZ+jhvV+Al4FPB8IYNe63edM97UewntpyteJKpTsDTypbf0zqAJxAm9s27Z5vX7hoK9Bn++npcDSCZR/GnAP8DCwQ8v6dan+cZbAgYO+Pk1fpzH2f2X7n61hvJ8y0xA87B/gufVNd0eHH35PpfqX4Apgg3H28xRgZV3+qW3bnlT/0Ehgi6aPPUzXaZxjvKk+xqIO24YltDR2nZh4CN6tPvZ/d9i2Rb1tKfV077PlOnXZ/yZUT5xWAhsO6/3U1uZJ/WU8mWs9TPdSU9dpnH2eUO/ztLb1mzNkoaWJ68TEQ/Ah9bHO7rCt6702Ez5N309Uv6FKqifMaw37/WSf4OH3qnr57cxc07ohM38LfA9YH3jFOPt5BbAe8L26Xut+1lA9iWk9XpPH7od+tPXRevmEX/PXNoyIQyLihIg4PCJmwnVp1/h1iogDIuL4iDgmIvaMiDldiu5WLy9p35BVN5ybgWdThZhBm+776SBgDvDVzHywS5lhuJ+aMJlrPUz3Uj+M97PpTyLi3fW99O6IeHG/GjZAcyLiLfU5vzciXjVG396u9xNwGdU/VheM8bNtNnlXvfxyZq7uUmZo7if7BA+/F9TLm7tsvwV4NVU/3+9McT/U+2n62P3Qj7YeUi87/aCEqo/Wl1tXRMRPgLdm5nWTPGbTpuM6tb8McU9EHJ6ZF0zi2FvWn9t6PPZ0me776dB6+Y9jlBmG+6kJk7nWw3QvTauIWBt4W/2128+mv6g/rfWWAAdl5l3T17qBegbwL23r7oiIgzPzv9vWd72fMvOxiLgD2JrqH1U/a7ylM0T9QulbgNVU/cy7GZr7ySfBw29uvVzWZfvo+g2nYT9NHbsfprWtEXEEsAfVywRf6VDkVKqXJzal+hXuS6n6Jr4EWBwRz5zMcadBk9fpG1R9FDej+i3DfOBTdd3zImKPaTz2dJu2tkbELlR/6f40M6/oUmxY7qcmzPafTdPtZOBFwEWZ+a22bSuBvwG2p3rBdyNgF6oX63YFvhMRG/SvqX1zFrA7VRDegOpX/P9I9ev8iyPiJW3lvZ8q/5vqHC/Jzi/rDt39ZAiWpigi3kD1lvGvqF48ebS9TGYem5lXZOa9mbk8M6/OzP2p3jbeBDiur43ug8z8bGZ+MzN/kZmrMvOmzDwBOJbqZ8+nBtzEmWr0143/1K1AifeTJi4ijqL683Yj8Nb27Zl5T2Z+ODN/mJkP1p/LqJ6s/wB4HvDOvja6DzLzo5m5ODN/nZkrM/OnmXkY1T8u1wNOGmwLZ6zRn00df0M1jPeTIXj4jf4LdG6X7aPrH5yG/TR17H6YlrZGxL5Uv+6/B9g124aQ68EZ9XLnCdabLv34f/olqr6J20bEU/t87KZM1/20MfBG4Hc88Ve1vZhp91MTZvvPpmlR/3bq81SjZLwqM+/vtW5mPsbjv+6eTffSeLr9+fF+itgaWED1QtxFE6k7k+8nQ/Dwu6lebtll+/PrZbe+cVPZT1PH7ofG2xoR+wNfBX4N7JKZN41TpZPf1MuZ8iuiaf9/mtUYyqMvX7aed9H3U230hbjzx3ghbiwz7X5qwmz/2dS4iHgfcBrVuKyvysyuE0CMYTbeS+Ppds5d76e6z/VzqP5hP9GHIMOklxfixjIj7ydD8PC7tF6+OtpmBKqfsu1E1U9nvFmkvk/19Gmntqdz1Pt9ddvxmjx2PzTa1oh4M9U4rr+kCsC3jFOlm9E32mfKD89p/38aES+g6iv2W+Delk2L62V7X2EiYguqv4DuZGZcq+m6TqMvxHXtCjGOmXY/NWEy13qY7qVGRcRfU03g8GOqAHzPJHc1G++l8XQ75673E9WTzfWBK3LmTL7SqKhm93wr1QtxXx6neDcz834a9Bhtfqb+YeIDyc8H5nfYj5Nl9HadDqL6YXA78Owejvti4Mld1t9bH/tNg74+TV4nqicjG3fY96Y8Prj8P7VtG2uCg68y8yY4aOR+atn+yrredbPpfmpp366MMV4p1cxu84HnNnCth+peavA6faiue3WnP38dym9H29jL9frdqcapTmDBoK9Jk9cJ2IoO43dTvRR3S73PE9q2PY3qSeZQTJbR1P3UUuatdBkDf9jvp6gbqCHWYUrRnwEvpxpf82aqm+6+lvIJkJnRtp/2aZOvovqBsQ9Vn9cFmXlbW50JHXuQmrhOEfEq4L+o/jL9CtVUrO0ezMzPtdRZSDVKwuV1+YepfujsQfWX9ZnAu3OG/GFs6Dq9nap/3Xep/rFwP/As4DVU/eeuBv4i237l32Gq27uofoDuQDUe7IyZ6rapP3ct2/+FavihozLztDGOu5AhuZ/qPvP71l+fAfwl1f1web3u3sw8ri67OdWEGHdm5uZt+5nwz5khu5f2ZYrXKSIOoppGezVVV4hOIxkszcyFLXWWUHUNuYLHp9t+MY+Pi/uhzPz45M+sWQ1dp5OoXha8jOq3Ab+lmpDltVSh9iLg9Zn5SIdjX0AV5s6l+pn2OqqRXC4A/vdM+DMHzf25a9nf5cD/ovoH5aIxjruEIbqfAJ8Ez5YP8KdUw76MAI9Q/eH+HLBRh7JJl5m8gI2pXqa4s97PCFXY26yJYw/6M9XrBLx9dP0Yn6VtdfYF/i9wK/BQy3VdRNu0kzPl08B12obqL+TrgPuoBuu/n+qH8JG0TTvdVveFVE/r7qUKeDdTzWG/3qCvS9PXqWXbRlTdkTrOEDes9xPVW/Y9/Vnh8dmmlnbZ14R/zgzLvdTEdephH0nbLIPAO4BvUs2gtry+RncB5wGvHPR1mabrtAtVV7YbqV5ke5TqKe9/Uo2n3HUWQaquNxcBD9R/Xq8DjqZt5rRBfxr+c7dVvf3n453nsN1PmT4JliRJUoF8MU6SJEnFMQRLkiSpOIZgSZIkFccQLEmSpOIYgiVJklQcQ7AkSZKKYwiWJElScQzBkiRJKo4hWJIkScUxBEuSJKk4hmBJkiQVxxAsSZKk4hiCJUmSVBxDsCRJkopjCJYkSVJxDMGSJEkqjiFYkiRJxfn/ARg5HKybqGQ2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "execution_count": 82,
     "metadata": {
      "image/png": {
       "height": 352,
       "width": 352
      },
      "needs_background": "light"
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sns.displot(sims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "15. Find the p-value for your data. What do you conclude about the migraine treatments? Use $\\alpha$=0.05 for NHST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0006"
      ]
     },
     "execution_count": 85,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(sims >= 1.5)\n",
    "pval = 6/10000\n",
    "pval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## Post Hoc Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Having obtained a significant result from the omnibus test, we can now try to track down the source of the significance. Which groups are actually different from which?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "16. Perform pairwise two-group comparisons and record the p-value for each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7fe751cd13d0>"
      ]
     },
     "execution_count": 90,
     "metadata": {
     },
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAALACAYAAABy7a/RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAABYlAAAWJQFJUiTwAAA37klEQVR4nO3df7xmZV0v/M9XUMxBRgUSlToIgfpkHkNNxVLUJ19qZqb4kjopoWKaSCJ06uAvNC17QEnE5DloTIk9aHjUgyl1jjghjWmChh5TEBx/BeqI7OGHoMD1/LHWzt1m75m9Z+7Z9+y53u/X635dc1/r+q513TfDns+sWWtd1VoLAAD05E7TngAAAKw0IRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtC8CpRVe+pqvdMex4AALuC3ac9AZbsgYceeuihSX5z2hMBANiJ1LYUORMMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALqz6kJwVf1pVX2sqr5RVT+oqmur6rNV9dqq2nve2AOqqm3hde4WjnNUVX26qm6oqpmqWl9VT9vC+N2q6viqumzOvD5SVYdN8vMDALD9dp/2BLbB8UkuTfK/knwnyZokj0pycpIXVdWjWmvfmFfzL0k+uMC+vrDQAarq1CQnJPlmkrOS3CXJkUnOr6qXtdbOmDe+kpyb5IgkX05yRpJ7JXlOkouq6lmttQ8t+5MCALBDrMYQvFdr7eb5nVX1xiQnJflvSX533ubPtdZOXsrOxzO3JyS5MskjWmvfH/tPSXJJklOr6sOttY1zyo7MEIA3JHni7Pyq6swkFyc5q6oubK1dv+RPCbAT27RpU2ZmZpY8fu3atdlnn3124IwAlmfVheCFAvDofRlC8MHbeYgXj+0bZwPweNyNVfX2JK9OcnSS186pecnYvmru/Fpr/1xV703y3Awh+eztnBvA1G3atCkHHnRQrt+8eck1d99rr1x15ZWCMLDTWHUheAt+dWwvW2Dbfavqd5LsneR7ST7ZWltoXJI8YWwvWGDbRzOE4CdkDMFVddckhyW5KcknFql57lgjBAOr3szMTK7fvDmHH3961uy931bH3/i9a7L+tOMyMzMjBAM7jVUbgqvqxCR7Jlmb5OFJfjFDAH7TAsN/eXzNrV+f5KjW2tfn9K1Jcr8kN7TWrl5gP1eM7SFz+g5KsluSq1prty6xZlFVdckimx64lHqAlbJm7/2y5777T3saANtk1YbgJCcmufec9xck+e3W2nfn9N2U5I8y3BR31dj3kAw30T0+yceq6qGttRvHbWvHdrEL3Wb77zGnb1tqAACYolUbgltr+yVJVd07w+UIb0ry2ap6Wmvt0nHMd5K8Zl7pRVX1pAw3rD0yyQuTvHXFJr4VrbWHLdQ/niE+dIWnAwCwS1p1zwmer7X27dbaB5I8KcM1v3+1hJpbk7xzfPvYOZtmz9quzcJm+6/bzhoAAKZo1YfgWa21ryX5YpKfraql3Hkxe9nEmjn7uDHJt5LsWVX3WaBm9skTl8/puzLJbUkOrKqFzqwvVAMAwBTtMiF4dN+xvW0JYx81tlfN679wbJ+8QM1T5o2ZfWTbhiR3S/JLS6kBAGC6VlUIrqpDquoOlx1U1Z3GxTJ+MsmGOQtcHFpVd/iMVfXEDCvPJck58zafObavrKp7zqk5IMlLk9ySOz7q7B1j+4bxkWmzNY/IsGrcd5O8f0kfEgCAHW613Rj31CR/UlUXJ/lqhmf+3jvJ45IcmOSaJMfMGf+WJAdX1YYMSyAnw9MhZp8F/OrW2oa5B2itbaiqtyR5RZLLquq8DMsmPyfDUsgvm7daXDIsmfzMDAtifLaqzs9wffJzMjw+7ZjW2tKfKg8AwA612kLw/07yMxmeCfzzGR47dmOG623fneT01tq1c8a/O8mvJ3lEhssS7pzk2xlWlzujtbbQ4hZprZ1QVZ/PcOb3RUluT3JpklNaax9eYHyrqt/IcFnE85O8LMnNSS5K8ob5QRsAgOlaVSG4tfaFJMcuY/y7krxrG4+1Lsm6ZYy/Nclp4wsAgJ3YqromGAAAJkEIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO6suhBcVX9aVR+rqm9U1Q+q6tqq+mxVvbaq9l6k5rCq+sg49gdVdVlVvbyqdtvCcZ5WVeuraqaqbqiqT1XVUVuZ21FV9elx/MxY/7Tt/cwAAEzWqgvBSY5PsibJ/0ry1iTvSXJrkpOTXFZVPzV3cFX9WpKLkjw2yQeSnJHkLklOS3LuQgeoqmOTnJ/kwUnOSXJWkvsmWVdVpy5Sc2qSdUnuM44/J8nPJTl/3B8AADuJ3ac9gW2wV2vt5vmdVfXGJCcl+W9Jfnfs2ytDIL0tyeGttc+M/a9OcmGSI6rqyNbauXP2c0CSU5Ncm+ThrbWNY//rk/xzkhOq6v2ttU/OqTksyQlJrkzyiNba98f+U5JckuTUqvrw7L4AAJiuVXcmeKEAPHrf2B48p++IJPsmOXc2AM/Zx6vGty+Zt5/nJ9kjyRlzQ+sYbP94fPvieTWz7984G4DHmo1J3j7u7+hFPxQAACtq1YXgLfjVsb1sTt8TxvaCBcZflOSmJIdV1R5LrPnovDHbUwMAwJSsxsshkiRVdWKSPZOsTfLwJL+YIQC/ac6wB4zt5fPrW2u3VtVXk/xskgOT/OsSaq6uqhuT7F9Vd2ut3VRVa5LcL8kNrbWrF5jqFWN7yBI/1yWLbHrgUuoBANi6VRuCk5yY5N5z3l+Q5Ldba9+d07d2bGcW2cds/z2WWbNmHHfTNh4DAIApWrUhuLW2X5JU1b2THJbhDPBnq+pprbVLpzq57dBae9hC/eMZ4kNXeDoAALukVX9NcGvt2621DyR5UpK9k/zVnM2zZ2HX3qHwP/Zftw01M/Pa5RwDAIApWvUheFZr7WtJvpjkZ6tqn7H7y2N7h+txq2r3JPfP8Izhq+Zs2lLNfTJcCvHN1tpN43FvTPKtJHuO2+ebfVrFHa4xBgBgOnaZEDy679jeNrYXju2TFxj72CR3S7KhtXbLnP4t1Txl3pjtqQEAYEpWVQiuqkOq6g6XHVTVncbFMn4yQ6idfVbveUk2JTmyqh4+Z/xdk7xhfPuOebs7O8ktSY4dF86YrblnhsU4kuTMeTWz7185jputOSDJS8f9nb3EjwkAwA622m6Me2qSP6mqi5N8Ncn3Mjwh4nEZHnN2TZJjZge31jZX1TEZwvD6qjo3w0pwT8/wKLTzkrx37gFaa1+tqt9PcnqSz1TVe5P8MMPCG/snefPc1eLGmg1V9ZYkr8iwdPN5GZZmfk6SeyV5mdXiAAB2HqstBP/vJD+T4ZnAP5/hsWM3Zrje9t1JTm+tXTu3oLX2wap6XJJXJnlWkrsm+UqGwHp6a63NP0hr7W1VtTHDY9iel+GM+ReTvKq19pcLTay1dkJVfT7Dmd8XJbk9yaVJTmmtfXj7PjYAAJO0qkJwa+0LSY7dhrp/zHAWeTk15yc5f5k165KsW04NAAArb1VdEwwAAJMgBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3VlUIrqq9q+qFVfWBqvpKVf2gqmaq6uKqekFV3Wne+AOqqm3hde4WjnVUVX26qm4Yj7G+qp62hfG7VdXxVXXZOK9rq+ojVXXYJL8DAAC23+7TnsAyPTvJO5JcneTjSb6e5N5JnpnknUmeUlXPbq21eXX/kuSDC+zvCwsdpKpOTXJCkm8mOSvJXZIcmeT8qnpZa+2MeeMryblJjkjy5SRnJLlXkuckuaiqntVa+9CyPy0AADvEagvBlyd5epK/ba3dPttZVScl+XSSZ2UIxO+fV/e51trJSznAeOb2hCRXJnlEa+37Y/8pSS5JcmpVfbi1tnFO2ZEZAvCGJE9srd081pyZ5OIkZ1XVha2165f3cQEA2BFW1eUQrbULW2vnzw3AY/81Sc4c3x6+nYd58di+cTYAj8fYmOTtSfZIcvS8mpeM7atmA/BY889J3ptk3wwhGQCAncCqCsFb8aOxvXWBbfetqt+pqpPG9iFb2M8TxvaCBbZ9dN6YVNVdkxyW5KYkn1hKDQAA07XaLodYUFXtnuR549uFwusvj6+5NeuTHNVa+/qcvjVJ7pfkhtba1Qvs54qxPWRO30FJdktyVWttoQC+UM2iquqSRTY9cCn1AABs3a5yJvhNSR6c5COttb+b039Tkj9K8rAk9xxfj8twU93hST42Bt9Za8d2ZpHjzPbfYztrAACYolV/JriqjstwI9uXkjx37rbW2neSvGZeyUVV9aQMN6w9MskLk7x1Baa6JK21hy3UP54hPnSFpwMAsEta1WeCq+rYDAH2i0ke31q7dil142UL7xzfPnbOptmztmuzsNn+67azBgCAKVq1IbiqXp7kbRme9fv48QkRy/Hdsf33yyFaazcm+VaSPavqPgvUHDy2l8/puzLJbUkOHK9NXkoNAABTtCpDcFX9QZLTknwuQwD+zjbs5lFje9W8/gvH9skL1Dxl3piMj0TbkORuSX5pKTUAAEzXqgvBVfXqDDfCXZJhYYpNWxh76PyllMf+JyY5fnx7zrzNs88bfmVV3XNOzQFJXprkliRnz6t5x9i+YXxk2mzNIzKsGvfd3HEBDwAApmRV3RhXVUcleX2Gyw8+keS4YcXi/2Bja23d+Ou3JDm4qjZkWAI5SR6SHz+z99WttQ1zi1trG6rqLUlekeSyqjovw7LJz8mwFPLL5q0WlwxLJj8zw4IYn62q85PsPdbsluSY1trmbf3cAABM1qoKwUnuP7a7JXn5ImP+Icm68dfvTvLrSR6R4bKEOyf5dpL3JTmjtbbQ4hZprZ1QVZ/PcOb3RUluT3JpklNaax9eYHyrqt/IcFnE85O8LMnNSS5K8ob5QRsAgOlaVSG4tXZykpOXMf5dSd61jcdalx+H6aWMvzXDdcqnbcvxAABYOavummAAANheQjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdGf3aU8AAHZ2mzZtyszMzJLHr127Nvvss88OnBGwvYRgANiCTZs25cCDDsr1mzcvuebue+2Vq668UhCGnZgQDABbMDMzk+s3b87hx5+eNXvvt9XxN37vmqw/7bjMzMwIwbATE4IBYAnW7L1f9tx3/2lPA5gQN8YBANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSvGAcyzadOmzMzMLKtm7dq1lsgFWEWEYIA5Nm3alAMPOijXb968rLq777VXrrrySkEYYJUQggHmmJmZyfWbN+fw40/Pmr33W1LNjd+7JutPOy4zMzNCMMAqIQQDLGDN3vtlz333n/Y0ANhB3BgHAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3VlVIbiq9q6qF1bVB6rqK1X1g6qaqaqLq+oFVbXg56mqw6rqI1V17VhzWVW9vKp228KxnlZV68f931BVn6qqo7Yyv6Oq6tPj+Jmx/mnb+7kBAJisVRWCkzw7yVlJHpnkU0n+LMn7kzw4yTuTvK+qam5BVf1akouSPDbJB5KckeQuSU5Lcu5CB6mqY5OcP+73nPGY902yrqpOXaTm1CTrktxnHH9Okp9Lcv64PwAAdhK7T3sCy3R5kqcn+dvW2u2znVV1UpJPJ3lWkmdmCMapqr0yBNLbkhzeWvvM2P/qJBcmOaKqjmytnTtnXwckOTXJtUke3lrbOPa/Psk/Jzmhqt7fWvvknJrDkpyQ5Mokj2itfX/sPyXJJUlOraoPz+4LAIDpWlVngltrF7bWzp8bgMf+a5KcOb49fM6mI5Lsm+Tc2QA8jr85yavGty+Zd5jnJ9kjyRlzQ+sYbP94fPvieTWz7984G4DHmo1J3j7u7+itf0IAAFbCqgrBW/Gjsb11Tt8TxvaCBcZflOSmJIdV1R5LrPnovDHbUwMAwJSstsshFlRVuyd53vh2bhB9wNhePr+mtXZrVX01yc8mOTDJvy6h5uqqujHJ/lV1t9baTVW1Jsn9ktzQWrt6geldMbaHLPGzXLLIpgcupR4AgK3bVc4EvynDTWwfaa393Zz+tWM7s0jdbP89tqFm7bx2OccAAGCKJnomuKp+Osl1rbXNWxhz9yT3bK19fULHPC7DTWlfSvLcSexzmlprD1uofzxDfOgKTwcAYJc06TPBX03ye1sZc9w4bruNjx57a5IvJnl8a+3aeUPmn7Wdb7b/um2omZnXLucYAABM0aRDcI2vHa6qXp7kbUm+kCEAX7PAsC+P7R2uxx2vI75/hhvprlpizX2SrEnyzdbaTUnSWrsxybeS7Dlun+/gsb3DNcYAAEzHNK4J3i/Jjduzg6r6gwyLXXwuQwD+ziJDLxzbJy+w7bFJ7pZkQ2vtliXWPGXemO2pAQBgSrb7muCqet68rocu0JckuyX56SS/leTz23G8Vyd5fYZFKJ60wCUQc52X5E+THFlVb5uzWMZdk7xhHPOOeTVnJ/mvSY6tqrPnLJZxzyQnjWPOnFdzZobrkV9ZVR+cs1jGAUlemuSWcb8AAOwEJnFj3Lokbfx1S/Jr42u+2cskbkryum05UFUdlSEA35bkE0mOm7dKcpJsbK2tS5LW2uaqOiZDGF5fVedmWAnu6RkehXZekvfOLW6tfbWqfj/J6Uk+U1XvTfLDDAtv7J/kzXNXixtrNlTVW5K8IsllVXVehqWZn5PkXkleZrU4AICdxyRC8OxKaJXkL5J8MMmHFhh3W5LvJflka+26bTzW/cd2tyQvX2TMP2QI5kmS1toHq+pxSV6ZYVnluyb5SobAenprrc3fQWvtbVW1McmJGZ4/fKcMN9+9qrX2lwsdtLV2QlV9PsOZ3xcluT3JpUlOaa19eFmfEgCAHWq7Q/DcUDieqf1ga+2vtne/ixzr5CQnb0PdPyZ56jJrzk9y/jJr1mVOAAcAYOc00ecEt9YeP8n9AQDAjrCrrBgHAABLNvEQXFWPq6oPV9V3qupHVXXbAq9bJ31cAABYqkkvm/wrGW6M2y3J1zMsPCHwAgCwU5loCM5w09qPkvxKa+3vJ7xvAACYiElfDvHgJO8VgAEA2JlNOgTfkGExCgAA2GlNOgR/LMmjJ7xPAACYqEmH4D9IclBVvaoWWM8YAAB2BpO+Me61Sf5PktcleX5VfS7JdQuMa621F0z42AAAsCSTDsG/PefXB4yvhbQkQjAAAFMx6RB8/wnvDwAAJm6iIbi19rVJ7g8AAHaEiS+bDAAAO7tJL5v800sd21r7+iSPDQAASzXpa4I3ZrjpbWvaDjg2AAAsyaSD6F9l4RB8jyQPTfKfkqxP4tphAACmZtI3xv32Ytuq6k5JXp3kxUmOmuRxAQBgOVbsxrjW2u2ttddluGTiTSt1XAAAmG8aT4fYkORJUzguAAAkmU4IvleSNVM4LgAAJFnhEFxV/3eS5yT5wkoeFwAA5pr0c4Iv3MJxfirJ7HOEXz/J4wIAwHJM+hFphy/S35J8P8nfJTm1tbZYWAYAgB1u0o9IswwzAAA7PaEVAIDu7NCli6vq7hlWi5tprW3ekccCAIClmviZ4Kravar+sKq+kuS6DItjfL+qvjL279DgDQAAWzPpp0PcJckFSR6X4Wa4byS5Osl9khyQ5I1JnlxVT2qt/XCSxwYAgKWa9JngV2R4QsTfJnlQa+2A1tqjW2sHJHlAkvOT/NI4DgAApmLSIfg3MyyE8YzW2hVzN7TWrkzyzCT/J8l/mfBxAQBgySYdgn8myUdba7cvtHHs/2iSgyZ8XAAAWLJJh+AfJtlzK2PWJPnRhI8LAABLNukQfFmSI6pq34U2VtU+SY5I8i8TPi4AACzZpEPwGUn2TfLpqnpBVR1YVT9RVfevqqOTfGrcfsaEjwsAAEs26WWT31dVD03yh0n++wJDKsn/01p73ySPCwAAyzHxhStaaydV1f9M8oIkP59kbZKZJJ9N8hettU9O+pgAALAcO2T1ttbaPyX5px2xbwAA2F7bfU1wVd2lqj5dVR+rqjtvZdzHquqftjQOAAB2tEncGPdbSR6W5M2ttUUffTYuk3xKkl+IxTIAAJiiSYTgZya5qrX2ka0NbK1dkOSKJM+ewHEBAGCbTCIE/3yS9csYf1GSh07guAAAsE0mEYL3SfLtZYz/dpK9J3BcAADYJpMIwT/I1pdKnmvPJDdP4LgAALBNJhGCv5Hk4csY//AkX5/AcQEAYJtMIgSvT/LoqtpqEK6qhyU5LMnHJ3BcAADYJpMIwWckaUn+pqoetNigqnpgkr9JcluSP5/AcQEAYJts94pxrbUvV9Xrk5yc5LNVdV6SC5N8cxxyvyRPTPKsJHskeU1r7cvbe1wAANhWE1k2ubX2+qq6Nclrk/xmkt+YN6SS/CjJK1trfzKJYwIAwLaaSAhOktbaH1fVe5I8P8ljktxn3HR1kouTnN1a+9qkjgcAANtqYiE4ScaQ+9pJ7hMAACZtEjfGAQDAqiIEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO6suBFfVEVX1tqr6RFVtrqpWVecsMvaAcftir3O3cJyjqurTVXVDVc1U1fqqetoWxu9WVcdX1WVV9YOquraqPlJVh03icwMAMDm7T3sC2+BVSf5zkhuSfDPJA5dQ8y9JPrhA/xcWGlxVpyY5Ydz/WUnukuTIJOdX1ctaa2fMG19Jzk1yRJIvJzkjyb2SPCfJRVX1rNbah5YwTwAAVsBqDMHHZwinX0nyuCQfX0LN51prJy9l5+OZ2xOSXJnkEa2174/9pyS5JMmpVfXh1trGOWVHZgjAG5I8sbV281hzZpKLk5xVVRe21q5fyhwAANixVt3lEK21j7fWrmittR10iBeP7RtnA/B43I1J3p5kjyRHz6t5ydi+ajYAjzX/nOS9SfbNEJIBANgJrLoQvI3uW1W/U1Unje1DtjD2CWN7wQLbPjpvTKrqrkkOS3JTkk8spQYAgOlajZdDbItfHl//rqrWJzmqtfb1OX1rktwvyQ2ttasX2M8VY3vInL6DkuyW5KrW2q1LrFlUVV2yyKalXPsMAMAS7Opngm9K8kdJHpbknuNr9jriw5N8bAy+s9aO7cwi+5vtv8d21gAAMEW79Jng1tp3krxmXvdFVfWkDDesPTLJC5O8daXntpjW2sMW6h/PEB+6wtMBANgl7epnghc0XrbwzvHtY+dsmj1ruzYLm+2/bjtrAACYoi5D8Oi7Y/vvl0O01m5M8q0ke1bVfRaoOXhsL5/Td2WS25IcWFULnVlfqAYAgCnqOQQ/amyvmtd/4dg+eYGap8wbk/GRaBuS3C3JLy2lBgCA6dqlQ3BVHVpVd/iMVfXEDItuJMn8JZfPHNtXVtU959QckOSlSW5Jcva8mneM7RvGR6bN1jwiw6px303y/m38GAAATNiquzGuqp6R5Bnj2/3G9tFVtW789abW2onjr9+S5OCq2pBhlbkkeUh+/MzeV7fWNszdf2ttQ1W9JckrklxWVedlWDb5ORmWQn7ZvNXikmHJ5GdmWBDjs1V1fpK9x5rdkhzTWtu8rZ8ZAIDJWnUhOMlDkxw1r+/A8ZUkX0syG4LfneTXkzwiw2UJd07y7STvS3JGa22hxS3SWjuhqj6f4czvi5LcnuTSJKe01j68wPhWVb+R4bKI5yd5WZKbk1yU5A3zgzYAANO16kJwa+3kJCcvcey7krxrG4+zLsm6ZYy/Nclp4wsAgJ3YLn1NMAAALEQIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO6suhBcVUdU1duq6hNVtbmqWlWds5Waw6rqI1V1bVX9oKouq6qXV9VuW6h5WlWtr6qZqrqhqj5VVUdt5ThHVdWnx/EzY/3TtvWzAgCwY6y6EJzkVUmOTfLQJN/a2uCq+rUkFyV5bJIPJDkjyV2SnJbk3EVqjk1yfpIHJzknyVlJ7ptkXVWdukjNqUnWJbnPOP6cJD+X5PxxfwAA7CRWYwg+PskhSfZK8pItDayqvTIE0tuSHN5ae0Fr7fczBOhPJjmiqo6cV3NAklOTXJvk4a21l7bWjk/ykCRXJjmhqh49r+awJCeM2x/SWju+tfbSJA8b93PquF8AAHYCqy4Et9Y+3lq7orXWljD8iCT7Jjm3tfaZOfu4OcMZ5eSOQfr5SfZIckZrbeOcmu8n+ePx7Yvn1cy+f+M4brZmY5K3j/s7egnzBQBgBay6ELxMTxjbCxbYdlGSm5IcVlV7LLHmo/PGbE8NAABTsvu0J7CDPWBsL5+/obV2a1V9NcnPJjkwyb8uoebqqroxyf5VdbfW2k1VtSbJ/ZLc0Fq7eoE5XDG2hyxlwlV1ySKbHriUegAAtm5XPxO8dmxnFtk+23+PbahZO69dzjEAAJiiXf1M8KrTWnvYQv3jGeJDV3g6AAC7pF39TPD8s7bzzfZftw01M/Pa5RwDAIAp2tVD8JfH9g7X41bV7knun+TWJFctseY+SdYk+WZr7aYkaa3dmOF5xXuO2+c7eGzvcI0xAADTsauH4AvH9skLbHtskrsl2dBau2WJNU+ZN2Z7agAAmJJdPQSfl2RTkiOr6uGznVV11yRvGN++Y17N2UluSXLs3AUuquqeSU4a3545r2b2/SvHcbM1ByR56bi/s7fngwAAMDmr7sa4qnpGkmeMb/cb20dX1brx15taaycmSWttc1UdkyEMr6+qczOs4Pb0DI9COy/Je+fuv7X21ar6/SSnJ/lMVb03yQ8zLLyxf5I3t9Y+Oa9mQ1W9JckrklxWVedlWJr5OUnuleRlcxfeAABgulZdCM6w5PFR8/oOHF9J8rUkJ85uaK19sKoel+SVSZ6V5K5JvpIhsJ6+0MpzrbW3VdXGcT/Py3DG/ItJXtVa+8uFJtVaO6GqPp/hzO+Lktye5NIkp7TWPrxNnxQAgB1i1YXg1trJSU5eZs0/JnnqMmvOT3L+MmvWJVm3nBoAAFbern5NMAAA3IEQDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALqz+7QnAOz8Nm3alJmZmSWPX7t2bfbZZ58dOCMA2D5CMLBFmzZtyoEHHZTrN29ecs3d99orV115pSAMwE5LCAa2aGZmJtdv3pzDjz89a/beb6vjb/zeNVl/2nGZmZkRggHYaQnBwJKs2Xu/7Lnv/tOeBgBMhBAMAGw39w6w2gjBAMB2ce8Aq5EQDABsF/cOsBoJwQDARLh3gNXEYhkAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAutNFCK6qjVXVFnlds0jNYVX1kaq6tqp+UFWXVdXLq2q3LRznaVW1vqpmquqGqvpUVR214z4ZAADbYvdpT2AFzST5swX6b5jfUVW/luT9SW5O8t4k1yb51SSnJXlMkmcvUHNskrcl+V6Sc5L8MMkRSdZV1c+11k6cyKcAAGC79RSCr2utnby1QVW1V5KzktyW5PDW2mfG/lcnuTDJEVV1ZGvt3Dk1ByQ5NUNYfnhrbePY//ok/5zkhKp6f2vtkxP9RAAAbJMuLodYpiOS7Jvk3NkAnCSttZuTvGp8+5J5Nc9PskeSM2YD8Fjz/SR/PL598Y6aMAAAy9PTmeA9quq3kvx0khuTXJbkotbabfPGPWFsL1hgHxcluSnJYVW1R2vtliXUfHTemC2qqksW2fTApdQDALB1PYXg/ZK8e17fV6vq6NbaP8zpe8DYXj5/B621W6vqq0l+NsmBSf51CTVXV9WNSfavqru11m7ang8BAMD26yUEn53kE0n+T5LrMwTYY5O8KMlHq+rRrbV/GceuHduZRfY123+POX1LqVkzjttiCG6tPWyh/vEM8aFbqgUAYGm6CMGttdfN6/pCkhdX1Q1JTkhycpJfX+l5AQAwHb3fGHfm2D52Tt/s2dy1Wdhs/3XbULPYmWIAAFZQ7yH4u2O7Zk7fl8f2kPmDq2r3JPdPcmuSq5ZYc59x/990PTAAwM6h9xD8qLGdG2gvHNsnLzD+sUnulmTDnCdDbK3mKfPGAAAwZbt8CK6qB1XVmgX6D0hyxvj2nDmbzkuyKcmRVfXwOePvmuQN49t3zNvd2UluSXLsuN/ZmnsmOWl8e2YAANgp9HBj3HMyrNh2UZKvZXg6xEFJfiXJXZN8JMNqb0mS1trmqjomQxheX1XnZlgJ7ukZHoV2XoallDOn5qtV9ftJTk/ymap6b368bPL+Sd5stTgAgJ1HDyH44xnC688neUyG63OvS3JxhucGv7u11uYWtNY+WFWPS/LKJM/KEJa/kuQVSU6fP36seVtVbUxyYpLnZTjL/sUkr2qt/eUO+WQAAGyTXT4Ejwth/MNWB96x7h+TPHWZNecnOX+5xwIAYGXt8tcEAwDAfEIwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0Z/dpTwB2Bps2bcrMzMySx69duzb77LPPDpwRALAjCcF0b9OmTTnwoINy/ebNS665+1575aorrxSEAWCVEoLp3szMTK7fvDmHH3961uy931bH3/i9a7L+tOMyMzMjBAPAKiUEw2jN3vtlz333n/Y0AIAVIAQDACyTe0lWPyEYAGAZ3EuyaxCCAQCWwb0kuwYhGABgG7iXZHWzWAYAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QvAEVdX+VfUXVfVvVXVLVW2sqj+rqntOe24AAPzY7tOewK6iqg5KsiHJTyb5UJIvJfmFJL+X5MlV9ZjW2vemOEUAAEbOBE/On2cIwMe11p7RWvvD1toTkpyW5AFJ3jjV2QEA8O+E4AkYzwI/KcnGJG+ft/m1SW5M8tyqWrPCUwMAYAEuh5iMx4/t37fWbp+7obV2fVX9Y4aQ/KgkH1vpyS3Xpk2bMjMzs+Txa9euzT777LMDZwQATNNys0Gy8+eDaq1New6rXlWdkuTEJCe21t68wPYzkrw0ye+21t6xlX1dssim//wTP/ETuz3oQQ/a7vluya233povfOELuf3227c+eHSnO90phxxySHbbbbcdOLMd54c//GGuuOKK7PmT++dOu995q+Nvv/VHueE738zBBx+cu9zlLisww+nq7ftZ7udNVv9nXi6/J7ZstX/ebdHbd9Tb573tttty+eWXLysbJEM+ePCDH5zdd9+x51wvvfTSv26t/Zfl1gnBE1BV/z3JMUmOaa29c4Htb0xyUpKTWmt/spV9LRaCH5zkhgyXXOxIDxzbL+3g46wGvouB7+HHfBcD38PA9/BjvouB7+HHVvK7+NK2hGCXQ+xkWmsPm+bxZ0P4tOexM/BdDHwPP+a7GPgeBr6HH/NdDHwPP7Yavgs3xk3G7EUyaxfZPtt/3Y6fCgAAWyMET8aXx/aQRbYfPLaXr8BcAADYCiF4Mj4+tk+qqv/wnVbV3ZM8JslNSf5ppScGAMAdCcET0Fq7MsnfJzkgw1Mg5npdkjVJ3t1au3GFpwYAwALcGDc5v5th2eTTq+qJSf41ySMzPEP48iSvnOLcAACYwyPSJqiqfirJ65M8OcneSa5O8oEkr2utfX+acwMA4MeEYAAAuuOaYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCWbaqemdVtfH1M9Oez0qoqp+qqj+vqk9V1TVVdUtV/VtVfaKqjq6qO097jiulqg6uqj+oqgur6htV9cOq+nZVfaiqHj/t+a2UqrpzVf1eVZ1dVZ8bv4dWVS+c9tx2lKrav6r+Yvy9f0tVbayqP6uqe057biulqo6oqreN/+9vHv+bnzPtea20qtq7ql5YVR+oqq9U1Q+qaqaqLq6qF1RVN/miqv60qj42/jz8QVVdW1WfrarXVtXe057fNFXVb83JCzvdz0bPCWZZqupXk/zPJDck2TPJwa21r0x3VjteVR2e5ENJPpXkqiTXZlgQ5SlJfirJx5M8qbV265SmuGKq6twkz0nyxSQXZ/guHpDk6Ul2S/J7rbXTpzfDlVFV90gyuwjOt5P8MMPvhWNaa++c1rx2lKo6KMOqmD+Z4f+FLyX5hQyrYn45yWNaa9+b3gxXRlV9Lsl/zvAz8JtJHpjkPa2135rmvFZaVb04yTsyLAr18SRfT3LvJM9MsjbJ+5M8u3UQMqrqh0kuzfAz8TtJ1iR5VJKHJ/m3JI9qrX1jejOcjnEBsc9n+HNhz+yEPxstm8ySVdW+Sc5K8t4k+yV53HRntKI2JLlna+32uZ3jGeC/zxAEnpnkfVOY20q7IMmfttY+O7ezqh6X5H8lOaWq/qa1dvVUZrdybkry1CSfa61dXVUnJ3ntdKe0Q/15hgB8XGvtbbOdVfWWJMcneWOSF09pbivp+Azh9ysZfgZ+fLrTmZrLM/zF92/n/lysqpOSfDrJszL8THz/dKa3ovZqrd08v7Oq3pjkpCT/LcnvrvispqiqKsnZSb6X5H8kOXG6M1pYN/9cwUT897F96VRnMQWttR/OD8Bj/4+SfHB8e/CKTmpKWmvr5gfgsf8fkqxPcpckh630vFba+Hviox2E/dmzwE9KsjHJ2+dtfm2SG5M8t6rWrPDUVlxr7eOttSt6OMO5Ja21C1tr58//udhauybJmePbw1d8YlOwUAAezZ4U6eLPhnmOS/KEJEdn+PmwUxKCWZKq+u0kz0jyOz38k+dSVdVuGc4GJsll05zLTuJHY7vLXxbSmdlrvf9+gdBzfZJ/THK3DP8EDH4ODH51bLv6s6GqHpTkTUne2lq7aNrz2RKXQ7BVVfWfkrw1yTmttQ9Nez7TVFX7JDk2SSXZN8kvJ/mZJH/dWjt/mnObtvH3yRMzXCawU//gY9keMLaXL7L9igxnig9J8rEVmRE7paraPcnzxrcXTHMuK62qTsxw7evaDNcD/2KGAPymac5rJY3//d+d4Rrxk6Y8na0Sgtmi8Q7fv8xwE8hxU57OzmCf/MfrPluSU7MK/mffkapqjyTvSbJHkv/aWvv+VkpYXdaO7cwi22f777Hjp8JO7k1JHpzkI621v5v2ZFbYiRluDpx1QZLfbq19d0rzmYbXJPn5JL/YWvvBtCezNS6H6MD4GKO2jNfcx/0cn+Hmj2NWe7DZzu8hSdJa+1JrrTL8BfI/Zfh+XpTkoqq61wp/pG02ie9izr52y/A3/8dkuGny1JX6HNtrkt8D9K6qjktyQoYnhzx3ytNZca21/cY/H/bLcFPggUk+W1WHTndmK6OqHpnhhNCbW2ufnPZ8lsKZ4D5cmWSxC/cX8m9JUlWHZLjj++zW2kd2xMRW2DZ9Dwtprd2W4Z973lpV307y/yV5fYZLJVaDiXwXYwA+J8mzM9wE8lur7Iahif2e2MXNnuldu8j22f7rdvxU2BlV1bEZLpv7YpInttaunfKUpqa19u0kH6iqSzNcQvRXGc6O77LGyyD+KsPnffWUp7NkQnAHWmtP3MbS/yvDP28fXVVHLzLmiuFJKPn11toHt/E4K2I7voet+ejYHr6D9j9xk/guxsfDvSdDAP7rJM8b/3KwauzA3xO7mi+P7SGLbJ+9+32xa4bZhVXVy5OcluQLGQLwd6Y7o51Da+1rVfXFJA+tqn1aa5umPacdaM/8+OfDzWMumO+sqjorww1zL1+piW2JEMyWbEzyrkW2/UqGf/L5mySbx7G9ut/YdnMndFXdJcOZ31/L8Lf/oxd6hBy7jNln4T6pqu4077mwd89wKcxNSf5pGpNjeqrqDzJcB/y5JL+8iwe9bXHfsV1VJwi2wS1ZPC8cmuE64Ysz/IV6p7lUQghmUa21zyVZcJnDqlqfIQSf1MmKcYcm+Zf5Zzqras8M/wSYJH+74hObgvEmuP+R4dFw70ryIgF419Zau7Kq/j7DEyBemuRtcza/LsMKWf9va22nfR4ok1dVr85wGdglGVbM7O4SiPGywW+31mbm9d8pyR9lWGBmw2q/p2ZrxpvgFssLJ2cIwX9pxThYnV6T5DFVtSHDtcA3ZVgi9ykZ7ojfkORPpja7lXVmhgC8Kcm3krxmgX/6Wt9aW7/C81pxVfWHGZbNTZKHju3RVfWL468v3tl+6G+H383w+/z0qnpikn9N8sgMzxC+PMkrpzi3FVNVz8jwzPRkOBGQJI+uqnXjrze11nbK1bEmqaqOyhCAb0vyiSTHLfBzYGNrbd0KT22lPTXJn1TVxUm+mmGFtHtnuKH8wCTXJDlmetNjS4RgWJqzMjwm7hcyXPt7tyTfz3AG5H1J/qK11svlEPcf230y/OVgMet3/FSm7sm54/Lhh+U/rpi3S4Tg8WzwwzMEnydn+MP/6gz/EvK6Xf1M1xwPTXLUvL4Dx1eSfC076RKxEzb7c2C3JC9fZMw/JFm3EpOZov+d4Vnxv5jhbOc9MqyQdnmGp+ac3uMZ8tWiVteN3AAAsP08JxgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7vz/Z37ToTPds4YAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "execution_count": 90,
     "metadata": {
      "image/png": {
       "height": 352,
       "width": 352
      },
      "needs_background": "light"
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pop = np.concatenate([A,B])\n",
    "\n",
    "store = np.zeros(10000)\n",
    "for i in range(10000):\n",
    "    Asim = np.random.choice(pop, len(migraines))\n",
    "    Bsim = np.random.choice(pop, len(migraines))\n",
    "    Amed = np.median(Asim)\n",
    "    Bmed = np.median(Bsim)\n",
    "    difMed = Amed - Bmed\n",
    "    store[i] = difMed\n",
    "\n",
    "sns.displot(store)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0407"
      ]
     },
     "execution_count": 100,
     "metadata": {
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "meds = migraines.median()\n",
    "meds\n",
    "diffMed = -3\n",
    "\n",
    "(np.sum(store >= 3) + np.sum(store <= -3)) / 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "17. Use the list of p-values from your two-group comparisons and the Benjamini-Hochberg method below to determine which two-group comparisons have a significant difference at $\\alpha$=0.05. See more details in \"Controlling the False Positive Rate with Multiple Comparisons Corrections\" PDF on Week 6 on CCLE site. \\\n",
    "\\\n",
    "`import statsmodels.stats.multitest as smm` \\\n",
    "`smm.multipletests(pvals, alpha=my_alpha, method='fdr_bh')`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-95-d55528c377ca>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-95-d55528c377ca>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    np.sum(store => 3)\u001b[0m\n\u001b[0m                  ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "18. Write a sentence or two describing your findings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (system-wide)",
   "language": "python",
   "metadata": {
    "cocalc": {
     "description": "Python 3 programming language",
     "priority": 100,
     "url": "https://www.python.org/"
    }
   },
   "name": "python3",
   "resource_dir": "/ext/jupyter/kernels/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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}