NOTES 04.04 - HYPOTHESIS TESTING - COMPARING SAMPLE MEANS.ipynb

{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# 04.04 - COMPARING SAMPLES"]}, {"cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["endpoint https://m5knaekxo6.execute-api.us-west-2.amazonaws.com/dev-v0001/rlxmooc\n"]}, {"data": {"text/html": ["<h3>See <a href='https://m5knaekxo6.execute-api.us-west-2.amazonaws.com/dev-v0001/rlxmooc/web/login' target='_blank'>my courses and progress</a></h2>"], "text/plain": ["<IPython.core.display.HTML object>"]}, "execution_count": 26, "metadata": {}, "output_type": "execute_result"}], "source": ["!wget --no-cache -O init.py -q https://raw.githubusercontent.com/rramosp/20201.xai4eng/master/content/init.py\n", "import init; init.init(force_download=False); init.get_weblink()"]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": ["from scipy import stats\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from progressbar import progressbar as pbar\n", "%matplotlib inline"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Comparing populations. A two-sided test\n", "\n", "We have some exam scores for two groups of 10 students each"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### STEP 1: Define the model you want to challenge (the NULL  Hypothesis $H_0$)\n", "\n", "We define the $H_0$ as a model where there is no statistical difference between the two groups of students."]}, {"cell_type": "markdown", "metadata": {}, "source": ["### STEP 2: Define your REAL WORLD sample and your test statistic\n", "\n", "- **SAMPLE**: we sample two groups of 10 students exam scores each. We **KNOW** there is a difference (`mean_diff`) because we create the use case in such way. We want to see if the test detects it.\n", "\n", "    - Can it detect it with larger/smaller group sizes?\n", "    - Can it detect it with larger/smaller mean differences?"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["group A  [  96.08 101.76  89.45  88.81 106.66 105.94  95.99 108.45  94.67 114.32 ]  mean 100.21\n", "group B  [ 129.23 127.05 115.34 105.17 117.75 121.96 103.69 120.42 118.77  89.56 ]  mean 114.89\n"]}], "source": ["real_world_sample_size = 10\n", "mu, sigma = 100, 10\n", "\n", "def generate_AB(mu=mu, sigma=sigma, real_world_sample_size=real_world_sample_size, mean_diff=5):\n", "    A = np.random.normal(loc=mu, scale=sigma, size=real_world_sample_size)\n", "    B = np.random.normal(loc=mu+mean_diff, scale=sigma, size=real_world_sample_size)\n", "    return A, B\n", "\n", "A,B = generate_AB()\n", "print (\"group A  [\", \" \".join([\"%6.2f\"%i for i in A]), \"]  mean %6.2f\"%np.mean(A))\n", "print (\"group B  [\", \" \".join([\"%6.2f\"%i for i in B]), \"]  mean %6.2f\"%np.mean(B))\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["if we look at these two sets of scores, how sure can we be that their respective populations are different (different $\\mu$)?\n", "\n", "\n", "#### Define a test statistic to measure the mean difference of the samples\n", "\n", "$$\\text{ttest}(A,B) = \\frac{\\bar{A}-\\bar{B}}{\\sqrt{\\frac{S_A^2}{N_A}+\\frac{S_B^2}{N_B}}}$$\n", "\n", "where:\n", "\n", "- $\\bar{A}$, $\\bar{B}$ is the mean of the two real world samples\n", "- $S_A$, $S_B$ is the standard deviation of the two real world samples\n", "- $N_A$, $N_B$ is the size (number of elementso) of the two real world samples\n", "\n", "The **denominator** factor is included so that we can compare the means with different standard deviations.\n", "\n", "\n", "### STEP 3: We simulate $H_0$, and understand its distribution of $\\text{ttest}$ \n", "\n", "individual student scores **NEED NOT** come from a normal disitribution. Regardless their origin, the distribution of the test statistic will be normal."]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100% (3000 of 3000) |####################| Elapsed Time: 0:00:01 Time:  0:00:01\n", "100% (3000 of 3000) |####################| Elapsed Time: 0:00:01 Time:  0:00:01\n"]}], "source": ["ttest = lambda a,b: (np.mean(a)-np.mean(b))/np.sqrt(np.std(a)**2/len(a)+np.std(b)**2/len(b))\n", "\n", "n = 3000\n", "\n", "a = np.r_[[stats.norm(loc=mu, scale=sigma).rvs(real_world_sample_size) for _ in pbar(range(n))]]\n", "b = np.r_[[stats.norm(loc=mu, scale=sigma).rvs(real_world_sample_size) for _ in pbar(range(n))]]\n", "\n", "t = np.r_[[ttest(i,j) for i,j in zip(a,b)]]"]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3gV5bn38e9NQIMaxQNSCCj4Fk8gRI1QtwKxgqJUwIqiKEVbpPrCVrFq8VCxotbag6d6gFqlVRStdGu0WEXrQrSoQAUBLYpANcB+VQ6FCAkJ3O8faxIXyaxkBTLMCv4+15WLNTPPzPyyWMmdmWfmGXN3REREamoWdwAREclOKhAiIhJKBUJEREKpQIiISCgVCBERCdU87gCN5aCDDvKOHTvGtv+vvvqKvffee5fvd8mSJQAcccQRocvjylWXbMwE2ZkrGzNB+lz1fR6j1NTeqzilZpo3b96X7t46tKG77xZfxx9/vMfp9ddfj2W/ffr08T59+qRdHleuumRjJvfszJWNmdzT56rv8xilpvZexSk1EzDX0/xe1SkmEREJpQIhIiKhVCBERCTUbtNJLfJNUFFRQUlJCWVlZbHm2G+//fjwww9rzR8/fjxA6LKopcsUt2zJlZubS/v27WnRokXG66hAiDQhJSUl5OXl0bFjR8wsthwbN24kLy+v1vxmzZInJeK4iildprhlQy53Z82aNZSUlNCpU6eM19MpJpEmpKysjAMPPDDW4iBNj5lx4IEHNvjIUwVCpIlRcZAdsSOfGxUIEREJpQIhIhn74osvOPnkk+nZsyfPPfdc9fxBgwaxatWq6ulbbrmF/Px8CgoKKCgoYPr06bW2tWLFCp588sldknvkyJF88MEHteZPnjyZMWPGNGhbHTt25Msvv6w1/8wzz2T9+vV1rjt58uTt3qeo3HLLLfz617/e6e2ok1q+Me6e8VGdy/PLyqvbjO13+K6I1OQ89dRTXHbZZfTr14+hQ4cyePBgXnjhBY499ljatWtXPdQGwNixY7nmmmvSbquqQAwbNizy3I888kjk+wgrgjVNnjyZrl270q5du4y3W1lZSfPm8fyqjvQIwsz6m9kSM1tqZuPqaHeOmbmZFabMuz5Yb4mZnR5lThHJTIsWLdi0aRPl5eXk5ORQWVnJPffcw3XXXdfgbY0bN45Zs2ZRUFDA3XffzYABA3j//fcBOPbYY7n11lsBuPnmm/n973+Pu3PttdfStWtXjjnmGJ5++ula2/zqq68YMGAA3bt3p2vXrtVtioqKmDt3LgCPPfYYhx9+OD169OCtt96qXveLL77gnHPO4YQTTuCEE06oXrZmzRpOO+00unTpwsiRI/E0T+GsOrJYsWIFRx11FJdeeildunRh0KBBbN68mWeffZa5c+dy4YUXUlBQwObNm5k3bx59+vTh+OOP5/TTT2f16tXVea+66ioKCwu5/fbbOfTQQ9m2bVv199ihQwcqKir4/e9/zwknnED37t0555xz2LRpU4P/H+oSWYEwsxzgAeAM4GjgAjM7OqRdHnAl8E7KvKOB84EuQH/gwWB7IhKjYcOG8fzzzzN48GBuuOEGHnzwQYYPH85ee+1Vq+3vfvc7unXrxg9/+EPWrVtXa/mdd95Jr169mD9/PmPHjqVXr17MmjWL//znPzRv3rz6F/SsWbPo3bs3f/nLX5g/fz4LFizg1Vdf5dprr63+hVrlb3/7G+3atWPBggUsWrSI/v37b7d89erVjB8/nrfeeos333xzu9NOV155JWPHjmXOnDlMmzaNkSNHAvDzn/+ck08+mcWLF3P22Wfz6aef1vs+ffzxx4wePZrFixfTqlUrpk2bxpAhQygsLGTKlCnMnz+f5s2b89///d88++yzzJs3jx/+8IfceOON1dvYsmULc+fOZfz48RQUFDBz5kwAXnzxRU4//XRatGjB97//febMmcOCBQs46qij+MMf/lBvtoaI8giiB7DU3Ze5+xZgKjAopN0E4JdA6vVXg4Cp7l7u7suBpcH2RCRG++23H3/961+ZOXMmxx13HC+88AJDhgzh0ksvZciQIbz33nsAXH755XzyySfMnz+ftm3b8pOf/KTebffq1Ys33niDt956iwEDBlBaWsqmTZtYvnw5RxxxBG+++SYXXHABOTk5tGnThj59+jBnzpzttnHMMccwY8YMfvrTnzJr1iz222+/7Za/8847FBUV0bp1a/bYYw+GDh1avezVV19lzJgxFBQUMHDgQDZs2EBpaSlvvPEGF110EQADBgxg//33r/d76dSpEwUFBQAUFBSwYsWKWm2WLFnCokWL6NevHwUFBdx2222UlJRUL0/NNnTo0OqjoalTp1YvW7RoEb169eKYY45hypQpLF68uN5sDRHlia184LOU6RKgZ2oDMzsO6ODufzWza2us+3aNdfNr7sDMRgGjANq0aUMikWic5DugtLQ0lv1XdYql23dcueoSV6b8svI6l7fYVk5+2XIAEonoOxIzUfO92m+//di4cWN8gQJbt27lZz/7GWPHjuWxxx7j+OOPZ/DgwZxzzjlMmjSJvfbaq/p0xwUXXMB5551XK/emTZuorKysnn/kkUfy7rvv0r59e0455RRWrVrF/fffT/fu3dm4cSNbtmyhrKysun1FRQWbN2+unt66dStt27Zl5syZvPLKK1x//fX06dOHcePGsXXrVr766is2b95MRUVF9TplZWVs2bKFjRs3snXrVmbMmEFubm51Rndn27ZtlJaWVq/j7pSWlrLnnntu9/1UzS8tLaVFixbV7c2Mr776qnofVa9LS0s58sgjee2117bbTlW7qtcAp5xyCtdffz3//ve/mTt3LieccAIbN25kxIgRPPnkk9UFYtasWWzcuJHy8vLtMlQpKysjkUhk/DMYWye1mTUDfgtcvKPbcPdJwCSAwsJCLyoqapRsOyKRSBDH/lu1agWQdt9x5apLXJnq76Rezsrc5F2m5xVlRyd1zffqww8/jP2uXID33nuPzz//nDPPPJN7772Xvffem3333ZctW7bQvHlzSktLadu2LQAzZsygW7dutXK3adOGzZs3bzf/0EMP5fnnn+e2227jq6++4pprruGaa64hLy+PU089lYkTJ/LjH/+YtWvXMnv2bO65557q9Tdu3MjGjRtp06YNl156KW3btuWRRx4hLy+PnJwc9t57b4qKihg3bhxbtmxh33335YUXXqB79+7k5eVx+umnM3nyZK69Nvm36vz58ykoKKCoqIji4mJuuukmXnrpJdavX88+++xT6/sxM/bZZx8geUd51fJmzZqx5557kpeXR6tWrdi2bRt5eXkcd9xxrF27lkWLFnHiiSdSUVHBRx99RJcuXarzVm0jLy+PHj16cNNNN3HWWWdV/9yXlpby7W9/m9zcXKZNm0Z+fj55eXnsueee1ftMlZuby7HHHpvxz2CUBWIl0CFlun0wr0oe0BVIBDdwfAsoNrOBGawrIjGaMGECv/zlL4HkEcLgwYO58847+fGPfwzAddddx/z58zEzOnbsyMSJE2tto1u3buTk5NC9e3cuvvji6n6I1157jZYtW9KrVy9KSkro1asXAGeffTazZ8+me/fumBl33XUX3/rWt7bb5sKFC7n22mtp1qwZLVq04KGHHtpuedu2bbnllls48cQTadWqVfVpIID77ruP0aNH061bNyorK+nduzcPP/ww48eP54ILLqBLly7813/9F4cccsgOv28XX3wxl112GS1btmT27Nk8++yzXHHFFfznP/+hsrKSq666ii5duoSuO3ToUM4999zt/vKfMGECPXv2pHXr1vTs2bPRjy4tXY/8Tm/YrDnwEXAqyV/uc4Bh7h56kszMEsA17j7XzLoAT5Lsd2gHvAZ0dvet6fZXWFjoVVcpxCGuv4qr9pnucFFHEF9ryBFEtlzmGnYEcdRRR8UXKJBufKE4nyiXDWMehcmmXFWfn9TPlZnNc/fCsPaRHUG4e6WZjQFeBnKAR919sZndSvIJRsV1rLvYzJ4BPgAqgdF1FQcREWl8kfZBuPt0YHqNeTenaVtUY/p24PbIwok0cVdddRXz589v1G0WFBRwzz33NOo2penSUBsiIhJKQ22INFFx/aW/efNmzjjjDGbOnMnChQu5/PLL2bBhAzk5OVxyySWceeaZtdapujooPz95tfqYMWMYOXIkS5YsYdiwYVRUVDBx4kROPPFEKisr6d+/P8XFxaE34NXl6aef5vbbb2fr1q1873vfq+5IHzt2LK+//jqQvLz2888/Dx03qaioiNWrV9OyZUsAXnnlFQ4++GDuv/9+Jk6cyCGHHMJzzz3HHnvswZtvvsm0adO4++67geSd2MOHD+dvf/tbgzJnMxUIEWmQRx99lLPOOoucnBz22msv/vSnP9G5c2dWrVpFQUEBJ598cuh6Q4cO5Xe/+9128yZOnMi9995Lx44dufLKK5k2bRoPPfQQF110UYOLw5o1a7j22muZN28erVu3ZsSIEbz22muceuqp1b/EAe6///7qG/rCTJkyhcLCwlrz3n//fe644w5efvllvve97zFhwgSeeuqp6jatW7embdu2vPXWW5x00kkNyp6tdIpJRBpkypQpDBgwAIDDDz+czp07A9CuXTsOOOAA1q5dm/G2qsZ22rRpEy1atGD9+vW88MIL/OAHP2hwrhUrVtC5c2dat24NQN++fZk2bVqtdk899RQXXHBBg7bt7lRUVFTnfOKJJzjjjDM44IADtms3ePBgpkyZ0uDs2UoFQkQytmXLFpYtW8ahhx5aa9m7775LRUVF2vsEpk2bRrdu3RgyZAiffZYcZGH06NHccccdjBgxghtuuIEJEyZwww03VD+6tCEOO+wwlixZwooVK6isrOS5556r3k+Vf//73yxfvpzvfve7abdzySWXUFBQwIQJE6oH5hszZgzf+c53+PTTTznppJN47LHHGD16dK11CwsLmTVrVoOzZysVCBHJ2Jdffll9F2+q1atXM3z4cO64447QX+5nnXUWK1as4P3336dfv36MGDECgEMOOYREIsHs2bPZa6+9KCkp4aijjmL48OEMHTqUjz6q+96VVPvvvz8PPfQQQ4cOpVevXnTs2JGcnO3H+Jw6dSpDhgypNb/KlClTWLhwIbNmzWLWrFk8/vjjAAwfPpz33nuPJ554grvvvpsrrriCl156iSFDhjB27NjqkVYPPvjgXfK8h11FBUJEMtayZctazzXesGEDAwYM4Pbbb9/uzuRUBx54YPXYRSNHjmTevHm12tx4443cdttt3HfffYwcOZK77rqLn//859u1+eyzz6ofQvTwww/X2sZZZ53FO++8w+zZszniiCM4/PDtb3icOnVqnaeXqjrR8/LyGDZsGO++++52y1etWsW7777L4MGD+c1vfsPTTz9Nq1atqsdTKisrq+7g3h2ok1pEMrb//vuzdetWysrKyMvLY8uWLZx99tn84Ac/YMiQIds9MCjV6tWrq8dmKi4urnU3+MyZM2nXrh2dO3dm06ZNNGvWjGbNmtV6vkGHDh3qvPfj888/5+CDD2bdunU8+OCDPPPMM9XL/vWvf7Fu3TpOPPHE0HUrKytZv349Bx10EBUVFbz44ov07dt3uzY/+9nPqp9TsXnzZsxsu5wfffQRXbt2TZuvqVGBEJEGOe2005g9ezYDBw7kmWee4Y033mDNmjVMnjyZ8vJyfvGLX3DEEUdw8803U1hYyMCBA7nvvvsoLi6mefPmHHDAAUyePLl6e+7ObbfdVj2c9ahRo7jwwguprKysNZZSfa688koWLFgAJB80lHoEMXXqVM4//3yCsd+qFRQUMH/+fMrLyzn99NOpqKhg69at9O3bl0svvbS6XdWVT8cddxyQfDbGMcccQ4cOHaofmPT6669Xd+DvDiIbi2lX01hMidDlGovpaxqLqXH885//5K677mLq1Km1ln3Tx2Lq3bs3zz///HbPjMiGXFUaOhaT+iBEpEGOO+44evfuXf3MAkn64osvuPrqqzN6oFBToQIhIg02fPjwtFcCfVO1bt2awYMHxx2jUalAiDQxu8tpYdm1duRzowIh0oTk5uayZs0aFQlpEHdnzZo12z1ONRO6ikmkCWnfvj0lJSV88cUXseYoKysL/WXzv//7vwDVN47tSukyxS1bcuXm5tK+ffsGraMCIdKEtGjRgk6dOsUdg0QiwbHHHltr/uWXX169fFdLlylu2ZorE5GeYjKz/ma2xMyWmtm4kOWXmdlCM5tvZm+a2dHB/I5mtjmYP9/Mat8yKSIikYrsCMLMcoAHgH5ACTDHzIrd/YOUZk+6+8NB+4HAb4H+wbJP3D38vn2RFPXd3yAiOybKI4gewFJ3X+buW4CpwKDUBu6+IWVyb0A9byIiWSLKPoh8IHWs3RKgZ81GZjYauBrYA0gdg7eTmb0HbABucvdaY+ia2ShgFECbNm1iOe9ZpbS0NJb9Vz0VK92+48pVl8bOlF9W3ijbabGtnPyy5QAkEtkxImc2/v9B+lz1fR6j1NTeqzhlmin2Tmp3fwB4wMyGATcBI4DVwCHuvsbMjgeeM7MuNY44cPdJwCRIDrUR55AScQ0fUTX0crp9fxOG2misU0ypQ22cV5SdQ21ki3S56vs8RqmpvVdxyjRTlKeYVgIdUqbbB/PSmQoMBnD3cndfE7yeB3wCZMdPrIjIN0SURxBzgM5m1olkYTgfGJbawMw6u/vHweQA4ONgfmtgrbtvNbPDgM7Asgizimwn06OSbBnUTyQKkRUId680szHAy0AO8Ki7LzazW4G57l4MjDGzvkAFsI7k6SWA3sCtZlYBbAMuc/fMH3QrIiI7LdI+CHefDkyvMe/mlNdXpllvGlD7aeMiIrLLaCwmEREJpQIhIiKhVCBERCSUCoSIiIRSgRARkVAqECIiEkoFQkREQqlAiIhIqNgH6xNpyjQkh+zOdAQhIiKhVCBERCSUCoSIiIRSgRARkVAqECIiEkoFQkREQqlAiIhIqEgLhJn1N7MlZrbUzMaFLL/MzBaa2Xwze9PMjk5Zdn2w3hIzOz3KnCIiUltkN8qZWQ7wANAPKAHmmFmxu3+Q0uxJd384aD8Q+C3QPygU5wNdgHbAq2Z2uLtvjSqvZJ9Mb0ITkWhEeQTRA1jq7svcfQswFRiU2sDdN6RM7g148HoQMNXdy919ObA02J6IiOwiUQ61kQ98ljJdAvSs2cjMRgNXA3sA301Z9+0a6+ZHE1NERMLEPhaTuz8APGBmw4CbgBGZrmtmo4BRAG3atCGRSESSMROlpaWx7H/9+vUAafcdV666ZJopv6w8+jApWmwrJ79seSTbTiRW7dB62fj/B+lz1fd5jFJTe6/ilGmmKAvESqBDynT7YF46U4GHGrKuu08CJgEUFhZ6UVHRTsTdOYlEgjj236pVK4C0+44rV10yzbSr+yDyy5azMrdTJNs+r2jHBuvLxv8/SJ+rvs9jlJraexWnTDNF2QcxB+hsZp3MbA+Snc7FqQ3MrHPK5ADg4+B1MXC+me1pZp2AzsC7EWYVEZEaIjuCcPdKMxsDvAzkAI+6+2IzuxWY6+7FwBgz6wtUAOsITi8F7Z4BPgAqgdG6gklEZNeKtA/C3acD02vMuznl9ZV1rHs7cHt06UREpC66k1pEREKpQIiISCgVCBERCaUCISIioVQgREQkVOx3Uot8E2R609/Yfjt2Q51IFHQEISIioVQgREQklAqEiIiEUoEQEZFQKhAiIhJKBUJEREKpQIiISCgVCBERCZVRgTCzv5jZADNTQRER+YbI9E7qB4FLgPvM7M/AY+6+JLpYsjv7fGP5Ln+cqIg0XEZHBO7+qrtfCBwHrABeNbN/mNklZtYiyoAiIhKPjE8ZmdmBwMXASOA94F6SBWNGHev0N7MlZrbUzMaFLL/azD4ws/fN7DUzOzRl2VYzmx98FddcV0REopXRKSYz+x/gCOBx4Cx3Xx0setrM5qZZJwd4AOgHlABzzKzY3T9IafYeUOjum8zscuAuYGiwbLO7FzT4OxIRkUaRaR/E74PnS1czsz3dvdzdC9Os0wNY6u7LgvZTgUFAdYFw99dT2r8NXJRxchERiVSmBeI2YHqNebNJnmJKJx/4LGW6BOhZR/sfAS+lTOcGRyeVwJ3u/lzNFcxsFDAKoE2bNiQSiTo2H63S0tJY9r9+/XqAtPuOK1ddWmwrJ79sedwxasmGXInEqu2ms/H/D9Lnqu/zGKWm9l7FKdNMdRYIM/sWyV/0Lc3sWMCCRfsCe+1kxtT9XAQUAn1SZh/q7ivN7DDg72a20N0/SV3P3ScBkwAKCwu9qKiosSI1WCKRII79t2rVCiDtvuPKVZdnXniZlbmd4o5RS37Z8thznVe0/fMgsvH/D9Lnqu/zGKWm9l7FKdNM9R1BnE6yY7o98NuU+RuBG+pZdyXQIWW6fTBvO2bWF7gR6OPu5VXz3X1l8O8yM0sAxwKf1FxfRESiUWeBcPc/An80s3PcfVoDtz0H6GxmnUgWhvOBYakNgqOSiUB/d/88Zf7+wCZ3Lzezg4CTSHZgi4jILlLfKaaL3P0JoKOZXV1zubv/NmS1qmWVZjYGeBnIAR5198Vmdisw192LgV8B+wB/NjOAT919IHAUMNHMtpG8FPfOGlc/iYhIxOo7xbR38O8+O7Lx4Mqn6TXm3Zzyum+a9f4BHLMj+xQRkcZR3ymmicG/P981cUREJFvUd4rpvrqWu/sVjRtHRESyRX2nmObtkhQiIpJ1MrmKSUR2kZqj3OaXhY98O7bf4bXmiTS2+k4x3ePuV5nZC4DXXB5ccSQiIruh+k4xPR78++uog4iISHap7xTTvODfmWa2B3AkySOJJe6+ZRfkExGRmGQ63PcA4GGSQ10Y0MnMfuzuL9W9poiINFWZjub6G+AUd18KYGb/B/gr24++KiIiu5FMnyi3sao4BJaRHLBPRER2U/VdxfT94OVcM5sOPEOyD+JckoPxiYjIbqq+U0xnpbz+f3z9vIYvgJaRJBIRkaxQ31VMl+yqICIikl0yvYopl+QjQbsAuVXz3f2HEeUSEZGYZdpJ/TjwLZJPmJtJ8ulw6qQWEdmNZVogvu3uPwO+CsZnGgD0jC6WiIjELdMCURH8u97MugL7AQfXt5KZ9TezJWa21MzGhSy/2sw+MLP3zew1Mzs0ZdkIM/s4+BqRYU4REWkkmRaIScFzon8GFAMfAL+sawUzywEeAM4AjgYuMLOjazR7Dyh0927AswTPnTazA4DxJI9SegDjg/2LiMguklGBcPdH3H2du89098Pc/eCqp83VoQew1N2XBeM2TQUG1dju6+6+KZh8m2TfBiT7Oma4+1p3XwfMAPpn+k2JiMjOy/QqpgOBW4CTSN4oNwuY4O5r6lgtH/gsZbqEuvstfsTXQ3eErZsfkmsUMAqgTZs2JBKJur6NSJWWlsay//Xr1wOk3XdcuerSYls5+WXL445RSzbmSpcpkVgVQ5qvpftc1fd5jFI2ftYhO3NlminTsZimAm8A5wTTFwJPA313JFxNZnYRUMjXN+JlxN0nAZMACgsLvaioqDHi7JBEIkEc+2/VqhVA2n3Hlasuz7zwMitzO8Udo5b8suVZlytdpvOK4n1gULrPVX2fxyhl42cdsjNXppkyLRBt3X1CyvRtZja0nnVWAh1SptsH87ZjZn2BG4E+7l6esm5RjXUTGWaVmIQ9+SxMrUNBEclKmXZSv2Jm55tZs+DrPODletaZA3Q2s07BsyTOJ9nBXc3MjgUmAgPd/fOURS8Dp5nZ/kHn9GkZ7E9ERBpRfYP1bSTZ52DAVcATwaJmQClwTbp13b3SzMaQ/MWeAzzq7ovN7FZgrrsXA78C9gH+bGYAn7r7QHdfa2YT+HpAwFvdfe2OfpMiItJw9Y3FlLczG3f36cD0GvNuTnmdtg/D3R8FHt2Z/YuIyI7LtA8CMxsI9A4mE+7+YjSRREQkG2TUB2FmdwJXkrxB7gPgSjP7RZTBREQkXpkeQZwJFLj7NgAz+yPJu6CvjyqYiIjEK9OrmABapbzer7GDiIhIdsn0COIO4D0ze53kFU29gVqD74mIyO6j3gJhZs2AbcB3gBOC2T919/+NMpiIiMSr3gLh7tvM7Dp3f4YaN7qJiMjuK9M+iFfN7Boz62BmB1R9RZpMRERilWkfxFCSd1T/3xrzD2vcOCIiki0yLRBHkywOJ/P1cN8PRxVKRETil2mB+COwAbgvmB4WzDsvilAiIhK/TAtEV3dPfVzo62b2QRSBRKR+mQ6tPrZfvM+NkKYt007qf5rZd6omzKwnMDeaSCIikg0yPYI4HviHmX0aTB8CLDGzhYC7e7dI0omISGwyLRD9I00hWS3T0xkisnvJqEC4+7+jDiIiItkl4+dB7Agz6w/cS/KJco+4+501lvcG7gG6Aee7+7Mpy7YCC4PJT919YJRZRXZH6syWnRFZgTCzHOABoB9QAswxs2J3T7366VPgYsIfXbrZ3QuiyiciInWL8giiB7DU3ZcBmNlUYBDJBw4B4O4rgmXbIswhIiI7oCHPg2iofOCzlOmSYF6mcs1srpm9bWaDGzeaiIjUJ9I+iJ10qLuvNLPDgL+b2UJ3/yS1gZmNAkYBtGnThkQiEUPMpNLS0lj2v379eoC0+26MXPll5Tu1fk0ttpWTX7a8UbfZGLIx167KlEisalD7dJ+r+j6PUYrrZ7A+2Zgr00xRFoiVQIeU6fbBvIy4+8rg32VmlgCOBT6p0WYSMAmgsLDQi4qKdi7xTkgkEsSx/1atkg/6S7fvxsjV2Je55pctZ2Vup0bdZmPIxly7KtN5RQ3rpE73uarv8xiluH4G65ONuTLNFOUppjlAZzPrZGZ7AOeT4fMkzGx/M9szeH0QcBIpfRciIhK9yAqEu1cCY4CXgQ+BZ9x9sZndamYDAczsBDMrAc4FJprZ4mD1o4C5ZrYAeB24s8bVTyIiErFI+yDcfTowvca8m1NezyF56qnmev8Ajokym4iI1C3KU0wiItKEqUCIiEgoFQgREQmlAiEiIqFUIEREJJQKhIiIhFKBEBGRUCoQIiISSgVCRERCqUCIiEgoFQgREUxJX1UAAA1zSURBVAmlAiEiIqFUIEREJJQKhIiIhFKBEBGRUCoQIiISKtIHBolI09CQ546P7dew51dL0xXpEYSZ9TezJWa21MzGhSzvbWb/NLNKMxtSY9kIM/s4+BoRZU4REaktsgJhZjnAA8AZwNHABWZ2dI1mnwIXA0/WWPcAYDzQE+gBjDez/aPKKiIitUV5BNEDWOruy9x9CzAVGJTawN1XuPv7wLYa654OzHD3te6+DpgB9I8wq4iI1BBlH0Q+8FnKdAnJI4IdXTe/ZiMzGwWMAmjTpg2JRGKHgjaG0tLSWPa/fv16gLT7boxc+WXlO7V+TS22lZNftrxRt9kYsjFXNmZKJFal/VzV93mMUlw/g/XJxlyZZmrSndTuPgmYBFBYWOhFRUWxZUkkEsSx/1atWgGk3XdduTLumMzdgWB1yC9bzsrcTo270UaQjbmyMdN5RYen/VzV93mMUlw/g/XJxlyZZoryFNNKoEPKdPtgXtTriohII4iyQMwBOptZJzPbAzgfKM5w3ZeB08xs/6Bz+rRgnoiI7CKRFQh3rwTGkPzF/iHwjLsvNrNbzWwggJmdYGYlwLnARDNbHKy7FphAssjMAW4N5omIyC4SaR+Eu08HpteYd3PK6zkkTx+Frfso8GiU+UREJD0NtSEiIqFUIEREJJQKhIiIhFKBEBGRUCoQIiISSgVCRERCqUCIiEioJj0Wk4jsenfP+Ij8svLQsbxK1m2ubqMHCzV9OoIQEZFQKhAiIhJKBUJEREKpQIiISCgVCBERCaUCISIioVQgREQklAqEiIiEirRAmFl/M1tiZkvNbFzI8j3N7Olg+Ttm1jGY39HMNpvZ/ODr4ShziohIbZHdSW1mOcADQD+gBJhjZsXu/kFKsx8B69z922Z2PvBLYGiw7BN3L4gqn4iI1C3KI4gewFJ3X+buW4CpwKAabQYBfwxePwucamYWYSYREclQlGMx5QOfpUyXAD3TtXH3SjP7D3BgsKyTmb0HbABucvdZNXdgZqOAUQBt2rQhkUg06jfQEKWlpbHsf/369QBp911Xrvyy8ohS1a3FtnLyy5bHsu+6ZGOubMwE6XPtuS05FlN+2XISiVW7NFNcP4P1ycZcmWbK1sH6VgOHuPsaMzseeM7Murj7htRG7j4JmARQWFjoRUVFuz5pIJFIEMf+W7VqBZB233XlChtsbVfIL1vOytxOsey7LtmYKxszQfpc5c1aArAytxPnFe3awfri+hmsTzbmyjRTlKeYVgIdUqbbB/NC25hZc2A/YI27l7v7GgB3nwd8AmhoSBGRXSjKI4g5QGcz60SyEJwPDKvRphgYAcwGhgB/d3c3s9bAWnffamaHAZ2BZRFmzXrp/tpPHV45TLphmUVE6hNZgQj6FMYALwM5wKPuvtjMbgXmunsx8AfgcTNbCqwlWUQAegO3mlkFsA24zN3XRpVVRERqi7QPwt2nA9NrzLs55XUZcG7IetOAaVFmExGRuulOahERCaUCISIioVQgREQklAqEiIiEytYb5b4xdAmqiGQrFQgRiUSmf/yM7ad7YLOVTjGJiEgoFQgREQmlAiEiIqHUByEisVJfRfbSEYSIiIRSgRARkVA6xRQR3d8gIk2djiBERCSUjiAaSA/mEYlHfT9fVT+D6sxuPDqCEBGRUCoQIiISKtJTTGbWH7iX5CNHH3H3O2ss3xP4E3A8sAYY6u4rgmXXAz8CtgJXuPvLUWbV6SERke1FViDMLAd4AOgHlABzzKzY3T9IafYjYJ27f9vMzgd+CQw1s6NJPp+6C9AOeNXMDnf3rVHlFZHdQ2P/sfdN7tOI8giiB7DU3ZcBmNlUYBCQWiAGAbcEr58FfmdmFsyf6u7lwHIzWxpsb3aEeUVEatnZgrMrLmCJqohFWSDygc9SpkuAnunauHulmf0HODCY/3aNdfNr7sDMRgGjgslSM1vSONF3yEHAl3Ht/OrTjki3KNZcaWRjJsjOXNmYCerJVcfnMUpN8r1qDFc3fJXUTIema9SkL3N190nApLhzAJjZXHcvjDtHTdmYKxszQXbmysZMkJ25sjETZGeuTDNFeRXTSqBDynT7YF5oGzNrDuxHsrM6k3VFRCRCURaIOUBnM+tkZnuQ7HQurtGmGBgRvB4C/N3dPZh/vpntaWadgM7AuxFmFRGRGiI7xRT0KYwBXiZ5meuj7r7YzG4F5rp7MfAH4PGgE3otySJC0O4Zkh3alcDoJnAFU1ac6gqRjbmyMRNkZ65szATZmSsbM0F25sookyX/YBcREdme7qQWEZFQKhAiIhJKBaIRmdkEM3vfzOab2Stm1i7uTABm9isz+1eQ7X/MrFUWZDrXzBab2TYzi/USQDPrb2ZLzGypmY2LM0sVM3vUzD43s0VxZ6liZh3M7HUz+yD4v7sy7kwAZpZrZu+a2YIg18/jzlTFzHLM7D0zezHuLFXMbIWZLQx+T82tq60KROP6lbt3c/cC4EXg5rgDBWYAXd29G/ARcH3MeQAWAd8H3ogzRMqQMGcARwMXBEO9xG0y0D/uEDVUAj9x96OB7wCjs+S9Kge+6+7dgQKgv5l9J+ZMVa4EPow7RIhT3L2gvnshVCAakbtvSJncG8iKKwDc/RV3rwwm3yZ5X0ms3P1Dd4/zzvcq1UPCuPsWoGpImFi5+xskr+zLGu6+2t3/GbzeSPIXX60RDnY1TyoNJlsEX7H/7JlZe2AA8EjcWXaUCkQjM7Pbzewz4EKy5wgi1Q+Bl+IOkUXChoSJ/ZdetjOzjsCxwDvxJkkKTuXMBz4HZrh7NuS6B7gO2BZ3kBoceMXM5gXDFaWlAtFAZvaqmS0K+RoE4O43unsHYAowJltyBW1uJHmaYEq2ZJKmx8z2AaYBV9U4ao6Nu28NTu22B3qYWdc485jZ94DP3X1enDnSONndjyN5WnW0mfVO17BJj8UUB3fvm2HTKcB0YHyEcarVl8vMLga+B5zqu+jmlwa8V3HSsC4NYGYtSBaHKe7+l7jz1OTu683sdZL9N3F28J8EDDSzM4FcYF8ze8LdL4oxEwDuvjL493Mz+x+Sp1lD+wJ1BNGIzKxzyuQg4F9xZUkVPLjpOmCgu2+KO0+WyWRIGAGCofj/AHzo7r+NO08VM2tddWWembUk+QyaWH/23P16d2/v7h1Jfqb+ng3Fwcz2NrO8qtfAadRRSFUgGtedwSmU90m+8VlxGSDwOyAPmBFc2vZw3IHM7GwzKwFOBP5qZpE+MTCdoPO+akiYD4Fn3H1xHFlSmdlTJJ9/coSZlZjZj+LORPKv4uHAd4PP0fzgL+S4tQVeD37u5pDsg8iay0qzTBvgTTNbQHJ8u7+6+9/SNdZQGyIiEkpHECIiEkoFQkREQqlAiIhIKBUIEREJpQIhIiKhVCCkSTGzW8zsmuD1rWbWN3jdKxjJc76ZtQxGsF1sZr+KN3HDmdlTwci7Y+POIt9supNamix3Tx3r6kLgF+7+BEAwxswBmT6q1syapwxoGBsz+xZwgrt/O+4sIjqCkKxnZjea2Udm9iZwRMr8yWY2xMxGAucBE8xsipkVA/sA88xsaHCn7TQzmxN8nRSsf4uZPW5mb5F8Nnpd7R41s4SZLTOzK1Iy/CD4a3+BmT0ezAvdTo3vKdfMHrPkuPzvmdkpwaJXgPzgSKhXjXUmm9lDZvZ2kKMoyPWhmU1OaXeamc02s3+a2Z+DsZMws5uDPIvMbFJwZzTB9/VLSz5T4aOq/ZpZl2De/OB7TB0pQL4J3F1f+sraL+B4YCGwF7AvsBS4Jlg2GRhS83UwXZry+kmSA5QBHEJyqAiAW4B5QMsM2v0D2BM4CFhDckjpLiSfr3FQ0O6AurZT4/v6CfBo8PpI4FOSY/Z0BBaleS8mkxyO3EgO5bIBOIbkH3rzSD4L4SCS4+rsHazzU+Dm1HzB68eBs4LXCeA3weszgVeD1/cDFwav96h6n/T1zfnSKSbJdr2A//FgDKng6KCh+gJHB38wQ3LgtH2C18XuvjmDdn9193Kg3Mw+JzlkwXeBP7v7lwDuvrau7fjXzywAOJnkL2Dc/V9m9m/gcJK/9Ovygru7mS0E/p+7LwQws8Uki0t7kg8+eivY/x4kh+wAOMXMriNZbA8AFgMvBMuqBt6bF2yHYL0bLflcg7+4+8f1ZJPdjAqEfBM0A77j7mWpM4NfoF9l2K48ZdZW6v7ZCd1OI6nKsa1Gpm1Bpq0kxyK6IHUlM8sFHgQK3f0zM7uF5BFLze1Wf2/u/qSZvUPyoTfTzezH7v73Rv5+JIupD0Ky3RvA4ODKpDzgrB3YxivAf1dNmFnBTrar8nfgXDM7MGh/QAO2M4tkxzpmdjjJU1GN8YS9t4GTzOzbwbb3DrZfVQy+DI6KhtS3ITM7DFjm7vcBzwPdGiGfNCEqEJLVPPmIy6eBBSSfhDdnBzZzBVAYdLR+AFy2k+2qsi0GbgdmWnJ0zKohsDPZzoNAs+BU0dPAxcEprJ3i7l8AFwNPWXJ009nAke6+Hvg9yaGdXyaz9/E8YJEln9TWFfjTzuaTpkWjuYqISCgdQYiISCgVCBERCaUCISIioVQgREQklAqEiIiEUoEQEZFQKhAiIhLq/wOKWp46BI+1vwAAAABJRU5ErkJggg==\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["# A,B = generate_AB()\n", "\n", "plt.hist( t, density=True, alpha=.5, bins=30);\n", "plt.axvline(np.percentile(t,2.5), color=\"black\", label=\"%5 two sided interval\\n(2.5% - 97.5%) \")\n", "plt.axvline(np.percentile(t,97.5), color=\"black\")\n", "plt.grid(); plt.legend(); \n", "plt.xlabel(\"difference of means\"); plt.ylabel(\"probability\")\n", "plt.show()\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### STEP 5: How rare is our real world sample w.r.t. $H_0$?\n", "\n", "we measure how much probability mass is left outside the red dashed lines (our real world sample)\n", "\n", "- small $p_{value}$ means our real world sample is rare\n", "- large $p_{value}$ means our real world sample is quite common (expected)\n", "\n", "For a confidence interval $\\alpha=0.05$, if $p_{value}<\\alpha$ we reject $H_0$ and consider $A$ and $B$ actually come from different populations.\n", "\n", "\n", "observe we display $\\text{ttest}(A,B)$ and $\\text{ttest}(B,A)$ (red dashed lines) since we want to check if they are different in any direction"]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["group A  [  96.08 101.76  89.45  88.81 106.66 105.94  95.99 108.45  94.67 114.32 ]  mean 100.21\n", "group B  [ 129.23 127.05 115.34 105.17 117.75 121.96 103.69 120.42 118.77  89.56 ]  mean 114.89\n"]}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfXwU1fX48c8hBMJDKgqIGp4VUHkKEASrQFCeKiKKKIIi1CK1hZ+CVYtPaAWttlqRCgq2QKsIIqgEpSJQoqAoBI0oKIKAEqRfEUQSSQKB8/tjJusm2U0msJNd4nm/XvPa3Zk7MyebvXt35s6cK6qKMcYYU1yVaAdgjDEmNlkDYYwxJiRrIIwxxoRkDYQxxpiQrIEwxhgTUtVoBxAp9erV06ZNm0Zug1u2OI+tWnkq/uOPP1KrVq3I7d+jLW6crcLEGa24ShOLMUFsxhWLMUH4uMr6PPqpXO9VOev3iYjF/2FwTBs2bPhOVeuHLKiqlWLq1KmTRlSPHs7k0apVqyK7f4969OihPUqJM1pxlSYWY1KNzbhiMSbV8HGV9Xn0U7neq3LW7xMRi//D4JiADA3zvWqnmIwxxoRUaU4xRdx990U7AmOMX6x+e2INRDi9ekU7AmOMX6x+e2INRDiZmc5jcnJ04zDGdeTIEbKyssjLy4t2KJxyyil89tlnJeY/8MADACGX+S1cTCEdPuw8VqvmX0CucsXlo4SEBBo2bEh8fLzndayBCGfcOOcxPT2qYRhTKCsri8TERJo2bYqIRDWW7OxsEhMTS8yvUsXp1ozGVUzhYgqpAq9iKldcPlFV9u3bR1ZWFs2aNfO8nnVSG3OSyMvLo27dulFvHMzJR0SoW7duuY8+rYEw5iRijYM5Xsfz2bEGwhhjTEjWQBhjPNm7dy8XX3wxbdq04fXXXw/MHzhwIN98803g9d///neSkpJITk4mOTmZpUuXltjWzp07efHFFysk7lGjRrF58+YS8+e88gpjx44t17aaNm3Kd999V2L+ZZddxoEDB0pdd86cOUXeJ788+OCDPP744xHZlnVSh/PII9GOwETYk8u/KHV5Ul5+oMz43i0rIqSTyrx587jlllsYNGgQffv2ZejQoSxZsoQOHTpw1llnFSk7fvx47rjjjrDbKmwghg0b5nfY/OMf/yg5MykJTj0Vdu+OyD5CNYLFzZkzhzZt2pR4r0pTUFBA1arR+5r29QhCRPqJyBYR2SYiE0opd7WIqIikBM27211vi4j09TPOkH75S2cyxgAQHx/PoUOHyM/PJy4ujoKCAqZMmcJdd91V7m1NmDCB1atXk5yczJNPPkn//v3ZuHEjAB06dOChhx4CYOLEiTz33HOoKnfeeSdt2rShbdu2vPTSSyW2+eOPP9K/f3/at29PmzZtAmVSU1PJyMgAYPbs2bRs2ZILLrmEdzdsCKy7d+9err76ajp37kznzp159913Adi3bx99+vShdevWjBo1Cg0zAmfhkcXOnTs577zzuPnmm2ndujV9+vQhNzeXhQsXkpGRwfXXX09ycjK5ubls2LCBHj160KlTJ/r27cuePXsC8Y4bN46UlBQefvhhmjRpwrFjxwJ/Y6NGjThy5AjPPfccnTt3pn379lx99dUcOnSo3P+HsvjWQIhIHDAN+BVwPjBURM4PUS4RuA34IGje+cB1QGugHzDd3V7Fee89ZzLGADBs2DAWL15M7969+cMf/sD06dMZPnw4NWvWLFH26aefpl27dtx00018//33JZY/+uijdOvWjczMTMaPH0+3bt1YvXo1P/zwA1WrVg18Qa9evZru3bvzyiuvkJmZyccff8yKFSu48847A1+ohVasWMFZZ53Fxx9/zKeffkq/fv2KLN+zZw8PPPAA7777LmvefJPNn3wSWHbbbbcxfvx41q9fz6JFixg1ahQAf/rTn7j44ovZtGkTV111FV9//XWZ79PWrVsZM2YMmzZtok6dOixevJjBgweTkpLC3LlzyczMpGrVqvy///f/WLhwIRs2bOCmm27i3nvvDWzj8OHDZGRk8MADD5CcnMzbb78NwOuvv07fvn2Jj49n0KBBrF+/no8//pjzzjuPf/7zn2XGVl5+HkFcAGxT1e2qehiYDwwMUW4S8BgQfP3VQGC+quar6g5gm7u9inPPPc5kjAGcG77eeOMNMjIyaN++PUuWLGHw4MHcfPPNDB48mLVr1wIwdOhQvvzySzIzMznzzDP5wx/+UOa2u3XrxjvvvMO7775L//79ycnJ4dChQ+zYsYNWrVqxZs0ahg4dSlxcHA0aNKBHjx6sX7++yDbOP/98li9fzh//+EdWr17NKaecUmT5Bx98QGpqKvXr16fa3r0MufTSwLIVK1YwduxYkpOTueKKKzh48CA5OTm888473HDDDQD079+fU089tcy/pVmzZiS7N9h26tQpZKOyZcsWPv30U3r37k1ycjKTJ08mKysrsHzIkCFFnhceDc2fPz+w7NNPP6Vbt260bduWuXPnsmnTpjJjKy8/T24lAbuCXmcBXYILiEhHoJGqviEidxZb9/1i6yYV34GIjAZGAzRo0ID0CN7Ulux2OGV63GZOTk5E9+9VYcdYuH1HK67SRCumpLz8UpfHH8snKW8HAOnp/ncmehH8Xp1yyilkZ2dHNyDXY489xvjx45k9ezadOnXiyiuv5Prrr2fq1KnUqVMncLpj6NChXHvttSXiPnToEAUFBYH55557LuvWraNhw4b07NmTb775hr///e+0b9+e7OxsDh8+TF5eXqD8kSNHyM3NLbLd5s2b8/bbb/PWW29x991306NHDyZMmMDRo0f58ccfyc3N5ciRI2RnZ1OjoICjx45x+PBhsrOzOXr0KMuXLychISGwPVXl2LFj5OTkBPajquTk5FC9evUif0/h/JycHOLj4wPlCwoKAvssjCM7O5ucnBzOPfdcVq5cWWQ7heUKnwP07NmTu+++m6+++oqMjAw6d+5MdnY2I0aM4MUXXww0EKtXryY7O5v8/PwiMQTLy8sjPT3dcx2MWu+HiFQB/gaMPN5tqOpMYCZASkqKpqamRiQ2AOrUAZzzgV6kp6d7LhtJdcqIM1pxlSZaMZXdSb2D3QnOXabXpsZGJ3Xwe/XZZ59F/Y5ccE6h7Nmzh8suu4ynnnqKWrVq8Ytf/IIjR45QtWpVvv32W1q3bg3A8uXLadeuXYm4GzRoQG5ubpH5TZo0YfHixUyePJkff/yRO+64gzvuuIPExEQuvfRSZsyYwW9/+1v279/P2rVrmTJlSpH1v/jiCxo3bszNN9/MmWeeyT/+8Q8SExOJi4ujVq1apKamMmHCBA4fPkyCKq8sX077rl1JTEykb9++zJkzhzvvdH6nZmZmkpycTGpqKmlpadx333385z//4cCBA9SuXbvE3yMi1K5dG3DuJi9cXr16dXJyckhMTKROnTocO3aMxMREOnbsyP79+/n000+58MILOXLkCF988QWtW7cOxFu4jcTERC644ALuu+8+BgwYEKjzOTk5nHPOOSQkJLBo0SKSkpJITEykevXqVK9ePeRnJSEhgQ4dOniug342ELuBRkGvG7rzCiUCbYB09waOM4A0EbnCw7rGmCi59957uf/++wHnCOHKK6/k0UcfDXQsP/744/z+979HRGjatCkzZswosY127doRFxdH+/btGTlyZKAfYuXKldSoUYNu3bqRlZVFt27dALjqqqtYu3Yt7du3R0T4y1/+whlnnFFkm5s2bWLw4MFUqVKF+Ph4nnnmmSLLzzzzTB588EEuvPBC6iQkkHzuuYFlU6dOZcyYMbRr146CggK6d+/Os88+ywMPPMDQoUNp3bo1v/zlL2ncuPFxv28jR47klltuoUaNGqxdu5aFCxdy66238sMPP1BQUMC4ceMCDWtxQ4YM4Zprrinyq3/SpEl06dKF+vXr06VLF1+OLiVcr/wJb1ikKvAFcCnOl/t6YJiqhjxRJiLpwB2qmiEirYEXcfodzgJWAi1U9Wi4/aWkpGjhlQoRUdi6ejwVEq1fxYX7DHe4aEcQPynPEUSsXOZa/AjivPPOi25ArnD5haI5opzlYipb4Wco+HMlIhtUNSVUed+OIFS1QETGAsuAOGCWqm4SkYdwRjBKK2XdTSKyANgMFABjSmscfDFlSoXuzhhTgRo1KruM8bcPQlWXAkuLzZsYpmxqsdcPAw/7FlxZLM23iWHjxo0jszAlfYQkJycz5efywyjEpbmmJEu1Ec6KFc5kjKl8Dh50JlMqS7URzuTJzqONPGViULR+6efm5tKvXz8WL15MZmYmv/vd7zh48CBxcXHce++9gev/gxVeHZSU5FypPnbsWEaNGsWWLVsYNmwYR44cYcaMGVx44YUUFBTQr18/0tLSQt6AV5pFixbxt7/9jaNHj3L55Zfz2GOPAU7aj1WrVgHO5bXffvstBz5w78v9xS8C66emprJnzx5q1KgBwFtvvcXpp5/O3//+d2bMmEHjxo157bXXqFatGmvWrGHRokU8+eSTgHMn9vDhw3nzzTfL94bGOGsgjDGezZo1i0GDBhEXF0fNmjX597//TYsWLfjmm2/o1KkTaWlp/CLoS7fQkCFDePrpp4vMmzFjBk899RRNmzbltttuY9GiRTzzzDPccMMN5W4c9u3bx/3338+HH35I/fr1GTFiBCtXruTSSy8NfImDk0jwo48+CruduXPnkpKSUmLexo0beeSRR1i2bBmXX345kyZNYt68eYEy9evX58wzz+Tdd9/loosuKlfsscxOMRljPJs7dy4DBzoJEVq2bEmLFi0AOOusszj99NPZv3+/520V5nY6dOgQ8fHxHDhwgCVLlnDjjTeWO67t27dz9tlnU79+fQB69erFokWLSpSbN28eQ4cOLde2VZUjR44E4nzhhRf41a9+xWmnnVak3JVXXsncuXPLHXssswbCGOPJ4cOH2b59O02bNi2xbN26dRw+fDjsfQKLFi2iXbt2DB48mF27nAQLY8aM4ZFHHmHEiBHcc889TJo0iXvuuScwbGl5nHPOOWzdupWdO3dSUFDAa6+9FthPoa+++oodO3ZwySWXhN3Or3/9a5KTk5k0aVIgMd/YsWPp2rUrX3/9NRdddBGzZ89mzJgxJdZNSUlh9erV5Y49llkDYYzx5LvvvgvcxRtsz549DB8+nNmzZ4f8ch8wYAA7d+5k48aN9O7dmxEjRgDQuHFj0tPTWbt2LTVr1iQrK4vzzjuP4cOHM2TIEL74ovT7VoKdeuqpPPnkkwwZMoRu3brRtGlT4uKK5vecP38+gwcPLjG/0Ny5c/nkk09YvXo1q1ev5vnnnwdg+PDhfPTRR7zwwgs8+eST3HrrrfznP/9h8ODBjB8/PpBp9fTTT6+Q8R4qkjUQ4cyY4UzGGABq1KhRYkzjgwcP0r9/fx5++GG6du0acr26desGcheNGjWKDUFptgvde++9TJ48malTpzJq1Cj+8pe/8Kc//alImV27dgUGIXr22WdLbONXv/oVH3zwAWvXrqVVq1a0bFn0Zsf58+f/dHqpSRNnClLYiZ6YmMiwYcNYt25dkeXffPMN69at48orr+SJJ57gpZdeok6dOoF8Snl5eYEO7srCOqnDicKdoMbEslNPPZWjR48GGonDhw9z1VVXceONNzJ48OCw6+3Zs4czzzwTgLS0tBJ3g7/99tucddZZtGjRgkOHDlGlShWqVKlSYnyDRo0alXrvx969e0lMTOT7779n+vTpLFiwILDs888/5/vvv+fCCy90ZgQl5QMnqd6BAweoV68eR44c4fXXX6dXsSsY77///kA6kdzcXESkSJxffPEFbdq0CRvfycgaiHCWLHEeBwyIbhzGxJA+ffqwZs0aunTpwoIFC3jnnXfYt28fc+bMAZzhLs877zwmTpxISkoKV1xxBVOnTiUtLY2qVaty2mmnBcqC0wE8efLkQDrr0aNHc/3111NQUFAil1JZ7rrrrsDQohMnTixyBDF//nyuu+463Lxv4GZBTk5NJTMzk/z8fPr27cuRI0c4evQovXr14uabbw6sX3jlU8eOHQFnbIy2bdvSqFGjwIBJq1aton///uWKOdb5loupolkupvSQyy0X008sF9OJ+/DDD3nyySeZPn265WIqpnv37ixevLjEmBEncy4m64MwxnjWsWNHevbsGRizwDj27t3L7bff7mlAoZOJNRDGmHK56aabwl4J9HNVv359rrzyymiHEXHWQBhjjAnJGghjjDEh2VVM4bg3yRhjKqFmzaIdwUnBGohwbEARYyqvatWiHcFJwddTTCLST0S2iMg2EZkQYvktIvKJiGSKyBoROd+d31REct35mSJS8rZJv730kjMZYyKmadOmfPfdd75tf+TIkSxcuLDE/PT0dC6//PKfZuzf70wVrHbt2uVeJzMzExEpNZV4XFwcycnJtG/fno4dO/Lee++dSJgBvh1BiEgcMA3oDWQB60UkTVU3BxV7UVWfdctfAfwN6Ocu+1JVozesW+FNOkOGRC0E401Z9zeYyFNVjh07dlyJ9fxSrktv9+51HotlZI1F8+bN4+KLL2bevHn069cvZJkaNWoE7jJftmwZd999N2+//fYJ79vP/+4FwDZV3a6qh4H5wMDgAqoaPKRTLaBy3LVnTEVITS05TZ/uLDt0KPTywruYv/uu5LIy7Ny5k1atWnHjjTfSpUsXdu3axV//+lc6d+5Mu3bteOCBBwJlr7zySjp16kTr1q2ZOXNmqdt9+eWXuf322wF46qmnaN68OeCk8C4cW2HlypV06NCBtm3bctNNN5Gfnw84RyR//OMf6dixI6+++mqR7b755puce+65dOzYkVdeeSXkvjdt2sQFF1xAcnIy7dq1Y+vWraXGX7t2be68805at25Nr169WLduHampqTRv3py0tDTAGSBp4MCBpKam0qJFC/785z+H3He49y6YqvLyyy8zZ84cli9fXiIXVigHDx6M2P0YfvZBJAHB+XazgC7FC4nIGOB2oBoQnIe3mYh8BBwE7lPVEnl0RWQ0MBqgQYMGYe8mPh7J7q34mR63mZOTE9H9e3XAjTPcvqMVV2kiHVNSXn5EthN/LJ+kvB0ApKfHRlbO4PfqlFNOITs7O7CsRohfzAV5eRzJzoZDh0IuP5KXR0F2NpKTQ0Kx5blB2w4Xy9atW5k+fTpTp04lPT2dzZs3s3LlSlSVIUOG0Lx5c1JSUnjqqac47bTTyM3NJTU1lT59+lC3bl1UlZycnEDyPoAOHTrw5z//mezsbFatWkWdOnXYsmULq1atomvXruzdu5cRI0aQlpZGixYtGD16NE8++SRjxoxBValduzZvv/02R48eZcWKFeTm5rJ3715GjRrFkiVLOPvssxk5ciQFBQWB969GQQEAU6dOZfTo0QwZMoTDhw9z9OhRsrOzw8b/448/0rVrVyZOnMiwYcOYMGECr7zyCp9//jm33HILPXv2JC8vjw8++ID333+fmjVr0qNHD/r27RtI05Gdnc3KlStLvHdvvvlmicGG3n//fRo3bszpp5/ORRddxMKFCwPjcRT53+Xm0q5dO/Ly8vi///s/lixZUuSzUigvL4/09HTPdTDqndSqOg2YJiLDgPuAEcAeoLGq7hORTsBrItK62BEHqjoTmAlOqo2Ipm9w0xp73Wa00kfUKSPOn0OqjUidYgpOtXFtamym2iiSsiHE2ANVgQSAxMSwywmzvKxkELVr16ZJkyZceumlZGdns2bNGlatWkX37t0BpwHJysqia9euzJ49O/CLfvfu3fzvf/+jadOmiAi1a9cu8nckJiaSm5sL/JQ6/MMPPyQjI4NBgwbxzTff0Lx588AX7KhRo5g2bRoTJkxARLjxxhtJTEwkOzub+Ph4atSowe7du2nevDkdOnQAnL6JmTNn/rTfqs470aNHDx5++GH27dvHoEGDAgMgPfHEEyHjr1atGoMGDUJE6NChA9WrV+e0004LjBeRmJhIQkICffr0CYybccUVV/DRRx/Ro0ePwN8b6r3bvXt3iZQcixcv5vrrrycxMZHhw4fz73//mxtuuKHE/6ZGjRps3LgRgLVr1zJq1Cg+/fTTn3JPuRISEujQoYPnOuhnA7EbCL4UqKE7L5z5wDMAqpoP5LvPN4jIl0BLIILJlowx5VWrVq3Ac1Xl7rvv5re//W1g3pYtW/jggw9YsWJFYJyH1NTUMk+N/PKXv2T27Nm0atWKbt26MWvWLNauXcsTTzzBzp07PcdUXsOGDaNLly688cYbXHbZZcyYMYMqVaqEjT8+Pj7wpVulSpXAkVCVKlUocI9KgBJfzMVfh3rvijt69CiLFi1i8eLFPPzww6gq+/btIzs7m1tvvZWPPvqIs846i6VLlxZZ78ILL+S7775j7969nH766cf93oC/DcR6oIWINMNpGK4DhgUXEJEWqrrVfdkf2OrOrw/sV9WjItIcaAFs9zHWkkJcCWF+PrwelcRKUr9o6Nu3L/fffz/XX389tWvXZvfu3ezbt4+cnBxOPfVUatasyeeff877779f5ra6devGxIkTmThxIh06dGDVqlXUqFGDU045hVatWrFz5062bdvGOeecw/PPPx/4NR7Oueeey86dO/nyyy85++yzi4wfDUBQP0fz5s259dZb+frrr9m4cSPNmjUrd/zFLV++nP3791OjRg3eeOONIhlsIfR7Fx8fX+QLfeXKlbRr145ly5YF5o0YMYJXX32V2bNnh933559/ztGjR6lbt2654y7OtwZCVQtEZCywDIgDZqnqJhF5CMhQ1TRgrIj0Ao4A3+OcXgLoDjwkIkeAY8Atqlqx16TVq1ehuzPmZNOnTx8+++yzwBgLtWvX5qGHHqJbt24sWbKE8847j1atWoUdSChYt27d2LVrF927dycuLo5GjRpx7rnnAs5pkdmzZ3PNNddQUFBA586dueWWW0rdXkJCAjNnzqR///7UrFmTbt26FT0nHx8PwIIFC3j++eeJj4/njDPO4J577qFWrVo8++yz5Yq/uAsuuICrr76arKwsrrnmGlJSiiZLDfXevfDCC0UaiHnz5nHVVVcVWe/qq6/mmWeeKTFud25uLsnJzkWfqsq//vWviOTLsnTf4RS2+CNHeipu6b69Oxn6ILzy+wgi1tJ9FwqXwvqkSfddeC+GDz8E58yZQ0ZGBk8//XT54/KZpfuOlDlzfmokjDGVy759zmRKFfWrmIwxpjIZOXIkIz2eeYh1dgRhjDEmJGsgjDHGhGQNhDHGmJCsDyKcYjefGGMqkXPOiXYEJwU7gginZk1nMsZETMyk+46Lc6YKVt50302bNqVt27YkJyfTtm1bFi9eHLZsQUEB9evXZ8KEEiMrHDdrIMKZPv2nzJjGmCIK033HknKl+/72W2c6CaxatYrMzEwWLlzIrbfeGrbc8uXLadmyJS+//DKRur/NGohwFixwJmNK8eTyLzxNvrB038ef7vv77+H772M+3XewstJ4z5s3j9tuu43GjRuzdu3aMrfnhTUQxhjPtm7dyu9//3vWrVvHli1b2Lp1K+vWrSMzM5MNGzawfv16AGbNmsWGDRvIyMhg6tSp7CvlprRu3bqx2s0su3r1aurWrcvu3btZvXo13bt3Jy8vj5EjR/LSSy/xySefUFBQwDOFA3oBdevW5cMPP2Tw4MGBeXl5edx8880sWbKEDRs28L///S/kvp999lluu+02MjMzycjIoGHDhqXG/+OPP3LJJZewadMmEhMTue+++1i+fDmvvvoqEydODGx33bp1LFq0iI0bN/Laa69RPMvDW2+9VeK9e+edd0LG2LNnT9q0aUOPHj2YPHlyyDJ5eXmsWLGCAQMGMHTo0JK5p46TdVIbc7IqLZ9/zZqlL69Xr/TlYTRp0oSuXbuSnZ3NW2+9xVtvvRVIqZ2Tk8NXX31F586dmTp1auAX/a5du9i6dWvY5HFnnHEGOTk5ZGdns2vXLoYNG8Y777zD6tWrGTRoEFu2bKFZs2a0bOmkNRkxYgTTpk1j3LhxAAwJMerj559/TrNmzQLpu2+44YaQRzIXXnghDz/8MFlZWUXSfYeLv1q1aoFR3dq2bUv16tWJj4+nbdu2RbLO9u7dO/D3DhgwgDVr1hTJxxTqvdu6dWsg/XewVatWUa9ePb788ksuvfRSUlNTS/RlvP766/Ts2ZMaNWpw9dVXM2nSJKZMmXLC+ZisgTDGeGbpvisu3XdxZ599Ng0aNGDz5s2sX7+e5557DoClS5cyb9481qxZExiDYt++ffz3v/+ld+/e5XtDirFTTMaY49K3b19mzZpFTk4OwAmn+3788cfp3r17IN139erVS6T7Bsqd7hsIe8olON33wIED2bhxIz/88EPE0n3n5ubyxhtvlBgpLtR7920ZnebffvstO3bsoEmTJowZM4bMzEwyMzOpXbs2q1ev5uuvv2bnzp3s3LmTadOmReQ0kx1BhBNjw3QaE2tO6nTfbrbZBY8+GrPpvgv17NmTuLg4jhw5wqOPPkqDBg2KLH/11Ve55JJLigzjOnDgQO666y7y8/OLzC8vS/cdIZbu27vKlO7bq+NNC27pvr2LlbTalu775+Dxx53JGFP5/O9/zmRK5WsDISL9RGSLiGwTkRK394nILSLyiYhkisgaETk/aNnd7npbRKSvn3GG9PrrzmSMqXx++MGZfDBy5MjA0cPJzrc+CBGJA6YBvYEsYL2IpKnq5qBiL6rqs275K4C/Af3chuI6oDVwFrBCRFqqajlulTQnO99uMDuJqWqJK2KM8eJ4uhP8PIK4ANimqttV9TAwHxgYXEBVDwa9rAUU/gUDgfmqmq+qO4Bt7vaM+dlKSEhg3759EUujYH4+VJV9+/aRkJBQrvX8vIopCdgV9DoL6FK8kIiMAW4HqgGXBK0bfG1ZljvPmJ+thg0bkpWVxd69e6MdCnl5eSG/bArvWI5GnqZwMYVU2P9QAXGWKy4fJSQkBO4U9yrql7mq6jRgmogMA+4DRnhdV0RGA6MBGjRoEPZKnuPRNjcXgE88bjMnJyei+/fqwIEDQPirmKIVV2m8xpSUl+9/MEHij+WTlLfDl22np39zXOvF4v8PnKeW1n4AABvMSURBVLhCZSYdP348AFOmTKnokMLGFEpbN+PpJ4895mdIQPni8ttXX30FeP9c+dlA7AYaBb1u6M4LZz5QmGDF07qqOhOYCc5lrhG9nPODDwDwusVoXU5ap04dgLD7Ppkvc63oPgg/L3O9NvXEL3ONJeHiKuvz6KdyvVflrN8nIhb/h15j8rMPYj3QQkSaiUg1nE7ntOACItIi6GV/YKv7PA24TkSqi0gzoAWwzsdYjTHGFOPbEYSqFojIWGAZEAfMUtVNIvIQkKGqacBYEekFHAG+xz295JZbAGwGCoAxFX4F06RJzuP991fobo0xFcDqtye+9kGo6lJgabF5E4Oe31bKug8DD/sXXRlWrnQe7QNkTOVj9dsTu5PaGGNMSNZAGGOMCckaCGOMMSFF/T6ImBVm9CtjTCVg9dsTayDCWbQo2hEYY/xi9dsTayCMqQBeb/o73nEjjPGD9UGEc/fdzmSMqXysfntiRxDhrF0b7QiMMX6x+u2JHUEYY4wJyRoIY4wxIVkDYYwxJiTrgwinnANrGGNOIla/PbEGIpwXXoh2BMYYv1j99sROMRljjAnJUwMhIq+ISH8R+fk0KOPGOZMxpvKx+u2J11NM04FfA1NF5GVgtqpu8S+sGJCZGe0IKq1vs/MrfDhRY4qw+u2JpyMCVV2hqtcDHYGdwAoReU9Efi0i8X4GaIwxJjo8nzISkbrASGAU8BHwFE6DsbyUdfqJyBYR2SYiE0Isv11ENovIRhFZKSJNgpYdFZFMd0orvq4xxhh/eTrFJCKvAq2A54EBqrrHXfSSiGSEWScOmAb0BrKA9SKSpqqbg4p9BKSo6iER+R3wF2CIuyxXVZPL/RcZY4yJCK99EM+540sHiEh1Vc1X1ZQw61wAbFPV7W75+cBAINBAqOqqoPLvAzd4jtxvLS2rpjGVltVvT0RVyy4k8qGqdixrXrHlg4F+qjrKfT0c6KKqY8OUfxr4n6pOdl8XAJlAAfCoqr4WYp3RwGiABg0adJo/f36Zf4tfcnJyqF27doXvd5x7JcaUKVNCLo9WXKX5/oeDHKlSPdphlBB/LD/qcZ2eWHT/sfj/g/BxlfV59NPJ9l5FU3BMPXv23BDuh36pRxAicgaQBNQQkQ6AuIt+AdSMVLAicgOQAvQImt1EVXeLSHPgvyLyiap+Gbyeqs4EZgKkpKRoampqpEIqt/T0dKKx/zp16gCE3Xe04irNgiXL2J3QLNphlJCUtyPqcV2bWvSXbSz+/yB8XGV9Hv10sr1X0eQ1prJOMfXF6ZhuCPwtaH42cE8Z6+4GGgW9bujOK0JEegH3Aj1UNb9wvqrudh+3i0g60AH4svj6vhk92nmcObPCdmmMqSBWvz0ptYFQ1X8B/xKRq1W1vGP0rQdaiEgznIbhOmBYcAH3qGQGzqmob4PmnwocUtV8EakHXITTgV1xvrDr9I2ptKx+e1LWKaYbVPUFoKmI3F58uar+LcRqhcsKRGQssAyIA2ap6iYReQjIUNU04K9AbeBlEQH4WlWvAM4DZojIMZxLcR8tdvWTMcYYn5V1iqmW+3hcPSzulU9Li82bGPS8V5j13gPaHs8+jTHGREZZp5hmuI9/qphwjDHGxIqyTjFNLW25qt4a2XBiSLLdo2dMpWX125OyTjFtqJAoYlEUruM2xlQQq9+eeLmKyRhTQYpnuU3KC535dnxvuxPY+K+sU0xTVHWciCwBStxy7V5xVDnd4Gb9sJGnjKl8rH57UtYppufdx8f9DiTmZGVFOwJjjF+sfntS1immDe7j2yJSDTgX50hii6oeroD4jDHGRInXdN/9gWdxUl0I0ExEfquq//EzOGOMMdHjNd33E0BPVd0GICJnA28A1kAYY0wl5bWByC5sHFzbcRL2VV4XXhjtCIwxfrH67UlZVzENcp9miMhSYAFOH8Q1OMn4Kq8//znaERhj/GL125OyjiAGBD3/P34ar2EvUMOXiIwxxsSEsq5i+nVFBRJzrr7aeVxU3iznxpiYZ/XbE69XMSUAvwFaAwmF81X1Jp/iir59+6IdgTHGL1a/PanisdzzwBk4I8y9jTM6XOXupDbGmJ85rw3EOap6P/Cjm5+pP9DFv7CMMcZEm9cG4oj7eEBE2gCnAKeXtZKI9BORLSKyTUQmhFh+u4hsFpGNIrJSRJoELRshIlvdaYTHOI0xxkSI1/sgZrrjRN8PpOGMMHd/aSuISBwwDegNZAHrRSSt2NChHwEpqnpIRH6HM+70EBE5DXgASMG5rHaDu+735fjbTsyll1bYrowxFczqtyeeGghV/Yf79G2gucdtXwBsU9XtACIyHxgIBBoIVV0VVP59wE2xSF9guarud9ddDvQD5nnc94m7v9T2zxhzMrP67YnXq5jqAg8CF+H8ol8NTFLV0i4FSAJ2Bb3OovR+i9/wU+qOUOsmhYhrNDAaoEGDBqSnp5f2Z/gqJycnKvs/cOAAQNh9Ryuu0sQfyycpb0e0wyghFuMKF1N6+jdRiOYn4T5XZX0e/RSLn3WIzbi8xuT1FNN84B3AvXiY64GXgF7HE1xxInIDzumkHmWVDaaqM4GZACkpKZqamhqJcBy/+pXz+B9v6abS09OJ6P49qlOnDkDYfUcrrtIsWLKM3QnNoh1GCUl5O2IurnAxXZsa3QGDwn2uyvo8+qlcn/Vy1u8TEYt10GtMXhuIM1V1UtDrySIypIx1dgONgl43dOcVISK9gHuBHqqaH7RuarF10z3GGhm5uRW6u8og1MhnoZQ4FDSmoln99sTrVUxvich1IlLFna4FlpWxznqghYg0c8eSuA6ngztARDoAM4ArVPXboEXLgD4icqrbOd7Hw/6MMcZEUFnJ+rJx+hwEGAcUjs9XBcgB7gi3rqoWiMhYnC/2OGCWqm4SkYeADFVNA/6Kc0XUyyIC8LWqXqGq+0VkEj8lBHyosMPaGGNMxSgrF1PiiWxcVZcCS4vNmxj0PGwfhqrOAmadyP6NMcYcP699EIjIFUB392W6qr7uT0gx4vLLox2BMcYvVr898XqZ66NAZ2CuO+s2EblIVe/2LbJouyPs2TNjzMnO6rcnXo8gLgOSVfUYgIj8C+cu6MrbQBhjzM+c16uYAOoEPT8l0oHEnNRUZzLGVD5Wvz3xegTxCPCRiKzCuaKpO1Ai+Z4xxpjKo8wGQkSqAMeArjj9EAB/VNX/+RmYMcaY6CqzgVDVYyJyl6ouoNiNbsYYYyovr30QK0TkDhFpJCKnFU6+RmaMMSaqvPZBDMG5o/r3xeZ7Tf198rn22mhHYIzxi9VvT7w2EOfjNA4X81O672f9Ciom/L54W2iMqTSsfnvitYH4F3AQmOq+HubOq7zN8KFDzmPNmtGNwxgTeVa/PfHaQLRR1fODXq8Skc1hS1cGl13mPMbYQB/GgPfU6uN7R3fciJhl9dsTr53UH4pI18IXItIFyPAnJGOMMbHA6xFEJ+A9Efnafd0Y2CIinwCqqu18ic4YY0zUeG0g+vkahYlpXk9nGGMqF08NhKp+5XcgxhhjYovn8SCOh4j0A57CGVHuH6r6aLHl3YEpQDvgOlVdGLTsKPCJ+/JrVb3Cz1hLGDmyQndnjB+sMzsMq9+e+NZAiEgcMA3oDWQB60UkTVWDr376GhhJ6KFLc1U12a/4ymQfIGMqL6vfnvh5BHEBsE1VtwOIyHxgIBBoIFR1p7vsmI9xHJ/vvnMe69WLbhzGmMiz+u1JecaDKK8kYFfQ6yx3nlcJIpIhIu+LyJWRDc2DwYOdyRhT+Vj99sTXPogT1ERVd4tIc+C/IvKJqn4ZXEBERgOjARo0aEB6BG96ST5wAIBMj9vMycmJ6P69OuDGGW7fkYgrKS//hNYvLv5YPkl5OyK6zUiIxbgqKqb09G/KVT7c56qsz6OfyvNZL2/9PhHR+m4ojdeY/GwgdgONgl43dOd5oqq73cftIpIOdAC+LFZmJjATICUlRVMjOUJUHWcAPa/bTE9P91w2kuqUEWck4or0Za5JeTvYndAsotuMhFiMq6Jiuja1fJ3U4T5XZX0e/VSuz3oFxhmt74bSeI3Jz1NM64EWItJMRKoB1+FxPAkROVVEqrvP6wEXEdR3YYwxxn++NRCqWgCMBZYBnwELVHWTiDwkIlcAiEhnEckCrgFmiMgmd/XzgAwR+RhYBTxa7OonY4wxPvO1D0JVlwJLi82bGPR8Pc6pp+LrvQe09TO2Mv3ud1HdvTHGR1a/PYnlTuroGjIk2hEYY/xi9dsTP/sgTm67djmTMabysfrtiR1BhDN8uPMYY5enGWMiwOq3J3YEYYwxJiRrIIwxxoRkDYQxxpiQrIEwxhgTknVSh/OHP0Q7AmOMX6x+e2INRDgDBkQ7AmOMX6x+e2KnmMLZssWZjDGVj9VvT+wIIpzf/tZ5tOukjal8rH57YkcQxhhjQrIGwhhjTEjWQBhjjAnJGghjjDEhWSd1OPfdF+0IjDF+sfrtiTUQ4fTqFe0IjKkw5Rl3fHzv8o1fHZOsfnvi6ykmEeknIltEZJuITAixvLuIfCgiBSIyuNiyESKy1Z1G+BlnSJmZzmSMqXysfnvi2xGEiMQB04DeQBawXkTSio0t/TUwErij2LqnAQ8AKYACG9x1v/cr3hLGjXMe7TppYyofq9+e+HkEcQGwTVW3q+phYD4wMLiAqu5U1Y3AsWLr9gWWq+p+t1FYDvTzMVZjjDHF+NkHkQQEj+mXBXQ5gXWTihcSkdHAaIAGDRqQHsFfA8kHDgCQ6XGbOTk5Ed2/VwfcOMPtOxJxJeXln9D6xcUfyycpb0dEtxkJsRhXLMaUnv5N2M9VWZ9HP5Xns17e+n0iovXdUBqvMZ3UndSqOhOYCZCSkqKpqamR23idOgB43WZ6errnspFUp4w4S4vLc8dkwnEEVoqkvB3sTmgW2Y1GQCzGFYsxXZvaMuznqqzPo5/KVQcrMM5ofTeUxmtMfp5i2g00Cnrd0J3n97rGGGMiwM8jiPVACxFphvPlfh0wzOO6y4BHRORU93Uf4O7Ih1iKRx6p0N0ZYyqQ1W9PfGsgVLVARMbifNnHAbNUdZOIPARkqGqaiHQGXgVOBQaIyJ9UtbWq7heRSTiNDMBDqrrfr1hD+uUvK3R3xpgKZPXbE1/7IFR1KbC02LyJQc/X45w+CrXuLGCWn/GV6r33nEf7IBlT+Vj99uSk7qT21T33OI8xdvWBMSYCrH57Ysn6jDHGhGQNhDHGmJCsgTDGGBOSNRDGGGNCsk7qcKZMiXYExhi/WP32xBqIcJKTox2BMcYvVr89sQYinBUrnEcbWMSYIp5c/gVJefkhc3llfZ8bKBPTAwtZ/fbEGohwJk92Hu0DZEzlY/XbE+ukNsYYE5I1EMYYY0KyBsIYY0xI1kAYY4wJyTqpw5kxI9oRGGP8YvXbE2sgwmnVKtoRGGP8YvXbEzvFFM6SJc5kjKl8rH574msDISL9RGSLiGwTkQkhllcXkZfc5R+ISFN3flMRyRWRTHd61s84Q3riCWcyxlQ+Vr898e0Uk4jEAdOA3kAWsF5E0lR1c1Cx3wDfq+o5InId8BgwxF32para/fDGGBMlfh5BXABsU9XtqnoYmA8MLFZmIPAv9/lC4FIRER9jMsYY45GfndRJwK6g11lAl3BlVLVARH4A6rrLmonIR8BB4D5VXV18ByIyGhgN0KBBA9IjOHxg8oEDAGR63GZOTk5E9+/VATfOcPsuLa6kvHyfoipd/LF8kvJ2RGXfpYnFuGIxJggfV/VjTi6mpLwdpKd/U6ExlacOlrd+n4hofTeUxmtMsXoV0x6gsaruE5FOwGsi0lpVDwYXUtWZwEyAlJQUTU1NjVwEdeoA4HWb6enpnstGUp0y4iwtrlDJ1ipCUt4Odic0i8q+SxOLccViTBA+rvwqNQDYndCMa1MrNllfuepgOev3iYjWd0NpvMbkZwOxG2gU9LqhOy9UmSwRqQqcAuxTVQXyAVR1g4h8CbQEMnyMt6jnn6+wXRljKpjVb0/8bCDWAy1EpBlOQ3AdMKxYmTRgBLAWGAz8V1VVROoD+1X1qIg0B1oA232MtaRGjcouU4HC/doPTq8cSri0zMb8rMVY/Y5VvjUQbp/CWGAZEAfMUtVNIvIQkKGqacA/gedFZBuwH6cRAegOPCQiR4BjwC2qut+vWEN66SXncciQ0ssZY04+Vr898bUPQlWXAkuLzZsY9DwPuCbEeouARX7GVqZnnnEe7QNkTOVj9dsTu5PaGGNMSNZAGGOMCckaCGOMMSFZA2GMMSakWL1RLvoWLqyQ3dglqMZEQQXV75OdNRDh1KsX7QiMOal5/fEzvnfF3nENWP32yE4xhTNnjjMZYyofq9+eWAMRjn2AjKm8rH57Yg2EMcaYkKwPwhgTVTHdV/EzZ0cQxhhjQrIGwhhjTEh2iimcpUvLLlMKu7/BmBh2gvX758IaiHBq1ox2BMYYv1j99sQaiHCmT3cef//7IrNtYB5joqOs+lVYBz11Zoep36Yo64MIZ8ECZzLGVD5Wvz2xBsIYY0xIvp5iEpF+wFM4Q47+Q1UfLba8OvBvoBOwDxiiqjvdZXcDvwGOAreq6jI/Yy1++DrYHet5oZ02Msb8TPnWQIhIHDAN6A1kAetFJE1VNwcV+w3wvaqeIyLXAY8BQ0TkfJzxqVsDZwErRKSlqh71K15jTOXgpS+wPD8Af8436Pl5BHEBsE1VtwOIyHxgIBDcQAwEHnSfLwSeFhFx589X1Xxgh4hsc7e31sd4jTGmhBO9+KQiLmDxqxETVfVnwyKDgX6qOsp9PRzooqpjg8p86pbJcl9/CXTBaTTeV9UX3Pn/BP6jqguL7WM0MNp92QrY4ssf40094Lso7j+cWIwrFmOC2IwrFmOC2IwrFmOC2IwrOKYmqlo/VKGT+jJXVZ0JzIx2HAAikqGqKdGOo7hYjCsWY4LYjCsWY4LYjCsWY4LYjMtrTH5exbQbaBT0uqE7L2QZEakKnILTWe1lXWOMMT7ys4FYD7QQkWYiUg2n0zmtWJk0YIT7fDDwX3XOeaUB14lIdRFpBrQA1vkYqzHGmGJ8O8WkqgUiMhZYhnOZ6yxV3SQiDwEZqpoG/BN43u2E3o/TiOCWW4DToV0AjDkJrmCKiVNdIcRiXLEYE8RmXLEYE8RmXLEYE8RmXJ5i8q2T2hhjzMnN7qQ2xhgTkjUQxhhjQrIGIoJEZJKIbBSRTBF5S0TOinZMACLyVxH53I3tVRGpEwMxXSMim0TkmIhE9RJAEeknIltEZJuITIhmLIVEZJaIfOveKxQTRKSRiKwSkc3u/+62aMcEICIJIrJORD524/pTtGMqJCJxIvKRiLwe7VgKichOEfnE/Z7KKK2sNRCR9VdVbaeqycDrwMRoB+RaDrRR1XbAF8DdUY4H4FNgEPBONIMISgnzK+B8YKib6iXa5gD9oh1EMQXAH1T1fKArMCZG3qt84BJVbQ8kA/1EpGuUYyp0G/BZtIMIoaeqJpd1L4Q1EBGkqgeDXtYCYuIKAFV9S1UL3Jfv49xXElWq+pmqRvPO90KBlDCqehgoTAkTVar6Ds6VfTFDVfeo6ofu82ycL76k6EYF6shxX8a7U9Trnog0BPoD/4h2LMfLGogIE5GHRWQXcD2xcwQR7CbgP9EOIoYkAbuCXmcRA196sU5EmgIdgA+iG4nDPZWTCXwLLFfVWIhrCnAXcCzagRSjwFsissFNVxSWNRDlJCIrROTTENNAAFW9V1UbAXOBsaVvreLicsvci3OaYG6sxGROPiJSG1gEjCt21Bw1qnrUPbXbELhARNpEMx4RuRz4VlU3RDOOMC5W1Y44p1XHiEj3cAVP6lxM0aCqvTwWnQssBR7wMZyAsuISkZHA5cClWkE3v5TjvYomS+tSDiISj9M4zFXVV6IdT3GqekBEVuH030Szg/8i4AoRuQxIAH4hIi+o6g1RjAkAVd3tPn4rIq/inGYN2RdoRxARJCItgl4OBD6PVizB3IGb7gKuUNVD0Y4nxnhJCWMANxX/P4HPVPVv0Y6nkIjUL7wyT0Rq4IxBE9W6p6p3q2pDVW2K85n6byw0DiJSS0QSC58DfSilIbUGIrIedU+hbMR542PiMkDgaSARWO5e2vZstAMSkatEJAu4EHhDRHwdMTAct/O+MCXMZ8ACVd0UjViCicg8nPFPWolIloj8Jtox4fwqHg5c4n6OMt1fyNF2JrDKrXfrcfogYuay0hjTAFgjIh/j5Ld7Q1XfDFfYUm0YY4wJyY4gjDHGhGQNhDHGmJCsgTDGGBOSNRDGGGNCsgbCGGNMSNZAmJOKiDwoIne4zx8SkV7u825uJs9MEanhZrDdJCJ/jW7E5Sci89zMu+OjHYv5ebM7qc1JS1WDc11dD/xZVV8AcHPMnOZ1qFoRqRqU0DBqROQMoLOqnhPtWIyxIwgT80TkXhH5QkTWAK2C5s8RkcEiMgq4FpgkInNFJA2oDWwQkSHunbaLRGS9O13krv+giDwvIu/ijI1eWrlZIpIuIttF5NagGG50f+1/LCLPu/NCbqfY35QgIrPFycv/kYj0dBe9BSS5R0Ldiq0zR0SeEZH33ThS3bg+E5E5QeX6iMhaEflQRF52cychIhPdeD4VkZnundG4f9dj4oyp8EXhfkWktTsv0/0bgzMFmJ8DVbXJppidgE7AJ0BN4BfANuAOd9kcYHDx5+7rnKDnL+IkKANojJMqAuBBYANQw0O594DqQD1gH05K6dY442vUc8udVtp2iv1dfwBmuc/PBb7GydnTFPg0zHsxBycdueCkcjkItMX5obcBZyyEejh5dWq56/wRmBgcn/v8eWCA+zwdeMJ9fhmwwn3+d+B693m1wvfJpp/PZKeYTKzrBryqbg4p9+igvHoB57s/mMFJnFbbfZ6mqrkeyr2hqvlAvoh8i5Oy4BLgZVX9DkBV95e2Hf1pzAKAi3G+gFHVz0XkK6Alzpd+aZaoqorIJ8D/qeonACKyCadxaYgz8NG77v6r4aTsAOgpInfhNLanAZuAJe6ywsR7G9zt4K53rzjjGryiqlvLiM1UMtZAmJ+DKkBXVc0Lnul+gf7osVx+0KyjlF53Qm4nQgrjOFYspmNuTEdxchENDV5JRBKA6UCKqu4SkQdxjliKbzfwt6nqiyLyAc6gN0tF5Leq+t8I/z0mhlkfhIl17wBXulcmJQIDjmMbbwH/r/CFiCSfYLlC/wWuEZG6bvnTyrGd1Tgd64hIS5xTUZEYYe994CIROcfddi13+4WNwXfuUdHgsjYkIs2B7ao6FVgMtItAfOYkYg2EiWnqDHH5EvAxzkh4649jM7cCKW5H62bglhMsVxjbJuBh4G1xsmMWpsD2sp3pQBX3VNFLwEj3FNYJUdW9wEhgnjjZTdcC56rqAeA5nNTOy/D2Pl4LfCrOSG1tgH+faHzm5GLZXI0xxoRkRxDGGGNCsgbCGGNMSNZAGGOMCckaCGOMMSFZA2GMMSYkayCMMcaEZA2EMcaYkP4/XraY8oqbeUYAAAAASUVORK5CYII=\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["print (\"group A  [\", \" \".join([\"%6.2f\"%i for i in A]), \"]  mean %6.2f\"%np.mean(A))\n", "print (\"group B  [\", \" \".join([\"%6.2f\"%i for i in B]), \"]  mean %6.2f\"%np.mean(B))\n", "\n", "\n", "plt.hist( t, density=True, alpha=.5, bins=30);\n", "plt.axvline(np.percentile(t,2.5), color=\"black\", label=\"%5 two sided interval\\n(2.5% - 97.5%) \")\n", "plt.axvline(np.percentile(t,97.5), color=\"black\")\n", "plt.axvline(ttest(A,B), color=\"red\", ls=\"--\", label=\"real world sample A-B\")\n", "plt.axvline(ttest(B,A), color=\"red\", ls=\"--\", label=\"real world sample B-A\")\n", "plt.grid(); plt.legend(); \n", "plt.xlabel(\"difference of means\"); plt.ylabel(\"probability\")\n", "plt.show()"]}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["empirical (simulated) p-value 0.0013\n"]}], "source": ["#A,B = generate_AB()\n", "\n", "k = np.mean( ttest(A,B)<t)\n", "k = (k if k<0.5 else 1-k)*2\n", "print (\"empirical (simulated) p-value %.4f\"%k)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["#### which corresponds to the value obtained with the analytical formula\n", "\n", "As implemented in [`scipy.stats.ttest_ind`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html). With larger simulations, the result gets closer to the analytical formula.\n", "\n"]}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [{"data": {"text/plain": ["0.005726493310989277"]}, "execution_count": 8, "metadata": {}, "output_type": "execute_result"}], "source": ["stats.ttest_ind(A, B).pvalue"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Understanding accuracy of large and small simulations\n", "\n", "observe the relation between the empirical (simulated) and analytical $p_{value}$ is much tighter (more accurate) when the simulation is larger."]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": ["\n", "    \n", "def test_simulations(n, mu=mu, sigma=sigma, real_world_sample_size=real_world_sample_size):\n", "    a = np.r_[[stats.norm(loc=mu, scale=sigma).rvs(real_world_sample_size) for _ in pbar(range(n))]]\n", "    b = np.r_[[stats.norm(loc=mu, scale=sigma).rvs(real_world_sample_size) for _ in pbar(range(n))]]\n", "\n", "    t = np.r_[[ttest(A, B) for A,B in zip(a,b)]]\n", "    r = []\n", "    for _  in range(1000):\n", "        A,B = generate_AB(mu=mu, sigma=sigma, real_world_sample_size=real_world_sample_size)\n", "        ts = ttest(A,B)\n", "        \n", "        k = np.mean( ts<t)\n", "        k = (k if k<0.5 else 1-k)*2\n", "        r.append([k, stats.ttest_ind(A, B).pvalue])\n", "    r = np.r_[r]\n", "    \n", "    plt.scatter(r[:,0], r[:,1], s=2)\n", "    plt.grid(); plt.xlabel(\"empirical p-value\"); plt.ylabel(\"analytical p-value\")\n", "    \n", "    em, an = r[:,0], r[:,1]\n", "    \n", "    print (\"p_value mean error simulated/analytic %.3f\"%np.mean(np.abs(em-an)))\n", "    print (\"p_value mean error when p_value<0.05  %.3f\"%np.mean(np.abs(em[an<0.05]-an[an<0.05])))"]}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100% (100 of 100) |######################| Elapsed Time: 0:00:00 Time:  0:00:00\n", "100% (100 of 100) |######################| Elapsed Time: 0:00:00 Time:  0:00:00\n"]}, {"name": "stdout", "output_type": "stream", "text": ["p_value mean error simulated/analytic 0.165\n", "p_value mean error when p_value<0.05  0.004\n"]}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5xV1X338c8PGJBhFEcxeOUSxSB4Q4i00YcMjVZURPNoEiWJ0WqsT0pjGmMdjZpEeV6S8KQNbeyTGE1i0tqJ1cYiTiWSSpAkIlARw2CAcItIIugIDNcZ+PWPs8949uGcmX1mzj7X7/v1mtc6Z+3L+S0Hz2/22muvZe6OiIhUrz7FDkBERIpLiUBEpMopEYiIVDklAhGRKqdEICJS5foVO4BcDRkyxEeMGNGjY3fv3s2gQYPyG1CJU5urg9pcHXrT5uXLl2939+MybSu7RDBixAiWLVvWo2MXLlxIQ0NDfgMqcWpzdVCbq0Nv2mxmm7JtU9eQiEiVUyIQEalySgQiIlVOiUBEpMopEYiIVLnYEoGZfd/M3jKz32TZbmb2D2a2zsxWmtl5ccUiIiLZxXlF8ENgShfbLwVGBT+3AP8/xlhERCSL2BKBuy8C3ulilyuBH3nCS8DRZnZCXPGIiBTa8k2tXP/oEpZvai12KF2yONcjMLMRwDx3PzPDtnnALHdfHLz/OXCnux/2tJiZ3ULiqoGhQ4eOb2pq6lE8bW1t1NXV9ejYcqU2Vwe1ubj+sHMf23cdYMiR/Tn+qCM66zdu382u/R0cOaAfI4b0/ino3rR58uTJy919QsaN7h7bDzAC+E2WbfOAC1Pe/xyY0N05x48f7z31wgsv9PjYcqU2Vwe1uTCWbXzHP/3IS75s4zuh+lPvetaH3znPT73r2Uj7/8tLm/zcr833f3lpU06f35s2A8s8y/dqMUcNbQFOSXl/clAnIlKSGp98lUVrt9P45Kuh+psvHEm/PsbNF44M1Y8fXs+PbprI+OH1ofoHm1fTuqedB5tXxx5zFMWca2guMMPMmoCJwA5331rEeEREuvSHnftDZVLjZWfQeNkZkc8zuLaGXfs7GFxbk9f4eirO4aP/Cvwa+ICZvWFmN5nZrWZ2a7BLM7AeWAd8D/hcXLGIiOTDXZedQX1tDXfl8KWfycB+fUJlscV2ReDu13Wz3YG/iuvzRUTy7QPHH8lZJw3mA8cf2bsTmYXLFMs3tTJnwRpuu+j0w7qU4lIa6UhEpIRkG/Y5Z8EaFq3dzpwFa3p1/kH9+4bKVA88s4pFa7fzwDOrQvWPL9nM6q07eXzJ5l59diZKBCIiabJ9Gd920elMGjWE2y46vVfnv/eKsUwaNYR7rxh72LbdBw6GyqTZ81+n45Aze/7rvfrsTMpuYRoRkdhl6bpJjgLqra7Os7f9YKhMuuOS0exYv4I7Lhnd689PpysCEZE0904dk/iLfeqYWM7/+JLNjLv/Zxm7eXbsORAqk6ZPHMYZJxzF9InD8h6PEoGISJps4/9zle1ew+z5r9O6pz1jN89dl40JRibFk4QyUSIQkYrX1V/g+ZDrzeU7LhlNfW1Nxm6e6ROH8cp9fx7LX/7Z6B6BiFS8mfNa2NN+kJnzWnr1BfuFpld4esWbXHXuiXzr2nGd9ckvfCDU95+8qZx+c3n6xGEF/aLvjq4IRKQKeFrZM0+veDNUJmUbTZSvLqa4KRGISMW7Z+pY6mtruGfq4cM1c3HuyYNDZVK5fOFno0QgIhUvX08Er/ljW6isFLpHICIVr/Gplax9q42tO/bx/Bc/3OPzHDx0KFRWCl0RiEjFmNW8mtPubmZW2vTOf9ixN1R2J9soo2HHDgqVlUKJQEQqxvdeXE/HIed7L64P1ec6Nn/mvFW07mln5rzwFBOzrj6bSaOGMOvqs/MWcylQ15CIVIyjBvajdU8HRw0Mf7XlOlxzf8ehUJmUrykmSo2uCESk7GSbibPjoIfKKOfJ1AWUPDziacqeEoGIlJ0Hm1voOOQ82NwSqs9XF5CllZVOXUMiUnaOru0P7A/K9+T+xG7mr/y/nPR+Hlm84bA1iCuVEoGIlJ1jB/UPlT11z9QxzJ7/+mFz/uS6BnG5UyIQkbJz1MCaUNlTpTbnT7HoHoGIlJ3kBG/JsqfinpW0XCgRiEjJyja981Xnnhgqe2rmvJbgZnFL9ztXMCUCESlZjU+tZNHa7TQ+tTJU/61rx3HWSYNDU0H3RKVOGZErJQIRKVlbWveEynyr1CkjcqVEICIl65i6AaEy3268YCT1tTXceEF1DBPNRolARErWwJq+oTLfnvvNVlr3tPPcb7bGcv5yoeGjIlKyBvXvGyrzLdtSktVGVwQiUrI+/sFh1NfW8PEPaqx/nJQIRKTosq0jkK+um2zDUJOLzs9ZsKZX5y93SgQiklfZvnS72vbI4g10HHIeWbwhVJ9tUfhcPzvbF36u569UukcgInmV/NIFDpu7P9u2my8cmXGSt1zn/39gXgsrfv8uO/d18PRfXdBZn+1eQKWuL5CrWBOBmU0B5gB9gUfcfVba9mHAY8DRwT6N7t4cZ0wiEq+ubsBm25a3Sd7cw2VAX/hdi61ryMz6Ag8BlwJjgOvMLH2S8HuAJ9x9HHAt8E9xxSMixZf8Qh4/vD6W8997xVgmjRrCvVeMjeX8lSrOewTnA+vcfb27HwCagCvT9nHgqOD1YODNGOMRkQLQDdjyY+7xrMVmZtcAU9z95uD9p4GJ7j4jZZ8TgJ8B9cAg4CJ3X57hXLcAtwAMHTp0fFNTU49iamtro66urkfHliu1uTqUUpvf2X2AP+7cx9CjjuCYXq4X0JVMbf7dtjb2HDhIbf++nHpcafz3yKfe/J4nT5683N0nZNpW7JvF1wE/dPdvmtmfAj82szPdPTQDlLs/DDwMMGHCBG9oaOjRhy1cuJCeHluu1ObqUEptvuqhX7Li9+2ce8qg0A3bfMvU5r/9vwt4a9d+3nfkAF7+WEPG48pZXL/nOLuGtgCnpLw/OahLdRPwBIC7/xo4AhgSY0wiErcsN2wLYfuu/aFSookzESwFRpnZSDPrT+Jm8Ny0fTYDHwEwszNIJIJtMcYkIjEr5tPAJx59RKiUaGJLBO7eAcwA5gOrSYwOWmVm95vZtGC324HPmtmrwL8CN3hcNy1EpCCKOZHbnOvOY9KoIcy57ryCf3Y5i/UeQfBMQHNa3X0pr1uA+DoRRaTgCjGR2/JNrWzcvpvlm1pDQ1H1vEDPaIoJEemRbNM59ORZgVzXDp6zYA279ndoiGqeKBGISI88MK+FRWu380Ae1vv92txVtO5p52tzV0Xa/7aLTufIAf2qfo6gfCn28FERKVd5HB3UHqwZ3B5x7eDxw+vZtWFQbE8oVxtdEYhIl7J12+RzOodp55wYKqWwlAhEpEuz579O6552Zs9/PVTfk3sB2dYdeGf3gVAphaVEICJA9pu/d1wymvraGu64ZHSvzgP5W3dA8kv3CEQEgAeeWcWKN3awc287T8+4sLN++sRhTJ8Y/eGwW3+8jG1tB1i9dSdL77k4tC1f6w5IfikRiAgAuw8cDJU9ta3tQKhMlbd1BySv1DUkIgAMGtAvVPbUpFFDQqWUPiUCEQHg4xNOScwRNOGU7nfuwtZ394ZKKX1KBCJVJttw0HzNEbTl3X2hUkpfpERgZgPN7ANxByMi8cs2HDRfI3dOCmb+PEkzgJaNbhOBmV0BrACeC96fa2bp00mLSJn4xIRT6NfH+ERaF1CuzwVkGyY665pzmDRqCLOuOSdvMUu8otwV+iqJ9YcXArj7CjMb2dUBIlK6WrbupOOQ07J1Z6/Ok224qYaClp8oXUPt7r4jrU5rBoiUqSlnnkB9bQ1TzjyhV+fJ13BTKb4oiWCVmU0H+prZKDP7R+BXMcclIr2UrevmiaWbad3TzhNLo035nO08g/r3DZVSvqIkgr8GxgL7SawithP4QpxBiUjvPfDMqsQ00c+kTe1sFi67cVvTKyxau53bml4J1edz0jkprm7vEbj7HuDLwY+IlIm3gwnc3k6byO3eqWOYs2BN5NFBbwbPA7yZ9lyA7gVUjm4TgZm9QIZ7Au7+Z7FEJCJ58dbO/aEyKdcv8BOPHsgbrXs58eiBeY1PSkeUUUNfSnl9BHA10BFPOCKSL337AAeDMoJZzas7J4RLnQ/ocw2nMXv+63yu4bR4ApWi6/afiLsvT/n5pbt/EWiIPzQR6Y17po6lvraGe6ZG68P/7qL1dBxyvrtofag+X08cS+mK0jV0TMrbPsB4YHBsEYlITpZvamXj9t0s39Qaehgs1+mjs0neS9BaAZUrykXjcmBZUP4auB24Kc6gRCS6B+a1sGt/R68Xkb/y3BNDZVJPViKT8hJl1JCeIhYpYbv3tYfKntJykdUrayIws//d1YHu/u/5D0dEcjXoiJpQmbR8U2vnMNHUv+az1asLqHp1dUVwRRfbHFAiECmgbF/g904dw2tLf8W9U8eE9r+t6RXeaN3L+u27WXzne6O95yxYw6K12wFCw0j1XED1ypoI3P3GQgYiIl1rfPJV1m7bzdZ39/L87Q2d9eOH17Nrw6DD+vC3tO4NlUn6y1/SRVqTzswuJzHNROcE4+5+f1xBicjhcl3wxdPKJP3lL+mirEfwHeATJOYcMuBjwPCY4xKRNCfVDwyV3Uk+SBb1gTKpXlH+iXzI3a8HWt39a8CfArqmFCmwj4x+H/36GB8Z/b5QfepzBKkGBBlggDKBdCPKv5BkB+MeMzsRaAd6N5G5iOTsB7/cQMch5we/3BCqb3xqJbv2d9D41MpQfa5PFkv1inKPYJ6ZHQ3MBv6bRJfj96Kc3MymAHOAvsAj7j4rwz4fJ7EKmgOvuvv0aKGLVJf9Bz1UJv1hx14YGpQp8vVksVS+KA+UPRC8fMrM5gFHZFix7DBm1hd4CLgYeANYamZz3b0lZZ9RwF3ABe7eambvy3w2EcnmrsvGsGP9Cu66bEz3O4tkEOVm8Uozu9vMTnX3/VGSQOB8YJ27r3f3A0ATcGXaPp8FHnL3VgB3fyuX4EUk8Zf/GSccpb/+pcfMvevlh81sOIlRQ58ADgE/AZ5w9y7XuTOza4Ap7n5z8P7TwER3n5Gyz9PAGuACEt1HX3X35zKc6xbgFoChQ4eOb2pqitzAVG1tbdTV1fXo2HKlNleO17a89zfYWSeF532s1DZ3RW3OzeTJk5e7+4RM26J0DW0CvgF8I+jKuRf4Ookv7t7qB4wiMa31ycAiMzvL3d9Ni+Fh4GGACRMmeENDQ48+bOHChfT02HKlNleOGxqf7Xy98ZMNoW2V2uauqM35E/WBstSrgoPA30Y4bAtwSsr7k4O6VG8AS9y9HdhgZmtIJIalUeISEZHei3KPYAnw02Dfj7n7+e7+zQjnXgqMMrORZtYfuBaYm7bP0wSL3JjZEBLPJ6xHpIpd/+gSRjQ+y/WPLil2KFIlojxHcL27nxcM/dwT9cTu3gHMAOYDq0ncV1hlZveb2bRgt/nA22bWArwA3OHub+fWBJHKkpwQLlkmTRo1JFSK5EuUewS/TXnbDJwX9eTu3hwck1p3X8prB74Y/IhUlS80vcLTK97kqnNP5FvXjut2f80PJHHJ9dlziyUKkQr2haZXGNH4LF9oeiVU//SKN0NlkuYIkkLL9Z9apCeKReQ92b7wrwqWhLwqbWnIK84+MVSKxC3qqKHzgAsBN7Pz3P2/4w1LpHLU9e9L24GD1PUPj7j+1rXjMnYJ/WLNtlApErcoo4buAx4DjgWGAD8ws3viDkykUuztOBgqu/Ph048LlSJxi3JF8EngHHffB2Bms4AVwMw4AxOpFAcPhcvu6IpACi3KPYI3SVmZDBjA4Q+GiUiOHl+ymXH3/4zHl4Rna7njktHU19ZwxyWjixSZVJsoVwQ7gFVm9jyJqaIvBl42s38AcPfPxxifSMV6sLmFXfsP8mBzS2jCOE0fLYUWJRH8NPhJWhhPKCLlK9dnAgCOHzyQXW+1cfzgaEtPisQlygNljxUiEJFysHxTK3MWrOG2i05n/PD6zvrUIaLpieDWSe/nkcUbuPnCkaH6WVef3XkukWLSIysiObj1n5exaO12bv3nZaH6bM8EAFw89ng+dOqxXDz2+FD9+OH1/OimiaGEIlIMkZ4jEJGEbbsOhMqkbM8EAMxZsKZz3iBNEyGlSIlAJGbJrh91AUmpypoIzOwZEqOEMnL3adm2iZS7bPcCeiLZBSRSqrq6Ivh/BYtCpMRk686p6Wu0H3Rq+mr+RakcWROBu/+ikIGIlJJjBvUPlUk1fYJE0EeJQCpHlLmGRpnZk2bWYmbrkz+FCE6kWLLNGHrP1LHU19Zwz9SxxQhLJBZRbhb/APgK8PfAZOBGNOxUKsSs5tWdY/wbLzuj2/311K9Uoihf6APd/eeAufsmd/8qcHm8YYkUxiOLN9BxyHlk8YZQ/XF1/UOlSCWLkgj2m1kfYK2ZzTCzjwJ1McclUhBTzz4hVCZtazsQKkUqWZREcBtQC3weGA98CvhMnEGJ5Nus5tWcdnczs5pXh+pXvbkzVIpUo24Tgbsvdfc2d3/D3W9096vd/aVCBCeSq8eXbGb11p2HTe388Ivr6TjkPPxieJzDH3bsDZVJk0YNCZUilSzKqKHnzezolPf1ZjY/3rBEupZtLv8Hm1voOOQ82NwSqq/p0ydUJt112Rjqa2u467Ixofof3TSRjbMu14NgUhWidA0Ncfd3k2/cvRV4X3whiXTvweYWWve0H/aFn81XpiWGfX5lWnjY5/SJw3jlvj/XSCCpalGGjx4ys2HuvhnAzIbTxdQTIoUwuLY/u/bvZXBteFTPrv0HQ2WShn2KZBclEXwZWGxmvwAM+F/ALbFGJdKNgTV9Q2XSuScPBt4OShGJIsrN4ueA84CfAE3AeHfXPQIpiOWbWrn+0SUs39Qaqh/Uv2+oTNr0zp5QKSLdy5oIzGx0UJ4HDCOxiP2bwLCgTiR2ycnf5ixYE6o/amBNqEy645LR9OtjWvhdJAdddQ19kUQX0DczbHPgz2KJSCTFzr3toTIpOTNoskyaPnEYC/eup0H3A0Qi62r20eR9gEvdfV/qNjM7ItaoRAIr3tgRKpOuOvfEzsXiRaR3otws/hWJewTd1YkUTFdLQ4pIbrpaoex44CRgoJmNIzFiCOAoElNOiORNrrOAikj+dDVq6BISq5SdDPwdiXsF3wT+Brg7ysnNbIqZ/dbM1plZYxf7XW1mbmYTooculeR7wRQQ33tRS12IFFpX9wgeAx4zs6vd/alcT2xmfYGHgIuBN4ClZjbX3VvS9juSxMR2S3L9DKkcx9YN4K1d+zm2bkCovr62htY97dTX1mQ5UkR6K8oUE7PNbLaZ5Xq9fj6wzt3Xu/sBEs8gXJlhvweArwP7MmyTKuHBw+qe9tB63YB+oVJE8s/cu54tIviL/VreW5ns+0CTu3c5b6+ZXQNMcfebg/efBia6+4yUfc4DvuzuV5vZQuBL7r4sw7luIXiaeejQoeObmpqitzBFW1sbdXXVtZRCubT5tS3vjQo666T3ngpueXMnB93pa8aYE4+KdK5yaXM+qc3VoTdtnjx58nJ3z9z97u6Rf4APA1uA3cBjwGld7HsN8EjK+08D30553wdYCIwI3i8EJnQXw/jx472nXnjhhR4fW67Kpc3D75zX+ZPqwWdb/NS7nvUHn22JfK5yaXM+qc3VoTdtBpZ5lu/VKNNQ9zWzaWb2U+BbJG4Yvx94Bmju4tAtwCkp708O6pKOBM4EFprZRuBPgLm6YSypWrbupOOQ07JVC8eIxCVKx+ta4AVgtrv/KqX+STOb1MVxS4FRZjaSRAK4Fpie3OjuO4DOVT+66hqSyldb04c97YeorQn/bTLlzBN4bcsOppx5QpYjRaS3otwsPtvdb0pLAgC4++ezHeTuHcAMYD6wGnjC3VeZ2f1mNq3HEUtF2tt+KFQmPbF0M6172nli6eZMh4lIHnT1QNk/Eqw7YGaHbe8qCaTs00xa95G735dl34buzieVq28f6DiUKEOS//Yy/BsUkfzoqmtIXTRSMH3N6CAxOijVvVPHMGfBGm676PQiRSZS+bp7oEykIPr368P+gwfp3y98STB+eL3WDRaJWbc3i83sOOBOYAzQOeuou2saasmbbEtMikj8otws/hcSN3tHAl8DNpIYESSSN30tXIpI4URJBMe6+6NAu7v/wt3/Ai1KIz101bcXM6LxWa769uJQ/YCgS2hAvyj/JEUkn6L8X5dcGmqrmV0eTEl9TIwxSQXLttDMPVPHUl9bwz1TxxYjLJGqFuWBsplmNhi4HfhHEusR/E2sUUnZW76ptXO0z/jh9d3uP33iMKZreUmRoug2Ebj7vODlDmByvOFIpUguOg9o1I9IiYs6auizwIjU/YN7BVLlsv3lnxz3r/H/IqUvStfQfwAvAgsAje2TkAeeWcWKN3awc287T8+4sLNe4/9FykeURFDr7nfGHomUrC77+zUFhEjZizJqaJ6ZXRZ7JFKykv39cxasOWzbvVPHMGnUEO6dOiZUv3xTK9c/uoTlm1pD9ZZWikjxRUkEt5FIBnvNbKeZ7TIzTQ5fRcaccBT9+hhjTjh8hbBkF1D6lULjUytZtHY7jU+tDNX/5aT306+P8ZeT3h9rzCISXZRRQ0cWIhApXT9Z9ns6Djk/WfZ7Gi+LtnT1urfaQmVS42VnRD6HiBRGpMc4zazezM43s0nJn7gDk9JxxyWjqa+t4Y5LRkc+xtNKESldUYaP3kyie+hkYAWJJSV/jaaZqBofOP5IzjppMB84XheHIpUo6j2CDwKb3H0yMA54N9aopKQ8MK+FRWu388C8lmKHIiIxiJII9rn7PgAzG+DurwMfiDcsKSnu4TLFrObVnHZ3M7OaV4fq+1i4FJHSFSURvGFmRwNPA8+b2X8Am+INS0rJ9rb9oTLVI4s30HHIeWTxhlD9zKvOor62hplXnVWQGEWk56KMGvpo8PKrZvYCMBh4LtaopKS88e6+UJnq5gtH8sjiDdx84chQvSaREykfUZ4s7uTuv4grECm+XGcMBQ0HFakEOSUCqWyNT61k7VttbN2xj+e/+OFihyMiBaLloKTT5nd2h0oRqQ5KBPIet3ApIlVBiUA6DTu2NlSKSHVQIpBON14wkvraGm68YGT3O4tIxVAikE5PLN1M6552nli6udihiEgBKRHIe7TIjEhVUiKQTtkWmRGRyqbnCKTTj3+9kUVrt3PMoP6RHygTkfIX6xWBmU0xs9+a2Toza8yw/Ytm1mJmK83s52Y2PM54pGtPr3gzVCb1SStFpLLE9v+2mfUFHgIuBcYA15lZep/DK8AEdz8beBL4RlzxSM/N/GgwgdxHNYGcSCWKs2vofGCdu68HMLMm4Eqgc1J7d38hZf+XgE/FGI/0kCaQE6ls5hnmmM/Lic2uAaa4+83B+08DE919Rpb9vw38wd1nZth2C3ALwNChQ8c3NTX1KKa2tjbq6up6dGy5yqXNr23Z0fn6rJMGxxVS7PR7rg5qc24mT5683N0nZNpWEjeLzexTwAQg40xn7v4w8DDAhAkTvKGhoUefs3DhQnp6bLnK1OZZzas7p45OnTn0hsZnO19v/GT4mHKi33N1UJvzJ877f1uAU1LenxzUhZjZRcCXgWnufvjKJ9Jj1z+6hNe27OD6R5eE6rMtJnPuyYNDpYhUhzivCJYCo8xsJIkEcC0wPXUHMxsHfJdEF9JbMcZSlRat3c4Hz0qUqbItJvP0jAsLGZ6IlIjYEoG7d5jZDGA+0Bf4vruvMrP7gWXuPheYDdQB/2aJp1k3u/u0uGKqNgP6WqhM0mIyIpIq1nsE7t4MNKfV3Zfy+qI4P7/a7T/ooVJEJBM9I1SCHl+ymXH3/4zHl2jyNxGJnxJBEc1qXs1pdzczq3l1qP7B5tW07mnnwbT6XI06blCoFBHJRImgiL6zaD0dh5zvLFofqj/+qAGhsqeev72Bs04azPO3N/TqPCJS2ZQIStCsa85h0qghzLrmnGKHIiJVoCQeKJOw8cPr+dFNE4sdhohUCV0RiIhUOSUCEZEqp0RQAF9oeoURjc/yhaZXih2KiMhhlAgKINuCLyIipUCJoIiuOvfEUCkiUgwaNVRE37p2HN+6dlyxwxCRKqcrAhGRKqdEICJS5ZQIRESqnBKBiEiVUyIQEalySgQiIlVOiSCPrn90CSManz1ssXgRkVKmRJBHyUXi0xeLFxEpZUoEBXB0bb9QKSJSSpQICuDRz5zPpFFDePQz5xc7FBGRw+hP1ALQQjMiUsp0RSAiUuWUCEREqpwSgYhIlVMiEBGpckoEIiJVTqOGsvjgzOfZ1naA4+r6s/Sei0PbLv7mQtZu282o4wbx/O0NxQlQRCRPdEWQxba2A6Ey1dptu0OliEg5UyLogVHHDQqVSQP79QmVIiLloOq7hh5fspnZ81/njktGM33isEjHZOsOWj3z0jxGJiJSGLH+6WpmU8zst2a2zswaM2wfYGY/CbYvMbMRccaTyd0/fY3WPe3c/dPXCv3RIiIlIbZEYGZ9gYeAS4ExwHVmNiZtt5uAVnc/Dfh74OtxxXPmfc/x2pYdnHnfc5H275NWiohUqji7hs4H1rn7egAzawKuBFpS9rkS+Grw+kng22Zm7u75DqbtwMFQ2Z31sy7PdwgiIiXJYvjOTZzY7BpgirvfHLz/NDDR3Wek7PObYJ83gve/C/bZnnauW4BbAIYOHTq+qakp53he27KDoQPhj3vhrJMGd9b/ZssOHDDgzJT6StHW1kZdXV2xwygotbk6qM25mTx58nJ3n5BpW1ncLHb3h4GHASZMmOANDQ05n+OGxme5/awOvvlaPzZ+8r3jcz9TeVm4cCE9+e9VztTm6qA250+cXeBbgFNS3p8c1GXcx8z6AYOBt+MIZuOsyznrpMFsVJePiEhInIlgKTDKzEaaWX/gWmBu2j5zgc8Er68B/iuO+wMiIpJdbF1D7t5hZjOA+UBf4PvuvsrM7geWuftc4FHgx2a2DniHRLIQEZECivUegbs3A81pdfelvN4HfCzOGEREpGsaJi8iUuWUCEREqpwSgYhIlVMiEOVMPSsAAAbPSURBVBGpcrE9WRwXM9sGbOrh4UOA7d3uVVnU5uqgNleH3rR5uLsfl2lD2SWC3jCzZdkesa5UanN1UJurQ1xtVteQiEiVUyIQEaly1ZYIHi52AEWgNlcHtbk6xNLmqrpHICIih6u2KwIREUmjRCAiUuUqMhGY2RQz+62ZrTOzxgzbB5jZT4LtS8xsROGjzK8IbZ5kZv9tZh3B6nFlL0Kbv2hmLWa20sx+bmbDixFnPkVo861m9pqZrTCzxRnWCS873bU5Zb+rzczNrKyHlEb4Hd9gZtuC3/EKM7u51x/q7hX1Q2LK698B7wf6A68CY9L2+RzwneD1tcBPih13Ado8Ajgb+BFwTbFjLlCbJwO1wev/UyW/56NSXk8Dnit23HG3OdjvSGAR8BIwodhxx/w7vgH4dj4/txKvCM4H1rn7enc/ADQBV6btcyXwWPD6SeAjZmYFjDHfum2zu29095XAoWIEGIMobX7B3fcEb18isUpeOYvS5p0pbwcB5T4aJMr/zwAPAF8H9hUyuBhEbW9eVWIiOAn4fcr7N4K6jPu4ewewAzi2INHFI0qbK02ubb4J+M9YI4pfpDab2V+Z2e+AbwCfL1Bscem2zWZ2HnCKuz9byMBiEvXf9dVBl+eTZnZKhu05qcREIBJiZp8CJgCzix1LIbj7Q+5+KnAncE+x44mTmfUB/g64vdixFNAzwAh3Pxt4nvd6N3qsEhPBFiA1Q54c1GXcx8z6AYOBtwsSXTyitLnSRGqzmV0EfBmY5u77CxRbXHL9PTcBV8UaUfy6a/ORwJnAQjPbCPwJMLeMbxh3+zt297dT/i0/Aozv7YdWYiJYCowys5Fm1p/EzeC5afvMBT4TvL4G+C8P7sKUqShtrjTdttnMxgHfJZEE3ipCjPkWpc2jUt5eDqwtYHxx6LLN7r7D3Ye4+wh3H0HiXtA0d19WnHB7Lcrv+ISUt9OA1b3+1GLfJY/pzvtlwBoSd9+/HNTdT+IfCMARwL8B64CXgfcXO+YCtPmDJPobd5O4+llV7JgL0OYFwB+BFcHP3GLHXIA2zwFWBe19ARhb7JjjbnPavgsp41FDEX/HDwa/41eD3/Ho3n6mppgQEalyldg1JCIiOVAiEBGpckoEIiJVTolARKTKKRGIiFQ5JQKpWmb2q55s6+acXzWzL/U8qpw/b6OZDSnU50llUiKQquXuH0qvC540z7hNpFIpEUjJM7NPmdnLwdzr3zWzvkF9m5nNNrNVZrbAzM43s4Vmtt7MpgX73GBm/xHUrzWzr6Scty0oG8zsRTObC7Skbgte3xnM8f+qmc0K6j5rZkuDuqfMrLabNvzQzL5jZsvMbI2ZTc2wz61mNjvl/Q1m9u3g9dNmtjxo6y0Zjh1hZr9Jef8lM/tq8PpUM3suOP5FMxsd5b+7VA8lAilpZnYG8AngAnc/FzgIfDLYPIjE9CBjgV3ATOBi4KMknsRMOh+4msR6DB/LMg/NecBt7n562udfSmIa4Inufg6JGT0B/t3dPxjUrSYxu2l3RgSxXA58x8yOSNv+VBB70idIzBcE8BfuPp7E5HmfN7NcZst9GPjr4PgvAf+Uw7FSBfoVOwCRbnyExKRaS4MlIwYCyXmDDgDPBa9fA/a7e7uZvUbiSzfpeXd/G8DM/h24EEifi+Zld9+Q4fMvAn7gwboG7v5OUH+mmc0EjgbqgPkR2vKEux8C1prZemA0iakgCM69Lbia+RMScwSNBn4ZbP68mSWTxCnAKCJMlGhmdcCHgH9LWXJjQIRYpYooEUipM+Axd78rw7Z2f2+OlEPAfgB3P5Ts6w+kz6OSaV6V3TnG9UPgKnd/1cxuABoiHHNYHGb2A2Ac8Ka7X0biCuDjwOvAT93dzayBREL6U3ffY2YLScyXlaqD8BV+cnsf4N3gakokI3UNSan7OXCNmb0PwMyOsdzXHr44OG4giWmZf9ndASmeB25M3gMws2OC+iOBrWZWw3tdVd35mJn1MbNTSSxF+Ft3v9Hdzw2SAMBPSXRFXcd73UKDgdYgCYwmMdVyuj8C7zOzY81sADAVOlcs22BmHwviNzM7J3rzpRooEUhJc/cWEour/MzMVpL4Yj6h66MO8zKJ/veVwFOewxTF7v4ciWmAl5nZChJ97AD3AktIJJXXI55ucxDLfwK3uvthyyq6eyuJew7D3f3loPo5oJ+ZrQZmkZhqOf24dhL3RV4m8d8oNaZPAjeZ2askZq2MfelDKS+afVQqWtBtM8HdZxQ5jh8C89z9yWLGIZKJrghERKqcrghERKqcrghERKqcEoGISJVTIhARqXJKBCIiVU6JQESkyv0PpWcj1C18z0kAAAAASUVORK5CYII=\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["test_simulations(n=100, mu=mu, sigma=sigma)"]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100% (1000 of 1000) |####################| Elapsed Time: 0:00:00 Time:  0:00:00\n", "100% (1000 of 1000) |####################| Elapsed Time: 0:00:00 Time:  0:00:00\n"]}, {"name": "stdout", "output_type": "stream", "text": ["p_value mean error simulated/analytic 0.172\n", "p_value mean error when p_value<0.05  0.016\n"]}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3zU9Z3v8dcnFy4hCKlBULm7sYjaoqBsLUthqwvFa4/21NpWbWWpe5atbltOaStaL6vZ5bS7nN2etRbd6jmHg11dKdocqfQ0Uus2YlYUuRQQTUSwgEYgXHP5nD9+M+NkMkl+CfllMjPv5+PBg5nfZebzdXA+872buyMiIvmrINMBiIhIZikRiIjkOSUCEZE8p0QgIpLnlAhERPJcUaYD6K7y8nIfP358j+49fPgwQ4YM6d2A+jmVOT+ozPnhZMpcW1u7391HpDuXdYlg/PjxvPzyyz26t7q6mlmzZvVuQP2cypwfVOb8cDJlNrO6js6paUhEJM8pEYiI5DklAhGRPKdEICKS5yJLBGb2iJntNbPXOzhvZvbfzWyHmb1mZhdGFYuIiHQsyhrBT4G5nZz/DFAR+7MA+OcIYxERkQ5ElgjcfR3wfieXXA085oHfAcPN7PSo4hERkfQyOY/gTODtpOe7Ysf2pF5oZgsIag2MHDmS6urqHr1hY2Njj+/NVipzflCZc8u7B4+x/9AJyocOYNQpgxLHoypzVkwoc/eHgIcApk2b5j2dUKEJKPlBZc4PuVbmFTX1LF2zlUVzJnHnL1+nubWQooJWdtw/K3FNVGXO5Kihd4AxSc9Hx46JiOSdpWu20nCkiaVrtjJ/xgSKCoz5MyZQW9fAjQ/XUFvXENl7ZzIRrAZujI0e+mPggLu3axYSEclVK2rqueCeX7Kipp5FcyZRVlLMojmTWDzvHHbcP4/F885h2dptrNu+n3uf3sRb+w9HkhAiaxoys/8DzALKzWwXcBdQDODuDwJVwDxgB3AE+EpUsYiI9CeVVVtY/sKbFBocb3GWrtnKK3f+GTdMH9vmutq6Bg4ea2bK6GFgxqHjzSxbu43Hbpneq/FElgjc/QtdnHfgL6N6fxGR/qi2roEH1+0EoNVI1AJSr1m2dhsHjzaxYdcBZlaUc9ulZ7Nx/YvcNuPsXo9JM4tFRPrIipp6Pv/jf088v+rjZ3D+mcMA2vQDxJuDMEskganjyhhfPoSp48p6Pa6sGDUkIpJtausauPeZzeDOkivPZeq4Mpau2Upzq1NUYDz+tU8kvvA3vnOAhiNNADx2y3RuuzT41R9PAFFTjUBEpJckj/BZtnYbG97+gA27DrBs7TYAPnV2sC/MFR87nanjyrjt0rOZWVHOojmTEr/8AaaOK+OxW6b3SRIA1QhERHpNokmH4Nf8wWPN4J74gn9+2742f8e/8IF2HcV9STUCEZFeMve80ykrKWbuecEv/iVXTOaUwcWJ88lDRPsT1QhERHoo3gQUb8t/9vU9NBxp4tnX9/DRUUO5+ZEaDh1v4eDRJlYtnMEN08dm9Jd/R5QIRES6KT4PYNQpA9n1wTGgfSfvsrXbOHS8JbjBLFOhhqJEICISUrwG8Nsd+2lx2HPgWNpOXoj1ERxtAjOWXDE5k2F3SYlARCSkeGfwmLLB7DlwjPkzJrB43jlpr506roxVC2f0cYQ9o0QgIhJSX4/v7ysaNSQi0oHUlT/7enx/X1EiEBHpwL3PbA5W/nxmc6ZDiZQSgYhIBw4fa2rzd1xf7BHQl9RHICKSIr5b2JCBwVfkkEHFwdpBT28KhoK6s2HXAYBeXxI6E5QIRESSrKip57tPbQSguaU12AsAuPfpTYkv/yljhrcZNprtlAhERJIsXbM18fjT54xk9YbdtAKjywYnNohZcsXknOowViIQkbyTbonouEVzJiU2kV+6ZiutseMHjpzghW//aWYCjpg6i0Ukr9TWNTD/0fXtloiO++iooZx/5jA+Omooi+ZMoqS4gJLiQr4zr3/PDj4ZqhGISF5ZtnYbDUeaGDqwiLNGDGnTzr+ipp47f/46za0OBB3B/XGRuN6mRCAieaWz2cH3PbOZ5lanwMiZjuAwlAhEJK8kLwwXFx8u2tQSrBY6qKggpzqDu6I+AhHJWZVVW/ij71ZRWbWl0+vue2ZTsGewGWUlxdxxxbl9FGH/oBqBiOSc+C/8A0ebaHVY/sKbHa4SWlm1hSNNrYnnr9z5Z30VZr+hRCAiOSO+X8DLbzVwpKmF4kKjwGH+jAkdXv/gup2J58UF/XsDmagoEYhITqita0hsDTmwMPhCH3/qEJ77xqc6vH7+o+sTzwcWWt41CcUpEYhIVovXAg4ea05sDTn21CGcPmxQhyN/4kmg4UgTZSXFLL/porzqHE6lRCAiWenIiRZufLiGg0eb2LDrAFNGDwu9BMTiJ16l4UgTJcUFeZ8EQIlARLJQbV0Dde8dZt32Y20WgOvqCz1ee3gntuF8YUF+DRPtiBKBiGSVeLPOV89yykqKu7UAXHzP4YoRQ9h/+ASL5kyKONrsoEQgIv1e/Jf8bZeenVgioqjAutWsU1vXwMFjzUwZPazdQnP5LtJEYGZzgWVAIbDc3StTzo8FHgWGx65Z7O5VUcYkItkluWMXPlz6Ydypjd36Ml+2dhsb3v6AmRXlSgIpIksEZlYI/Ai4DNgFrDez1e6evPnnHcDP3P2fzWwyUAWMjyomEcke8Ulh5aUDE6N74v0Aj90ynerq6k7vT11qOnmNIWkryhrBxcAOd98JYGYrgauB5ETgwCmxx8OA3RHGIyJZorauIbEKaHNLa+jO4LjUVUSXrd3GY7dMz4ltJaMQZSI4E3g76fkuIPVT+D7wSzP7K2AIcGmE8YhIlli2dhvNrU5RgfGdeZNDLwUd70t4pf4DmludQoPzzxymWkAXzN2jeWGz64C57j4/9vzLwHR3X5h0zTdiMfzAzD4BPAyc5+6tKa+1AFgAMHLkyKkrV67sUUyNjY2Ulpb26N5spTLnh1wo85ETLew9eIzTThkEkHhcMqAw7fWpZX7/8Al2f3AMxxlUVEhzaysjTxnER4YM6JP4+8LJfM6zZ8+udfdp6c5FWSN4BxiT9Hx07FiyW4C5AO7+72Y2CCgH9iZf5O4PAQ8BTJs2zWfNmtWjgKqrq+npvdlKZc4P2V7mFTX13PnL12luhZkVpaGacKqrq9k9eCIPVG1h1CkDefegc+h4IUUFxuNf+0ROdghH9TlHuQz1eqDCzCaY2QDgemB1yjX1wKcBzOwcYBCwL8KYRKSfSe4PKCqw0M04R060cOfPX+fQ8Wa27zvMqFMGUlZSzD1Xn5eTSSBKkdUI3L3ZzBYCawiGhj7i7pvM7B7gZXdfDXwT+ImZ/TVBx/HNHlVblYj0K8lrBMWTQHe+xPcePEZzKxQaTCwfQuV1H1cC6KFI5xHE5gRUpRy7M+nxZuCTUcYgIv1PsFLoSxw63kzFaaXdHhUEcNopg5hZUdrt+6Q9zSwWkT5VWbWlzR4AQwYUtusTqKzawvIX3mT+jAlcdu6oxKzi5C/8kjT3Sc8oEYhIn1lRU98mCcSXe0iWvFnM8hfeZPOeg6zbvh9AX/wRUSIQkT4R7xSOu3XmxDbbR9bWNXDv05vY9ofGxLF4jQA0IzhKSgQiErnkmb7xTuHUSWLL1m5jw64DiedTRg9LJArVBKKlRCAikbvvmU00tzoFRtox/vGVQStGDAEzhgwobNdkJNFRIhCRPhDsITyoqDDtCJ97n9nMhrc/YMqY4az6Sw0k7GtRTigTEQHgjismU1ZSzB1XTE4cq61r4MaHa6ita4D49CFNI8oI1QhEpNfFl5BeNGcSN0wfm/iTbPGTr7F9byN7Dhyj8tqPJYaISt9TjUBEet3dT2+i4UgTdz+9qcNr3j1wNPF3fI8BTQzLDCUCEelVtXUNHG8OFhBuam7t8LrvzAuai74zb3KH10jfUNOQiPSaFTX1LFm1MfF8wcyJHV6brrlIMkOJQER6xYqaer771IdJIHkegPRvoRKBmQ0Gxrr77yOOR0SyTLxj+HhTS+JYxWmlmgeQRbpMBGZ2JfDfgAHABDObAtzj7ldFHZyI9E8rauq5e/UmmlpbKS4s4HhzK8WFwVyBa6acwT9cf0GGI5TuCNNZ/H2Cjeg/AHD3DcCECGMSkX5u6ZqtHG9ppdWDDuGykmKKC4JE8Pw27S2VbcIkgiZ3P5ByTLM+RPLYojmTGFhYgAFnlA1m+U0XcccV51JWUsyiOZMyHZ50U5hEsMnMbgAKzazCzP4ReDHiuESkn6ita+CyHz7P+Xc9y4qaeiAY8fP7v/kMf1JRzq6Goyxbu40bpo/llTv/TCOBslCYRPBXwLnAceD/AAeB26MMSkT6j/gM4EPHW3igavOHy0IQLA0d311MsleXicDdj7j799z9InefFnt8rC+CE5G+taKmngvu+WXil/+Kmnq27/1wf4BhJQNYt30/y9ZuA9CM4BwRZtTQr0nTJ+DufxpJRCKSMXc/vYnjza3c/fQmbpg+lqVrtrY5X146kInlQ1QDyDFh5hF8K+nxIOBaoDmacEQkk5pSloZYNGdSYkOZspJillwxWb/+c1CXicDda1MO/dbMXoooHhHpQ7V1Ddz7zGZwZ8mV57Jg5sTEpvEQdAp/dNTQtJvHS+4I0zT0kaSnBcBUYFhkEYlI5GrrGli2dhsHjzWz4e0PgGCryMdumd5uWYh4P4DkrjBNQ7UEfQRG0CT0JnBLlEGJSO+rrNrCT36zkzOGD+bUIQPYsOsAU0YPY8qY4eDObZeenUgQ+vWfX8I0DWkWsUiWq6zawoPrdgLwdsNRTi0dmBj2Gf/Cr61rYP6j62k40gRow/h80mEiMLP/1NmN7v5vvR+OiEThod/sTDweUzY4bafvvc9spuFIE0MHFmlUUJ7prEZwZSfnHFAiEMkCt698hdbYAPCBhQX85tsdjPyO7Rd81oghahbKMx0mAnf/Sl8GIiK9r7augVUbdiee33VVx0tDL7nyXO0bnKfC7kdwOcEyE4Pix9z9nqiCEpGTU1vXwK3/82X2NZ5IHLtmyhlp1wFK3mhe/QL5Kczw0QeBEmA2sBy4DtA8ApF+JvkL/dnX97RJAk/+xSVpm3tq6xoSE8aWrtmqBePyVJhF5y5x9xuBBne/G/gEoLqjSD+zdM1WGo40sXTNVm679GxGlA4AgppAR0lg/qPraW51igpMy0fnsTBNQ0djfx8xszOA94DTw7y4mc0FlgGFwHJ3r0xzzX8m2PzGgVfd/YYwry2Sz2rrGrj36U3sPnCU9xpP8DefKGTRnMmJGsHUcWWsv+Oydvck9wHEh4qWlRSz/KaL1EGcx8IkgmfMbDiwFPgPgi/sn3R1k5kVAj8CLgN2AevNbLW7b066pgL4DvBJd28ws9N6UAaRvBHfIvJ4S2ub4/sPnWDhVWM77QMoLx3YZiVRJQGJCzOh7N7YwyfN7BlgUJody9K5GNjh7jsBzGwlcDWwOemaPwd+5O4Nsffa253gRfJFfFZwq7ddCvi0oQN4r/EE5UMHtLsnXmvY+M4BWhyaW1rb7R2gGcQCYO6d7zppZq8BK4HH3f2N0C9sdh0w193nx55/GZju7guTrlkFbAM+SdB89H13fzbNay0AFgCMHDly6sqVK8OG0UZjYyOlpaU9ujdbqczZ7d2Dx9h/6ASeshK8AeVDBzLqlGAgX7oyv7X/MIeON8euN84YPoiPDGmfMLJVLn3OYZ1MmWfPnl3r7tPSnQvTNHQl8HngZ2bWCjwO/Mzd63sUTfv3rwBmAaOBdWZ2vrt/kHyRuz8EPAQwbdo0nzVrVo/erLq6mp7em61U5uxVWbWFB1/cSfAbCQoNTh82iGVfuLDdr/h0ZY7XCDDLyeWjc+Vz7o6oyhymaagO+Dvg72Jt+kuAvyX+r7Nj7wBjkp6Pjh1Ltguocfcm4E0z20aQGNaHC18kd/0kaVmIgUXG7++b1+U9qYvGrVo4I8oQJUeEnVA2jqBW8HmgBfivIW5bD1SY2QSCBHA9kDoiaBXwBeBfzKycYFjqTkTy3Iqa+sSyEAbcdeV5nV5/5EQLNz5c02ZZaU0Ok7C6nEdgZjXAU7FrP+fuF7v7D7q6z92bgYXAGmALQXPSJjO7x8yuil22BnjPzDYDvwYWuft7PSyLSFarrWtIbAy/dM1WHCgqMJ74i0vSjgZKvn7vwWOs274f3LWZvHRbmBrBje7+ewAzG9WdF3f3KqAq5didSY8d+Ebsj0jeSl4mGoItIpPnBKSzbO224MsfuH7MIGZWlGoUkPRImD6C3yc9rQIujC4ckfwS3yoy3pwDHw7p7Gq5h/iv/tsuPZtDb76qpiDpsVB9BEkskihE8lDqRjAAt86cGPoXffIWktVvRhKi5InuJoIuZxSLSOeS9wuObwRz1oghLLnyXDXrSEaEHTV0ITADcDO70N3/I9qwRHLXvU9vYsOuA1SMGNJuu0iRTAgzauhO4FHgVKCcYKjnHVEHJpKzLGhhHTKomMduma4kIBkXZhnqLwIXuftd7n4X8MfAl6MNSyS3JA/1XHLFZGZWlLPkismZDksECNc0tJtgZ7JjsecDaT9DWETSSLdT2GO3TNcIH+lXwiSCA8AmM3uOYOHDy4CXzOy/A7j71yOMTyRrraipZ8mqjbQkrReniV7SH4VJBE/F/sRVRxOKSPaLjwj6yJABbTaNh453ChPJtDATyh7ti0BEsl3q7OC4ihFDqLzu40oC0m91dx6BiKSorNrCj9ftbLNjwDVTzuCZ1/Ywf8YEFs87J+19K2rqeaBqM6OGDaby2o8pUUjGKBGInKTUJHDrzIksnncO/3D9BZ3et3TNVg4db+HQ3kaWrd2mDmTJmDDDR0UkpHgS6Exl1Rb+6LtVnH/mMIYOLKTitFJ1IktGdVgjMLOngQ73sXT3qzo6J5LL4pvBL5oziRumj+VrMyey/IU3O20Ggg87kn+7Yz8tDi++8R477u96sxmRqHXWNPTf+iwKkSyxoqae7z61EQiadm6YPpbF887pshaQvMDc6LLBvHvgGPNnTOiLkEW61GEicPfn+zIQkf6qsmoLP/nNTk4tHcB7SRPDFs2ZFPo1lq3dRsORJspKill2/QXqGJZ+pcvO4tg+xQ8AkwlmGAPg7hMjjEsko25f+QqrNuzmmiln8PSru2lx2HsoSAJFBcY9V5/X5X4ByfsHJ+8doCQg/U2YUUP/AtwF/D0wG/gK6mSWHBefDLZqw25Glw1mV8NRiguNc08/pdPlopO//JN3ENOyEtKfhUkEg939V2Zm7l4HfN/MaoE7u7pRJJskdwKPGDqAfYdOMGLoAJZdf0Hiy72jBBC/t7x0INv3NgJtdxAT6c/CJILjZlYAbDezhQQLzpVGG5ZI36ms2sLyF96ksMA43tzK0jVbWX7TRW2+/Lv6Nb90zVYajjTR3NLaZo8B1QIkG4Rp4rkNKAG+DkwFvgTcFGVQIn0lvixEc6vT1NxKWUlxYsP47uwV8PlpYygqML44fZz2GJCsE2atofWxh40E/QMiWS/dpvELQkwG68jmPQdpbnU27znYWyGK9Jkwo4aeAz7n7h/EnpcBK919TtTBiURhRU09d/78dZpbP5wvGWZGcGfUHyDZLEwfQXk8CQC4e4OZnRZhTCKRSK0FGPDx0cN6tGl88uigqePK1B8gWS1MImg1s7HuXg9gZuPoZOkJkf7osh9Us33f4TbHDFi1cEaPXi91aKhINguTCL4HvGBmzxP8v/MnwIJIoxLpJfFhnQ1HmhLHRgwdwPuNJ/jzP+n5nEg1BUkuCdNZ/KyZXUiwaT3A7e6+P9qwRE5OfFmIlpS6a8WIITz3zVkn/fpqCpJc0tnqo5PcfWssCUCwiT3A2FhT0X9EH55I971/+AQP/qbtTmEFBjsfuDxDEYn0b53VCL5B0AT0gzTnHPjTSCIS6aF4B+7Fg44BhQAMHVjIkRMtXTYDpS4tLZJPOlt9NN4P8Bl3P5Z8zswGpblFpM/FZwVfctapvPjGezS3OhedH7QHFRhsvHtuqNd5oGoLh44380DVFiUCyTthZha/GPKYSJ/7cWxW8Lrt+2ludYoKjOElAygqMBZ0ozN40ICCNn+L5JPO+ghGAWcCg83sAoIRQwCnECw50SUzmwssI6inL3f3yg6uuxZ4ArjI3V8OH77ko+Qx/Ml9wfHlIc44upMd988K/XoraurZF1ti+tiJlt4NViQLdNZHMAe4GRgN/DDp+EHgu129sJkVAj8CLgN2AevNbLW7b065bijBekY13Ypc8k5tXQO3r3yFtxuOJo5dM+WMxL4B8c3iq6t3dvQSaS1dszXx+DvzJvdOsCJZpLM+gkeBR83sWnd/sgevfTGww913ApjZSuBqYHPKdfcCfwss6sF7SJ648eGaxASuuPis3ngC6I7kWsWiOZPUUSx5zdw7nyRsZjuBJ4FH3H1L6Bc2uw6Y6+7zY8+/DEx394VJ11wIfM/drzWzauBb6ZqGzGwBsUlsI0eOnLpy5cqwYbTR2NhIaWl+raCd7WV+9+Ax9h863m4q+4ihAxl1SvoxC2HK/Ma+Ro6caKFkQCFnjcje/z5x2f4594TK3D2zZ8+udfdp6c6FmVn8ceB64OHYvgSPECw6d1LLLMZe64cEzU+dcveHgIcApk2b5rNmzerRe1ZXV9PTe7NVtpd5/OJfkPzPtLAAfva1SzpdG6izMsdrAnsOGNv3NjJlzHBWfe6TvRx138v2z7knVObe0+UQCXc/5O4/cfdLgG8TbFu5x8weNbM/6uTWd4AxSc9Hx47FDQXOA6rN7C2CmcurzSxtxpL8UlvXwGU/fL7NsYrTSrtMAl2JrxE0ZGARMyvKWXKF+gREwixDXQhcTrAXwXiCCWb/m2DNoSqgo8VW1gMVZjaBIAFcD9wQP+nuB4DypPeppoOmIck/i598LbHlI8DMivIeLemQOlFMm8iLtBemaWg78Gtgqbsnzx94wsxmdnSTuzfHtrZcQzB89BF332Rm9wAvu/vqkwlccldl1ZY2SeD+z57f407c+IJzS9ds5YbpY7VGkEgaYRLBx9y9Md0Jd/96Zze6exVBrSH5WNpN7919VohYJAclj+B5btO7PLjuw+GfJcUF3UoCK2rqObDnILtr6rlh+tg2I4JEJL3OJpT9I7F9B8ys3fmukoBIWH/xv2rZe+g4W989xPuHTySODyw07rji3G691tI1W/nqWZ6oAcT/iEjHOqsRqK1e+sTeQ8cTf986cyLLX3iT+TMmhNo6MnWnsEVzJnFg5wbVAES6oasJZSKRuX3lK6zasLvNscXzzulWAjh4rDmx9eRjt0znhuljqT66k1mqBYiEFmbU0AiCYaOTgcQMHnfXMtRyUlKTwK0zwy8SFx8GOmX0MGZWlGunMJGTEKaz+H8DjxMMIb0VuAnYF2VQktviQzqThR0ZFK8JzD3vdEDDQEV6Q5hEcKq7P2xmt7n788DzZrY+6sAk99TWNbSZH2AEoxFGlw0O3aGrTeNFel+YRBDf9XuPmV1OsGXlR6ILSXJNfPOY04cNarNy6NdmTmTznoNpm3VSO4HjtGm8SO8LkwjuM7NhwDeBfyTYj+CvI41KcsqP1+3EgbcbjlJxWinvHjjKd+ZN7rQW0NEvf00IE+l9XSYCd38m9vAAMDvacCRXdLR5zHPf+FSo+/XLX6TvhB019OcE6wwlrnf3r0YXlmS7L/3kdxxtbmX9m+8n+gLaT0vsmH75i/SdMBu0/hwYBqwFfpH0RyStyqotHG1uBeBocytfmzmRogLja90YHioifSdMH0GJu3878kgkJ6TuJFYxYkjoSWIikhlhagTPmNm8yCORrFdb19AmCdw6cyLPfXNWu2tufLiG2rqGPo5ORDoSJhHcRpAMjprZQTM7ZGYntTuZ5KabHqlJPC4rKUrUApK//OOjgZat3ZapMEUkRZhRQ0P7IhDJTpVVW9osGx23/KaLE4+Th4JqNJBI/xOmjwAzKwMqaLvW0LqogpL+L/7r/rc79rc7N2X0sA4ngWk0kEj/E2b46HyC5qHRwAaCvYX/HdCic3liRU09D1RtZtSwwXx60mk8/vLblJcOZPveRkaXDWZXbLbwlNHDOGVwcbtf+/ryF+nfwtQIbgMuAn7n7rPNbBJwf7RhSaYFX/5bGFZSzO4PjtLqcGhvI2/uP0xzq9Pc0ppY9TP+6z95xJC++EWyR5hEcMzdj5kZZjbQ3bea2Ucjj0wyYkVNPUt+vpGWYBoAh443J86VFBdy4yfG8fjLb7NoziQ+Ompom/WA1P4vkp3CJIJdZjYcWAU8Z2YNQF20YUmmLF2zNZEEIFgZdHBRAUMGFbPkislMHVeWGA2UWgNQE5BIdgozauizsYffN7NfE8wyfjbSqKRPxZeHfvfAUT59zkiefm03La0ws6K80y921QBEckOoUUNxsf0IJEfU1jVw79ObeGPf4UQT0PPb9vHG/ZeHul81AJHc0K1EILkh3qRTWECiGaikuJDCArrc9L2jfQJEJHspEeSJ2roG7n1mM4ePNbF932EgSAJDBxZy1mlDE+3/XdEOYSK5R4kgx23/QyM3L/4FAwuN4y3e5lxZSTHLb7qoW7/s1S8gknuUCHLcseYWoCiRBEqKCzh71CmhawCp1C8gknuUCHJUvB/gm+d/eCx1ApiICIRbfVSyzIqa+jbLQUOwO1h8rH8qLQ0tkt9UI8gxqRvDQJDtF3SyO5g6gEXymxJBDkndGOb+z57PGUd3svOLszq9Tx3AIvlNTUM55OakjWEAbpg+NtR98Q5g9R2I5KdIawRmNhdYBhQCy929MuX8N4D5QDOwD/iqu2sdo26orNrC8hfeZMKpJRw63pI4PrOiPINRiUg2iaxGYGaFwI+AzwCTgS+Y2eSUy14Bprn7x4AngL+LKp5cFN8drLnVE5PEINgXQG39IhJWlDWCi4Ed7r4TwMxWAlcDm+MXuPuvk67/HfClCOPJGStq6lm6ZisfHGlqd+6aKWfwD9dfkIGoRCRbmbt3fVVPXtjsOmCuu8+PPf8yMN3dF3Zw/T8B77r7fWnOLQAWAIwcOXLqypUrexRTY2MjpaWlPbq3v//tY1UAAAxHSURBVHj34DH2HTre7viIoQMZdcqgdsdzoczdpTLnB5W5e2bPnl3r7tPSnesXo4bM7EvANOBT6c67+0PAQwDTpk3zWbNm9eh9qqur6em9mZS80NstD75Ii3/4sXU1SSxby3wyVOb8oDL3nigTwTvAmKTno2PH2jCzS4HvAZ9y9/Y/dfPY7StfYdWG3ZQOKKTxRAt7PjjK6cPb7hGc3BeglUFFpCeiTATrgQozm0CQAK4Hbki+wMwuAH5M0IS0N8JYskb8y3zueaezasNuABpPBKOB3j14nJ9+9eIOv+w1MUxEeiKyRODuzWa2EFhDMHz0EXffZGb3AC+7+2pgKVAK/KuZAdS7+1VRxdTfXfNPL7Bh1wEAana+nzheOqCQ4qICFs2Z1G7Rt+RagCaGiUhPRNpH4O5VQFXKsTuTHl8a5ftnk9tXvpJIAgBNra1MGTMc3Fly5bkdNvWk1gJUExCR7uoXncX5Lt36QAv+ZGJik/hUqgWISG9SIugHUpPArTM7TgKgWoCI9C4lggyqrNrCj9ftbHNsyuhhnSYB0CJxItK7lAgyYEVNPXevfr3d1pEDC41VC2d0eb92CROR3qREkAHfe2ojqfO5x5QN1tIQIpIRSgR97MaHa9okgdHDB7HsCxdqApiIZIwSQR+66L7n2Nd4IvG8YsQQnvvmrMwFJCKCNqbpM5f9oLpNEphZUa4kICL9ghJBH7h95Stt9gsoLtASECLSfygRRCy+cFyylV+7hNq6Bm58uIbauoYMRSYiElAiiNCKmvp2SWD08EFMHVeWmBS2bO22DEUnIhJQZ3EE4ltIpkqeMaxJYSLSXygR9LIpd6/hg6PN7Y7PrChvM2NYk8JEpL9QIuhFMyp/1S4JzKwo1xe+iPRrSgS9IN3qocMHF7HhrjkZikhEJDx1Fp+k2rqGdklg9PBBSgIikjVUIzgJ6YaG3v/Z87lh+tgMRSQi0n1KBD2QLgEAjCgdoCQgIllHiaCbZlT+il0fHGt3fMroYaGWkBYR6W+UCLopNQlcM+UMLR8tIllNncUn4a3Ky/nyJ8ZrqQgRyWpKBCdJS0WISLZT01BIHXUQa6kIEcl2SgSdWFFTz3ef2pj2XFGsLqWlIkQk2ykRpHHZD6rb7B+QSjuLiUguUSJI0lUCAI0SEpHck/eJoLaugWv/+cUOz2vNIBHJdXmdCDqaHDZ8cBGNx1uYP2NCm6WjRURyUV4mghU19SxZtZEWb39OawWJSL7Ju0QwcfEvaE1zXE1AIpKv8mpC2cZ3DigJiIikiDQRmNlcM/u9me0ws8Vpzg80s8dj52vMbHxUsZx357Ptjt06cyJvVV6uJCAieS2yRGBmhcCPgM8Ak4EvmNnklMtuARrc/Y+Avwf+Nqp4Gk+0tHl+/2fPV0ewiAjR1gguBna4+053PwGsBK5OueZq4NHY4yeAT5uZRRgTAKUDCtUhLCISY+5phs70xgubXQfMdff5sedfBqa7+8Kka16PXbMr9vyN2DX7U15rAbAAYOTIkVNXrlzZ7Xg2vnOAkYPhD0fh/DOH9bRYWaexsZHS0tJMh9GnVOb8oDJ3z+zZs2vdfVq6c1kxasjdHwIeApg2bZrPmjWr26/x1cW/4K/Pb+bvNxax84vdvz9bVVdX05P/XtlMZc4PKnPvibJp6B1gTNLz0bFjaa8xsyJgGPBeFMHsrLyc888cxs7Ky6N4eRGRrBVlIlgPVJjZBDMbAFwPrE65ZjVwU+zxdcD/86jaqkREJK3ImobcvdnMFgJrgELgEXffZGb3AC+7+2rgYeB/mtkO4H2CZCEiIn0o0j4Cd68CqlKO3Zn0+BjwuShjEBGRzuXVzGIREWlPiUBEJM8pEYiI5DklAhGRPBfZzOKomNk+oK6Ht5cD+7u8KreozPlBZc4PJ1Pmce4+It2JrEsEJ8PMXu5oinWuUpnzg8qcH6Iqs5qGRETynBKBiEiey7dE8FCmA8gAlTk/qMz5IZIy51UfgYiItJdvNQIREUmhRCAikudyMhGY2Vwz+72Z7TCzxWnODzSzx2Pna8xsfN9H2btClHmmmf2HmTXHdo/LeiHK/A0z22xmr5nZr8xsXCbi7E0hynyrmW00sw1m9kKafcKzTldlTrruWjNzM8vqIaUhPuObzWxf7DPeYGbzT/pN3T2n/hAsef0GMBEYALwKTE655r8AD8YeXw88num4+6DM44GPAY8B12U65j4q82ygJPb4L/Lkcz4l6fFVwLOZjjvqMseuGwqsA34HTMt03BF/xjcD/9Sb75uLNYKLgR3uvtPdTwArgatTrrkaeDT2+Ang02ZmfRhjb+uyzO7+lru/BrRmIsAIhCnzr939SOzp7wh2yctmYcp8MOnpECDbR4OE+f8Z4F7gb4FjfRlcBMKWt1flYiI4E3g76fmu2LG017h7M3AAOLVPootGmDLnmu6W+Rbg/0YaUfRCldnM/tLM3gD+Dvh6H8UWlS7LbGYXAmPc/Rd9GVhEwv67vjbW5PmEmY1Jc75bcjERiLRhZl8CpgFLMx1LX3D3H7n7WcC3gTsyHU+UzKwA+CHwzUzH0oeeBsa7+8eA5/iwdaPHcjERvAMkZ8jRsWNprzGzImAY8F6fRBeNMGXONaHKbGaXAt8DrnL3430UW1S6+zmvBK6JNKLodVXmocB5QLWZvQX8MbA6izuMu/yM3f29pH/Ly4GpJ/umuZgI1gMVZjbBzAYQdAavTrlmNXBT7PF1wP/zWC9MlgpT5lzTZZnN7ALgxwRJYG8GYuxtYcpckfT0cmB7H8YXhU7L7O4H3L3c3ce7+3iCvqCr3P3lzIR70sJ8xqcnPb0K2HLS75rpXvKIet7nAdsIet+/Fzt2D8E/EIBBwL8CO4CXgImZjrkPynwRQXvjYYLaz6ZMx9wHZV4L/AHYEPuzOtMx90GZlwGbYuX9NXBupmOOuswp11aTxaOGQn7GD8Q+41djn/Gkk31PLTEhIpLncrFpSEREukGJQEQkzykRiIjkOSUCEZE8p0QgIpLnlAgkb5nZiz0518Vrft/MvtXzqLr9fm+ZWXlfvZ/kJiUCyVvufknqsdhM87TnRHKVEoH0e2b2JTN7Kbb2+o/NrDB2vNHMlprZJjNba2YXm1m1me00s6ti19xsZj+PHd9uZnclvW5j7O9ZZvYbM1sNbE4+F3v87dga/6+aWWXs2J+b2frYsSfNrKSLMvzUzB40s5fNbJuZXZHmmlvNbGnS85vN7J9ij1eZWW2srAvS3DvezF5Pev4tM/t+7PFZZvZs7P7fmNmkMP/dJX8oEUi/ZmbnAJ8HPunuU4AW4Iux00MIlgc5FzgE3AdcBnyWYCZm3MXAtQT7MXyug3VoLgRuc/ezU97/MwTLAE93948TrOgJ8G/uflHs2BaC1U27Mj4Wy+XAg2Y2KOX8k7HY4z5PsF4QwFfdfSrB4nlfN7PurJb7EPBXsfu/BfyPbtwreaAo0wGIdOHTBItqrY9tGTEYiK8bdAJ4NvZ4I3Dc3ZvMbCPBl27cc+7+HoCZ/RswA0hdi+Yld38zzftfCvyLx/Y1cPf3Y8fPM7P7gOFAKbAmRFl+5u6twHYz2wlMIlgKgthr74vVZv6YYI2gScBvY6e/bmbxJDEGqCDEQolmVgpcAvxr0pYbA0PEKnlEiUD6OwMedffvpDnX5B+ukdIKHAdw99Z4W39M6joq6dZVOdzNuH4KXOPur5rZzcCsEPe0i8PM/gW4ANjt7vMIagD/GdgKPOXubmazCBLSJ9z9iJlVE6yXlayZtjX8+PkC4INYbUokLTUNSX/3K+A6MzsNwMw+Yt3fe/iy2H2DCZZl/m1XNyR5DvhKvA/AzD4SOz4U2GNmxXzYVNWVz5lZgZmdRbAV4e/d/SvuPiWWBACeImiK+gIfNgsNAxpiSWASwVLLqf4AnGZmp5rZQOAKSOxY9qaZfS4Wv5nZx8MXX/KBEoH0a+6+mWBzlV+a2WsEX8ynd35XOy8RtL+/Bjzp3Vii2N2fJVgG+GUz20DQxg6wBKghSCpbQ75cfSyW/wvc6u7ttlV09waCPodx7v5S7PCzQJGZbQEqCZZaTr2viaBf5CWC/0bJMX0RuMXMXiVYtTLyrQ8lu2j1UclpsWabae6+MMNx/BR4xt2fyGQcIumoRiAikudUIxARyXOqEYiI5DklAhGRPKdEICKS55QIRETynBKBiEie+/+MTq4L3b6OrAAAAABJRU5ErkJggg==\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["test_simulations(n=1000, mu=mu, sigma=sigma)"]}, {"cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [{"name": "stderr", "output_type": "stream", "text": ["100% (6000 of 6000) |####################| Elapsed Time: 0:00:02 Time:  0:00:02\n", "100% (6000 of 6000) |####################| Elapsed Time: 0:00:02 Time:  0:00:02\n"]}, {"name": "stdout", "output_type": "stream", "text": ["p_value mean error simulated/analytic 0.168\n", "p_value mean error when p_value<0.05  0.009\n"]}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3xU9Z3/8dcnVxPCJTaIKCDahiKiokTR1qWwxZW13rpqtdaqrdTarb/Sbeuv3q3gVlq23WW3/VWpWrW7abS6UnSpVLpNrdVGSEWRi0DRRLwUaIdLuIRcPr8/zsw4CYEMISeTmfN+Ph48mDnnzMzny0A+nO/l8zV3R0REoisv0wGIiEhmKRGIiEScEoGISMQpEYiIRJwSgYhIxBVkOoCDVVFR4aNHj+7Ra3fu3MmAAQN6N6B+Tm2Ojii2W21OX319/RZ3H9rVuaxLBKNHj2bZsmU9em1tbS1Tpkzp3YD6ObU5OqLYbrU5fWbWsL9z6hoSEYk4JQIRkYhTIhARiTglAhGRiFMiEBGJuNASgZk9aGabzOy1/Zw3M/t3M1tvZq+a2alhxSIiIvsX5h3BQ8D0A5z/e6Ay/us64EchxiIiIvsRWiJw9+eAvx7gkguBRzzwB2CImQ0PKx4RkWxUXdfIKbN+RXVdY2ifYWHuR2Bmo4Gn3X18F+eeBua4+/Px578Gvunu+6wWM7PrCO4aGDZs2MSampoexdPU1ERZWVmPXput1OboiGK7o9Dm1e9up7XdKcgzjh8+qMdtnjp1ar27V3V1LitWFrv7fGA+QFVVlfd0JaFWIUZDFNsM0Wx3FNr8Tl0jcxev4cZzxjJl0qhQ2pzJWUNvAyNTno+IHxMRkbgrJo3i5Tv+jg8fOZCrHqhj1962Xv+MTCaChcBV8dlDZwDb3P3dDMYjIpIx9Q0xrnqgjvqGWJfn5y1Zy3PrtrBp+55e/+zQuobM7GfAFKDCzDYCdwKFAO5+L7AIOBdYD+wCPhdWLCIi/Vl9Q4wZDy8ltqsFgEeunUR1SpfQFZNGMXPaGACOGNTU658fWiJw9093c96BL4f1+SIi/V19Q4x5S9ayfXcLsV0tlJcWMn38cK56oI5lb8bY1dLGPYtWc8WkUUw8ppxHrp1EbW1tr8eRFYPFIiK5JPG//YoBRazbvJMJI4cwubKCwwcUccuTKwAojHfcHzmoOPR4lAhERPpQfUOMO37xGq3tTmubM7myItntc/GPXkheV5ifz5kfLE+eC5MSgYhIH0mMBSTWBdx87vFcMWlU8nhCcYFx23njuGLSqD6JS4lARCREqYO+z7z2bnIs4LKqkdyzaBU/+f0bDCjKTx6//+rTmHhMeZ/GqEQgIhKS1G6guYvXcP/VpwEwbvgg7n1uAwA7NjUlxwhmThvT50kAlAhEREJRXdeYTAIFecaN54xNzvw5ZdavktdVHlHG7eeNy0gCSFAiEBHpRYkpoS83xpJJ4NEvngnAVQ/UMXPaGG48Z2yHNQKZpkQgItJLUu8CKo8oo6CpmcuqRibXCizfuA0IFoz1hwSQoEQgInKIuroLmHPxSQDJFcOp4wD9jRKBiMghStQBStwFfGzMUGY8vJSKsuLkbKBMjwMciBKBiMghSvwvf+a0Mbz+3o7k6uDWtvaMzgZKlxKBiEgveHfbHq55sI629veP3Xxu3y0KOxRKBCIih6Bz5dDSwjzKSwv7zYygdCgRiIj0UOosodLCfPLzsucuIJUSgYhIGjrvDwAwd/Ga5Cyhn86Y1K/HAQ4kkzuUiYhkhUSpiNiuFuYuXpM8fuM5YykvLWTWheOzNgmA7ghERParviHG7KdW8qfNTR1KRSRcMWlU1nUDdUWJQEQkRedqoYnVwACHDyjiw0cOzGB04VDXkIhI3FdrXuaWJ1cku4DGDR9Evr1/ftOOZuYtWZu5AEOiRCAiQnAnsGD5O8nnl1WN5P7n36DNobggj3yDkeUl/bJExKFS15CIRF51XWNyNTBAfh48uuyt5LjAneefkBNjAfujRCAikZVaLC6htDCPo8tLWbepKWM7hvU1JQIRiaxEsbjCvGAgYOjAIu69sip5rr/XCOotGiMQkcgaN3wQAC3tDsDmHXuZt2RtciexKCQB0B2BiETYf9U1AGBASbxLKBcHg7ujRCAikVJd18g9i1Zz5KBiBpcWsaN5Nx86ooxnv/axTIeWMUoEIhIZ9Q0xbluwgnaHHZtbmTBiMMdVDIjkXUAqJQIRySmJmUCpA73BXcAq2tohPhxAHnD7+SdEZhzgQJQIRCRndN4b4JFrJzFn0WrufW5D8prUctFKAgElAhHJGfOWrE3uETxu+CBOmfUrtu1uSZ4vLczP6nLRYVEiEJGsltoVlOjrnz5+eHLDmOJ8o7XdGT6khHmXn6Ik0IVQE4GZTQfmAfnA/e4+p9P5UcDDwJD4NTe5+6IwYxKR3JJYFAZBV9C44YOS5SLyDe68YHxOl4foDaElAjPLB34InA1sBJaa2UJ3X5Vy2W3AY+7+IzMbBywCRocVk4jkjvqGGG9u2cn08ccBMHPaGOobYh3GA048erCSQBrCXFl8OrDe3Te4+16gBriw0zUODIo/Hgy8g4hIN+obYlzzYB07mlt5bGkjj1w7CYAZDy9NXlOcb9x+/gmZCjGrmLuH88ZmlwDT3X1G/PlngUnufkPKNcOBXwHlwABgmrvXd/Fe1wHXAQwbNmxiTU1Nj2JqamqirKysR6/NVmpzdESh3X/duZc/b99DQV4ee1rbGFYCsb15tLsnj+WbUVyYx/DBJZQW5Wc65F7X0+956tSp9e5e1dW5TA8Wfxp4yN2/Z2ZnAj81s/Hu3p56kbvPB+YDVFVV+ZQpU3r0YbW1tfT0tdlKbY6OXG93fUOMGfe9SGt7PgOLC/jg0CGcX7SN76/Io7mtneKCPCYdOzTnC8WF8T2H2TX0NjAy5fmI+LFU1wKPAbj7i8BhQEWIMYlIFkqsD0jsD/CZSaPAjKB3OZBvRKpQXG8KMxEsBSrN7FgzKwIuBxZ2uqYR+DiAmR1PkAg2hxiTiGSh1PUBj37xTH69ZhPL39rKrr1tjDq8hPLSQm47T+MBPRVa15C7t5rZDcBigqmhD7r7SjObBSxz94XA14Efm9k/EaT2azysQQsRyVozp41h++6W+F0AvLdtNxBUDZ1zycm6CzhEoY4RxNcELOp07I6Ux6uAj4YZg4hkt/qGGDc98SobNjfR5sHdwc3njmPu4jUcNaRQSaAXaGMaEenX5i1Zy7pNQRIoyDNmThvDFZNG8fIdf8fhA4oyHV5OyPSsIRER6htizH5qJZhx+3lBMbhExdAhpUWMKC9h2669KhQXEiUCEcm42U+vYvnGbcnHgw4r4OXGGDua29jRvJvJlRU8cu3fZjjK3KVEICIZkygYt3NPUCF0YHE+W3bsYflbexg6MOj2OXJwSeQ3jgmbEoGI9LlEAti+p5Xlb21lwsghTK6sYOa0MVzz4EsA7Nnbzoq7pmc40mjQYLGI9KnqukYuu+/FoGKoO5MrKzjj2MNZ8fY2Xn9vBzefezzlpYXcfO7xmQ41MnRHICJ9au7iNckVwomicJ+69wXaHO5ZtIoVd01XxdA+pjsCEQlVfUOMqx6oo74hBsCN54ylvLSQGWcdy+ynV3HNg3W0xZeRHjm4JIORRpcSgYiEKrFxzLwlawGSawBWvbud5W9tZUdzGwOL85kwcghzLj4pw9FGk7qGRCRUiRk/iY1jrv/PZWzesZd8g5HlJXxgQBG3n3+C1gdkUFqJwMxKgFHu/nrI8YhIlkvdQ3jiMeVMPKacmdPGMPvpVfxp0w52NLcB0ObQ1NzK776p9QGZ1m0iMLPzgX8BioBjzWwCMMvdLwg7OBHJHskpobtbkovDZk4bs8+xwjxoaQ9+v/GcsZkMWeLSuSP4FsG2k7UA7r7czI4NMSYRyUKJsYDEmoBxwwdx2X0v0truTBg5hAkjh4C7uoH6oXQSQYu7b7N4+dc4lYoWibiv1rzMguXvUJgPd11wYoexgInHlHPKrF+9P030PNUI6s/SmTW00syuAPLNrNLM/gN4IeS4RKSfW7D8HQBa2oK1Ac+ufI/fr9/CzJqXqW+IJaeJzrpwvJJAP5dOIvg/wAlAM/AzYDvw1TCDEpH+qbqukVNm/YrqukaOGFgMBJvD3HjOWO5//g3aHDbGdjNvydrkNFEtDuv/uk0E7r7L3W9199PcvSr+eE9fBCci/Ud1XSO3PLmC2K4W5i5ew4+unMiEkUM4ecRgPnzkQGacdSz5BiPKVSQu26Qza+g3dDEm4O6a8yUSIXMXr0k+Tsz2Wfvedna1tHPTE6/y7Nc+xk2qD5SV0hks/kbK48OAi4HWcMIRkf4ksTnMkYNLuKxqJI8ue4vLqkby2NJGVr+3nebW4P+IiT2EJTt1mwjcvb7Tod+b2UshxSMi/UR9Q4zbnlxBO7BjUxPDBx/G/VefxoyHlxLb1dLh2pvPHZeZIKVXpNM1dHjK0zxgIjA4tIhEpF+Yt2Qt7fHHBmzYsjNZJXRgcT5DSovYGt8+UgPC2S2drqF6gjECI+gSegO4NsygRCTzZk4bw7vb9tD4l520tDsbY0H3T0Ge8dDnJ2lKaA5Jp2tIq4hFImp3SxvN8RrReQYDivK1gXwO2m8iMLN/ONAL3f2/ez8cEekPqusauX3BiuQ+AQB3X3SiuoBy1IHuCM4/wDkHlAhEckx9Q4zZT69ixcatHZLA9ZOPUxLIYftNBO7+ub4MRET6Vudy0dV1jdzxi9dobQ8yQL7BcUPLmHPxSeoKynHp7kfwCYIyE4cljrn7rLCCEpHwzX5qJcs3buPlxhgPfX5Sci/hfIMTjx6sKqER0m2JCTO7F7iMoOaQAZcCx4Qcl4iEqLqukVfi+wPsaG5j9tOrkkXiZl90IgtuOEtJIELSuSP4iLufZGavuvtdZvY94JdhByYivSvRFTRu+CDufW5Dx5PuXDFplMYBIiqdRJBYO77LzI4C/gIMDy8kEelNiQSw+t3tbG7ay3PrtiTPFRfkMaq8hNvPPyGDEUqmpZMInjazIcBc4I8EM4Z+nM6bm9l0YB6QD9zv7nO6uOZTBLugOfCKu1+RXugi0p36hliXJSEAvv1JTQeVQDoLymbHHz5hZk8Dh7n7tu5eZ2b5wA+Bs4GNwFIzW+juq1KuqQRuBj7q7jEzO6InjRCRjqrrGpm7eA2F+Xn7JAED/llJQFKkU2voVaAGeNTd/0SwQU06TgfWu/uG+PvUABcCq1Ku+QLwQ3ePAbj7poOIXUQ6SawDeHXjVtq72FC2cugA5lxysgaCpQNzP/D2w2Z2DMGsocuAduBR4DF3b+zmdZcA0919Rvz5Z4FJ7n5DyjULgLXARwm6j77l7s908V7XAdcBDBs2bGJNTU3aDUzV1NREWVlZj16brdTm3LdrbxvvbtvD4IJWtrXms2tvW5fXHT2khMMHFPVxdOGK2ncNPW/z1KlT6929qqtz6XQNNQDfBb4b78q5HfgOwQ/uQ1UAVAJTgBHAc2Z2ortv7RTDfGA+QFVVlU+ZMqVHH1ZbW0tPX5ut1Obcd9UDdTy3bg9fPxEe21jKO1t3c3hZEUcPKeWMYw/n0WVvceM5Y/mHHOwKitp3DeG0Od0FZal3BW3A/03jZW8DI1Oej4gfS7URqHP3FuANM1tLkBiWphOXiARVQrfvaeWwghjvbN1Nm0Nrm7Pgyx8F0K5h0q10xgjqgELgMeDSRJ9/GpYClWZ2LEECuBzoPCNoAfBp4CdmVgGMAdJ9f5HISowF7GxuBXcGHFZIXh60eVAmOrGVpEg60rkjuMrdXwcwsyPTfWN3bzWzG4DFBN1ID7r7SjObBSxz94Xxc39nZqsI7jRudPe/HHQrRCJm3pK1LH+rQw8q5x9hTK6sSNYOEklXOmMEr6c8XQScmu6bu/ui+GtSj92R8tiBr8V/iUg36htizKx5mXe2dtwjeMLIIQwfvIdHLp2Uocgkm3Vba6gTCyUKEelWdV0jl933Ihtju2n34B9jQZ5x/eTjWPDlj1Ja1BvzNySK0hosTpHWimIR6R2J8hDTxw9Plog2gt3CvvA3x2kgWHpFurOGTgXOAtzMTnX3P4YblogAzH56Fcvf2pqsD1SQZ8y6cLxWBUuvSmfW0B0EpacTO5L9xMx+7u53hxqZSIRV1zVy11MraW5tTx4ryDMe/eKZGgiWXpfOHcFngJPdfQ+Amc0BlgNKBCIhqK5r5JYnV3Q4lm8w68LxSgISinQSwTsEO5PtiT8vZt+FYSJyCJLrAva0sGHLzuRxA44uL2He5acoCUho0kkE24CVZvYsQanos4GXzOzfAdz9KyHGJxIJM3/2RzZu3ZN8rrEA6UvpJIIn478SasMJRSR6qusauWfRanY0tyaPTRg5hNvPG6c7AOkz6Swoe7gvAhGJijmLVnP/829w5KDiDncBANdP1pRQ6XsHu45ARA5BdV1jcr/g1CSgfQIkk5QIRPpIV7OBACZXVvDItSoNIZmjRCDSR+56amWH57oLkP5iv4nAzJ4imCXUJXe/IJSIRHJIfUOM2U+tBDP2piwO08bx0p8c6I7gX/osCpEcVN8Q47P317GrJdg6ckR5Ce9t28OMs45VEpB+Zb+JwN1/25eBiOSSOYtWJweFIVgZrEVh0l+lU2uoErgHGEewwhgAdz8uxLhEslKiWujz67ckjw0szufmc7UuQPqvdAaLfwLcCfwrMBX4HAe/j4FIzqtviHHNg3XsaG6jMN9ob3NGlJfw/Df/NtOhiRxQOomgxN1/bWbm7g3At8ysHrijuxeK5LLE//5nThvDT198kwXL30meG/2BAQwffBgzp43JXIAiaUonETSbWR6wLr4H8dtAWbhhifR/85as5bl1W/j9+i20pcyvU4kIyTbpdPHMBEqBrwATgSuBq8MMSiQbjBs+CKBDErhowlEs+PJHlQQkq6RTa2hp/GETwfiASGSlbh15//NvJI9XHlHGnItPUgKQrJTOrKFngUvdfWv8eTlQ4+7nhB2cSH+Sui5g2Zt/pbXdVS5ackI6YwQViSQA4O4xMzsixJhE+p3qukZufXJFylJ7Y3JlBTOnjdFdgGS9dBJBu5mNcvdGADM7hgOUnhDJJfUNMa7/z2Vs3rG3w/HbzhunuwDJGekkgluB583stwQ75/0NcF2oUYlkWH1DjJsef4UNW3Z2GAyGYM8AJQHJJekMFj9jZqcCZ8QPfdXdtxzoNSLZqrqukbmL11AxoIh1m4O9g43gFrgw37jrAo0HSO45UPXRse6+Jp4EINjEHmBUvKvoj+GHJ9I3EruG5ecZza3ttLY5lUMH8N72Zm4+93j98JecdqA7gq8RdAF9r4tzDmjdvGS9xHTQxKKw9nanvLSQG88Zqx/+EhkHqj6aGAf4e3fvsLGqmR3WxUtEskpqhdChA4uI7WxhxlnHas9giZx0BotfAE5N45hIv1ffEGPmz/7IO9v20J4yCHz04BKW3np25gITyaADjREcCRwNlJjZKQRjZgCDCEpOiGSd2U+t7LBpPAQbxtx+/gkZikgk8w50R3AOcA0wAvh+yvHtwC3pvLmZTQfmAfnA/e4+Zz/XXQw8Dpzm7svSeW+RdNU3xJj99Creju1ic9P76wEMePxLH9GCMIm8A40RPAw8bGYXu/sTB/vGZpYP/BA4G9gILDWzhe6+qtN1AwkK29Ud7GeIdKe+Icbl971IS3vHxQDaLEbkfemMEcw1szOAB9199UG89+nAenffAGBmNcCFwKpO180GvgPceBDvLXJAifUAe1raOiSBfIMv/M1xGhAWSWHuB64WEf8f++W8vzPZgwRF57Z387pLgOnuPiP+/LPAJHe/IeWaU4Fb3f1iM6sFvtFV15CZXUd8NfOwYcMm1tTUpN/CFE1NTZSVRWsrhai2+c1tbfvUQRlSWsTI8pKMxNQXovpdq83pmTp1ar27V3V1Lp2VxTuAHwM/NrOPAdXAv5rZ48Bsd19/0BEB8c1uvk8wDtFdDPOB+QBVVVU+ZcqUnnwktbW19PS12SqKba5Z+Az/ssI6HJtcWcEjn5mUoYj6RhS/a7W5d6RThjof+ATBHcFoggVm/0VQc2gRsL+9+N4GRqY8HxE/ljAQGA/UmhnAkcBCM7tAA8bSU/UNMTbvaCbxVzvPYMM9n8hsUCL9XDpjBOuA3wBz3f2FlOOPm9nkA7xuKVBpZscSJIDLgSsSJ919G1CReH6griGR7ny15mUWLH+HoQOLuGp0cKy4II87NS1UpFvpJIKT3L2pqxPu/pX9vcjdW+N7HC8mmD76oLuvNLNZwDJ3X9ijiEXiEj/8E0XhADbv2MvA4gKe0LRQkbQdaEHZfxD/9xXvuungQEkg5ZpFBN1Hqcfu2M+1U7p7P5FUC5YHdRBTB4UvmnAUoyu2KQmIHIQD3RGoi0b6neq6Ru5ZtJojBxVTkAet7cHx6ye/PyW0trY2cwGKZKHuFpSJZFx9Q4zZT61k5942Nmxuos1hx+ZWKo8oY0tTsyqFihyidGYNDQW+CYwDklVH3V1lqCV09Q0xrpj/Is1tHReFHVcxgDkXn6QuIJFekM5g8X8BjxJMIb0euBrYHGZQIomB4MJ8oyUlCUwYMZjbzz9BCUCkF6WTCD7g7g+Y2Ux3/y3wWzNbGnZgEm2JgeCWNqc4P4/W9naVhhAJSTqJoCX++7tm9gmCLSsPDy8kiaL6hhgza17m3a27Of/ko5J3AoX5RvV1Z+gOQCRE6SSCu81sMPB14D8I9iP4p1CjkkhJ3SkM3r8bKC8t5P6rT1MSEAlZOrWGno4/3AZMDTcciZqrHqjjuXVb9jmuJCDSd9KdNfQFgjpDyevd/fPhhSVR0TkJTK4Mqo7MnDZGSUCkj6TTNfQL4HfAEqAt3HAkKi76wfMs37itw7Fvf/JErQcQyYB0EkGpu38z9EgkEuobYtz0xKus29SxfJWSgEjmpJMInjazc+N1g0R6rPOgcIKSgEhmpZMIZgK3mFkzwVRSA9zdB4UameSMOYtWc99zG/bZMazyiDKtDhbpB9KZNTSwLwKR3LO/GUEQVAn9t8tP6eOIRKQr6dwRYGblQCUdaw09F1ZQkt3mLFrN/Oc20N7FuSMGFvGjK6t0FyDSj6QzfXQGQffQCGA5cAbwIqCic5IUlIdexZGDS9iwualDEhhSWsADV5+uH/4i/VS6YwSnAX9w96lmNhb4drhhSbaob4hx0+OvsG7zTgB2bGqiON9ob3OKCow7zx+vgWCRfi6dRLDH3feYGWZW7O5rzOzDoUcm/Vp9Q4x5S9byxpadvBXbnTyeZ9Dc5kyurOCRaydlMEIRSVc6iWCjmQ0BFgDPmlkMaAg3LOnPqusauX3BCtoc8lKODx1YxD9N+zDPvPYuM6eNyVh8InJw0pk19Mn4w2+Z2W+AwcAzoUYl/VJ1XSN3LXytwyYxR5WXsDF+R9Da5lwxaZS6gkSyTFqzhhLi+xFIxNQ3xPhqzcsduoAASgvz+ccpHwJg7uI13HjO2EyEJyKH6KASgUTP2d+rTQ4EpxoRvxN45rV3eeTaSboLEMliSgTSpeq6Ru5+ehW7WjrWGRxRXsK8+EKweUvWaixAJAcoEcg+EvsFpzLgnzvVBNKsIJHcoEQgHcxZtHqfJDCivISKsmI+fKSqjYjkorzuL5GoqG+I7VMdtDAPNsZ2s/ytrcxbsjZDkYlImHRHEHFBaYjVgLOj+f3xgOL8PEZ9oJR1m5oYWJzPB48YqPEAkRylRBBRidpAqT/8E5740keYeEx5cvWwto0UyW1KBBHUXXnoxA/9iceUa0BYJAKUCCIgMRW0rd0ZdXhJh3UBhfnQ0gaVQwfw7NenZC5IEcmYUAeLzWy6mb1uZuvN7KYuzn/NzFaZ2atm9mszOybMeKKmuq6Rle9s55YnV7CrpY3mtnbWbd7JkNIg/0+urODM4yoAGD6kJJOhikgGhXZHYGb5wA+Bs4GNwFIzW+juq1IuexmocvddZvYl4LvAZWHFFCWJ/YG/fuL7dYGK8/MYdXgJcy45Odn9U98QA9BAsEiEhdk1dDqw3t03AJhZDXAhkEwE7v6blOv/AFwZYjyR0dUm8ddPPo6bzj1+n2s1DiAi5t55S/FeemOzS4Dp7j4j/vyzwCR3v2E/1/8AeM/d7+7i3HXAdQDDhg2bWFNT06OYmpqaKCsr69Frs8V72/eweUdz8vmwEigoLuHwAUUZjKpvReF77koU2602p2/q1Kn17l7V1bl+MVhsZlcCVcDHujrv7vOB+QBVVVU+ZcqUHn1ObW0tPX1tNqiua+SWZ1aQ+rV+qwrm1cON5xwXmcJwuf49708U2602944wE8HbwMiU5yPixzows2nArcDH3L2583npXn1DjC/95zI27djb4fhFE45id8tbxHYFZaKjkghE5OCEOWtoKVBpZseaWRFwObAw9QIzOwW4D7jA3TeFGEvOqq5r5NIfvdAhCRhBfaDPnjmaYYMOo7y0UHsFiMh+hXZH4O6tZnYDsBjIBx5095VmNgtY5u4LgblAGfBzMwNodPcLwoop11z0g+dZvnFbh2MTRgzmtXe2szG2m3lL1vL5Dxbx8h1TMhOgiGSFUMcI3H0RsKjTsTtSHk8L8/Nz2VdrXu6QBEoL87ntvHHMXbyG1nanIM+YOW0MO954JYNRikg26BeDxZK+OYtWM/+5DbR3Or5q9vTk48S2kROPKaf2jb6NT0SyjxJBFulqfQAEawQStHm8iBwsJYIs0dWuYftbJCYicjCUCPq5/U0NBVj17vYMRCQiuUaJoB+rrmvkjl+8Rmt7x9XfF004ir/u3Kv6QCLSK5QI+qH6hhjXPvwSW3e1AsG6ACdIAP92+SkZjU1Eco8SQT9T3xDjU/e9QFvKtKB//uSJGgAWkdAoEfQjXe0cNrmyQklAREKlRNAP1DfEmPmzP7Jx654OxydXVqhEtIiETokgw/a3NkBJQET6ihJBBnXVFTRiyGE8f9PHMxSRiESREkGGnDXn1/t0BQ0szmfep0/NUEQiElVKBH2sq7sAgG9rZpCIZIgSQR+prmtk7uI1xHa1dDhelG/87Lozk5vJi4j0NSWCPtDVvuVXD9QAAAnXSURBVAGgBWIi0j8oEYSouq6RuxaupLmtc9HoYFaQkoCI9AdKBCHZ37TQCSMGs+CGszIQkYhI15QIell9Q4zrf7qMzU37VgtVEhCR/kiJoBed/b1a1m3euc9xJQAR6c+UCHrB/rqBNBgsItlAieAQ7G9NACgJiEj2UCLogf3dAYASgIhkHyWCg7S/uwDtHywi2UqJIE3VdY3c+uQKvNPxgcX5PPT5SVoZLCJZS4mgGwcaB1ClUBHJBUoE+1HfEONT975AW+dbAKBy6ACe/fqUPo9JRCQMSgRdmHDXYrbubu3ynMYCRCTXKBF0cvxtv2R3a8faQPkGj13/EY0DiEhOUiKIq2+IcemPXqBzebihZUUsve3sjMQkItIXlAiAD93yP7TuWyBUawJEJBIinQj2NyOorCif12ZNz0BEIiJ9L9REYGbTgXlAPnC/u8/pdL4YeASYCPwFuMzd3wwzJoDT7n62y+qgoBlBIhI9oSUCM8sHfgicDWwElprZQndflXLZtUDM3T9kZpcD3wEuCyumN7fsZPRN/9PlOY0FiEhU5YX43qcD6919g7vvBWqACztdcyHwcPzx48DHzczCCGb8Hc+wo7nrKaFPfOkjSgIiElnm3sWKqd54Y7NLgOnuPiP+/LPAJHe/IeWa1+LXbIw//1P8mi2d3us64DqAYcOGTaypqTnoeFa8vY1hJfDn3R2Pn3j04IN+r2zS1NREWVlZpsPoU1FsM0Sz3Wpz+qZOnVrv7lVdncuKwWJ3nw/MB6iqqvIpU6Yc9HvccMczfPHDzXxvRUGkSkPU1tbSkz+vbBbFNkM02602944wE8HbwMiU5yPix7q6ZqOZFQCDCQaNe91rs6ZTW1vLm5+ZEsbbi4hkrTDHCJYClWZ2rJkVAZcDCztdsxC4Ov74EuB/Pay+KhER6VJodwTu3mpmNwCLCaaPPujuK81sFrDM3RcCDwA/NbP1wF8JkoWIiPShUMcI3H0RsKjTsTtSHu8BLg0zBhERObAwu4ZERCQLKBGIiEScEoGISMQpEYiIRFxoK4vDYmabgYYevrwC6HoD4tylNkdHFNutNqfvGHcf2tWJrEsEh8LMlu1viXWuUpujI4rtVpt7h7qGREQiTolARCTiopYI5mc6gAxQm6Mjiu1Wm3tBpMYIRERkX1G7IxARkU6UCEREIi4nE4GZTTez181svZnd1MX5YjN7NH6+zsxG932UvSuNNk82sz+aWWt897isl0abv2Zmq8zsVTP7tZkdk4k4e1Mabb7ezFaY2XIze97MxmUizt7UXZtTrrvYzNzMsn46aRrf8zVmtjn+PS83sxmH9IHunlO/CEpe/wk4DigCXgHGdbrmH4F7448vBx7NdNx90ObRwEnAI8AlmY65j9o8FSiNP/5SRL7nQSmPLwCeyXTcYbc5ft1A4DngD0BVpuPug+/5GuAHvfWZuXhHcDqw3t03uPteoAa4sNM1FwIPxx8/DnzczKwPY+xt3bbZ3d9091eB9kwEGIJ02vwbd98Vf/oHgl3yslk6bd6e8nQAkO2zQdL59wwwG/gOsKcvgwtJum3uNbmYCI4G3kp5vjF+rMtr3L0V2AZ8oE+iC0c6bc41B9vma4FfhhpR+NJqs5l92cz+BHwX+EofxRaWbttsZqcCI939f/oysBCl+3f74ni35+NmNrKL82nLxUQg0oGZXQlUAXMzHUtfcPcfuvsHgW8Ct2U6njCZWR7wfeDrmY6ljz0FjHb3k4Bneb+Ho0dyMRG8DaRmxxHxY11eY2YFwGDgL30SXTjSaXOuSavNZjYNuBW4wN2b+yi2sBzs91wDXBRqROHrrs0DgfFArZm9CZwBLMzyAeNuv2d3/0vK3+f7gYmH8oG5mAiWApVmdqyZFREMBi/sdM1C4Or440uA//X4CEyWSqfNuabbNpvZKcB9BElgUwZi7G3ptLky5ekngHV9GF8YDthmd9/m7hXuPtrdRxOMBV3g7ssyE26vSOd7Hp7y9AJg9SF9YqZHyEMadT8XWEsw8n5r/Ngsgr8gAIcBPwfWAy8Bx2U65j5o82kEfY07Ce5+VmY65j5o8xLgz8Dy+K+FmY65D9o8D1gZb+9vgBMyHXPYbe50bS1ZPmsoze/5nvj3/Er8ex57KJ+nEhMiIhGXi11DIiJyEJQIREQiTolARCTilAhERCJOiUBEJOKUCCSyzOyFnpzr5j2/ZWbf6HlUB/15b5pZRV99nuQmJQKJLHf/SOdj8ZXmXZ4TyVVKBNLvmdmVZvZSvO76fWaWHz/eZGZzzWylmS0xs9PNrNbMNpjZBfFrrjGzX8SPrzOzO1Petyn++xQz+52ZLQRWpZ6LP/5mvMb/K2Y2J37sC2a2NH7sCTMr7aYND5nZvWa2zMzWmtl5XVxzvZnNTXl+jZn9IP54gZnVx9t6XRevHW1mr6U8/4aZfSv++INm9kz89b8zs7Hp/LlLdCgRSL9mZscDlwEfdfcJQBvwmfjpAQTlQU4AdgB3A2cDnyRYhZlwOnAxwX4Ml+6nDs2pwEx3H9Pp8/+eoATwJHc/maCiJ8B/u/tp8WOrCaqbdmd0PJZPAPea2WGdzj8Rjz3hMoJ6QQCfd/eJBMXzvmJmB1Mtdz7wf+Kv/wbw/w7itRIBBZkOQKQbHycoqLU0vmVECZCoG7QXeCb+eAXQ7O4tZraC4IduwrPu/hcAM/tv4Cygcy2al9z9jS4+fxrwE4/va+Duf40fH29mdwNDgDJgcRpteczd24F1ZrYBGEtQCoL4e2+O382cQVAjaCzw+/jpr5hZIkmMBCpJo1CimZUBHwF+nrLlRnEasUqEKBFIf2fAw+5+cxfnWvz9GintQDOAu7cn+vrjOtdR6aquys6DjOsh4CJ3f8XMrgGmpPGafeIws58ApwDvuPu5BHcAnwLWAE+6u5vZFIKEdKa77zKzWoJ6Wala6XiHnzifB2yN302JdEldQ9Lf/Rq4xMyOADCzw+3g9x4+O/66EoKyzL/v7gUpngU+lxgDMLPD48cHAu+aWSHvd1V151IzyzOzDxJsQ/i6u3/O3SfEkwDAkwRdUZ/m/W6hwUAsngTGEpRa7uzPwBFm9gEzKwbOg+SOZW+Y2aXx+M3MTk6/+RIFSgTSr7n7KoLNVX5lZq8S/GAefuBX7eMlgv73V4En/CBKFLv7MwQlgJeZ2XKCPnaA24E6gqSyJs23a4zH8kvgenffZ1tFd48RjDkc4+4vxQ8/AxSY2WpgDkGp5c6vayEYF3mJ4M8oNabPANea2SsEFStD3fZQso+qj0pOi3fbVLn7DRmO4yHgaXd/PJNxiHRFdwQiIhGnOwIRkYjTHYGISMQpEYiIRJwSgYhIxCkRiIhEnBKBiEjE/X+oJ6DP3hpwdAAAAABJRU5ErkJggg==\n", "text/plain": ["<Figure size 432x288 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["test_simulations(n=6000, mu=mu, sigma=sigma)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["\n", "\n", "\n", "For more info:\n", "\n", "- Sampling distribution of the sample mean [Khan Academy](https://www.khanacademy.org/math/ap-statistics/sampling-distribution-ap/sampling-distribution-mean/v/sampling-distribution-of-the-sample-mean)\n", "\n", "- Central limit theorem [Wikipedia](https://en.wikipedia.org/wiki/Central_limit_theorem)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Refs\n", "\n", "- [scipy.stats.ttest_ind](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html)\n", "- [Wikipedia T-test](https://en.wikipedia.org/wiki/T-test#Independent_two-sample_t-test)"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": []}], "metadata": {"kernelspec": {"display_name": "p38", "language": "python", "name": "p38"}, "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.3"}}, "nbformat": 4, "nbformat_minor": 4}