why and when we need multiple comparison correction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_context(\"poster\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'r = 0.0882, p = 0.7114')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFkCAYAAAAe8OFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X2cZFVh5//PvbeququnpqcZmSecEURnzgAjMTSQGB+ARBEGiBqI0TVRf8ZoEn+bn3mxSyLZbDZuzDPr5kHjRje6ribZCKv+JCgYZUCMirS6OCBnEEQGYYYZhpmemu7qqrr37B/3Vk/19FN1d1XdutXf94umuqtuV52q6frWuefRc84hIiLd56ddABGR1UoBLCKSEgWwiEhKFMAiIilRAIuIpEQBLCKSklzaBZDeZox5D/AW4r+VTwC/b62dNXZxvuOMMWcAHwU2E3/g/4m19hPJ77wO+H0gAo4Av2KtfcQYUwQ+AFwMeMA3gHdZayc7+mRnP6cAuAm4Inlef26t/dAcx90MvLDpqucDd1lrf7bpmLcBr7PWXnPK7w4AtwL/zVp78ym3FYCvAJ+y1v55e56V9BLVgGVexpjdwOuBUWAXcBnw80s87g+Bb1hrf4w4yP7GGLM5CdlPAD9nrX0x8DngL5Pf+R3iwDs/+SoC7+nEc1zEO4EdxM/pIuDdxpiLTz3IWnudtfbFyfP4FeAo8C4AY8x6Y8yHgL8g/jCZZox5CfA14KXzPP5/Bc5u03ORHqQacB8yxlxK/IY/AZSAi6y1U023/yswdMqvfdVa+65Trnsd8PfW2hPJ730U+EXgn5ZwXACsM8Z4yWPWiWu8AXEgrUvuowRUku/vBh6z1kbJ/X0bOK+F510H/hi4ElgD3Git/d+nHDMC7Jnj1z9lrX3fHM/rb621deBZY8w/Js/r3nkevwD8D+Dd1tr9ydWvB54E/h1wzSm/8hvAbwM3znFfv0T82vzznE9W+oICuH/tAs621v7w1BustT/V4n1sA77U9PMTwNYlHvce4tPonwc2ANdba58GMMb8KvCvxphniAP5pUn57mjckTHmTODdwDtaKG8ATFhrR40x5wN3G2O+Yq091DjAWnsUeHEL99V4Xvubfn6CuEY+n18GnrTWfrrp8T6UPI+3nnqwtfaNyW0zAtgY8yLg/wMuIW6KkT6lAO5f++cKX1hSDdgHmtt7PSCc4y4XOu6TwJ9aa//GGLMd2GOM+TowCfxH4Nyk3fc3gFuMMS9utDEbY0aBTwN/ba29dZHn2/DXANba+40x3wVeAdzS9NyXUgNu9fk3/CatfVDMyxizDvg48CZr7QljzEruTnqcArh/lee7YQk14MeBM5p+PoO4FtjSccaY04GXAT+TPO7DxpgvEocixKH/SPL9B4D3A88BDhtj3gB8EPh/rbV/32J5IW7iaPA5JTCXWANu9fljjPlx4vfTXS2XdG6vBk4D/j4J3+cBrzLGDFtr/+MK71t6jDrhZCGfBd5kjFmT9Na/FfjMEo57hjiwrgNIAvkVxKMavgVcYozZlNzHa4EfWGsPG2OuIe6Qu3yJ4Qvw5uSxLgB2srJA/CzwNmNMLqk5v4G5nz/EzQVfnmuEyFJYa//JWntWU6fe/w+8X+Hbn1QDlnlZaz+XtEfeCxSIA+njMN1+e6G19u3zHZcMQ/tZ4K+MMb9L3Pn2R9baryT38WfETRJV4mFor0ke+s+JT/c/0nQK/lVr7buaH3eeYr/UGPMO4srFL1hrn13BS/A3wAuA/5M8r/9mrb0rKft7k9eoEYzbgcdW8FiyCnlajlKyxBizFvgra+1b57jNARustYe7XjCRZVAThGTNi4mHmolknmrAIiIpUQ1YRCQlCmARkZQogEVEUqJhaC0aGxv7r8QdQN8ZHR19d9rlEZHsUwC37sXAJaVS6RLiefqpGBsbA2B0dDStIqwqer27KyOvt7f4Ia1RE4SISEoUwCIiKVEAi4ikRAEsIpISBbCISEo0CkKkg5xz7DtY5uax/dy17xAT1ZChQsClOzZw3YXb2L6xhOe1rVNdMkYBLNIh9TDipjsst+09QOSgVAgYKeYJI8fnHzjI7Q8eZPeuzVx/uSEX6GR0NVIAi3SAc46b7rDcev9TDBfz+E213FzgMVL0iZzj1vufAuCGK3aqJrwK6WNXpAP2HSxz294Ds8K3me95DBfz3Lb3AA8/Pe8OUtLHFMAiHXDL2H4ix7zh2+B7HpGDW8bm3GpO+pwCWKQD9uw7RKkQtHRsqRBwp326wyWSXqQAFumAiWpI4LfWphv4HpPVhXa7l36lABbpgKFCQBi1tttMGDmKLdaWpb8ogEU64NIdGyi3WKstV0MuMxs7XCLpRQpgkQ64dnQbvgfRInsuRpHD9+Da0a1dKpn0EgWwSAfs2FRi967NjE/W5g3hKHKMV2rs3rWZ7RtLXS6h9AJNxBDpAM/zuP5yAzBjJlzge4SRo1wN8T24+vwtXH+50SSMVUoBLNIhucDnhit28roLtnLzffvZs+8Q5ck6xULA7l2buW50K9s3rU27mJIiBbBIB3mex45Na7nxqnO58aq0SyO9Rm3AIiIpUQCLiKREASwikhIFsIhIShTAIiIpUQCLiKREASwikhIFsIhIShTAIiIpUQCLiKREASwikhIFsIhIShTAIiIp6ZvV0IwxAfBrwFuAc4AAeBT4R+DPrLWVFIsnIjJLX9SAk/D9LPBXwE7g68Ae4AzgvcAeY8xQagUUEZlDXwQw8HbgKuB+YKe19pXW2iuB7cDXgJ8AfjfF8omIzNIvAfzW5PLd1tofNa601h4mbpYAeEO3CyUispB+CeDDwEPAvXPcti+5PKN7xRERWVxfdMJZa69Z4OaLkssnulEWEZFW9UsNeE7GGI+4Ew7gljTLIiJyKs85l3YZOsYY80fAbwMHgRdZaw8t977Gxsb2AJe0qWgiklGjo6Neu+6rb2vAxpj3EofvFPD6lYSviEgn9EUbcDNjTA74APAOoAL8nLX27nbdf6lUwhjTrrtbsrGxMQBGR0dTK8Nqote7u1bb691XAWyMKQGfAq4AjgKvaWf4ioi0U98EsDHmNOCLwCiwH9htrd2bbqlERObXFwFsjCkAtxGH74PAq621GnYmIj2tLwKYeKjZTxLXfC9Vh5uIZEHmA9gYsx74jeTHQ8D75+sks9b+YrfKJSKymMwHMHAxUEy+vyD5mo8CWER6RuYD2Fr7BaBtA6NFRLqlbydiiIj0OgWwiEhKFMAiIilRAIuIpEQBLCKSEgWwiEhKFMAiIilRAIuIpEQBLCKSEgWwiEhKFMAiIilRAIuIpEQBLCKSEgWwiEhKMr8cpYj0Pucc+w6WuXlsP3ftO8RENWSoEHDpjg1cd+E2tm8s4Xmrb1VZBbCIdFQ9jLjpDsttew8QOSgVAkaKecLI8fkHDnL7gwfZvWsz118+9042/UwBLCId45zjpjsst97/FMPFPH5TLTcXeIwUfSLnuPX+pwD4mQ1uVdWE1QYsIh2z72CZ2/YemBW+zXzPY7iY57a9B9g/Xu9yCdOlABaRjrllbD+RY97wbfA9j8jBnh9WulSy3qAAFpGO2bPvEKVC0NKxpULAt56a6nCJeosCWEQ6ZqIaEvittekGvsdUPepwiXqLAlhEOmaoEBBGrqVjw8gxkFtdkbS6nq2IdNWlOzZQroYtHVuuhlywZaDDJWqLDe26IwWwiHTMtaPb8D2I3MK14Chy+B5ceuZgl0q2fPd8//Dz2nVfCmAR6Zgdm0rs3rWZ8cnavCEcRY7xSo3duzazbbh3pyZM1UIOHJvkwSfH17frPhXAItIxnudx/eWGq8/fQnmqztHJGvUwwjlHPYw4OlmjXK1z9flbuP5y05OTMKr1iMPlKQ6OV5istdac0qre/bgRkb6QC3xuuGInr7tgKzfft589+w5RnqxTLATs3rWZ60a3sn3T2rSLOUstjBifrHFiqk7kHC32JS6JAlhEOs7zPHZsWsuNV53LjVelXZqF1cOI8UqdE1N1wijqSPA2KIBFRIiHwR2v1Dhe6XzwNiiARWRViyJHearG+GSdMHKEi4zYaCcFsIisSs45ylN1xidr1CPX8oSRdlIAi8iqU56qMT5RpxZFqQRvgwJYRFaN8lSN45N1qmFEFDnSi96YAlhE+ppzjhPVek8Fb4MCWET6UqNzrVwJqUW9FbwNCmAR6Sv1MOL49Dje7o5qWCoFsIj0hWo95HilzkQ1JIzcogsA9QIFsIhkVhQ5Jmp1TkyFTNXCjk0Z7hQFsIhkTqUacqIa13Yj53qyfbcVCmARyYRaGHFiqs7EVEg9ijJX252LAlhEelbWmxgWowAWkZ6Ry+UIgqBvmhgWowAWkZ4QRo7J0KMyFbL2eAXXZ7XduSiARSR15akaxyZqHJusUa+Hqa7P0E0KYBFJzVQt5OhkjUotDt1qrZ52kbpKASwiXRdGjmMTNcpTNULnyMCciY5QAItIVzWaG2phd3ad6GUKYBHpilObG0QBLCIdVq1HjFdqTEzVV3Vzw1wUwCJ9yDnHvoNlbh7bz137DjFRDRkqBFy6YwPXXbiN7RtLeJ7X0TJM1UKOT9WZ6OC27lmnABbpM/Uw4qY7LLftPUDkoFQIGCnmCSPH5x84yO0PHmT3rs1cf7khF/htf/yJap3jlXpfzlxrNwWwSB9xznHTHZZb73+K4WIev6mWmws8Roo+kXPcev9TANxwxc621IRn7TqhpoaWtP/jT0RSs+9gmdv2HpgVvs18z2O4mOe2vQd4+Onyih4vihzHJ2s8dbTCM+XqdAebwrc1CmCRPnLL2H4ix7zh2+B7HpGDW8aeWNbj1MOIoxNVnjw6yZGJKpVk9ppyd2kUwCJ9ZM++Q5QKQUvHlgoBd9qnl3T/1XrEkXKVp45VODpRoxqmu6171qkNWKSPTFRDRor5lo4NfI/y5OJTfxvtu/26JGSaFMAiHdbNIWFDhYAwcuSCxe8vjBzFBWrLWdhVOOsUwCId1O0hYZfu2MDnHzjISHHx+ypXQ3bv2jzr+ihyHK/UOF7p/V2Fs05twCId0jwkrDSQY6SYJxf4eJ5HLvAZKeYpDeS49f6nuOkOi2tD0F07ug3fY9EdgaPI4Xtw7ejW6eviBXKqPHVs8mT7rsK3oxTAIh3S7SFhADs2ldi9azPjk7V5QziKHOOVGj97/ha2byxRDyOePVHlqaNx8E7VFbzdogAW6ZBuDQlr5nke119uuPr8LZSn6hydrFEPI5xz8dCxyXj5xzdevI1fv+yFHDkRj2g4NllVjTcFCmCRDun0kLD55AKfG67YyYfffCFXnreJ0MVtukOFgDdetI2/fOOL+fkLn8eh41Mcr2hZyDSpE06kQzoxJKxVnuexY9NabrzqXG68avZQsvJUTbPVeoACWKRD2jkkbLnCZERDuRIvBamhZL1FTRAiHXLpjg2Uq2FLx5arIZeZjW177HoYceREPFW4ecaawre3KIBFOmQlQ8KWK4xcPKLhWIXxyWrSvqvY7VUKYJEOWcqQsN27NrN9Y2lFj1eeqnHg2CTjk+pYywq1AYt0SGNIGDBjJlzge4SRo1wN8T24+vwtXH+5WfZ05MZea1M1rUiWNQpgkQ5qDAl73QVbufm+/ezZd4jyZJ1iIWD3rs1cN7qV7ZvWLuu+tdda9imARTrs1CFhK6W91vqHAlikzdq9+lk9jJio1qnUI6q1KAldBW8/UACLtFE7Vz+rVENOVOtMVEOiKA5dZW5/UQAvUTWMdwTI5zzyvk8u8Ah8r+NbfEvva8eGmFO1kIlqyGQ1XoPXqabb1xTASxQPG6ri4ZH8Bx4EnofveQSBN/19zvfwA4+c55EPfHxfId3PlrL62e0PHOS60W28YGOJqVrIiWpIpRpSjyI1L6wiCuBliN8cjubzwXryg1fjZDATd8A0vs8HPgP5gMG8z0AuHo4k/aOV1c8C32PdQJ61gxH/fP+TXDe6TaG7iimA22xWNjeNDapHIVP1kOOVOJQLOZ9CLqCQNGf4flJ7VjBn0lyrn/ke5HyfwPcoFgIGcnEzRK0eH3/1j52RUmmlFyiAO8w5x2PPTPCFvQf45mNHmKyFFPMBF5+1nt3nb+bM9UP4nj+jOcPHI5/zyQc++SBuygh8n5zvMTAw0JadE6T9GqufBb5HMR8wlEy6aHwoR85RDSMmp+ocn6px5EQt7SJLyhTAHRRGjr+751HuevgwkYOhfMDwYI4wgru/f5h7HjnMJdtP520vO/uU5ghHLYzwkkT2mpo0Dk3Es6eePl4h8DwKgZ/UpH11BC6gGxtjnjEySDEfMJALcC7+95+qhdQjRxg5qvV4OFnkoB45BvPtX/1MsmXJAWyM+Vvg49baezpQnr7hXBy+d9pDlAZz+DT1iPswPJAjwnGnPQTA219+9owAcDRaL2a2NZcnKniex4lKPak1xwHtex6FnM9gzmcwHyiQm3R6Y8xaGFENI37q7Ofwr48+gwdUahETU/HY3blM1EIu2X76Cp+ZZN1yasBvB37ZGPM48D+BT1prbXuLlX2PPTPBXQ8fnhW+zXw8SoM57nr4MK86bzNnPWdNS/ftGuNBHbgkoEPiLWcmq3EY+54Xd/jlfAbycVvzatSOoWG+PzOUq/WQqXoUfzWtv/DjzzuNz3znSXKBN++/OcRNEb4Hr55jR2JZXZazGtp/Ah4GzgR+B3jQGPN1Y8y7jDH6SE/cntS2FnojQnx75OLjV8rB9OltNYwoV2o8c6LKgWMVnjgywaHjUxybqHKiWqdSDanWQ+phFC/S3aftyksZGvbF7x3k+4dmboyZzxeI/IBjE1UOHZ/iiSMTHDhW4Ui5yvFKvIFlo4lh62lFXnL2esqVOtE8UyYi5yhP1blk++mcuX6o7c9XsmXJNWBr7XuB9xpjLgR+CXg9cDFwEfBfjDFfIK4Zf85aO9XOwmbJvY8dYajFNr6hfMA3fnCEd17ygraWIW7GcNPD5hrtyo0a3ow4ah46Rzy5xPe9pDbNjEvPj+8jwMP3IUh6+XvRfEPDcn48NruQ8ykEPjkfjk/V+cJ3D7DppwaohY6pWsShiZAoinh2ooqb2Ro0i+d5vO1lZwPMaPcPfAijuNnB9+Ays4G3vexsNRHJ8jvhrLX3AfcZY34TeBVxGL8GuAa4Ghg3xnwK+IS19u52FDZLJmshw4OtvbyBDxOV1nZOWInpduVFa7sOwqaA9maGtdd0ReP6wPfIBX4SavGwul6YfHL/j47xvNOK0x8QuSAeTQLxy9D4gKpFDvDYY5/mmh/bMv1aTVTiOkSrY3QD3+PtLz+bV523mS989wD3PnaEiUrIYD7gku2nc8WuzZzZYlOT9L8Vj4Kw1kbA7cDtxpgB4KeJA/jNwC9zsr3474APW2tXfq6dAcV8QBjFHW6LCSN6skfcNX3jTr3llECqJ738nncyoD3iMbC5nIeHhyNeMjGKHPVkF4hGIOaDkwGeX0ZHWEMUOSr1kKmmRWuG8sF0bdO5uNOsMSqhFkZU6/FXhON4pb7iCRGe53HWc9bwq5e+gF+lvWc10l/aNgzNGOMBLwWuBHYDjY/5KeL24t8H3mOM+QNr7R+263F71cVnrefu7x9meGDxl7hfesTnGrlRj0K8+uzjGqr1aDqsG7MGvWQa9/TxDprP1l3S+dioyDdC3/NI1k842VF5fLJOuRKST+6vHkXUw7kXtenVD0LpXysOYGPMTwJvJG4L3kj8XnLA3cDHgU8BzwH+H+DfAf/ZGONZa9+30sfuZa/etZl7HjlMhFv1PeILVShn5PV004ij1qYWmRc9d92q+yCU7FhWABtjzgP+DfAG4CxONgU+TNwB9z+ttT9s+pXjwO8ZYx4CPgn8OtDXAXzWc4a4ZPvpc44Dbmj0iF9mNqhHvEP0QSi9bDkTMe4Hzkt+9IBngf9FPDnj64v8+leTy+GlPm6LZXslcCNwPlAAxoA/ttbe3onHW4h6xHuDPgilly2nBrwLqAGfJ25i+Jy1ttVJ7UPAPwHfWMbjLsgY81bgo8Rtzl8GAuAy4AvGmHdaa/+23Y+5GPWIp08fhNLLlhPAvwH8g7X2maX+orX2IeJmi7YyxmwBPgQcA15mrd2bXH8R8C/AXxhj/tla+6N2P/Zi1COeLuccjx+ZoDEj+PhEjWfCKfKBx8iaPK944QaufJE+CCUdy5mI8dedKMgK/VtgAPijRvgCWGu/aYz5U+APgHcAv5dS+SQFcy2GtOa0YLrmC/EY7K2nqdlB0tEvq6FdkVx+Zo7bPk0cwFeiAF41lrMYEjBr6VAvrHPuhjxrNp/gzPVDaqJYJZxzTNZCjk7UODYZf40nl6etGWjb42Q+gJPxx+cCEfC9OQ7Zl9x2XjL8rT8XPZAZlroY0mU7N3HnQwdnLR06MVHj/xys8tBn9s6zdKhkQbUeTYfo0cmToXqsKWCPVWaGbS2cOypu3H1O28qV+QAGTiNufjhkra2eeqO1tm6MOUw8RnktMN7l8kkKlrQYUuT4L1+0HJ2ozQps3/cY8mGgEMy7dKh0Vxg5ypU6xyZrHJ2scmyyngRnNbmszwrXyXYNLG+zfgjgRu/JxALHTCaXJVYYwBMTEzz+wAMruYu2eKAHytDL7n7oWQLPMVlZfIBOWI94+OAUW0o+U5X6nMdMVSrknOOLDzzJWYUyZ6zth7dO+pxzTIWOctVxohZRrjrK1Ygv/+A+yrWIcjXiRNVRriWX1YiJ2twzGbOoH/6KGiteL7hQ1SmX0ueqYcRQvrV/7olavCDPYrVaz/NwznHvk1O81vTDW6f9apHjRDWaFagnqhHlmpsVqCeqEfV+SdNl6Ie/osYCrsUFjhlMLk+s9MGGhoY47/nPW+ndLFuj5nveeectcuTqNvzNb+KfsqbEfKrlMrnAozg4+09oshKfPDVuyxccj4y7VfH6N5/qN76ONnVGNbelNq6bqPbmqf5CSgM51hXzDBfjy5FigeFijpFinuFinnVNXyNDeb7xg2fb9tj9EMDjxCF8ujEmZ62dcQ5pjMkBpwMVa+3RNAoo3beUxZDCCNYVe2vp0HZr9Oofm6xN9+yPN3dGzWhHrU/fnrXKaSHnx8E5mGfdUH46REemA7bAukbQDhUYHswtaxuqdsl8AFtrnTHmQeJF4XcAD55yiCHe+eO73S6bpGepa0CsGWhtFbReWTGtFkaze/HnqKU2/zxfr36v8j2SmmkcoOtOqZGODCXXJWG7bjDPYD5beyFmPoATXyAO4NcyO4Bfm1ze1tUSSaqWsgbECzeWeHaytdn0nVgxLXLxOsTzBepcX1k81V8zECSn+CeDdGYNNc+RA09QKvhc8KJzWDuYm3cbqX7RLwH8UeAG4LeMMbdba8cAkm2TbiAeBfHBFMsnXbaUNSAu27mJ//S5B9qyYlrzqf6M0/uJU4ZINY1JPV6prXgR+G7LBx4jxUJSG82xbqjp++T6Rjtqo+baykL7D0wdBOKa72rQFwFsrX3MGHM98AHga8aYLxGPePhp4uf4Zmvt02mWUbqv1cWQnHPz1pYbWxZN1kLKU3XO3bKW7zx+lLvsoXlrp1k81W+cxg8Pnjy1b/6KwzQJ2aE8g7lsner3qr4IYABr7QeTrY9uAF5OvCraPcD7rLVfSrVwkprmxZDe4c6ePtUfn6yx/8gke380Hodm0ra7/8gkkUu2Tpre0BQaQ8m/+sgRvvrIkdSeTysap/rNYTpXsDZO/Uur4FS/V/VNAANYa28Fbk27HNJ5zjkqtWhW59P0qf0cc/jHM3iqD3FNPvA9Ai8e7B54HmefvoZLd27ktKG4V7+5HXUle+pJd/VVAEt21cJoxvjS6e9nBWk9mX6azVP95l77GTXURoAO5rjTPs23Hz/K2sEcgefNOtWPcBwqT3H4eIXXvvgMNQVkmAJY2q65V3/OIK3MvG58ssaJjPbqz3lqP8/pfiu9+j84fIK9T46zbijf0iJCrzpvM2dpLePMUgDLghY71Z/RATWR3VP9fODNGG/a+JoqH6VU8Nj5/OfNGNi/rkOn+ktaRMjFx7/zEi30n1UK4FWmHs4M01OHSx2dqE6f6jeW6Ks2tpPICN+DtYNNY0ynO6JmDpc6OeW0MO8A/ump3zs3dqXs9z52hKEWJ3oM5QO+8YMjCuAMUwBnWORmz9U/1nTKP16Zff2Jqeyd6g8VgrlnPyU/N075G22rpYFcZtfsnayFDA/297RoOUkB3COcc1SaFo0+td208fXk4aOcqDoqd/1rZk/155pe2jwjat3QzFlSq6lXv5iPt0zKtfCUe2VatCyfArhD6mHEeCU+pW+eAdWY/dR8qj9eia9b2ql+b9R8BnI+69cUOL1UYF3TKlLxgP1kiFRSQ13oVF9iS1lEqBPToqW7FMBLdGKqzv3fP7zoXP0snuqvKQSzlt9bl/Te7/3RMR4+VMbDYygfUMh5OAeT9QjPg11nDGu7njZY6iJCC02Llt6nAF6ih58+we/dtT/tYixqvl79xnXNbaaN2+Y61XfO8ZGvPMrjRyZ4zprCrFDIB/6szS1Vw12+pSwidJnZwJnrtaNzlimAM8AD1g7mGBkqkHNV1uR9tm1aP91OOh2sTR1SxXzQliBc6uaWGpe6MktZROhtL1v+h51zbtYO0MV8wMVnreeKF23WDtBdogBOQTEfxGE5Rw11ZqdUPESqNHiyV//kjhimK2XVuNTua3URoeUKI8ff3fPorB2gwwju/v5h7nnkcN/sAN3rHzQK4BXK+d6MXvv5evgbq0gND+YptNLF3SN6bVxqr7+h2qV5EaFfpX2vp3Nx+M7VxJHzYXgg1zdNSln4oFEAL9HZG4b4wL85c3o8artO9XtVL41LzcIbqtetlialrHzQZKcq1iOGB/Ocs2WYLeuKDBVyfR2+cHJcais6OS61+Q01VAgYHsiR8z084o03hwdyDBUC7rSH+Lt7HsW5jA2Q7pLlNCll0VI/aH54ZKLLJWyUQfqGc44fHD7B3+x5hLd97Ju88cNf520f+yYf2vMIjz1zYlmhdPFZ65motVarnaiF/MTz1y/5MVqRlTdUr1tOk1IWZeWDRgHcJ8IoHi72u5/dy1e+fxjf8xhOVt+6+/uH+Q+lkL4cAAAUBUlEQVSf2ctHvvIo4RKnzr1612Z8L14CcSGdHpealTdUr5ushbQ6sTDwodLih2+vycoHjQK4D3Ty9LwxLrVcqc8bwo1xqZdsP71j41Kz8obqdb3SpNRpWfmgUQD3gU6enjfGpV5mNjBRDRmfqlOPHA5HPXKMT9WZqIUrHpe6mKy8oXpdrzQpdVpWPmgUwH2g06fnjXGpf/DaXbzihadPL7geJZtZvu+1u/iVV7ygoyMPsvKG6nW90qTUaVn5oNEwtD7QjbG6nRqX2iotUtMeq2Wqc1bW1FANuA+shtPz1VJz67ReaVLqtF7pu1iMasB9YDWsIbtaam7d0Ompzr2gW2tqrJQCuA+shtPz+d5QAzmPnO9TjRxDOZ9XnrORX7joeXF7tAdeEtTNIz/mrA85iOuBje/7W9pNSt2QhQ8aBXAfyEp710oFvsevvPxsrti1hS9/72m+d2AcgMG8z4VnrudV523i7NNLeJ6H78XHN2o2YeSSMdCOU0fhRS4+/W4cU48cUfL9muIgzjlyvtcUznPH82oJ7yzp9Q8aBXAf6NfTcy+uxOJ5cT028D1ygc+5ZxR40XNHGCoEDBZaa04JfG9ZozSeHfJxwOZ1RZxz1KKIWj2O4KSSjQOiKA7m0DnCMA7vaIHx1o1bXJLa7pQbFOKrgwK4D2SlvWsuvncyYPHioXK+5+H7kAt8cr5HPuczEPjkUtgbrlqtAkyvYDdAAAOt/W6U1KZjTU0gLv4KcURhHNqNAHfEtfVaPSJyJ4PZozm0V09TSb9TAPeJLLR3NXhJ22zgewwVguldjJubDPqB73sUljk22jlHNQlh34tfl0YNPIoctTBpNgkddRdNh3GjRp21zVpXKwVwH+nV9i4vqdl6SbtsIRd3ng0VcqnUarPA8zwG5hitMsDs65yLAzmMIuqRo1Z3VGoh9SiaEcqqNfceBbB0hO8RNyUkQTKY8xnI+xRy2RsC1+s8z6OQ8zh1WH8tjKiFUVxLTtqk4+aOk23Uzs1u2tBKnt2jAJYVizvK4mYFz4Oc71MsBAwVAgo5bUOflnzgz7nRarN6GMVtzlFEPQnqWj2uSZ8M5DiR3fT/VItuFwWwLJmX1G49D/K+P91Zlgs8CoE/56mz9KZc4JMLZjdtNJo16lFEPYyoRydryGHk4uBOqs8zQlpNHEuiAJaWBZ6Hn3ScDRUCBvt8O6bVrNGsUVhgtYJ4lEc03RlYa4yjTmrV0YzRGmqDnosCWBYVePEIhdJgjrWDee23JkBjlEdAYZ7b49Ea0XTzRmOSSyOg4WTtuRHOjdEeq4UCWOblJ+25a5Pg9RW8sgS+7zHgB3OO3GjUlOvOEYYR9RDqLqJUHCB0jsD3VkWnoAJYZvGTGu+agRzriovXeJ1z7DtY5uax/dy17xAT1ZChQsClOzZw3YXb2L6xpKYKmSEe9x3MmtNysBjPPNy0dpBqGFIL4/HQtSiansDi3HyTwbNHASzAyRlpgbe0poZ6GHHTHZbbkkXhS4WAkWKeMHJ8/oGD3P7gQXbv2sz1lxuN+ZVFNWYeDhYCBptqzo3mjGo9YrIeMlWLZo/SyGD7sgJ4FWuMZsgnw8aK+YCBfOvDxpxz3HSH5db7n2K4mMdv+r1c4DFS9Imc49b7nwLghit2qiYsyzLdnJEPWEt+5uSTMP6+FsVjnkMXzTm+uRc7ABXAq5DnxR1rxUKO0kBAsbC8P4N9B8vctvfArPBt5nsew8U8t+09wOsu2MqOTWtXUnQRYP7JJxAHbmNFu8bU7UYbc2NSSnP7ssOlFs4K4FWksaLYQD5uJljpeN1bxvbHe9EtUqv1vXgvulvGnuA9u89Z0WOKLMbzPPKBB3OMb25oLHhUCyOqyYzBWuhm1JajLgSzAngVaNR4B/IBawdzDC2zxnuqPfsOUWpxOchSIeBO+7QCWHpC4HsEp7QzN2rOYXhy2dFqGFINo5Odf21OYwVwH2sO3uHB3LKbGuYzUQ0ZKeZbOjbwPcqT9bY+vkg7NWrO+YAZwRxGjqmk468aRgS+1+L+3ItTAPeZwJu56tjagVzLi5Yv1VAhIIwcuWDxjrUwchQ7VA6RToqXTc0xlMw4+eWXPf/xdt23xgX1Cd+LF19ZN5Rny0iR5542xIa1Ax0LX4BLd2ygXG1th+VyNeQys7FjZRHpoofbdUcK4IyLVx/zWDuYZ8u6QUaGCouugNUu145ui7eKX6RhLIriveiuHd3alXKJZIUCOKMazQzFfMCm4UGeUxro+kSHHZtK7N61mfHJ2rwhHEWO8UqN3bs2s31jqavlE+l1agPOEM+DwYECge+zdjBPsRCwps0da0srj8f1lxuAGTPhAt8jjBzlarwX3dXnb+H6y40mYYicQgGcAfEW6z5DhYDTigEDgcfppRZ3huywXOBzwxU7ed0FW7n5vv3s2XeI8mSdYiFg967NXDe6le2afCEyJwVwD2tMnFgzkGNkqEDgezzuQqrVto2CaQvP89ixaS03XnUuN16VdmlEskMB3KMCP95dYv2awowZa1HUW+ErIsunAO4xjTV41xXzlAZzfdduqqUrRU5SAPeIRnPD0ECOkWK+L5du1NKVIjMpgHtAo7lhZCjf9unCvUJLV4rM1p/v9oxoNDcMF+MF0Ps5cLR0pXRbFpq7FMApaN7yZ3gwtypOt7V0pXRTVpq7+v+d30MC72RTwxkjRdavKayK8IXlLV0pshzNzV2lpj4Vz/PIBT4jxTylgRy33v8UN91hU92FeXW8+1M0vQh6zue0NXHwNsb0riYT1bDl5xz4HpMtLvIjcqqlNnc9/HS5yyVsKkdqj9znTgZvwPqhAmeMFBkuFlbt1u6NpStboaUrZSWW09yVFrUBt5lHvIFgIfBZW8yxptB/Y3mX49IdG/j8AwcZKS7+mV+uhuzetXnR47LQySLdl6WdWhTAbdLYfaKQ81lbzKe6SE4vunZ0G7c/eJDIuQVrJq0uXZmVThbpvizt1KK/zBXwvZNjeIcH82xYO8DmdUWF7xzauXRlljpZpPuy1NylAF6GXLIO7/BggQ1rB3juaUWeUxro20kU7dBYuvLq87dQnqpzdLJGPYzijRDDiKOTNcrVektLV2apk0W6L0s7tSiAlygf+GwZKbJlpMj6UoEhtfG2rLF05YfffCFXnreJ0DmOTdYInWP3rs185M0X8ltXnrNok0GWOlmk+7K0U4uqbEsU+F7XtvzpR+1YujJLnSzSfY3mrrmmvTc0mruuPn9Lqju1KIAlc7LUySLdl6WdWhTAkjmNTpZcsPgbJ+1OFklHVnZqUQBL5nRiTLH0nyzs1KLGTMmcLHWyiCxEASyZ084xxSJpUhOEZE6WOllEFqIAlkzKSieLyEIUwJJZWehkEVmI2oBFRFKiABYRSYkCWEQkJQpgEZGUKIBFRFKiABYRSYkCWEQkJQpgEZGUKIBFRFKiABYRSYkCWEQkJQpgEZGUKIBFRFKiABYRSYkCWEQkJQpgEZGUZH5BdmNMEbgeeD3wQsABDwEfAz5grY3SK52IyPwyXQM2xpSAu4H/DDwXuAv4OrAT+EvgfxtjgvRKKCIyv0wHMPA7wIXAl4AXWGuvtNb+DHAe8DDwGuDtKZZPRGReWQ/gtySX77DWHm1caa19DPj3yY9v6HahRERakdk24KT54WFgv7X20TkO2ZdcntG9UomItC6zAWytLQOXLHDIRcnlE10ojkjXOefYd7DMzWP7uWvfISaqIUOFgEt3bOC6C7exfWMJz/PSLqYswHPOpV2GtjPGDADfBbYD77LWfnCl9zk2NraHhQNfpGvCyPEPD5T52v4KDijmwPcgcjBZB8+Dl2wd5I3nlQh8hXA7jY6Otu0F7akasDHmk8BoC4d+2lr7nnnuwwP+O3H4fi/5XqRvOBeH71f3T7Im7+E31XIDD0oFiJzjq49PAvCmXaoJ96qeCmDgTMC0cNyWua5Mhpx9BHgT8Czwc9baqfYVD0qlEsa0UsTOGBsbA2B0tJXPKVmpXny97YHjfOvO+9iwrjQjfE81NOT41qE6v7bNsGPT2i6WcPl68fXupJ4KYGvty5b7u8aYNcA/ANcAR4DLrbUPtatsIr3ilrH9RI4Fwxfi2yMHt4w9wXt2n9Ol0slSZH0YGgDGmE3EkzCuIe50e4W1dizdUol0xp59hygVWptfVCoE3Gmf7nCJZLkyH8DGmDOBrxG3HX8XeIm19oF0SyXSORPVsOWOtcD3mKyGHS6RLFemA9gYsx74F+D5xFOSX26t1bAz6WtDhYAwam30Uhg5ii3WlqX7Mh3AwAeJF+D5NnCFtfZYyuUR6bhLd2yg3GKttlwNucxs7HCJZLl6qhNuKYwx5xCvgAZwHPjwPKMTDllrf7NrBRPpsGtHt3H7gweJnFuwIy6KHL4H145u7WLpZCkyG8DEkyIaf32vWOC4HwIKYOkbOzaV2L1rM7fe/xTDxfycIRxFjvFKjavP38L2jaUUSimtyGwAW2s/BHwo7XKIdJvneVx/eXy2d9veA0QuHu0Q+B5h5ChXQ3wPrj5/C9dfbjQJo4dlNoBFVrNc4HPDFTt53QVbufm+/ezZd4jyZJ1iIWD3rs1cN7qV7RmZfLGaKYBFMsrzPHZsWsuNV53LjVelXRpZjqyPghARySwFsIhIShTAIiIp6cv1gDthbGzsCeC5vu8zNDSUWjnK5TIQr8omnafXu7uy8HqXy+W7gO+Mjo6+e6X3pU641pUAoiia/iNJUy+UYTXR691dPf56t21jBgVw635AvOZEGfh+ymURkXR9px13oiYIEZGUqBNORCQlCmARkZQogEVEUqIAFhFJiQJYRCQlCmARkZQogEVEUqIAFhFJiQJYRCQlCmARkZQogEVEUqIAFhFJiQJYRCQlWo4yo4wxReB64PXACwEHPAR8DPiAtTZKr3T9wRjzSuBG4HygAIwBf2ytvT3VgvUZY0wA/BrwFuAcIAAeBf4R+DNrbSXF4nWUlqPMIGNMCbgTuBA4AtxLHBA/CQwBnwWutdaGqRUy44wxbwU+CkwBXyYOhcuAPPBOa+3fple6/pGE72eBq4jX2v4GUCP+Wx5Jfv5pa+1EaoXsIDVBZNPvEIfvl4AXWGuvtNb+DHAe8DDwGuDtKZYv04wxW4APAceAC621u621rwZeCowDf2GMeW6aZewjbycO3/uBndbaV1prrwS2A18DfgL43RTL11EK4Gx6S3L5Dmvt0caV1trHgH+f/PiGbheqj/xbYAB4v7V2b+NKa+03gT8FBoF3pFS2fvPW5PLd1tofNa601h4mbpaAPv5bVgBnTNL88DBwr7X20TkO2ZdcntG9UvWdK5LLz8xx26eTyyu7VJZ+d5i47+LeOW7r+79ldcJljLW2zMKbAl6UXD7RheL0HWOMB5wLRMD35jhkX3LbecYYz1qrTpQVsNZes8DNff+3rBpwHzHGDAD/IfnxljTLkmGnETc/PGOtrZ56o7W2TlxrGwLWdrlsq0byQfje5Me+/VtWDbgHGGM+CYy2cOinrbXvmec+POC/E3defC/5XpZuTXK5UK/7ZHJZIu6Uk/b7Q+IzvYPAn6Vclo5RAPeGMwHTwnFb5royGcrzEeBNwLPAz1lrp9pXvFWlMX56oaYF75RLaSNjzHuB3yYeAvh6a+2hlIvUMQrgHmCtfdlyf9cYswb4B+Aa4jHBl1trH2pX2VahcnJZXOCYweTyRIfLsqoYY3LAB4hHmFSIKxJ3p1uqzlIbcIYZYzYBdxGH7xPAK6y1Y+mWKvPGiUP49CQQZkiuOx2oNA8BlJVJRvd8jjh8jwKvttZ+Pt1SdZ4COKOMMWcSD1QfBb4LvMRa+0C6pcq+ZFTDg8Qz33bMcYghft98t5vl6mfGmNOAPcTD//YDL+/3mm+DAjiDjDHrgX8Bng/cTfwH27dDdVLwheTytXPc1rjuti6Vpa8ZYwrEr+Uo8QffTzVPful3CuBs+iDxAjzfBq6w1h5LuTz95qPEbZC/ZYyZHp1ijLkQuIF4FMQHUypbv3kv8boP+4FLV1tFQovxZIwx5hzgAeIe+LuJ/3Dncsha+5tdK1ifMcb8OnGHUI14zQ0P+Gnijus3W2s/kWLx+kJyJvcEcYfnt5h74gsA1tpf7Fa5ukmjILLnEk4Of3rFAsf9EFAAL5O19oPGmMeJa7wvJx4SdQ/wPmvtl1ItXP+4mJOjTS5IvubTlwGsGrCISErUBiwikhIFsIhIShTAIiIpUQCLiKREASwikhIFsIhIShTAIiIpUQCLiKREASwikhIFsIhIShTAIiIpUQCLiKREq6GJzMEY8yLgPqAA/La19k9Ouf2XgI8Tr5J2kbVWO2TIkqkGLDKHJFB/P/nx94wxz2/cZozZCvxl8uONCl9ZLgWwyPz+BPgm8Zq1HwQwxnjEO2aMAF8G3p9a6STztB6wyAKSHUi+DQwAbwDWE4fxs8D5q20LHWkvBbDIIowxNxDXhp8CSsBa4A3W2v+VasEk89QEIbK4Pwe+BmwhDt9PKnylHRTAIouw1kac3KoeQJ1u0hZqghBZhDFmJyfbgSEeevbj1tqH0iuV9APVgEUWYIwJgI8Bg8CHgY8k3/+P5DaRZVMAiyzsBuAniDvgfiv5+QDxluo3pFgu6QNqghCZhzFmF/FsuAHgOmvtLcn1vwD8I1AFLrDWPpBeKSXLFMAiczDG5ICvA6PAZ621rz3l9tuAK4Ex4CettfXul1KyTk0QInO7kTh8x4F3zXH7rwMnkmPe08VySR9RDVhEJCWqAYuIpEQBLCKSEgWwiEhKFMAiIilRAIuIpEQBLCKSEgWwiEhKFMAiIilRAIuIpEQBLCKSEgWwiEhK/i/blMevL2712QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def generate_date(means = [0,0],cov = [[1,0],[0,1]],size = 20,):\n",
    "    return np.random.multivariate_normal(means,cov,size = size)\n",
    "# test\n",
    "temp = generate_date()\n",
    "df = pd.DataFrame(temp,columns = ['x','y'])\n",
    "g = sns.lmplot('x','y',data = df)\n",
    "r,p = stats.pearsonr(df['x'],df['y'])\n",
    "g.fig.suptitle(f'r = {r:.4f}, p = {p:.4f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = dict(r=[],p=[])\n",
    "for ii in range(int(1e4)):\n",
    "    temp = generate_date()\n",
    "    df = pd.DataFrame(temp,columns = ['x','y'])\n",
    "    r,p = stats.pearsonr(df['x'],df['y'])\n",
    "    results['r'].append(r)\n",
    "    results['p'].append(p)\n",
    "results = pd.DataFrame(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "level of false alarm = 0.054\n"
     ]
    }
   ],
   "source": [
    "level = np.sum(results['p'] < 0.05) / results.shape[0]\n",
    "print(f'level of false alarm = {level:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels as sms\n",
    "class MCPConverter(object):\n",
    "    \"\"\"\n",
    "    input: array of p-values.\n",
    "    * convert p-value into adjusted p-value (or q-value)\n",
    "    \"\"\"\n",
    "    def __init__(self, pvals, zscores=None):\n",
    "        self.pvals = pvals\n",
    "        self.zscores = zscores\n",
    "        self.len = len(pvals)\n",
    "        if zscores is not None:\n",
    "            srted = np.array(sorted(zip(pvals.copy(), zscores.copy())))\n",
    "            self.sorted_pvals = srted[:, 0]\n",
    "            self.sorted_zscores = srted[:, 1]\n",
    "        else:\n",
    "            self.sorted_pvals = np.array(sorted(pvals.copy()))\n",
    "        self.order = sorted(range(len(pvals)), key=lambda x: pvals[x])\n",
    "    \n",
    "    def adjust(self, method=\"holm\"):\n",
    "        \"\"\"\n",
    "        methods = [\"bonferroni\", \"holm\", \"bh\", \"lfdr\"]\n",
    "         (local FDR method needs 'statsmodels' package)\n",
    "        \"\"\"\n",
    "        if method is \"bonferroni\":\n",
    "            return [np.min([1, i]) for i in self.sorted_pvals * self.len]\n",
    "        elif method is \"holm\":\n",
    "            return [np.min([1, i]) for i in (self.sorted_pvals * (self.len - np.arange(1, self.len+1) + 1))]\n",
    "        elif method is \"bh\":\n",
    "            p_times_m_i = self.sorted_pvals * self.len / np.arange(1, self.len+1)\n",
    "            return [np.min([p, p_times_m_i[i+1]]) if i < self.len-1 else p for i, p in enumerate(p_times_m_i)]\n",
    "        elif method is \"lfdr\":\n",
    "            if self.zscores is None:\n",
    "                raise ValueError(\"Z-scores were not provided.\")\n",
    "            return sms.stats.multitest.local_fdr(abs(self.sorted_zscores))\n",
    "        else:\n",
    "            raise ValueError(\"invalid method entered: '{}'\".format(method))\n",
    "            \n",
    "    def adjust_many(self, methods=[\"bonferroni\", \"holm\", \"bh\", \"lfdr\"]):\n",
    "        if self.zscores is not None:\n",
    "            df = pd.DataFrame(np.c_[self.sorted_pvals, self.sorted_zscores], columns=[\"p_values\", \"z_scores\"])\n",
    "            for method in methods:\n",
    "                df[method] = self.adjust(method)\n",
    "        else:\n",
    "            df = pd.DataFrame(self.sorted_pvals, columns=[\"p_values\"])\n",
    "            for method in methods:\n",
    "                if method is not \"lfdr\":\n",
    "                    df[method] = self.adjust(method)\n",
    "        return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "level of false alarm = 0.0000000000\n"
     ]
    }
   ],
   "source": [
    "df_pval = results.sort_values(['p'])\n",
    "pvals = df_pval['p'].values\n",
    "converter = MCPConverter(pvals = pvals)\n",
    "d = converter.adjust_many()\n",
    "df_pval['p_corrected'] = d['bonferroni'].values\n",
    "level = np.sum(df_pval['p_corrected'] < 0.05) / df_pval.shape[0]\n",
    "print(f'level of false alarm = {level:.10f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2a5221dfba8>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAEXCAYAAAAqfto4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8FVX+//FXkE4ERYIooCJyDy0IhKJIC01QAghio1hQbOD+XL+sYEPYyAqPVVFWXAsKsv5+2IC1t0WaC0IQKSIH5WuhF6kBQkgyvz/m3kvKTb8hN5P38/HIY2DmMzNnzp17P1POnIlyHAcREREvqFDaBRAREQkXJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfEMJTUREfGMiqVdAK9Ys2bNWqARkAz8XMrFEREpKy4DooFf4uLi2hR3YUpq4dMIqOX/q1/KZRERKWsahWMhSmrhkwzUqlChAtWrVy/8zMnJAERHR4e5WN6k+ioc1VfhqL4Kpzj1dfz4cTIyMsD9DS02JbXw+RmoX716dYwxhZ55zZo1AEWatzxSfRWO6qtwVF+FU5z6stYGkmJYbtsUOakZY24D3gC6WGuXh5juAyYBnYHzcAv8CjDTWpsRIv5CYCLQG7gA+B34FzDNWnsyRPw5wATgOqAhsAd4H5hkrT0SIr4a8H+A4binuYeAj4EnrLW7Crn5IiISgYrU+tEYcyUwI4/plwOrgZuA34DPcBPPDODNEPENgG+B0ZxONjWBycBnxphK2eJrAkuAvwAZwEf+4Z+BFcaYWtniKwH/BqYAZwOfAAeAO4HvjDEXFaoCREQkIhU6qRljBgOf47ZWCTU9Cjdx1QRGWGs7W2sHAz5gPTDMGDMk22wzgQbA49battba63FbxHwFdAceyBafCLQCXgWaW2uH+pc/F2jun57ZWNwzwI+BJtba6621LXGTXD3gxUJVgoiIRKQCJzVjTANjzJu4l/jOwr3cF0pv3ISz2Fr7r8BIa+0+4D7/f4NJyrgXYfsDW3GTTCD+GDAKSMdNSoH4c3DPsI4ADwUuZVpr0/zLPwiMMsbU8MdH4Z7BOcCYbJcyHwcs0N8YE5aWNyIiUnoKc6aWCIwAkoArgM25xPX1Dxdmn2Ct/QbYC3Q2xpztH301EAV8mP1em7X2d+A74GJjTHP/6K5ANWCRtfZotvhk3LO7akA3/+hY3Cb266y1v2aLzwA+8P+3Xy7bIyIiZURhktpm4Fago7V2Qx5xLfzDjblMt/71BpJUfvGB5Bl7huJFRKSMKnDrR2vt0wUMvcA/zK1FYWD8+REaXyzJycnB5q1FUZx5yyPVV+GovgqnrNZXXFxcqawvEuqrJPp+rOEfHs9l+gn/MNDQJNLiRUSkjCqJh68D98WcXKZHZRtGWnyxREdH5/sAouM47N27lwMHDgTHpaenExUVRYUK6mO6IPw9EKi+Ckj1VTh51Vft2rWpW7cuUVFh+ckoUWeqiI7/17UoZ4iZHr4Oi5JIaoHSVctlelX/8FiExpe4I0eOcODAAaKiojjvvPM4++yzSUtLIyoqiho1auS/AOHYMffjUn0VjOqrcDLXV0ZGBqmpqRw9epQ//viDAwcOULVqVWrVqpXPUqQ0lMRh207/sF4u07Pf44q0+BJ38OBBAM4//3xiYmKoWrVqmTjqEymPKlSoQNWqVYmJieH8891b74HvsESekkhqgVaGzbNP8D8z1hT32bNN+cX7NfMPAy0uSzq+xKWkpABQs2bNM7VKEQmDwHc28B2WyFMSSe0z/3BQiGmdgBhgeaZnzALxA4wxWcrj776qDfCbtTaQBJfiNu7oFXjAOlN8NNAL95LjMgBr7Y+4XXW1NcY0zBZfARiAe7/t80JuZ5E5/gvQZ5111plapYiEQeAeW+A7LJGnJJLaEuAHoLcx5q7ASGNMDG53WADPBMZba3/BTWwGt6/HQHwN4DXc3ksyxx8D5gDnAjONMRX98RVxu7s6B3gl24PZL/mXMytbIvwrbvdaC6y1W4u32SLidbpNEPnC3lDEWpthjLkD+A/wijFmFO59re64iehVa+2H2Wa7H/gGeNQYMxD3Ae1OuPe7PsVNSpk9CsQDI3F7J/kOaAtcCqzF7e0/s+dwu+LqDfxkjPkvbhJtifs2gDHF3GwREYkAJdK+11q7CuiI209kE6AP7iXAe4B7Q8T/L9ABmI17efJa3D4cJwCD/f06Zo4/gJv0XgAqAQm4TfenAfH+7rIyx6fidsf1V9zn1frjPpf2CnCFXj0jIuINRT5Ts9Z2z2f6JuD6QixvG3B7IeIPAH/y/xUk/jjwhP9PRMoQx3F06U8KRE9iSkSaMWMGxhhmzpyZf3ARpaam8re//Y0uXbrQsmVLOnfuzMaNuXURKqXh6NGjJCYm8sEHH+QfXEQzZ87EGMOMGbm+IlLKECW1SBcVVTb+yqBXXnmF2bNnc+zYMbp3705sbCwNGzbMf0Y5Y6ZNm8bcuXNJT08v7aJIGVESPYqIlAnr168H4IknnmDQoFBPoEhpC3RXJVJQOlOTcis1NRWACy64IJ9IESkrlNQk4n300UcMGDCA2NhY4uPjefrppzly5EjI2HXr1nH//fdzxRVXEBsby9VXX81zzz2XpcPU+fPnY4xhxYoVAIwcOTLHPZWDBw8ybdo0rr76alq2bEmHDh0YNWoUy5Yty7HOESNGYIzhp59+YtiwYbRs2ZKuXbtmiS1IuQpiz549TJkyhd69e9OqVSt69uzJhAkT2L59e47YnTt3MnHiRHr06EGHDh3o2bMnY8eODZ6hZtajRw86duzIjz/+yMCBA2nZsiU9e/Zk8+bNwfubn332GRMmTKB169Z07NiRl146/aTN7t27mThxIvHx8cH7k+PHj2fbtm0ht+PIkSM8//zzXHPNNVx++eV069aNBx54AGttMMYYw3vvvQfAhAkTMMbw7bffBqenpKTwz3/+k4SEBFq1akX79u258847WbVqVch1Hj58mKlTp9KjRw9atWrFwIED+eSTTwpW8VJmKKlJRPvwww956KGHyMjIID4+noyMDN544w1uu+02Dh8+nCV2/vz53HzzzSxatIiGDRsSHx/PyZMn+ec//8nNN9/MoUOHALjoootISEggJiYGgE6dOpGQkBB8u8Lvv//OwIEDmTVrFikpKfTo0SOYBO+8807+8Y9/hCzrmDFj2LZtG927d6dChQq0aNGiUOXKz+bNmxk8eDBz5syhQoUKdO/enRo1ajB//nyGDBnCb7/9Foxdt24dAwYMYN68eVSqVInu3btTv359vvjiC2666aZgssgsNTWV0aNHk5KSQteuXalYsSKNGzcOTn/uuef49NNP6dSpE3Xq1KFJkyYAbNq0iUGDBjFv3jyqVKlCfHw8MTExLFiwgMGDB+dIort27WLo0KHMnDmT5ORkunXrRr169fj888+5/vrrWbt2LQAJCQlcdNFFALRp04aEhATq1KkDuEnxlltu4bnnnuPAgQN06tSJZs2asWLFCkaOHMm8efOyrPPgwYMMGzaM119/Hcdx6N69O47j8OCDD/Lxxx8XqP6ljHAcR39h+EtKSlqclJTkbN682cnPpk2bnE2bNmUZl5yc7CQnJ+cMdt/qEPl/YfbCCy84Pp/P8fl8zvTp052MjAzHcRwnJSXFue+++xyfz+dMmDAhGP/zzz87LVq0cOLi4pykpKTg+NTUVOfxxx93fD6f8+c//znLOm699VbH5/M5K1euDI7LyMhwrrvuOsfn8zmTJk1yUlNTg9PWrVvndOjQwfH5fM6SJUuC44cPH+74fD6nR48ezpEjRxzHcZz09PQilyuU9PR0Z8CAAY7P53NmzJgRrA/HcZwZM2Y4Pp/PGTVqlOM4jnPixAmnc+fOjs/nc15++WUnIyMjuH8tXrzYiY2NdVq0aJFlX42Pj3d8Pp8zdOjQ4DYHtiHwWTRr1sz54YcfspTp5MmTTs+ePR2fz+fMnTs3S5kXLFjgGGOc+Ph45+TJk8Hxd999t+Pz+ZxHH300S/3Onz/f8fl8Tt++fYPjHnnkEcfn8znvv/9+lmX/z//8j+Pz+Zxx48Y5J06cCI7/4YcfnI4dOzotWrRwtmzZEhw/adIkx+fzOX/605+ylOXll18O7mcvvPBCcHyu30cn9Pc3UpWFr//mzZudpKQkJykpabETht9inalJRLv00ksZO3Zs8BmlKlWqkJiYSOXKlfn44485ftx99+ucOXM4deoUDzzwQJZ3OlWqVInHHnuM888/n08++YQ9e/bkub7Vq1fzww8/0LhxYx599FEqVaoUnNaqVSvGjx8PwKxZs3LMe91113H22WcDp/sIDFe51q5dy+bNm2nZsiVjxozJ8szWvffeS9OmTUlLSyM1NZVPP/2UvXv30rlzZ0aPHp0ltlu3bowePZpTp04xZ86cHOu5+eabg9uc/V1icXFxNG9+ul/wChUq8OWXX7Jt2zZ69+7N8OHDs8QPGjSIPn36sGPHDr744gvAvXz69ddfExMTwxNPPJGlfq+77jo6d+5MrVq12L9/f651sWfPHj7++GPq1q3L5MmTqVq1anBa8+bNGTt2LKdOnWLu3LmAewa6YMECqlatyuTJk6lcuXIwfvTo0Vx++eW5rkvKHiU1iWj9+vXL8eN67rnn0qpVK1JSUoLPlQXutXTs2DHHMipXrkyHDh3IyMggKSkpz/WtXr0agD59+oTscLpv376cddZZfPfddzmamTdt2jRHfLjKFbhPFB8fn2PaWWedxb///W9mz55N5cqVg9vQt2/fkMu65pprgNPbmt825DUtr+0D6NKlS5byB+K7dOmSJbkEzJo1i3nz5gUvM4ayevVq0tPTad26dZaEFtC5c+cs69ywYQPHjx+ndevWId+M0bNnz1zXJWWPmvRLRKtfv37I8XXr1gVg7969gNtQAWDAgAF5Lm/Xrrx7RAssL7f1VqtWjdq1a7Nv3z4OHz5M7dq1g9NCvTQyXOXat28fAPXq5fZawNMC29CgQYOQ0wPjQ50N5fXiy1DTAuVOTEwkMTEx13kD9VCY7chNYJ1ffPFFnm+ZD6wzUB+Bd6Fll9tnLWWTkppEtFBH4nD61R8VK7q7cOCsqX///nl2p3TxxRfnub7AcvNaRuDZqexnGqHmCVe50tLS8pyeWWAbcpNb+SHnJcf8pgWW1alTJ84777xc573ssssAwvIQdWCdPp8vz6RW0G619Aoob1FSk4gWOMrOLnC0Hjjir1u3Ljt27GDcuHHFOgsInAHm1hQ9OTmZgwcPUrVqVaKjowu0vHCUK9BSM3D2kd3ixYs5fvw4nTp1Cm5DqGb+cHrb8kpChS3XoEGDGDhwYL7xgcuKuW3H6tWr2bVrFx07dsz1zCqwzlatWvHUU0/lu85Ave/cuTPk9MDZo3iD7qlJRPvmm29yjNuzZw8bN24kOjo62HChXbt2ACxZsiTkckaNGsWNN94Y8hmtzNq3bw/Al19+GfKs4vPPPycjI4MOHToUqPzhKlfbtm0BWL58eY5pjuOQmJjIQw89RFpaWnAbPvvssxyxAJ9++ilAgbchL4HtW7p0acjp06dPZ+DAgbzzzjvA6e1YsWJFyLPPF154gXHjxgUbzoQ62wqsc+XKlZw8eTLH9CVLltC3b1+efPJJAFq2bEnNmjX5/vvvQ15yze2zkbJJSU0i2rJly3j77beD/09OTuYvf/kLaWlpXH/99cFLaCNGjKBChQo8++yzWRpdOI7DP/7xD5YvX8727dvzbAgB7g998+bN2bp1K0899RSnTp0KTtu4cSPTpk0DYNiwYQUqf7jKdeWVV9KoUSPWrl3L66+/nmXazJkz2bZtW/D5sX79+lG3bl2WL1/OK6+8kuVy5NKlS3nttdeoVKkSN954Y4G2IS/XXnstMTExfPTRR7z11ltZpi1btoxZs2ZhrSU2NhaASy65hKuuuopdu3YxderULAcOCxcuZNWqVTRq1CgYX6VKFcDt2DjgoosuIj4+nu3btzNx4kROnDgRnLZ9+3YmTZrEL7/8QqNGjQC3pektt9zCqVOnePjhh4MtZgHeeeedkAcKUnbp8qNEtNatW/PEE0/wzjvvUL9+fZKSkvjjjz9o06YNd999dzAuNjaWhx9+mKeffprhw4fTvHlz6tevz5YtW/j111+pWrUqzz//fMj7SJlFRUXx7LPPcuutt/LWW2+xaNEiWrVqxcGDB1mzZg3p6encd999dO/evUDlD1e5AonxtttuY+rUqcyfP59LL72UrVu38vPPP1OnTh2mTJkCuI1Znn/+eUaPHs0zzzzD+++/T5MmTdi9ezcbNmygYsWKTJw4kWbNmhVoG/JSrVo1pk+fzt13383kyZOZM2cOTZo0Yf/+/Xz//fcAjB8/Psu6EhMTGTZsGG+++SaLFi2iRYsW7Nixg40bNwaXFzhDC9xrfPHFF1mzZg233norcXFxJCYmMmLECBYsWMCSJUuIjY0lPT2dVatWkZqamuMRg3vvvZekpCSWL19O7969iYuLY+fOnWzYsIHWrVsHyypln87UIl3pP1ZdsL8SMnLkSCZNmkRycjKLFi2iRo0ajB07lpkzZwaP4gNuu+023nzzTeLj49m5cyeLFy8mIyOD6667joULFwYvW+WnUaNGLFiwgNtvv51KlSqxaNEitm7dSteuXZk9ezZ/+lOBXuEX9nI1b96c+fPnM3ToUI4ePcqiRYs4fPgwgwcP5r333styD6pt27YsWLCAG264gZMnT7J48WJ2795N//79mTdvHjfccEOhtiEv7dq1Y+HChQwdOpTU1FSWLFnCzp07g/V1++1ZX5N44YUX8v7773PbbbcBsGjRIrZv307fvn159913s5y13nDDDQwYMIC0tDSWLVvGTz/9BLj35t59913GjBnDeeedx8qVK9m4cSPNmjUjMTGR6dOnZ2kAUrVqVWbNmsWf/vQnoqOj+frrrzl69ChPPvkkt9xyS9jqQkpfVH4tpaRg1qxZsxjoFh0dnWeLLIAff/wRIMvR67FjxwCoUaNGSRXRU1RfhaP6Kpy86ivU9zdSnam3QhUnjVhrA32gLomLi+te3LLoTE1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU1ERDxDSU3KDb2R4sxSfUtpUFITz9i9ezdjxoyhQ4cOtGrVij59+pCRkcGvv/7KqFGj2LFjR2kXsdzYsGFDWN/ZFkrv3r0xxrB9+/YSXY+ULXrzdYQ7U+9DKq5IOCgfN24cq1aton79+lx55ZXUqlWLChUqcPfdd/Prr7+WdvHKlZtvvplTp06VdjGkHFJSE8/YsGEDAG+99RYXXHBBcHxGRkZpFancUp1LadHlR/GM1NRUgCwJTUTKFyU1iViHDh1i+vTpDBo0iLZt29KyZUu6du3KuHHjslxOHD9+PMYY0tPTATDGYIzh22+/xRjD77//DkDPnj0xxmRZx+7du5k4cSLx8fG0bNmSzp07M378eLZt25ajPMYYBg8ezIoVK+jTpw+xsbH069eP/fv357st//3vf7nnnnvo1KkTbdq0YeDAgcyZMyeYiDNbsmQJo0aNon379sTGxnL11Vfz97//ncOHD2eJC2zf1KlTeeONN7jiiito3bo199xzDwA9evSgY8eO/Pjjj9x000107NiRnj17snnz5uAyli5dyh133EH79u1p1aoVCQkJzJo1K2S5ANavX8+DDz5Ily5daN26Nddeey0zZswgOTkZgPnz5+f4LHr06JFlGVu3bmXcuHF07tyZli1b0qNHDxITE/njjz9CrnPdunXcfffddOzYkbi4OMaMGRPy8xEBJTWJUPv372fIkCG89NJLHD9+nE6dOtGxY0dOnjzJBx98wMiRI9m9ezcAbdq0ISEhgSj/DciEhAQSEhKoU6cOCQkJVK9eHYBevXqRkJAQXMemTZsYNGgQ8+bNo0qVKsTHxxMTE8OCBQsYPHgw69evz1GuvXv3ct9991GtWjWuuuoqatasSZ06dfLclpdffpk77riDpUuX0rhxYzp16sSePXuYMmUKDz74YJZLdX//+98ZPXo0K1asoGnTpsTHx3PixAleffVVBg8eHLJRxKJFi5g6dSrNmjWjZcuWXHzxxcFpqampjB49mpSUFK666ioqVqxI48aNAZg5cyZ33XUXq1atokmTJnTt2pX9+/czbdo07rzzzhyJbeHChdxyyy188sknXHjhhXTp0oXk5GT+8Y9/cMcdd3Dy5EkuuuiiHJ9Fr169gstYvnw5Q4YM4YMPPuC8886jR48eVK5cmblz5zJkyJAcyWrRokUMGzaMxYsX07hxY6688kpWr17NTTfdxJEjR/KsdymnHMfRXxj+kpKSFiclJTmbN2928rNp0yZn06ZNWcYlJyc7ycnJOWLdJhiR/xduTz75pOPz+ZwpU6Y4GRkZwfFHjx51brrpJsfn8znTp0/PMk+zZs0cn8+XY1m9evVyfD6fs23btuC4kydPOj179nR8Pp8zd+7cLPELFixwjDFOfHy8c/LkyeB4n8/n+Hw+Z8yYMcEypaen57kd69evd5o2bep06NDB2bBhQ3D8kSNHnIEDBzo+n8/5+OOPHcdxnK+++srx+XzOlVdemWX/OHnypPPoo486Pp/PueGGG4LjV65cGSzTnDlzguMDZYqPj3d8Pp8zdOhQ5+DBg05ycnJw2jfffOP4fD6ne/fuzpYtW4LzHjt2zLnnnnscn8/nPPPMM8HxO3fudFq3bu20aNHCWbJkSXB8SkqKM2rUKMfn8zkvv/xycHyoz+KPP/5w2rdv7zRv3tz54osvguMzMjKcmTNnOj6fz7nxxhuD45OTk51OnTo5TZs2dT7//PPg+MOHDzs33nhjcNszf67hktv30XFCf38jVVn4/m/evNlJSkpykpKSFjth+C3WmZpEpHPPPZcuXbowduzY4FE/QHR0NP379wcInqkVxZdffsm2bdvo3bs3w4cPzzJt0KBB9OnThx07dvDFF1/kmHfEiBHBMlWokPdX6O233yYjI4OxY8fSsmXL4Pizzz6bhx56iEaNGrFz504AZs+eDcAjjzxCs2bNgrGVK1fmySef5JJLLuH7778nKSkpyzoqV67MTTfdFPx/9jLdfPPNVKpUKcu0WbNmAfDYY4/RpEmTYGz16tV56qmnqFq1Km+99VbwbG3hwoUcP36cESNG0LVr12B8lSpVGD9+PA0bNsz18mHAe++9x+HDhxk+fDi9e/cOjo+KiuLee+8lNjaWtWvXsnbtWgC++uor9u/fzzXXXEOfPn2C8TVr1uSpp57Kc11SfimpSUR64IEHeO2114iOjg6OO3DgACtWrAj+qBenyfi3334LQMeOHUNO79KlCwCrVq3KMa1p06YFXk9g/vj4+JDr+Oyzz7jzzjtJS0tj7dq1VKxYMcsPfkDFihWDP+zZy3TppZdSuXLlXMuQvbzp6enBOgy1/bVr16Z58+YkJyezadOmfLfjsssu46uvvmLChAm5lgHyr/OrrroKgNWrV2cZBj6LzBo3bswll1yS5/qkfFKTfolYv//+O//617/47rvv+OWXX4KNEaLC8PDerl27AEhMTCQxMTHXuOxngxUqVKBmzZoFXs++ffsAqFevXp5xhw4d4tSpU9SrV48qVaqEjGnQoAFAjoYptWrVynPZ2acfOnSIlJQUAOLi4vKcd9euXbRu3Tq4HcVpWRqo83vvvbdAcXv37gXg/PPPDxnXoEEDPX8oOSipSUT68MMPefjhh0lPT+eSSy6ha9euNG7cmNjYWHbv3s0TTzxRrOUHGmd06tSJ8847L9e4yy67LMv/C5tQ09LSChTnOE6+yw+0KMx+VpbfJdDs0wPLqVatWpZGHKHExMQABd+OvATW27Nnz2DjnVAKeiZ81llnFbtM4j1KahJxjh07xsSJE6lQoQIvvfQS3bp1yzJ97ty5xV5H4Md60KBBDBw4sNjLy2s9O3bsYM+ePVx44YVZpqWlpfH222/TqFEj2rdvT6VKldi3bx8nT54MebYWaBmYVxIuiHPOOYdKlSqRlpbG1KlTC5QcYmJi+OWXX9i9ezcNGzbMMf3dd9+ldu3a9OjRI9fEXLduXX799VfuuOMO2rVrl+86A2e3gXuO2QXOHkUy0z01iThbt27l2LFjtGjRIkdCA/jmm2+AgvdaEepHNvCjunTp0pDzTJ8+nYEDB/LOO+8UtNghtWnTJtf1fPfdd0yePJnZs2dTqVIl2rRpQ1paGl9++WWO2LS0NL766isg93tSBVW5cmUuv/xyTp06xYoVK3JMT01NZfDgwdxyyy3BRwjatm2b63bs2LGDxx57jGeffTZY10Wp8wkTJjBkyBD+85//AHDllVcCBLc7s927d2OtzXdbpfxRUpOIEzhC37JlS5bnltLT03nppZf4+uuvAXJ9QDi7wFlP4J4cwLXXXktMTAwfffQRb731Vpb4ZcuWMWvWLKy1xMbGFmtbbr75ZqKiopgxYwY///xzcPzhw4d5+umnARgwYAAAt956KwBTpkzhxx9/DMaeOnWKSZMm8fvvvxMbG0urVq2KVabM65o4cSJbtmwJjk9LS+Ovf/0rP/zwA8ePHw/ex7v++uuDz5MFGnAApKSkMHny5CzbAacvkR49ejQ47sYbb6RatWq8/vrrORL3/PnzWbBgAVu2bOHyyy8H3IfHGzZsyKJFi7IcXBw/fpxHHnkkeDlTJDNdfpSIU7duXa655ho++eQTEhIS6NChAxUrVmT9+vXs27ePyy67jJ9//rlAPXkAXHzxxWzZsoUHHniApk2bMmXKFKKjo5k+fTp33303kydPZs6cOTRp0oT9+/fz/fffA25PJZmb1hdFu3btGDNmDDNmzGDQoEF06NCBypUrs3btWg4dOkRCQkLwEYVevXpxxx138PrrrzNkyBDatWvHOeecw7p169i9ezcNGjTg2WefLVZ5Avr06cN+mfr2AAAXh0lEQVStt97KnDlzGDx4MC1btqROnTps3LiRXbt2Ubt27SzratiwIRMnTuTxxx9n5MiRxMXFUatWLdavX8/evXvp0KEDo0aNCsZfcsklbNq0ieHDh3PZZZfxzDPPUK9ePf72t78xbtw4xowZg8/n45JLLuG3337DWkuFChWYNm1a8GH2KlWqMG3aNO666y4ef/xx3n77berXr8+aNWtITU2lUaNG/PLLL2GpD/EOnalFuNJ/rLpgf+E2ZcoUxowZQ7169Vi5ciXr16+nQYMGPPHEEyxYsICaNWuyYcMGDhw4kO+y/vKXvxAXF8fu3btZuXJl8JJau3btWLhwIUOHDiU1NZUlS5awc+dOunbtyuzZs7n99tvDsi1jxozhpZdeIi4ujnXr1rF8+XLq1q3LhAkTmDp1apbYhx9+mBdffJEOHTqwadMmlixZQnR0NPfffz8LFizgoosuCkuZwH0e7sUXX6R9+/Zs3bqVpUuXUrVqVUaMGMHChQu59NJLs8Rff/31zJ07l27duvHTTz+xZMkSqlWrxv3338+rr75KxYqnj5EnTZpEs2bN2Lp1K//97385dOgQAP369ePdd9+lf//+HDx4kK+//pqjR49y9dVX8+6779KvX78s62zbti3vvPMO11xzDbt27WLZsmX4fD7mzp2bb4tSKZ+inJL4RSqH1qxZsxjoFh0dnaN/wewCl5YynwUcO3YMgBo1apRUET1F9VU4qq/Cyau+Qn1/I9WZenVVcdKItTZwa2BJXFxc9+KWRWdqIiLiGUpqIiLiGUpqIiLiGUpqIiLiGUpqIiLiGUpqIiLiGUpqIiIFpEegIp+SWikI9Iunbn5EypZAf6PheP2RlAwltVJQtWpVAI4cOVLKJRGRwgh8ZwPfYYk86vuxFJx77rmcOHGCPXv2kJaWxtlnn63LGiIRyHEcHMchNTWVo0ePBvsbPffcc0u5ZJIbJbVSULNmTVJSUjhw4AD79+9n//79pKenExUVle8LH8UVuAyk+ioY1Vfh5FVftWvXLtTbz+XMUlIrBVFRUZx//vnUqFGDI0eOcPz4cY4fPw6Q5xuB5bSUlBRA9VVQqq/CyVxfUVFRVKxYkerVq1OzZk2io6NLuXSSFyW1UhQdHR38gqxZswYoG52kRgLVV+GovgpH9VV26VqEiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4hpKaiIh4RsWSWrAxZgTwZh4hT1lrH8sU3w6YCLQHooEfgOettf83l+X7gElAZ+A84GfgFWCmtTYjRPyF/uX3Bi4Afgf+BUyz1p4s9AaKiEjEKbGkBrTxD78E9oaY/n3gH8aY3sDHuGeOS4DjQE/gLWNMC2vto5lnNMZcDiwFagLfAKuBeGAGcAUwPFt8A2AF0ABYC3wHXAVMBnoYY/pYa08VZ2NFRKT0nYmkdru1dkduQcaYarhnTAC9rbVf+8c3BhYDjxhj5ltr1/jHR+GeAdYERlhr/+UfHwN8BQwzxiyw1r6faTUzcRPa49baRH98DWAh0At4AHim+JssIiKlqSTvqbUG9uSV0PxGAHWBtwIJDcBauxUY7//vA5niewOtgMWBhOaP3wfclz3eGGOA/sBWYEqm+GPAKCAdGFuoLRMRkYhUIknNGNMIOAdYU4Dwvv7hwhDTPsRNOv0KEm+t/Qb3UmdnY8zZ/tFXA1HAh9nvtVlrf8e9FHmxMaZ5AcoqIiIRrKTO1AKXHvcYY2YYY342xqQYY6wx5nFjTNVMsS38w43ZF2KtPQLsBGKMMefnFx+YDXe7Akkqv/jN/mFs7psjIiJlQUkntduBYbgtGb/Fva81GfiP/14auC0RAXblsqzA+EBSK+l4EREpo0qqoUggqb0D3OG/f4Ux5hLcy4adgETgIaCGP/ZELssKjI/2DwPxx0sovliSk5NZs6YgV11DK8685ZHqq3BUX4VTVusrLi6uVNYbCfVVUmdq1+Ne9hsRSGgA1tpfgdsABxhtjKmEe8/MsdY6uSwrKtswcF+spOJFRKSMKpEzNWttCrApl2nfG2O2Aw0BH3AMOMcYU9U/X3aB+2+B5JjsH1YLERuO+GKJjo7GbXBZOIEjnNI6wiprVF+Fo/oqHNVX0RSlvqy1JCcn5x9YQKXVTdZu/7A6bkMQgHq5xGa/J1bS8SIiUkaFPakZY842xrxijHnPGJPbmWAj/3AHp1sl5mhSb4ypCVwI7LPW7vGPzis+CmiKe0lzU37xfs38ww25TBcRkTKiJM7UkoHrgCFAt+wTjTF9gTrABmvtTuAz/6RBIZaVAJwFfJJpXF7xnYAYYLm19mi2+AHGmCzba4y5CLdRy2/W2pCXS0VEpOwIe1LzN/h41f/fGf6OhIFg11cz/f9N9A/fx31g+jZjzDWZYi8FnsZt4PFsplUswX1EoLcx5q5M8TGZlh3s8spa+wtuYjO4jxME4msAr+EmTXWRJSLiASXVpP+vQBfcHvStMWa5f3w8UAV41lr7DrgPWPuT0/vAR8aYJcBR3A6NqwOPWmvXBxZsrc0wxtwB/Ad4xRgzCve+WXfgXOBVa+2H2cpzP27Hx48aYwbiPqDdCfd+2qfAS2HefhERKQUl0lDEWnsCNymNB37FTWadgJXAEGvtQ9niP8C9VPkF7uXAbsB64AZr7RSysdauAjriJsImQB/gN+Ae4N4Q8f8LdABm416evBY4CEwABltr04q5ySIiEgFKrJd+a20qMNX/V5D4/3K6X8eCxG/CfR6uoPHbcHs4ERERj9Kbr0VExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDOU1ERExDMqlnYBRIoiLi6utIsgIhFISS1C6EdaSpL2LykvlNSk9EVFndn1Oc6ZXZ+UOToIKLuU1CJNaf3An+n1Zl63nDk6gBCPU1ITkZJXmgdNOlAsV5TUpPzSj52I56hJv4iIeIaSmoiIeIaSmoiIeIaSmoiIeIaSmoiIeIaSmoiIeIaSmoiIeIaSmoiIeIaSmoiIeIaSmoiIeEa56ybLGNMLeARoBVQG1gBPW2s/L9WCiYhIsZWrMzVjzG3Al0AnYBWwArgK+MwYM7oUiyYiImFQbpKaMeYC4J/AYaCdtfYaa+3VuEntCPC8MaZ+aZZRRESKp9wkNWAsUAV4zlq7MTDSWrsamAZUBXS2JiJShpWnpNbXP1wYYtoC/7DfGSqLiIiUgHKR1IwxUUBzIAP4MUTIFv+0Fv5YEREpg6KccvDyQGNMbeAPYJ+1tm4uMXuAukAta+2Rwq5jzZo124Ei35OLjo4u6qwiIhEhOTm5OLPviIuLa1DcMpSXJv01/MPjecSc8A+jcRuOFFaxslIxdwYRkbIuLEf25SWpZfiHeZ2WRmUbFtYvQCMgGfi5iMsQESlvLsNNaL+EY2HlJakFToOq5RFT1T88VpQVxMXFtSnKfCIiEj7loqEI7uXEZKCOMSZHIvePqwOkWGsPnenCiYhIeJSLpGatdYBNwFmAL0SIwa2LDWeyXCIiEl7lIqn5feYfDgoxLTDukzNUFhERKQHlKam9AaQADxtj4gIjjTHtgL/gtn6cWUplExGRMCgXz6kFGGPuA14ETgH/wW3p2AO3wcxIa+2/SrF4IiJSTOUqqQEYY/rjnpm1BU4C64CnrLX/KdWCiYhIsZW7pCYiIt5Vnu6piYiIxympiYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZ5SXXvpLlDGmF/AI0AqoDKwBnrbWfl6IZfiASUBn4Dzc19e8Asy01maEiL8QmAj0Bi4Afgf+BUyz1p4s1gaVsDDVVz/g/wDtcV9bsQv4FEi01m7PFlsRt0PrKrksboe1ttgvJywpxa0vY0xD3P0jN99Yaztnm6fM7l9QvDozxiwGuhVgNZOstU9mmm8rcGke8ZWstWkFWG6pMcbchtv7Uhdr7fJCzFeo/cUYcw4wAbgOaAjsAd7HrdOivM8ySM+pFVOmneAksAi30+R4oBJwt7X2lQIs43JgKVAT+AbY61/GOcBb1trh2eIbACuABsBa4H+Bq4B6wGKgj7X2VPG3LvzCVF/jgb/hvidvFe4Xog1wEbAP6Gqt3ZwpvhXuQ/ZbgZUhFnnAWvtA0beq5ISpvgYA/wbWE7rTbmut/Wum+DK7f0Hx68wY8yjQLJfJ0cBA/7+vs9Yu9M9TCziI+939Kpd5R4Y6QI0UxpgrgS9wt7HASa2w+4sxpiawDPeAw+Luk3G476PcBHSy1h4u6nboTK0YjDEXAP8EDgOdrbUb/ePb4+7YzxtjPrbW7shjGVHAm7gJbUSgqy5jTIx/GcOMMQuste9nmm0m7g70uLU20R9fA1gI9AIeAJ4J68aGQZjqqzmQiHvm1cdau8I/vhIwHbgP9wftykyzBd5194a19qnwblXJCUd9+QW2f5q19q0CrLpM7l8QnjrLax8xxrzp/+ezgYTm1xq3270PrbV3FXMzzjhjzGBgNkV7+3Rh95dE3IT2KnCPtTbDfzXldWCEf/rYom2J7qkV11jcS1rPBb48ANba1cA03BePjs5nGb1xP+DFmfuetNbuw/2BBnenAMAYY4D+uGcdUzLFHwNGAekUY4coYeGorxG4R97PBhKafxmncC9H7gOuMMZcnGmewI/6mmJvwZkVjvqCQmx/Gd+/IHx1loMx5hbc/W8j7qWzzMrkPmaMaeBP1O/jfq/2FHL+Qu0v/suOd+K+4/KhwJmr/7Lsfbhnu6P8SbFIlNSKp69/uDDEtAX+Yb+iLsNaG7gU2dkYc7Z/9NWcPiLMyBb/O/AdcLH/jCbShKO+UnEvoy3NPsGf2AKvhL8w06TAD853BStmxAhHfYG7/cnAlgLEluX9C8JXZ1kYY6I5fbZxj7U2NVtImUxquGdFI4Ak4Apgc97hORR2f+kKVAMWWWuPZotPxj2brkbB7mmGpKRWRP7Lhs1x7+v8GCJki39aC39sblr4hxtzmW5xP6fATpFffGCnjM1jnWdcuOrLWjvRWnt5qA6o/Ud3gXranmm9rYHdwABjzCpjzFFjzD5jzP/zH2lGnHDVlzGmNu69xi3An40x64wxx40xO40xr/hv8GdWJvcvCOt3MpRHce8Rve0/2MyuDe5ZiTHGLDXGHDLGHDTGfGiM6VDIdZ1Jm4FbgY7W2qK8JLmw+0uJ719KakV3Lu5ljj9CHLUFTqf3A9WBs7NPz+QC/3BXLtMD488vYnykCFd95eVh3HsCq6212/zjLsW9X1kPeBn3nXpf+4c3AauNMVcVcX0lKVz1FTiDaIt7eWgv7vZXBO4C1mRL7GV1/4IS2sf8BwYPAA4wOcT0KrgNS84C5uJe4vzav67+wHJjzPWF3ZgzwVr7tLX2zWI0YIm43y8ltaILXPM9nkfMCf8wr5uv+S0n+zIKGx8pwlVfIRljrsFtwp2B+2qhgMCP+g4gzlrb1Vo7ALel1TO4P25vG2OqFnadJSxc9RXY/h8AY63tba29Fnf7/x9uss/ceKSs7l9QcvvYvbiJ8ENr7aYQ02NxDxKOAr2stR2stddZa5sAD+K2upxtjKlXiHWWFRH3+6WkVnSBI5u8nomIyjYsynKyL6Ow8ZEiXPWVgzHmWk7f6H7EWrs40+T3cS+/dbDWfh8Y6T9q/wvuPZD6wKDCrPMMCFd9PYd7ttrdWhu43xi4kX8n/mRvjLmigOuN1P0LSmAfM8acBdzv/++0UDHW2iTcM5DY7JfFrbXTce/v1QBuK8g6y5iI+/1Sk/6iS/YPq+UREzj6P1aM5WRfRmHjI0W46isLY8wduJcVKwKTrbVTM0+31jrAtlDz+psSf4L7jEwcMK+g6z0DwlJf1tp0TjeeyT7tuDFmEW5DgTjcZ/jK6v4FJbOPdcVNWL/kci8NAGvt7jyW8SHuQVNcAddZlkTc75fO1IruCO4HVMf/jEUW/nF1gBRr7aE8lrPTP8zt0kT2a9CFjY8U4aqvzPP8FZiFe4b2oLV2YhHKFfgxql6EeUtS2OsrF9m3v6zuX1AydTbYP3y7GOWK1H0sHCLu90tJrYj8ZwCbcH9QfSFCDG795teiKNAKKEcTaX8Lraa4rao25RfvF+gJoSgtmUpMGOsLY0yUMeY14DHcXiNu8l/mCRV7vzHmbX+3SaE08g+35zK9VISrvowxE40x7xljcmtNln37y+T+BeHdxzK5xj+cn1uAMeZGY8z/9T/HFkpE7mNhUtj9pcT3LyW14vnMPwx1PyYw7pNiLKMTEAMsz/RMRyB+gDEmy+dnjLkIt2HAb7nc0C5t4agvcBt4jMI9Mr/aWvtOHrGXAjfgNlvOwt84ZKj/v18UYL1nWjjqqxUwBLcOsjDG1AX6AKdwW+tlXmdZ3L8gfPsYxpjzcPef47jdP+WmLnAzboOS7MuIAgLd3BW4b9MypLD7y1LcxiC9sj9g7X8WsBfu2fayohZISa143sBtGv6wMSZ4vdwY0w63EcIJ3C5kAuMbG2Oa+vuJC1iC2zKttzHmrkyxMZnmDXYx47/Z/xnuUefkTPE1gNdwj1IjsgsjwlBfxpi+uC3K0oBrrbVL8lnnLNwz3WHGmCGZllMJmAFcDHxqrY3Eh2bDsX+97B8+lPnRBf8PyOu4jzu8FrgnVMb3LwhPnQW09w/X5tMR8TzcA6zOxpgHMy07CngC96HmDcAHRdukyGCMuchfV3UC4wq7v/gbKM3BffxiZuAysX/4Im5/t69kfzC7MNShcTEZY+7D/TBOAf/BbbXTA7fhwsjMXV8ZY37F/RG93Vo7O9P4Dv55o4Fvca87d8f94F+11mbp1scYcylux8f1cE/nLe5Z3QW4PdUPyOdLWGqKW1/GmJVAR9xWe4vzWNVT1tof/fM8gNsvZBSwGrcH8Y64/dVtBrpZa/eGaRPDKkz71zPAn3Fbnn2D+/xUF9z7S8uAvtba45niy+z+BeGpM/+0e3ET4Fxr7ch81jkE9xGJSrgHqZtxm/r7cO+pdbXW/lT8rStZ5vQbCnJ0aJxpWva3ExRqf/E/9/df3ET4v7i9jrTFPStei1tXgQYlhaYztWKy1s4EEnBbjnXBPbpbDvTO/OXJZxmrcH9k3wea4F4S+g24hxCXNKy1/wt0wO2ANAa4FrfPtAnA4Ej+wSlOfRljqnP66Lk+MCyPv+DDm9baF3D72Pwct377415SegpoH6kJDcK2fz2Ee/nxG9zLQX1xb8T/BeiZOaH548vs/gXhqTO/GP8w33th1u1w/Erce2/nAwNwHwSfAbQqCwmtqAq7v1hrD+AmvRdwDwIScA+4pgHxxUlooDM1ERHxEJ2piYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZyipiYiIZ/x/QFwc3HDqh88AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(df_pval['p'],color = 'red',label = 'before corrected')\n",
    "plt.hist(df_pval['p_corrected'],color = 'blue',label = 'after corrected')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}