Two Sample Hypothesis Testing by Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two-Sample Hypothesis Testing by Simulation\n",
    "(c) 2018 by Jordan Ryda  \n",
    "__[Framework](http://allendowney.blogspot.co.uk/2016/06/there-is-still-only-one-test.html)__ inspired by Allen Downey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set(color_codes=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_stat(group_a, group_b):\n",
    "    return(np.mean(group_a) - np.mean(group_b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_from_groups(group_a, group_b, replace=True):\n",
    "    merged_group = np.append(group_a, group_b)\n",
    "    \n",
    "    sampled_merged_group = np.random.choice(merged_group, len(merged_group), replace = replace)\n",
    "    sim_group_a = sampled_merged_group[:len(group_a)]\n",
    "    sim_group_b = sampled_merged_group[len(group_a):]\n",
    "    \n",
    "    return(sim_group_a, sim_group_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 tail | H1: A != B [DEFAULT]\n",
    "# 1 tail | H1: A > B\n",
    "\n",
    "def run_hyp_test(group_a, group_b, num_sims, replace=True, alpha=.01, tails=2):\n",
    "    \n",
    "    if tails not in [1,2]:\n",
    "        print(\"tails can only be 1 or 2\")\n",
    "        return\n",
    "\n",
    "    delta_dist_h0 = []\n",
    "    delta_star = test_stat(group_a, group_b)\n",
    "\n",
    "    # Run Simulations\n",
    "    for i in range(num_sims):\n",
    "        delta_sim = test_stat(*sample_from_groups(group_a, group_b, replace))\n",
    "        delta_dist_h0.append(delta_sim)\n",
    "\n",
    "    # Produce plot of delta distribution with observed delta as dashed line\n",
    "    sns.distplot(delta_dist_h0, fit=sp.stats.norm, kde=False)\n",
    "    plt.axvline(delta_star, color='b', linestyle='dashed', linewidth=2)\n",
    "    \n",
    "    # Comparing against absolute value of real statistic for ease\n",
    "    if tails == 1:\n",
    "        H1 = \"A > B\"\n",
    "        p_val = np.mean(delta_dist_h0 > delta_star)\n",
    "    elif tails == 2:\n",
    "        H1 = \"A = B\"\n",
    "        p_val = np.mean(delta_dist_h0 > abs(delta_star))\n",
    "    \n",
    "    print(\"p-value\", str(round(p_val*100, 3)) + '%')\n",
    "    \n",
    "    # Result of test.  Divide by tails, assuming symmetry\n",
    "    if p_val < alpha/tails:\n",
    "        print(\"Reject H0 | \" + H1)\n",
    "    else:\n",
    "        print(\"Don't Reject H0 | A = B\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate and Plot Groups to Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x112399518>"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD3CAYAAAANMK+RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADwRJREFUeJzt3V+MXGd5x/Hvrtf2YnkcbaUBtVFCSmmf9KKKE5c6iASv1LghRMgtbaqoIqUghFL5IlajEgU5iivBBVJwVSCUyBAMCESFA1y0cvANSVzj4pLGEhHhSUMaUrmqtLHseB3XadbeXszZdmtm/8x4zszsu9+PZOnMOTNzHr/7zm/eeeecMyOzs7NIkla+0UEXIEnqDQNdkgphoEtSIQx0SSqEgS5JhRgb5M6npqaLPcRmYmIDp06dG3QZA2c7tNgOLbZDy+W2Q7PZGGm33hF6TcbG1gy6hKFgO7TYDi22Q0td7bCsEXpE/Atwprr5b8AjwN8AM8ChzPyriBgFPg9cB7wOfCQzX+h9yZKkdpYM9IgYB0Yyc3LeuuPAHwIvAv8QEdcDvwqMZ+Y7I+JG4NPAjlqqliT9gpGlzhSNiK3AV4Gf03oD2AM8kpm/WW2/B1gH/DJwLDO/Wa0/kZlXLvbcMzMXZv0IJkkdazuHvpwpl3PAQ8AXgV8HDgKn522fBt4GbAJenbf+QkSMZebMQk9c8pcjzWaDqanpQZcxcLZDi+3QYju0XG47NJuNtuuXE+jPAy9k5izwfES8CvzSvO0NWgG/oVqeM7pYmEuSems5R7l8mNZ8OBHxK7SC+7WI+LWIGAFuBQ4DR4D3Vve7EfhxLRVLktpazgj9S8D+iPhHYJZWwF8Evg6soXWUyw8j4p+B7RHxA1rzOx+qqWZJUhtLBnpm/jfwJ2023XjJ/S4Cd/eoLklShzyxSJIKMdBT/6VeeeL4ia4eN7l50SNrpRXFEbokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCM0U1dLo961Na7RyhS1IhDHRJKoSBLkmFcA5dq1o38/VeoVHDyhG6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwjNFpQ51c3bpHduvraES6f9zhC5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqxLJOLIqINwNPA9uBGWA/MAs8C+zMzIsR8SBwe7V9V2Yeq6ViSVJbS47QI2It8AjwX9WqvcDuzLwZGAF2RMQNwDZgK3An8HA95UqSFrKcKZeHgC8A/1Hd3gI8WS0fBG4BbgIOZeZsZr4MjEVEs9fFSpIWtuiUS0T8GTCVmd+LiPur1SOZOVstTwNXAJuAk/MeOrd+arHnn5jYwNjYmm7qXhGazcagSxgKnbZDY+N4TZUMlv2hxXZoqaMdlppD/zAwGxG3AJuBrwJvnre9AZwGzlTLl65f1KlT5zoqdiVpNhtMTU0PuoyB66Ydps+er6mawbI/+LqYc7ntsNCbwaJTLpn57szclpmTwHHgT4GDETFZ3eU24DBwBLg1IkYj4mpgNDNf6bpaSVLHurl87r3AvohYBzwHHMjMCxFxGDhK601iZw9rlCQtw7IDvRqlz9nWZvseYM9lVyRJ6oonFklSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQnRzPXRJHXr86Esd/xLT5OYr6ylGxXKELkmFcISuWnUzMpXUHUfoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKseRP0EXEGmAfEMAscDdwHthf3X4W2JmZFyPiQeB2YAbYlZnHaqpbknSJ5YzQ3weQme8CdgOfBPYCuzPzZmAE2BERNwDbgK3AncDDtVQsSWpryUDPzO8CH61uvhU4DWwBnqzWHQRuAW4CDmXmbGa+DIxFRLP3JUuS2llyygUgM2ci4ivAHwB/BGzPzNlq8zRwBbAJODnvYXPrpxZ63omJDYyNremm7hWh2WwMuoTBe+EkjY3jg65iKHTaDqX2n1L/X52qox2WFegAmfnBiLgP+CHwpnmbGrRG7Weq5UvXL+jUqXPLr3SFaTYbTE1ND7qMoTB99vygSxi4xsbxjtuhxP7j66LlctthoTeDJadcIuKuiLi/unkOuAj8KCImq3W3AYeBI8CtETEaEVcDo5n5StcVS5I6spwR+reBL0fEU8BaYBfwHLAvItZVywcy80JEHAaO0nqj2FlTzZKkNpYM9Mx8DfjjNpu2tbnvHmDPZVclSeqYJxZJUiEMdEkqhIEuSYUw0CWpEAa6JBVi2ScWSU8cP9HxYzxLVOofR+iSVAhH6NKQ6uYTEcDk5it7XIlWCkfoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhTDQJakQBrokFcJAl6RCGOiSVAgDXZIKYaBLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklSIsUEXoP574viJQZcgqQaO0CWpEIuO0CNiLfAocA2wHvgE8BNgPzALPAvszMyLEfEgcDswA+zKzGP1lS1JutRSI/QPACcz82bgPcDngL3A7mrdCLAjIm4AtgFbgTuBh+srWZLUzlJz6N8CDlTLI7RG31uAJ6t1B4HfAxI4lJmzwMsRMRYRzcycWuzJJyY2MDa2puvih12z2Rh0CW01No4Xvb9h1a92GNZ+N2fY6+uXOtph0UDPzLMAEdGgFey7gYeq4AaYBq4ANgEn5z10bv2igX7q1Lnuql4Bms0GU1PTgy6jremz5/u2r8bG8b7ub1j1sx2Gtd/BcL8u+uly22GhN4MlvxSNiKuA7wNfy8xvABfnbW4Ap4Ez1fKl6yVJfbJooEfEW4BDwH2Z+Wi1+pmImKyWbwMOA0eAWyNiNCKuBkYz85WaapYktbHUHPrHgQnggYh4oFp3D/CZiFgHPAccyMwLEXEYOErrTWJnXQVLktpbag79HloBfqltbe67B9jTk6okSR3zxCJJKoSBLkmFMNAlqRBenEsqTDcXX5vcfGUNlajfHKFLUiEMdEkqhIEuSYUw0CWpEAa6JBXCQJekQhjoklQIA12SCmGgS1IhDHRJKoSBLkmFMNAlqRAGuiQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhfA3RbWqvXj+xx0/5m3jv1VDJdLlM9Al+cPShXDKRZIK4Qh9hetmZCWpTI7QJakQBrokFcIpFxWhm6NVpNI4QpekQhjoklQIA12SCuEcuoaO8+FSd5YV6BGxFfhUZk5GxNuB/cAs8CywMzMvRsSDwO3ADLArM4/VVLMkqY0lAz0iPgbcBbxWrdoL7M7MJyLiC8COiPg5sA3YClwFPAa8o56SpcHq5hPEetby+vk3OnqM14xRp5Yzh/4z4P3zbm8BnqyWDwK3ADcBhzJzNjNfBsYiotnTSiVJi1pyhJ6Zj0XENfNWjWTmbLU8DVwBbAJOzrvP3PqpxZ57YmIDY2NrOip4JWk2G7Xvo7FxvPZ9XK5Oa1zP2poqGaz14539v4b9b9tt/+7H62IlqKMduvlS9OK85QZwGjhTLV+6flGnTp3rYvcrQ7PZYGpquvb9TJ89X/s+LscJftrxVEOJ1o93PuUyzXD/bbvp3/16XQy7y22Hhd4Mujls8ZmImKyWbwMOA0eAWyNiNCKuBkYz85VuCpUkdaebEfq9wL6IWAc8BxzIzAsRcRg4SutNYmcPa5QkLcOyAj0zXwJurJafp3VEy6X32QPs6V1pkqROeKaoJBXCQJekQhjoklQIA12SCuHFuYbIsP8+aFenvHd4Mo2k7jlCl6RCGOiSVAinXKQh1e114b1K4+rlCF2SCmGgS1IhDHRJKoSBLkmFMNAlqRAe5SKpK92cCHfH9mtrqERzHKFLUiEMdEkqhIEuSYVwDn0V6vYMRK0M3fx9Pbu0DI7QJakQBrokFcJAl6RCGOiSVAgDXZIK4VEuNXn86EtMnz0/6DIkrSIGuqS+6WagM7n5ynqKKZBTLpJUCANdkgrhlIskzy4thCN0SSqEI/Rl6Oa6z42N4zVU8ou8LoukOY7QJakQjtAlDbVuPiF3a6UfImmgDxGnT7SSdNNfr9v4jhoq0RynXCSpEAa6JBWip1MuETEKfB64Dngd+EhmvtDLfUiS2uv1HPrvA+OZ+c6IuBH4NLCjx/voWj+/XPnp2eO8fv6Nvu1PWgm6eV14AtPy9TrQbwIeB8jMf4qI3+7x8/+vfoZzN1/+rB9fW0MlkurUr1y5Y/u1tTxvrwN9E/DqvNsXImIsM2fa3bnZbIx0u6O6GqS9fu5L0mrQbDZ6/py9/lL0DDC/ytGFwlyS1Fu9DvQjwHsBqjl0D6yWpD7p9ZTLd4DtEfEDYAT4UI+fX5K0gJHZ2dlB1yBJ6gFPLJKkQhjoklQIA12SCuHVFnskIrYCn8rMyYh4O7AfmAWeBXZm5sVB1tcvl7TD9cDfA/9abf7bzPy7wVXXHxGxFngUuAZYD3wC+AmrrE8s0A7/zirrExGxBtgHBK2//93AeWroD47QeyAiPgZ8EZj7maK9wO7MvJnW0T5Dc/mDOrVphy3A3sycrP4V/cKd5wPAyerv/x7gc6zOPtGuHVZjn3gfQGa+C9gNfJKa+oOB3hs/A94/7/YW4Mlq+SBwS98rGox27XB7RDwVEV+KiN6fGjecvgU8UC2PADOszj6xUDusqj6Rmd8FPlrdfCtwmpr6g4HeA5n5GDD/ikMjmTl3POg0cEX/q+q/Nu1wDPjLzHw38CLw4EAK67PMPJuZ01VYHaA1Klt1fWKBdlitfWImIr4CfBb4OjX1BwO9HvPnwhq03pFXo+9k5tNzy8D1gyymnyLiKuD7wNcy8xus0j7Rph1WbZ/IzA8Cv0FrPv1N8zb1rD8Y6PV4JiImq+XbgMMDrGWQvhcRv1Mt/y7w9GJ3LkVEvAU4BNyXmY9Wq1ddn1igHVZdn4iIuyLi/urmOVpv7j+qoz94lEs97gX2RcQ64DlaHzdXoz8HPhsRbwD/yf/NI5bu48AE8EBEzM0h3wN8ZpX1iXbt8BfAX6+yPvFt4MsR8RSwFthFqw/0PCM89V+SCuGUiyQVwkCXpEIY6JJUCANdkgphoEtSIQx0SSqEgS5JhfgfVGb3XccenQ0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x1 = np.random.binomial(49, .357, 4000)\n",
    "x2 = np.random.binomial(33, .53, 1000)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "sns.distplot(x1, bins=np.arange(x1.min(), x1.max() + 1), ax=ax, kde=False)\n",
    "sns.distplot(x2, bins=np.arange(x2.min(), x2.max() + 1), ax=ax, kde=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Hypothesis Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value 31.6%\n",
      "Don't Reject H0 | A = B\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD3CAYAAAAALt/WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VNed//H3jLpQQeChY0SRDhICRDNCNFGNbWxciHuL45L17i9eezdx7E2yySZrx+tN8SZxNq6xib2OjQ02xZgqECoU0QQShyo6QoAaCPX5/SEJi6I+M2fK9/U8eizmztz70bHmO1fnnnuOxW63I4QQwnNZTQcQQgjROVLIhRDCw0khF0IIDyeFXAghPJwUciGE8HD+rj5gYWGZWw+TiYoKpaio3HQMozy9DZ749VoA3vvx9A7vw9PbwFGkHdynDWy2cEtz2+SM/Cr+/n6mIxgnbSBt0EjawTPaQAq5EEJ4OCnkQgjh4aSQCyGEh5NCLoQQHs7lo1aEcLbOjFYRwhPJGbkQQng4KeRCCOHhpJALr/OL97fwi/e3mI4hhMtIH7nwOkcKykxHEMKlpJAL0YrUHSeu+3hKYl8XJxHi+qRrRQghPJwUciGE8HBSyIUQwsO12keulPID3gYUYAe+r7Xe3WT788CTQGHDQ89orbUTsgohhLiOtlzsvB1Aaz1RKZUC/Ccwr8n2McCjWutsx8cTov2mjOxjOoIQLmWx21tf50Ep5a+1rlFKPQZM11o/1mRbHrAH6AUs01q/2tK+ampq7Z4wv6/wbisy8695bM6EaOx2O1VVVQQGBmKxWJp9buPzhXChZheWaNPww4Yi/gFwFzD/qs2fAH8CSoFFSqm5Wuulze3LHVbaaInNFk5hoW+PQ/aFNii7UHH5++LzZ9iU9jV//M/dHD58iLq6OiIiIhgyJJZZs+YQcMMwAoOCr9mHt7cR+MbvQmvcpQ1stvBmt7V5HLnW+jGl1IvAJqVUvNb6olLKAvxea10CoJRaBowCmi3kQjhb/ulSAKJ7RbT4vKqqSlYv+Rub076mrq6WwMBA4uLiCQnpwrlzZ9i+PZtt27YSEdmdOXd/l2GJyZfP0oVwJ2252PkI0K+hy6QcqGv4AogAdiul4oCLwHTgPSdlFaJN/uOvW4GWZ0E8e+YEH7/1KmfPnKC7rTdTb/4O//DYfIIazrxttnBycvaxePHnfPrp//Hp+//N2OTZ3Dr/Sfz9A1zycwjRVm05I/8CeF8ptQEIAP4ZuEspFaa1fksp9TKwDqgE1mitlzsvrhCdl5u7h3d+9xLlF8tImjqXWbc/TEBg0OUi3qhXr958//v/hG1wEn9/73W2Zqzk3NlTPPT0vxEYGGQovRDXarWQa60vAve2sH0BsMCRoYRwlKtvry84eYR33/g3KisuMe/+ZxmTPKvVfXS39ebJ519l4Qe/ZW/OZj5++1UeevplZ0UWot3khiDhM0qKzvLhn/+DiksXufvhH7SpiDcKDAzi3u/+K0MTxnFI7+SrT/5MW0Z8CeEKMmmW8Ak1NdV8+v5/U1ZynpvvfJyR46a2ex/+/gF857v/ynv/8xN2bknlld/9mYnT513zPJlMS7ianJELn7DqqwUcy9cMHzOZ5Gl3dHg/AQGBPPDkjwmPiGLVVws4cfSAA1MK0TFSyIXXO7Qvh8zUJdzQsy933P8PnR5CGBHZjbsf+Wfq6mpZ+MFvqaqsaP1FQjiRFHLhdX72+Fh+9vhYACorL7H44z9isVi5++EfEBQU4pBjDFYjSJ4+j3OFp1i7/P8csk8hOkoKufA60b0iLt8MlPr13yk+f4ZJM+6k34BYhx5nxq0P0O2GXmSmLuXUsUMO3bcQ7SGFXHit/PzDZKYuJap7T1LmNDuC9rLUHScuf63IzL/8fXMCAoOYe+8z2O11fPXp/1JXV9fsc4VwJhm1IrzOX7/eC9jZtepN6upqueXuJwhw0g08Q4YmkjB6Eru3bSQnO61Do2GE6Cw5IxdeZ8POk2zYeYqtWzcxeGgiKmGcU4836/ZH8PcPYPXSv1FdVenUYwlxPVLIhdeyWCzcPO9Rp090FdW9B0kpcykpOkvWBpmhQrieFHLhtWbOvJlefQe65FiTZ95DcEgX0tcsorz8okuOKUQjKeTCq9TW1l7+/oknnnbZcUNCu5A8fR7lF8v47LNPXHZcIUAKufAy69evvfx93779XHrsCVPnEtolnL///WPKy917ARXhXaSQC69ht9v58MP3jR0/KDiEpKlzuXChjGXLvjKWQ/geKeTCa2zZsolDhw4QWFfKgJ7NL4vlTOMmzSEgIJAFHy1gTfbRVseiC+EIUsiF1/j0048B+Oe7Y/n37zp3yGFzuoRFMCppBsXnC8ndkWEkg/A9UsiFV8jPP8TmzVmMGJHI0KHxRrMkp9yOxWJl45rFMme5cAkp5MIrNI4UuffeBwwngW623sSPHM+p44c4vH+36TjCB0ghFx6vuLiIb75ZTp8+fZk4cQpP/HotT/x6besvdKLk6XcCkL52sdEcwjdIIRce78svv6Cqqor58+/Hz8/PdBwA+kfHcuOgOPbnbuNc4SnTcYSXa3XSLKWUH/A2oAA78H2t9e4m228HfgbUAO9prd92UlYhrlFdXc2iRQsJCg4lvN+YK0aImB4tctPkWzh6KI+tGSu5Z9ZYo1mEd2vL7Ie3A2itJyqlUoD/BOYBKKUCgN8B44CLQLpS6iutdYFz4gpxpfT0NM6fP0fS1LkEBTtm0QhHiR+RRGiXCLZnrWX1lgfxDwi45jmyvqdwhFa7VrTWi4HGe50HAMVNNscBB7TWRVrrKmAjMMXhKYVoxtKl9X3QY5NnGU5yLf+AAEaNn0b5xVJyd2WZjiO8WJvmI9da1yilPgDuAuY32RQBlDT5dxkQ2dK+oqJC8fd3j37M5thsZm4mcSee0AYnTpxgy5ZNJCYmMjjm2tV/wsOCO7X/zr4eYOqs20lf+yXbs1Yxceq1Hzae0M6ekNHZ3L0N2rywhNb6MaXUi8AmpVS81voiUAo0/QnDufKM/RpFRe49B4XNFk5hYZnpGEa5exs09n2vWfYxdrud6PhJlF24dgHk6z3WVuFhwZ16faPgLt0ZGDucgzqHQwcPYOt55fwv7tzO4P6/C67gLm3Q0odJq10rSqlHlFIvNfyzHKhr+ALIA2KUUt2UUoHUd6tkdi6uEK2rra1le9ZagoJDSRg18YptScN6kjSsp6Fk1xqXfDMAW9NXGk4ivFVbhh9+AYxSSm0AvgH+GbhLKfW01roaeKHh8UzqR63IxBLC6Q7kbaO05BwjxkwmMOjKLpDY/l2J7d/VULJrDR1xE13CI9mxeR3V1VWm4wgv1GrXSkMXSrMr12qtlwBLHBlKiNZkZ6wCYIwbXuS8mr9/AIk3TSN9zWL27d7KsFHJpiMJLyM3BAmPU1ZaxL7cbHr3G0Sf/oOv2b7vWDH7jrV4qcblEselALBjS6rRHMI7SSEXHicneyN1dXWMTppx3e1ZewrI2uNetzL07DOA3v0GsT93GxfK3OtDRng+KeTC4+zckorV6nfNRU53l3hTCnV1tezettF0FOFlpJALj5Kff4hTxw8xJG4UXcJbvGXB7QwfMxmr1cqOzammowgvI4VceJSVK1cAMHLcVMNJ2i8svCtD4kZz8thBzpw6ajqO8CJSyIXHqKurY9WqFQQFhTA0wcwKQJ2VeFMKIBc9hWNJIRceY9euHRQUnCY+cQIBgUGm43SIShhHcEgou7ZsoK6u1nQc4SWkkAuPsXLl1wCMbBjK54kCAgIZNmoipSXnyD+QazqO8BJSyIVHqKysJDV1DTZbD6KHDGvxuY/OUTw6R7koWfsNHz0ZQEavCIeRQi48wubNWVy4cIGZM2djtXr2r230kHjCIqLI3ZlJTU2N6TjCC3j2O0L4jHXrVgMwbdpMw0k6z2r1Y1hiMuUXy8jO3mI6jvACUsiF26usrCA9PY3evfugVFyrz1+akc/SjHznB+uE4aMnAbBmjcyIKDpPCrlwe5s3Z3HpUjnTps3EYrG0+vzzpZWcL610QbKO6xcdS2SUjbS0VCor3TurcH9SyIXb86ZulUZWq5WE0RO5ePEimzfLMnCic6SQC7fW2K3Sp09fYmPddyRKRwwfVd+9snbtKsNJhKeTQi7c2qZNmVy6dKnN3SqepHf/QfTt25/09A1cunTJdBzhwaSQC7e2bt0awLu6VRpZLBamT59JRUUFmzZlmI4jPJgUcuG2KioqyMhIo2/ffsTExJqO4xRTp04HYP36dYaTCE/W6lJvQpiSlZXRoW6VmH6eM71tTEwsvXv3ISNjI5WVlQQFeeYcMsIsOSMXbis1tWPdKhMSejEhoZczIjmcxWJh6tTpXLpUztatm03HER5KCrlwS992q/RnyJAY03GcaurUaQCsX7/WcBLhqVrsWlFKBQDvAdFAEPArrfVXTbY/DzwJFDY89IzWWjsnqvAlWVnpVFRUMH36jHaPVjlXUgFA98hgZ0RzuLi4YdhsPUhPT6OmpgZ/f+nxFO3T2hn5w8A5rfVkYA7wx6u2jwEe1VqnNHxJERcO0ZnRKssyj7As84ijIzmN1WplypQUyspK2b4923Qc4YFaK+SfAT9t+N4CXD1V2xjgJaXURqXUS44OJ3xTZWUlWVkZ9OnTl8GDvbtbpdG3o1eke0W0X4t/w2mtLwAopcKBhcBPrnrKJ8CfgFJgkVJqrtZ6aUv7jIoKxd/fr+OJXcBmCzcdwTiTbbB+/TYuXSrn/vvvo0ePiGu2h4e1rcukrc9z1uvborGdp0+fRLdu3UhP38Arr/wSPz/3eY/I+8H926DVzjilVH9gEfCm1vrjJo9bgN9rrUsa/r0MGAW0WMiLiso7FdjZbLZwCgvLTMcwynQbLF1avxLQ2LHJ181RdqGiTftp6/OuJzwsuFOvb6umP9/EiVNYsmQxa9duJDFxtNOP3Ramfxfcgbu0QUsfJi12rSilegIrgRe11u9dtTkC2K2UCmso6tMB6eATnVJbW0t6ehpRUd2Ij08wHcelvh29IjcHifZprY/8ZSAK+KlSKrXh6yGl1NMNZ+IvA+uANGCP1nq5k/MKL5ebu5uiovNMmjTFrboXXGHUqLGEhYWzYcM66urqTMcRHqS1PvLngOda2L4AWODoUMJ3bdiQCsCkSVPNBnGR1B0nrvj3kPgx7Nicyt69uT73F4noOLkhSLgNu93Oxo3rCQkJZfTosR3ez20TBnDbhAEOTOY68SMnADJ6RbSPFHLhNg4dOsiJE8dJSkru1Jwj3SODPeZmoKsNHppIYFAwGzakYrfbTccRHkIKuXAbGzeuB2DyZN/oVrmegIBAYuJGc+LEcQ4dOmg6jvAQci+wMK6xn3j5N6vw8/OnpsvAa/qO2yNz92kAj5k462pxI5PYsyODtLRUBg8eYjqO8AByRi7cQvH5M5w6foiBMcMJDunSqX3tP17C/uMlDkrmerHxo/H39yctLdV0FOEhpJALt7A3p34K17iR4w0nMS84pAtjxoxj//59nDzZ8b9MhO+QQi7cQu7OTQAMTRhnOIl7mDw5Bfj2uoEQLZFCLowrv1jKkYO59IuOJTyym+k4bmHSpClYLJbL4+qFaIkUcmGc3r0Vu72OuBFJpqO4jW7dupOQMIKcnJ0UFZ03HUe4OSnkwri8XfXdKnEjbjKcxL1MnpyC3W4nPX2D6SjCzUkhF0ZdunSJA3t3YOvVnxt69HXIPrtFBNEtwvMXMZ4yJQWADRukn1y0TMaRC6O2bMmiprrKoWfjc5OjHbYvk+oX1hhCdvZmLl68QJcuYaYjCTclZ+TCqLS0+rNN6R+/vsmTU6iuriYrK9N0FOHGpJALY2pqasjI2EhEZHf69B9sOo5bauxekZuDREukkAtjdu7cTllZKUNH3ITFYnHYfj9coflwhXesAz54cAy9e/chKyuDqqoq03GEm5JCLoxpPMuUbpXmWSwWJk9Oobz8ItnZW0zHEW5KCrkwwm63k5a2nvDwCKKHxJuO49ake0W0Rgq5MELrPAoLz5CcPBE/Pxk81ZJhw4YTFdWNjRs3UFtbazqOcENSyIURjaNVJk1KMRvEA/j5+TFp0hSKi4vYvXuX6TjCDUkhF0Zs3LiewMAgbrpJ+sfbonESLeleEdcjhVy43LFjRzl8+BDjxo0nJCTEdByPMHr0WEJDu5CWtl6WgBPXaLFzUikVALwHRANBwK+01l812X478DOgBnhPa/2286IKb+HsJd2ShvV0yn5NCgwMZMKEiaxZs5KDB/czZEis6UjCjbR2Rv4wcE5rPRmYA/yxcUNDkf8dMBuYCjytlPK+d5BwuA0bUrFarSQnT3bK/mP7dyW2f1en7Nukxg8+mdpWXK214QKfAQsbvrdQf+bdKA44oLUuAlBKbQSmNLymWVFRofj7+3UsrYvYbOGmIxjnrDY4c+YMe/bkMG7cOGJi+gMQHuaeK96bzHW99p8792ZeeeUXZGam8eKL/2I0i69x9zZosZBrrS8AKKXCqS/oP2myOQJoujBiGRDZ2gGLisrbn9KFbLZwCgvLTMcwyplt8NVXywFISpp8+RhlFyoceox9x4oBOnVWHh4W7PBc7dFc+48ZM47MzHR27txLnz6OmS2yJfJ+cJ82aOnDpNWLnUqp/sA6YIHW+uMmm0qBpnsOB4o7mFH4iMZuAWf1jwNk7Skga0+B0/ZvUo/oRADe/+RLUnecIHWHrOkpWr/Y2RNYCfyT1nrNVZvzgBilVDfgAvXdKv/tlJTCK3ydsY/s7K307jeIvafr2HtailBzmivQQxPG8ZXFSt6uLCZOn+fiVMJdtdZH/jIQBfxUKfXThsfeBrpord9SSr0AfEP9mf17Wmt5Z4pm7duzlbq6WuJGytjxjuoSHsmNg4Zy9FAeF0qLCYvwvou6ov1a6yN/Dniuhe1LgCWODiW8U+OSbvEjxhtO4tniRoznyMFc9uZsZuzE2abjCDcgNwQJl6isrGB/3ja623pj69XfdByP1jhbZO6uLMNJhLuQQi5cYuvWLVRXVRI3Ismhc4/7oqjuPejVdyCH9+VQcemi6TjCDUghFy7x7dzj0q3iCHEjxlNbW8O+3G2mowg3IIVcOF1NTQ3p6RsIj4ii74AYpx/v0TmKR+copx/HpPiGC8Z5O6V7RUghFy6Qk7OTkpISho4Yj9Uqv3KO0KP3jXS7oRf787ZRWVlpOo4wTN5VwukabwKSbhXHsVgsxI0YT1VlBdnZm03HEYbJ0izCqex2Oxs3ricsLJyBMQkuOebSjHwA5iZHu+R4psSNSCJ97Zd8ung5VaGDrtmekuj8W/iFe5AzcuFU+/btpaDgNMnJk1y2pNv50krOl3p/d0O/6FjCIqLYm7NFloDzcVLIhVM0zgOy4LP6+8W69R9uOJH3sVqtDB0+jvKLpRw9lGc6jjBICrlwqrxdm/APCGTI0FGmo3ilxpuDGu+aFb5JCrlwmrNnTlB4+hhD1EgCg9xzznFPNzAmgeCQUPJ2bZIl4HyYFHLhNHk7688SZZIs5/H3DyA2fiwlRYWcOn7IdBxhiBRy4TR5uzZhtVpRCWNNR/FqjcM6pXvFd0khF05RUnSW40f2MWBwPKFdIlx67Jh+kcT0a3WxKq8xJH4U/v4B5MkkWj5LCrlwij07MgEYljjR5ceekNCLCQm9XH5cU4KCQhg8NJEzp45x7sxJ03GEAVLIhVPs2ZGBxWK9PCeIcK7G7pVc6V7xSVLIhcMVFBRw7PBeoofEG1nB5lxJBedKzC2cbIJKGIelYQk44XukkAuHW79+LWCmWwVgWeYRlmUeMXJsU7qERTBgcDzH8/dRWnLedBzhYlLIhcOlpq6RbhUD4kfWd6/s3SWTaPkaKeTCoQoKCti9e5exbhVfNnR44zBE6V7xNW2axUgpNR54TWudctXjzwNPAoUNDz2jtdYOTSg8iuluFV/WtZuNPv0Hc3j/bsovlpmOI1yo1UKulPoR8AhwvcUBxwCPaq2zHR1MeKbU1DVYrdKtYsqwxGROHjtI3q5N3DpxqOk4wkXa0rVyELi7mW1jgJeUUhuVUi85LpbwRI3dKiNHjpJuFUMSRk8CICc7zXAS4UqtnpFrrT9XSkU3s/kT4E9AKbBIKTVXa720pf1FRYXi7+/X7qCuZLOFm45gXEfaYPnyLwC4/fbbCA8zP0lWZzO4w8/QXuFhNxI9OI7D+3cDFdhstk7vU94P7t8GHZ7pXyllAX6vtS5p+PcyYBTQYiEvKirv6CFdwmYLp7DQt/sXO9oGS5cux2q1Mnr0BHYdNTeO+7YJAwAou9DxDOFhwZ16vUnxicnkH8zj88+/4p577u3UvuT94D5t0NKHSWdGrUQAu5VSYQ1FfTogfeU+qmm3Srdu3Y1m6R4ZTPdIzzubdpRhiclYLFbWrl1lOopwkXYXcqXUg0qppxvOxF8G1gFpwB6t9XJHBxSeoXG0yrRpMwwnEeGR3YiOGUZOzk4KCk6bjiNcoE1dK1rrfCCp4fuPmzy+AFjglGTCo6xZsxKr1cqUKdNMRyFzd33x8qWJs642fPQkDu/LYe3aVTzwwCOm4wgnkxuCRKek7jjB56u2kpe3h4GxI9h1tILUHSeMZtp/vIT9x0uMZjAtfmQSfn5+rFkj3Su+QAq56LRdWzcAMHLsFMNJRKPQLhGMG5fEvn17OXbMt+ad8UVSyEWn2O12crLT8A8IvHyLuHAPM2fOBpCzch8ghVx0yqnjhzh75gRDE8YRHBJqOo5oYuLEyQQGBrF69UpZmNnLSSEXndLYrTJ8zGTDScTVunQJY8KEiRw9ms+hQwdMxxFOJIVcdFhtbS052zYSHNKFmLjRpuOIq6TuOEGvQWMAePejhaTuOGH8QrRwDinkosN27txOWcl5ho1Kxj8gwHScy7pFBNEtIsh0DLcQO2wMQcGh7Nq6gbq6OtNxhJNIIRcdtnLl1wCMGONeo1XmJkczNznadAy3EBAYxLDEZEqLzzXMvyK8kRRy0SEVFRWkpq4lMuoGBgyONx1HtCBxfP1NWju3rDOcRDiLFHLRIRs2pFJefpHEm6ZhtcqvkTu7ceBQorr3JHdHFpWVl0zHEU4g70DRIStW1E9ymTguxWyQ6/hwhebDFbJQVSOr1crIcSlUVVWQt1OWgfNGUshFuxUUFJCdvYXhw0fSvUcf03FEGyTelALA9s3SveKNpJCLdlu5cjl2u505c24zHUW0UbcbenHjoDjy9++WGRG9kBRy0S52u52vv15KYGAQ06bNNB1HtMOom6Zht9svjzYS3kMKuWiXPXtyOH78GFOmTCUsLMx0HNEO9eP9A1m+fKncsu9lpJCLNmm8K/C9v30KQJ+YCXKXoIcJDunCsMRkTpw4xvbtspiXN5FCLtqssvISu7I3EBl1A4PUcNNxRAeMSZ4FwNKliw0nEY7U4cWXhe/JyU6jqrKCidPvxGr1Mx2nWUnDepqO4LYGDIrjxhujWb9+HSUlxURGdjUdSTiAnJGLNtuasQqLxcroJPdelzO2f1di+0uBuh6LxcLcufOorq5m5coVpuMIB5FCLtrk5LGDnDx6gNhhY4iMusF0HNEJc+bcir+/P0uXLpaLnl5CCrlok60ZKwEYO3G24SSt23esmH3Hik3HcFtdu0YxefJUDh8+xJ49OabjCAdoUyFXSo1XSqVe5/HblVJblFKZSqmnHJ5OuIXy8nJ2ba2/yBkTN8p0nFZl7Skga0+B6Rhube7cOwFYskQuenqDVgu5UupHwDtA8FWPBwC/A2YDU4GnlVJylckLrVmzkqrKCkYnzXTri5yi7caMGUefPn1Zs2YVJSXy14una8uolYPA3cCCqx6PAw5orYsAlFIbgSnAZy3tLCoqFH9/9y4GNlu46QjGNbaB3W5n2bLFWCxWpsy4lfCw4FZe6T46m9WTftb2aPx/+8gjD/Paa6+RmvoNTz75ZKvP92Xu3gatFnKt9edKqejrbIoASpr8uwyIbG1/RUXlbQ5ngs0WTmFhmekYRjVtg127dpCbm0v8yCT8AsMpu1BhOF3bdSZreFiwR/2s7dH4/3bKlNm88cYbfPTRx8ydOx8/v2tPsOT94D5t0NKHSWcudpYCTfccDsjfaF5m4cJPAEiaOtdwEuFo4eHhzJ59C6dPnyIjY6PpOKITOlPI84AYpVQ3pVQg9d0qmY6JJdxBQcFp0tLWM2RIjKwC5KXuuus7AHzxxaeGk4jOaHchV0o9qJR6WmtdDbwAfEN9AX9Pay2Tb3iRRYsWUltby/z592OxWEzHEU4wePAQRo0aQ3b2FvLzD5uOIzqoTbfoa63zgaSG7z9u8vgSYIlTkgmjKioqWLJkMZGRXZkxYzaZeWdNR2qzR+co0xE8yt1338v27dl88cVnvPDCj0zHER0gNwSJ61q1agVlZaXcccddBAUFmY4jnGjixMn07NmLFSuWylBEDyWFXFyjrq6Ozz77P/z8/LjzzntMxxFO5u/vz733PkhFRQWLFi00HUd0gBRycY3U1FTy8w8zc+ZsbLYepuO029KMfJZm5JuO4VFuu+0OwsMj+PzzT6mo8M5hl95MCrm4wrrtx/mv3/4BgJjRcy4vKOFJzpdWcr600nQMjxIaGspdd82npKSYr79eajqOaCcp5OIK+Qf2cOTgXoYmjKNH7xtNxxEudM899xIYGMgnn/yNmpoa03FEO0ghF1dIW/0FAJNnSd+4r4mK6sYtt8zl1KmTrF+/znQc0Q5SyMVl+/ZpDuRtZ8jQEfQfKEP4fEljF1r0iJlYLFb+8s47rNt+3HQs0Uay1Ju47KOPPgBgxq33Gk4inKmlax7dbb1JGJVMzraN7M3ZDLPjXJhMdJSckQsADh06SGrqGnr3G4QaNtp0HGFQyi33YbFYWff1J9TV1ZmOI9pACrkA4N13/4LdbmfGbQ96/O34Mf0iienX6kScohm2nv0YMWYyp0/ks3r1atNxRBtI14pg795c0tJSGTZsODHxnn82PiGhl+kIHi9lzr3kbEvj1f/6Ld//4W+xWq8850tJ7GsombgeOSP3YY0XuP7rd28AMG7adzz+bFw4RvcefRg5LoXTJ46wZ0eG6Thqi7cgAAAPSUlEQVSiFVLIfdyRg7kcyNvOwNjhDIodbjqOQ5wrqeBcidyd2FlTZ8/HarWybvkn1NbWmo4jWiCF3IfZ7XZWL/0IgBm3PWg4jeMsyzzCsswjpmN4vG623oyffDNnz5wgO2OV6TiiBVLIfVjeriyOHMxFJYzjxoFDTccRbuiWOx8hMCiYtcv/j4pLF03HEc2QQu6jKisr+WbxB1itftw87zHTcYSbCo+MYsqs+ZRfLGXDys9NxxHNkELuoz777BOKzhUwfsqt3NBTRiCI5k1ImUtklI3M1CUUnSswHUdchxRyH3T27FkWLHif0C4RpMy5z3Qc4eYCAoOYdccj1NbWsPLLD03HEdchhdwHvfPOn7l0qZwZtz1ASGgX03GEBxg+ehL9BsSyZ0cGB/buMB1HXEUKuY/ZuXM7y5cvYdCgIYyeMMt0HOEhLBYLc+99GovFytJP/0JlpQzvdCet3tmplLICbwIjgUrgSa31gSbb3wAmAWUND83TWpc4IavooMZJkqqrq3jztV9gsViYPu97+Pn5GU7mHLdNGGA6glfq038wE1LmkrHuKz788D2eeupZ05FEg7ackd8JBGutJwA/Bn5z1fYxwM1a65SGLynibmr9N59x7sxJxk+51auHG3aPDKZ7ZLDpGF5p2q33Exll4+OPF3D48EHTcUSDthTyScAKAK11FjC2cUPD2XoM8JZSKl0p9YRTUopOO30in42rFxEZZWPG3IdMxxEeKigohLnfeYra2lpef/1VuePTTbRl0qwIoOlZdq1Syl9rXQN0Af4A/BbwA9YppbZqrXc1t7OoqFD8/d37T3qbLdx0BIcKDfZnyd/fpK6ulvse/wE3dI9q9TXhYZ57Rrsu+xgA08b079R+PLkNHOnqdhibNJlzR+ewYsUKli5dyJNPPmkomeu4e01oSyEvBZr+FNaGIg5QDryhtS4HUEqtpb4vvdlCXlRU3sGormGzhVNYWNb6Ez3IVws/4Fj+fkaOS6HfwATKLrR8oSo8LLjV57iz3MPnARirbB3eh6e3gaM01w7PPvsCW7Zs5Y033mDo0JEo5b1dde5SE1r6MGlL10o6cCuAUioJyGmyLRZIV0r5KaUCqO+G2dbxqMLRtm3bStqqz4nq3pPb5nv/mZNwja5du/LSSz+jpqaGX/7yp1RUyIeeSW0p5IuACqVUBvA74Hml1AtKqTu01nnAAiALWA98qLXe47y4oj2Ki4v51a/+HYvFyvzHXiA4RMaMC8e56aYk5s+/n6NHj/Dmm2+YjuPTWu1a0VrXAd+/6uG9Tba/Drzu4Fyik+x2O6+99ivOni1k5u0P0z861nQk4YWeeeYfyc7ewuLFnzNy5ChmzJhtOpJPkhuCvNSHH75HevoGRo8ey6QZd5mOI7xM46IkmXlnmfvAcwQFhfDKK//BgQP7TEfzSVLIvUjjm+vNDxby7rt/ITLKxsx7/vGaZbqEcCRbr/7c/chzVFdX8fLLP6S4uNh0JJ8j73Avc+r4YT5f8AaBgcE89PRLhIV3NR3J5bpFBNEtIsh0DJ8SN2I8KXPu4/TpU/z85y9TU1PT+ouEw0gh9yKlxef4+O1Xqa6q5O5HnqNX34GmIxkxNzmaucnRpmP4nJQ59zJp0hS2bdvKr3/9S+rq6kxH8hlSyL1EcXERH/zp55QUFTJj7kPEj0wyHUn4GKvVyk9/+kvi4xNYufJr/vjH32O3203H8glSyL3AhQsX+Jd/+QGFBcdJnnYHU2bdYzqS8FEhISG89tpviY4eyMKFn7BgwfumI/kEKeQerrz8Ii+++Dz792vGJM/i5jsfx2KxmI5l1IcrNB+u0KZj+KTUHSfYfvgi85/4N7p2s/HOO//L3/72V9OxvF5bbtEXbqq4uIgf/vA5tN7LzJmzmTT3aZ8v4sI9RHTtzmP/+Av++sef8dZbb5J36DQzb3/4it/PlERZYtBR5IzcQy1as53vfu9xtN7L6KQZTLztaaxW956MTPiW7rbefO+5V+hu603a6i9Y9tlb1NXJbInOIIXcA+Xl5fLO71/iXOEpJs+8m3kP/KPXLhIhPFvXbja+99wr9OwTzeaNK/j47VepuHTRdCyvI4XcwyxZsph/+qenuFBaxM13Ps6sOx6R7hTh1sIiuvLED37J4KGJ7NuTzVu/eZGzBSdMx/IqUsg9REVFBa+99itef/0VgoNDePiZnzBx+jzTsYRok5DQsMu/s2fPnOAvv/kRa9asNB3La8jFTjeXuuMERw7lsfijP3Cu8BS9+w3i/u/9iKjuPU1HE6Jd/Pz8uPnOx+nVdyBf/f3P/OIXP2Hjxg08//wPiYiINB3Po0khd2OXLl1ixaL3yUxdAkDytDuYcduDBATK7ectSRomH3LubOS4qfQbEMPnf3uDNWtWsnnLFm6950niEydgsVhkNEsHWFx951VhYZlb3+rlDquB1NXV8c03y3n77T9z9mwh3W29ufOh/8eAQXEuOb6sjiNt0MiZ7VBbW8vGNYtI/frv1NbWMDAmgVvveZL7bk12yvE6yh1qQkOOZi+GSSG/isn/aXa7nczMjbz77lvs368JDAxiwrQ7mDzrHgJdeBYuRUzaoJEr2uHsmROsWPQ++/ZkY7FYmTPnVh5++HH697/RqcdtKynk1yGF/Fo1NTWsXbuajz76gMOHDwIwe/YtPPXUs+Sdcv0scp5exPYdq59GNbZ/x2d+9PQ2cBRXtsO+3Gy+WfwBhaePYbVamT59Fvfd95Dx9UA9oZBLH7lBx48fY+nSL1mxYinnz5/Hz8+PWbPm8OCDjzB4cAwAeadkmFZ7Ze0pADpXyIXrxcaPYcjQUVhL9/HBB++yevU3rF79DXFxw5g3725SUmYQGhpqOqZbkkLuYl+sziZ3Zya5O7M4nl+/mkpwSBfuuec+vvOd++nTRy70CN9ltVpJSZnBlCnT2LQpky+//JzMzHTy8vbwm9+8RuywMSSMmsSQ+FEEBYXIhdEG0rVyFUf/GVVaWkJOzi62bt1MdvZm8vMPA2CxWBkYk8DopBnEjRjvViNRPL1boXHCrEfnqA7vw9PbwFHcoR2Kz59hW9YacrZt5NyZkwD4+fnTf+BQZk+fSmLiKGJiFEFBznkPSdeKD7Hb7Zw/f46jR4+g9V727s1lx64czp89ffk5AYFBxMSPIX5kEkOH30SXsAiDiYXwDF279WD6rQ8w7Zb7KTiZz+7tGezP3Ub+gd28dWA3AP7+/gwZEkNcXAJxcfFERw/kxhsHEBraxXB612i1kCulrMCbwEigEnhSa32gyfangGeAGuBXWuulTspqVFVVFefOneXcubOcPXv28venT58kVx/k3JmTVFZeuuI1IaFhDB6aSL8BMQxWI+kXHYu/f4Chn0AIz2axWOjVdyC9+g5k5tyHuFhWQlDlMXbv3kVu7h4OHNjH3r15LFr07Wtsth7ceOMA+vTpR48ePbDZvv3q2jWKsLAw/P09/3y2LT/BnUCw1nqCUioJ+A0wD0Ap1Qv4ATAWCAY2KqVWaa0rnRG2sPAM1dXV1NbWUldXR11d43/t1NXVUlvb9LGmX/XbqqurqKyspKqqiqqq+u/1kUJqaqqpqa6iurqKupoqLpSVEuRXw4ULZVy8eJELF8q4cOFCs7n8/QPoZuvNDT360N3Wh979BtHnxsFEde8p86AI4SRdwiMhPJLElAQSU6CmuppTJw5z6thBCgtOcLbgOGfPnCQ7ewvZ2Vua3U9wSBdCQsMI7RJOSGg4QcHBBAQGE92nG8HBIXTvHkltrYXg4BCCgoLw9/fH398fP7/G//pdfqzx8frH/LBarVgsVqxWK4GBgfTq1dspbdGWQj4JWAGgtc5SSo1tsu0mIL2hcFcqpQ4AI4DmW62DPvvsE/7wh986erctCgsLIywsnJ49e6FUFN2733DF1w033EDPnr3IPVkjK9ULYZh/QAD9o2PpHx17xeNVlRUUny+ktPgcpcVnKSk+R2nJOS5dLKO84etS+QUKTh6hpqb68ut2OCHjD3/4MrfffqfD99uWQh4BlDT5d61Syl9rXXOdbWVAi5MmtNRh35Jnn32KZ599qiMvdbrhw00nEE19Z5bZccdCuFpbTiNLgfCmr2ko4tfbFg4UOyibEEKINmhLIU8HbgVo6CPPabJtMzBZKRWslIoE4oDdDk8phBCiWa2OI28yamUEYAG+S31hP6C1/qph1MrT1H8ovKK1/ty5kYUQQjTl8huChBBCOJYMtRBCCA8nhVwIITycFHIhhPBwnn9vaicopUKAvwE9qB8D/5jWuvA6zwsFMoAfa61XuDal87WlHZRSr1N/c5g/8JbW+m2XB3UCmYKiTW3wPHB/wz+Xa61/4fqUztVaGzR5zjLgS631/7o+ZfN8/Yz8H4AcrfVk4EPgJ80870+AN18VbrEdlFLTgCFa6wnUF/MXlVJRro/pFJenoAB+TP0UFMAVU1BMBG4GXlVKuc80lY7TUhsMAh4CkoEkYLZSaoSRlM7VbBs08SvALX/vfb2QX55+APgamHn1E5RS/0r92fhOF+ZytdbaIRN4ouF7O+AHVOMdrpiCgvp5gxpdnoJCa10CNE5B4W1aaoNjwBytda3W2g4EAN44v29LbYBSaj5Qx7fvE7fiM10rSqnvAc9f9XAB304xcM30AkqpGUCM1voZpdRE56d0vo60g9a6AqhQSgUAH1DftdL8LGKexaFTUHioZttAa10NnFVKWYDXge1a631GUjpXs22glEoAHgTmAz8zkq4VPlPItdbvAu82fUwp9QXfTjFwvekFvgcMUEqlAkOB0Uqp01prZ8yn4xIdbAcaulIWAqla61edndOFZAqKltsApVQw8B71H2TPujibq7TUBo8CfYG1QDRQpZTKd6frZT5TyJvROP3AZuAWIK3pRq31g43fK6X+CnziyUW8BS22Q8PF0DXAb7TWH7k+nlOlA7cDnzYzBcV/NhSyILx3Copm26DhTPxLYK3W+jVD+Vyh2TbQWv+o8Xul1M+B0+5UxEEK+Z+BD5RSG4Eq6v98Qin1X8BCrfVmk+FcqMV2oP5i3yDgqYZRHADf1VofNhHWwRYBs5RSGTRMQaGUeoFvp6D4H+o/2KzAvzV0M3mbZtuA+ushU4EgpdQtDc9/SWudaSaq07T4e2A2WuvkFn0hhPBwvj5qRQghPJ4UciGE8HBSyIUQwsNJIRdCCA8nhVwIITycFHIhhPBwUsiFEMLD/X++jcKPNgEEbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Configurables\n",
    "num_sims = 5000\n",
    "replace = True  # False is permutation, True is sampling with replacement\n",
    "alpha = .01 # Significance level - we care only about order of magnitude really\n",
    "\n",
    "run_hyp_test(x1, x2, num_sims, replace, alpha)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}