{
"cells": [
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from jupyterthemes import get_themes\n",
"from jupyterthemes.stylefx import set_nb_theme\n",
"themes = get_themes()\n",
"set_nb_theme(themes[1])"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ethen 2017-04-23 14:48:05 \n",
"\n",
"CPython 3.5.2\n",
"IPython 5.3.0\n",
"\n",
"numpy 1.12.1\n",
"pandas 0.19.2\n",
"matplotlib 2.0.0\n",
"tqdm 4.11.2\n"
]
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# 1. magic for inline plot\n",
"# 2. magic to print version\n",
"# 3. magic so that the notebook will reload external python modules\n",
"%matplotlib inline\n",
"%load_ext watermark\n",
"%load_ext autoreload \n",
"%autoreload 2\n",
"\n",
"from tqdm import trange\n",
"from scipy.stats import bernoulli, ttest_ind\n",
"\n",
"%watermark -a 'Ethen' -d -t -v -p numpy,pandas,matplotlib,tqdm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The multiple hypothesis testing problem\n",
"\n",
"This is simply walking through the content at following link. [Blog: The multiple hypothesis testing problem](http://www.claudiobellei.com/2016/10/30/multiple-testing-part1/).\n",
"\n",
"Suppose a particular event has a chance of 1% of happening at the first attempt. Now, if we make $N$ attempts what is the probability that this event will have happened at least once among the $N$ attempts? The answer is $1 − 0.99^N$. For $N = 5$, the rare event already has an almost 5% chance of happening at least once in the five attempts. This is the rationale behind why multiple testings can complicate things in statistical inference.\n",
"\n",
"Suppose we formulate a hypothesis test by defining a null hypothesis $H_0$ and alternative hypothesis $H_1$. We then set a type-I error level $\\alpha$, which means that if the null hypothesis $H_0$ were true, we would incorrectly reject the null with probability $\\alpha$. Then given $N$ tests the probability of rejecting the null in any of the tests can be written as\n",
"\n",
"\\begin{align}\n",
"P(\\mathrm{rejecting\\ the\\ null\\ in \\ any \\ of \\ the \\ tests}) = P(r_1 \\lor r_2 \\lor \\dots \\lor r_n)\n",
"\\end{align}\n",
"\n",
"In which $r_j$ denotes the event \"the null is rejected at the j-th test\". To compute the probability above we can convert it to:\n",
"\n",
"\\begin{align}\n",
"P(r_1\\lor r_2\\lor\\dots\\lor r_n) = 1 - P(r_1^* \\land r_2^* \\land\\dots\\land r_n^* ) = 1 - \\prod_{j=1}^N P(r_j^*)\n",
"\\end{align}\n",
"\n",
"Where $r_j^*$ denotes the event \"the null is NOT rejected at the j-th test\".\n",
"\n",
"Suppose that we do a test where we fix the type-I error $\\alpha$ to 0.05 (0.05 is simply the conventionally used number). By definition, if we do one test only we will reject the null 5% of the times if the null is actually true. What if we make 2 tests? What are the chances of committing a type-I error then? The error will be:\n",
"\n",
"\\begin{align}\n",
"P(r_1 \\lor r_2) = 1 - P(r_1^* \\land r_2^*) = 1-P(r_1^*)P(r_2^*) = 1-0.95 \\times 0.95 = 0.0975\n",
"\\end{align}\n",
"\n",
"So the consequence of doing $n$ multiple test is that the effective type-I error will become higher than expected, or to be more explicit, it will become:\n",
"\n",
"$$\\mathrm{Type \\ I \\ error} = 1-(1-\\alpha)^N$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To prevent this from happening, we can:\n",
"\n",
"- In many cases, it is not even necessary to do multiple tests\n",
"- If multiple testing is unavoidable (for example, we are testing multiple hypothesis in a single test because we have multiple groups), then we can just correct the type-I error as an effective type-I error $α_{effective}$. In order to recover the original type-I error (thereby “controlling” the multiple testing), we must ensure that:\n",
"\n",
"\\begin{align}\n",
"1-(1-\\alpha_\\mathrm{effective})^N &= \\alpha \\ \\Rightarrow \\ \\alpha_\\mathrm{effective} = 1 - (1-\\alpha)^{1/N} \\approx \\frac{\\alpha}{N}\n",
"\\end{align}\n",
"\n",
"For the last part, we're assuming $\\alpha\\ll 1$ and using taylor expansion we can obtain $(1-x)^m \\approx 1-mx$. [StackExchange](http://math.stackexchange.com/questions/210110/whats-the-name-of-the-approximation-1xn-approx-1-xn). And this goes under the name of **Bonferroni correction**."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example\n",
"\n",
"We will draw samples from two Bernoulli distributions A and B, each with a probability $p=0.5$ of success. Each hypothesis test looks like:\n",
"\n",
"\\begin{eqnarray}\n",
"H_0 &:& \\mu_B -\\mu_A = 0 \\nonumber \\\\\n",
"H_1 &:& \\mu_B -\\mu_A \\neq 0 \\nonumber\n",
"\\end{eqnarray}\n",
"\n",
"Where $\\mu_A$ and $\\mu_B$ are the two sample means. By our definition, the null $H_0$ is true as we are going to set both $\\mu_A$ and $\\mu_B$ to 0.5. The section below computes the probability of committing a type-I error as a function of the number of independent t-tests, assuming $\\alpha=0.05$. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# set parameters\n",
"alpha = 0.05 # Type I error rate\n",
"p1 = 0.5 # Probability for population 1\n",
"p2 = 0.5 # Probability for population 2\n",
"\n",
"# simulation parameters\n",
"n_samples = 500 # number of samples each test will use\n",
"n_tests = 20 # max number of tests\n",
"n_sims = 300 # number of simulations from which to draw average\n",
"\n",
"def run_exp(n_tests, n_samples, p1, p2, alpha):\n",
" \"\"\"\n",
" simulations without Bonferroni correction, if it\n",
" committed a Type I error in any one of the test\n",
" return 1, else 0\n",
" \"\"\"\n",
" for i in range(n_tests):\n",
" testA = bernoulli(p1).rvs(n_samples)\n",
" testB = bernoulli(p2).rvs(n_samples)\n",
" _, pvalue = ttest_ind(testA, testB)\n",
" if pvalue < alpha:\n",
" return 1\n",
" \n",
" return 0\n",
"\n",
"def run_exp_corrected(n_tests, n_samples, p1, p2, alpha):\n",
" \"\"\"simulations with Bonferroni correction\"\"\"\n",
" for i in range(n_tests):\n",
" testA = bernoulli(p1).rvs(n_samples)\n",
" testB = bernoulli(p2).rvs(n_samples)\n",
" _, pvalue = ttest_ind(testA, testB)\n",
" if pvalue < alpha / n_tests:\n",
" return 1\n",
" \n",
" return 0"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 20/20 [02:56<00:00, 14.28s/it]\n"
]
}
],
"source": [
"p_rejects = []\n",
"p_rejects_corrected = []\n",
"for n_test in trange(1, n_tests + 1):\n",
" # results without using Bonferroni correction\n",
" p_reject = [ run_exp(n_test, n_samples, p1, p2, alpha) \n",
" for _ in range(n_sims) ]\n",
" p_reject_avg = np.mean(p_reject)\n",
" p_rejects.append(p_reject_avg)\n",
" \n",
" # results using Bonferroni correction\n",
" p_reject_corrected = [ run_exp_corrected(n_test, n_samples, p1, p2, alpha) \n",
" for _ in range(n_sims) ]\n",
" p_reject_corrected_avg = np.mean(p_reject_corrected)\n",
" p_rejects_corrected.append(p_reject_corrected_avg)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAAF9CAYAAACAkIXdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl01NX9//HnnSQwZCEbAyQhCRJ2kUUREKsg9ue+W6tY\n0dpaK6nVbqHa2lKtXSRdtCq0KG60im1xKaK1/aqgqJGCsoNoEEJIWBJIwoQM2e7vjwlD9kxgkpkk\nr8c5nMznzp2ZN+Qc8s69n/t+G2stIiIiItL1OIIdgIiIiIicGCVyIiIiIl2UEjkRERGRLkqJnIiI\niEgXpUROREREpItSIiciIiLSRSmRExEREemilMiJiIiIdFFK5ERERES6KCVyIiIiIl1UeLAD6GjG\nmMuBy2NiYr41fPjwYIcjIiIi0qa1a9cWWWtdbc0zPaXX6sSJE+2aNWuCHYaIiIhIm4wxa621E9ua\np61VERERkS5KiZyIiIhIF6VETkRERKSLUiInIiIi0kV1+1OrbamqqiI/Px+PxxPsUKQdwsLCiIuL\no1+/fjgc+n1ERER6ph6fyOXn5xMTE8PgwYMxxgQ7HPGDtZaqqir27dtHfn4+aWlpwQ5JREQkKHr8\nUobH4yExMVFJXBdijKFXr16kpKRQXl4e7HBERESCJmQSOWNMgjHmZWNMuTFmlzHmxhbm/dkY4673\n56gx5vBJfvbJvFyCRFuqIiLS04XS1urjQCUwABgPLDfGrLfWbq4/yVp7B3DHsWtjzDNAbSfGKSIi\nIhISQmJJwxgTBVwL/Mxa67bWrgJeBWb5+bpnOz5KERERkdASEokcMByottZurze2Hji1jdddCxwA\n3m3uSWPM7caYz4wxB/Ly8gITaR23283cuXNxuVw4HA5cLhdz587F7XYH9HOOqa6u7pD3FRERka4r\nVBK5aKCs0VgZENPG624BnrMtNIy11i601g6z1roCebLR7XYzZcoU5s2bR1FREdZaioqKmDdvHlOm\nTAlYMjd48GAeeughxo4dS1RUFA8++CAZGRnExMQwevRoXn75Zd/c9PR01q5dC8Df/vY3jDFs3uzd\nlV60aBFXXXVVQGISERGR0BEqiZwb6NtoLBZo8RCDMSYNmA4813FhNS87O5vc3Nwmtec8Hg+5ublk\nZ2cH7LNeeOEFli9fTklJCSNGjOC9996jtLSUuXPnctNNN1FYWAjAtGnTWLFiBQArV65kyJAhvPvu\nu77radOmBSwmERGRnqSzd+HaI1QSue1AuDFmWL2xccDmFuaD9/659621Ozo0smbMnz+/xQLCHo+H\nBQsWBOyz7rrrLlJTU+nTpw/XXXcdycnJOBwOrr/+eoYNG8bq1asBbyK3cuVKAN577z3uvfde37US\nORER6SyhnPSciM7ahTtRIZHIWWvLgZeAB4wxUcaYLwFXAItbednNwDOdEF4TxcXFJ/V8e6Smpvoe\nP/fcc4wfP564uDji4uLYtGkTRUVFgDeRe++99ygsLKSmpoavfvWrvP/+++zcuZPS0lLGjx8fsJhE\nRESa09lJT2ckjZ25C3ciQiKRq5MJ9AH2A88Ds621m40xaXX14nw3uRljzgIGAf8IRqCJiYkn9Xx7\nHKtxt2vXLr71rW/x2GOPUVxcTElJCWPGjOHY7YFDhw4lMjKSRx99lHPPPZe+ffsycOBAFi5cyJe+\n9CXVXBMR6cE6a5WsM5OezkoaO3MX7kSEzE93a+1Ba+1V1tooa22atfb5uvE8a220tTav3twP6+ad\nVCHgE5WZmYnT6Wz2OafTyezZswP+meXl5RhjcLlcADz99NNs2rSpwZxp06bx2GOP+bZRp0+f3uBa\nRER6ns5cJevMpKezksbO3IU7ESGTyHUlWVlZZGRkNEnmnE4nGRkZZGVlBfwzR48ezQ9/+EPOOuss\nBgwYwMaNGzn77LMbzJk2bRqHDx/m3HPPbfZaRER6ns5cJevMpKcjk0ZrLaUVVWzfd5j+Y88l6rT/\nR+zUG0i48Du4rv05/a661zc3kLtwJ8K0ULmj25k4caJds2ZNk/GtW7cyatSodr+f2+0mOzubBQsW\nUFxcTGJiIrNnzyYrK4vo6OhAhCx+ONHvn4hIT+FyuXz3U7f0/P79+7vcZzkcDlrLYRwOBzU1NU3G\na2otxe6j7C3zsLfUc/xr/cdlHo5UNn3tMbUeN7sfuQGn08mcOXO4//77A/J3qs8Ys9ZaO7GteaHU\noqtLiY6O5v777++Qb56IiEigdOYqWWZmJvPmzWt2pSzQtx4lJiY2TRrDwgmPTiQsph/xyYNZ+G4u\ne0uPsreswpes7T98lOrak1vEcjij6RMTx5C0lA7ZhWsPJXIiIiLdWLMJT6PnAyUrK4ulS5c22coN\nxK1HNbWW/Yc9FJRUsKfEw7m3P8B7azdBVALhfV2ExfQjLDK2wWt+/fq2E/48Z4SDpNg+DOjbm36R\n4ezcuo6PV71N2b48oh3V/PB7d/PjOT8K+i5cyCRyxpgEYBFwAVAE3HvswEMzc4cAfwKmAUeBp6y1\nczorVhERka6iM1fJoqOjycnJafetR9ZayjzVFJRUUFjqTdQKSirq/fFud9Y0WElLI3L8iXVtio+M\nYEBfJwNjnSTFOhnQt/7XPgzs66Rvn3Bf5QivM4FvndDndaSQuUfOGPMC3sMX3wTGA8uBqdbazY3m\n9QK2Ao8DfwFqgOHW2g2tvX+g75GT0KDvn4hI646dWm1plSwnJ6fDV5Uqq2vZV+ZhT73kbE+Jh8LS\n44ma++jJ9xS3tTVEOaoZmtKPlPioFpM1Z0RYAP5WHatL3SNnjIkCrgXGWGvdwCpjzKt4uzfc02j6\n14ECa+0f6o21msSJiIj0VCe6StYenqoa8g9VkH/oCPmHKthd93XPIW+idsB9lECsG/WL7kVSbB+S\n45wkx/UhJa4PyXF9SIr1XveL7k2Yw7T9Rt1ISCRywHCg2lq7vd7Yery9VBubAuw0xryBd51zE/Bd\na+3GDo9SREQkgI5VQJg/f74vwcrMzAx4BYSTPaBXWV1LQUlFvSTtCLsPehO33YcqOHD46EnH6Ixw\nHE/OYusStDhng2StK6ykdbZQSeSigbJGY2VATDNzBwHn4W3h9RZwN/CqMWaktbay/kRjzO1AFhB3\nrJCuiIhIKGhuy/NYod6lS5d2ypbnMdU1tRSWenwrafkHG66s7S3znNSKmjEwIMZJUv2VtLpVtGN/\n4iMjGt2TJv4IlUTODfRtNBYLNNe5oQJYZa19A8AY8zvgPmAU3lU8H2vtQmAheO+RC3DMIiIiJ8yf\nQr2BKnFlreWA+yh5xUfYVXzEtw16LFErLG18kKB9HAaSYvuQmtCHQfGRpMZHMii+Dynx3qRtQF8n\nvcLVg6AjhEoitx0IN8YMs9Z+Vjc2DtjczNwNwNnNjIuIiHQZ/nQmaE8iV1NrKSipIO+gN1nbVVzO\nruIj7CwuJ+/gkVYL3LbFGBjY18mg+D6+JG1QQqTvemCsk4gwJWrBEBKJnLW23BjzEvCAMeY2YALe\nrdOpzUz/K/BDY8yXgXeAu/CWK9naWfGKiIicrBMp1Hu02nuo4FiSVj9h233oCFU1J76q5orp3SBR\nS62XqCXFOekdrvvTQlFIJHJ1MoGngP1AMTDbWrvZGJMGbAFGW2vzrLWfGmNuAv4M9Ac+Bq5ofH+c\ntO3ee+9lwIABfO973wtqHJMmTeLpp5/m1FNPbTA2efJk3n77be644w6++93vBjFCEZHAa6lQr4lw\nEh6XRGL6CP68MrdB0lZQWnHC96rF9A4nvV8kaQmRdUlaJKnx3q3QQfF9dJCgiwqZOnIdravWkXvs\nscd45pln2LhxIzNnzuSZZ54JyPseOHCA8ePH8/nnn9OnT58O/ay2/P3vf+fFF19k6dKlgPdm38mT\nJ7Nu3ToqKyuZNWsWr7/+erOvDfXvn4hIczxVNcx5IJvnXnkTG9OfiPgUwuOTiIhLIiw6/oTft190\nL9ISIhmcGEVaYsOvOkzQtXSpOnLSsuTkZO677z7efPNNKioqAva+zzzzDJdccokvievIz2rLFVdc\nwR133MHevXsZOHAgGzZs4IYbbiAmJoYtW7YwePDgTotFRCRQJUFqay17yzzsOFDOjiI3Ow6Uk3vA\nzRdF5ewpqcDa04i79LR2xWYMJMf28SZr/SJJS4hicGIkaYmRpCdGEd1bP9Z7Gn3HQ9w111wDwJo1\na8jPzw/Y+77xxht84xvf6JTPqq6u5je/+Q2LFi3i8OHDPProo+Tn51NVVcVPf/pTnE4nZ5xxBm++\n+Sa33HILGzZs4PTTTwdg/fr1jBs3LmCxiIi05kRKgriPVrPjgLsuYSv3Pf6iqJyKqvYfMHBQ603Q\nXNGkJ3gTtPS6RE1boNKYErkeauPGjYwYMaJTPuu+++5jzZo1rF+/nnfffZc5c+bgcDjIycnxzRk1\nahTr16/3xXb55ZcD3kTu6quv7pQ4RURaLAlytJKdRW7u/u1CzrrwqgYJ2/4TKIbrMDAoPpIhrihO\n6RfFEFc0p9QlbMlxfXpcdwI5cUrkGhl8z/JO+6ydv7200z6rsZKSEmJimqu3HFhlZWU8/PDDbNmy\nhdjYWCZPnsy2bdv41a9+1eDzY2JiKCwsBGDRokW+8d/+9rcdHqOIhL7O6oAw/y9PUNM3maghaUT0\nSycicRDhCSlExCVhwiN4qxreWu5/kYS4yAiG1CVqQ1xRvsfpiZE6BSoBoUSuG5g+fTorV65s9rmz\nzz6bVatWNRmPj4/n8OHm6i0H9vPefvtthg8fzpAhQwCorKwkNja2ySnUw4cPExcXd0LxiEj31hEd\nEDxVNXy+381n+w+zfZ+bz/Z5H0fe8heiTPvqoUWEGdITo5pN2BKierXrvUTaS4lcN7BixYp2v2bs\n2LFs376dM888s0M/r6CggOTkZN/1woULSUlJabIauHXrVm666aZ2xyIiwdUZK2Un0wHhWML2+X43\n2/fVJW37D5N38EizZTxMK0lc9eFijHs/s666sC5Ri2JIv2gGxfchXMVwJUiUyDUSzO3O5lRXV1Nd\nXU1NTQ01NTV4PB7Cw8MJDz+5b90ll1zCypUr+drXvtahnzVo0CDWrVtHYWEheXl5LF68GLfbTWVl\nJb16eX9T9Xg8rF27lmefffak/k4i0rk6q1eoPx0Q7r3v5+Qe8K6sbd93mM/2u/lsnzdha0/nKYOl\n6lABlQfyqCraRVVxPlUH86k6uIfeDsucOXO4/+r2nTQV6UhK5ELcgw8+2OA3zb/+9a/MnTuXX/zi\nFyf1vjfffDPjx4+noqLCV4KkIz7roosu4oILLmDUqFEkJCTw0ksvMWfOHGbMmOHbgl22bBnTp09v\nsHInIqGvs3qFNuhwEBZBROIgevVLJyIxlYh+aUS40hn983+3K2FzGEhPjGJY/2iGDYhm+IAYhvWP\nYUAkTPvS1CZ/L6fTSUZGBllZWSf99xEJpJApCGyMSQAWARfgbbl1r7X2+Wbmfb1uXv1CZ5dZa1e0\n9v5dtSBwR/rJT35C//79g97ZYfLkySxatIgxY8a0+7U9+fsn0pLOOhjgcrma7UxQ//n9+/ef8Psf\nOHyUrYVl3HDHD6mM6k+v/qcQkZiKcfh/SMAYSE+IZNiAGIbXJWxD+0eT4YpusYzHsX+/BQsW+P79\nZs+eHfB/P5HW+FsQOJQSuRcAB/BNYDywHJhqrd3caN7XgdustV9qz/srkeue9P0Taai57U44vqIU\nqO1OAIfDQWs/QxwOBzU1bddRq66p5YuicrYUlrGlsIythYfZWljGgXaU9TiWsA3tfzxhGzag9YRN\nJJR1qc4Oxpgo4FpgjLXWDawyxrwKzALuCWpwIiJdSGdtd0LLvULrP99YaUUV2wrL2FqXsG0pLGP7\nvsMcra71+3OrDhVSVbQLW1JAQnglf318Hqel91fCJj1SSCRywHCg2lq7vd7YemB6C/MnGGOKgIPA\nYuA31trqjg1RRCT0+XMwIFCJXGZmJvPmzWv285zOPtx0x/f496ZCttStsG0tLCP/kP/t/5wRDkYM\n7MvopL5kJPRi7Vv/YulTj1K8d493u1jbnSKhsbVqjDkH+Ie1dmC9sW8BX7PWTm80dwhggV3AqcCL\nwGJr7W+aed/bgSwgzuVy9WvuXg1tzXVt+v6JNBSo7U5/+LZxd+ZR2zeJXv1PoVf/U+g9MINerlMg\nwun3ew3s62RUUgyjk/syKsn7Z3BilDocSI/VpbZWATfQt9FYLNCkYq21dke9y43GmAfwJmtNEjlr\n7UJgIXjvkQtYtCIiIepEtjvbo6qmlk/3HmbjnlI25JeQetvjlO8tw+JfHbWIMMPQ/jHepC3Ju9o2\nMqmvCueKnKBQSeS2A+HGmGHW2s/qxsYBm1t5zTEW0K9sIiK0td3pZPbs2X6/V02tZccBNxvyvUnb\n+vxSthSWUdnkfrbmk7j4yAjvCtvA46tsQ/tH0ytcxXNFAiUkEjlrbbkx5iXgAWPMbcAE4ApgauO5\nxpiLgY+ttfuMMSOBnwH/6NSARURCVFZWFkuXLm13HTRrLXkHj7A+v5SNdUnb5j2llFf6tw17Sr8o\nRif39a2yjUrqy4C+vTFGv2eLdKSQSOTqZAJPAfuBYmC2tXazMSYN2AKMttbmAecDzxhjooF9wF+B\nXwcpZhGRkBIdHU1OTk6rddCstRSWenwrbd5t0lJKK6r8+oyUuD6MHRTL2EFxjB0Uy5iUWGL7RHTw\n30xEmhMShx06g+rIhY7Bgwfz5JNP8uUvf/mk30vfP5G2FbmPsiG/pC5x8/4pcvtXo80V05txg2I5\nLSWOsamxnJYSS7/o3h0csYh0tcMO0kWsWLGCm266ifz8/GCHItKldFa3hcrqWjYXlLJ21yHW7jrE\n+t0lFJQ2X46ksbjICE5LiWXcoDhOG+T9qu1RkdCmRO5Eud2QnQ3z50NxMSQmQmYmZGVBCNU0qq6u\nbtL0vrkxEek4Hdlc/lB5JR/nHWLNrkOs3XmI9fklfhXXjeoVxpiUWMalxvmSt9SEPkraRLoYHR06\nEW43TJkC8+ZBURFY6/06b5533O0O2Eft3r2ba665BpfLRWJiInfeeSe1tbU8+OCDpKen079/f26+\n+WZKS0sB2LlzJ8YYFi1aRFpaGjNmzGh2DCAnJ4epU6cSFxfHuHHjWLFihe9zDx48yK233kpycjLx\n8fFcddVVlJeXc/HFF1NQUEB0dDTR0dEUFBRQW1vLb3/7WzIyMkhMTOSrX/0qBw8e9L3X4sWLSU9P\nJzExkV/96lcB+7cR6Sr86bbgD2stuQfc/H3Nbn78zw2c//sVTPjlf/nms2tYsCKX1TsPNpvE9Q53\ncHpaHF+fOpg/fHUc//eDc9n4iwt58dtn8ZNLRnH5uGTSEiOVxIl0QVqWORHZ2ZCbC42P93s83vHs\nbAhA5fSamhouu+wyZsyYweLFiwkLC2PNmjU888wzPPPMM7zzzju+RO7OO+9k8eLFvteuXLmSrVu3\n4nA42LdvX5OxPXv2cOmll7J48WIuuugi3nrrLa699lq2bduGy+Vi1qxZREdHs3nzZqKjo/nggw+I\niorijTfeaLK1+sgjj/DKK6+wcuVKXC4Xd911F9/5znd44YUX2LJlC7Nnz+b1119n8uTJ3HvvvdqW\nlR7nRLsteKpq2LinlDU7vdukH+cd4mB5ZZufl5rQh4npCZyRHs+EtDiGD4ghIky/t4t0RzrscCI3\ny7tc3hW41p5vpotEe3344YdcccUVFBYWNtgKPf/887n22mvJzMwE4NNPP2XMmDFUVFSQn5/PKaec\nQm5uLkOGDAG8q3SNxx566CE2bdrUIPm78MILufHGG7ngggtISUmhuLiY+Pj4BjE1d4/cqFGjeOyx\nxzj//PMBKCwsJC0tjYqKCn7961+zZcsWlixZAkB5eTnx8fG8/vrrOuwgPYa/3RYOHD5ad2/bQdbu\nOsSmPWVU1rS+TRruMJyaEsvE9HgmpsdzRno8/fv631FBREKTDjt0pOLik3veT7t37yY9Pb3J/WwF\nBQWkp6f7rtPT06murvatvAGkpqY2eb/6Y7t27eIf//gHy5Yt841VVVVx3nnnsXv3bhISEpokcS3Z\ntWsXV199NQ7H8d/4w8LC2LdvHwUFBQ0+Nyoq6qQry4t0Nc13WzBE9Euld8ooYjNOZ3r2O+wsPtLm\ne8VFRnBGWjyn1yVu41Lj1CxepAdTInciEhNbX5ELUKKSmppKXl5ek8MJycnJ7Nq1y3edl5dHeHg4\nAwYM8K2UNXevS/2x1NRUZs2axRNPPNFkXmFhIQcPHqSkpIS4uLgW36P+ez311FOcffbZTZ5LSkpi\n69atvusjR45QHKBEV6SryMzMZF7276iNG4QzbSzO1DH0ShlJmPP4AYeWkrgh/aI4Iz2eiYO9q21D\n+kXjUP9REamjmyZORGYmOFvYunA6oR0tcFozadIkkpKSuOeeeygvL8fj8fD+++8zc+ZM/vjHP/LF\nF1/gdrv5yU9+wvXXX9+uk6g33XQTy5Yt480336SmpgaPx8OKFSvIz88nKSmJiy++mMzMTA4dOkRV\nVRXvvvsuAAMGDKC4uNh3uALgjjvu4Kc//akvuTxw4ACvvvoqAF/5yld47bXXWLVqFZWVlfz85z+n\ntrbtE3UiXV1trWVrYRmLVn3BrrSLGTD7OZJu/gPx079On4yJDZK4Y3qFO5iYHs+3pw3hiZsnsva+\nL/P2j6aTfd04rj8zjaH9Y5TEiUgDWpE7EVlZsHRp0wMPTidkZHifD4CwsDCWLVvGXXfdRVpaGsYY\nbrzxRh5++GEKCgo499xz8Xg8XHjhhTz66KPteu/U1FReffVV5syZw8yZMwkLC2PSpEksWLAA8J40\n/f73v8/IkSOprKzkvPPO49xzz2XkyJHMnDmTIUOGUFNTw5YtW7j77rux1nLBBRdQUFBA//79uf76\n67nyyis59dRTefzxx7nxxhspLy/nBz/4AYMGDQrIv49IKLHWsrP4CB/kFvFBbjE5ucUU1z+YENH0\nlz8nVZw9IonJGf04Iz2BMSl96R2ubVIR8Z8OO5zozfLH6sgtWHC8jtzs2SFXR66702EHCabC0go+\n+LyYD3KL+TC3qM3CuwP7Opk6NJGzhiRy5uAE0lXyQ0RaoMMOHS062ltiJABlRkSkayh2HyVnx0He\nzy3iw9xivigqb3V+QlQvzspIZGpGIlMz+jFYiZuIBFjIJHLGmARgEXABUATca619vo3XvAXMACKs\ntdUdH6WIdCdttc0q81SxesdBPsgt5oPcIrbtPdzq+8X0DmfykATOyujH1IxERgzQPW0i0rFCJpED\nHgcqgQHAeGC5MWa9tXZzc5ONMV8DIjoxPhHpRpprm1VccpiHX3iD5ze5GTX9SjYVHKa2lbtPeoc7\nOHNwAmdlJHL20H6MSe5LuArvikgnColEzhgTBVwLjLHWuoFVxphXgVnAPc3MjwXmAjcDH3ZmrCLS\nPRxrm1UTM5C+48+kzymn0zt5JCY8gipgw56mq2/hDsP41DimDvWuuE1Ii9PhBBEJqpBI5IDhQLW1\ndnu9sfXA9Bbm/xpYAOzt4LhEpJupqKzhg9winlznJvHW+YT37d/iXGNgTHIsUzMSOSvDe0Ahqneo\n/LcpIhI6iVw0UNZorAyIaTzRGDMROBu4G2i1joUx5nYgC4hzuVwtzrPW6gbkLkj16MRf+YeO8M62\n/by1bT8f5hZztLqWiFEzmp1bWbQLz64NVO7eyJ6PVxAbqTs4RCR0hUoi5wb6NhqLBRrsbRhjHMB8\n4G5rbXVbyZe1diGwELzlR5qb43Q6fTc5K5nrGqy1VFVVsW/fPqKiooIdjoSg6ppa1u46xNuf7ued\nbfvZvs/d4twajxvPFx9Tkfs/KnZ+Qm15CQAul0tJnIiEvFBJ5LYD4caYYdbaz+rGxgGNDzr0BSYC\nL9YlXcduTsk3xlxnrX2vvR88aNAg8vPzOXDgwAmGLsEQHh5ObGws/fr1C3YoEiIOlleycvt+3t52\ngJWf7qfM0/JB9mH9o+lV/BnvLZnP4S/Wg224uut0OpkdoA4tIiIdKWQKAhtjlgAWuA2YACwHptY/\ntWq82duAei9LBVbj3WI9YK2tV0a9oZYKAotI12StZUthGe9s28/b2/bzye4SWvrvrFe4g7OGJHL+\nqP6cN6I/qQmRzZ5aBW8Sl5GRQU5ODtEq7i0iQdIVCwJnAk8B+4FiYLa1drMxJg3YAoy21uZR74CD\nMeZYz5t9qiMn0v0dqazm/c+LeXvbflZ8up/CVjopJMU6OW9kf2aM6M/UoYlE9mr43110dDQ5OTlk\nZ2ezYMEC3y0Ws2fP9tWRExEJdSGzItfRtCIn0jXlFR/h7W37ePvTA+TsKKayuvlDLg4DE9LimTHS\nu+o2KilG972KSJfVFVfkREQAWPfFPn7+5KusKwJHfEqL8/o6w5k2oj8zRrqYNrw/CVG9OjFKEZHg\nUyInIiHh8/1ulm8o5F/r8sktOgKk4IhvOm+oK5Ivj05ixsj+nJ4Wp04KItKjKZETkaD5oqic5RsK\neG1DYYt9TGurjuLJ20BF7v+ozd/ALbO/wT0X39/JkYqIhCbdIycinSqv+AivbSxg+YZCNhc0rgPu\nVVvloSL3fxzZtoqK3DXY6qO+51wuF/v37++scEVEgkL3yIlIyMg/dITlGwpZvrGQDfmlzc7pHe7g\nvBH9+etvvs+Rz1djq442O6+4uLgjQxUR6VKUyIlIhygoqeD1jYW8tqGQdbtLmp3TK8zBucNdXD4u\nifNHDSC6dzhLf7CV8haSOIDExMSOCllEpMtRIiciAbOvzMPrGwtZvqGQNbsONTsnIsxwzjAXl41N\n4sujB9DX2bANVmZmJvPmzWtQpPcYdVwQEWlI98iJyEk5cPgob2zyrrz9b+fBZrsrhDsMZw/tx6Vj\nk7hw9MBWe5iq44KIiO6RE5EAc7vdZGdnM3/+fA4dqcI14f+Rfs417Kcvtc0kb2EOw9SMRC49LYkL\nTx1IvJ813tRxQUTEf1qRE+nC6idXxxKezMzMgCc8brebyWedzR4S6T1qOs7B4zGOsCbzHAYmn5LI\nZeOSuOgVHm6RAAAgAElEQVTUgSRG9w5YDCIiPYm/K3JK5ES6qM7agsw94Oa7Dy9h0+FIHJGxTZ63\ntpaBDjd3XnEWF44ZSP8YZzPvIiIi7aGtVZFuLjs7u0kSB+DxeMjNzSU7O5v77z+xwrmeqhpe31jI\nktW7Wb3zIJCEI7LRnPzNHNm2iiOfvk9FnzBm/Ua13UREOptW5ES6KJfLRVFRUavPt7dw7paCMpb8\nL4+XP9nDYU91k+eryw7g3vhf3Bv+S03ZAd+4w+GgpqamXZ8lIiIt04qcSDfXVmFcfwvnuo9W8691\nBSz5X16zxXrDHIajO9ZQ/L9leL74GGxtkzmq7SYiEhxK5ES6qMTExFZX5FpLrqy1rNtdwpLVu1m2\noYAjlU1X09ITI7n+zFS+cvogHv/9R8x7dVOzSZxqu4mIBE+biZwxJgx4E7jUWttyuXUR6VQnUji3\n5EglL3+yhyWrd/PpvqZN6nuFObhwzEBmnpnKlCGJOBwGgKysLJYuXdriwYqsrKwA/s1ERMRfft0j\nZ4zZBYyw1jb9idFF6B456W78PbVqrSVnx0Fe/F8er2/aS2V101W1Yf2juWFSGtdMSGmx3tuxUieq\n7SYi0vECWn7EGPN14Bzg59baPScfXudTIifdUWvJVYWNYOnH+bz4v918UVTe5LV9IsK4bGwSN0xK\n4/S0OIwxQfgbiIhIcwKdyB37Fb7+ZANYa23TqqAnwBiTACwCLgCKgHuttc83M+8G4H4gCfAAbwDf\ntdaWtfb+SuSkJ6iptbz32QGWrN7N/23dR3UzLRdOS4nlhkmpXDEumRhny62yREQkeAJ9anXYScbj\nj8eBSmAAMB5YboxZb63d3GjeB8A0a+1eY0w08BfgQeCuTohRJCTtK/Pwwuo8/rEmnz0lFU2ej+kd\nzpUTkrnhzDTGpDQt6isiIl2TX4mctTb32GNjjMtae6C1+e1ljIkCrgXGWGvdwCpjzKvALOCeRrHk\nNXp5DTA0kPGIdBWf73ez8N1cXv5kD1U1TVffJqbHc8OkNC49LYk+vQKyeC4iIiHEr0SubuXrT8CN\nQC9jzFHgeeB71tqmR9/abzhQba3dXm9sPTC9hXi+BCwH+gJHgKsDEINIl/FJ3iH+vDKX/2zZR+O7\nI+IjI7j29EHcMCmVof1jghOgiIh0Cn+3Vv8EJACnA7uAdLzbmX8Cbg1AHNFA43vcyoBmfwpZa1cB\nscaYFOBbwM7m5hljbgeygDiXyxWAMEWCx1rLyu0H+PPKXHJ2HGzy/Bnp8Xx96mAuOHUAvcO1+iYi\n0hP4m8hdAmRYa48dfdtijLkF+DxAcbjxrq7VFwu0utpnrd1jjPk3sARvktn4+YXAQvAedghMqCKd\nq7qmluUbC/nzyh1sLWx6puf8kf25Y3oGZw5OCEJ0IiISTP4mch4gEahfwyABCFSB4O1AuDFmmLX2\ns7qxcUDjgw7NCQcyAhSHSMjwVNXwjzW7WfjeDnYfbHiAIcxhuHJcMt+elsGIgdo+FRHpqfxN5J4C\n/mOM+R3Ht1Z/gLdcyEmz1pYbY14CHjDG3AZMAK4Apjaea4z5GvCetTbPGJMO/Ap4KxBxiISC0iNV\nLM7ZydPv76S4vLLBc30iwrj+zFRuO+cUBsVHBilCEREJFf4mcr8E9gI3AclAAfAI8EQAY8nEmzDu\nB4qB2dbazcaYNGALMLruxOpo4CFjTDxwCHgduDeAcYgExd5SD4tW7eD5j/Iob9T7NC4yglvOGswt\nUweT0ELnBRER6XnaLAhc12v1p8BDXbnXqgoCS6hqrYRIcqyT284Zwg2TUons5e/vXSIi0tUFrCCw\ntbbGGHM33lU5EQmQ1kqIDB8QzbfPzeCK8clEhDmCE6CIiIQ8f3/F/yveMh8LOzAWkW7jWA/U+fPn\n+3qgZmZm8qMf/Yi1BRUsWJHLR180LSEyMT2eO6ZlMGNkfxwO9T4VEZHW+ZvIjQcyjTFzgN3U67lq\nrZ3REYGJdFVut5spU6aQm5uLx+MBoKj4IH965X0W7x9EbWxyk9eohIiIiJwIfxO55+r+iEgbsrOz\nfUmcCe9F1Glfpu+ka4iIG0htvXkqISIiIierzUSu7rBDCl38sINIZ5k/fz6eyipizriC2LO+SlhU\nXIPnbdVRbp02QiVERETkpPl72OEudNhBpE3WWsrjhpB81TeISBzU4LmaijIOr32N8nXL+cXvS4IU\noYiIdCf+bq3+DR12EGnVtr1lPPjaVvpf+/MG49Vl+ylb/QruDW9iq46ivr8iIhIoOuwgcpKK3Ef5\n/X+28+L/8qitV0ak9mg5pR8soWztMqipBsDpdDJ79uwgRSoiIt2NDjuInCBPVQ1Pv7+Tx9/5HPfR\nat+4w4Djiw8p+M8TVBza7xt3Op1kZGSQlZUVjHBFRKQb8iuRs9YGpKeqSHdgreWNTXv5zRtbmzSz\nP2dYP+67dDQp0dPIzg5jwYIFvjpys2fPJisri+jo6CBFLiIi3U2bLbp8E425FZgJuKy1E4wx5wAD\nrLX/7MgAA0UtuiQQNuSX8OBrW1m9s2Ex3wxXFPddOprpI1wYo0K+IiJycvxt0eVX7x9jzP14m9o/\nB5xSN1yAmtVLF+J2u5k7dy4ulwuHw4HL5WLu3Lm43e42X7u31MMP/r6OKx57v0ESFxcZwf1XnMq/\nv3cu543sryROREQ6lV8rcsaYPOAMa+0BY8wha2288f7EOmitje/wKANAK3I9W3PdFuD4fWs5OTnN\nbnlWVNbwl3dz+cvKHVRU1fjGwx2GW6YO5q4Zw4iNjOiUv4OIiPQc/q7I+XvYIRwoq3t8LPOLBtpe\nyhAJAfW7LdTn8XjIzc0lOzub+++/3zdeW2t5df0eHnrjU/aWNXzNl0cN4CeXjGSIS/e6iYhIcPm7\nIvcUUA78ANhnrU0wxvweiLLW3tHBMQaEVuR6NpfLRVFRUavP79/vPWG6dtdBHli2hfX5pQ3mjBwY\nw88vG83Uof06NFYREZFAr8h9H1gMlAK9jTFlwErgphMPUaTzFBcXt/n87oNH+O2/t7F8Q2GD5/pF\n9+ZHFwznuomphDl0D5yIiIQOf8uPlAJXGGOSgXRgt7U2v0MjEwmgxMTEFlfkTK8+DJxxC+f/YSWV\n1cfb2vcKd3Dbl04h87yhRPf293ceERGRztOun07W2gK8p1VFupTMzEzmzZvX8B454yD6tC8Td+4s\nwqLiGyRxl41N4scXjSQ1QU3tRUQkdGmZQXqErKwsli5d6jvw4EwbS/yM2+g1YEiDeeMGxfKzy0Yz\ncXBCkCIVERHxnxI56RGio6PJycnh/nkP88JnlvD00xs8P7Cvkx9fPIIrx6Xg0H1wIiLSRSiRkx7j\ng11u/i98EuHpVb6xPhFh3DEtg9vPHUKfXmFBjE5ERKT9lMhJt3ekspoHl2/l+Y/yGoxfe/ogsi4c\nwcBYZ5AiExEROTmtJnLGmFqOFwBu8jRgrbVaxpCQtWlPKXct+YQdB8p9Y0mxTv54/XimDEkMYmQi\nIiInr60VuWGdEoVIgNXWWp54bwe/+8+nVNUc/13k0tOS+PXVp6mtloiIdAutJnLW2tzOCkQkUI41\nuP8g93gR4MheYdx/xal85YxBamwvIiLdhu6Rk27l35v2cs9LGyg5cvxAw7jUOB65fjyD+0UFMTIR\nEZHAUyIn3cKRymoeWLaFJf/b7RtzGPjOeUO56/xhRIQ5ghidiIhIx1AiJ13ehvwSvrdkHTuKjh9o\nSInrwx+vH8+kU1TYV0REuq92JXLGmCRgkLX2fx0Uj4jfamotC9/dwe//8ynVtccPNFw+LpkHrxpD\nbB8daBARke7Nr/0mY8wgY8y7wBfAO3Vj1xpj/hKoQIwxCcaYl40x5caYXcaYG1uYd4sxZq0xpswY\nk2+MmWeM0cpiD1NQUsHXnszhoX9v8yVx0b3D+cNXx/GnG8YriRMRkR7B3xuH/gL8HxANHLuL/C3g\nwgDG8jhQCQwAvgYsMMac2sy8SOB7QD9gMnA+8KMAxiEh7vWNhVz8yHvk7DjoG5uQFsfrd53DNafr\nVKqIiPQc/q5kTQEut9bWGmMsgLW2xBgTF4ggjDFRwLXAGGutG1hljHkVmAXcU3+utXZBvcs9xpi/\nAecFIg4JbeVHq/nFvzbzj7X5vjGHgTtnDOOuGUMJ14EGERHpYfxN5PYDQ4DPjw0YY0YCu1t8RfsM\nB6qttdvrja0Hpvvx2nOBzQGKQ0LUut0lfG/JJ+wsPuIbS4nrwyM3jGfiYB1oEBGRnsnfRO73wDJj\nzK+BcGPMdcB9QHaA4ogGyhqNlQExrb3IGPMNYCJwWwvP3w5kAXEulysAYUpnq6m1/HllLn/87/YG\nBxquHJ/ML68aQ1+n7oUTEZGey69Ezlr7pDHmEPBtoLDu6y+ttf8MUBxuoG+jsVjgcEsvMMZcBfwG\n+LK1tqi5OdbahcBCgIkTJ7bUM1ZC1J6SCr7/4jpWf3H8Xrjo3uE8eNUYrpqQEsTIREREQoPfpz2t\ntUuBpR0Ux3a8K33DrLWf1Y2No4UtU2PMRcATwKXW2o0dFJME0bL1Bfzk5Y0c9lT7xs5Ij+fh68eT\nmhAZxMhERERCh9+JnDHmZmAmkAwUAEustc8GIghrbbkx5iXgAWPMbcAE4ApgajNxzAD+BlxtrV0d\niM+X0OE+Ws3PX93ESx/v8Y05DNx1/jDuPE8HGkREROrzK5EzxvwG+ArwJ2AXkAb81Bgz0lp7b4Bi\nyQSewnuwohiYba3dbIxJA7YAo621ecDP8G67vl6vzMR71tqLAxSHBMnHeYf43pJ15B08fqAhNaEP\nD18/njPSdaBBRESkMWNt27eOGWP2A2dYa3fXG0sD1lpru8QpgokTJ9o1a9YEOwxphrWW+Sty+cN/\nt1NT70DDNRNSuP/KU4nRgQYREelhjDFrrbUT25rn79aqGyhtNFZaNy5ywmprLXP/tZnFObt8YzG9\nw3nw6jFcOV4HGkRERFrjbyL3B+CfdVus+UAqMAf4Xd3KHAB1W58ifqmuqeX7Sz5m2cZ9vrGavZ9y\nQVoF5w89O4iRiYiIdA3+bq3W+vFe1lobdvIhdQxtrYaWo9U1ZC5ew1ufHq8cU75lBUXL/4izVwQZ\nGRnk5OQQHR0dxChFRESCw9+tVX+PAEb48afXiYUqPc2Rympue7ZhEnd43b8peu0PUFuDx+MhNzeX\n7OxA1ZsWERHpnvxN5C4Caq21Na396chApXsoraji5kWree+z40lc6UdLOfjmY2CPL/x6PB4WLFjQ\n3FuIiIhIHX8TuXlAoTHmYWPMGR0ZkHRfxe6j3PhEDmt2HfKNlby7mJIVTzc/v7i4s0ITERHpkvxK\n5Ky1pwKXAhZ4zRiz2RhzjzEmtUOjk25jb6mHr/7lQzYXHG+pezTnb5R++GKLr0lMTOyM0ERERLos\nv8vkW2vXWmu/D6QAP8Tb5eELY8zbxpjrTb3qvCL17Sou5yt//oDcA+WAt1PDvGvHcvu0YTidzmZf\n43Q6mT17dmeGKSIi0uX4dWrVN9mYdOCmuj+9gOeAPLxdGXZaa6/riCADQadWg2P7vsPc9ORH7D98\nFIBwh+HhG8Zz2dhk3G43U6ZMITc3F4/H43uN0+nUqVUREenRAnpq1RjzbWPMKuATIB34lrU2w1p7\nv7X2aWA6oBZZ0sDG/FKu/8uHviSud7iDJ26eyGVjkwGIjo4mJyeHOXPm4HK5cDgcuFwu5syZoyRO\nRETED/7Wkfs38CzwsrXW08KcS6y1rwc4voDRilzn+mhHMd98dg3uo9UARPcO58lbJjJliO57ExER\naUtAWnQZY5Zbay+11l7U1huFchInnWvFp/v59uK1HK32lhOJi4zg2VsnMS41LsiRiYiIdC9tteg6\np1OikG7j9Y2F3L3kE6pqvCu9rpje/PWbkxkxMCbIkYmIiHQ//vZaFWnTP9bs5sdLN1Bbt1ufEteH\nv902mcH9ooIbmIiISDfVViLnNMY819oEa+3NAYxHuqhn3v+CXyzb4rse0i+Kv942meS4PkGMSkRE\npHtrK5GzQG5nBCJdk7WWx9/5nN/9Z7tvbHRSX5775iT6RfcOYmQiIiLdX1uJ3FFr7f2dEol0OdZa\nfvvGNv7y7g7f2OlpcTx96yRi+0QEMTIREZGeoa1ETt0apFm1tZafvbqJv32U5xs7e2giC2dNJKq3\nbr0UERHpDG39xP1rp0QhXUpVTS1Z/1jPK+sKfGP/b/QAHp05AWdEWBAjExER6VlaTeSstWp2KQ14\nqmr47guf8N8t+3xjV45P5nfXjSMizO/WvSIiIhIA2gMTv5Ufreb2xWt4//Ni39iNk9N48MoxOBza\nhRcREelsSuTEL6UVVdz69Go+zivxjX172hDuuWgkxiiJExERCYYW98KMMS/We3xr54QjoajIfZQb\nFuY0SOKyLhyhJE5ERCTIWrup6UJz/Kf0I50RjISegpIKvvqXD9laWOYb+8Xlo/nOeUOVxImIiARZ\na1ur7wEfGmO200qHB3V26L52FpXztSc/Yk9JBQAOAw9dO5brJqYGOTIRERGB1hO564CvAOmow0OP\nc9hTxdefXu1L4iLCDI/cMIFLTksKcmQiIiJyTIuJnLXWQ10dOWNMhDo89BzWWn72yiZ2Fh8BoHe4\ngz/POoPzRvQPcmQiIiJSn1+nVq21vzDGDANmAinAHuAFa+1nHRmcBMfSj/c0KPY77ytjlcSJiIiE\nIL8quBpjLgfWAiOBg8AIYI0x5ooOjE2CIPeAm5+9ssl3fd0Zg7hyfEoQIxIREZGW+FtH7tfAldba\nd44NGGOmA48B/+qAuCQIPFU1fPf5T6ioqgFgiCuK+688NchRiYiISEv87ak0CO8p1vpW1Y1LN/Hb\nN7axpa7MSK9wB4/NPJ3IXqoZLSIiEqr8TeTWAT9sNPaDunHpBv67ZR/PfLDTd/3TS0YxOrlv8AIS\nERGRNvm73DIbWGaMuRvYDaQCR4DLOyow6TyFpRVk/XO97/r/jR7AzWelBzEiERER8Ye/p1a3GWNG\nAVOAZKAA+MhaW9WRwUnHq6m13L1kHSVHvN/KpFgn2V8Zq64NIiIiXYDfN0BZa6vx3hcn3cijb3/G\n6i8OAt7ODY/cMIG4yF5BjkpERET84e89ctIN5ewo5k9vHS8FePf5w5l0SkIQIxIREZH2UCLXQx0q\nr+R7S9ZRa73Xk09J4M4ZQ4MblIiIiLSLErkeyFpL1j83sLfMA0B8ZAQP3zCeMIfuixMREelK/O3s\n8EdjzPiODkY6x7Mf7OT/tu7zXWd/ZRxJsX2CGJGIiIicCH9X5MKAN40xm4wxPzbGqBBwF7W5oJRf\nv77Nd33r2YP58ugBTea53W7mzp2Ly+XC4XDgcrmYO3cubre7M8MVERGRVhhrrX8TjQkDLga+BlwG\nfAQ8B7xkrQ35n+4TJ060a9asCXYYQVV+tJrLH13FjqJyAE5N7stLmVPpHR7WYJ7b7WbKlCnk5ubi\n8Xh8406nk4yMDHJycoiOju7U2EVERHoSY8xaa+3Etub5fY+ctbbGWvuatXYm3npyLuAZYK8x5klj\njDqrh7i5/9rsS+Iie4Xx6MwJTZI4gOzs7CZJHIDH4yE3N5fs7OxOiVdERERa53ciZ4zpa4z5pjHm\nHeBdvCty5wCjADfwRseEKIHwyid7+OfafN/1L68cwxBX86tq8+fPb5LEHePxeFiwYEGHxCgiIiLt\n41dBYGPMP4EL8SZwfwZesdYerff8D4DSDolQTtrOonJ++vJG3/XVE1K49oyWb3MsLi5u9f3ael5E\nREQ6h7+dHXKAO621e5t70lpba4xpese8BF1ldS3ffeETyitrABicGMkvrxrT6msSExMpKipq9XkR\nEREJPn+3Vs9pLokzxrx07LG19kjAopKAmffvbWzc410sjQgzPDrzdKJ7t56/Z2Zm4nQ6m33O6XQy\ne/bsgMcpIiIi7edvIndeC+PTAxSHdIB3Pt3Pk6u+8F3/+KKRnDYots3XZWVlkZGR0SSZO3ZqNSsr\nK+CxioiISPu1ujRjjHmg7mGveo+PGQLs6pCo5KTtK/Pww7+v913PGNmfb37pFL9eGx0dTU5ODtnZ\n2SxYsIDi4mISExOZPXs2WVlZKj0iIiISIlqtI2eMebru4deAv9V7ygL7gEXW2s87LrzA6Ul15Gpq\nLbMWfcQHud5DCf1jevPG3eeQGN07yJGJiIiIP/ytI9fqipy19ta6N/vAWvtEoIKTjvXnlbm+JM4Y\nePiG8UriREREuqEWEzljzGBr7c66y7eMMUOam2et3dERgcmJWbvrIH/473bf9Z3nDWVqRr8gRiQi\nIiIdpbUVuY1ATN3jz/Fup5pGcyzePqwSAkqPVHHXC+uoqfVul09Mj+fu84cFOSoRERHpKC0mctba\nmHqP/e4AIcFhreWelzawp6QCgL7OcB6ZOYHwMH3rREREuiv9lO8m/vZRHm9sOl7qb95XxpES1yeI\nEYmIiEhHa+0euffwbp22ylp7bkAjknbbtreMX762xXd905Q0LhozMIgRiYiISGdo7R65JzstCjlh\nFZU13Pn8JxytrgVg5MAY7rt0dJCjEhERkc7Q2j1yz3ZmIMaYBGARcAFQBNxrrX2+mXljgN8DZwCJ\n1trGBzB6lAde28zn+90AOCMcPHbjBJwROn8iIiLSE7S2tTrLWru47vE3WppnrX0qQLE8DlQCA4Dx\nwHJjzHpr7eZG86qAvwPzgVcC9Nld0msbCnhh9W7f9f1XnMrQ/jGtvEJERES6k9a2VmcCi+sez2ph\njgVOOpEzxkQB1wJjrLVuYJUx5tW6z72nwQda+ynwqTFm6Ml+ble2++AR7l260Xd92dgkvjoxNYgR\niYiISGdrbWv1knqPz+vgOIYD1dba7fXG1gPTO/hzu6Sqmlq++8InHD5aDUBqQh9+fc1pGNOjd5lF\nRER6nFZbdNVnjIkDLgWSgQJgubW2JEBxRANljcbKOF6Q+IQYY24HbgdIS0s7mbcKKb//z3bW7fb+\n04c7DH+6YQJ9nRFBjkpEREQ6m1915IwxM4CdwF3AmcB3gZ3GmPMDFIcb6NtoLBY4fDJvaq1daK2d\naK2d6HK5TuatQsYHuUX8eWWu7/pHF45gQlp8ECMSERGRYPF3Re4x4HZr7d+PDRhjrsN7QGFkAOLY\nDoQbY4ZZaz+rGxsHND7o0OM98n+f+R6fM6wft5/TbAtcERER6QH87eyQDCxtNPYyEJCqs9bacuAl\n4AFjTJQx5kvAFRw/bOFjvJxAr7prpzGmdyDiCHWf7j3MR18cBCDMYXjo2rE4HLovTkREpKfyN5Fb\nDHyn0dhs4LkAxpIJ9AH2A88Ds621m40xacYYtzHm2E1u6UAFx1frKoBPAxhHyFqcs9P3+MJTB5Cs\nFlwiIiI9mr8tuhzAHcaYOcAeIAVvvbecQAVirT0IXNXMeB7ewxDHrncCPW4Z6rCnipc/3uO7njVl\ncPCCERERkZDQnhZdT3RkINK6lz/ZQ3llDQDD+kczZUhCkCMSERGRYAuZFl3SMmstz324y3c966x0\n1YwTERGRdtWRGwBMAvpRb2szgC26pAU5Ow76+qlG9Qrj6gkpQY5IREREQoFfiZwx5irgr8BnwKl4\nDxqMAVYRgBZd0rr6hxyuPj2FGBX/FREREfw/tfogcKu1dgJQXvf1dmBth0UWIMaYy40xC0tLS4Md\nygnZW+rhzc37fNc3nzU4eMGIiIhISPE3kUuz1v6j0dizwM0BjifgrLXLrLW3x8bGBjuUE/L86jxq\nar2HhyefksDwASfVtUxERES6EX8Tuf1198iBtzXXWUAGENYxYQlAVU0tL6zO811rNU5ERETq8zeR\newL4Ut3jPwLvAOuB+R0RlHi9uXkvBw4fBaB/TG8uOHVAG68QERGRnsSvww7W2ofqPX7OGLMCiLLW\nbu2owIQGJUdmTkojIszfvFtERER6gvaUHwkDpuDtu1pAALs6SFOf7j3M6rq+quEOw42T09p4hYiI\niPQ0/pYfGQu8AjiBfGAQ4DHGXG2tXd+B8fVYDfuqDmRAX2fwghEREZGQ5O9e3VPA40CKtXYS3l6r\nj6Each2icV/Vm6akBzEaERERCVX+JnLDgYettRag7usjwLCOCqwne+nj431Vhw9QX1URERFpnr+J\n3OvAFY3GLgeWBzYcsdayOKdeX9Up6qsqIiIizWvxHjljzGLA1l2GAUuMMWuB3UAqcAbwaodH2MN8\nuKO4QV/Vq9RXVURERFrQ2mGHzxtdb6r3eAvwZuDDCTxjzOXA5UOHDg12KH5ZXK/kyDWnD1JfVRER\nEWlRi4mctfb+zgyko1hrlwHLJk6c+K1gx9KWwtIK/rPleF/VWWfpkIOIiIi0rD115Kbj7a2aAuwB\nFltr3+mguHqkF1bv9vVVnTJEfVVFRESkdX4ddjDG3Ab8HdgLvAQUAi8YY0J+laurqKxu2Fd11pTB\nwQtGREREugR/V+TmAP+vfvFfY8yLwFK8fVjlJNXvqzqgr/qqioiISNv8LT+SiPeAQ32fAipwFiD1\nS4601FfV7XYzd+5cXC4XDocDl8vF3LlzcbvdnRmqiIiIhAh/E7lVwB+MMZEAxpgoIBv4oKMC60m2\n7S1r0Fd15qSmfVXdbjdTpkxh3rx5FBUVYa2lqKiIefPmMWXKFCVzIiIiPZC/idwdwFig1BizDygB\nxgHf7qjAepL6JUda6quanZ1Nbm4uHo+nwbjH4yE3N5fs7OwOj1NERERCS5uJnPG2FegDnA+cgrej\nwynW2mnW2oIOjq/bK/NU8fInx/uqtlRyZP78+U2SuGM8Hg8LFizokPhEREQkdLWZyNX1Vd0I1Fpr\n8621q621+R0fWs/w8sd7OFKvr+rkU5q/7bC4uLjV92nreRERCUFuN8ydCy4XOBzer3Pnese74ud0\nVyH872e8eVobk8z/b+/Oo+Qq6zSOf59OmLTQMUBoggRIJGERGBEIENBxARFUUGdwFpaARxYhB3Vc\n0kC66KgAABaUSURBVHAcJQIOzCQenaMIDiNrJOIGIjqiIAYBaTRRGY2szS5LaJBOGgkk4Td/3Nvp\nSqWrujtdy723ns85ddJ171u33nrrvbd+ee+76Hbg5Ii4t/5Zqo9Zs2bF0qVLm52NDUQEh335l+uX\n5Drv/Xsy56DpQ6bt7Oykt7e34rE6OztZsWJFPbJpZmb10N8Ps2dDTw+U3nFpb4cZM6C7Gzo68vM+\nRdWk8pO0LCJmDZdupH3klgA3Svq8pJMkfXjgMaZctrg7ewbXVe2YMJ6/33eHimnnzp1Le/vGfecA\n2tvbOf300+uSRzMzq5OFCzcODiB53tOT7M/T+xRVxstvpC1ylVZwiIg4pLZZqo8stsidtmgZNy5/\nGoATDprGue/fq2LagVGr5QMe2tvbmTFjBt3d3XT4f1RmZvnR2QlV7rTQ2Qm1uNPSqPcpqiaVX01b\n5CLiHRUemQ/iJB0l6ZK+vr5mZ2UDT/W9xE33lKyrOrv6uqodHR10d3fT1dW1wTxyXV1dDuLMzPJo\nuL7Nter73Kj3KaqMl1/VFrl03rjPAnsBvwUuiIiXG5S3mspai9yXfnYfX7nlQSBZV/WaUw9qco7M\nzKyh3CKXDzlvkfsayXQj9wIfBL5Yg7y1vFfWvsriXz++/vkJFQY4mJlZgc2dm3SYH0p7O9Sq73Oj\n3qeoMl5+wwVyRwDviogu4N3AkfXPUvH9dPnT9PYPrqt62B5eV9XMrOXMm5eMeiwPEgZGQ86bl6/3\nKaqMl99wgdwWEfEUQEQ8Dkyqf5aKr3Qlh0rrqpqZWcF1dCRTV3R1bTg/WVdXbae0aNT7FFXGy2+4\nPnJ/Bd4LKN30A+D9Jc+JiFvqmcFayUofuXufXskR/3UbkKyr+quzDmHbIZbkMjMzs9ZVqz5yK4DL\ngEvTx3Nlz78xxny2nKtK11XdazsHcY2U4Zm5zXLL55VZU41oHrkiyEKL3MrVa5h9/s/XL8n17VNn\nc+DOk5uap5bhmc3Nas/nlVnd1HplB6uBa5c9scG6qgdUWFfV6iDjM3Ob5ZLPK7OmcyDXIBHBou7B\n26pzDpqOpCqvsJq66KKNf2wGrF4NF1/c2PyYFYHPK7OmcyDXIL/qeY6eZ18E0nVV95na5By1mIzP\nzG2WSz6vzJrOgVyDlE45cvS+U+mYML6JucmQRnWUnjxMX8Th9pvZxnxemTWdA7kGKF9X9fhh1lVt\nGQMdpRcsSJY/iUj+XbAg2V7LYC7jM3Ob5ZLPK7OmcyDXAIvveox1ryajgw/aeTK7TJnY5BxlRCM7\nSmd8Zm6zXPJ5ZdZ0DuTq7JW1r/KtDdZVdWvceo3sKJ3xmbnNcsnnlVnTFX4eOUlHAUfNnDnzlAce\neKDh7//Du5/kY9/6HZCsq3r7mYd4Sa4BbW3J7dRq+9eta1x+zMzMMsLzyKUi4oaIOHXSpOYsE7vo\nzkfW/33sAdMcxJVyR2kzM7MxcVRRR/c8tZLfPPIXIFlX9ZgDdhzbAYu2FI47SudL0epfo7n8zKwO\nCn9rdUAzluj6zHV/YPFdjwFw5Btfx4XH7rvpByviUjhF/ExF5e9qbFx+ZjZKvrXaZCtXr+EHv/vz\n+udzxjrlSBGXwnFH6fwoYv1rJJefmdWJW+Tq5PI7HuacG/4EwG5TJnLjv/7d2Jbk6uxM5lirtn/F\nik0/vlk1rn9j4/Izs1Fyi1wTbbyu6rSxr6vqpXCsmVz/xsblZ2Z14kBujPr7+5k/fz6dnZ20tbXR\n2dnJR87+Eg/Vel1Vj/C0ZnL9GxuXn1l9eBCRA7mx6O/vZ/bs2SxYsIDe3l4igt7eXn54T9/6NEfv\nO5UtarGuqkd4WjO5/o2Ny8+s9hq5zGOGOZAbg4ULF9LT08Pqkg7M4yZuw4SdB29pz6nVSg5eCsea\nyfVvbFx+ZrXnQUSAA7kxueiiizYI4gA63nQEahsHwLon/8TMbWu0rqpHeFozuf6NjcvPrPYaucxj\nhnnU6hi0tbWxQfm1jWeH0y9nXMdWAPRefwH999xe0/c0MzMzCr/Mo0etNsDksg7Km+928Pogbu2q\nXjZ//sFmZMvMzKz4PIgIcCA3JnPnzqW9pM/LxH3fu/7v1X+8mdNP+0gzsmVmZlZ8HkQENDCQk7S1\npOskvSjpUUnHVkn7CUlPS1op6TJJE0r2LZG0WlJ/+rivMZ9gY/PmzWPGjBm0t7ezWefrad9hTwBi\n3Vq27X+Aee7AbGZmVh8eRAQ0tkXua8ArwBTgOOBiSXuWJ5J0OHAWcCgwDdgZOKcs2RkR0ZE+dqtv\ntivr6Oigu7ubrq4uOg8+ev3212/Wx69vvZkOd2A2MzOrDw8iAho02EHSFsBfgL0i4v5021XAkxFx\nVlnaxcAjEfGZ9PkhwOKI2C59vgT4ZkR8YzR5qOcSXX0vrWH2+T/npTVJp8rvfOQgDnj91nV5LzMz\nMyu+rA122BVYOxDEpe4GNmqRS7fdXZZuiqTSXosXSOqVdIekt9c8t6N07W+fWB/E7b7dRPafvlWT\nc2RmZmatoFGBXAewsmzbSmCoSdY6gL6ydJSkPZPkdutU4BLgBkkzhnpTSadKekDSs4899tim5n1Y\nW23+N0yfvDkAx8+uwbqqZmZmZiNQg7WjRqQfeG3ZtknAqhGknZT+uwogIu4q2XelpGOA9wBfLT9Q\nRFxCEuwxa9asut1D/sA+U3nf3ttz24O9zJrm1jgzMzNrjEa1yN0PjJe0S8m2vYHlQ6Rdnu4rTfdM\nRDxX4dgBNL0JrK1NvG3Xztqsq5oFXojYzLLE1ySzITVsZQdJ15AEXScD+wA/Bg6OiOVl6Y4ArgAO\nAZ4CrgO6I+IsSVsCBwK3AmuBfyZpcdunrP/dRuo52KFwBhYiLl/DbmBIdwuNBjKzDPA1yVpQ1gY7\nAMwFXgOsABYDp0fEckk7pfPB7QQQETcCC4BfAI8CDwPz02NsBnwBeBboBT4KfGC4IM5GyQsRm1mW\n+JpkVpHXWrWNdXZCb2/1/StWNC4/ZtbafE2yFpTFFjnLi+cqdUcc4X4zs1ryNcmsIgdytjEvRGxm\nWeJrkllFDuRsY16I2MyyxNcks4ocyNnGvBCxtQpPaZEPvibli8+rhnIgZxvzQsTWCgamtFiwIOlI\nH5H8u2BBst0/Otnha1J++LxqOI9aNbPWNH9+8uNSPqUFJC09XV1wzjmNz5dZnvm8qpmRjlp1IGdm\nrclTWpjVns+rmvH0I2Zm1XhKC7Pa83nVcA7kzKw1eUoLs9rzedVwDuTMrDV5Sguz2vN51XAO5Mys\nNXlKC7Pa83nVcA7kxsrz5ZjlU5GntPB1yZqlyOdVRnnU6lgMzJfT07PhUOuB/3m40ppZo/m6ZFYI\nHrXaCAsXbnyxhOR5T0+y38yskXxdMmspbpEbC8+XY2ZZ4+uSWSG4RS4l6ShJl/T19dX+4J4vx8yy\nxtcls5ZS+EAuIm6IiFMnTZpU+4N7vhwzyxpfl8xaSuEDubryfDlmljW+Lpm1FAdyY+H5cswsa3xd\nMmspDuTGwvPlWCWex8uaxdelsfP5azniUatmteZ5vMzyy+evZYRHrZo1i+fxMssvn7+WM26RM6s1\nz+Nlll8+fy0j3CJn1iyex8ssv3z+Ws44kDOrNc/jZZZfPn8tZxzImdWa5/Eyyy+fv5YzDuTMas3z\neJnll89fyxkHcma15nm8zPLL56/ljEetmpmZmWWMR62amZmZFZwDOTMzM7OcciBnZmZmllMO5MzM\nzMxyyoGcmZmZWU45kDMzMzPLKQdyZmZmZjnlQM7MzMwspxzImZmZmeWUAzkzMzOznHIgZ2ZmZpZT\nDuTMzMzMcsqBnJmZmVlOOZAzMzMzy6nCB3KSjpJ0SV9fX7OzYmZmZlZThQ/kIuKGiDh10qRJzc6K\nmZmZWU0VPpAzMzMzKyoHcmZmZmY55UDOzMzMLKccyJmZmZnllAM5MzMzs5xyIGeWZ/39MH8+dHZC\nW1vy7/z5yfY8K+rnMjOrMUVEs/PQELNmzYqlS5c2OxtmtdPfD7NnQ08PrF49uL29HWbMgO5u6Oho\nXv42VVE/l5nZKEhaFhGzhkvnFjmzvFq4cONgB5LnPT3J/jwq6ucyM6sDt8iZ5VVnJ/T2Vt+/YkXj\n8lMrRf1cZmajMNIWuZYJ5CQ9Czza7Hw02TZAlV/IlpL7stgP9hsuzTJYNoJDZaosavi5RitT5dBk\nLotBLotBLotEo8phWkR0DpeoZQI5A0lLRxLdtwKXxSCXRcLlMMhlMchlMchlkchaObiPnJmZmVlO\nOZAzMzMzyykHcq3lkmZnIENcFoNcFgmXwyCXxSCXxSCXRSJT5eA+cmZmZmY55RY5MzMzs5xyIGdm\nZmaWUw7kCkTSBEmXSnpU0ipJv5f07gppPyRpnaT+ksfbG5zlupG0RNLqks92X5W0n5D0tKSVki6T\nNKGRea2nsu+3P/3Ov1ohbaHqhKQzJC2V9LKkK8r2HSrpXkl/lfQLSdOqHGdrSddJejE9t46te+Zr\nrFJZSJot6SZJz0t6VtJ3Jb2uynFGfF5lVZWymC4pyur/56ocp8j14riycvhrWjZDzvGY93ox3G9n\n1q8XDuSKZTzwOPA2YBLwWeA7kqZXSH9nRHSUPJY0JJeNc0bJZ9ttqASSDgfOAg4FpgE7A+c0MI91\nVfr9AtsBLwHfrfKSItWJJ4EvAJeVbpS0DXAt8Dlga2Ap8O0qx/ka8AowBTgOuFjSnvXIcB0NWRbA\nViQdt6eT1P9VwOXDHGvY8yrjKpXFgC1LPt95VY5T2HoREVeXXTvmAg8Bv61yrDzXi4q/nXm4Xoyv\n58GtsSLiReDzJZt+JOlhkpnyH2lGnnLgRODSiFgOIOlcYDFJcFc0RwMrgNuanZFGiIhrASTNAnYo\n2fUPwPKI+G66//NAr6TdI+Le0mNI2oKk3PaKiH7gdknXA3PIUR2pVBYR8ZPSdJIuBG5tbO4aq0q9\nGLGi14shnAhcFQUdHTnMb+dkMn69cItcgUmaAuwKLK+QZB9JvZLul/Q5SUUL7C9IP98dVW4R7gnc\nXfL8bmCKpMl1z13jjeRiXPQ6AWXfeXoRfzDdXm5XYG1E3F+y7e4KaYvgrVS+XgwYyXmVZ49KekLS\n5WlrzFBapl6ktxHfClw1TNLC1Iuy387MXy8cyBWUpM2Aq4Ery//XkPolsBewLcn/II4B5jUuh3V3\nJslt0qkkt45ukDRjiHQdQF/J85XpvxPrm73GSi/GbwOurJKs6HViQPl3Dsn3PtR33sFgnRguba5J\neiNwNtW/85GeV3nUC+xPcot5P5Lv+OoKaVumXgAnALdFxMNV0hSmXgzx25n564UDuQKS1AYsIrlP\nf8ZQaSLioYh4OCJejYg/AOcCH2xgNusqIu6KiFUR8XJEXAncAbxniKT9wGtLnk9K/11V7zw22Bzg\n9moX46LXiRLl3zkk3/tQ3/lo0uaWpJnAT4CPR0TFW++jOK9yJyL6I2JpRKyNiGdIrp3vkjTUj3BL\n1IvUCVT/D2Bh6kWF387MXy8cyBWMJAGXknS0PDoi1ozwpQGobhlrvkqfbzmwd8nzvYFnIuK5huSq\ncYa9GA+hqHVig+887dcyg6FvKd4PjJe0S8m2vSukzaW0tfZm4LyIWDTKlxe1jkDy2WDo38nC1wsA\nSW8Gtge+N8qX5q5eVPntzPz1woFc8VwMvAE4KiJeqpRI0rvTfgBI2p1kRM71jclifUnaUtLhktol\njZd0HEkfjxuHSH4VcJKkPSRtRVIOVzQwu3Un6WCSWx7VRqsWrk6k3307MA4YN1AfgOuAvSQdne6f\nD9w9VBeEtD/MtcC5kraQ9BbgfST/a8+NSmUhaSpwC3BhRHx9mGOM5rzKrCplcaCk3SS1pX1kvwIs\niYjy22qFrxclSU4Evh8RFVuUilIvqPzbmf3rRUT4UZAHSd+OAFaTNPEOPI4Ddkr/3ilN+0XgGeBF\nkmHl5wKbNfsz1KgcOoHfkDRnvwB0A4el+zYoh3TbJ9OyWEky9cKEZn+GGpfHfwOLhthe6DpBMgot\nyh6fT/e9E7iXZDqWJcD0ktd9BvhJyfOtgR+k5fIYcGyzP1utyoLkRynKrhf9Q5VFtfMqT48qZXEM\n8HD6PT9F8p+87VqxXqT72tPv+dAhXleoekGV3850f6avF15r1czMzCynfGvVzMzMLKccyJmZmZnl\nlAM5MzMzs5xyIGdmZmaWUw7kzMzMzHLKgZyZmZlZTjmQM7OWIukKSV9o0nsrXYz9L5J+3Yw8mFmx\nOJAzs6aS9IikFenSNwPbTpa0pInZqpe3AIcBO0TEAeU7JX1I0u21eKO0XN9Zi2OZWXY5kDOzLBgH\nfLzZmRgtSeNG+ZJpwCORLOVjZjZmDuTMLAsWAp+WtGX5DknTJUXpGpCSlkg6Of37Q5LukPRlSS9I\nekjSwen2x9PWvhPLDruNpJskrZJ0a7pw/MCxd0/3PS/pPkn/VLLvCkkXS/pfSS8C7xgiv9tL+mH6\n+gclnZJuPwn4BnCQpH5J55S97g3A10v2v5BunyDpi5Iek/SMpK9Lek26bxtJP0o/9/OSbkvXCl1E\nsgTbDemxutK1ML8p6bk0/W8G1tY1s/xyIGdmWbCUZA3DT2/i6w8E/g+YDCwGrgH2B2YCxwMXSuoo\nSX8ccB6wDfB74GqA9PbuTekxtgX+BbhI0h4lrz0W+HdgIjDUbdBrgCeA7YEPAudLOiQiLgVOA+6M\niI6ImF/6ooi4p2z/QFD7H8CuwJvSzzMVODvd96n0vTqBKSTrPkZEzCFZ5/Go9FgLSBZAnwTsmJbT\naSRrR5pZjjmQM7OsOBv4qKTOTXjtwxFxeUSsA75NEqycGxEvR8TPgFdIgqABP46IX0bEy8C/kbSC\n7QgcSXLr8/KIWBsRvwO+D/xjyWuvj4g7IuLViFhdmon0GG8GzoyI1RHxe5JWuBM24TMhScCpwCci\n4vmIWAWcTxJgAqwBXgdMi4g1EXFbVF5Aew1JADczItZFxLKIWLkp+TKz7HAgZ2aZEBF/BH4EnLUJ\nL3+m5O+X0uOVbyttkXu85H37gedJWtCmAQemtx5fSG9vHgdsN9Rrh7A9MBBwDXiUpBVtU3QCmwPL\nSvJzY7odklvSDwI/S28pVyu7RcBPgWskPSlpgaTNNjFfZpYRDuTMLEvmA6ewYeAzMDBg85JtpYHV\npthx4I/0luvWwJMkQdqtEbFlyaMjIk4veW2lFi/SY2wtaWLJtp2AP48wX+XH7iUJQvcsyc+kiOgA\niIhVEfGpiNgZeB/wSUmHDnWstMXunIjYAziYpPVxk1oKzSw7HMiZWWZExIMkt0Y/VrLtWZJA6HhJ\n4yR9GJgxxrd6j6S3SPobkr5y3RHxOEmL4K6S5kjaLH3snw5EGEn+Hwd+BVyQDi54I3AS8M0R5usZ\nYIc0X0TEq8D/AF+WtC2ApKmSDk//PlLSzPQWbB+wDni15Fg7DxxY0jsk/W060nYlya3WgbRmllMO\n5Mwsa84FtijbdgowD3gO2JMkWBqLxSStf88D+5EMiCC9Jfoukj5oTwJPA/8JTBjFsY8Bpqevvw6Y\nHxE3j/C1twDLgacl9abbziS5fdotaSVwM7Bbum+X9Hk/cCdwUUT8It13AfDZ9Jbsp0laMb9HEsTd\nA9xKcrvVzHJMlfvFmpmZmVmWuUXOzMzMLKccyJmZmZnllAM5MzMzs5xyIGdmZmaWUw7kzMzMzHLK\ngZyZmZlZTjmQMzMzM8spB3JmZmZmOeVAzszMzCyn/h/DTmXPvwok1gAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# change default figure size and font size\n",
"plt.rcParams['figure.figsize'] = 10, 6\n",
"plt.rcParams['font.size'] = 12\n",
"\n",
"# plot results\n",
"n = np.arange(1, n_tests + 1)\n",
"plt.semilogy(n, p_rejects, 'ko', markersize = 8, label = 'raw')\n",
"plt.semilogy(n, 1 - (1 - alpha) ** n, linewidth = 3, label = r'$1-(1-\\alpha)^n$')\n",
"plt.semilogy(n, p_rejects_corrected, 'ro', markersize = 8, label = 'corrected')\n",
"ticks = [0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7]\n",
"plt.yticks(ticks, ticks)\n",
"plt.ylim([0, 0.8])\n",
"plt.xlabel('Number of tests')\n",
"plt.ylabel('Probability of Type I error')\n",
"plt.legend(loc = 'best')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From the plot we can see, without the correction, there's a rapid increase of the type-I error as the number of tests grows. Whereas applying the Bonferroni correction does succeed in controlling the error at the 5%."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Reference\n",
"\n",
"- [Blog: The multiple hypothesis testing problem](http://www.claudiobellei.com/2016/10/30/multiple-testing-part1/)"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
},
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "81px",
"width": "252px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": false,
"toc_section_display": "block",
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 1
}