p-values.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demonstrate that the p-value is uniformly distributed under $H_0$\n",
    "\n",
    "In this demonstration, we are testing whether the sampled values come from the standard normal distribution (spoiler: they do!). That is:\n",
    "\n",
    "$$H_0: \\text{the sampled values come from } N(\\mu=0,\\sigma=1)$$\n",
    "$$H_A: \\text{the sampled values do not come from } N(\\mu=0,\\sigma=1)$$\n",
    "##### Sample from the standard normal distribution \n",
    "We will start by sampling 1000 times from the standard normal distribution ($\\mu=0, \\sigma=1$) using `scipy.stats.norm.rvs`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAD9FJREFUeJzt3X+sX3V9x/HnS/w5JQFGwQ7Qi0s1otGqHWO6LChzImyiLk5dlM6xVTfwx+ISqizTZGOrMXMbxunqYKJR/DF1dJYgWEHmFPViGoRVRoNVKqStv5CNTFN474/vuePafNr77e0993zv7fORfPP9ns/38z3nfdLmvvI5Pz4nVYUkSft6yNAFSJImkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUtNDhy7gUBx77LE1NTU1dBmStKTcdNNN36uqFXP1W9IBMTU1xfT09NBlSNKSkuTb4/TzEJMkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKlpSd9JLU2yqfWbB9nujg1nD7JdLT+OICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJauotIJKclOS6JNuS3JrkjV37MUmuTXJ79350154klyTZnuTmJM/sqzZJ0tz6HEHsBd5cVU8GTgPOT3IKsB7YUlWrgC3dMsALgVXdax3w3h5rkyTNobdHjlbV3cDd3ed7k2wDTgDOAU7vul0OXA9c2LV/sKoKuDHJUUlWduuRNCYfdaqFsijnIJJMAc8AvgIcP/NHv3s/rut2AnDnrJ/t7NokSQPoPSCSPAb4JPCmqvrxgbo22qqxvnVJppNM79mzZ6HKlCTto9eASPIwRuHw4ar6VNe8K8nK7vuVwO6ufSdw0qyfnwjcte86q2pjVa2pqjUrVqzor3hJOsz1eRVTgEuBbVX1rllfbQLWdp/XAlfOaj+3u5rpNOAezz9I0nB6O0kNPAd4NfCNJFu7trcCG4CPJzkP+A7wsu67q4CzgO3AfcBreqxNkjSHPq9i+iLt8woAZzT6F3B+X/VIkg6Od1JLkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJaurtmdTSJJhav3noEqQlyxGEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ19RYQSS5LsjvJLbPa3p7ku0m2dq+zZn33liTbk9yW5AV91SVJGk+fI4gPAGc22v+2qlZ3r6sAkpwCvAJ4Svebf0hyRI+1SZLm0FtAVNUNwA/G7H4O8NGq+klVfQvYDpzaV22SpLkNcQ7igiQ3d4egju7aTgDunNVnZ9cmSRrIYgfEe4FfBFYDdwN/07Wn0bdaK0iyLsl0kuk9e/b0U6UkaXEDoqp2VdX9VfUA8H4ePIy0EzhpVtcTgbv2s46NVbWmqtasWLGi34Il6TC2qAGRZOWsxZcAM1c4bQJekeQRSU4GVgFfXczaJEk/66F9rTjJFcDpwLFJdgJvA05PsprR4aMdwGsBqurWJB8H/hPYC5xfVff3VZskaW69BURVvbLRfOkB+l8MXNxXPZKkg+Od1JKkJgNCktQ0Z0AkOWYxCpEkTZZxRhBfSfKJJGclad2vIElahsYJiCcCG4FXA9uT/FWSJ/ZbliRpaHMGRI1c212V9AfAWuCrSb6Q5Fd6r1CSNIg5L3NN8vPAqxiNIHYBr2d0Y9tq4BPAyX0WKEkaxjj3QXwZ+BDw4qraOat9Osn7+ilLkjS0cQLiSVXVnDivqt6xwPVIkibEOCepr0ly1MxCkqOTfLbHmiRJE2CcgFhRVT+aWaiqHwLH9VeSJGkSjBMQ9yd53MxCksezn2c1SJKWj3HOQVwEfDHJF7rlXwPW9VeSJGkSzBkQVXV1kmcCpzF68tufVNX3eq9MkjSocaf7fgTwg67/KUmoqhv6K0uSNLRxbpR7B/By4Fbgga65AANCkpaxcUYQL2Z0L8RP+i5GkjQ5xrmK6Q7gYX0XIkmaLOOMIO4DtibZAvz/KKKq3tBbVZKkwY0TEJu6lyTt19T6zYNte8eGswfb9nI2zmWulyd5FPC4qrptEWqSJE2AcR45+lvAVuDqbnl1EkcUkrTMjXOS+u3AqcCPAKpqKz4DQpKWvXECYm9V3bNPm3MxSdIyN85J6luS/C5wRJJVwBuAL/VbliRpaOOMIF4PPIXRJa5XAD8G3tRnUZKk4Y1zFdN9jGZ0vaj/ciRJk2KcuZiuo3HOoaqe10tFkqSJMM45iD+d9fmRwG8De/spR5I0KcY5xHTTPk3/MevhQZKkZWqcQ0zHzFp8CPAs4LG9VaRlachpGCTNzziHmG5idA4ijA4tfQs4r8+iJEnDG+cQk3dNS9JhaJxDTC890PdV9amFK0eSNCnGOcR0HvBs4PPd8nOB64F7GB16MiAkaRkaJyAKOKWq7gZIshJ4T1W9ptfKJEmDGmeqjamZcOjsAp7YUz2SpAkxTkBcn+SzSX4vyVpgM3DdXD9KclmS3UlumdV2TJJrk9zevR/dtSfJJUm2J7k5yTPnvUeSpAUxZ0BU1QXA+4CnA6uBjVX1+jHW/QHgzH3a1gNbqmoVsKVbBnghsKp7rQPeO07xkqT+jHMOAuDrwL1V9bkkP5fkyKq690A/qKobkkzt03wOcHr3+XJGJ7sv7No/WFUF3JjkqCQr9zm0JUlaROM8cvQPgX8B/rFrOgH413lu7/iZP/rd+3Gz1nnnrH47uzZJ0kDGOQdxPvAcRs+BoKpu58E/7AsljbbmU+uSrEsynWR6z549C1yGJGnGOAHxk6r66cxCkocy/0eO7uouk525XHZ3174TOGlWvxOBu1orqKqNVbWmqtasWLFinmVIkuYyTkB8IclbgUcleT7wCeDf5rm9TcDa7vNa4MpZ7ed2VzOdBtzj+QdJGtY4AbEe2AN8A3gtcBXwZ3P9KMkVwJeBJyXZmeQ8YAPw/CS3A8/vlunWeQewHXg/8McHuR+SpAV2wKuYkhwBXF5Vr2L0h3tsVfXK/Xx1RqNvMTrXIUmaEAccQVTV/cCKJA9fpHokSRNinPsgdjB6itwm4H9mGqvqXX0VJUka3n5HEEk+1H18OfCZru+Rs16SpGXsQCOIZyV5PPAd4N2LVI8kaUIcKCDeB1wNnAxMz2oPo/sgntBjXZKkge33EFNVXVJVTwb+uaqeMOt1clUZDpK0zI0zm+sfLUYhkqTJMs6NcpKkw5ABIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElS04GeSS1JS8LU+s2DbHfHhrMH2e5icQQhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkpkFmc02yA7gXuB/YW1VrkhwDfAyYAnYAv1NVPxyiPknSsCOI51bV6qpa0y2vB7ZU1SpgS7csSRrIJD0P4hzg9O7z5cD1wIVDFbMcDTVnvqSlaagRRAHXJLkpybqu7fiquhugez+u9cMk65JMJ5nes2fPIpUrSYefoUYQz6mqu5IcB1yb5Jvj/rCqNgIbAdasWVN9FShJh7tBRhBVdVf3vhv4NHAqsCvJSoDuffcQtUmSRhY9IJI8OsmRM5+B3wBuATYBa7tua4ErF7s2SdKDhjjEdDzw6SQz2/9IVV2d5GvAx5OcB3wHeNkAtUmSOoseEFV1B/D0Rvv3gTMWux5JUpt3UkuSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqWvRnUkvScjG1fvNg296x4ezet+EIQpLUZEBIkpoMCElSkwEhSWryJPUAhjyxJUnjcgQhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDVN3GR9Sc4E/h44AvinqtrQx3acME+SDmyiRhBJjgDeA7wQOAV4ZZJThq1Kkg5PExUQwKnA9qq6o6p+CnwUOGfgmiTpsDRpAXECcOes5Z1dmyRpkU3aOYg02upnOiTrgHXd4n8nua33qkaOBb63SNtaDO7P5Ftu++T+LKC845B+/vhxOk1aQOwETpq1fCJw1+wOVbUR2LiYRQEkma6qNYu93b64P5Nvue2T+7P0TNohpq8Bq5KcnOThwCuATQPXJEmHpYkaQVTV3iQXAJ9ldJnrZVV168BlSdJhaaICAqCqrgKuGrqOhkU/rNUz92fyLbd9cn+WmFTV3L0kSYedSTsHIUmaEAbEQUjyF0luTrI1yTVJfmHomg5Fkncm+Wa3T59OctTQNR2KJC9LcmuSB5Is2atLkpyZ5LYk25OsH7qeQ5XksiS7k9wydC2HKslJSa5Lsq37v/bGoWvqkwFxcN5ZVU+rqtXAZ4A/H7qgQ3Qt8NSqehrwX8BbBq7nUN0CvBS4YehC5muZTjfzAeDMoYtYIHuBN1fVk4HTgPOXwb/PfhkQB6Gqfjxr8dHscxPfUlNV11TV3m7xRkb3nSxZVbWtqhbrxsm+LLvpZqrqBuAHQ9exEKrq7qr6evf5XmAby3i2h4m7imnSJbkYOBe4B3juwOUspN8HPjZ0EWpON/PLA9WiA0gyBTwD+MqwlfTHgNhHks8Bj218dVFVXVlVFwEXJXkLcAHwtkUt8CDNtT9dn4sYDZ0/vJi1zcc4+7PEzTndjIaX5DHAJ4E37XNkYVkxIPZRVb8+ZtePAJuZ8ICYa3+SrAV+EzijlsA1zwfx77NUzTndjIaV5GGMwuHDVfWpoevpk+cgDkKSVbMWXwR8c6haFkL3cKYLgRdV1X1D1yPA6WYmWpIAlwLbqupdQ9fTN2+UOwhJPgk8CXgA+Dbwuqr67rBVzV+S7cAjgO93TTdW1esGLOmQJHkJ8G5gBfAjYGtVvWDYqg5ekrOAv+PB6WYuHrikQ5LkCuB0RrOf7gLeVlWXDlrUPCX5VeDfgW8w+jsA8NZuBohlx4CQJDV5iEmS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZENICSvJL3fM1Hpnk0d0zA546dF3SfHijnLTAkvwl8EjgUcDOqvrrgUuS5sWAkBZYN4fS14D/BZ5dVfcPXJI0Lx5ikhbeMcBjgCMZjSSkJckRhLTAkmxi9CS4k4GVVXXBwCVJ8+LzIKQFlORcYG9VfaR7vvSXkjyvqj4/dG3SwXIEIUlq8hyEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU3/B9Sr5ziBb/ptAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 1000 samples from standard normal distribution\n",
    "sample = stats.norm.rvs(loc = 0, scale = 1, size = 1000)\n",
    "plt.hist(sample)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('frequency')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Compute p-values\n",
    "For each of the sampled values, we will compute the p-value, or the probability of observing a value _at least_ as extreme as the observed value under $H_0$. To do this, we can compute the cumulative density function (CDF) for the standard normal distribution at each sampled value. The CDF, evaluated at some value $x$, gives the probability that the specified distribution will take on a value less than or equal to $x$. This is the equivalent of a one-sided statistical test for the lower-tail.\n",
    "\n",
    "Note: For the purposes of this tutorial, we could also evaluate the survival function, or $1-CDF_{\\mu,\\sigma}(x)$, because this is equivalent to a one-sided statistical test for the upper tail.\n",
    "\n",
    "In this tutorial, will use the [scipy implementations](https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.stats.norm.html) of the CDF and the survival function - or $\\left(1-CDF \\right)$ - for the normal distribution: `scipy.stats.norm.cdf` and `scipy.stats.norm.sf`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAElCAYAAADqeCmyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGfNJREFUeJzt3XmYZXV95/H3R5pFFEXsxshmu+ACPKikVRIzkRESARUwowhBbQ2RmHHcE0XjjEzGBaMTjePaEQQVVEQj7ZaIyOISMI0iCsiIgHQLQrGK4oZ+88c5Lbfb2vreW11dv36/nqeeOves39+t7k+d+p17fidVhSSpXXeb7wIkSXPLoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBr01ekv2SrNlIxzonyV/200cl+cIY931Jkv366eOSfHiM+35NkvePa39qy6L5LkDaVFXVKcApM62X5CRgTVW9dob97TmOuvpfFh+uql0G9v3GcexbbfKMXppjSTyh0rwy6DW0JFcneXWSS5PckuQDSbaZYt1jk5y+3rx/SvKOfvp5SS5LcnuSK5P81TTHrSQPGXh9UpLXD7x+SpKLktya5GtJ9p5mX3+S5LtJbkvyTiADy56b5Cv9dJK8LckN/boXJ9kryTHAUcArk/wkyacH3ptXJbkY+GmSRf28AwYOv02Sj/Vt/kaSR87UxiT3AD4P7NQf7ydJdlq/KyjJIX1X0a19d9QjBpZdneRv+jbc1tcw6c9NbTDoNaqjgCcBDwYeCkzVffER4OAk9wJIsgVwOHBqv/wG4CnAvYDnAW9Lss+GFtNvcyLwV8B9gfcBK5NsPcm6i4FP9DUvBr4PPH6KXf8p8Md0bdweeCZwU1WtoOve+YequmdVPXVgmyOBJwPbV9Wdk+zzUODjwA5078Onkmw5Xfuq6qfAQcC1/fHuWVXXrteuh9K93y8FlgCfAz6dZKuB1Q4HDgQeCOwNPHe642phM+g1qndW1eqquhl4A124/Y6q+gHwDeCwftYTgTuq6vx++Wer6vvVORf4AvBfhqjn+cD7quqCqvp1VZ0M/ALYd5J1DwYurarTq+pXwNuBH02x318B2wEPB1JVl1XVdTPU8o7+vfnZFMsvHDj2PwLbTFHnhnom8NmqOrPf91uBuwN/uF5t1/Y/t08DjxrDcbWJMug1qtUD0z8AdgJI8vmBroWj+uWnctcvgj/nrrN5khyU5PwkNye5lS6EFw9RzwOAV/RdFrf2+9p1bV3r2Wmw/upG+Fs9yXpU1ZeAdwLvAq5PsmLtXyfTmHRfky2vqt8Aa6aoc0PtRPezGNz3amDngXUGf6HdAdxzDMfVJsqg16h2HZjeDbgWoKoOGuhaWPvJlY8D+yXZBXgafdD33SqfoDvzvF9VbU/X3RAmdwew7cDr3xuYXg28oaq2H/jatqo+Msl+rhusP0nWa886quodVfX7wJ50XTh/u3bRVJtMta/e4LHvBuxC//4xfRtn2u+1dL/w1u57bbt+OMN2apRBr1G9MMkuSXYAXgN8bKoVq2oCOAf4AHBVVV3WL9oK2BqYAO5MchBdn/hULgL+PMkWSQ4EnjCw7J+BFyR5XH8B9R5Jnpxku0n281lgzyR/1n8y5sWsG6i/leQx/T63BH4K/Bz4db/4euBB09Q7ld8fOPZL6bqYzp9FG68H7pvk3lPs9zTgyUn27+t9Rb/vrw1Roxpg0GtUp9L1p1/Zf71++tU5FTiAgW6bqrqdLmRPA26h69ZZOc0+XgI8FbiV7mLwpwb2tYqun/6d/b6uYIoLjVV1I/AM4HjgJmB34KtTHPNedL9EbqHrFrmJ7i8QgBOAPfquok9Nsf1kzqDrT78FeDbwZ32f+kxt/C7dxdYr+2Ou091TVZcDzwL+H3Bjv5+nVtUvN6A2NSQ+eETDSnI18JdV9cX5rkXS1Dyjl6TGGfSS1Di7biSpcZ7RS1LjDHrNqwwM3TvJsqGHJx5l27mUZGk/js2CHugsyXuT/M9+epN8r3WXBf2PTQvfuIbu3VBJCti9qq4YcT/HAQ+pqmeNpbB5sqGfoKqqF8xtRRonz+ilzchC/0tCwzHoN2NJHtEPYXtr34VyyMCyk5K8K8ln+2F0L0jy4IHlD09yZj82zeVJDp/mOIuTfKY/zs1Jvtzf8r92yNwD+um798e9JcmlwGPW289OST6RZCLJVUlePLBs2m3X2895/eS3+rF4ntnPn3J443RDDv+wfy8u7+86PZDubuBn9vv51qze+N9t08r+fbkiyfP7+dsk+Vm6ETZJ8tokd+au0T9fn+Tt/fTWSd6a5Jok1/fdKnfvl+2XZE1f/4/o7kpev4YP0Q1f8em+Ha/s5388yY/SDWV8XpI9B7ZZZ2hobeKqyq/N8AvYku6u0dfQDUHwROB24GH98pOAm4HH0nXxnQJ8tF92D7oxZZ7XL9uH7g7MPac41puA9/bH3JJuVMq1n/i6Gjignz4e+DLdsL27At+he3ITdCclFwL/q6/3QXR34j5ppm2nqKnoulzWvt6HbqjkxwFbAMv72rYGHta3d6d+3aXAg/vp4+ie9jTb931pf+xF/etzgXfTjVz5KLphIPbvl50H/Ld++gt0wygfNLDsaf302+nuJN6BboTNTwNv6pftB9wJvLlvy92nqOu3P4eBeX/R72/r/hgXDSw7CXj9wDGmfK/9mv8vz+g3X/vSjVh4fFX9srrRGT/DusMMf7Kqvl7dWOqncNdQtk8Brq6qD1TVnVX1DbpByZ4+xbF+BdwfeEBV/aqqvlx9QqzncLoByW6uqtXAOwaWPQZYUlV/39d7Jd2QBEfMYtvZmG5441/Thd0eSbasqqur6vsbuP/fkWRX4I+AV1XVz6vqIuD9dMMhQPdL4Al9d8vefZuekO4hIY8Bvpwkfe0v69t+O/BG7npfAH4DvK6qflFTD5n8O6rqxKq6vap+QfcL7ZHTjK+jTZhBv/naCVhd3RC2a/2A2Q1l+wDgcVl3KOCjgN9LslvuGp74J/36b6H76+EL6Z4edex0Na1Xz1oPoHuq0uAxXwPcbxbbzsaUwxtXd8H2pXRhd0OSj64/vsyQdgLWhvNg3Wt/BufSnS3vA3wbOJNucLN9gSuqG6tnCd0olxcO1P2v/fy1Jqrq52tfZPIhpNeRbjC145N8P8mP6c74YbihozXPvDCz+boW2DXJ3QbCfjfg/89i29XAuVX1J1MsX2ds8z7IXkEXpHsCZyf5j6o6a73t1g4bfMlAPYPHvKqqdp/imNNtOxtrhzd+w2QLq+pU4NS+j/x9dF0hz2bmIYOncy2wQ5LtBsJ+N+4aTvhrdN1GT6N7vy9NshvdU6vO7de5EfgZXbfZVMMQr1NjVR000zp0A8sdSjcA3dXAvekGX5tq6Ghtwjyj33xdQDfc7iuTbJnus+xPBT46i20/Azw0ybP7bbdMN4zvIyZbub/I+ZC+m+HHdF0hv55k1dOAVye5T7ox6180sOzrwI/7i4p3788490rymFlsO5n1hxaecnjjJA9L8sR04+b/nC5YB4coXpr+4nLf3uOSnDPD8em7mL4GvKm/+Lo3cDRdNxlVdQfddYkXclewf43uMYnn9uv8pq/9bUl27I+/c5InzXT8Gd6P7ei6rm6i+4vhjRu4P21CDPrNVHVD1h5C9/zRG+kuCD6nuiFwZ9r2drrx4o+gOyv9EXdd7JvM7sAXgZ8A/w68u6rOmWS9/03XdXEV3cXHDw0c89d0v4ge1S+/ka4/+94zbTuF44CT++6Ow2v64Y23prvYe2Pf1h3puo2ge5gKwE1JvtFP78rUwx2v70i6C7TXAv9C15d+5sDyc+kuYH994PV2dBdj13pVX+/5fTfLF+n+EtgQbwJe278ffwN8kO79/CFwKXeNk68FyLFupDFLchHdJ2dumu9aJDDoJal5dt1IUuMMeklqnEEvSY3bJD5Hv3jx4lq6dOl8lyFJC8qFF154Y1UtmWm9TSLoly5dyqpVq+a7DElaUJLM6g5wu24kqXEGvSQ1zqCXpMYZ9JLUuBmDPsmJSW5I8p2BeW9J8t0kFyf5lyTbDyx7df+knMuHGFhJkjRmszmjPwk4cL15ZwJ7VdXedMPavhogyR50A13t2W/z7iRbjK1aSdIGmzHoq+o8ukfKDc77Qv/UIehGtdulnz6U7nFzv6iqq+hG1HvsGOuVJG2gcfTR/wXw+X56Z9Z9ys8a1n1i0W8lOSbJqiSrJiYmxlCGJGkyIwV9kr+je/DwKWtnTbLapMNjVtWKqlpWVcuWLJnxxi5J0pCGvjM2yXK6h0TvP/Cg5zV0D11Yaxe6BypI0iZr6bGfnbdjX338k+f8GEOd0Sc5kO6pNof0jztbayVwRJKtkzyQ7slCX59sH5KkjWPGM/okH6F7Ev3iJGuA19F9ymZr4MzuMaCcX1UvqKpLkpxG9+ixO4EX9o+AkyTNkxmDvqqOnGT2CdOs/wbgDaMUJUkaH++MlaTGGfSS1DiDXpIaZ9BLUuM2iSdMSdp0tP6Z8s2RZ/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxnnDlKRNxnzerNUyz+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4/x4pTSD+frIn2Oza1w8o5ekxi34M3qfhiNJ0/OMXpIaZ9BLUuMMeklq3IxBn+TEJDck+c7AvB2SnJnke/33+/Tzk+QdSa5IcnGSfeayeEnSzGZzRn8ScOB6844Fzqqq3YGz+tcABwG791/HAO8ZT5mSpGHNGPRVdR5w83qzDwVO7qdPBg4bmP/B6pwPbJ/k/uMqVpK04Ybto79fVV0H0H/fsZ+/M7B6YL01/bzfkeSYJKuSrJqYmBiyDEnSTMZ9MTaTzKvJVqyqFVW1rKqWLVmyZMxlSJLWGjbor1/bJdN/v6GfvwbYdWC9XYBrhy9PkjSqYYN+JbC8n14OnDEw/zn9p2/2BW5b28UjSZofMw6BkOQjwH7A4iRrgNcBxwOnJTkauAZ4Rr/654CDgSuAO4DnzUHNkqQNMGPQV9WRUyzaf5J1C3jhqEVJksbHO2MlqXELfvTK+eQ45ZIWAs/oJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcd4ZuwDN1x2588m7gaXheUYvSY3zjF4Lwub4V4w0Lp7RS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhrnDVPSJsqbxDQuntFLUuMMeklq3EhBn+RlSS5J8p0kH0myTZIHJrkgyfeSfCzJVuMqVpK04YYO+iQ7Ay8GllXVXsAWwBHAm4G3VdXuwC3A0eMoVJI0nFG7bhYBd0+yCNgWuA54InB6v/xk4LARjyFJGsHQQV9VPwTeClxDF/C3ARcCt1bVnf1qa4CdRy1SkjS8Ubpu7gMcCjwQ2Am4B3DQJKvWFNsfk2RVklUTExPDliFJmsEoXTcHAFdV1URV/Qr4JPCHwPZ9Vw7ALsC1k21cVSuqallVLVuyZMkIZUiSpjNK0F8D7Jtk2yQB9gcuBc4Gnt6vsxw4Y7QSJUmjGKWP/gK6i67fAL7d72sF8Crg5UmuAO4LnDCGOiVJQxppCISqeh3wuvVmXwk8dpT9SpLGxztjJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNW6koE+yfZLTk3w3yWVJ/iDJDknOTPK9/vt9xlWsJGnDjXpG/0/Av1bVw4FHApcBxwJnVdXuwFn9a0nSPBk66JPcC/hj4ASAqvplVd0KHAqc3K92MnDYqEVKkoY3yhn9g4AJ4ANJvpnk/UnuAdyvqq4D6L/vONnGSY5JsirJqomJiRHKkCRNZ5SgXwTsA7ynqh4N/JQN6KapqhVVtayqli1ZsmSEMiRJ0xkl6NcAa6rqgv716XTBf32S+wP0328YrURJ0iiGDvqq+hGwOsnD+ln7A5cCK4Hl/bzlwBkjVShJGsmiEbd/EXBKkq2AK4Hn0f3yOC3J0cA1wDNGPIYkaQQjBX1VXQQsm2TR/qPsV5I0Pt4ZK0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEjB32SLZJ8M8ln+tcPTHJBku8l+ViSrUYvU5I0rHGc0b8EuGzg9ZuBt1XV7sAtwNFjOIYkaUgjBX2SXYAnA+/vXwd4InB6v8rJwGGjHEOSNJpRz+jfDrwS+E3/+r7ArVV1Z/96DbDzZBsmOSbJqiSrJiYmRixDkjSVoYM+yVOAG6rqwsHZk6xak21fVSuqallVLVuyZMmwZUiSZrBohG0fDxyS5GBgG+BedGf42ydZ1J/V7wJcO3qZkqRhDX1GX1WvrqpdqmopcATwpao6CjgbeHq/2nLgjJGrlCQNbS4+R/8q4OVJrqDrsz9hDo4hSZqlUbpufquqzgHO6aevBB47jv1KkkbnnbGS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaN3TQJ9k1ydlJLktySZKX9PN3SHJmku/13+8zvnIlSRtqlDP6O4FXVNUjgH2BFybZAzgWOKuqdgfO6l9LkubJ0EFfVddV1Tf66duBy4CdgUOBk/vVTgYOG7VISdLwxtJHn2Qp8GjgAuB+VXUddL8MgB2n2OaYJKuSrJqYmBhHGZKkSYwc9EnuCXwCeGlV/Xi221XViqpaVlXLlixZMmoZkqQpjBT0SbakC/lTquqT/ezrk9y/X35/4IbRSpQkjWKUT90EOAG4rKr+cWDRSmB5P70cOGP48iRJo1o0wraPB54NfDvJRf281wDHA6clORq4BnjGaCVKkkYxdNBX1VeATLF4/2H3K0kaL++MlaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1Lg5C/okBya5PMkVSY6dq+NIkqY3J0GfZAvgXcBBwB7AkUn2mItjSZKmN1dn9I8FrqiqK6vql8BHgUPn6FiSpGksmqP97gysHni9Bnjc4ApJjgGO6V/+JMnlQx5rMXDjkNsuVJtjm2HzbLdtblzeDAzf5gfMZqW5CvpMMq/WeVG1Algx8oGSVVW1bNT9LCSbY5th82y3bd48zHWb56rrZg2w68DrXYBr5+hYkqRpzFXQ/wewe5IHJtkKOAJYOUfHkiRNY066bqrqziT/A/g3YAvgxKq6ZC6OxRi6fxagzbHNsHm22zZvHua0zamqmdeSJC1Y3hkrSY0z6CWpcQsm6GcaUiHJ1kk+1i+/IMnSjV/leM2izS9PcmmSi5OclWRWn6ndlM126IwkT09SSRb8x/Bm0+Ykh/c/60uSnLqxa5wLs/j3vVuSs5N8s/83fvB81DkuSU5MckOS70yxPEne0b8fFyfZZ2wHr6pN/ovugu73gQcBWwHfAvZYb53/Dry3nz4C+Nh8170R2vxfgW376b/eHNrcr7cdcB5wPrBsvuveCD/n3YFvAvfpX+8433VvpHavAP66n94DuHq+6x6xzX8M7AN8Z4rlBwOfp7sPaV/ggnEde6Gc0c9mSIVDgZP76dOB/ZNMduPWQjFjm6vq7Kq6o395Pt39CgvZbIfO+D/APwA/35jFzZHZtPn5wLuq6haAqrphI9c4F2bT7gLu1U/fmwV+L05VnQfcPM0qhwIfrM75wPZJ7j+OYy+UoJ9sSIWdp1qnqu4EbgPuu1GqmxuzafOgo+nOBhayGduc5NHArlX1mY1Z2Byazc/5ocBDk3w1yflJDtxo1c2d2bT7OOBZSdYAnwNetHFKmzcb+n9+1uZqCIRxm3FIhVmus5DMuj1JngUsA54wpxXNvWnbnORuwNuA526sgjaC2fycF9F13+xH91fbl5PsVVW3znFtc2k27T4SOKmq/m+SPwA+1Lf7N3Nf3ryYswxbKGf0sxlS4bfrJFlE96fedH8mbepmNYxEkgOAvwMOqapfbKTa5spMbd4O2As4J8nVdP2YKxf4BdnZ/ts+o6p+VVVXAZfTBf9CNpt2Hw2cBlBV/w5sQzf4V6vmbOiYhRL0sxlSYSWwvJ9+OvCl6q9wLFAztrnvxngfXci30G87bZur6raqWlxVS6tqKd11iUOqatX8lDsWs/m3/Sm6C+8kWUzXlXPlRq1y/GbT7muA/QGSPIIu6Cc2apUb10rgOf2nb/YFbquq68ax4wXRdVNTDKmQ5O+BVVW1EjiB7k+7K+jO5I+Yv4pHN8s2vwW4J/Dx/rrzNVV1yLwVPaJZtrkps2zzvwF/muRS4NfA31bVTfNX9ehm2e5XAP+c5GV0XRjPXcgnb0k+Qtf9tri/7vA6YEuAqnov3XWIg4ErgDuA543t2Av4fZMkzcJC6bqRJA3JoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+ilSSRZmuS7SU7uxwY/Pcm2812XNAyDXpraw4AVVbU38GO6Zx5IC45BL01tdVV9tZ/+MPBH81mMNCyDXpra+uODOF6IFiSDXprabv046NCNjf6V+SxGGpZBL03tMmB5kouBHYD3zHM90lAWxDDF0jz5TVW9YL6LkEblGb0kNc7x6CWpcZ7RS1LjDHpJapxBL0mNM+glqXEGvSQ17j8BMbrHOwVryZUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute p-values\n",
    "pvals = stats.norm.cdf(sample, loc = 0, scale = 1)\n",
    "plt.hist(pvals)\n",
    "plt.xlabel('p')\n",
    "plt.title('p-value distribution\\none-sided test, lower-tail')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### QQ-plot of p-values\n",
    "\n",
    "Here, you will plot your own QQ-plot and visualize what a QQ-plot should look like when the observed values follow the theoretical distribution.\n",
    "\n",
    "The QQ-plot is a plot of one set of quantiles against another. When testing for concordance/discordance of one distribution against another (e.g. observed vs. theoretical), you will plot the quantiles of the theoretical distribution on the x-axis and the quantiles of the observed distribution on the y-axis.\n",
    "\n",
    "What's our theoretical distribution? Under $H_0$, the p-value is **uniformly distributed**. Therefore, for $N$ observed values, the quantiles of the theoretical distribution can be computed by:\n",
    "\n",
    "$$\\frac{i}{N} \\forall i\\in[0,...,N]$$\n",
    "\n",
    "To obtain the quantiles of our observed distribution, we will simply sort the values, as the p-values are essentially percentiles already!\n",
    "\n",
    "But wait - we have to make one more adjustment! Because people are usually interested in small p-values, it's more comprehensible to plot the p-values on a log-scale. Standard practice is to use $-\\log_{10}p$. Thus, we are actually plotting the sorted values of $-\\log_{10}p_{\\text{observed}}$ on the y-axis, and the sorted values of $-\\log_{10}(p_{\\text{expected}})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEdCAYAAAD5KpvoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYU9XWx/HvYhiKCAKCSkdREUVEHSt6BRt2FBvKVbGhKCr3KooFaXqx3NdesPeCBRErgoKFiyjNAoigojQVRJoMZWC9f5wzGEIykwmTZDLz+zzPPCQ5+5ysQyBr9t7n7GXujoiISCyVMh2AiIiUXUoSIiISl5KEiIjEpSQhIiJxKUmIiEhcShIiIhKXkoSIiMSlJCEiInEpSUhGmdkcMzsyTe/V0symmNkKM7syHe8ZJ460nXMqmNnTZnZLpuOQ9FCSkI3MrJuZfWNmq8zsVzN7yMy2KWmbFMa3pV+u1wJj3b2mu99XWnHJljGzOmbmZrZj1Ov3m9ljmYpLAkoSAoCZXQ3cDvQGtgEOBJoDH5hZbqJtyrhmwLRMB1FRmVnlOJvaAiuBOVGvtwampjImKZ6ShGBmtYABwBXu/r67r3P3OcAZwI7A2Ym0KeL4c8zsejObbmZ/mtlTZlYtTttWZjbWzJaa2TQzOyl8/TmgKfCWma00s2tLuP9HQAfggXD/Xbcwzj5m9lrUa/ea2X0R238Ih7amm9kpRfz9uJntHPF843COmTU0s9fNbJGZ/RQ9TGZm15nZ/PB9ZprZEUm8xxwzu8bMvjazZWY2NPK8zWxvM5scvsdQoFrUsePGGB77OjP7GvgrTqJoC0z3zReS2wOYEu/vTdJDSUIADib4jz8s8kV3Xwm8BxydYJuidAU6Ai2AXYGbohuEvZG3gA+A7YArgBfMrKW7nwP8Apzo7lu7+x0l3P9w4FOgZ7j/98nGGXoJOC5MnphZDkHCfDHc/gNwKEGPawDwvJk1iPeXE4uZVQrP5yugEXAE0MvMOobbWwI9gf3cvWYY95ySvEeEM4BjCBJ+G6Bb+B5VgOHAc0Bd4FXg1ERjDJ0FHA/UdveCGO+9N/Bt1LlvD9QDvk7yfKSUKEkIBP8ZF8f5D7wQqJ9gm6I84O5z3X0JcCvBF0e0A4Gtgdvcfa27fwS8HadtLFu6f6Jx4u4/A5OBk8OXDgdWufvn4fZX3X2Bu29w96HALGD/EsQBsB9Q390HhufzI/AY0CXcvh6oCuxuZrnuPsfdfyjhexS6L4x3CcGXftvw9QOBXOCesPf4GvBlCWIsPPZcd8+P895tgbPMbHHhD/AdMCv8JQQzu9XMPjWz18xsqyTPUZKgJCEAi4F6cYYCGgCLEmljZl3DoZyVZvZeVJu5EY9/BhrGOE5DYK67b4hq2yjB89jS/ePGGefcXuTvJHI2f/ciMLNzzWxqOOy1lGB8vV4J4oBgDqVh4THC49wAbA/g7rOBXkB/4Hcze9nMYv29JuLXiMerCJItBOc/P2oo6OdEYwxF/p1uwsyqAq0IenBtI34eJpyPMLPWQAt3PxQYDVxQ8tOTZClJCMB4YA3QOfJFM6sBHAt8nEgbd38hHMrZ2t2PjXqPJhGPmwILYsSxAGgSDmFEtp0fPi6u+Elx+yciZpxxzu1VoL2ZNQZOIUwSZtaM4LfpnsC27l6bYDjF4rznKiDyt+Mdwj/nAj+5e+2In5ruflxhQ3d/0d0PIfiydoILC0ryHsVZCDQys8jYm0Y8LjZGiv7cWhN8D41y93mFP8Au/D0fcSjBkCbhn4ckGLuUAiUJwd2XEYyb329mx5hZrpk1J/gSXAy8kEibYt7mcjNrbGZ1CX7THBqjzQTgL+Da8PjtgROBl8PtvwE7FfEexe2fiETiBMDdFwFjgacIvihnhJtqEHwxLgIws/MJvgzjmUpwcUCOmR0DHBa+/gWwPJz4rR5ub21m+4XHbWlmh4e/ja8G8gmGoEryHsUZDxQAV5pZZTPrzKbDZkXGmIC9CSatV0a9vh9/X9lUB1gWPl5GMDciaaIkIQCEE8E3AP8FVgA/EfzmeaS7/5VomyK8SDCh/GP4s9nNWO6+FjiJoGeyGHgIONfdvwubDAZuCoc1rkli/0QUG2eM9kcSMdTk7tOB/yP4gv0N2BMYV8QxriJIZksJhl2Gh8dZH77eluDvejHwOMFkOATzEbeFr/9KMFl/Q0neozjh32lngonsP4Ezibh4IYEYi9OWTec4MLP6BD2jwiTxZ8TxtgGWJHhsKQWm8qUSi5ldQNBzaOfuvyTbJmw3B7jI3UenItbSki1xVjRmtidwvbufbWbdgarufn+m46oo4t3cIhWcuz9pZusILn2NmQASaSOypdz9GzP72cw+BX4Hzs10TBWJkoTE5e7PlUYbkS3l7tdnOoaKSsNNIiISlyauRUQkLiUJERGJK+vnJOrVq+fNmzfPdBgiIlll0qRJi929uOV0sj9JNG/enIkTJ2Y6DBGRrGJmPxffSsNNIiJSBCUJERGJS0lCRETiUpIQEZG4lCRERCQuJQkREYkr6y+BFRGpSIZPmc+dI2eyYGk+DWtXp3fHlpy8d0mKL5aMkoSISJYYPmU+1w/7hvx1QW2p+UvzuX7YNwApSxQabhIRyRJ3jpy5MUEUyl+3njtHzkzZeypJiIhkiQVL80v0emlQkhARyRINa1cv0eulQUlCRCRL9O7Ykuq5OZu8Vj03h94dW6bsPTVxLSKSJQonp3V1k4iIxHTy3o1SmhSiabhJRETiUpIQEZG4lCRERCSutCUJM6tmZl+Y2VdmNs3MBsRoU9XMhprZbDObYGbN0xWfiIhsLp09iTXA4e6+F9AWOMbMDoxqcyHwp7vvDNwN3J7G+EREJErakoQHVoZPc8Mfj2rWCXgmfPwacISZWZpCFBGRKGmdkzCzHDObCvwOjHL3CVFNGgFzAdy9AFgGbJvOGEVEyrS//oI+feDZZ9PydmlNEu6+3t3bAo2B/c2sdVSTWL2G6N4GZtbdzCaa2cRFixalIlQRkbLFHYYNg1at4Pbb4euv0/K2Gbm6yd2XAmOBY6I2zQOaAJhZZWAbYEmM/R919zx3z6tfv36KoxURybAffoDjj4dTT4XateGzz+C//03LW6fz6qb6ZlY7fFwdOBL4LqrZCOC88PFpwEfuvllPQkSkQli9GgYMgD32gE8/hbvugsmToV27tIWQzmU5GgDPmFkOQXJ6xd3fNrOBwER3HwE8ATxnZrMJehBd0hifiEjZ8d57cMUVQS+iSxf4v/+Dhg3THkbakoS7fw3sHeP1myMerwZOT1dMIiJlzty50KtXMP/QsiWMHg1HHJGxcHTHtYhIWbB2bTAhvdtuQS/iP/+Br77KaIIArQIrIpJ5Y8fCZZfBjBnQqRPccw80b57pqAD1JEREMufXX+Gf/4QOHSA/H956C4YPLzMJApQkRETSr6AA7r8/mHN49VXo2xemT4cTTsh0ZJvRcJOISDp9/jn06AFTp8LRRwfJYtddMx1VXOpJiIikwx9/wMUXw0EHwaJF8Mor8P77ZTpBgJKEiEhqbdgAjz8eJIOnnoKrrw4mqE8/HbJg/VINN4mIpMqUKcFVS59/DoceCg8+CHvumemoSkQ9CRGR0rZsGVx5JeTlwY8/wjPPwMcfZ12CAPUkRERKjzu8+GIwpPT770Ev4pZbgkX5spSShIhIaZg+HS6/PLgxbr/94J13YN99Mx3VFtNwk4jIlli5Eq67DvbaK1hGY8gQGD++XCQIUE9CRCQ57vDGG8FifHPnwvnnB2svlbMaN+pJiIiU1OzZcNxxQRGgOnWCIkBPPlnuEgQoSYiIJG71aujfH1q3hnHj4O67YdKktBYBSjcNN4mIJKKMFAFKN/UkRESK8ssv0LlzMLxUuXJQBOillypEggAlCRGR2AqLALVqFayxVEaKAKWbhptERKIlUARo+JT53DlyJguW5tOwdnV6d2zJyXs3yki4qaQkISIVSpFf7gsXwjXXBHdN77gjvP02HH98zGNcP+wb8tetB2D+0nyuH/YNQLlLFBpuEpEKo/DLff7SfJy/v9zf/PJnuO++oL70a68FRYCmTYuZIADuHDlzY4IolL9uPXeOnJmGs0gv9SREpMKI9eW+25xptOp0BSycHRQBeuAB2GWXIo+zYGl+iV7PZmnrSZhZEzMbY2YzzGyamV0Vo017M1tmZlPDn5vTFZ+IlH+RX+J1Vi3jtvfu443nr6Hmij+DMqLvv19sggBoWLt6iV7PZunsSRQAV7v7ZDOrCUwys1HuPj2q3afuXvYKvYpI1mtYuzoL/vyLM7/6gOs+foaaa/7ikf0789px5zPqtMS/dnp3bLnJnARA9dwcendsmYqwMyptScLdFwILw8crzGwG0AiIThIiIilxS7N11H2gN3vNn8mExnvQ9+gezG3YgsEnlazOQ+HktK5uShEzaw7sDUyIsfkgM/sKWABc4+7T0hiaiJRHy5ZB3750ePBBVteuy8DT+/DUju1oWGcrBif55X7y3o3KZVKIlvYkYWZbA68Dvdx9edTmyUAzd19pZscBw4HNBgjNrDvQHaBp06YpjlhEspY7vPBCcFnrokXQowfVbrmFm2vXRhOeiUnrJbBmlkuQIF5w92HR2919ubuvDB+/C+SaWb0Y7R519zx3z6tfDlddFJFSMH06dOgA55wDzZrBF18EVy5lcZW4TEjn1U0GPAHMcPe74rTZIWyHme0fxvdHumIUkXIgsgjQ11/DI4+UqyJA6ZbO4aZ2wDnAN2Y2NXztBqApgLsPAU4DephZAZAPdHF3T2OMIpJFNrl7eptq3J37A/vfOwjmzYMLLoDbbiuXNR7SKZ1XN30GWDFtHgAeSE9EIpLNIpfGaPbnAga88gj7/zSJZbu0Yptx4+DggzMdYrmgO65FJCvdOXImG1atoteE1+jx+WuszanMgCMu5sMOp/OJEkSpUZIQkay066RPeHH0IzRb+itvtjqMWztcwO81t8VWrM10aOWKkoSIZJdffoFevXjqjTf4oW5jzupyK+Ob7bVxc3lcGiOTlCREJDusXRvUlB44ENyZ1rMPZ9U8mOUb/r5Is7wujZFJWipcRMq+MWOgbVvo0ydYqXXGDPa4fzADT9+HRrWrY0Cj2tUZ3HnPCnEXdDqpJyEiZVcxRYAqytIYmaSehIiUPQUFmxYBuvnmIosASeqoJyEiZcv48UF96alToWNHuP/+hGo8SGqoJyEiZcPixXDRRXDwweQv+JUbz76ZHffqSbvX5zJ8yvxMR1dhqSchIpm1YQM88UQwKb18ObPOvZQuOxzNH1YF+LsONaD5hwxQT0JEMmfy5GD5jO7doXVrmDqVbq1O35ggCuWvW8+dI2dmKMiKTUlCRNJv6VK44grYbz/46Sd49lkYOxb22GOTOtSR4r0uqaUkISLp4w7PPx9ctfTQQ8EE9cyZQc2HoEpA3DumdSd1ZihJiEh6TJu2aRGgL78MrlyKKgLUu2NLqufmbPKa7qTOHCUJEUmtlSvh2muDO6YjiwDts0/M5ifv3YjBnffUndRlRImvbjKzrsCX7v59CuIRkfLCHYYNg169giJAF14YFAGqt1lF4s3oTuqyI5lLYBcBD5lZFWAx8L279yndsEQkq82eDT17wsiRQRnRoUNVBChLlXi4yd0/ACa4+z+A84CtSz0qEclO+fnQr19wOev//gf33AMTJypBZLFkb6arZWb7AN8CNUoxHhHJVu++G1zW+uOPcPbZ8N//QoMGmY5KtlCySeLfQA+gJzCy9MIRkazzyy9w1VUwfHhwaeuHH8LhhwNBHeo7R85kwdJ8GtauTu+OLTl570ZxX5eyJ6kk4e7rgPtKORYRySZr18Jdd8GgQRRscB7reBF373k89b+A3nWCtZauH/YN+evWA38vrzHx5yW8Pmn+Zq+Dlt0oi8zdS76T2YNADXfvZmZHh/MUGZGXl+cTJ07M1NuLVExjxgQ3wn33HQs6dOScPc7ihxp/X7VUPTeHarmV+HPVus12zTFjfYzvnUa1qzOuz+EpDVv+ZmaT3D2vuHbJ3iexFvgxfJzQp2pmTcxsjJnNMLNpZnZVjDZmZveZ2Wwz+zqc9xCRsmLhwmC+4fDDg57EO+9w+tHXbpIgIFhrKVaCAGImCNCyG2VVsnMSq4BtzCwXaJrgPgXA1e4+2cxqApPMbJS7T49ocyywS/hzAPBw+KeIZFJBATz4IPTty/rVa3jm8HO4o+3JbPt1NeaX8Ms9Xk9Cy26UTcn2JAYDPwAPAi8msoO7L3T3yeHjFcAMIHoAshPwrAc+B2qbmS6PEMmk8eMhLw969eK31vtw/MUPMXC/M1mdW5X5S/OxOLvVrp4bc3mNsw5oomU3skiyPYl7gQnApe6+oaQ7m1lzYO/wGJEaAXMjns8LX1sYtX93oDtA06aJdmREpEQWL+bni3rS7M2hLKhZjwe79uPdXQ7iz/yCTZo5YOGfharn5tD/pD0AYl7FlNesrq5uyhJJTVwDmNmhQDfgTXcfUYL9tgY+Bm5192FR294BBrv7Z+HzD4Fr3X1SvONp4lqklG3YAI8/ztprr8NWrOCJvE7c1+4sVlUpejioUe3q+tLPIolOXCfVkzCzfwCtgXzgFiChJBHOYbwOvBCdIELzgCYRzxsDC5KJUURKbsyL77Ndn3+zx9wZTGnSmptO7cGs+s2K3U9XJpVfyQ43HQJ8QvBlvyyRHczMgCeAGe5+V5xmI4CeZvYywYT1MndfGKetiJSWpUv58ZJe/OPV51iyVS16nXA1w3dvv7HGQ1E0n1C+JZskCgqHhEqgHXAO8I2ZTQ1fu4Hw6ih3HwK8CxwHzCa4gur8JOMTkQQMnzyPr25/kMvfHkKz/OU8t/dx3HXoP1leLf6SbLWr51KjamUNLVUQySaJTma2Chjl7gkVng2TSpG/lngwQXJ5kjGJSAIKl8TYatZ3DBr1MCfP/ZapDXblvNP7M22HnYvct3BCWkmh4kg2SXQmuDrpFDPb2d0vKsWYRKQURa6TtE31XHzlCi775EUumPgmK6tsRZ+OPRm619G4xb4iPseMDe7qNVRQySaJm/l7WY6OpRmQiGy5wsRQeB+DA7hz0JQx9P3ocRquWMzLbY7m9sPO48+ttol7nOq5OaoKV8ElmyTWAr+FjzuglWBFMi5mYgj/bL5kPgNGP8JhP01m+nY70rPTdUxu1CrmcdRzkEjpXJZDRFJg+JT59B8xjaX5f6+VVJggqq5bw2Wfv8alE15lbU4u/Y/oznP7HM/6Sjkxj6Weg0RLNkmMAXYFhgAvlF44IpKIeL2GSB1++JIBo4bQdNlvDN/9MG7tcCGLtq67WbvC/Rup5yAxJJskjiaoJ9Ex/Hmn1CISkbiK6jUUarTsd27+8FE6zvqc2XUbc1aX/zC+WZuN23MrGVtXq8zSVes0pCTFSjZJ1AauA64FLiy9cEQknuFT5m9SxCda7vp1XPTlcK4c9zJucNth3Xhiv06sy8lVb0GSlmySGAjs5u4zzazEC/yJSMndOXJm3ARx0M9fMeiDh9l5yTze3/UgBh5xMQtrbafEIFss2fKl88xsefi4T+mGJCKFIuceYqm/cgk3ffQEnWZ8zM+1d6Dbaf0Y22I/6myVy90n6qY32XLJLvD3ErDMzKoDT7j7J6UblkjFlcikdM6G9Zw7+W3+/enzVFlfwD3tzuLhA06jXv3a3KNeg5SiZIebZrj7QAAze4hgsT8R2ULR8w6xEsQ+82Zwy6iH2P33nxi7477cdeLlXNDtaGYqMUgKlDhJmNljwG5mthT4Glha6lGJVDDFDSsB1F21jOvGPs2Z34xiQc16XHLyDRzbtwcj9mmcxkiloilxknD3i82sGsHaTfsD25Z6VCIVRKxLWqOZb6DLVx9w3cdPU2NtPkMOOJX7Du5Cne3q8ogShKRYssNNNwCtgL8I6l2LSDESmWuI1vrX2dzywUO0Xfg9nzdpTd+jgiJAquEg6ZJskqjr7qebWRXgbrS8t0iREplriFRr9Uqu/vQ5zpn8Ln/U2IarTriat/fowHp0SaukV7JJYo2Z7QN8A9QoxXhEypVEhpM24c4p08Zww5gnqZu/nGf2PYEnj+zGp7d15t7UhioSU7JJ4kbgivBnaOmFI1I+lDg5ALsumsOgUUM4YO63TGnQkm6n9+fHJi0Z3HnPFEYqUrRkb6ZbDdxZyrGIZK1Erk6KZ6u1+Vw17qWNRYCuP6YnL7c5moZ1ajBYw0qSYQknCTM7FuhPsG7TV8Bd7v55iuISyRrFrakUlzvHzhzHzR8+RoOVfzBsn2OoducdDD58T10NImVGSXoSDwH/BKYD+wL/NbMH3f2llEQmUsYlM6RUqPmS+QwcNYR/zJnCdzu04Pv7n6Rzt5NSEKXIlilJkvjN3ceFj0eb2XhgAqAkIRXOTcO/4fnPfynxflXXreGKCa/S44th5FSvBvfey26XXcZulZOdHhRJrZL8y5xjZrcAA919LbAOWJHozmb2JHAC8Lu7t46xvT3wJvBT+NKwwqU/RMqKLek9dPjhSwZ9+AiN//wVunaFO++EBg1SEKVI6SlJknCgM3Cxmc0iKFv6gpnt4u6zEtj/aeAB4Nki2nzq7ieUICaRtEl27qHRst95bPJz7P7FGGjVCl5/ETp0SFGUIqUr4STh7mcBhEtytAb2AtoAj5vZTu7epJj9PzGz5smHKpJZRdVziCV3/TqunDKCHuNeonKlSnDbbfCvf0GVKimMUqR0JTMQegrQCVgPVAKGlOLk9UFm9hWwALjG3aeV0nFFtliil7fm5hhPNVnBIff0g5kzoXNnuPtuaNo0xRGKlL5kkkR7d+9S+MTMHqR0Jq8nA83cfaWZHQcMB3aJ1dDMugPdAZrqP56kwfAp8xNab2nXguU8OXUojf/zJuy0E7z7Lhx7bDpCFEmJSknsU9XMjjezNmZ2PFC9NAJx9+XuvjJ8/C6Qa2b14rR91N3z3D2vfv36pfH2IkUa8Na0uAmizla53Htqa+bsMIsPHr2ExmPfh3794NtvlSAk6yXTk7iMYAJ7T2Au0LM0AjGzHQgus3Uz258ggf1RGscWSVYwWf01+evil3KfckQNuOgk+PprOOYYuP9+2HnnNEYpkjrJ1JNYBTxf0v3CkqftgXpmNg/oB+SGxxwCnAb0MLMCIB/o4u6JrKYskhLF3QtRd9UyBv7vObj9fWjSBF5/HU45BczSGKVIam3xHTxmdp27315cu8Kro4rY/gDBJbIiGVdUgjDfwFlfjeTaj5+hVsFquO466NsXamhBZCl/kilf+krkU6AtUGySEMkWw6fM54U4CSKyCNDEZm3Ie/cl2H33NEcokj7J9CSWu/tFhU/M7OFSjEck42JNUkcXAep1wtW0H9CLvN1VPlTKt2SSxK1Rz28sjUBEMqXIyWl3Ok/7iOvHPLWxCNBdh/6TToe14mTVl5YKoCRLhf874nHhw2XAJGBJ6YYlkh5FzT1sVgTojAHMabIrt56yp2o8SIVRkp5EXvjzVvj8eOBL4FIze9Xd7yjt4ERSKV6CqLFm1cYiQCuq1uC6Y67glTZH0fWg5rxzsqrEScVSkiSxLbBP4Q1vZtYPeA34B0FvQklCskLc4SV3jps5jr5hEaCX2hzNHYedx59bbUPt6rncogQhFVBJkkRTYG3E83UEy2jkm9ma0g1LJDXi9R4iiwBN224nLjv5eqY02g0ILuHrf9IeaY5UpGwoSZJ4EfjczN4k+H9zAvCSmdUgqFYnUmbF6z1UXbeGy8e/wiVfvM6anCr0O/ISnt/7ONZXytnYpuuBTTUHIRVWSZYKH2Rm7wKHECSJS919Yri5ayqCEykNw6fM59+vTGVD1HWth8/+gv6jH6Hpst94Y/f2/KfDhSzaus7G7XW2yqXfiXsoQUiFVtJLYAuADQSLYZa8NJdImg2fMp9/DZ26yX0PjZb9Tr8PH+XoWZ8za9smdDnrP3zetA0A/zywqeYeRCKU5BLYq4CLgdcJehLPm9mj7n5/qoIT2RJdHxvPuB/+vjq7SsE6LvryDa7431DcYHD7bjyZ14l1ObmAEoRILCXpSVwIHODufwGY2e3AeEBJQsqU4VPm0/vVqUROPxw8ZyqDRg2hxZJ5vLfrwQw64iIW1Npu43YlCJHYSpIkjKAaXaH14WsiGRfvqqXtVvzBTWOe4KQZnzCndgO6ndafsS3yNm6vUSVHN8eJFKEkSeIpYIKZvRE+Pxl4ovRDEklMUctp5GxYz3mT3uZfnz1PlfUF3N3ubIYceBprKv9dX1q9B5HileTqprvM7GOgHUEP4nx3n5KyyERiSKQI0L7zpnPLBw/RatEcxuy0L/2OvJRf6jTYpI0ShEhiSnR1k7tPIri7WiStYs0zRKu7ahl9xj7FGd+MZn7N+lxyyg2M3OWgTYoAVa1cidtPbaPhJZEEFZskzGwFseu/G+DuXqvUoxKJUFyFuEob1tPl6w+49uNnqLE2n4cPOI37Du5CfpVqm7Rr16IuL1x8UKrDFSlXik0S7l4zHYGIxDJ8yvwiE0RQBOhB2i6cxfime9L3qB7Mrtd0s3YaXhJJTjKV6U5097eKbymyZYZPmU+voVNjbqu1eiXXfPIc/5wSFAG68sRrGNHqsM3qS+uuaZEtk2zRISUJSam4Q0xhEaAbxjxJnfwVG4sAraga1JfWnINI6UomSejeCEmp6DulC+26aA6DPniYA+ZNY3LDlpx3xkCmbd8CUI9BJFWSSRKxJrFFSsVNw7/ZLEFEFgFaXm1rrj3mSl5tcyRVcitzj3oNIimVTJJIipk9SbC8+O/u3jrGdgPuBY4DVgHd3H1yuuKTzIo5vBRVBOjFvTpyx2HnsbR6LU1Ei6RJ2pIE8DTwAPBsnO3HAruEPwcAD4d/SjkWb+5hxyXzGRAWAfp2+xabFAFq16KuEoRImiSTJH5L5o3c/RMza15Ek07As+7uBMWNaptZA3dfmMz7SdkW78qlautWc9n4VzcWAbo5LAK0ISwCpHsdRNKrxEnC3Y9KRSBAI2BuxPN54WtKEuVMvN7D4bO/YMDoR2iy7DeG7dGBwe0v2KQIkIaYRNIvncNNxYl11VTMSXIz6w5Jjy93AAARFElEQVR0B2jadPMbp6TsinVzXONlv9Fv9KMcNXsC32/bdJMiQIWUIEQyoywliXlAk4jnjYEFsRq6+6PAowB5eXm62ipLRA8xFVcEqJDmIEQypywliRFATzN7mWDCepnmI8qHWPc9tJszlYGjHqbFkvm8u+vBDDriYhbWqr9JGwO6qgchklHpvAT2JaA9UM/M5gH9gFwAdx8CvEtw+etsgktgz09XbJI6B9w6it9WrN34fLsVf9D3o8c58btPmVO7AeedPoCPd9p3k300OS1SdqQtSbj7WcVsd+DyNIUjadD1sfEbE0TOhvV0m/QWvT57IW4RoF22q8Gof7fPULQiEktZGm6SciTyzunIIkAf7ZRH/yMv2awIUPXcSkoQImWQkoSUusJLXOuuWsb1Y57i9G/jFwEqNLhzmxhHEpFMU5KQUjV8ynxe/N9PdP1qJL0/eZYaa/N56MDTuP+gzYsAAVSuZPz39L20/pJIGaUkIVsssrTongtn8cYHD7HXr7P4X9M29D2qBz/Ua7LZPlrSWyQ7KElI0iLvnK61eiU3f/IsXae8x+IateMWAQK458y2Sg4iWUJJQpJy1F1jmfX7X+DOqd9+xPVjgyJAT+97Incf2nVjEaBo7VrUVYIQySJKElIikTfGtVw0h0EfPMT+86YzqeFunHvGIKZvv1PcfXfZrobufxDJMkoSkpDIJTVqrFlFr3Evcv7EEZsUAXKrFHd/3QMhkp2UJCSuzZbzduf47z6j70ePscPKJZsUASpKrao5ShAiWUpJQjYTq9ZDdBGgHiffsLEIUFHUgxDJbkoSsonotZY2KQJUuepmRYBi0dpLIuWHkoQAEVcrRThi9gT6j36UJst+4/U9OjC4wwUsrlEnzhHUaxApj5QkKrhYVeKiiwCdedZgJjQterluFQUSKZ+UJCqw6DoPVQrWcfEXw+g5/hU2mPGf9ufzZF4nCnLi/zPRjXEi5ZuSRAUUa2I6kSJAm7TXvINIhaAkUcFET0xvv2IxN330BCd+9yk/1WnAuacP4JOoIkCFlBhEKh4liQqkTb/3Wb5mPQCV1xdw3qS3+Ne4F8ldX8Bdh3TlkQNO3aQIUKHta1Zhwo1HpTtcESkDlCQqiMgEkTdvGoM+eLjIIkAAlQzuOkNzDiIVmZJEObfbje+yer0DsO1fS+kz9mlO/3Y082rVp/spN/LBLgdqpVYRiUtJopyKnJyutGE9Z381kt4fP0P1dWuKLgJkMHvw8ekOV0TKKCWJcihyaGnPhbO4JYEiQKC5BxHZXFqThJkdA9wL5ACPu/ttUdu7AXcC88OXHnD3x9MZYzaL7D3UWr2S3psUAerNiFb/iDm0VKtqDl8POCbd4YpIFkhbkjCzHOBB4ChgHvClmY1w9+lRTYe6e890xVVebEwQJSwCpAQhIkVJZ09if2C2u/8IYGYvA52A6CQhJVS47lJJiwBpeElEipPOJNEImBvxfB5wQIx2p5rZP4DvgX+5+9wYbYS/k0ONNau4MaIIUO9jr+S1PeMXATLgbl29JCIJSGeS2HwwHDzq+VvAS+6+xswuBZ4BDt/sQGbdge4ATZs2Le04y7yNay5FFAHabuWfvLxXR+447NwiiwDNuU1XLolI4tKZJOYBkZfVNAYWRDZw9z8inj4G3B7rQO7+KPAoQF5eXnSiKdea93kHgJ3+mMeAUUM49OepfLt9Cy495UamNmwZd79qOcZ3tx6XrjBFpJxIZ5L4EtjFzHYkuHqpC3B2ZAMza+DuC8OnJwEz0hhfmVa45lK1dau5fPyrXDLhdVbnVqXvUZfyQttjiywCpGW8RSRZaUsS7l5gZj2BkQSXwD7p7tPMbCAw0d1HAFea2UlAAbAE6Jau+Mqywt5DSYsAgYaXRGTLmHt2j9bk5eX5xIkTMx1GyjTv8w6Nl/5Kvw8f5ajZX/D9tk3pe3SPYosAgRKEiMRnZpPcPa+4drrjugwqrBZXpWAdl38xjCvGD2W9VeLW9hfwVN5JRRYBApURFZHSoyRRxhQmiMgiQO+0bMegwy/m11r1itxX9R5EpLQpSZQRhXdMb79iMQ989AQnJFAEqJDumhaRVFGSKAOa93mHyusLuDAsAlR5w3r+75CuPBqnCFAkLektIqmkJJFBhVct5c2bxi0jH2K3xT/zYYv96H/kJcytvUOR+2pJDRFJByWJNCtcSgOCIkDXj32K0779kHm16nNx55sYtfMBMVdqjaTeg4iki5JEmhROSEPJigBFMuAnXdYqImmkJJEGkQmizcLvGfTBw+z16yzGNWvDzUf14IdtYxcBiqTeg4hkgpJEihXOO9RavZJrP36Gs6e+z6Kt63DFib15K04RoGhKECKSKUoSKVI492C+gVO//Yg+Y5+iTv4Knso7ibsP6crKqlsVub+GlkSkLFCSSIHC3kNkEaCJjVpxzpk9mLFd/CJAhbSchoiUFUoSpajwhrgaa1bxr89eoNukt1iWQBGgQkoOIlLWKEmUkuZ93gF3TvjuU2766HG2W/knL7XtyB3/OI9l1WvG3U91HkSkLFOS2EKxigB9s30LLjnlRr4qoggQqM6DiJR9ShJJiLykNZkiQKAEISLZQUmihAp7DgBHzppA/9GP0Hj577ze+nAGtz9fRYBEpFxRkkjQbje+y+r1QYGmxkt/pf/oRzjyhy+ZWa8pZ5x9G180aV3sMZQgRCTbKEkUI3KtpSoF67g4iSJAmpwWkWylJBFH5LwDwCE/TWHA6CG0WDKft1sewi2HX1RsESBQ70FEspuSRAyR8w7br1hM34giQOecMZBPd9yn2GNoYlpEygMliQiRyaHy+gLOS6IIkJKDiJQnShJsmhwA9pv7LYM+eLhERYA0rCQi5VFak4SZHQPcC+QAj7v7bVHbqwLPAvsCfwBnuvucVMUTnRySKQLUrkVdXrj4oFSFKCKSUWlLEmaWAzwIHAXMA740sxHuPj2i2YXAn+6+s5l1AW4HzkxFPJEJIroI0IMHns4DB51ZbBEg9R5EpLxLZ09if2C2u/8IYGYvA52AyCTRCegfPn4NeMDMzN29NAM54NZRGx8nUwRIyUFEKop0JolGwNyI5/OAA+K1cfcCM1sGbAssLs1AfluxFoBzJr/NgFGPJFwESMlBRCqadCaJWN++0T2ERNpgZt2B7gBNmzZNOqBxzdryZN5J3FNMESBVhhORiiqdSWIeEDmO0xhYEKfNPDOrDGwDLIk+kLs/CjwKkJeXl/RQ1I/bNuaWIy4uso16DyJSkaUzSXwJ7GJmOwLzgS7A2VFtRgDnAeOB04CPSns+AmD7mlU2DjnFomU0REQCRZdKK0XuXgD0BEYCM4BX3H2amQ00s5PCZk8A25rZbODfQJ9UxDLhxqPYvubmN8W1a1GXObcdrwQhIhKyFPyinlZ5eXk+ceLETIchIpJVzGySu+cV1y5tPQkREck+ShIiIhKXkoSIiMSlJCEiInEpSYiISFxZf3WTmS0Cft6CQ9SjlJf9yKDyci7l5Tyg/JyLzqPs2dJzaebu9YtrlPVJYkuZ2cRELgPLBuXlXMrLeUD5ORedR9mTrnPRcJOIiMSlJCEiInEpSYQLBZYT5eVcyst5QPk5F51H2ZOWc6nwcxIiIhKfehIiIhKXkoSIiMRVYZKEmR1jZjPNbLaZbbYEuZlVNbOh4fYJZtY8/VEWL4Hz6GZmi8xsavhzUSbiLI6ZPWlmv5vZt3G2m5ndF57n12a2T7pjTFQC59LezJZFfCY3pzvGRJhZEzMbY2YzzGyamV0Vo02Z/1wSPI9s+UyqmdkXZvZVeC4DYrRJ7XeXu5f7HyAH+AHYCagCfAXsHtXmMmBI+LgLMDTTcSd5Ht2ABzIdawLn8g9gH+DbONuPA94jKGl7IDAh0zFvwbm0B97OdJwJnEcDYJ/wcU3g+xj/vsr855LgeWTLZ2LA1uHjXGACcGBUm5R+d1WUnsT+wGx3/9Hd1wIvA52i2nQCngkfvwYcYWaxam5nUiLnkRXc/RNilKaN0Al41gOfA7XNrEF6oiuZBM4lK7j7QnefHD5eQVAcLLq4e5n/XBI8j6wQ/j2vDJ/mhj/RVxul9LuroiSJRsDciOfz2PwfzcY2HlTRWwZsm5boEpfIeQCcGg4FvGZmTWJszwaJnmu2OCgcMnjPzPbIdDDFCYcs9ib4zTVSVn0uRZwHZMlnYmY5ZjYV+B0Y5e5xP5NUfHdVlCQRK6tGZ+NE2mRaIjG+BTR39zbAaP7+DSPbZMPnkajJBOvk7AXcDwzPcDxFMrOtgdeBXu6+PHpzjF3K5OdSzHlkzWfi7uvdvS3QGNjfzFpHNUnpZ1JRksQ8IPI36sbAgnhtzKwysA1lbwih2PNw9z/cfU349DFg3zTFVtoS+cyygrsvLxwycPd3gVwzq5fhsGIys1yCL9YX3H1YjCZZ8bkUdx7Z9JkUcvelwFjgmKhNKf3uqihJ4ktgFzPb0cyqEEzujIhqMwI4L3x8GvCRhzNBZUix5xE1PnwSwXhsNhoBnBteTXMgsMzdF2Y6qGSY2Q6FY8Rmtj/B/7s/MhvV5sIYnwBmuPtdcZqV+c8lkfPIos+kvpnVDh9XB44EvotqltLvrsqldaCyzN0LzKwnMJLgCqEn3X2amQ0EJrr7CIJ/VM+Z2WyCLNwlcxHHluB5XGlmJwEFBOfRLWMBF8HMXiK4wqSemc0D+hFMyuHuQ4B3Ca6kmQ2sAs7PTKTFS+BcTgN6mFkBkA90KYO/gAC0A84BvgnHwAFuAJpCVn0uiZxHtnwmDYBnzCyHIJG94u5vp/O7S8tyiIhIXBVluElERJKgJCEiInEpSYiISFxKEiIiEpeShIiIxKUkISIicSlJSIVjZiuLbyUioCQhIiJFUJKQCs3M/m1m34Y/vSJe72tm35nZKDN7ycyuKcExXw6LwEwws5/N7PjURC+SehViWQ6RWMxsX4JlJQ4gWElzgpl9TLDkyakES0xXJlgxdFIJDr0XMNzdzzSzQ4C7gHdKM/ZIZlbH3f9M1fGlYlOSkHLHzEYDO8TYdKO7vxnx/BDgDXf/K9xvGHAoQQ/7TXfPD19/K/xzJ+BGYBt3Py18rQbwELCWYIXOYUA9oLDM5HSgTmmeXwx3U8waXWbWEdjO3Z9LcSxSzihJSLnj7kcm2DRe9a6Yr7v7j8CFZvZaxMudgdfc/S0zG0pQKnOWu68Ot+9DUGa2cPnqAcBWBOVn+wF3AH2BgcD1BKVBRwItCRahy4na58pwHwd+Dt9vt3A47N6otr0IejHLCXpLlyf49yKykeYkpCL7BDjZzLYKewSnAJ8CnwEnhkXotwaKmlNozN+V2tYTDDU1DfetQfClfXe4vTtQHVhKULd4EfAL8H8EX/55wEvufj1BFbJto/cBehD0cq529/uAxcDz7v7fOG2fcfcbgKrAzOT/qqSiUk9CKix3n2xmTwNfhC897u5TAMxsBEEP4GdgIkFJyFjmESSKqQS/dO0FvEAw9FQL+I+7jwvb7g1cXlgUKkxAOwEF7r7SzPYD/he23cbdF5lZ9D5PAQ9HvH+bMM5Yx38aeMTMagKLy+hS2FLGaalwkRjMbOvwi3srgh5Hd4KEcStwFEFCGRz2Fh4AVhP0QC4BLnb3zX5rN7MTgbMJeh6fEAxVDQDOICgodQlBz6AywRDWh1H7fETQIziJoG7AYIJ5lVOB24Cdo9rWAToS1H1Y6e7XlubfkVQMShIiMZjZi8DuQDWCIZvBCe43H2ji7huSeM/n3P2cku4nkkpKEiIiEpcmrkVEJC4lCRERiUtJQkRE4lKSEBGRuJQkREQkLiUJERGJS0lCRETiUpIQEZG4lCRERCSu/wfWTdd01sYYwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sorted_pvals = sorted(-np.log10(pvals))\n",
    "\n",
    "plt.scatter(x = sorted(-np.log10(np.arange(1, pvals.shape[0] + 1)/pvals.shape[0])), y = sorted_pvals)\n",
    "plt.plot(sorted(-np.log10(np.arange(1, pvals.shape[0] + 1)/pvals.shape[0])), \n",
    "         sorted(-np.log10(np.arange(1, pvals.shape[0] + 1)/pvals.shape[0])), 'r-')\n",
    "plt.title(r'QQ-plot of p-values under $H_0$')\n",
    "plt.xlabel(r'$-\\log_{10}p_{expected}$')\n",
    "plt.ylabel(r'$-\\log_{10}p_{observed}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What does the p-value distribution look like when $H_0$ is rejected?\n",
    "\n",
    "Let's sample from a normal distribution with $\\mu=1$ this time, but test whether the mean of the sampled values is **greater** than that of the standard normal distribution. This is a **one-sided statistical test for the upper tail**, and we'll set a significance cutoff of $\\alpha=0.05$.\n",
    "\n",
    "$$H_0: \\mu_{sample}= 0$$\n",
    "$$H_A: \\mu_{sample}>0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEA9JREFUeJzt3X+s3XV9x/HnS2DClAUYBRt+eHGpRjSzamdQEoPgNgQn6MLERW0cW80C/lhctipLdMncMMvc4uZ0dTKrUxxOGZ0wsFaUOH+2rkNYZTbaSUdD61TEkWEK7/1xvndeu0/v/fb2nvs9t/f5SE7u+X7u95zvC9L01c/3Z6oKSZIO9KihA0iSJpMFIUlqsiAkSU0WhCSpyYKQJDVZEJKkJgtCktRkQUiSmiwISVLT0UMHOBwnn3xyTU1NDR1DkpaUbdu2fbuqVsy13pIuiKmpKbZu3Tp0DElaUpL8R5/13MUkSWqyICRJTRaEJKnJgpAkNVkQkqQmC0KS1GRBSJKaLAhJUpMFIUlqWtJXUkuTbGr9TYNsd9c1Fw+yXR15nEFIkposCElSkwUhSWqyICRJTRaEJKnJgpAkNVkQkqQmC0KS1GRBSJKaLAhJUpMFIUlqsiAkSU0WhCSpyYKQJDVZEJKkJgtCktRkQUiSmiwISVKTBSFJarIgJElNFoQkqcmCkCQ1HT10AGmcptbfNHQEaclyBiFJarIgJElNFoQkqcmCkCQ1ja0gkpyR5LYkO5LcleR13fhJSTYn+Xr388RuPEnekWRnkjuSPGNc2SRJcxvnDGI/8IaqejJwDnBlkrOB9cCWqloFbOmWAV4ArOpe64B3jTGbJGkOYyuIqtpTVV/p3j8A7ABOAy4BNnarbQQu7d5fAry/Rr4AnJBk5bjySZJmtyjHIJJMAU8HvgicWlV7YFQiwCndaqcB98z42O5u7MDvWpdka5Kt+/btG2dsSVrWxl4QSR4LfBR4fVV9f7ZVG2P1/waqNlTVmqpas2LFioWKKUk6wFgLIskxjMrhg1X1sW74vuldR93Pvd34buCMGR8/Hbh3nPkkSQc3zrOYArwX2FFVb5/xq03A2u79WuDGGeOv7M5mOge4f3pXlCRp8Y3zXkznAq8Avppkezf2JuAa4PokVwDfAi7rfnczcBGwE3gQeNUYs0mS5jC2gqiqz9I+rgBwQWP9Aq4cVx5J0qHxSmpJUpMFIUlqsiAkSU0WhCSpySfKSUeYoZ6it+uaiwfZrsbHGYQkqcmCkCQ1WRCSpCYLQpLUZEFIkposCElSkwUhSWqyICRJTRaEJKnJgpAkNVkQkqQmC0KS1GRBSJKaLAhJUpMFIUlqsiAkSU0WhCSpyYKQJDVZEJKkJgtCktRkQUiSmiwISVKTBSFJarIgJElNFoQkqcmCkCQ1HT10AC0PU+tvGjqCpEPkDEKS1GRBSJKaLAhJUtPYCiLJtUn2Jrlzxthbkvxnku3d66IZv3tjkp1J7k7yi+PKJUnqZ5wziPcBFzbG/7SqVnevmwGSnA1cDjyl+8xfJjlqjNkkSXMYW0FU1e3Ad3qufgnw4ap6qKq+CewEnjWubJKkuc1ZEElOWuBtXpXkjm4X1Ind2GnAPTPW2d2NtfKsS7I1ydZ9+/YtcDRJ0rQ+M4gvJvlIkouS5DC39y7gZ4DVwB7gT7rx1vdW6wuqakNVramqNStWrDjMOJKkg+lTEE8ENgCvAHYm+cMkT5zPxqrqvqp6uKoeAd7Dj3Yj7QbOmLHq6cC989mGJGlhzFkQNbK5ql4G/DqwFvhSks8kefahbCzJyhmLLwamz3DaBFye5NFJzgJWAV86lO+WJC2sOW+1keSngZczmkHcB7yG0V/oq4GPAGcd5HPXAecBJyfZDbwZOC/Jaka7j3YBrwaoqruSXA/8G7AfuLKqHj6c/zBJ0uHpcy+mzwMfAC6tqt0zxrcmeffBPtTNOA703lnWfyvw1h55JEmLoE9BPKmqDnbA+G0LnEeSNCH6HKT+RJITpheSnJjk1jFmkiRNgD4FsaKqvje9UFXfBU4ZXyRJ0iToUxAPJzlzeiHJ4znINQqSpCNHn2MQVwOfTfKZbvm5wLrxRZIkTYI5C6KqbknyDOAcRlc8/1ZVfXvsySRJg+r7yNFHM7rx3tHA2Ummb8YnSTpC9blQ7m3AS4G7gEe64QIsCEk6gvWZQVzK6FqIh8YdRpI0OfqcxfQN4JhxB5EkTZY+M4gHge1JtgD/N4uoqteOLZUkaXB9CmJT95IkLSN9TnPdmOQ44MyqunsRMkmSJkCfR47+ErAduKVbXp3EGYUkHeH6HKR+C6Mnv30PoKq2c5BnQEiSjhx9CmJ/Vd1/wJj3YpKkI1yfg9R3JvlV4Kgkq4DXAp8bbyxJ0tD6zCBeAzyF0Smu1wHfB14/zlCSpOH1OYvpQUZ3dL16/HEkSZOiz72YbqNxzKGqzh9LIknSROhzDOK3Z7w/FvhlYP944kiSJkWfXUzbDhj65xkPD5IkHaH67GI6acbio4BnAo8bWyJJ0kTos4tpG6NjEGG0a+mbwBXjDCVJGl6fXUxeNS1Jy1CfXUwvme33VfWxhYsjSZoUfXYxXQE8B/hUt/w84NPA/Yx2PVkQknQE6lMQBZxdVXsAkqwE3llVrxprMknSoPrcamNquhw69wFPHFMeSdKE6DOD+HSSWxndh6mAy4HbxppKkjS4PmcxXZXkxcBzu6ENVXXDeGNJkobWZwYB8BXggar6ZJKfTHJ8VT0wzmCSpGH1eeTobwB/D/xVN3Qa8A/jDCVJGl6fg9RXAucyeg4EVfV14JRxhpIkDa9PQTxUVT+cXkhyND0eOZrk2iR7k9w5Y+ykJJuTfL37eWI3niTvSLIzyR1JnjGf/xhJ0sLpUxCfSfIm4LgkPw98BPjHHp97H3DhAWPrgS1VtQrY0i0DvABY1b3WAe/q8f2SpDHqUxDrgX3AV4FXAzcDvzfXh6rqduA7BwxfAmzs3m8ELp0x/v4a+QJwQndBniRpILOexZTkKGBjVb0ceM8CbO/U6YvuqmpPkuljGacB98xYb3c3tgdJ0iBmnUFU1cPAiiQ/MeYcaW2+uWKyLsnWJFv37ds35liStHz1uQ5iF6OnyG0C/nt6sKrePo/t3ZdkZTd7WAns7cZ3A2fMWO904N7WF1TVBmADwJo1a+Y8WC5Jmp+DziCSfKB7+1Lg4926x894zccmYG33fi1w44zxV3ZnM50D3H/A/Z8kSYtsthnEM5M8HvgW8OeH+sVJrgPOA05Osht4M3ANcH2SK7rvvaxb/WbgImAn8CDgnWIlaWCzFcS7gVuAs4CtM8bD6PjAE2b74qp62UF+dUFj3WJ0QZ4kaUIcdBdTVb2jqp4M/E1VPWHG66yqmrUcJElL35zXQVTVby5GEEnSZOlzoZwkaRmyICRJTRaEJKnJgpAkNVkQkqSmvo8claRZTa2/abBt77rm4sG2fSRzBiFJarIgJElNFoQkqcmCkCQ1WRCSpCYLQpLUZEFIkposCElSkwUhSWqyICRJTRaEJKnJgpAkNVkQkqQmC0KS1GRBSJKaLAhJUpMFIUlq8olyy8iQT/yStPQ4g5AkNVkQkqQmC0KS1GRBSJKaLAhJUpMFIUlqsiAkSU0WhCSpyYKQJDUNciV1kl3AA8DDwP6qWpPkJODvgClgF/ArVfXdIfJJkoadQTyvqlZX1ZpueT2wpapWAVu6ZUnSQCZpF9MlwMbu/Ubg0gGzSNKyN1RBFPCJJNuSrOvGTq2qPQDdz1MGyiZJYri7uZ5bVfcmOQXYnORrfT/YFco6gDPPPHNc+SRp2RtkBlFV93Y/9wI3AM8C7kuyEqD7ufcgn91QVWuqas2KFSsWK7IkLTuLXhBJHpPk+On3wC8AdwKbgLXdamuBGxc7myTpR4bYxXQqcEOS6e1/qKpuSfJl4PokVwDfAi4bIJskqbPoBVFV3wCe1hj/L+CCxc4jSWqbpNNcJUkTxIKQJDVZEJKkJgtCktRkQUiSmiwISVKTBSFJarIgJElNFoQkqcmCkCQ1WRCSpCYLQpLUZEFIkposCElSkwUhSWqyICRJTRaEJKnJgpAkNQ3xTGpJWlBT628aZLu7rrl4kO0uFmcQkqQmC0KS1GRBSJKaLAhJUpMFIUlqsiAkSU2e5jqAoU7Jk6RD4QxCktRkQUiSmiwISVKTBSFJarIgJElNnsUkSfM05BmJi3GjQGcQkqQmC0KS1LRsdzF5sZokzW7iZhBJLkxyd5KdSdYPnUeSlquJKogkRwHvBF4AnA28LMnZw6aSpOVpogoCeBaws6q+UVU/BD4MXDJwJklaliatIE4D7pmxvLsbkyQtskk7SJ3GWP3YCsk6YF23+IMkd4891Y87Gfj2Im9zoZh9GEs5Oyzt/Eds9rztsL778X1WmrSC2A2cMWP5dODemStU1QZgw2KGminJ1qpaM9T2D4fZh7GUs8PSzm/2wzNpu5i+DKxKclaSnwAuBzYNnEmSlqWJmkFU1f4kVwG3AkcB11bVXQPHkqRlaaIKAqCqbgZuHjrHLAbbvbUAzD6MpZwdlnZ+sx+GVNXca0mSlp1JOwYhSZoQFsQ8JPnjJF9LckeSG5KcMHSmvpJcluSuJI8kWRJndyzV268kuTbJ3iR3Dp3lUCU5I8ltSXZ0f15eN3SmvpIcm+RLSf61y/77Q2c6VEmOSvIvST4+ZA4LYn42A0+tqp8F/h1448B5DsWdwEuA24cO0scSv/3K+4ALhw4xT/uBN1TVk4FzgCuX0P/3h4Dzq+ppwGrgwiTnDJzpUL0O2DF0CAtiHqrqE1W1v1v8AqPrNZaEqtpRVYt9ceHhWLK3X6mq24HvDJ1jPqpqT1V9pXv/AKO/rJbEXQ1q5Afd4jHda8kcbE1yOnAx8NdDZ7EgDt+vAf80dIgjmLdfGViSKeDpwBeHTdJft4tmO7AX2FxVSyY78GfA7wCPDB1k4k5znRRJPgk8rvGrq6vqxm6dqxlNxT+4mNnm0if7EjLn7Vc0PkkeC3wUeH1VfX/oPH1V1cPA6u744A1JnlpVE38sKMkLgb1VtS3JeUPnsSAOoqqeP9vvk6wFXghcUBN2rvBc2ZeYOW+/ovFIcgyjcvhgVX1s6DzzUVXfS/JpRseCJr4ggHOBFyW5CDgW+Kkkf1tVLx8ijLuY5iHJhcDvAi+qqgeHznOE8/YrA0gS4L3Ajqp6+9B5DkWSFdNnFiY5Dng+8LVhU/VTVW+sqtOraorRn/VPDVUOYEHM118AxwObk2xP8u6hA/WV5MVJdgPPBm5KcuvQmWbTnQwwffuVHcD1S+X2K0muAz4PPCnJ7iRXDJ3pEJwLvAI4v/szvr37V+1SsBK4LckdjP6BsbmqBj1ddKnySmpJUpMzCElSkwUhSWqyICRJTRaEJKnJgpAkNVkQkqQmC0KS1GRBSAsoyc91zwk5NsljuucRPHXoXNJ8eKGctMCS/AGj++gcB+yuqj8aOJI0LxaEtMC6e0Z9Gfgf4DndnUWlJcddTNLCOwl4LKP7dR07cBZp3pxBSAssySZGT747C1hZVVcNHEmaF58HIS2gJK8E9lfVh7rnaX8uyflV9amhs0mHyhmEJKnJYxCSpCYLQpLUZEFIkposCElSkwUhSWqyICRJTRaEJKnJgpAkNf0vF2ADHDLB8zgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 1000 samples \n",
    "sample2 = stats.norm.rvs(loc = 1, scale = 1, size = 1000)\n",
    "plt.hist(sample2)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('frequency')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Compute p-values\n",
    "\n",
    "This time we have to use the survival function because we are testing for how likely we are to observe a value _greater_ than the sampled value under the standard normal distribution - by definition, an upper tail test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAElCAYAAADqeCmyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHIlJREFUeJzt3X2cXVV97/HP1xACyEOADBSSwCAE5eGlkQ4P99LbUsJVCEqwBQxFibxooxZbuVI1cPu6xlbuxbaaXqqisSBReYpgIQK2IpAitgkdMKRA4BogkiEpGcgDIAVJ+N0/9ho5mZwzZ895mJNZ+b5fr/Oavddee+/fOgO/2Vnn7N9WRGBmZvl6S6cDMDOz9nKiNzPLnBO9mVnmnOjNzDLnRG9mljknejOzzDnR23ZP0kmS+kboXIsl/WFaPk/Sj1p47EclnZSW50r6bguPfZmkv2/V8SwvO3U6ALPtVURcB1xXr5+ka4G+iPjzOsc7qhVxpT8W342ISRXH/t+tOLblyVf0Zm0myRdU1lFO9NYwSaskXSrpMUkbJH1L0i41+s6RdPOgtv8r6cq0fIGkFZJekvSUpI8Ocd6QdFjF+rWSvlCx/j5JyyRtlPQvkt45xLH+u6THJW2S9BVAFds+Iun+tCxJ8yStS32XSzpa0mzgPOAzkl6W9IOK9+azkpYDv5S0U2o7peL0u0i6KY35IUnvqjdGSW8FfggcmM73sqQDB08FSTojTRVtTNNRR1RsWyXpz9IYNqUYqv7eLA9O9Nas84D3AocChwO1pi9uAKZL2hNA0hjgHOD6tH0d8D5gT+ACYJ6kY4YbTNrnGuCjwL7AN4BFksZV6TsBuCXFPAF4EjixxqHfA/w2xRjHAx8EXoiI+RTTO38VEbtHxPsr9jkXOB0YHxGbqxxzBvA9YB+K9+FWSWOHGl9E/BI4DViTzrd7RKwZNK7DKd7vi4Eu4E7gB5J2ruh2DnAqcAjwTuAjQ53XRjcnemvWVyJidUSsBy6nSG7biIhfAA8BZ6amk4FXImJJ2n5HRDwZhX8GfgT8twbi+SPgGxGxNCK2RMQC4DXghCp9pwOPRcTNEfE68LfAf9Q47uvAHsA7AEXEiohYWyeWK9N78581tj9Yce4vA7vUiHO4PgjcERF3pWP/DbAr8F8HxbYm/d5+AExtwXltO+VEb81aXbH8C+BAAEk/rJhaOC9tv543/xD8AW9ezSPpNElLJK2XtJEiCU9oIJ6DgUvSlMXGdKzJA3ENcmBl/FFU+FtdpR8RcQ/wFeCrwHOS5g/862QIVY9VbXtEvAH01YhzuA6k+F1UHns1MLGiT+UftFeA3VtwXttOOdFbsyZXLB8ErAGIiNMqphYGvrnyPeAkSZOAD5ASfZpWuYXiynP/iBhPMd0gqnsF2K1i/TcqllcDl0fE+IrXbhFxQ5XjrK2MX5IGjWcrEXFlRPwmcBTFFM6nBzbV2qXWsZLKc78FmER6/xh6jPWOu4biD97AsQfG9Wyd/SxTTvTWrIskTZK0D3AZcFOtjhHRDywGvgU8HREr0qadgXFAP7BZ0mkUc+K1LAP+QNIYSacCv1Ox7ZvAxyQdnz5Afauk0yXtUeU4dwBHSfq99M2YP2XrhPprko5NxxwL/BJ4FdiSNj8HvG2IeGv5zYpzX0wxxbSkxBifA/aVtFeN4y4ETpc0LcV7STr2vzQQo2XAid6adT3FfPpT6fWFobtzPXAKFdM2EfESRZJdCGygmNZZNMQxPgm8H9hI8WHwrRXH6qWYp/9KOtZKanzQGBHPA2cDVwAvAFOAn9Y4554Uf0Q2UEyLvEDxLxCAq4Ej01TRrTX2r+Y2ivn0DcCHgd9Lc+r1xvg4xYetT6VzbjXdExFPAB8C/g54Ph3n/RHxq2HEZhmRHzxijZK0CvjDiPhxp2Mxs9p8RW9mljknejOzzHnqxswsc76iNzPLnBO9dYwqyvZW2dZwaeJm9rXy0k1xs9Lyr+sC2fbHVfWsY1pVtne4JAUwJSJWNnmcucBhEfGhlgTWYcN9XyLitDaHZC3iK3qzHYxcNnmH40S/g5J0RCpfuzFNoZxRse1aSV+VdEcqobtU0qEV298h6a5Ul+YJSecMcZ4Jkm5P51kv6Sfpdv+BcrmnpOVd03k3SHoMOHbQcQ6UdIukfklPS/rTim1D7jvoOPelxYdTHZ4PpvaapY1VlBt+Nr0XT6Q7Tk+luBP4g+k4D5d4z7eZUhr0HsyVdLNqly4esix0nTFsUza53vsiae/0u+tP57tdRfmKgX1+/TQu285FhF872AsYS3HH6GUU5QdOBl4C3p62XwusB46jmN67DrgxbXsrRT2ZC9K2Yyjuvjyqxrn+D/D1dM6xFBUpB77ttQo4JS1fAfyEomTvZOARiqc2QXFB8iDwv1K8b6O4C/e99fatEVNQTLkMrB9DUSb5eGAMMCvFNg54exrvgalvN3BoWp5L8aSnsu/7SYPjGvQezKWoknlWeq/+DHgaGFvR95E0xn0o7uL9Qr0xVOy7LO27a8n3ZV/g9ylq7uxBUavo1ortiylumIPi7uP7O/3ftl/VX76i3zGdQFGt8IqI+FUUlRlvZ+sSw9+PiAeiqKN+HW+WsX0fsCoivhURmyPiIYqCZGfVONfrwAHAwRHxekT8JFJmGOQcimJk6yNiNXBlxbZjga6I+IsU71MU5Qhmlti3jKFKG2+hSPhHShobEasi4slhHn846pUurlUWukx55nplk7cSES9ExC0R8UoUZSouZ+uaOzZKONHvmA4EVkdRvnbALyhXxvZg4HhtXQb4POA3JB2kN0sTv5z6/zXFvx5+pOLJUXOGimlQPAMOpniiUuU5LwP2L7FvGTVLG0fxweTFFFfb6yTdOLi2TIvVK11ctSz0UGOosW9dknaT9A1Jv5D0InAfMF7FQ2NsFHGi3zGtASYPzJUnB1GujO1q4J9j6zLAu0fExyPimXizNPHuUBQsi4hLIuJtFMW1PiVpWpXjblUyOMVTec6nB51zj4iYXmLfMoYsbRwR10fEb1Ek0wC+mPYb7t2Gv6Si9HBKmF2D+gxVunir7VSUha43hgbjvYRi6ur4iNiT4glbULt8tG2nnOh3TEspks5nJI1V8V329wM3ltj3duBwSR9O+45VUcL3iGqd0weEh0kS8CLFVMiWKl0XApemDwAnAX9Sse0B4MX0YeKuKkr3Hi3p2BL7VjO4rHDN0saS3i7pZBU1818F/pOtyxN3V/7BTB+oLq5x3v9H8ZzY01WUD/5zimmhSkOVLobaZaGHU5657PuyRxrvxnS+zw3jWLYdcaLfAUVRrvYMimePPg98DTg/ivK39fZ9iaJW/EyKq8n/oLjC3eaZrMkU4MfAy8C/Al+LiMVV+n2eYiriaYqyx9+pOOcWij9EU9P254G/B/aqt28Nc4EFaYrjnBi6tPE4ig97n09j3Y8iwULx4STAC5IeSsuTqVHqOCI2AX+cYn+W4o/t4Bu7hipdDDXKQtcZQ1lzqXhfKB6tuCvF2JcA/zjM49l2wrVuzFpI0jJgWkS80MC+cxniBiy5LLQ1yDdOmLVQRPgh27bd8dSNmVnmPHVjZpY5X9GbmWVuu5ijnzBhQnR3d3c6DDOzUeXBBx98PiIG34uxje0i0Xd3d9Pb29vpMMzMRhVJpe4C99SNmVnmnOjNzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5raLO2Ob0T3njo6de9UVp3fs3GZmZZW+ok+Pb/uZpNvT+iGSlkr6uaSbJO2c2sel9ZVpe3d7QjczszKGM3XzSWBFxfoXgXkRMYXi0WUXpvYLgQ0RcRgwjzcfpGxmZh1QKtGnBy6fTvGsS9KDnk8Gbk5dFgBnpuUZaZ20fVrqb2ZmHVD2iv5vgc8Ab6T1fYGNEbE5rfcBE9PyRGA1QNq+KfXfiqTZknol9fb39zcYvpmZ1VM30Ut6H7AuIh6sbK7SNUpse7MhYn5E9ERET1dX3XLKZmbWoDLfujkROEPSdGAXYE+KK/zxknZKV+2TgDWpfx8wGeiTtBOwF7C+5ZGbmVkpda/oI+LSiJgUEd3ATOCeiDgPuBc4K3WbBdyWlhelddL2e8IPpjUz65hmbpj6LPApSSsp5uCvTu1XA/um9k8Bc5oL0czMmjGsG6YiYjGwOC0/BRxXpc+rwNktiM3MzFrAJRDMzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWebKPBx8F0kPSHpY0qOSPp/ar5X0tKRl6TU1tUvSlZJWSlou6Zh2D8LMzGor84Sp14CTI+JlSWOB+yX9MG37dETcPKj/acCU9DoeuCr9NDOzDijzcPCIiJfT6tj0Guph3zOAb6f9lgDjJR3QfKhmZtaIUnP0ksZIWgasA+6KiKVp0+VpemaepHGpbSKwumL3vtQ2+JizJfVK6u3v729iCGZmNpRSiT4itkTEVGAScJyko4FLgXcAxwL7AJ9N3VXtEFWOOT8ieiKip6urq6HgzcysvmF96yYiNgKLgVMjYm2annkN+BZwXOrWB0yu2G0SsKYFsZqZWQPKfOumS9L4tLwrcArw+MC8uyQBZwKPpF0WAeenb9+cAGyKiLVtid7MzOoq862bA4AFksZQ/GFYGBG3S7pHUhfFVM0y4GOp/53AdGAl8ApwQevDNjOzsuom+ohYDry7SvvJNfoHcFHzoZmZWSv4zlgzs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDLnRG9mlrkyjxLcRdIDkh6W9Kikz6f2QyQtlfRzSTdJ2jm1j0vrK9P27vYOwczMhlLmiv414OSIeBcwFTg1PQv2i8C8iJgCbAAuTP0vBDZExGHAvNTPzMw6pG6ij8LLaXVsegVwMnBzal9A8YBwgBlpnbR9WnqAuJmZdUCpOXpJYyQtA9YBdwFPAhsjYnPq0gdMTMsTgdUAafsmYN9WBm1mZuWVSvQRsSUipgKTgOOAI6p1Sz+rXb3H4AZJsyX1Surt7+8vG6+ZmQ3TsL51ExEbgcXACcB4STulTZOANWm5D5gMkLbvBayvcqz5EdETET1dXV2NRW9mZnWV+dZNl6TxaXlX4BRgBXAvcFbqNgu4LS0vSuuk7fdExDZX9GZmNjJ2qt+FA4AFksZQ/GFYGBG3S3oMuFHSF4CfAVen/lcD35G0kuJKfmYb4jYzs5LqJvqIWA68u0r7UxTz9YPbXwXObkl0ZmbWNN8Za2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnmnOjNzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpa5Ms+MnSzpXkkrJD0q6ZOpfa6kZyUtS6/pFftcKmmlpCckvbedAzAzs6GVeWbsZuCSiHhI0h7Ag5LuStvmRcTfVHaWdCTFc2KPAg4Efizp8IjY0srAzcysnLpX9BGxNiIeSssvASuAiUPsMgO4MSJei4ingZVUebasmZmNjGHN0UvqpnhQ+NLU9AlJyyVdI2nv1DYRWF2xWx9V/jBImi2pV1Jvf3//sAM3M7NySid6SbsDtwAXR8SLwFXAocBUYC3wpYGuVXaPbRoi5kdET0T0dHV1DTtwMzMrp1SilzSWIslfFxHfB4iI5yJiS0S8AXyTN6dn+oDJFbtPAta0LmQzMxuOMt+6EXA1sCIivlzRfkBFtw8Aj6TlRcBMSeMkHQJMAR5oXchmZjYcZb51cyLwYeDfJS1LbZcB50qaSjEtswr4KEBEPCppIfAYxTd2LvI3bszMOqduoo+I+6k+737nEPtcDlzeRFxmZtYivjPWzCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWWuzKMEJ0u6V9IKSY9K+mRq30fSXZJ+nn7undol6UpJKyUtl3RMuwdhZma1lbmi3wxcEhFHACcAF0k6EpgD3B0RU4C70zrAaRTPiZ0CzAauannUZmZWWt1EHxFrI+KhtPwSsAKYCMwAFqRuC4Az0/IM4NtRWAKMH/QgcTMzG0HDmqOX1A28G1gK7B8Ra6H4YwDsl7pNBFZX7NaX2gYfa7akXkm9/f39w4/czMxKKZ3oJe0O3AJcHBEvDtW1Slts0xAxPyJ6IqKnq6urbBhmZjZMpRK9pLEUSf66iPh+an5uYEom/VyX2vuAyRW7TwLWtCZcMzMbrjLfuhFwNbAiIr5csWkRMCstzwJuq2g/P3375gRg08AUj5mZjbydSvQ5Efgw8O+SlqW2y4ArgIWSLgSeAc5O2+4EpgMrgVeAC1oasZmZDUvdRB8R91N93h1gWpX+AVzUZFxmZtYivjPWzCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWuTL16K2G7jl3dOS8q644vSPnNbPRyVf0ZmaZc6I3M8tcmWfGXiNpnaRHKtrmSnpW0rL0ml6x7VJJKyU9Iem97QrczMzKKXNFfy1wapX2eRExNb3uBJB0JDATOCrt8zVJY1oVrJmZDV/dRB8R9wHrSx5vBnBjRLwWEU9TPCD8uCbiMzOzJjUzR/8JScvT1M7eqW0isLqiT19q24ak2ZJ6JfX29/c3EYaZmQ2l0UR/FXAoMBVYC3wptatK36h2gIiYHxE9EdHT1dXVYBhmZlZPQ4k+Ip6LiC0R8QbwTd6cnukDJld0nQSsaS5EMzNrRkOJXtIBFasfAAa+kbMImClpnKRDgCnAA82FaGZmzah7Z6ykG4CTgAmS+oDPASdJmkoxLbMK+ChARDwqaSHwGLAZuCgitrQndDMzK6Nuoo+Ic6s0Xz1E/8uBy5sJyszMWsd3xpqZZc6J3swsc070ZmaZc6I3M8ucE72ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXOid7MLHNO9GZmmXOiNzPLnBO9mVnm6lavtO1P95w7OnbuVVec3rFzm1ljfEVvZpY5J3ozs8zVTfSSrpG0TtIjFW37SLpL0s/Tz71TuyRdKWmlpOWSjmln8GZmVl+ZK/prgVMHtc0B7o6IKcDdaR3gNIrnxE4BZgNXtSZMMzNrVN1EHxH3AesHNc8AFqTlBcCZFe3fjsISYPygB4mbmdkIa3SOfv+IWAuQfu6X2icCqyv69aW2bUiaLalXUm9/f3+DYZiZWT2t/jBWVdqiWseImB8RPRHR09XV1eIwzMxsQKOJ/rmBKZn0c11q7wMmV/SbBKxpPDwzM2tWo4l+ETArLc8CbqtoPz99++YEYNPAFI+ZmXVG3TtjJd0AnARMkNQHfA64Algo6ULgGeDs1P1OYDqwEngFuKANMZuZ2TDUTfQRcW6NTdOq9A3gomaDMjOz1vGdsWZmmXOiNzPLnBO9mVnmXKbYhqVTJZJdHtmscb6iNzPLnK/obVTww1bMGucrejOzzDnRm5llzonezCxzTvRmZplzojczy5wTvZlZ5pzozcwy50RvZpY5J3ozs8w50ZuZZa6pEgiSVgEvAVuAzRHRI2kf4CagG1gFnBMRG5oL06xzXMjNRrtWXNH/bkRMjYietD4HuDsipgB3p3UzM+uQdkzdzAAWpOUFwJltOIeZmZXUbKIP4EeSHpQ0O7XtHxFrAdLP/artKGm2pF5Jvf39/U2GYWZmtTRbpvjEiFgjaT/gLkmPl90xIuYD8wF6enqiyTjMzKyGpq7oI2JN+rkO+AfgOOA5SQcApJ/rmg3SzMwa13Cil/RWSXsMLAPvAR4BFgGzUrdZwG3NBmlmZo1rZupmf+AfJA0c5/qI+EdJ/wYslHQh8AxwdvNhmplZoxpO9BHxFPCuKu0vANOaCcrMzFrHz4w12075Ri1rFZdAMDPLnBO9mVnmnOjNzDLnRG9mljknejOzzDnRm5llzonezCxz/h69mW2lU9/fB3+Hv12c6M1su+GbxNrDUzdmZpnzFb2Z7fByn67yFb2ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWWubYle0qmSnpC0UtKcdp3HzMyG1pZEL2kM8FXgNOBI4FxJR7bjXGZmNrR2XdEfB6yMiKci4lfAjcCMNp3LzMyG0K4bpiYCqyvW+4DjKztImg3MTqsvS3qiwXNNAJ5vcN/RzOPeceyIY4YdZNz64larwx3zwWU6tSvRq0pbbLUSMR+Y3/SJpN6I6Gn2OKONx73j2BHHDDvmuNs15nZN3fQBkyvWJwFr2nQuMzMbQrsS/b8BUyQdImlnYCawqE3nMjOzIbRl6iYiNkv6BPBPwBjgmoh4tB3nogXTP6OUx73j2BHHDDvmuNsyZkVE/V5mZjZq+c5YM7PMOdGbmWVu1CT6eiUVJI2TdFPavlRS98hH2Xolxv0pSY9JWi7pbkmlvle7PStbPkPSWZJCUhZfwSszbknnpN/3o5KuH+kYW63Ef98HSbpX0s/Sf+PTOxFnK0m6RtI6SY/U2C5JV6b3ZLmkY5o+aURs9y+KD3SfBN4G7Aw8DBw5qM8fA19PyzOBmzod9wiN+3eB3dLyx0f7uMuMOfXbA7gPWAL0dDruEfpdTwF+Buyd1vfrdNwjMOb5wMfT8pHAqk7H3YJx/zZwDPBIje3TgR9S3I90ArC02XOOliv6MiUVZgAL0vLNwDRJ1W7cGk3qjjsi7o2IV9LqEop7FkazsuUz/hL4K+DVkQyujcqM+4+Ar0bEBoCIWDfCMbZamTEHsGda3osM7seJiPuA9UN0mQF8OwpLgPGSDmjmnKMl0VcrqTCxVp+I2AxsAvYdkejap8y4K11IcSUwmtUds6R3A5Mj4vaRDKzNyvyuDwcOl/RTSUsknTpi0bVHmTHPBT4kqQ+4E/iTkQmto4b7/31do+Xh4HVLKpTsM9qUHpOkDwE9wO+0NaL2G3LMkt4CzAM+MlIBjZAyv+udKKZvTqL4l9tPJB0dERvbHFu7lBnzucC1EfElSf8F+E4a8xvtD69jWp7LRssVfZmSCr/uI2knin/mDfXPo9GgVCkJSacA/xM4IyJeG6HY2qXemPcAjgYWS1pFMYe5KIMPZMv+N35bRLweEU8DT1Ak/tGqzJgvBBYCRMS/ArtQFP7KWctLyIyWRF+mpMIiYFZaPgu4J9InG6NY3XGnaYxvUCT50T5nC3XGHBGbImJCRHRHRDfF5xJnRERvZ8JtmTL/jd9K8eE7kiZQTOU8NaJRtlaZMT8DTAOQdARFou8f0ShH3iLg/PTtmxOATRGxtpkDjoqpm6hRUkHSXwC9EbEIuJrin3UrKa7kZ3Yu4tYoOe6/BnYHvpc+e34mIs7oWNBNKjnm7JQc9z8B75H0GLAF+HREvNC5qJtTcsyXAN+U9D8opi8+Mtov4CTdQDH9NiF99vA5YCxARHyd4rOI6cBK4BXggqbPOcrfMzMzq2O0TN2YmVmDnOjNzDLnRG9mljknejOzzDnRm5llzonezCxzTvRmZplzojerQlK3pMclLUg1wW+WtFun4zJrhBO9WW1vB+ZHxDuBFymeeWA26jjRm9W2OiJ+mpa/C/xWJ4Mxa5QTvVltg+uDuF6IjUpO9Ga1HZRqoENRF/3+TgZj1ignerPaVgCzJC0H9gGu6nA8Zg0ZFWWKzTrkjYj4WKeDMGuWr+jNzDLnevRmZpnzFb2ZWeac6M3MMudEb2aWOSd6M7PMOdGbmWXu/wPw7mx1CaKNlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pvals2 = stats.norm.sf(sample2, loc = 0, scale = 1)\n",
    "plt.hist(pvals2)\n",
    "plt.xlabel('p')\n",
    "plt.title('p-value distribution\\none-sided test, upper tail')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### QQ-plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEbCAYAAADTZlM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVNX5x/HP47oIAooKFlDE3kURK8ZgxYrYscVe0ow/E2NNLFGxJBpjYlaw0wRREStibLERQVDsYgdUQKXJ0pbn98e5q+uyy87cnZk7M/f7fr325e7Mnb3P3cHznDnn3OeYuyMiIumzQtIBiIhIMpQARERSSglARCSllABERFJKCUBEJKWUAEREUkoJQMqGmXUxMzezFQtwrlZm9oiZzTaz+/N9viZi+dTM9kkyBilNSgCSF2Z2splNMrP5ZvaVmd1qZqvWO2ZLMxsVNaJzzewZM9ulQPFdbmaDmvErjgTWAtZw96NyFJZIQSkBSM6Z2e+B64DzgVWBXYAuwFNmVhkdsxHwEjAJ2ADoCIwExpjZTgmEna31gQ/cfUnSgYjE5u760lfOvoBVgHnA0fUebwNMB06Kfh4IPN7A6/8NPLuc3/8c0A/4HzAbeBhYPXquC+DAitHPHYFRwLfAZOCM6PH9gUXA4ijWNxo51xbR+WYBbwO9o8evqPf60xp47eXACGAYMBd4HejayHk6AtW11xE9tj0wE6gENgKeAb6JHhsMtKtz7KfAPtH3dwNX1XmuJzCl3rkeAGYAnwDn1HluJ2AcMAf4Grgx6X9P+srvlz4BSK7tBrQEHqz7oLvPA54A9ose2hdoaOx8OPAzM2u5nHP8AjiV0JgtAf7RyHFDgSnRcUcC15jZ3u7+JHANMMzd27h71/ovjD6pPAI8BawJ/BYYbGabuftl9V5/RyPnPzS6xtWBIcDI2k9Adbn7NOAV4Ig6Dx8HjHD3xYARkl5HQlJaj5BgsmJmK0TX9AbQCdgbONfMekWH3Azc7O6rEJLO8GzPIaVFCUByrT0w0xseGvkS6FDnuC8bOaaC0Gg2ZqC7v+Xu3wN/Ao42s4q6B5jZesDuwAXuvsDdJwK3AydmeB27ED61XOvui9z9GeBR4NgMXw8w3t1rG/EbCYmxsTmOIbW/28wM6Bs9hrtPdvcx7r7Q3WdEv+vnWcRRa0egg7tfGV3Tx8CA6FwQPtFsbGbt3X2eu78a4xxSQpQAJNdmAu0bWYmzDmHoofa4dRo5xoGZZlZlZvOir4vrHPNFne8/IwyTtK/3ezoC37r73HrHdsrwOjoCX7j70piv/0mc0e+ZAnQ0s+PrXNcT0SEjgF3NrCOwB+Fv8F8AM1vTzO4zs6lmNgcYxLLXm4n1o/PPqv0CLiZMZgOcBmwKvGdmr5nZwTHOISVECUBy7RVgIXB43QfNrDVwAPB89NDTQEOrZ44GXo16qGdHQyxt3P2aOsesV+f7zoSe68x6v2casLqZta137NTo+6bK4E4D1ouGTRp6fSZ+iDP6PesC09x9cJ3rOgDA3WcRhpuOJgz/DHX32hj7RfFuGw3PnEAYFmrI98DKdX5eu873XwCfuHu7Ol9t3f3AKIYP3f1YwpDXdcCI6H2TMqUEIDnl7rMJk6S3mNn+ZlZpZl0IY+G1E5hEx+xmZleb2epm1tbMfgucAlzWxGlOiJaQrgxcSRgrr6kXxxfAy0A/M2tpZtsSeri15/8a6FKvga9rLKEx/WN0DT2BQ4D7MvtLALCDmR0efRo6l5AYlzesMoQwv3FE9H2ttoTJ5llm1omwuqoxE4EDo7/p2tF5a/0PmGNmF0T3MVSY2dZmtiOAmZ1gZh2iTyuzotfUIGVLCUByzt2vJwwt/JWwAuYTQq90n2jcHnf/kDBG35WwimUW8BfgMHcf08QpBhJWu3xFGFc/p5HjjiWsDJoGPARcVud3105Af2NmrzdwDYuA3oRPLTOBW4FfuPt7TcRW18PAMcB3hLmHw6P5gMaMAjYBvnb3N+o8fgXQjbDq6THqTbDXM5Awyfsp4RPFsDrXVENIYtsR3pOZhHmR2vsz9gfeNrN5hAnhvu6+IJMLldJkP37KFMkPMzuV0Ij1cPfPGzlmXULv+LLlrKrBzJ4DBrn77fmINVfM7HJgY3c/IelYRBqT91vmRdz9TjNbTFgi2mACcPcpZnYAcKiZtYmWjYpIHikBSEG4+8AMjplEuDNYRApAQ0AiIimlSWARkZQq6BCQmX1KWBVSAyxx9+6FPL+IiPwoiTmAPd29/k07DWrfvr136dIlz+GIiJSP8ePHz3T3Dk0fWeSTwF26dGHcuHFJhyEiUjLM7LNMjy30HIATasKPN7MzGzrAzM40s3FmNm7GjBkNHSIiIjlQ6ATQw927Ee6u/LWZ7VH/AHfv7+7d3b17hw4ZfYoREZEYCpoAorrnuPt0wq35pbDzk4hIWSpYAjCz1rWVGaMKg/sBbxXq/CIi8lOFnAReC3go7HXBisCQaGcmERFJQMESQLT70DJb74mICIycMJUbRr/PtFnVdGzXivN7bUaf7bPZfyh7Rb0MVEQkDUZOmMpFD06ienHYfmHqrGouejCUxcpnElApCBGRhN0w+v0fGv9a1YtruGH0+3k9rxKAiEjCps2qzurxXFECEBFJWMd2rbJ6PFeUAEREEnZ+r81oVVnxk8daVVZwfq/N8npeTQKLiCSsdqJXq4BERFKoz/ad8t7g16chIBGRlFICEBFJKSUAEZGUUgIQEUkpJQARkZRSAhARSSklABGRlFICEBFJKSUAEZGUUgIQEUkpJQARkZRSAhARSSklABGRlFICEBFJKSUAEZGUUgIQEUkpbQgjIlLHyAlTC74zV1KUAEREIiMnTOWiBydRvbgGgKmzqrnowUkAZZkENAQkIhK5YfT7PzT+taoX13DD6PcTiii/lABERCLTZlVn9XipUwIQEYl0bNcqq8dLnRKAiEjk/F6b0aqy4iePtaqs4PxemyUUUX5pElhEykIuVu/UHq9VQCIiCcqmQc/l6p0+23cq2wa/PiUAESmoTBr2bBv05a3eSUtjHkfB5wDMrMLMJpjZo4U+t4gkq7ZhnzqrGufHhn3khKk/OS7b5ZhpW72TK0lMAv8OeDeB84pIAYycMJUe1z7DBhc+Ro9rn/lJ455pw55tg5621Tu5UtAEYGbrAgcBtxfyvCKSW4018k318DNt2LNt0NO2eidXCv0J4O/AH4GljR1gZmea2TgzGzdjxozCRSYiGVleI99UDz/Thj3bBr3P9p3od/g2dGrXCgM6tWtFv8O30fh/Ewo2CWxmBwPT3X28mfVs7Dh37w/0B+jevbsXKDwRydDyGvmmevjn99rsJ5O70HDDHmc5ZppW7+RKIVcB9QB6m9mBQEtgFTMb5O4nFDAGEWmm5TXyHdu1YmoDz9f28LNp2NWg51/BEoC7XwRcBBB9AviDGn+R0rO8Rj6THr4a9uKhUhAikpXljc9rLL60JHIjmLs/BzyXxLlFJHPLu2mrscfVwy8duhNYRBrU1N24auRLnxKAiADL9vbnL1qi8gplTglARBrs7TdG5RXKhyaBRaTBtf2NUXmF8qEEICIZ9+pVXqG8KAGISKO9+natKrWks4xpDkBEGr2B6/LeW6nBL2NKACKSuq0QJVACEBFAa/vTSHMAIiIppQQgIpJSSgAiIimlOQCRFFpekTdJDyUAkZRpqsibpIeGgERSpql9eyU9lABEUqapfXslPTQEJFKG6o/x77l5B559bwbTZlWzghk17su8RkXe0kcJQKRM1Db6U2dVY0BtEz91VjWDXv38h+MaavxV5C2dlABEykD9id1lm/hlVZix1F2rgFJMCUCkDGRTz7/WUnc+ufagPEUkpUAJQKTE1B3fX7VVJWbw3fzFWf8ejfmLEoBICak/1DOrOvuGHzTmL4ESgEiRam5Pv3YiuFO9VUAa85daSgAiRai5Pf1OauQlA0oAIkVm5ISp/H74Gw0u18xEp3ateOnCvXIclZSjrO8ENrPjzWzTfAQjkna1Pf+4jb/G9iUbcT4BzABuNbMWwEzgA3e/MLdhiaRL3Zu4stEumhuYNX+xxvYla1knAHd/ysx+7u6XmFlboF8e4hJJjfrj/ZloVVlBv8O3UWMvzRK3GNwqZtYNWAi0zmE8IqlSO97fVONvwGorV2KEMX41/pILcSeBzwN+CfwGGJ27cETSI9PxfvX2JV9iJQB3Xwz8I8exiKRKJuUbtJxT8ilWAjCzfwGt3f1kM9vP3Z/KcVwiZSnTyV71+qUQ4g4BLQK+jr7fC2gyAZhZS+AFYKXovCPc/bKY5xcpCY2VaF6eCjM1/lIQcRPAfGBVM6sEOmf4moXAXu4+L3rdi2b2hLu/GjMGkaIWp0Szev5SSHFXAfUDPgL+BQzJ5AUezIt+rIy+4t3tIlIC4pRoVuOfchMmwE03Fex0cRPAzcAS4Gx3fzTTF5lZhZlNBKYDY9x9bAPHnGlm48xs3IwZM2KGJ5KckROmst0VT2V9U1endq3U+KfR/Plw112w887QrRtceinMnFmQU8dKAO5+GvAuMMDMemfxuhp33w5YF9jJzLZu4Jj+7t7d3bt36NAhTngiiaht+M8dNjHr4m0q4ZBC77wDv/sddOwIp54Kc+fCzTfDlCnQvn1BQoi7CmgPYGugGrgKGJXN6919lpk9B+wPvBUnBpFiEudu3hUMlrqWeqbKwoXw4INQVQUvvACVlXDkkXD22fCzn4FZQcOJOwm8O2FFz2B3n53JC8ysA7A4avxbAfsA18U8v0hRGDlhKpePejurHr8a/BT66CPo3x/uvDMM72y4IVx3HZx8Mqy5ZmJhxU0AS9z9xSxfsw5wj5lVEIaehmczfyBSbEZOmMr597/B4qWZrWVQmeaUWbwYHnkk9PbHjIGKCujdO/T299kHVog7BZs7cRPAoWY2nzCR+34mL3D3N4HtY55PpGjE6fVXVpjG+NPiiy9gwAC4/Xb48ktYd1248ko47bQw3l9E4iaAwwmN+WFmtrG7n57DmESKUpyGH0IRt8sO2UpDPuWspgZGjw69/cceA3c44AC47bbw3xWLc++tuFH9mR9LQfTKZUAixSjbSd4KM/52dFc1+uXuq6/CuH7//vDZZ2E8/8IL4YwzoEuXpKNrUi5KQeyJKoJKmcvmpq7KCuOGI9X4ly13ePbZ0Nt/6CFYsgT22gtuuAEOPRRatEg6wowVshSESMmaluFNXRruKWPffAP33BOGdT74AFZbDc45B848EzYrzfmduAngWWBToAoYnLtwRIpTu5Ur+W5+42P/6vWXKXd4+eXQ6A8fHtbx77ZbuFv3yCOhVaukI2yWuAlgP8J+AL2ir8dyFpFIkRk5YSrzFixp9Hn1+svQ7NkwaFAY5nnrLWjbFk4/Hc46C7bZJunociZuAmgHXAD8ETgtd+GIFIdM6va3a1XJxMv2K2BUknfjx4dGf8iQUKOnW7ewpLNvX2jTJunoci5uArgS2Nzd3zezpbkMSCRJYbXPm1Qvbvqf9ewsl4NKkfr+e7jvvtDwjxsXhnWOOy7csNW9e9LR5VXcLSGnmNmc6PsLcxuSSDIuHTmJQa9+nvHxHduV9vhv6r31Vhjbv/demDMHttoKbrkFTjgB2rVLOrqCiFsMbigwO6rpc4e7v5DbsETyL+6NXbV0Z28JWrAAHngg9PZffDEs2TzqqNDb79Gj4MXYkhZ3COhdd78SwMxuJRSGEyl6zW30a7VrValJ31Ly4YfhZq277grLOTfeOKzbP/nkgpVeLkZZJwAzGwBsbmazgDeBWTmPSiQPsh3iaUxlhXF5761yEJHk1eLF8PDDobf/n/+EcgyHHhp6+3vtVRTF2JKWdQJw9zOiDd63B3YC1sh5VCI5lM3EblO05LMEfPZZWLlzxx2hVEPnznDVVWHTlXXWSTq6ohJ3COhiYAvge8L+wCJFqbm9/tYtKrj6MO3TW/RqauCJJ0Jv//HHw2MHHRR6+/vvH0oxyzLiJoDV3f0oM2sB3AT8OocxiTRbLnr9J+zSmav6lM9NP2Xpyy9DT79//1CGee214ZJLwk1b66+fdHRFL24CWGhm3YBJQOscxiPSbMcPeIWXPvo29us1zFPkli6FZ54Jvf2RI0Pvf5994KabwoYrlZVJR1gy4iaAS4DfRl/DcheOSDzNWd2jBr9EzJgBd98devuTJ8Maa8D//V8oxrbJJklHV5Li3gi2ALghx7GIxBK3x68hnhLgHtbrV1XBiBGwaBHsvjtcfjkccQS0bJl0hCUt4wRgZgcAlxPqAL0B3Ojur+YpLpGMxGn8NbFbAmbNgoEDQ8P/zjuwyiqhENtZZ4U7diUnsvkEcCtwAvAOsAPwVzP7l7sPzUtkIssRd5JXvf4i5h5q8VRVwdChUF0NO+4YJnmPOQZaa7ox17JJAF+7+0vR90+b2SvAWEAJQAoqztLOlVZcgeuO2Fa9/mI0b15o8Kuq4PXXYeWVQz2es86CHXZIOrqylk0C+NTMrgKudPdFwGJgbn7CEvmp5kzy9thodQafsWseopJmefPNUIxt4ECYOzfU2f/Xv+D442HVVZOOLhWySQAOHA6cYWYfEraCHGxmm7j7h3mJToT4k7xa3VOEqqvDZG5VVdhpa6WV4Oijww1bu+6aumJsScs4Abj7sQBRGYitga7AtsDtZrahu6+XnxAlrUZOmMr5908k23u51OMvQu+/H3r7d98N330Hm24Kf/sbnHRSWM4piYizDPQw4FCgBlgBqNJEsORKc6t1apK3iCxaFG7UqqqCZ58NxdgOPzz09nv2VG+/CMRJAD3dvW/tD2b2LzQRLDnQ3Lo9avyLxCef/FiMbfr0UJLhmmvglFNCqQYpGnESwEpmdhDwBbAeoG2RpFma2/BrhU8RWLIkFGGrqoInnwy9+4MPDr39/fZTMbYiFScB/IowGbwNIQn8JqcRSSrkoja/Acer15+sqVNDT3/AAJgyJZRb/tOfQjG29TQtWOzi7AcwHxiUh1gkJfa98Tk+nP597Nerx5+wpUvh6adDb3/UqFCMbb/94B//CL1+FWMrGXGLwf3AzC5w9+tyEYyUv+Y2/hrnT9D06WFLxf794eOPw1aKv/99KMa20UZJRycxxNkScnjdH4HtACUAaVJzGn/V70mIO7zwQujtP/BA2Gbx5z+Hq6+Gww4L6/ilZMX5BDDH3U+v/cHM/p3DeKQMjZwwlXOHTYz1WvX4E/Ldd3DvvaHhf+89aNcOfvWrUJ5hiy2Sjk5yJE4CuLrez5dk8iIzWw+4F1gbWAr0d/ebY5xfSkjcu3jV40+AO/zvf6HRv+8+WLAAdt45DPscfXSo0SNlJZty0OfV+b7229nAeCCT/8OXAL9399fNrC0w3szGuPs7WcQrJSJOr1+rehIydy4MGRIa/okTQ9XNk04Kvf3tt086OsmjbD4BdI++Hol+Pgh4DTjbzO539+uX92J3/xL4Mvp+rpm9C3QilJeWMpLtEs9N1mzNmPN65i8gadjEiaE8w6BBoSJn167w73/DcceF+vtS9rJJAGsA3dx9HoCZXQaMAPYgfApYbgKoy8y6ANsTyknXf+5M4EyAzp07ZxGeFINsh3zU+BfY/PkwfHjo7Y8dG3bUOuaYcMPWzjurPEPKZJMAOgOL6vy8GFjf3avNbGGmv8TM2gAPAOe6+5z6z7t7f6A/QPfu3T2L+CRBcYZ81mrbQo1/obz7bujt33NP2G1r883DJuq/+AWsvnrS0UlCskkAQ4BXzexhwnDtwcBQM2tNhsM4ZlZJaPwHu/uD2QYrxSfuHb3q+RfAwoXw0EOht//88+EGrSOOCL39PfZQb1+yKgf9FzN7HNidkADOdvdx0dPHN/V6CzPHdwDvuvuNcYKV4hF3aecKBjcevZ1W9+TTxx+Hm7XuvBNmzIANNoBrrw3F2NZcM+nopIhkuwx0CWEJpxOGgLLRAzgRmGRmtS3Hxe7+eJa/RxIW94Yu9frzaMkSePTR0NsfPToUXzvkkNDb33dfWGGFpCOUIpTNMtDfAWcQhnAMGGRm/d39lkxe7+4vRq+TEtWcAm7apCVPpkyB228PxdimTYNOneDyy+G002DddZOOTopcNp8ATgN2dvfvAczsOuAVIKMEIKUt7g1doLt5c66mBp56KvT2H3003MDVqxfceiscdFDYeEUkA9n8SzHCLmC1alCPPhUuHTkpVuO/VtsWjL1k3zxElFJffx3G9fv3h08/DeP5F1wAZ5wRxvlFspRNArgLGGtmD0U/9yFM6kqZGjlhKucNm0iWW/IC6vXnjDs891zo7T/0UCjGtueecN110KcPtGiRdIRSwrJZBXSjmT1PmMw14BR3n5C3yCQxWuFTBL79NqzZr6qCDz6A1VaD3/wmlF7efPOko5MykdVgobuPJ9z1K2Uq7li/evw54A6vvhoa/WHDwjr+XXcNieCoo6CVdl+V3GoyAZjZXMKyz2WeAtzdVTSkDOiGrgTNmQODB4eG/803oW1bOPXUUIyta9eko5My1mQCcPe2hQhEkrPtZU8yZ2FN0wfWo8a/mV5/PTT6Q4bA99+Hypu33RaKsbVpk3R0kgJxdgQ7xN0fafpIKXbNWdevOj4xff99GN6pqoLXXgvDOsceG3r7O+6o8gxSUHE3hFECKHFxe/2g8f5Y3n479O7vvRdmz4YttwybqJ94YthtSyQBcRKAuiglbuOLHmNJjDqrups3SwsXhn10q6rgv/8NSzaPPDKUZ9h9d/X2JXFxEoBKNJeg5uzLq7H+LE2eHG7WuusumDkTNtoIrr8eTj4ZOnRIOjqRH+ie8RRQGYcCWLwYRo0KwzxjxoRibIceGnr7e++tYmxSlJQAytzOV4/h67mLmj6wHjX8Gfr881CI7fbb4auvQgG2K68Mxdg6dkw6OpHlipMAvs55FJIXcRr/VVaq4M0r9s9TRGWipgaefDKM7T/+eLiB64ADQm//gANUjE1KRtb/Ut1d1b2KXNzlner1N+Grr+COO8L4/uefw1prwUUXwemnQ5cuSUcnkjV1VcpIcyZ6/36Mavg0aOlSePbZ0NsfOTJsvLL33vC3v4Ux/srKpCMUiU0JoAw0p+FfcQXjr0d1VeNf38yZoQbPbbfBhx+GjdN/97tQjG3TTZOOTiQnlABKXNztGTXW3wB3ePnl0Nu///6wjr9HD/jzn8P6/ZYtk45QJKeUAEpU3NU9oI1aljF7NgwcGBr+t9+GVVYJ4/pnnQXbaE5EypcSQIlpTv0eUOP/E+PGhUZ/6FCYPx922CEs5+zbF1q3Tjo6kbxTAighm1/yOAtq4t+IrVIOhGJsQ4eGhn/8eFh55VB986yzoHv3pKMTKSglgBLQ3IZfpRyASZPChO7AgaH+/lZbwT//CSecAKuumnR0IolQAihizR3uSX3Dv2ABjBgRevsvvQQrrRR21jr7bNhtNxVjk9RTAihScVf3AKxoMLnfQTmOqIR88MGPxdi+/RY22QT++lc46SRo3z7p6ESKhhJAkWnOmn5I8Tj/okXw8MNhmOc//wnlGPr0Cb39PfdUMTaRBigBFJHmVO1M7Z28n34airHdcQd8/TWsvz5cfXXYU3fttZOOTqSoKQEUibgTvam8oaumJhRhq6qCJ54IY/kHHRR6+716hVLMItIkJYAisMGFj8XaZSd1vf5p00JPf8AA+OILWGcduPTScNNW585JRydScpQAEhR3lc+n16Zognfp0jCmX1UVxvhramDffeHvf4dDDlExNpFmUAJISNx9eVPT+M+YAXffHSZ1P/oI1lgDzjsvFGPbeOOkoxMpC0oACehy4WNZvyYVY/3u8OKLobc/YkRY2bPHHmGHrSOOCOv4RSRnCpYAzOxO4GBgurtvXajzFpttL3sy69eUfa9/1iy4997Q8L/7brgz9+yzQ3mGLbdMOjqRslXITwB3A/8E7i3gOYvOnIU1GR9rwCfl2vi7w2uvhUb/vvuguhp22gnuvBOOOSbU6BGRvCpYAnD3F8ysS6HOV4yyGfppWWG8d/WBeYwmIfPmwZAhoeGfMCFU3TzxxNDb79Yt6ehEUkVzAAWQ7d29ZVnD5403woTuoEEwd26os3/rrXD88aH+vogUXNElADM7EzgToHMZrO3OZuOWshvyqa4OO2tVVcErr4RJ3GOOCeP7u+yiYmwiCSu6BODu/YH+AN27d49fA7kIZLPUs6xW+bz3Xujt33MPfPdd2EP3xhtDMbbVV086OhGJFF0CKBfZ3t1b8o3/okXw0EOht//cc6EY2+GHh95+z57q7YsUoUIuAx0K9ATam9kU4DJ3v6NQ5y+UONU8/37MdnmKpgA++SSUXr7zTpg+Hbp0gX794JRTYK21ko5ORJajkKuAji3UuZIS5wavHhutXnr1fJYsgcceC7390aND7/6QQ0Jvf7/9VHpZpERoCCgH4tb0OWGXzlzVZ5s8RJQnU6eGTdMHDAjfd+wIf/5zKMa27rpJRyciWVICaKa4ZZxL5u7epUvhqafCpO4jj4RibL16hf10Dz44jPWLSEnS/73NEGfIB0qk8Z8+PYzr9+8fxvk7dIA//CEUY9tww6SjE5EcUAKIqSwbf3d4/vkwtv/gg7B4cVjBc801cNhhKsYmUmaUAGKI0/iv1bYFYy/ZNw/R5MC33/5YjO3996FdO/j1r0Nvf4stko5ORPJECSBLcRr/ouz1u8PYsaHRHzYMFiwId+fefTccdZSKsYmkgBJAFrJt/Ityy8a5c2Hw4NDwv/EGtGkDJ58cirFtV8L3I4hI1pQAMpRt4190vf4JE8JKnsGDQ0XOrl1DEjjuOGjbNunoRCQBSgAZKNnGf/58GD48NPRjx0LLltC3b7hha6edVJ5BJOWUAJqQTeNfNNU833kn9PbvvTfstrXFFnDzzaHu/mqrJR2diBQJJYDlyLbnn2jjv3BhWLpZVQUvvACVlXDkkaG3/7OfqbcvIstQAmhEyQz7fPTRj8XYZs4MN2ldd12Y2F1zzWRiEpGSoATQgGwa/0TW9y9eDI8+Gnr7Tz0FFRXQu3fo7e+zj4qxiUhGlADqyXbf3oI2/l98EYqx3X47TJsWCrBdcQWcdhp0KrLlpiJS9JQUPHxnAAALV0lEQVQA6sh2wrcgm7bX1ISSy7fdFnr97rD//vDvf8OBB6oYm4jEptYjkk3jv6LB5H55HvP/6qsfi7F99lkYz7/gAjjjDNhgg/yeW0RSQQmAIhrzd4dnnw1j+w89FDZe2WsvuP566NMHWrTIz3lFJJVSnwCyafzztoHLN9+EDdRvuw0++CCs1T/nnFCMbbPNcn8+ERFSngCy7fnntPF3h1deCb394cPDOv7ddoNLLw3r91u1yt25REQakMoEcPyAV3jpo28zPn5FI3fDPrNn/1iMbdKkUIfntNNCMbZtt83NOUREMpC6BJDtFo6brNmaMef1bP6Jx48Pjf6QIaFGz/bbhwneY48NFTlFRAosVQmg4I3/99/DffeFhn/cuDCsc+yx4Yat7t1VnkFEEpWaBHD8gFeyavxXWakifuP/1ls/FmObMwe23BJuuQVOOCHstiUiUgRSkwCyGfOP1fNfsAAeeCD09l98MSzZPOqo0Nvv0UO9fREpOqlIANms9umx0eoMPmPXzH/5hx+Gsfy77grLOTfeGG64IRRja98++2BFRAqk7BPABlk0/hlv4bh4MYwaFXr7Tz8dirH16RN6+3vtpWJsIlISyjoBjJwwlUxH/TMq5/z55zBgQCjG9tVXsN568Je/wKmnQseOzYpVRKTQyjoBnDtsYkbHLbfxr6mBJ54Ik7qPPx5u4DrwwNDbP+CA0PsXESlBZZsANr/k8YyO67HR6g0/8eWXcMcdocf/+eew9tpw8cVw+umw/vo5jFREJBllmwAyWfK5ovHTCd+lS+GZZ8LY/sMPh2Js++wDN94YNlyprMxjxCIihVWWCSDTVT8/lHSeORPuvjsM80yeDGusAeeeG4qxbbJJ/gIVEUlQWSaATHza78CwXr+qCu6/HxYtgt13h8svhyOOgJYtkw5RRCSvCpoAzGx/4GagArjd3a/N9Tl2vnrMcp9fZcE8rp7/BmxzAbz9NqyySujpn3UWbL11rsMRESlaBUsAZlYB/AvYF5gCvGZmo9z9nVye5+u5ixp8fNsvP+D4CU/Q+90XaLVkIey4Y1jO2bcvtG6dyxBEREpCIT8B7ARMdvePAczsPuBQIKcJoK6VF1XT+53nOX7iE2zz9UfMr1yJt/Y8mB37XQQ77JCv04qIlIRCJoBOwBd1fp4C7JyPE602fzb/9+IQDnv7Gdouqua99utz6b6/5OGtejLppqPzcUoRkZJTyATQUDW0ZdZqmtmZwJkAnTt3jnWi6sqV2P+Dl3lqk10YvN2BvN5pczDj78dsF+v3iYiUo0ImgCnAenV+XheYVv8gd+8P9Afo3r175vWb61hQ2ZLdz76TRSv+dN1+RnV+RERSopBVy14DNjGzDcysBdAXGJXrk2yyZpjQrd/41z4uIiJBwRKAuy8BfgOMBt4Fhrv727k+z5jzei7T2OdsW0cRkTJS0PsA3P1xILMiPc2gxl5EpGkqXC8iklJKACIiKaUEICKSUkoAIiIppQQgIpJS5h7rXquCMLMZwGfN+BXtgZk5CidJ5XIdUD7XousoPuVyLc29jvXdvUMmBxZ1AmguMxvn7t2TjqO5yuU6oHyuRddRfMrlWgp5HRoCEhFJKSUAEZGUKvcE0D/pAHKkXK4DyudadB3Fp1yupWDXUdZzACIi0rhy/wQgIiKNUAIQEUmpkk8AZra/mb1vZpPN7MIGnl/JzIZFz481sy6FjzIzGVzLyWY2w8wmRl+nJxFnU8zsTjObbmZvNfK8mdk/out808y6FTrGTGRwHT3NbHad9+PPhY4xE2a2npk9a2bvmtnbZva7Bo4plfckk2sp+vfFzFqa2f/M7I3oOq5o4Jj8t13uXrJfQAXwEbAh0AJ4A9iy3jG/Aqqi7/sCw5KOuxnXcjLwz6RjzeBa9gC6AW818vyBwBOEbUJ3AcYmHXPM6+gJPJp0nBlcxzpAt+j7tsAHDfzbKpX3JJNrKfr3Jfo7t4m+rwTGArvUOybvbVepfwLYCZjs7h+7+yLgPuDQesccCtwTfT8C2NvMGtqfOGmZXEtJcPcXgG+Xc8ihwL0evAq0M7N1ChNd5jK4jpLg7l+6++vR93MJGzLV3x+1VN6TTK6l6EV/53nRj5XRV/0VOXlvu0o9AXQCvqjz8xSW/cfwwzEediWbDaxRkOiyk8m1ABwRfUQfYWbrNfB8Kcj0WkvBrtHH+CfMbKukg2lKNIywPaHHWVfJvSfLuRYogffFzCrMbCIwHRjj7o2+J/lqu0o9ATSUDetn0UyOKQaZxPkI0MXdtwWe5sfeQakplfekKa8T6q50BW4BRiYcz3KZWRvgAeBcd59T/+kGXlK070kT11IS74u717j7dsC6wE5mtnW9Q/L+npR6ApgC1O0FrwtMa+wYM1sRWJXi/Fjf5LW4+zfuvjD6cQCwQ4Fiy7VM3rei5+5zaj/Ge9jutNLM2iccVoPMrJLQYA529wcbOKRk3pOmrqWU3hcAd58FPAfsX++pvLddpZ4AXgM2MbMNzKwFYaJkVL1jRgEnRd8fCTzj0axKkWnyWuqNyfYmjH+WolHAL6KVJ7sAs939y6SDypaZrV07JmtmOxH+f/om2aiWFcV4B/Cuu9/YyGEl8Z5kci2l8L6YWQczaxd93wrYB3iv3mF5b7sKuil8rrn7EjP7DTCasIrmTnd/28yuBMa5+yjCP5aBZjaZkD37Jhdx4zK8lnPMrDewhHAtJycW8HKY2VDCSoz2ZjYFuIwwyYW7VwGPE1adTAbmA6ckE+nyZXAdRwK/NLMlQDXQt0g7Fz2AE4FJ0ZgzwMVAZyit94TMrqUU3pd1gHvMrIKQoIa7+6OFbrtUCkJEJKVKfQhIRERiUgIQEUkpJQARkZRSAhARSSklABGRlFICEBFJKSUAKRtmNq/po0SklhKAiEhKKQFIWTKz88zsrejr3DqP/8nM3jOzMWY21Mz+kMXvvC/aoGOsmX1mZgflJ3qRwijpUhAiDTGzHQilDHYmVFQca2bPE0psHEEoIbwioWrk+Cx+dVdgpLsfY2a7AzcCj+Uy9rrMbDV3/y5fv19ECUBKhpk9DazdwFOXuPvDdX7eHXjI3b+PXvcg8DPCJ96H3b06evyR6L8bApcAq7r7kdFjrYFbgUWESo0PAu2B2q373gFWy+X1NeAmmqj3ZGa9gDXdfWCeY5EypAQgJcPd98nw0MZ2TWrwcXf/GDjNzEbUefhwYIS7P2JmwwhbD37o7gui57sRtu2sLU98BbAyYTvPy4DrgT8BVwIXEbZbHA1sRihWVlHvNedEr3Hgs+h8m0dDVDfXO/ZcwqePOYRPOb/O8O8i8hOaA5By9ALQx8xWjnryhwH/BV4EDok25G4DLG8Mf11+3CGrhjD80zl6bWtCg3xT9PyZQCtgFmGf1xnA58DfCA17d2Cou19E2P1pjfqvAX5J+HTye3f/BzATGOTuf23k2Hvc/WJgJeD9+H8qSTN9ApCy4+6vm9ndwP+ih2539wkAZjaK0HP/DBhH2GavIVMISWAioaPUFRhMGA5aBbjG3V+Kjt0e+HXtZj1RctkQWOLu88xsR+Dl6NhV3X2GmdV/zV3Av+ucf9sozoZ+/93AbWbWFphZhKWOpUSoHLSkipm1iRrllQmfFM4kJIOrgX0JyaJf1Mv/J7CA8MnhLOAMd1+mt21mhwDHET4xvEAYProCOJqw0c9ZhB79ioRhpf/Ue80zhJ58b0Ld936EeYwjgGuBjesduxrQi1C3f567/zGXfyNJDyUASRUzGwJsCbQkDKP0y/B1U4H13H1pjHMOdPcTs32dSL4pAYiIpJQmgUVEUkoJQEQkpZQARERSSglARCSllABERFJKCUBEJKWUAEREUkoJQEQkpZQARERS6v8BHLo2f5WV34EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sorted_pvals2 = sorted(-np.log10(pvals2))\n",
    "\n",
    "plt.scatter(x = sorted(-np.log10(np.arange(1,pvals.shape[0]+1)/pvals.shape[0])), y = sorted_pvals2)\n",
    "\n",
    "plt.plot(sorted(-np.log10(np.arange(1,pvals2.shape[0]+1)/pvals2.shape[0])), \n",
    "         sorted(-np.log10(np.arange(1,pvals2.shape[0]+1)/pvals2.shape[0])), 'r-')\n",
    "\n",
    "plt.title('QQ-plot of p-values')\n",
    "plt.xlabel(r'$-\\log_{10}p_{expected}$')\n",
    "plt.ylabel(r'$-\\log_{10}p_{observed}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple Testing Correction\n",
    "\n",
    "Let's perform 1000 independent samples t-tests to test for differences in the means of two sampled groups, A and B, each of size $n=10$. \n",
    "\n",
    "Among these 1000 tests, let's predetermine that 100 of them will yield **true positives** - that is, the means of samples A and B come from two different normal distributions. For the rest of the tests, we will sample from the same normal distributions - these are **true negatives**.\n",
    "\n",
    "##### Let's start by generating the samples for each group:\n",
    "- for the **true positive** tests, we will sample 10 values from the standard normal distribution for group A, and 10 values from $N(\\mu=1,\\sigma=1)$ for group B\n",
    "- for the **true negative** tests, we will sample 10 values from the standard normal distribution for group A and for group B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ntests = 1000\n",
    "TPtests = 100\n",
    "TNtests = Ntests - TPtests\n",
    "sampleSize = 10\n",
    "\n",
    "# True positive cases\n",
    "tp_samplesA = stats.norm.rvs(size = (TPtests, sampleSize))\n",
    "tp_samplesB = stats.norm.rvs(loc = 1, scale = 0, size = (TPtests, sampleSize))\n",
    "\n",
    "# True negative cases\n",
    "tn_samplesA = stats.norm.rvs(size = (TNtests, sampleSize))\n",
    "tn_samplesB = stats.norm.rvs(size = (TNtests, sampleSize))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Independent samples t-test\n",
    "We will perform a [**two-sided, independent samples t-test**](https://en.wikipedia.org/wiki/Student%27s_t-test#Independent_two-sample_t-test) to determine whether the two samples being compared have the same means:\n",
    "\n",
    "$$H_0: \\mu_A = \\mu_B$$\n",
    "$$H_A: \\mu_A \\neq \\mu_B$$\n",
    "\n",
    "This test assumes that the populations have identical variances (which is why we kept $\\sigma=1$).\n",
    "\n",
    "We will use the scipy implementation for the t-test: [`scipy.stats.ttest_ind`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# concatenate the samples for groups A and B into the same arrays\n",
    "all_samplesA = np.concatenate((tp_samplesA, tn_samplesA), axis = 0)\n",
    "all_samplesB = np.concatenate((tp_samplesB, tn_samplesB), axis = 0)\n",
    "\n",
    "# perform the t-test (row-wise) \n",
    "res = stats.ttest_ind(all_samplesA, all_samplesB, axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function returns the test statistic, $t$, and corresponding $p$ for each test performed. Let's take a look at the p-value distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEGlJREFUeJzt3X+MZWV9x/H3R1ZtrLSgOxgCbAfMYorGrjohNlaLxVrABrRRy6YqWuKKlaY/TFPUpBgbE6pSElMLXcIGbBVBKbqpWCXUSjViHQS3i0BdcIWVDTuCBRss7cK3f8zZdlxmd+7ec+8M99n3K7mZc57znHO+z87sJ2eee+6ZVBWSpHY9ZaULkCSNl0EvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatyqlS4AYPXq1TU9Pb3SZUjSRLn55pt/WFVTS/V7UgT99PQ0s7OzK12GJE2UJN8fpJ9TN5LUOINekhpn0EtS4wx6SWrckkGfZFOSXUm2Lmi7Ksmt3Wt7klu79ukkP1mw7ZJxFi9JWtogd91cDvwV8PE9DVX123uWk1wIPLSg/11VtW5UBUqS+lky6KvqxiTTi21LEuCNwK+NtixJ0qj0naN/OXB/VX13QduxSW5J8pUkL+95fElST30/MLUeuHLB+k5gTVU9kOQlwGeTPL+qHt57xyQbgA0Aa9as6VmGJGlfhg76JKuA3wJesqetqh4FHu2Wb05yF3A88ISPvVbVRmAjwMzMTK+/UD593uf77D607Re8ZkXOK0kHos/UzauAO6pqx56GJFNJDumWjwPWAnf3K1GS1Mcgt1deCXwdeF6SHUnO7jadyU9P2wC8AtiS5NvAZ4BzqurBURYsSTowg9x1s34f7W9dpO0a4Jr+ZUmSRsVPxkpS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1bsmgT7Ipya4kWxe0vT/JD5Lc2r1OW7DtPUm2JbkzyW+Mq3BJ0mAGuaK/HDhlkfaLqmpd97oOIMkJwJnA87t9/jrJIaMqVpJ04JYM+qq6EXhwwOOdAXyqqh6tqu8B24ATe9QnSeqpzxz9uUm2dFM7h3dtRwH3Luizo2t7giQbkswmmZ2bm+tRhiRpf4YN+ouB5wLrgJ3AhV17Fulbix2gqjZW1UxVzUxNTQ1ZhiRpKUMFfVXdX1WPVdXjwKX8//TMDuCYBV2PBu7rV6IkqY+hgj7JkQtWXwfsuSNnM3BmkqcnORZYC/xrvxIlSX2sWqpDkiuBk4DVSXYA5wMnJVnH/LTMduAdAFV1W5Krge8Au4F3VdVj4yldkjSIJYO+qtYv0nzZfvp/EPhgn6IkSaPjJ2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4JYM+yaYku5JsXdD24SR3JNmS5Nokh3Xt00l+kuTW7nXJOIuXJC1tkCv6y4FT9mq7HnhBVb0Q+HfgPQu23VVV67rXOaMpU5I0rCWDvqpuBB7cq+1LVbW7W70JOHoMtUmSRmAUc/S/C3xhwfqxSW5J8pUkLx/B8SVJPazqs3OS9wG7gU90TTuBNVX1QJKXAJ9N8vyqeniRfTcAGwDWrFnTpwxJ0n4MfUWf5CzgN4HfqaoCqKpHq+qBbvlm4C7g+MX2r6qNVTVTVTNTU1PDliFJWsJQQZ/kFOBPgdOr6pEF7VNJDumWjwPWAnePolBJ0nCWnLpJciVwErA6yQ7gfObvsnk6cH0SgJu6O2xeAXwgyW7gMeCcqnpw0QNLkpbFkkFfVesXab5sH32vAa7pW5QkaXT8ZKwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcb0eaiYdDKbP+/yKnHf7Ba9ZkfOqPV7RS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxAwV9kk1JdiXZuqDtWUmuT/Ld7uvhXXuSfDTJtiRbkrx4XMVLkpY26BX95cApe7WdB9xQVWuBG7p1gFOBtd1rA3Bx/zIlScMaKOir6kbgwb2azwCu6JavAF67oP3jNe8m4LAkR46iWEnSgeszR/+cqtoJ0H09oms/Crh3Qb8dXZskaQWM483YLNJWT+iUbEgym2R2bm5uDGVIkqBf0N+/Z0qm+7qra98BHLOg39HAfXvvXFUbq2qmqmampqZ6lCFJ2p8+Qb8ZOKtbPgv43IL2t3R337wUeGjPFI8kafkN9Ddjk1wJnASsTrIDOB+4ALg6ydnAPcAbuu7XAacB24BHgLeNuGZJ0gEYKOirav0+Np28SN8C3tWnKEnS6AwU9JLUsunzPr9i595+wWvGfg4fgSBJjTPoJalxTt1MoNZ/zVzMSo5ZmnRe0UtS4wx6SWqcQS9JjTPoJalxBr0kNc67biQ9wUrd5bRSd3W1zqDXAfE2R2nyGPQ9GHoaJ3++NCrO0UtS4wx6SWqcQS9JjXOOXtKThu9LjIdX9JLUOINekhpn0EtS4wx6SWqcQS9JjRv6rpskzwOuWtB0HPBnwGHA24G5rv29VXXd0BVKknoZOuir6k5gHUCSQ4AfANcCbwMuqqqPjKRCSVIvo5q6ORm4q6q+P6LjSZJGZFRBfyZw5YL1c5NsSbIpyeEjOockaQi9gz7J04DTgU93TRcDz2V+WmcncOE+9tuQZDbJ7Nzc3GJdJEkjMIor+lOBb1XV/QBVdX9VPVZVjwOXAicutlNVbayqmaqamZqaGkEZkqTFjCLo17Ng2ibJkQu2vQ7YOoJzSJKG1OuhZkmeAfw68I4FzR9Ksg4oYPte2yRJy6xX0FfVI8Cz92p7c6+KJEkj5SdjJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhq3qu8BkmwHfgw8BuyuqpkkzwKuAqaB7cAbq+pHfc8lSTpwo7qif2VVrauqmW79POCGqloL3NCtS5JWwLimbs4AruiWrwBeO6bzSJKWMIqgL+BLSW5OsqFre05V7QTovh6x905JNiSZTTI7Nzc3gjIkSYvpPUcPvKyq7ktyBHB9kjsG2amqNgIbAWZmZmoEdUiSFtH7ir6q7uu+7gKuBU4E7k9yJED3dVff80iShtMr6JP8bJJD9ywDrwa2ApuBs7puZwGf63MeSdLw+k7dPAe4NsmeY32yqv4xyTeBq5OcDdwDvKHneSRJQ+oV9FV1N/BLi7Q/AJzc59iSpNHwk7GS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWrc0EGf5JgkX05ye5LbkvxB1/7+JD9Icmv3Om105UqSDtSqHvvuBt5dVd9Kcihwc5Lru20XVdVH+pcnSepr6KCvqp3Azm75x0luB44aVWGSpNEYyRx9kmngRcA3uqZzk2xJsinJ4aM4hyRpOL2DPskzgWuAP6yqh4GLgecC65i/4r9wH/ttSDKbZHZubq5vGZKkfegV9EmeynzIf6Kq/h6gqu6vqseq6nHgUuDExfatqo1VNVNVM1NTU33KkCTtR5+7bgJcBtxeVX+5oP3IBd1eB2wdvjxJUl997rp5GfBm4N+S3Nq1vRdYn2QdUMB24B29KpQk9dLnrpuvAllk03XDlyNJGjU/GStJjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuLEFfZJTktyZZFuS88Z1HknS/o0l6JMcAnwMOBU4AVif5IRxnEuStH/juqI/EdhWVXdX1X8DnwLOGNO5JEn7Ma6gPwq4d8H6jq5NkrTMVo3puFmkrX6qQ7IB2NCt/meSO3ucbzXwwx77T5qDbbzgmA8WB92Y8xe9xvwLg3QaV9DvAI5ZsH40cN/CDlW1Edg4ipMlma2qmVEcaxIcbOMFx3ywcMzjMa6pm28Ca5Mcm+RpwJnA5jGdS5K0H2O5oq+q3UnOBb4IHAJsqqrbxnEuSdL+jWvqhqq6DrhuXMffy0imgCbIwTZecMwHC8c8BqmqpXtJkiaWj0CQpMZNTNAv9UiFJE9PclW3/RtJppe/ytEaYMx/nOQ7SbYkuSHJQLdaPZkN+uiMJK9PUkkm/g6NQcac5I3d9/q2JJ9c7hpHbYCf7TVJvpzklu7n+7SVqHNUkmxKsivJ1n1sT5KPdv8eW5K8eKQFVNWT/sX8G7p3AccBTwO+DZywV5/fAy7pls8ErlrpupdhzK8EntEtv/NgGHPX71DgRuAmYGal616G7/Na4Bbg8G79iJWuexnGvBF4Z7d8ArB9pevuOeZXAC8Gtu5j+2nAF5j/DNJLgW+M8vyTckU/yCMVzgCu6JY/A5ycZLEPbk2KJcdcVV+uqke61ZuY/7zCJBv00Rl/DnwI+K/lLG5MBhnz24GPVdWPAKpq1zLXOGqDjLmAn+uWf569PoczaarqRuDB/XQ5A/h4zbsJOCzJkaM6/6QE/SCPVPi/PlW1G3gIePayVDceB/oYibOZvyKYZEuOOcmLgGOq6h+Ws7AxGuT7fDxwfJKvJbkpySnLVt14DDLm9wNvSrKD+bv3fn95SlsxY31szNhurxyxJR+pMGCfSTLweJK8CZgBfnWsFY3ffsec5CnARcBbl6ugZTDI93kV89M3JzH/W9u/JHlBVf3HmGsbl0HGvB64vKouTPLLwN92Y358/OWtiLHm16Rc0S/5SIWFfZKsYv7Xvf39qvRkN8iYSfIq4H3A6VX16DLVNi5LjflQ4AXAPyfZzvxc5uYJf0N20J/tz1XV/1TV94A7mQ/+STXImM8Grgaoqq8DP8P8c3BaNdD/92FNStAP8kiFzcBZ3fLrgX+q7l2OCbXkmLtpjL9hPuQnfd4WlhhzVT1UVaurarqqppl/X+L0qppdmXJHYpCf7c8y/8Y7SVYzP5Vz97JWOVqDjPke4GSAJL/IfNDPLWuVy2sz8Jbu7puXAg9V1c5RHXwipm5qH49USPIBYLaqNgOXMf/r3Tbmr+TPXLmK+xtwzB8Gngl8unvf+Z6qOn3Fiu5pwDE3ZcAxfxF4dZLvAI8Bf1JVD6xc1f0MOOZ3A5cm+SPmpzDeOskXbkmuZH7qbXX3vsP5wFMBquoS5t+HOA3YBjwCvG2k55/gfztJ0gAmZepGkjQkg16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXFpFkOskdSa7ong/+mSTPWOm6pGEY9NK+PQ/YWFUvBB5m/m8eSBPHoJf27d6q+lq3/HfAr6xkMdKwDHpp3/Z+PojPC9FEMuilfVvTPQsd5p+P/tWVLEYalkEv7dvtwFlJtgDPAi5e4XqkoUzEY4qlFfJ4VZ2z0kVIfXlFL0mN83n0ktQ4r+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4/4XVGBIepgqFfsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(res.pvalue)\n",
    "plt.xlabel('p')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a significance cutoff of $\\alpha = 0.05$, how many significant p-values do we observe?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "132 significant results\n"
     ]
    }
   ],
   "source": [
    "nSignificant = sum(res.pvalue < 0.05)\n",
    "print('%d significant results' % nSignificant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many of these are true positives? How many are false positives?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "86 TPs, 46 FPs\n"
     ]
    }
   ],
   "source": [
    "# number of TPs\n",
    "tp = sum(res.pvalue[:TPtests] < 0.05)\n",
    "fp = nSignificant - tp\n",
    "print('%d TPs, %d FPs' % (tp, fp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize that in the p-value distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEKCAYAAAAcgp5RAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEpxJREFUeJzt3X+QXWV9x/H3FxLcgiiQLAxlSTfMBDDBdGB2EKtjU+O0iJEwDNZQxGiZZhCxVp0WrH9E2jqDpdTWGUGCAWNL+U1JRrSWQRiq09AGkQ1JoKS4JFuoiUFjAaMg3/5xT+iS7mZv7rl3b+6T92smc+8595z7fJ/8+OTZ557z3MhMJEnlOqjbBUiSOsugl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBVuWrcLAJg5c2YODg5OfcNPPNF4POmkqW9bkmp6+OGHf5SZ/ZMdt18E/eDgIOvWrZv6hhcsaDw+8MDUty1JNUXE080c59SNJBXOoJekwhn0klS4/WKOXpI65aWXXmJ0dJRdu3Z1u5SW9fX1MTAwwPTp01s636CXVLTR0VEOP/xwBgcHiYhul7PPMpMdO3YwOjrK7NmzW3oPp24kFW3Xrl3MmDGjJ0MeICKYMWNGrZ9IDHpJxevVkN+tbv0GvSQVzqCXpML1/Iexg5ff0/K5tzy1gzNOmNHGaiRpfNdddx2PPvoo11xzzav75s2bx5133snJJ5/c0bZ7PuglaV/UGRyOZ+TK9zR13PDwMKeeeuqr27t27WLLli3MmTOnrfWMx6kbSZoC69ev57TTTnvN9oknnsjBBx/c8bYd0UvSFNiwYQPnnnvuq1fQPP/88yxatGhK2nZEL0kdtnXrVvr7+3n66acZGRlhZGSEJUuWMH/+/FePef/738/VV1/dkfYNeknqsOHhYebNm/eafRs3buTNb34zAKtXr2bRokUMDw93pH2DXpI6bP369cydO/c1+zZs2MD8+fPZtWsXt99+OxdeeCE7d+7sSPsGvSR12J5B/9xzz5GZHHPMMVx11VU8//zzXHzxxWzYsIGf/exnbW/fD2MlHVCavRyynW666abXbB911FFs27aNLVu2MDIywt133w3AFVdcwfDwMG95y1va2r5BL0ldMmvWLFauXPnq9vLlyzvSjlM3klQ4g16SCmfQS1LhDHpJKpxBL0mFmzToI+KGiNgWEY+N2XdVRDweEcMR8Y8RccSY1z4dEZsj4omI+J1OFS5Jak4zI/qvAmfuse9e4JTMnA/8B/BpgIiYCywB5lXnXBMRnV+aTZI0oUmDPjMfBJ7bY98/Z+bL1eZaYKB6vhi4JTN/npk/ADYDp7exXknSPmrHHP3vA9+snh8HbB3z2mi1T5LUJbWCPiI+A7wM7L6/d7yvKs8Jzl0WEesiYt327dvrlCFJ+73rrruOSy655DX75s2bx+OPP97xtlteAiEilgKLgIWZuTvMR4Hjxxw2ADwz3vmZuQJYATA0NDTufwaS1HaffWOb36+5FSd77qsEI+JM4DLg7Mx8ccxLa4AlEfG6iJgNzAH+rX6ZktTb9uuvEoyIm4EFwMyIGAWW07jK5nXAvdXXYq3NzIszc0NE3AZspDGl89HM/GWnipekXtHNrxKcNOgz8/xxdq8cZ9/u4z8HfK5OUZJUkt1fJTh2Pv7SSy/lhBNO4IUXXuCSSy7hkEMOYcGCBVxwwQVtb987YyWpw/b2VYJ33XUX5513Htdffz1r1qzpSPsGvSR12N6+SnB0dJTjj29cw9Kp+XqDXpI6bG9fJTgwMMDo6CgAr7zySkfa9xumJB1Ymrwcsp0m+ipBgHPPPZdLL72Ue+65h/e+970dad+gl6QuOuyww7jxxhs72oZTN5JUOINekgpn0EtS4Qx6SSqcQS+peP+37mJvqlu/QS+paH19fezYsaNnwz4z2bFjB319fS2/h5dXSira7huSevl7L/r6+hgYGJj8wAkY9JKKNn36dGbPnt3tMrrKqRtJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSrcpEEfETdExLaIeGzMvqMi4t6IeLJ6PLLaHxHxxYjYHBHDEXFaJ4uXJE2umRH9V4Ez99h3OXBfZs4B7qu2Ad4NzKl+LQOubU+ZkqRWTRr0mfkg8NweuxcDq6rnq4Bzxuz/WjasBY6IiGPbVawkad+1Okd/TGY+C1A9Hl3tPw7YOua40WqfJKlL2v1hbIyzb9xFoCNiWUSsi4h1vbx8qCTt71oN+h/unpKpHrdV+0eB48ccNwA8M94bZOaKzBzKzKH+/v4Wy5AkTabVoF8DLK2eLwVWj9n/werqmzOAnbuneCRJ3THpF49ExM3AAmBmRIwCy4Ergdsi4iJgC/C+6vBvAGcBm4EXgQ93oGZJ0j6YNOgz8/wJXlo4zrEJfLRuUZKk9vHOWEkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKVyvoI+ITEbEhIh6LiJsjoi8iZkfEQxHxZETcGhGHtKtYSdK+aznoI+I44A+Bocw8BTgYWAJ8HvhCZs4Bfgxc1I5CJUmtqTt1Mw34lYiYBhwKPAu8E7ijen0VcE7NNiRJNbQc9Jn5X8BfAVtoBPxO4GHgJ5n5cnXYKHDceOdHxLKIWBcR67Zv395qGZKkSdSZujkSWAzMBn4VOAx49ziH5njnZ+aKzBzKzKH+/v5Wy5AkTaLO1M27gB9k5vbMfAm4C/gN4IhqKgdgAHimZo2SpBrqBP0W4IyIODQiAlgIbATuB86rjlkKrK5XoiSpjjpz9A/R+ND1e8D66r1WAJcBn4yIzcAMYGUb6pQktWja5IdMLDOXA8v32P0UcHqd95UktY93xkpS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMLVCvqIOCIi7oiIxyNiU0S8NSKOioh7I+LJ6vHIdhUrSdp3dUf0fwv8U2aeDPw6sAm4HLgvM+cA91XbkqQumdbqiRHxBuAdwIcAMvMXwC8iYjGwoDpsFfAAcFmdIqWe9dk31jh3Z/vq0AGtzoj+BGA7cGNEPBIRX4mIw4BjMvNZgOrx6DbUKUlqUZ2gnwacBlybmacCL7AP0zQRsSwi1kXEuu3bt9coQ5K0N3WCfhQYzcyHqu07aAT/DyPiWIDqcdt4J2fmiswcysyh/v7+GmVIkvam5aDPzP8GtkbESdWuhcBGYA2wtNq3FFhdq0JJUi0tfxhb+RhwU0QcAjwFfJjGfx63RcRFwBbgfTXbkCTVUCvoM/P7wNA4Ly2s876SpPbxzlhJKlzdqRt1weDl99Q6f+TK97SpEkm9wBG9JBXOEX0PGun7vZrv4B2X0oHEEb0kFc6gl6TCGfSSVDjn6KXJ1FmBUtoPOKKXpMI5oldPqHPvgPcNaFJ1f2rbz787wKBXT6h3Sen+/Y9Q6jSnbiSpcI7oa3A6QVIvcEQvSYUz6CWpcAa9JBXOoJekwvX8h7G1Lrs76AXg7d1p20v+1EFeKKCxej7oJf1/PTsIqXPjUjdvWtrP63bqRpIK54he+2Y/H7moDVzErTgGvaZOtwLE4NIBrvbUTUQcHBGPRMTXq+3ZEfFQRDwZEbdGxCH1y5QktaodI/qPA5uAN1Tbnwe+kJm3RMSXgYuAa9vQTkesfWoHS1q8QmGkr83FSFIH1BrRR8QA8B7gK9V2AO8E7qgOWQWcU6cNSVI9dadu/gb4E+CVansG8JPMfLnaHgWOq9mGJKmGlqduImIRsC0zH46IBbt3j3NoTnD+MmAZwKxZs1otQyqXHyKrTeqM6N8GnB0RI8AtNKZs/gY4IiJ2/wcyADwz3smZuSIzhzJzqL+/v0YZkqS9aTnoM/PTmTmQmYPAEuDbmXkBcD9wXnXYUmB17SolSS3rxHX0lwG3RMRfAI8AKzvQhiS9llNdE2pL0GfmA8AD1fOngNPb8b6SpPq8M7ZLaq0uWPf6fUc+0gHFoO+SeqsLSlLzXL1Skgp3wI/ozzhok6NrSUVzRC9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqnEEvSYUz6CWpcAa9JBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVruWgj4jjI+L+iNgUERsi4uPV/qMi4t6IeLJ6PLJ95UqS9lWdEf3LwKcy803AGcBHI2IucDlwX2bOAe6rtiVJXdJy0Gfms5n5ver5/wCbgOOAxcCq6rBVwDl1i5Qkta4tc/QRMQicCjwEHJOZz0LjPwPg6AnOWRYR6yJi3fbt29tRhiRpHLWDPiJeD9wJ/FFm/rTZ8zJzRWYOZeZQf39/3TIkSROoFfQRMZ1GyN+UmXdVu38YEcdWrx8LbKtXoiSpjjpX3QSwEtiUmX895qU1wNLq+VJgdevlSZLqmlbj3LcBFwLrI+L71b4/Ba4EbouIi4AtwPvqlShJqqPloM/M7wAxwcsLW31fSVJ7eWesJBXOoJekwhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvSQVzqCXpMIZ9JJUOINekgpn0EtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TCGfSSVDiDXpIKZ9BLUuEMekkqXMeCPiLOjIgnImJzRFzeqXYkSXvXkaCPiIOBLwHvBuYC50fE3E60JUnau06N6E8HNmfmU5n5C+AWYHGH2pIk7UWngv44YOuY7dFqnyRpik3r0PvGOPvyNQdELAOWVZvPR8QTLbY1E/hRi+c2XPHTWqd3Qf0+9x77fGA48Pp8RdTp8681c1Cngn4UOH7M9gDwzNgDMnMFsKJuQxGxLjOH6r5PL7HPBwb7fGCYij53aurm34E5ETE7Ig4BlgBrOtSWJGkvOjKiz8yXI+JS4FvAwcANmbmhE21JkvauU1M3ZOY3gG906v3HqD3904Ps84HBPh8YOt7nyMzJj5Ik9SyXQJCkwvVM0E+2pEJEvC4ibq1efygiBqe+yvZqos+fjIiNETEcEfdFRFOXWu3Pml06IyLOi4iMiJ6/QqOZPkfE71Z/1hsi4h+musZ2a+Lv9qyIuD8iHqn+fp/VjTrbJSJuiIhtEfHYBK9HRHyx+v0YjojT2lpAZu73v2h8oPufwAnAIcCjwNw9jrkE+HL1fAlwa7frnoI+/xZwaPX8IwdCn6vjDgceBNYCQ92uewr+nOcAjwBHVttHd7vuKejzCuAj1fO5wEi3667Z53cApwGPTfD6WcA3adyDdAbwUDvb75URfTNLKiwGVlXP7wAWRsR4N271ikn7nJn3Z+aL1eZaGvcr9LJml874c+AvgV1TWVyHNNPnPwC+lJk/BsjMbVNcY7s10+cE3lA9fyN73IfTazLzQeC5vRyyGPhaNqwFjoiIY9vVfq8EfTNLKrx6TGa+DOwEZkxJdZ2xr8tIXERjRNDLJu1zRJwKHJ+ZX5/KwjqomT/nE4ETI+K7EbE2Is6csuo6o5k+fxb4QESM0rh672NTU1rXdHTZmI5dXtlmky6p0OQxvaTp/kTEB4Ah4Dc7WlHn7bXPEXEQ8AXgQ1NV0BRo5s95Go3pmwU0fmr7l4g4JTN/0uHaOqWZPp8PfDUzr46ItwJ/V/X5lc6X1xUdza9eGdFPuqTC2GMiYhqNH/f29qPS/q6ZPhMR7wI+A5ydmT+foto6ZbI+Hw6cAjwQESM05jLX9PgHss3+3V6dmS9l5g+AJ2gEf69qps8XAbcBZOa/An001sEpVVP/3lvVK0HfzJIKa4Cl1fPzgG9n9SlHj5q0z9U0xnU0Qr7X521hkj5n5s7MnJmZg5k5SONzibMzc113ym2LZv5u303jg3ciYiaNqZynprTK9mqmz1uAhQAR8SYaQb99SqucWmuAD1ZX35wB7MzMZ9v15j0xdZMTLKkQEX8GrMvMNcBKGj/ebaYxkl/SvYrra7LPVwGvB26vPnfekplnd63omprsc1Ga7PO3gN+OiI3AL4E/zswd3au6nib7/Cng+oj4BI0pjA/18sAtIm6mMfU2s/rcYTkwHSAzv0zjc4izgM3Ai8CH29p+D//eSZKa0CtTN5KkFhn0klQ4g16SCmfQS1LhDHpJKpxBL0mFM+glqXAGvTSOiBiMiMcjYlW1PvgdEXFot+uSWmHQSxM7CViRmfOBn9L4zgOp5xj00sS2ZuZ3q+d/D7y9m8VIrTLopYntuT6I64WoJxn00sRmVWuhQ2N99O90sxipVQa9NLFNwNKIGAaOAq7tcj1SS3pimWKpS17JzIu7XYRUlyN6SSqc69FLUuEc0UtS4Qx6SSqcQS9JhTPoJalwBr0kFc6gl6TC/S/k8eyS9sNdtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(res.pvalue, bins = np.arange(0,1.05,0.05))\n",
    "plt.hist(res.pvalue[TPtests:], bins = np.arange(0,1.05,0.05))\n",
    "plt.legend([r'$H_A$',r'$H_0$'])\n",
    "plt.xlabel('p')\n",
    "plt.axvline(x = 0.05, color = 'red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also visualize the QQ-plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEbCAYAAADTZlM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcHGWdx/HPL5MJCTkIkCgkZDJci8uRyxECGTBAkEMFDIggiOyqWdcTcZVLRVAMHovirsrGY1UICIswgqCYLCIkC4GEmSScgjCBTJAk4iQkhJy//aOqmU6ne6a6U13d1f19v179oo+nup6aDs+v6nme+j3m7oiISP3pV+kKiIhIZSgAiIjUKQUAEZE6pQAgIlKnFABEROqUAoCISJ1SAJCaYWbNZuZm1j+BfQ0ys7vMbI2Z/U+599dHXTrNbFol6yDppAAgZWFmF5jZUjN73cz+amY/NLPdcsocbGZ3ho3oa2Z2n5lNTqh+XzWzG3fiK84E3grs6e7vj6laIolSAJDYmdnngW8CXwB2AyYDzcAfzKwxLLM/MB9YCuwLjALagDlmdngFql2sscCf3X1LpSsiUjJ310OP2B7AMGAdcFbO+0OAlcCHw9c3APfk2f5HwB97+f77gZnAI8Aa4DfAHuFnzYAD/cPXo4A7gVeB54CPhe+fBGwCNod1XVxgX/8Y7q8beAI4NXz/ypztP5Jn268CtwG3AK8BjwHjC+xnFLAhcxzhexOB1UAjsD9wH/C38L3ZwPCssp3AtPD5z4GvZ302FVies69fA6uAF4DPZH12OLAQWAu8Alxb6X9PepT3oSsAidtRwEDg9uw33X0d8DvgXeFbJwD5+s5vBY42s4G97ON84J8JGrMtwPcLlLsZWB6WOxP4hpkd7+6/B74B3OLuQ9x9fO6G4ZXKXcAfgLcAnwZmm9lB7n5FzvY/LbD/08Jj3AO4CWjLXAFlc/cVwEPAGVlvfxC4zd03A0YQ9EYRBKUxBAGmKGbWLzymxcBo4HjgQjM7MSxyHXCduw8jCDq3FrsPSRcFAInbCGC15+8aeRkYmVXu5QJlGggazUJucPfH3X098GXgLDNryC5gZmOAVuBid3/D3TuAnwAfingckwmuWq5x903ufh/wW+CciNsDLHL3TCN+LUFgLDTGcVPmu83MgLPD93D359x9jrtvdPdV4Xe9s4h6ZLwDGOnuV4XH9Dzw43BfEFzRHGBmI9x9nbs/XMI+JEUUACRuq4ERBWbi7E3Q9ZApt3eBMg6sNrPrzWxd+Lgsq8xLWc+XEXSTjMj5nlHAq+7+Wk7Z0RGPYxTwkrtvK3H77eoZfs9yYJSZnZt1XL8Li9wGHGlmo4BjCP4GDwKY2VvM7Fdm1mVma4Eb2fF4oxgb7r878wAuIxjMBvgI8A/A02b2qJm9p4R9SIooAEjcHgI2AtOz3zSzwcDJwJ/Ct+YC+WbPnAU8HJ6hfjzsYhni7t/IKjMm63kTwZnr6pzvWQHsYWZDc8p2hc/7SoO7AhgTdpvk2z6KN+sZfs8+wAp3n511XCcDuHs3QXfTWQTdPze7e6aOM8P6jgu7Z84j6BbKZz2wa9brvbKevwS84O7Dsx5D3f2UsA7Puvs5BF1e3wRuC383qVEKABIrd19DMEj6H2Z2kpk1mlkzQV94ZgCTsMxRZna1me1hZkPN7NPAPwFX9LGb88IppLsCVxH0lW/NqcdLwP8BM81soJmNIzjDzez/FaA5p4HPtoCgMf1ieAxTgfcCv4r2lwDg7WY2PbwaupAgMPbWrXITwfjGGeHzjKEEg83dZjaaYHZVIR3AKeHfdK9wvxmPAGvN7OLwPoYGMzvUzN4BYGbnmdnI8GqlO9xmK1KzFAAkdu7+LYKuhe8QzIB5geCsdFrYb4+7P0vQRz+eYBZLN/A14H3uPqePXdxAMNvlrwT96p8pUO4cgplBK4A7gCuyvjszAP03M3sszzFsAk4luGpZDfwQON/dn+6jbtl+A3wA+DvB2MP0cDygkDuBA4FX3H1x1vtXApMIZj3dTc4Ae44bCAZ5OwmuKG7JOqatBEFsAsFvsppgXCRzf8ZJwBNmto5gQPhsd38jyoFKOlnPVaZIeZjZPxM0YlPc/cUCZfYhODu+opdZNZjZ/cCN7v6TctQ1Lmb2VeAAdz+v0nURKaTst8yLuPvPzGwzwRTRvAHA3Zeb2cnAaWY2JJw2KiJlpAAgiXD3GyKUWUpwZ7CIJEBdQCIidUqDwCIidUoBQESkTlX1GMCIESO8ubm50tUQEUmNRYsWrXb3kX2XrPIA0NzczMKFCytdDRGR1DCzZVHLqgtIRKROKQCIiNQpBQARkTqlACAiUqcUAERE6lRVzwISEakXbe1dfPveZ1jRvYFRwwfxhRMP4vSJxaw/VDwFABGRCmtr7+LS25eyYXOw/EJX9wYuvT1Ii1XOIKAuIBGRCvv2vc+82fhnbNi8lW/f+0xZ96sAICJSYSu6NxT1flwUAEREKmzU8EFFvR8XBQARkQr7wokHMaixYbv3BjU28IUTDyrrfjUILCJSYZmB3pqdBWRmB5G1QDWwH/AVd/9eUnUQEalWp08cXfYGP1diAcDdnwEmAJhZA9AF3JHU/kVEZHuVGgM4HviLu0dOWyoiIvGqVAA4G7i5QvsWEREqMAhsZgOAU4FLC3w+A5gB0NTUlGDNREQqk5KhUipxBXAy8Ji7v5LvQ3ef5e4t7t4ycmSkVc1ERGKRScnQ1b0BpyclQ1t7V6WrVhaVCADnoO4fEalClUrJUCmJBgAz2xU4Abg9yf2KiERRqZQMlZJoAHD31919T3dfk+R+RUSiqFRKhkpRKggRkVClUjJUilJBiIiEKpWSoVIUAEREslQiJUOlKACISE2op/n7cVEAEJHUq9SSimmnQWARSb16m78fFwUAEUm9epu/HxcFABFJvXqbvx8XBQARSb16m78fFw0Ci0jq1dv8/bgoAIhITain+ftxUReQiEidUgAQEalTCgAiInVKAUBEpE4pAIiI1CkFABGROqUAICJSp5JeE3i4md1mZk+b2VNmdmSS+xcRkR5J3wh2HfB7dz/TzAYAuya8fxERCSUWAMxsGHAMcAGAu28CNiW1fxER2V6SXUD7AauA/zazdjP7iZkNTnD/IiKSJckA0B+YBPzI3ScC64FLcguZ2QwzW2hmC1etWpVg9URE6kuSYwDLgeXuviB8fRt5AoC7zwJmAbS0tHhy1RORaqI1fssvsQDg7n81s5fM7CB3fwY4Hngyqf2LSHWI0rBrjd9kJH0fwKeB2Wa2BJgAfCPh/YtIBWUa9q7uDTg9DXtbe9d25bTGbzISDQDu3uHuLe4+zt1Pd/e/J7l/EamsqA271vhNhu4EFpHERG3YtcZvMhQARCQxURt2rfGbDAUAEYlVW3sXU665j30vuZsp19y3Xf9+1Ib99ImjmTn9MEYPH4QBo4cPYub0wzQAHDOtCSwiselr9k4xi7drjd/yUwAQkdj0NsibaczVsFcPdQGJSGw0eyddFABEJDaavZMuCgAiEhvN3kkXjQGISGyKGeSVylMAEJFYaZA3PRQARKRoytRZGxQARKQoytRZOxQARKSgfGf6Ueb6SzooAIhIXoXO9HMb/wzN9U8fTQMVkbwKnek3mOUtr7n+6aMAICJ5FTqj3+quuf41QgFARPIqdEafycypTJ3ppzEAEcnrCycetEOff+ZMX3P9a0OiAcDMOoHXgK3AFndvSXL/IhKd7uqtfZW4AjjW3VdXYL8iUiSd6dc2jQGIiNSppAOAA38ws0VmNiPhfYuISJaku4CmuPsKM3sLMMfMnnb3B7ILhIFhBkBTU1PC1RMRqR+JBgB3XxH+d6WZ3QEcDjyQU2YWMAugpaXFk6yfSD1Tgrf6U3QXkJmda2b/UMJ2g81saOY58C7g8WK/R0Til0n70NW9Aacn7UNbe1elqyZlVMoYwCrgh2b2gJndbmbXRNzurcA8M1sMPALc7e6/L2H/IhKz3hK8Se0qugvI3f9gZu9098vDM/qZEbd7Hhhf7P5EpPy0mHt9KnUW0DAzmwRsBAbHWB8RqQAt5l6fSg0AFwGtwPXAvfFVR0QqQYu516eSZgG5+2bg+zHXRUQSkm/Gz8zph2kWUJ0pKQCY2Q+Awe5+gZm9y93/EHO9RKRMCi30MnP6Ycy/5LgK106SVGoX0Cbg+fC5/sWIpIhm/EhGqQHgdWA3M2sEdLuuSIpoxo9klBoAZgJ/AX4A3BRfdUSk3DTjRzJKTQVxHbAA+Li7b4uxPiJSBtmDvrsNaqSxwdi8tSfTimb81KdSZwF9xMyOBn5sZr9x9ztjrpeI7ITcBn/9pi1vNvjdGzbT2M/YfddGul/frBk/dazUWUDHAIcCG4CvAwoAIlUid5ZP94bNO5TZvM3ZdUB/2r/yrqSrJ1Wk1C6gVoIsnrPdfU2M9RGRnZRvlk8+GvSVUgPAFnefF2tNRGSnZLp9uiI27Br0lVIDwGlm9jowx901eVikQrIbfSNYci8KDfoKlB4ApgMTgfeZ2QHu/tEY6yQivSjU6PfW+Df2M4YM7K9BX9lOqQHgK/SkgjgxzgqJyI5KafQzRqvBlwJKDQCbgFfC58eijKAiZZM7q6eYdVJHDx+k/D5SkFJBiFS5qLN6cqmfX/pS6hXAH4F/IFgPYHZ81RGRXMVM18x0D6nbR6IoNQC8i2A9gBPDx91RNzSzBmAh0OXu7ylx/yI1L9Pv31eXjxp9KVWpAWA4cDHwReAjRW77WeApYFiJ+xapebn9/rnU6EscSg0AVwFvc/dnzCxyMjgz2wd4N3A1wbKSIpKjrb2Lz9+6mK2e/9xfjb7EpdRkcMvNbG34/JIiNv0ewVXD0FL2K1LrMmf+hRp/A83qkdiUmgzuZmCNmQ0CfuruD0TY5j3ASndfZGZTeyk3A5gB0NSkCUZS27Ln9zeYFWz4M5S+QeJUahfQU+5+FYCZ/ZAgMVxfpgCnmtkpwEBgmJnd6O7nZRdy91nALICWlpZipjyLpEZbexdfvfOJ7TJ19tX4a1qnxK3oAGBmPwbeZmbdwBKgO8p27n4pcGn4HVOBf8tt/EXqQV8DvPk0mDFz+mHq95dYFR0A3P1jZjaQIBfQ4cCesddKpEb1NcCbz6DGBjX+UhaldgFdBvwjsJ5gfeCiuPv9wP0l7lsklfoa4M1HM36knEoNAHu4+/vNbADwXeCTMdZJpOYUe+avs35JQqkBYKOZTQKWAoNjrI9ITck32FtIP4NtrrN+SU6pAeBy4NPh45b4qiNSG4pp+BvM+PezxqvBl8SVeiPYG8C3Y66LSE34UttSZj/8YqS0zerqkUqKHADM7GTgqwR5gBYD17r7w2Wql0jqFHPWD5raKZVXzHoAPyTI3zOZ4Eat75jZOWWplUjKZGb4RG38BzU2qNtHKq6YLqBX3H1++HyumT0ELABujr9aIulSzKItu+/ayBXvPUSNv1RcMQGg08y+Dlzl7puAzcBr5amWSLpEWbRFDb9Um2ICgAPTgY+Z2bMES0HONrMD3f3ZstROJAXa2rvo10siNwPOndzE108/LNmKifQhcgBw93MAwjQQhwLjgXHAT8xsP3cfU54qilSvvu7u1Vm/VLNSpoG+DzgN2EowiHy9u2scQOpSob5/ze2XNCglAEx197MzL8zsB2ggWOpUV4G+/23uavyl6pUSAHYxs3cDLwFjAK1QIXUnM+e/EC3cImlQzH0AGZ8AdgdOIbgp7FOx1kikyvU1599AC7dIKpSyHsDrwI1lqItIKlx51xO9zvl3UPePpEIpVwDbMbOL46iISBq0tXfx99d7v9t3tLp/JCVKWRLy1uyXwATgm7HVSKQKRc3zo3V7JU1KGQRe6+4fzbwwsx/FWB+RqtPW3sUX/mcxm7f1nt9Tc/4lbUoJAFfnvL48ykbhDWQPALuE+73N3a8oYf8iiWlr7+Jzt3bQ10Jewwc10v6VdyVTKZGYFJMO+qKs55mna4BFwKsRvmIjcJy7rzOzRmCemf1OKaWlGhWT2nlQYwNfPfWQBGolEq9irgBawsdd4et3A48CHzez/3H3b/W2sbs7sC582Rg+oq+OLVJmxebzB+X0l3QrZhbQnsAkd/+8u3+eIBiMBI4BLojyBWbWYGYdwEpgjrsvKLK+ImWR6ecvpvFvbFC6B4nJqlVwyy3w0Y/C5Mn02ecYk2KuAJqATVmvNwNj3X2DmW2M8gXuvhWYYGbDgTvM7FB3fzy7jJnNAGYANDU1FVE9kdJ9+95n+hzkzWYG3z5Tjb+UaMMGmD8f5swJHu3twfu77QbHHQdr1wbPy6yYAHAT8LCZ/YZg+ud7gJvNbDDwZDE7dfduM7sfOAl4POezWQQrjtHS0qIuIim7tvaugjl98mlsMDX+Upxt26CjA+bODRr8efPgjTegsRGOPBK+9jU44QR4+9uhf0lLtZekmHTQXzOze4BWggDwcXdfGH58bl/bm9lIYHPY+A8CpqH7B6TCMl0/UWmqp0S2bFnQ2M+dC//7v7B6dfD+oYfCxz8eNPjHHANDhlSsisWGmi3ANoLB2+idpYG9gV+YWQPB2MOt7v7bIr9DJBbFDviepwVdpC/d3fDHP/Y0+s+G62TtvTeccgpMmxY89t67svXMUsw00M8CHwN+TXAFcKOZzXL3/4iyvbsvASaWVEuRGEW9sQt0xi+92LQJHn64px//0UeDrp7Bg2HqVPjkJ4MG/+CDg0GjKlTMFcBHgCPcfT2AmX0TeAiIFABEKq3Ys/7Rwwcx/5LjylwrSQ13eOKJnn78P/0J1q+Hfv3g8MPh8suDBn/yZBgwoNK1jaSYAGAEq4BlbA3fE6l6xZz1QzDQq5w+wooVQf99plvn5ZeD9w88ED784aAff+pUGD68otUsVTEB4L+BBWZ2R/j6dOCn8VdJJH5X3vVE5MZfUzzr2Lp1wZl95iz/iXDRnxEj4PjjgwZ/2jQYO7ay9YxJMbOArjWzPwFTCM78/8nd28tWM5EYBIu3LGHD5m2RymuKZ53ZsgUWLuw5w3/oIdi8GQYOhKOPhvPPDxr98eODrp4aU9QsIHdfRJD7R6SqFdvwgwZ864I7PPdcT4N/332wZk1w2TdxIlx0UXCGP2UKDKr9dR36DABm9hr5c/YYQYqfYbHXSqQEpeTyAZ3117xVq4KGPtPoL1sWvD92LLz//UGDf/zxQTdPnekzALj70CQqIrIzvtS2lBsffrHo7XTWX4P6SrNw8cVBo3/AAVU7PTMppawI9l53v6vvkiLJOPfHDzH/L1Eykvcw4LsfmKCGvxZUaZqFNCh1QRgFAKkKpTT+AOdOblLjn2aF0iwcckjVpFlIg1ICQH1fM0nFlTLAmzF4QANXv0/5+1OntzQLJ5/cMz2zitIspEEpAUAZOqViSu3r36V/P755xjg1/GmRnWZh7lx45JHt0yx84hNBo1/FaRbSQB1ikhqlNP79DD54hBK5VT13ePLJnoHbGkizkAYKAJIKxTb+U/bfg9kfO7KMNZKdVuNpFtKglADwSuy1EOmFGv8aUWdpFtKg6ADg7ieUoyIi+bS1dxXV+CtvfxUplGZhl13qIs1CGqgLSKpWW3sXn7ulI1JZze6pAoXSLABMmgSf+1zQ4NdJmoU0UACQqtTW3sVFt3b0OeWsn8G1Z+mGropZvXr7fnylWUgVBQCpShf/egl9ZW821Pgnrq80C1/8YnCWrzQLqZBYADCzMcAvgb0I1hWe5e7XJbV/SY9zf/wQG7f0fZOXUjkkYNs2WLy4p8HPl2Zh2jRoaVGahRRK8hfbAnze3R8zs6HAIjOb4+5PJlgHqXJfalsaKbXDeUrlUD7LlvXM1FGahZqWWABw95eBl8Pnr5nZU8BoQAFAgOh5fTTTJ2aZNAuZRl9pFupGRa7ZzKwZmAgsqMT+pfqo8U+Q0ixIKPEAYGZDgF8DF7r72jyfzwBmADQ1NSVcO6mEqI3/lP33UONfir7SLFx2WdDgK81C3Uk0AJhZI0HjP9vdb89Xxt1nAbMAWlpalHiuxp1w7f08u3J9n+V0d2+RXn65p0snX5qFadPg2GOVZqHOJTkLyICfAk+5+7VJ7VeqV9TGf1BjPzX+fVm3Dh54oOcsP5NmYc89g8ZeaRYkjySvAKYAHwKWmlnm9s7L3P2eBOsgVSJq42/AzOnjyl+htMmkWcic5SvNgpQgyVlA89BiMnWv2MRumusfUpoFKQPduSGJKHbpRq3ZS0+ahcxZvtIsSMwUAKSsSlnEpX8/4zvvH19/jf8bbwR32mbO8tvbgzN/pVmQMlEAkLKJ2s+f7cC3DGbORVPLU6Fq01eahauuUpoFKSv9q5Ky+FLbUjX++SjNglQRBQCJXbGLuEANN/5R0iwcfzyMGlXZekpdUgCQWBU72As1lt4hk2Yh0+Bnp1l45zuVZkGqigKAxKbYxr8mGv7sNAtz58L99yvNgqSGAoDEopjGP/UNv9IsSI1QAJCdVkzj/700zu1XmgWpUQoAslOiLuACQUK3VDT+SrMgdUIBQEpS7A1eB75lcPUmdOstzcLEiUqzIDVLAUCKUsqdvVXZ518ozUJTE5x5ZtDgH3ccjBxZ2XqKlJECgERSyl29UEWNf29pFo49VmkWpC4pAEif3nb5Pbyxtfi1eSq6glehNAv9+8NRR8GVVwYNvtIsSB3Tv3wpqJSbujIqcub/4os9DX6+NAvTpgU3YynNggigACAFlHrWDwme+Xd3BzdeZRp9pVkQKYoCgOxgZxr/sp75K82CSKwUAORNbe1dXHhLR98Fc5Qtf7/SLIiUVZKLwv8MeA+w0t0PTWq/Ek0p/f279O/HN88YF2/D31uahcwNWEqzIBKLJK8Afg78J/DLBPcpfSil4R+2SwNLrjwpngpkp1mYOxcefzx4P5NmIZNqQWkWRGKX5KLwD5hZc1L7k/xK7ebJ2Om8/Vu2wKJFPQO3uWkWPvShoNGfMEFpFkTKTGMAdWTcFb9n7catJW9fUiK3TJqFTLdOvjQL06ZBa6vSLIgkrOoCgJnNAGYANDU1Vbg2tWFn5vMDDGwwnr76lOgbrF4dNPSZs3ylWRCpSlUXANx9FjALoKWlpbS5iALsfHcPRGz8M2kWMmf5SrMgkgpVFwBk5+3sGX/GW4cOYMHlJ+z4QXaahblz4cEHlWZBJIWSnAZ6MzAVGGFmy4Er3P2nSe2/Xhxw6d1sieG6aYcbupRmQaTmJDkL6Jyk9lVv4jrj3+6Gru5uaGvbMc3CXnsFaRYyUzSVZkEktXR9nmJxNPxT9t8jWKhl0yZYsADu+C/4RIE0C9OmBWf86scXqQkKACmzM3l6sg0b0I8lH2gK+vDfc/X2aRbe8Q6lWRCpAwoAKbCz8/czRq57ldbODs5f/ywT/7wIvrYi+EBpFkTqkgJAlSp1Ba5sgza9wREvLeXozg6mdHbwttXhfHylWRARFACqzs6c7Tds28phf32O1s52Wjs7mNT1NAO2bWFT/wEMmHoMnPAJpVkQkTcpAFSBkht9d5r/voLWZYtp7WznqGVLGLYxuGp4/K3786sp0zn/yx9lgNIsiEgeCgAVUup8/d1fX8NRy5bQ2tnO0Z0d7LN2JQDLh43knoOmMK95Alve+U6uv+gUlHNbRHqjAJCgL7Ut5caHXyxqm122bKJl+ZO0dnbQ2tnOIa88Tz+ctQN25aGx47j+iDOY1zyBzt1HMeWAPYMpnSIiESgAlNkRV8/hldc2RS5vvo2DV75Aa2c7UzoXc/jyJxi4ZROb+zXw2Ki38d3WDzKveSJL9j6QwYMGxJeXX0TqjgJAGTRfcndR5UetXUnrCx20LuvgqGWLGfF6kC75mRFN3DT+JB7cdyILxhzK6wOCfvyS0jKLiORQAIhJMTdoDd24niOXLaG1s4MpyzrY/9UuAFYO3p0/7TuJ+c0TmDd2AiuH7vnmNju9EIuISA4FgBIVk4ahcetmJqx4htbODo7ubGf8y8/S4NtY3ziQBWMOZfaEU5jXPJ4/jxi7Q5qFHZKyiYjERAGgCJG7dtw5cPWLtC7roLWzgyNeepwhmzaw1fqxZK8D+cHk9zNv34m0jzqIzQ2NO2yuRl9EkqAA0Iti+vIzaRZalwV33e61Lrg6eGH3vbnjkGOZ1zyBh5rGsXZg/nTJnde8O5Y6i4hEpQCQY99L7iZKT36hNAuvDhrG/40dz4PNE5jfPIHlu7214HcUvdSiiEiMFACIdqbfb9tWxuVJs7CxoZFH9jmEtkOO5cHmCTz51v1wK5xmoeAqWyIiCavLABBpANedsd0vc3R4A9aRy5awW1aahZ+1nMq85ok8us/BbGzcpdev0rRNEalGdRMAoqReGL5hLVM6F+dNs/C7g6Ywf+x45jdP4NVdd+tzfzrTF5Fql2gAMLOTgOuABuAn7n5NOffXV9fOLls28fblT755lp8vzcL85gm8sPuoSKtgaSBXRNIkyUXhG4AfACcAy4FHzexOd38y7n0VavgzaRamdAbTM3tLs7C1X0Of+1GDLyJpluQVwOHAc+7+PICZ/Qo4DYg1AOQ2/sWmWeiNunVEpJYkGQBGAy9lvV4OHFGOHQ3fsJaLHpwdOc1CbzRVU0RqVZIBIF8n+g7DsmY2A5gB0NTUVNKOXm8cxCnPzGPJXgdy04STebB5Qt40C4Woa0dE6kGSAWA5MCbr9T7AitxC7j4LmAXQ0tJSwpIpsKl/I4d/8pdsi9CPn6FGX0TqTZIB4FHgQDPbF+gCzgY+WK6dRWn81eiLSD1LLAC4+xYz+xRwL8E00J+5+xNx76fzmncXnAWk/nwRkR6J3gfg7vcA95R7PzqzFxHpW+GkNSIiUtMUAERE6pQCgIhInVIAEBGpUwoAIiJ1ytxLutcqEWa2Cli2E18xAlgdU3UqqVaOA2rnWHQc1adWjmVnj2Osu4+MUrCqA8DOMrOF7t5S6XrsrFo5DqidY9FxVJ9aOZYkj0NdQCIidUoBQESkTtV6AJhV6QrEpFaOA2rnWHQc1adWjiWx46jpMQARESms1q8ARESkAAUAEZE6lfoAYGYnmdkzZvacmV2S5/NdzOyW8PMFZtacfC2jiXAsF5jZKjPrCB8frUQ9+2JmPzOzlWaKtF6pAAAGPklEQVT2eIHPzcy+Hx7nEjOblHQdo4hwHFPNbE3W7/GVpOsYhZmNMbM/mtlTZvaEmX02T5m0/CZRjqXqfxczG2hmj5jZ4vA4rsxTpvxtl7un9kGwrsBfgP2AAcBi4OCcMp8Arg+fnw3cUul678SxXAD8Z6XrGuFYjgEmAY8X+PwU4HcEy4ROBhZUus4lHsdU4LeVrmeE49gbmBQ+Hwr8Oc+/rbT8JlGOpep/l/DvPCR83ggsACbnlCl725X2K4DDgefc/Xl33wT8Cjgtp8xpwC/C57cBx5tFXBw4WVGOJRXc/QHg1V6KnAb80gMPA8PNbO9kahddhONIBXd/2d0fC5+/BjwFjM4plpbfJMqxVL3w77wufNkYPnJn5JS97Up7ABgNvJT1ejk7/mN4s4y7bwHWAHsmUrviRDkWgDPCS/TbzGxMns/TIOqxpsGR4WX878zskEpXpi9hN8JEgjPObKn7TXo5FkjB72JmDWbWAawE5rh7wd+kXG1X2gNAvmiYG0WjlKkGUep5F9Ds7uOAufScHaRNWn6TvjxGkHdlPPAfQFuF69MrMxsC/Bq40N3X5n6cZ5Oq/U36OJZU/C7uvtXdJwD7AIeb2aE5Rcr+m6Q9ACwHss+C9wFWFCpjZv2B3ajOy/o+j8Xd/+buG8OXPwbenlDd4hbld6t67r42cxnvwXKnjWY2osLVysvMGgkazNnufnueIqn5Tfo6ljT9LgDu3g3cD5yU81HZ2660B4BHgQPNbF8zG0AwUHJnTpk7gQ+Hz88E7vNwVKXK9HksOX2ypxL0f6bRncD54cyTycAad3+50pUqlpntlemTNbPDCf5/+ltla7WjsI4/BZ5y92sLFEvFbxLlWNLwu5jZSDMbHj4fBEwDns4pVva2K9FF4ePm7lvM7FPAvQSzaH7m7k+Y2VXAQne/k+Afyw1m9hxB9Dy7cjUuLOKxfMbMTgW2EBzLBRWrcC/M7GaCmRgjzGw5cAXBIBfufj1wD8Gsk+eA14F/qkxNexfhOM4E/tXMtgAbgLOr9ORiCvAhYGnY5wxwGdAE6fpNiHYsafhd9gZ+YWYNBAHqVnf/bdJtl1JBiIjUqbR3AYmISIkUAERE6pQCgIhInVIAEBGpUwoAIiJ1SgFARKROKQBIzTCzdX2XEpEMBQARkTqlACA1ycwuMrPHw8eFWe9/2cyeNrM5Znazmf1bEd/5q3CBjgVmtszM3l2e2oskI9WpIETyMbO3E6QyOIIgo+ICM/sTQYqNMwhSCPcnyBq5qIivHg+0ufsHzKwVuBa4O866ZzOz3d397+X6fhEFAEkNM5sL7JXno8vd/TdZr1uBO9x9fbjd7cDRBFe8v3H3DeH7d4X/3Q+4HNjN3c8M3xsM/BDYRJCp8XZgBJBZuu9JYPc4jy+P79JHviczOxF4i7vfUOa6SA1SAJDUcPdpEYsWWjUp7/vu/jzwETO7Levt6cBt7n6Xmd1CsPTgs+7+Rvj5JIJlOzPpia8EdiVYzvMK4FvAl4GrgEsJllu8FziIIFlZQ842nwm3cWBZuL+3hV1U1+WUvZDg6mMtwVXOJyP+XUS2ozEAqUUPAKeb2a7hmfz7gAeBecB7wwW5hwC99eHvQ88KWVsJun+awm0HEzTI3w0/nwEMAroJ1nldBbwI/DtBw94C3OzulxKs/rRn7jbAvxJcnXze3b8PrAZudPfvFCj7C3e/DNgFeKb0P5XUM10BSM1x98fM7OfAI+FbP3H3dgAzu5PgzH0ZsJBgmb18lhMEgQ6CE6XxwGyC7qBhwDfcfX5YdiLwycxiPWFw2Q/Y4u7rzOwdwP+FZXdz91VmlrvNfwM/ytr/uLCe+b7/58B/mdlQYHUVpjqWlFA6aKkrZjYkbJR3JbhSmEEQDK4GTiAIFjPDs/z/BN4guHL4F+Bj7r7D2baZvRf4IMEVwwME3UdXAmcRLPTzLwRn9P0JupX+N2eb+wjO5E8lyPs+k2Ac4wzgGuCAnLK7AycS5O1f5+5fjPNvJPVDAUDqipndBBwMDCToRpkZcbsuYIy7bythnze4+4eK3U6k3BQARETqlAaBRUTqlAKAiEidUgAQEalTCgAiInVKAUBEpE4pAIiI1CkFABGROqUAICJSpxQARETq1P8DvKo03CwEEtAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x = sorted(-np.log10(np.arange(1,Ntests+1)/Ntests)), \n",
    "            y = sorted(-np.log10(res.pvalue)))\n",
    "\n",
    "plt.plot(sorted(-np.log10(np.arange(1,Ntests+1)/Ntests)), \n",
    "         sorted(-np.log10(np.arange(1,Ntests+1)/Ntests)), 'r-')\n",
    "\n",
    "plt.title('QQ-plot of p-values')\n",
    "plt.xlabel(r'$-\\log_{10}p_{expected}$')\n",
    "plt.ylabel(r'$-\\log_{10}p_{observed}$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the number of FPs consistent with the false positive rate?\n",
    "\n",
    "$$ FPR = \\frac{FP}{FP+TN}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FPR = 0.051\n"
     ]
    }
   ],
   "source": [
    "fpr = fp/(fp + sum(res.pvalue[TPtests:] >= 0.05))\n",
    "print('FPR = %.3f' % fpr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Bonferroni correction\n",
    "This is a more stringent correction that adjusts the significance cutoff $\\alpha$ to account for the number of tests.\n",
    "$$\\alpha_{\\text{Bonferroni}} = \\frac{\\alpha}{N}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13 significant results after Bonferroni correction\n",
      "13 TPs after Bonferroni correction\n",
      "0 FPs after Bonferroni correction, FPR = 0.000\n"
     ]
    }
   ],
   "source": [
    "alphaBonf = 0.05/Ntests\n",
    "bonfSignificant = sum(res.pvalue < alphaBonf)\n",
    "print('%d significant results after Bonferroni correction' % bonfSignificant)\n",
    "\n",
    "bonfTP = sum(res.pvalue[:TPtests] < alphaBonf)\n",
    "print('%d TPs after Bonferroni correction' % bonfTP)\n",
    "\n",
    "bonfFP = sum(res.pvalue[TPtests:] < alphaBonf)\n",
    "bonfFPR = bonfFP/(bonfFP + sum(res.pvalue[TPtests:] >= alphaBonf))\n",
    "print('%d FPs after Bonferroni correction, FPR = %.3f' % (bonfFP, bonfFPR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How did the number of true positives change after Bonferroni correction? How about the number of false positives?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}