NSCI801_Intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# NSCI 801 - Quantitative Neuroscience\n",
    "## Introduction\n",
    "Gunnar Blohm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Outline\n",
    "\n",
    "* Why quantitative Neuroscience?\n",
    "* Course overview & materials\n",
    "* The research process\n",
    "* Study design\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Why Quantitative Neuroscience?\n",
    "\n",
    "* We want to quantify observations\n",
    "* But data is corrupted by noise\n",
    "* Certain things are not directly observable (latent)\n",
    "    + we need models!\n",
    "* Ultimately we want to identify **_causal_** relationships"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Why Quantitative Neuroscience?\n",
    "\n",
    "* We want to quantify observations\n",
    "    - questionnaires\n",
    "    - measurements\n",
    "    \n",
    "But: such observations are variable..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Why Quantitative Neuroscience?\n",
    "\n",
    "* We want to quantify observations\n",
    "* But data is corrupted by noise\n",
    "    - noise in the process / system\n",
    "    - noise due to the measurement\n",
    "    - noise due to A/D conversion\n",
    "    - noise due to post-processing\n",
    "\n",
    "Thus: we need ways to infer reality from noisy data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Why Quantitative Neuroscience?\n",
    "\n",
    "* We want to quantify observations\n",
    "* But data is corrupted by noise\n",
    "* Certain things are not directly observable (latent)\n",
    "    - e.g. we cannot measure your thought process (yet), only the outcome!\n",
    "    - e.g. we cannot measure inflammation, only the body's reaction\n",
    "    - we often want good \"measures\" of latent variables\n",
    "* Ultimately we want to identify **_causal_** relationships\n",
    "\n",
    "Solution: we need models that causally link latent variables to measurable quantities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Course overview & materials\n",
    "\n",
    "* [course web site](http://www.compneurosci.com/NSCI801.html)\n",
    "* we will use Google Colab - you need a Google account!\n",
    "* all lecture materials will be in Python & Markdown\n",
    "* slides / tutorials will be shared on [GitHub](https://github.com/BlohmLab/NSCI801-QuantNeuro)...\n",
    "* download code from Github into Colab: File>>Open Notebook>>Github..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### for those interested...\n",
    "\n",
    "* Jupyter Notebook has a \"slide\" option that produces HTML5 slides\n",
    "    - [read more about it](https://medium.com/@mjspeck/presenting-code-using-jupyter-notebook-slides-a8a3c3b59d67)\n",
    "* install Reveal.js - Jupyter/IPython Slideshow Extension (RISE)\n",
    "    - [install and use](https://rise.readthedocs.io/en/maint-5.5/index.html)\n",
    "    \n",
    "But you **don't** need any of this!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Course overview & materials\n",
    "\n",
    "#### Goals of the course:\n",
    "* hands-on skills in signal processing, basic and advanced statistics, data neuroscience (machine learning) and model fitting methods\n",
    "* gain intuitive understanding of these topics\n",
    "* introduction to scientific programming in Python\n",
    "* familiarization with open science framework approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Course overview & materials\n",
    "\n",
    "#### [Specific topics](http://www.compneurosci.com/NSCI801.html):\n",
    "* intro to Python & Colab\n",
    "* signal processing\n",
    "* statistics and hypothesis testing\n",
    "* models & data neuroscience\n",
    "* causality, reproducibility, Open Science"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "<img style=\"float: center; width:500px;\" src=\"stuff/research-process.png\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Research design:\n",
    "* what is power?\n",
    "* what is effect size?\n",
    "* how to determine sample size?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Research design:\n",
    "* what is power?\n",
    "\n",
    "    Power calculations tell us how many samples are required in order to avoid a type I (false positive) or a type II (false negative) error\n",
    "    \n",
    "    Typically in hypothesis testing, only type II errors are considered: For a type II error probability of $\\beta$, the corresponding statistical power is $1 − \\beta$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "<img style=\"float: right; width:300px;\" src=\"stuff/effect-size.png\">\n",
    "\n",
    "#### Research design:\n",
    "* what is effect size?\n",
    "\n",
    "    Quantification of the difference between two groups\n",
    "    \n",
    "    E.g. Cohen $d=\\frac{\\mu_1-\\mu_2}{\\sigma}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Let's play - effect size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVyTV/Y/8A8JYQdRQFRAQAQFUURIUBEVUBE3XLpo96lFvvPTdrrM1E470860022mznSdLna6L2jrUqkLqCgisgRlFdAEI4K4gCD7EpL7+4OCxaLJDYmBcN6v13lJkueE8wAews197jUDwEAIIcRkCYxdACGEEMOiRk8IISaOGj0hhJg4avSEEGLiqNETQoiJMzd2Af25evUqKioqjF0GIYQMGZ6enhg9enS/jw3KRl9RUQGxWGzsMgghZMiQSqW3fIyGbgghxMRRoyeEEBNHjZ4QQkwcNXpCCDFx1OgJIcTEadXoY2JiUFZWBplMhs2bN9/yuNDQUHR1dWHNmjXcuYQQQgyH3S4EAgGTy+XM29ubiUQilp+fz/z9/fs97vDhw2zv3r1szZo1XLk3h1Qq1XgMBQUFBcWNuF3f1DiPXiKRQC6XQ6FQAAASExMRFxeH0tLSPsc9/vjj2LFjR5/579rmEmIslrY2CJg3By6eHlCr1VCczIcirxBqlUpjrkhkjqioaQgL8wMAFBdfwP79J9HW1mHosgnhorHRu7m5obKysvd2VVUVwsLC+hwzbtw4rFq1ClFRUX0avTa5PeLj47FhwwYAgLOzM99ZEMJJIBRi/iP3YWHCo7CwturzWG1lFXa88k+czbz1BSirV8/Gu+9twLhxTn3ub2howXObv8AnnySDMWaQ2gnhpXGM3szM7Df33fwD/Pbbb2Pz5s1Qq9XcuT22bt0KsVgMsViM2tpaTWURojNLWxv836fvYemT/w9lGVl474EN+GNQOP4sicJXf/wLmEqNhE/excL/e/Q3uWZmZnjvvQT8uOPPuHSpHiuWvwxrq9UQmcchOuoFSKUyfPjRRvy448+wsBiUF56TYUjjT2JVVRU8PDx6b7u7u6O6urrPMaGhoUhMTATQ/Wp8yZIl6Orq0iqXkDvJ0tYGCR+/A/cpk/Hd8y/jZNL+3sc629pQkHwYp4+k464XN2PxxniYW1hg/7sf9R6zdesmPLp+Eba8tQvPPfcFVKobL26OHCnEkSOFePLJOPz7P49h1+6/YNXKf6Czs+uOniMh/bntAL9QKGTl5eXMy8ur9w3VgICAWx7/+eef974Zy5sLLd5UoKAYSDz8n9fZP/PSWWDUvNseZ2Zmxu56cTPbUpTJQpbHMgDsT39azdQsif397/dr/Dzr1y9iapbEPvjg90Y/Z4rhERr6puYniI2NZWfOnGFyuZw9//zzDABLSEhgCQkJvzn2143+VrkDLJiCQqeY/8j9bEtRJpv70FqtjhcIhez3/3ufvSE9ytY8uJgpu3az7xOf1frzvfnmI0zNktgDD0Qa/dwpTD8G3OgHWcEUFNwx2tuTvXnqGHv4369x5dmNGsn+kZbEapp2s9Kyj5itrZXWuUKhgB1Ne53V1X/PxowZafSvAYVpx+36Jl0ZS0yemZkZ7nppMzpaWrHjH//iym2uq4fb1Vw42Qnx3x2laGlp1zpXpVLjsfXvwsrKAm+/s4G3bEL0hho9MXkzli2GT0gwft7yPprr6rly/f098MCqYKQV12PsgjjYOztpTvoVufwS/vHKNtxzzxwsWDCdK5cQfaFGT0yaUCTC4k3xqDxdipzdP3Pnv/raQ2hubseG+1+FuUiERf1MudTkrbd2QqG4gtffeLjfKceEGBo1emLSZt+7GqPGjcXetz/kzp01azJWrpyJf/1zJ2SFZ5D5426ErVkB5/HuXM/T2dmFl178FiEhE3H33eHcdRAyUNToickyt7RE9GMP4WyWFLKsW1/leisv/e0+XL5cj7ff/gkAcPDjz9DVqcTCBP5X9d99l4bCQgX+9vf76VU9ueOo0ROTJVm5FPZOo3Dwo8+4c2fM8MGiRcH4z793o7W1e+2a5mv1yN6xB8GxC+E4xpXr+dRqNV57dTsmT3ZHXFz/y4AQYijU6IlJ6lnLpqKgGOdO5nPn/+nZNWhoaMHHHx/oc3/aV98DAOY9tI77OXfsOIHy8kt4dvMazQcTokfU6IlJmrpgPpzc3ZD62TfcuV5errjrrtn46MP9aGxs7fPY9ctXcGpfCsLWrIC1gz3X86pUamx5axdmzpyMOXMCuOsiRFfU6IlJCl+7BteqLuL00XTu3N//PhaMAe+9l9Tv48e+ToSljTXEcUu5n/uLLw6jrq4J/28jfy4huqJGT0zOmIkT4BMajBPbdoHdtKKqJlZWFnh0/ULs3p2F6uq6fo+pPiPD+fwizL5nFfcbq+3tnfji80NYs2Y2XF0duXIJ0RU1emJyZt+7GsqODkh1mDd/zz1z4OTkgA//u++2x2Uk/ggXr/HwnSm+7XH9+eij/RCJzPHYY4u4cwnRBTV6YlIsrK0Qsmwx8g8cRsv1Bu78+A0xKC2txJEjhbc9riDlCJrr6jHzrjjuzyGXX0JKSh7iNyyGQED/BYnh0U8ZMSlTo+fDys4W2Tv3cOf6+o5DeHgAPv/soMZjVUolTu5NxpTICNg6juD+XP/7NAXjx7sgMnIqdy4hvKjRE5MiXrkUtZVVUJwq4M595JFodHWp8M03R7U6Xrr7Z5iLRAhespD7c+3Zk436+mY88rsF3LmE8KJGT0zGyHFj4BsWCulPtx9f749AIMCDD0XhwIFTuHxZu4XPLp0tR2VJGcRxy7g/X0eHEt9/l4bVq2fBwcGGO58QHlo1+piYGJSVlUEmk2Hz5s2/eXzFihUoKChAXl4epFIpwsNvrOehUChQWFjY+xghhhKybDEA4OSe/RqO/K3IyKlwd3fGV18e5srL/Wkv3AMmwdXHm/tzfvHFYVhbW+Kuu2j9G2J4t13MXiAQMLlczry9vXu3A/T39+9zjK2tbe/HU6dOZaWlpb23FQoFc3Jy0tsC+hQUt4pndnzNNn7xoU65//vfE+x6wzZmZWXBlWfvNIr9qyCDxfy/x3T6vGfOfsz2H/i70b92FEM/BrTxiEQigVwuh0KhgFKpRGJiIuLi+s40aGlp6f3Y1tYWjDFNT0uIXo329sQ4v4koSOF7RQ4AlpYirFo9Czt3ZqK9vZMrt+laHc7l5iEoJpr78wLAjz8cR3R0EJycHHTKJ0QbGhu9m5sbKisre29XVVXBzc3tN8etXLkSpaWl2Lt3Lx599MbqfowxpKSkIDc3F/Hx8bf8PPHx8ZBKpZBKpXB2duY9DzLMBcVEQ61Wo/DgUe7c2NgQODraIfH7NJ0+d37yYbhO8MIYXx/u3O3bj8PcXIhVq2bq9LkJ0YbGRt/flX/9vWLfvXs3/P39sXLlSrzyyiu994eHhyMkJASxsbHYuHEjIiIi+v08W7duhVgshlgsRm1tLc85EIKgRVFQ5BWgsYb/Z+fetXNx9ep1HD7MP1MHAIoOHYVapcJ0HV7VFxQocPbsRdx9zxydPjch2tDY6KuqquDh4dF7293dHdXV1bc8Pj09HT4+PnBy6t5y7dKlSwCAmpoa7Nq1CxKJZKA1E9KH6wQvjPX1QUFyKneupaUIS5aEYPeuLKhUfMsl9Giuq4dcegpBi6J0yv9h+3FERU2DszMN3xDD0NjopVIpfH194eXlBZFIhLVr12LPnr4Xo/j43PiTNTg4GBYWFrh27RpsbGxgZ2cHALCxscGiRYtQXFys51Mgw13Qoqhfhm2OcOcuWDAd9vY22LUrc0A1FCQfxmhvT4z1m8id+8MPxyEUCrF69ewB1UDIrWhs9CqVCps2bUJycjJKS0uxfft2lJSUICEhAQkJCQCANWvWoLi4GHl5efjggw9w7733AgBcXV1x/Phx5OfnIycnB3v37kVycrJhz4gMO0Ex0Th3Mh9Ntde4c1etmomGhhakpt5+yQNNig6nQdXVpdPwTWHheZw5U0XDN8SgjD4t6Oag6ZUU2oarjzfbUpTJwteu4c4VCATsytVv2Dff/lEvtWz4+G323M/bdcp9+eX7mbJrNxs92tHoX1OKoRkDml5JyGAWtDBS52Gb8HB/uLiMwO4BDtv0KEhJhYunh47DNxkQCmn2DTEMavRkSAuIjEBFfhGarvW/dvztrFo1C+3tnThw4JReaik5ehwAMGU+/xBMUdF5lJdfwrLlNFmB6B81ejJkOYx2gUfAZJQcy9ApP27lTBw6VIDm5ja91NN0rQ4VhacRMFe3JQ32/ixFdHQQbGws9VIPIT2o0ZMhK2Bu9yyV07+8kuYRFOQNb29XvQ3b9ChJOw6PqQGwcxrJnZuUlAMrKwtERU3Ta02EUKMnQ1bA3HBcq6rGlXIFd+6qVbOgUqmwZ0+2XmsqScuAQCBAQAT/q/pjx06jsbEVy2n4hugZNXoyJJlbWsJ3phglafyv5gFg+QoJMjJKUVvbqNe6qs/IUH/pMvzn8Td6pbILycmnsHQZ//aEhNwONXoyJPlKQmBhbYWSNP7x+bFjRyE42Af79uYaoLLuV/WTZktgbmHBnbv3ZynGjXPCjBn86+YQcivU6MmQFDB/DtpbWlCem8edu3jxDADAvn0GavTHMmBpYwMf8Qzu3H37cqFWq7GMXtUTPaJGT4akgHnhOHsiByqlkjs3dkkoqqpqUVxcYYDKAHn2SXS0tiFAh+Gb2tpGZGWdwdJlNE5P9IcaPRlyxk3yhaPraJ2mVZqbC7Fw4XQc2H/SAJV16+rshCxbqlOjB4Cfk3IgFvti7NhReq6MDFfU6MmQEzB/DtRqNUrTT3Dnzpo1GSNG2GK/ARs90H3x1KhxY3Vao/7nn7u33Fy6NFTfZZFhiho9GXImh89E1ekyNF/TbhPvX4uNDYFS2YVDh/INUNkNZRlZAIDJs8O4c4uLK3DhQg0WxfCP8RPSH2r0ZEixsrfD+KkBOHNCt/nvsUtCcfx4CZqa9HM17K00XKnBJVk5JoXzN3oAOJiShwULgiAU0n9RMnD0U0SGFF9JCITm5jo1+nHjRiEoyNug4/O/djYzB94zgiCy4l/SICUlD46OdhCLfQ1QGRluqNGTIWVS+Ey0N7egopB/A5vFi0MAwODj8z3OZGRDZGmJCSHB3LmHDuVDrVYjhoZviB5QoydDyqTZYZBl50LdpeLONfS0ypuVn8yHsqNDp+Gb+vpmSKUyLFzE/0uCkJtp1ehjYmJQVlYGmUyGzZs3/+bxFStWoKCgAHl5eZBKpQgPD9c6lxBtOXt6YJTbWJ2GbQQCAaKjpyFZT0sSa6OrowPnTuZj0izd5sSnJOchLMwPjo62eq6MDEe33bVEIBAwuVzOvL29mUgkYvn5+czf37/PMba2tr0fT506lZWWlmqd21/QDlMU/UX4urvYlqJMNsp9HHduSMhEpmZJbO3auXe05nkPrWNbijKZo+to/vMND2BqlsRWr55t9K89xeCPAe0wJZFIIJfLoVAooFQqkZiYiLi4uD7HtLS09H5sa2sLxpjWuYRoa9LsMNReqEJdVTV3bnR0EADgyJGB7Q3Lq+evDz8dpllmZ59BQ0MLYmJo+IYMjMZG7+bmhsrKyt7bVVVVcHNz+81xK1euRGlpKfbu3YtHH32UKxcA4uPjIZVKIZVK4ezszH0ixLQJzc0xUTJD52mVkVHTcPr0BVy5cl3Pld3eZfk5NFypwaTZ/MM3XV0qpKYW0jg9GTCNjd7MzOw39/W8Yv+13bt3w9/fHytXrsQrr7zClQsAW7duhVgshlgsRm1trcbCyfDiNX0qLG1sdGr0IpE5IiKmIPVwgQEq0+xMZjb8ZklgJuCf+5CSfApeXq7w8+v/BRIh2tD4k1dVVQUPD4/e2+7u7qiuvvWfzunp6fDx8YGTkxN3LiG34jc7DCplF+Q5/FMjw8L8YGNjidTUOzts0+NMRjZsRjjAY8pk7tyUlO7VOWmaJRkIjY1eKpXC19cXXl5eEIlEWLt2Lfbs2dPnGB+fG+t5BAcHw8LCAteuXdMqlxBt+M0Uo6KwGB0trdy50dFBUKvVSEsrMkBlmsmyu5dD9g3jX3pYobgCubwaUdG0vSDRncZGr1KpsGnTJiQnJ6O0tBTbt29HSUkJEhISkJCQAABYs2YNiouLkZeXhw8++AD33nvvbXMJ4WFlbwf3gEm9DZNXZNQ0nDpVjuvXWzQfbAAt9ddRfUaGiWEhOuUfSS3E/PlTaTkEMiBGnxZ0c9D0Sopfx5TICLalKJNNCJnOnWtjY8naO3ayN9542KjnsOLZP7A3pEeZuYUFd+4998xhapbEJBI/o38vKAZvDGh6JSHGNlESgs62dlQUnubODQ8PgIWFyGjj8z1kWbkQWVnCMyiQO/fIke4hp6goGr4huqFGTwY937BQKPIKdNpNKjp6GpTKLhw/btwhw3Mn86Dq6oLvTP415mtqGlBYqEDUL9cCEMKLGj0Z1OycRmKsr49Os20AIDIqCFlZZ9Da2qHnyvh0tLSi8nQpfMN020wk9XAhwsP9YWkp0nNlZDigRk8GtYmh3dMKZVn8b8Q6OtpixowJRps/fzNZdi48pvjD0taGOzc1tQDW1paYNYt/iiYh1OjJoDZxZijaGptwsewsd+7cuYEQCoVGH5/vIc8+CaG5uU7LFqelFaOrS0Xj9EQn1OjJoOYrCUX5yTyoVfzLEkdHB6G1tQPZ2WcMUBm/8/lFUHZ0wFeHaZZNTW3IzZXROD3RCTV6MmiNHDsGzuPdIc/WdXx+GtLTT6Ozs0vPlemmq7MT5/OKBjROL5H4wd7eWs+VEVNHjZ4MWhMlv4zP6/BG7OjRjggM9MSRQTJs00OWnYtxk3xhN2okd25qagHMzYWIiJhigMqIKaNGTwatiZJQNF2rw2VZOXduz1j2YBmf7yHLlgIAfMT8a9ecOFGG9vbO3iWXCdEWNXoyaPmGheo8rTIqahquX2/GqVP8vyQMqarkDNqamnUavmlv70RGRiki6Q1ZwokaPRmUXLzGY4SrywDmz0/D0aPFUKvVeq5sYNQqFc7l5mGiRPd1b6ZPnwBnZwc9V0ZMGTV6Mij1NEKZDm/EenqOho/P2EE3Pt9DlnMSLp4ecBzjyp17+JdrAubPn6rvsogJo0ZPBiXfsFDUX7qMa5VV3Lk3xucHx4VSN+tdtliH5RByc2VobGylcXrChRo9GXTMzMwwUTxjQMsSX7lSj9OnL+i5Mv24LCtH07U6nYZvVCo10tKKaZyecKFGTwadsX4TYTvSUef581FR0wbdbJubyXNO6jyf/khqIfz83ODh4aLnqoipokZPBp2eDTrkUv5GP2mSO8aNcxq04/M95DknMWK0C1y8xnPn9ozTR0bSOD3RjlaNPiYmBmVlZZDJZNi8efNvHr/vvvtQUFCAgoICZGRkYNq0G39WKhQKFBYWIi8vD1KpVH+VE5PlKwnFVUUFGq7UcOdGRw/O+fM363mTWZfhm+LiCtTWNtLwDdGaxkYvEAjwwQcfIDY2FgEBAVi3bh38/f37HKNQKDBv3jwEBQXhlVdewSeffNLn8cjISAQHB0Ms5t8zkwwvAnMhJoROH9CyxBUVV3Hu3GU9V6Zf1yqrUH/psk7DN4wxHDlSiMhIavREOxobvUQigVwuh0KhgFKpRGJiIuLi4vock5mZievXrwMAsrKy4O7ubphqicnzmOIPK1tbnd6INTMzQ2Tk1EGzLLEm8pyTmCieATMzM+7cI6mFGD/eBT4+Yw1QGTE1Ghu9m5sbKisre29XVVXBzc3tlsevX78e+/fv773NGENKSgpyc3MRHx9/y7z4+HhIpVJIpVI4OztrWz8xMT1DGeXSU9y5QUHeGDXKftAP2/SQZZ+E7UhHjPH14c7tOUdatphoQ2Oj7+/VBmOs32Pnz5+P9evX9xnHDw8PR0hICGJjY7Fx40ZERET0m7t161aIxWKIxWLU1tZqWz8xMb6SUFSfkaHlegN3bs/c8iNHhkajL8/RfZz+7NmLuHjxGo3TE61obPRVVVXw8PDove3u7o7q6urfHDd16lR8+umniIuLQ11dXe/9ly5dAgDU1NRg165dkEgk+qibmCBzCwt4BU8d0Pz5srIqVFfXaT54ELh+5SpqKirhq+NyCKmphfSKnmhFY6OXSqXw9fWFl5cXRCIR1q5diz179vQ5xsPDAzt37sSDDz4ImUzWe7+NjQ3s7Ox6P160aBGKi4v1fArEVHhOmwKRpSXkOfzDNubmQsydO2XIjM/3kOecxITQYAiEQu7cI6kFGD3aEVOm8E/RJMOLxkavUqmwadMmJCcno7S0FNu3b0dJSQkSEhKQkJAAAHjxxRfh5OSE//73v32mUbq6uuL48ePIz89HTk4O9u7di+TkZMOeERmyJoaFdi/6dTKPO1cs9oWdnfWQGZ/vIc/OhbW9Hdz8J3Hn9pwrzb4h2mCDLaRSqdFroLjzsfGLD9kfvvufTrl/+cu9TM2S2KhR9kY/D56wcxrJthRlsqj1D+qULy/fynbsfN7o50Fh/Lhd36QrY8mgYGFtBc9pgZDn6D4+n5dXjrq6Jj1XZljN1+pxSVY+oGWL58+fCoGA/iuTW6OfDjIoeAcHQSgy12lZYisrC8ye7Y/Uw0Nr2KaHLDu3+/zNzblzU1MLMXKkHaZP9zZAZcRUUKMng8JEyQx0KZU4n8/frMPD/WFpKRq0yxJrUi49BQtrK4yfxr8XbM9UUpp9Q26HGj0ZFCZKQnGh8DQ629q5c6OipqGrS4X09BIDVGZ45bl5UKvVOk2zvHy5HiUlFxAZRevTk1ujRk+MzsrOFu4BkyDX4WpYoHt8PifnLJqb2/Rc2Z3R1tiEi6VnMHEAyxZHRATA3Jx/iiYZHqjRE6ObENI9j1yuw4VSDg42EIt9h9z8+ZvJs0/CMygQIitL7tzU1ELY2VlDLPY1QGXEFFCjJ0Y3MSwEyo4OVBSe5s6dO3cKhELhkJs/fzNZzkmYi0Twms4/1p6W1r0JOo3Tk1uhRk+MbqJ4Bs7nFaGrs5M7NyoqCG1tHcjMLDNAZXeO4lQBVMounaZZ1tU1oaBAQevekFuiRk+MytZxBNwm+0E2gPnzGRml6OhQ6rmyO6uzrQ0Xikt0XvfmSGohZs/2h5WVhZ4rI6aAGj0xKh/xDABAuQ7r2zg7OyAoyHvQbxuoLXnOSXgE+sPKzpY7NzW1EFZWFpg1a7IBKiNDHTV6YlQTJSHoaG3FhdP8UyN71ngZ6uPzPWTZuRAIhfCeMZ07Nz39NLq6VDROT/pFjZ4Y1URJCM6dKoC6S8WdGxU1DY2NrcjNlWk+eAioKCiGsqMDvmH8wzdNTW2QSmWYTxuGk35QoydGY+/sBNcJXpDrsOwB0D0+n5ZWDJVKrefKjKOrsxPn84owUaz7OL1E4gc7O2s9V0aGOmr0xGh6ZpjospCZu7sz/PzccHSI7CalLVlOLtz8/WAzwoE798iRQohE5pgzJ8AAlZGhjBo9MRpfSQhaGxtxsYx/6KVn28BDh4b2hVI363lTuudNah4nTpSho0NJ4/TkN6jRE6PxkczAudw8MDX/0Ev0giBcuVKP4uIKA1RmPBdOl6CjtRW+OiyH0HM9Ac2nJzfTqtHHxMSgrKwMMpmsz8bfPe677z4UFBSgoKAAGRkZmDZtmta5ZHgaOXYMnD3cddo2EOh+RX/4cOEtN6ofqtRdKpw7mT+g9emDgydg5Eg7PVdGhrrb7loiEAiYXC5n3t7eTCQSsfz8fObv79/nmFmzZjFHR0cGgC1evJhlZWVpndtf0A5Tph/iuCVsS1EmG+Prw50bEDCeqVkSe/TRhUY/D0PE/IfvY1uKMpmDizN3bnh4AFOzJBYXN9Po50FxZ2NAO0xJJBLI5XIoFAoolUokJiYiLi6uzzGZmZm4fv06ACArKwvu7u5a55LhaaIkFM119bgiP8edu2BBz/h8vr7LGhR6rhLWZZw+J+csWlraaZye9KGx0bu5uaGysrL3dlVVFdzc3G55/Pr167F//37u3Pj4eEilUkilUjg7O2t9AmRomhgWgvLcPJ2GXqIXTIdMVo0LF2oMUJnxVZ+Ro7WhUaflEJTKLqSnn6ZxetKHxkZvZmb2m/tu9Z9z/vz5WL9+fe9YPE/u1q1bIRaLIRaLUVtbq6ksMoSN9vaEo+tonM3M4c41Nxdi/vzAIb8s8e0wtRrluXmYqMOFU0D3OH1goCdGj3bUc2VkqNLY6KuqquDh4dF7293dHdXV1b85burUqfj0008RFxeHuro6rlwyvPjNkgCATo1eLPaFvb2NyQ7b9JDn5MLJ3Q2j3MZy5x45UgQAiKSrZMkvNDZ6qVQKX19feHl5QSQSYe3atdizZ0+fYzw8PLBz5048+OCDkMlkXLlk+PGbJUHthSrUXbzEnbtgwXSo1ereZmaqejZJ1+Uq2by8cly/3kzj9KSXxm3nVSoVNm3ahOTkZAiFQnz22WcoKSlBQkICAODjjz/Giy++CCcnJ/z3v/8FAHR1dUEsFt8ylwxfQnNzTJTMQO6e/TrlRy8IQl7eOdTVNem5ssHlSrkCTdfqMDEsBDm7f+bKVanUSEsrxvxIavTkBqNPC7o5aHql6caEkOlsS1EmC4yax51ra2vF2jt2sjfeeNjo53En4oE3/85ePLxHp9wnnljO1CyJeXi4GP08KO5MDGh6JSH65DdLArVKBbmUfyGziIgpsLAQmdyyB7ciy87FiNEucJ3gxZ3bs3RzdDS9qie0BAK5w/xminGhuATtTc3cuQsWBKGjQ4mMjFIDVDb4nM2UAgD8Zodx5xYXV+DSpTosWBis77LIEESNntwx1g728Aj0721gvKKig5CRUYK2tg49VzY41V+6jKuKCvjNEuuUn5KSh4ULp/c7zZkML9ToyR0zUTwDAqEQZ09kc+e6uIzA9OkTcHiYDNv0OJuZA5/QGRCKRNy5B1Pyer9uZHijRk/uGL9ZErQ3t6Ci6DR3bs9UQVOfP3+zMydyYGljDa/p/HPie97LWLSIhm+GO2r05I7xmy1BufSUTtsGLlgwHdevN+PkyXIDVDZ4lUtPQaXswiQdxumvXr/y+X4AACAASURBVL2OvLxyLFzEvwctMS3U6MkdMcp9HJw93HE2i/9qWABYFBOM1NRCqHVYu34o62htxfnCIvjNluiUfzAlD3PmBMDGxlLPlZGhhBo9uSN6lj04c4K/0fv7e8DDwwUH9uu2t+xQd/ZEDjwCJsN2JP/aNSkpebCwEGHevEADVEaGCmr05I7wmylG/aXLqDl/gTt38eLu5XqTk/P0XdaQ0PPLUZddpzIyStHa2kHj9MMcNXpicGYCAXzDQnWeVhmzOASnT19AZaVpLkusSVVJGVobGnUap+/oUOLYsWIspEY/rFGjJwbnHjAZNiMcdFqt0sbGEvPmBSIlWbctB00BU6txNks6oHH6gIDxcHenfR6GK2r0xOAm/dKgZNm53Lnz5gXC0lKEAweGb6MHgLMnsuHoOhpjfH24c1NSuoe8YmL4d6wipoEaPTE4/4jZuFBUgpb669y5MTEz0NragWPHig1Q2dBRejwLAOA/ZyZ3bs+w1+JY3TYyIUMfNXpiULaOIzB+2hSUHsvQKT9m8QwcPVqEjg6lnisbWhqv1uBi2VlMjpitU/6+vblYuHA6RCKNK5MTE0SNnhiU3+wwCAQClKZncud6e7ti0iR3JA/zYZsepemZ8A6eBis7W+7cffty4eBgg/BwfwNURgY7rRp9TEwMysrKIJPJeveD/bVJkybhxIkTaG9vxzPPPNPnMYVCgcLCQuTl5UEq1W3WBRm6AubORtO1OlSVlHHn9owpHzgwPOfP36ws/QSE5ua91yTwSE0tREeHEkuX6rZAGhnaNDZ6gUCADz74ALGxsQgICMC6devg79/3VUFdXR2eeOIJvPXWW/0+R2RkJIKDgyEW0w/ZcGImEGBS+EyUHc+65abwtxOzeAYUiiuQyWifYQCoKDyN1sZG+OswfNPS0o60tGLELqFx+uFIY6OXSCSQy+VQKBRQKpVITExEXFxcn2NqamqQm5sLpXJ4j6OSvsZPDYCt4wiUpZ/gzrWyssDChcHYt5f+CuyhVqlwJiMbk+fM1Gnp4f37chEQMB5eXq4GqI4MZhobvZubGyorK3tvV1VVwc3NTetPwBhDSkoKcnNzER8fr1uVZEjyj5gNVVcXzugwfz4ycipsbCzx88/U6H+tND0TDi7OGDfZlzt3377u6a2xNPtm2NHY6Pt75cDzZ3h4eDhCQkIQGxuLjRs3IiIiot/j4uPjIZVKIZVK4exMF3aYAv+I2ThfUIS2Rv6NvJcvl6C5uQ1HjxYZoLKh60zGL9MsdRi+kcmqIZNVY8lS/qUUyNCmsdFXVVXBw8Oj97a7uzuqq7UfM7106RKA7uGdXbt2QSLp/42krVu3QiwWQywWo7a2VuvnJ4OTg4sz3AMmoUyH2TYAsHSZGCkpecN+WuXNmuvqcaGoRKdGD3QP30RFTYOVlYWeKyODmcZGL5VK4evrCy8vL4hEIqxduxZ79uzR6sltbGxgZ2fX+/GiRYtQXDy8L3wZLiaHd1/YU6rD+HxQkDc8PFywl4Zt+lWafgLjp02BreMI7tx9+3JhbW2JyEj+jUzI0KWx0atUKmzatAnJyckoLS3F9u3bUVJSgoSEBCQkJAAAXF1dUVlZiaeffhp/+ctfUFlZCXt7e7i6uuL48ePIz89HTk4O9u7di+TkZIOfFDG+KZFzcP3yFVw6y79RyPLlEqjVauzdy79kwnBQkpYBgUAA/7nh3LlpacVobm7D8uW6rZtDhi422EIqlRq9BgrdQ2RlyV7POcJW/flpnfKzsrewE5n/Mvp5DOb466Gf2CNvv6FT7vYfnmOVVV8wMzMzo58Hhf7idn2Trowleuc3SwILaysUpx7jznV1dYRE4kfDNhoUpx7DpNlhEFnx7xz10+4suLk5ITR0ogEqI4MRNXqid4FRc9Ha2Ijyk/wbhfRcuZmUpNuWg8NFceoxWFhbwW8m/0WIe/dK0dWlwsqV/AukkaGJGj3RK4FQiCnz5qD02AmdNgFftlyCCxdqUFh4Xv/FmZDy3FNoa2xCYNQ87tzr11uQllaMOGr0wwY1eqJXXsHTYDvSUadhG1tbK8TEBGPPT1kGqMy0qLtUKE0/gYB54TAT8P83/ml3FgICxsPXd5wBqiODDTV6oleBkRFQdnTgTEY2d25sbAisrS2xYwf/lMzhqCj1GOxGjYTXdP6pkj/91P39iYvj356QDD3U6IleBUbNhSwrFx2trdy5q9fMxtWr15GeXmKAykzPmeNZ6OrsRGDUXO7cysoanDwpp+GbYYIaPdGbsX4+cHJ3Q3FqGneulZUFli4NxU+7s6BWqw1QnenpaG2FLDtXp0YPdA/fzJo1Ga6ujnqujAw21OiJ3gRGzYNarcbptOPcuYsWBcPe3gY//kjDNjyKU4/B2cMdY/3495L96adsCAQCrFhBwzemjho90ZugRVE4n1eI5mv13Lmr18xGXV0TjhwpNEBlpqs49RjUKhWCYqK5c4uKzkMmq8Zdd/NfYUuGFmr0RC9cfbwx1tcH+QcOceeKROZYsUKCPXty0KXDlMzhrLmuHvKck5ges0Cn/O3b0hEVNQ0uLvzr5pChgxo90YvpixdArVKh4GAqd25U1DQ4Otph5w7dNhAf7vIPHIKLpwfc/P24c7dtS4dQKMSaNbqthkmGBmr0RC+CFy9AuTRPp2GbNWtmo6mpFQcP5hugMtNXeCgNKmUXpi/mf1VfXFyBkpILuHetbm/okqGBGj0ZMLfJfnDxGo+8Awe5c0Uic6xaPQtJSVJae15HbY2NOJuVg6BF/OP0ALAtMR0REQEYO3aUnisjgwU1ejJg02MXQKXsQtGho9y5MTHBcHJywHff8ueSG/IPHIaT+ziMnxrAnbttWzoEAgHupjdlTRY1ejJgQYuicTZbitaGRu7c+x+IRE1NA1JS+BdAIzcUp6ahq7NTp+Gbs2cvIj//HA3fmDBq9GRAxk8NgJP7OBToMNvG3t4aK1ZI8MP24zTbZoDam1tQlpGFoJjofvd51mT7tnTMmjUZ48e7GKA6YmxaNfqYmBiUlZVBJpNh8+bNv3l80qRJOHHiBNrb2/HMM89w5ZKhbXrsQnR1dqJIh0XMVq+eDWtrS3zzzREDVDb85B84DEfX0fCeEcSdu21bOgDg3nsj9F0WGQQ0NnqBQIAPPvgAsbGxCAgIwLp16+Dv79/nmLq6OjzxxBN46623uHPJ0CUwF2LGkkUoOXYC7U3N3Pn33T8P5eWXkJV1xgDVDT+njxxDe0sLQlcs4c5VKK4gM7MMDz4UZYDKiLFpbPQSiQRyuRwKhQJKpRKJiYmIi4vrc0xNTQ1yc3OhVCq5c8nQ5T9nFuydRkG6ey937tixoxAdHYTvvuVfF4f0r7OtHQXJqQiKiYKFtRV3/hefH0JgoCdCQ30NUB0xJo2N3s3NDZWVlb23q6qq4ObmptWT8+TGx8dDKpVCKpXC2dlZq+cnxiVeuQxN1+pQlpHJnbt2bQQEAgG+pdk2eiX9aS+sbG0xdUEkd+62beloa+vA736n21W2ZPDS2Oj7e2OHMabVk/Pkbt26FWKxGGKxGLW1tVo9PzEe25GOCJgbjpNJB3TaSeqBByMhlcpw9uxFA1Q3fClOFaD2QhXEcfzDN42Nrdi5MxNr182FpaXIANURY9HY6KuqquDh4dF7293dHdXV1Vo9+UByyeA2Y2kMhCJz5Oz+mT93hg+Cg33w9Vf8yyUQzaQ/7YVvWChGuY3lzv3i80MYOdKONiQxMRobvVQqha+vL7y8vCASibB27Vrs2bNHqycfSC4Z3CQrl+JCUQmulCu4c+PjY9DW1oGvv6bZNoaQu2c/1Go1QpfHcuemphbiwoUaPPwIDd+YGqYpYmNj2ZkzZ5hcLmfPP/88A8ASEhJYQkICA8BcXV1ZZWUla2hoYPX19ayyspLZ29vfMldTSKVSrY6jME64TfZjW4oy2ex7V3Pn2tpasYbGbeyzz580+nmYciR88g57fv8OZmZmxp378sv3M2XXbjZu3CijnweF9qGhbxq/QM6CKYwcK597ir2Re5RZO9hz565fv4ipWRKbNWuy0c/DlGPG0kVsS1EmmygJ4c6dMGEMU7Mk9uc/323086DQPm7XN+nKWMLFwtoaoSuWoOhwGtoam7jz4zfEoLi4ApmZZQaojvQoPJSGlusNmH3vau7cc+cuIzW1APEbYiAQUIswBfRdJFxCli2Gtb0djn/3A3fu9OkTIJH4YesnyQaojPxaV0cHsnfuQWDUXDi6jubO/+D9vfDycsWyZWIDVEfuNGr0hEv4ujWoLClDRUExd27Pm7C05MGdcWLbTpgJBJh1zyru3D17snHhQg02Pb7MAJWRO40aPdHaREkIxvr6IEOHV/MODja4/4H5+OGHDNTX8y+XQPjVV19GSdpxzLwrDuYWFly5KpUaH324DwsWTIe/v4fmBDKoUaMnWgtfdxda6q8j78Bh7tzHHlsEBwcbvPdukgEqI7dy/LsfYTdqJIIW8a9h8+mnKWhv78TGjUsNUBm5k6jRE62MHDsGgZERyNqxB10dHVy55uZCPPGHFTh6tAgnT8oNVCHpjyxLiquKCoSvu4s7t7a2EYmJ6Xjo4Sg4ONgYoDpyp1CjJ1qZfW/3OO+JbTu5c+++ew7Gj3fBv7fs0ndZRAvHv/8RntOm6LT71PvvJcHOzhqPPKLbNoVkcKBGTzSytLXBzLtXojj1GK5fvsKd//QzK1FWVoW9e3MNUB3RJPenfWhtbETkow9y5546VY6MjBI8+VQczM2FBqiO3AnU6IlG4WvXwMbBAYc//ZI7d968QISETMR//r1b68XwiH51tLbi+Hc/YtqC+XD18ebOf/21H+Dl5Yr775+v/+LIHUGNntyWhbUV5j20DqXHM1FVwr9ByDN/XI2rV6/TujZGlv7NNnS0tiL6sYe4c/fty0VeXjme+/NddAHVEEXfNXJbYWviYDdqJA59/AV3blCQN5YtE+O/H+xFe3un/osjWmttaETm9t0Ijl0IJ3ft9pP4tdde3Y5Jk9xx112zDVAdMTRq9OSWhCIRIh+5H/KckzifX8id//IrD6C+vhnvvENTKgeDo19+B1VXFyLXP8Cdu3NnJkpKLuD5F+7RafNxYlzU6MktiVcuxQhXFxzayj82HxY2CcuXS/DWv3aioaHFANURXk2115Cz62eI47q/rzwYY3jj9R8xbZo3LYswBFGjJ/0yt7BA9GMPoaKgGLIsKXf+y6/cj5qaBrxLF0gNKkc+/wZmMMPChEe5c7//Pg3l5Zfwt7/fR6/qhxhq9KRfEQ/cg1HjxmLfux9x586bF4iFC4Pxxus/oKWl3QDVEV3VV19GxrYdCFu9HK4TvLhyVSo1XvzrtwgO9sEDD8w3SH3EMKjRk9+wHemI6McexumjxyHPOcmd//IrD6C6+ho+/HC/AaojA3Xo48/R0dKKZU9v4s5NTDyGnJyzePW1h2BtbWmA6oghaNXoY2JiUFZWBplMhs2bN/d7zDvvvAOZTIaCggIEBwf33q9QKFBYWIi8vDxIpfxDAOTOW/T79bCwtsLP/36fO3f16tmIiJiCf7yyjWbaDFIt1xtwaOuXCJgXDt+wUK5cxhj+9MfP4O7ujKeeijNQhcQQbrtriUAgYHK5nHl7ezORSMTy8/OZv79/n2NiY2PZvn37GAAWFhbGsrKyeh9TKBTMyclJbzulUBg2Rnt7sn/mpbPVL/yRO9fGxpKdr/iM5eW/y4RCgdHPheLWYW5hwV44sJM9vf1LZibg/17t2Pk8a2jcxlxdHY1+LhTdMaAdpiQSCeRyORQKBZRKJRITExEX1/c3eVxcHL766isAQHZ2NhwdHTFmzBhNT00GoWVPb4KyvQMpH/6PO/eFF+7B+PEu2LTxI6hUagNUR/Slq7MT+975EG7+fhCvWMKd/9zmL2BlZYGXX77fANURfdPY6N3c3FBZWdl7u6qqCm5ublofwxhDSkoKcnNzER8ff8vPEx8fD6lUCqlUCmdnZ+4TIQM3dcF8TJk/Bwc//hzNdfVcub6+4/DMH1fhyy8PIyOjxEAVEn3KP3AIirxCLP/j47BzGsmVK5NV4/33fkb8hsWYM4d/sTRyZ2ls9P1No7p5zZLbHRMeHo6QkBDExsZi48aNiIiI6PfzbN26FWKxGGKxGLW1tVoVT/TH2sEBq1/4IypLynDs60Tu/HffS0BbWyc2P/uF/osjBsEYw/aXXoOFtRVWPfc0d/6LL34LheIKPv3fE7Cy4tvYhNxZGht9VVUVPDxu7DDj7u6O6upqrY+5dOkSAKCmpga7du2CRCLRS+FEv+KefQK2I0Zg219fhVql4sp9+OFoxMTMwIt//RZXr143UIXEEK4qKpDy0WeYvngBAqPmcuW2tLRjQ/x78PNzw4svrjVQhURfbjvALxQKWXl5OfPy8up9MzYgIKDPMUuWLOnzZmx2djYDwGxsbJidnV3vxxkZGSwmJmZAbypQ6D8mhc9kW4oy2eLHN3Dnenu7sobGbSz1yGtMoMObehTGD4G5kD3z41fsxcN7mJW9HXf+1q2Ps07lbhYc7GP0cxnOoaFvan6C2NhYdubMGSaXy9nzzz/PALCEhASWkJDQe8z777/P5HI5KywsZCEhIQwA8/b2Zvn5+Sw/P58VFxf35g6wYAo9hs0IB/aXlF3s2Z++Z+YWFly5QqGApR9/k9XVf888PFyMfi4Uuod7wCT2r/zj7L7XX+LOHTHCllVd/IIVFL7HrK0tjX4uwzUG3OgHWcEUegozMzP22H+3sDdPHWMeU/y581944R6mZkls3bp5Rj8XioHHgg2PsC1FmWzm3Su5cxctCmZdqp/Y//73hNHPY7gGNXqKfqPnP/bse1dz50ZETGGdyt3sm2/559tTDM7o/cV/Mo25B0zizv/b3+5japbEHn10odHPZTgGNXqK34TvTDH7V0EGu/+Nv3HnTpgwhl2t+ZaVlH7IRoywNfq5UOgveobynt+/g1k72HPlCgQCdiD5ZdbatoNNnz7B6Ocy3IIaPUWfcPb0YH9P28f+tOtbZmFtxZXr4GDDik9/wGqvfcd8fMYa/Vwo9B/jpwawN08dYwmfvMOEIhFXrrOzA7tQ+TkrP/cpGzNmpNHPZTjFgK6MJabFYbQLEj55B4wxfP6Hzehs0351SaFQgO8Tn4Wv7zjcteZ1lJdfMmClxFguFJXgh7+9Ab9ZEtz32osw49g+sLa2EatWvgoXFwccSP47HB1tDVgp4WH030Q3B72iN0xYOziwP+36lr2aeYh7DFYoFLBvv/sjU7MkFh+veYosxdCP+Q/fx7YUZeq07lF0dBBra9/J0o+/STNx7lDQ0A0Fs7K3Y098s5W9kXuU+YQGc+UKhQL2feKzTM2S2J/+xP/GLcXQjWVPbWRbijLZ0id/z527Zs1s1qX6ie3d9zdq9ncgqNEP8xjh6sL+uPMb9uapY2xKZARXrrm5kCVu28zULIk988wqo58LxZ2PNX99lm0pymT3/O3PTCAUcuWuX7+Idal+YunH32QjR/JfjEWhfVCjH8bh6uPN/npwN3s18xCbKAnhynV2dmCpR15japbEnnoqzujnQmG8WLxpA9tSlMkeffefTGTF9+p8zZrZrK19Jys+/QFdWGfAoEY/TGNKZAR7JSOZvZSaxMZN8uXKnTbNi51TfMpa23aw+++fb/RzoTB+zL53NftXQQZ7MvEz5uTuxpU7b14gq7+eyKoufsHmzg00+rmYYlCjH2ZhbmHBVv35abalKJM9ue1zNnLcGK78Rx9dyJqaf2CVVV+w0FC+XxAUph09Lx5ezTzEpsdEc+UGBnqysjMfMWXXbvbSS+tocxo9BzX6YRQegQHs6R++ZFuKMtmKZ//ANQ/aw8OF7T/wd6ZmSexw6qs0D5qi3xg5bgx7/OtP2JaiTLb2H39ldqO0/zmxtbVin33+JFOzJJZ27A3m7+9h9PMxlaBGPwzCzmkku+fvz7MtRZnspdQk5j83XOtca2tL9qc/rWbXG7axhsZt7P/+L5aZmZkZ/ZwoBm8IzIVs8eMb2D9PpbN/ZKSwOffdxfVG7f33z2fX6r5nncrd7N//fow5ONgY/ZyGelCjN+GwHenIFm/awP5x4iD756l0tuzpTczSVrv/NCKROduwYTGruvgFU7MktifpRebpOdro50QxdGK0tyfb8PHbbEtRJnv2p+9ZyPJYJjDXruE7OTmwjz7ayLpUP7HLV75mmzffRQ1/AEGN3gTD1cebrfrz0+z1nCPsXwUZ7OF/v8ZGe3tqlTt27Cj20kvr2MXqL5maJbFj6W+yOXMCDFovhWlHYNQ89syOr9mWokz2woGdLOKBe5mt4witcmfM8OkdMqy/nshef/1hNmEC3/tKFLfvm2a/fDCoSKVSiMViY5cx6Di6jkZg9FyErlgCjyn+UCm7cPLnAzjy+Te4qqi4fa6jLVasCMOau8KxePEMCIUCHDhwCu+9m4Tk5FN36AyIqfOfG47oxx6Cd/A0dCmVKEnLwKmfD+DMiRx0trXdNjc42AfPbl6Du+8Oh0AgQHr6aXz9VSr27Mmhncu0cLu+SY1+ELN2cIBn0BT4hAZj8pxZGOc3EQBQVXIGuXv2IW//wVtu4m1hYQ6x2A9RUdMwP3Iq5swJgEhkjoqKq9i+LR0ff3wA585dvpOnQ4aRsX4+CF2xBCHLFsPeaRS6Ojtx7mQ+yo5nQZFfiIulZ6FSKvvNdXNzwgMPROKhh6Pg79+9RalUKsOB/Sdx7FgxcnLOoqnp9r80hqMBN/qYmBi88847EAqF+PTTT/Hmm2/+5ph33nkHS5YsQWtrKx555BHk5eVpnctTsCkSWVnCxXM8Rnt7wnWCF0Z7e2LcJF+M9vYEAKiUXTh3Kh9l6ZkoOZbR59W7lZUFJkwYAx+fMfDxGYvAwPEInuGDwEBPiETmUKvVyM9X4NDBPPz44wnk5sqMdZpkGBKYC+EdHAT/iNnwj5iFMRMnAACUHR2oLpPhcrkCV8+dxxVFBa6eO4/6S5f77FkcFOSNJUtCEbskFLNmTYJQKIRarcbp0xdQVFSBstJKlJVVobS0EnL5JXR09P/LYzgYUKMXCAQ4e/YsFi5ciKqqKkilUqxbtw6lpaW9x8TGxuLxxx/HkiVLEBYWhnfeeQczZ87UKpe34MHEzMwMQpEI5hYimFtYwFwkgtDC4pfbIlhYW8Pa3h5W9rbd/9p1/2s7cgQcR7vA0dUZo1ydYTfCHkIzQGDGYKZWoa22Bi1XL6HxYhVarlyCsqEO1pbmcHZ2gKurI1xGO8LV1RFjx46Eu7tzn5pqahpw8qQceafKIZXKkJZWjPr6ZiN9hQjpy97ZCV5BgfAMmgqPQH+M9vaEg7NT7+NqlQpN1+rQWFOLhqs1aKy5htaGRrQ3N0PQ1QFfNzv4+zgjcJIrJnq5wMN9VJ/nr6lpQHV1XXdcvIZr1xrR0NCK69db0NDQ8su/rWhpaUd7eyc6OpRob1f+8m/3bZVKfae/LHpxu75prilZIpFALpdDoVAAABITExEXF9enWcfFxeGrr74CAGRnZ8PR0RFjxoyBl5eXxlx9qri2C5ai7iVVzcxu3N/9sRl++Qe9D/V+bNb3+D55/Xz8S6LZrw7+1cO//fhXdwjAIDRjEAp67lQCqLvpTCwBeP0SNyiVXbh6tQFXrlzHlSvXUVxcAcW5yygvvwy5vBrl5ZdRV9cEQgarptprKDqchqLDab33WTs4wNXbE6MneGHkuDEY4eIMB1cXjHIbB6+gqbC2t4dQdKNVdQHIB5B/HTBvYBhpqcIoyy44Wqhha24FO28XBPuqEGGuhpU5gznnYuxqxqBmAOsJ4MZtdP+rZgADu3HMrx67lV8/1Pe4GzcaW7vgN+YuvoK1oLHRu7m5obKysvd2VVUVwsLCNB7j5uamVW6P+Ph4bNiwAQDg7Ozc7zGaVFRfh4VI2PtFZOyXbwS6vwvslwd67gdjvW9L41f3qX85pienNx8AU3c/h1qtBlOpoFKpoVarwFRqqFQqqFVqqNVqqFUqqJRd6OpUQtnZCZVSCWVHJzraO9HZ2YXOTiU6O7ugVKp+c7u5uQ1NTW1obGxFU1NbbzQ2tur0dSFkMGtrbMT5giKcLyi65TEW1lawsrODtb1d71/JIkvL7r+kLUQwt7Ds/tfSovevazOBAGYCM1iIzGFnawE7axFsbUSws7GAlaU5LERCiMwFsLQwh8hc0HvbXGgGgaD7xZ9QYAYzs+6PBTCDQIAbtwXdLx4FZjfu72F28wmYmfV5gdn3oRt3NDVrvz8ED42N/tdF9GA3/dq61THa5PbYunUrtm7dCqD7TxBdzJ36O53yCCGDW2dbOzrb2tFYU2vsUoYkjY2+qqoKHh4evbfd3d1RXV2t1TEWFhYacwkhhBjebSfhC4VCVl5ezry8vJhIJGL5+fksIKDvxTVLlixh+/btYwBYWFgYy87O1jq3v6ALpigoKCj44nZ9U+MrepVKhU2bNiE5ORlCoRCfffYZSkpKkJCQAAD4+OOPsW/fPixZsgRyuRytra343e9+d9tcQgghdw5dMEUIISbgdn2Tc+IRIYSQoYYaPSGEmDhq9IQQYuKo0RNCiIkblG/GXr16FRUVt192d7BxdnZGbe3wupiDznl4oHMeGjw9PTF69OhbPm70+Z+mEMNx7j+d8/AIOuehHzR0QwghJo4aPSGEmDghgL8ZuwhTcerU8NuSj855eKBzHtoG5ZuxhBBC9IeGbgghxMRRoyeEEBNHjd4AnnnmGTDG4OTkpPngIe6f//wnSktLUVBQgJ07d2LEiBHGLskgYmJiUFZWBplMhs2bNxu7HINzd3dHamoqSkpKUFxcjCeeeMLYJd0xAoEAp06dQlJSkrFL0Sujz/E0pXB3d2cHDhxg58+fZ05OTkavx9CxcOFCJhQKGQD2xhtvsDfeeMPoNek7BAIBk8vlCxQ6fAAAAk1JREFUzNvbu3dfBX9/f6PXZcgYM2YMCw4OZgCYnZ0dO3PmjMmfc0889dRT7Ntvv2VJSUlGr0VfQa/o9ew///kPnn322VtumWhqDh48CJVKBQDIysqCu7u7kSvSP4lE0rvJvVKp7N3k3pRdvnwZeXl5AIDm5maUlpbCzc3NyFUZnpubG5YuXYpPP/3U2KXoFTV6PVq+fDkuXryIwsJCY5diFI8++ij2799v7DL0rr9N7odD0+vh6emJ4OBgZGdnG7sUg3v77bfx7LPPQq1WG7sUvdK4wxTp6+DBgxgzZsxv7n/hhRfw/PPPY9GiRUaoyrBud8579uwBADz//PPo6urCt99+e6fLMzieTe5Nja2tLXbs2IEnn3wSTU1Nxi7HoJYuXYqrV6/i1KlTmDdvnrHL0Tujjx+ZQgQGBrIrV64whULBFAoFUyqVrKKigrm6uhq9NkPHQw89xE6cOMGsra2NXoshYubMmezAgQO9t5977jn23HPPGb0uQ4e5uTk7cOAAe+qpp4xey52I1157jVVWVjKFQsEuXbrEWlpa2Ndff230uvQURi/AJEOhUAyLN2NjYmLY6dOnmbOzs9FrMVTousn9UI8vv/yS/ec//zF6HcaIefPm0ZuxhPR4//33YW9vj4MHDyIvLw8ffvihsUvSu19vcl9aWort27eb/Cb34eHheOihhxAVFYW8vDzk5eUhNjbW2GURHdESCIQQYuLoFT0hhJg4avSEEGLiqNETQoiJo0ZPCCEmjho9IYSYOGr0hBBi4qjRE0KIifv/XtgSwnAmKTwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effect size d = 0.5\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import math\n",
    "\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "x = np.linspace(-5, 5, 200)\n",
    "mu1 = -.3\n",
    "sigma = 1\n",
    "plt.plot(x, stats.norm.pdf(x, mu1, sigma))\n",
    "mu2 = .2\n",
    "sigma = 1\n",
    "plt.plot(x, stats.norm.pdf(x, mu2, sigma))\n",
    "plt.show()\n",
    "print(\"Effect size d =\", abs((mu1-mu2)/sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Let's play - random samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAQx0lEQVR4nO3cX0zT97/H8RcMbn5ygYMJS8HqRdlGzu+icZURsuCN05oYTMYFMwvGJUDM1PiLFyQuxHPpLn4xLjJDGhYlw3jDYJjBpsZdELOR76DlTygKCTE0WJwYDKLLxPM9F+c3IqPQIi0tn/N8JCS03w/t2w/sueYL36ZJsgUA2PTSkz0AACA+CDoAGIKgA4AhCDoAGIKgA4AhMpL1xA8fPtT9+/eT9fQAsCk5nU5t27Yt4rGkBf3+/fvyeDzJenoA2JQsy1rxGKdcAMAQBB0ADEHQAcAQBB0ADEHQAcAQBB0ADEHQAcAQBB0ADEHQAcAQSbtSdD3+PfTLisdO/7N0AycBgNTBK3QAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMETUoBcUFOj27dsaGRnR8PCwTp48uWxNeXm5Zmdn5ff75ff71dDQkJBhAQAri/rmXAsLCzp9+rT8fr+ysrLU19enmzdvKhgMLlnX09OjgwcPJmxQAMDqor5CD4fD8vv9kqSnT58qGAzK4XAkfDAAwNqs6Ry60+mU2+1Wb2/vsmOlpaUKBALq6upScXFxxK+vqamRZVmyLEu5ubmvNzEAIKKY3w99y5Ytamtr06lTpzQ3N7fkWH9/v5xOp+bn5+X1etXR0aGioqJlj+Hz+eTz+SRJlmWtc3QAwKtieoWekZGhtrY2tba2qr29fdnxubk5zc/PS5K6u7uVmZmpnJyc+E4KAFhVTEFvbm5WMBjU+fPnIx7Py8tb/Nzj8Sg9PV0zMzPxmRAAEJOop1zKyspUXV2twcHBxV+OnjlzRtu3b5ckNTU1qbKyUseOHdPCwoKeP3+uqqqqxE4NAFgmatDv3LmjtLS0Vdc0NjaqsbExbkMBANaOK0UBwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBAEHQAMQdABwBBRg15QUKDbt29rZGREw8PDOnnyZMR1Fy5c0NjYmAYGBuR2u+M+KABgdRnRFiwsLOj06dPy+/3KyspSX1+fbt68qWAwuLjG6/XK5XLJ5XKppKREly5d0gcffJDQwQEAS0V9hR4Oh+X3+yVJT58+VTAYlMPhWLKmoqJCLS0tkqTe3l5lZ2crPz8/AeMCAFaypnPoTqdTbrdbvb29S+53OByanJxcvB0KhZZFX5JqampkWZYsy1Jubu5rjgwAiCTmoG/ZskVtbW06deqU5ubmlhxLS0tbtt627WX3+Xw+eTweeTwePXr06DXGBQCsJKagZ2RkqK2tTa2trWpvb192PBQKqbCwcPF2QUGBpqam4jclACCqmILe3NysYDCo8+fPRzze2dmp6upqSVJJSYmePHmicDgcvykBAFFF/SuXsrIyVVdXa3BwcPGXo2fOnNH27dslSU1NTerq6tKBAwc0Pj6uZ8+e6ejRo4mdGgCwTNSg37lzJ+I58r87fvx4XAYCALwerhQFAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwBEEHAEMQdAAwRNSgNzc3a3p6WkNDQxGPl5eXa3Z2Vn6/X36/Xw0NDXEfEgAQXUa0BZcvX9bFixfV0tKy4pqenh4dPHgwroMBANYm6iv0np4ePX78eCNmAQCsQ1zOoZeWlioQCKirq0vFxcUrrqupqZFlWbIsS7m5ufF4agDAf0Q95RJNf3+/nE6n5ufn5fV61dHRoaKioohrfT6ffD6fJMmyrPU+NQDgFet+hT43N6f5+XlJUnd3tzIzM5WTk7PuwQAAa7PuoOfl5S1+7vF4lJ6erpmZmfU+LABgjaKecrl69ar27Nmj3NxcTU5O6uzZs8rMzJQkNTU1qbKyUseOHdPCwoKeP3+uqqqqhA8NAFguatAPHz686vHGxkY1NjbGbSAAwOvhSlEAMARBBwBDEHQAMARBBwBDEHQAMMS6rxQ11b+Hfol4/+l/libl8RM9z2oSPetG/NuSuX+RpNo8MAOv0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEAQdAAxB0AHAEFGD3tzcrOnpaQ0NDa245sKFCxobG9PAwIDcbndcBwQAxCZq0C9fvqz9+/eveNzr9crlcsnlcqm2tlaXLl2K64AAgNhEDXpPT48eP3684vGKigq1tLRIknp7e5Wdna38/Pz4TQgAiMm6z6E7HA5NTk4u3g6FQnI4HBHX1tTUyLIsWZal3Nzc9T41AOAV6w56Wlrasvts24641ufzyePxyOPx6NGjR+t9agDAK9Yd9FAopMLCwsXbBQUFmpqaWu/DAgDWaN1B7+zsVHV1tSSppKRET548UTgcXvdgAIC1yYi24OrVq9qzZ49yc3M1OTmps2fPKjMzU5LU1NSkrq4uHThwQOPj43r27JmOHj2a8KEBAMtFDfrhw4ejPsjx48fjMgwA4PVxpSgAGIKgA4AhCDoAGCLqOfT/r/71X5H/Tv5f9vUlt9PTDm7EOHjF//ztexBJKn5fls4d28/XSjb63xfLnm+0eO5BvH6mYt2nRH3/eIUOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIYg6ABgCIIOAIaIKej79u3T6OioxsbGVF9fv+x4eXm5Zmdn5ff75ff71dDQEPdBAQCry4i2ID09XY2Njdq7d69CoZAsy1JnZ6eCweCSdT09PTp48GDCBgUArC7qK/Tdu3drfHxcExMTevHiha5du6aKioqNmA0AsAZRg+5wODQ5Obl4OxQKyeFwLFtXWlqqQCCgrq4uFRcXR3ysmpoaWZYly7KUm5u7jrEBAH8X9ZRLWlrasvts215yu7+/X06nU/Pz8/J6vero6FBRUdGyr/P5fPL5fJIky7Jed2YAQARRX6GHQiEVFhYu3i4oKNDU1NSSNXNzc5qfn5ckdXd3KzMzUzk5OXEeFQCwmqhBtyxLLpdLO3bsUGZmpqqqqtTZ2blkTV5e3uLnHo9H6enpmpmZif+0AIAVRT3l8vLlSx0/flw//fST3njjDX3zzTcaGRlRXV2dJKmpqUmVlZU6duyYFhYW9Pz5c1VVVSV8cADAUlGDLv3faZTu7u4l9zU1NS1+3tjYqMbGxvhOBgBYE64UBQBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMARBBwBDEHQAMERMQd+3b59GR0c1Njam+vr6iGsuXLigsbExDQwMyO12x3VIAEB0UYOenp6uxsZGeb1eFRcX65NPPtF77723ZI3X65XL5ZLL5VJtba0uXbqUsIEBAJFFDfru3bs1Pj6uiYkJvXjxQteuXVNFRcWSNRUVFWppaZEk9fb2Kjs7W/n5+YmZGAAQUZoke7UFH3/8sfbv36+amhpJ0qeffqqSkhKdOHFicc3169d17tw53blzR5J069Yt1dfXq6+vb8lj1dTUqLa2VpL0zjvv6O7du6sOl5ubq0ePHq35H7WRmDF+NsOczBgfzPj6nE6ntm3bFvFYRrQvTktLW3afbdtrXiNJPp9PPp8v2lMusixLHo8n5vXJwIzxsxnmZMb4YMbEiHrKJRQKqbCwcPF2QUGBpqam1rwGAJBYUYNuWZZcLpd27NihzMxMVVVVqbOzc8mazs5OVVdXS5JKSkr05MkThcPhxEwMAIjoDUn/vdoC27Y1Njam1tZWnThxQt9++62+++471dXV6f3331dfX5/GxsZUWlqqr776Sl6vV7W1tXrw4EFcBuzv74/L4yQSM8bPZpiTGeODGeMv6i9FAQCbA1eKAoAhCDoAGCKlgl5ZWanh4WG9fPlSu3btWnHdxMSEBgcH5ff7ZVnWBk4Y+4yxvF1ComzdulU3btzQvXv3dOPGDWVnZ0dcl4x93AxvIxFtxvLycs3Ozsrv98vv96uhoWHDZ2xubtb09LSGhoZWXJPsfYw2YyrsY0FBgW7fvq2RkRENDw/r5MmTEdcley/Xwk6Vj3fffdcuKiqyf/75Z3vXrl0rrpuYmLBzcnJSdsb09HR7fHzc3rlzp52ZmWkHAgH7vffe27AZv/zyS7u+vt6WZNfX19vnzp1LiX2MZV+8Xq/d1dVlS7JLSkrsX3/9dUO/v7HMWF5ebl+/fj0pP39/fXz44Ye22+22h4aGIh5P9j7GMmMq7GN+fr7tdrttSXZWVpZ99+7dlPuZXMtHSr1CHx0d1b1795I9xqpimTGWt0tIpIqKCl25ckWSdOXKFR06dGjDnns1m+FtJJL9vYtVT0+PHj9+vOLxZO+jFH3GVBAOh+X3+yVJT58+VTAYlMPhWLImFfYyVikV9FjZtq0bN27ot99+W3xLglTicDg0OTm5eDsUCi37IUmkvLy8xesAwuHwipcJb/Q+xrIvyd67WJ+/tLRUgUBAXV1dKi4u3rD5YpXsfYxVKu2j0+mU2+1Wb2/vkvs3y15KMVz6H283b96M+H+3L774YtkFSyspKyvTgwcP9NZbb+nmzZsaHR1VT09PyswY61shrMdqM8Yq0fv4d/F8G4lEieX5+/v75XQ6NT8/L6/Xq46ODhUVFW3UiDFJ9j7GIpX2ccuWLWpra9OpU6c0Nze35Nhm2Mu/bHjQ9+7du+7H+Ouipd9//13t7e3avXt3XEO03hk34q0QVptxenpa+fn5CofDys/P18OHDyOuS/Q+/t1meBuJWJ7/1f/gu7u79fXXXysnJ0czMzMbNmc0yd7HWKTKPmZkZKitrU2tra1qb29fdnwz7OVfNt0pl3/84x/Kyspa/Pyjjz7S8PBwkqdaKpa3S0ikzs5OHTlyRJJ05MgRff/998vWJGMfN8PbSMQyY15e3uLnHo9H6enpKRVzKfn7GItU2cfm5mYFg0GdP38+4vHNsJevSvpvZv/6OHTokD05OWn/8ccfdjgctn/88Udbkv3222/bP/zwgy3J3rlzpx0IBOxAIGAPDw/bZ86cSbkZ9Z/fjN+9e9ceHx/f8BnffPNN+9atW/a9e/fsW7du2Vu3bk2ZfYy0L3V1dXZdXd3imosXL9rj4+P24ODgqn/tlKwZP//8c3t4eNgOBAL2L7/8YpeWlm74jFevXrWnpqbsP//8056cnLQ/++yzlNvHaDOmwj6WlZXZtm3bAwMDtt/vt/1+v+31elNuL2P94NJ/ADDEpjvlAgCIjKADgCEIOgAYgqADgCEIOgAYgqADgCEIOgAY4n8BpjA8s1FKTjcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu1 = -1\n",
    "mu2 = 1\n",
    "sigma = 1\n",
    "N = 10 # number samples\n",
    "s1 = np.random.normal(mu1, sigma, N)\n",
    "s2 = np.random.normal(mu2, sigma, N)\n",
    "plt.hist(s1, 30, density=True)\n",
    "plt.hist(s2, 30, density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Research design:\n",
    "\n",
    "* how to determine sample size? (aka power calculations)\n",
    "    \n",
    "    - you essentially simulate your statistical analysis\n",
    "    - you need to make meaningful assumptions, e.g. group difference, variability, power\n",
    "    - you want to know how many samples you need so that you can reliably identify the hypothesized effect\n",
    "\n",
    "* many tools available, e.g. [G*Power](http://www.gpower.hhu.de/), [WebPower online](https://webpower.psychstat.org/wiki/), [powerandsamplesize.com](powerandsamplesize.com), ...\n",
    "* for Python: [StatsModels package](http://www.statsmodels.org/stable/index.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Research design - let's compute sample size\n",
    "\n",
    "This is for a repeated measures t-test..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUddr/8feUzEx6h/SEktAERIooIAIisLiigiviigWBVdFl19+Kq+vio89aEEV9RMVIWVxcQMUFZBUQAQstUoNJgEDakF5IbzNzfn+EDAkEkkCSk0zu13WdK1POmblnyPXh5Hu+RQMoCCGE6PC0ahcghBCiZUigCyGEg5BAF0IIByGBLoQQDkICXQghHIQEuhBCOAgJdCGuwl133UVqairFxcVcf/31apcjhJ0im2xtvSUlJSllZWVKcXGxkpmZqaxYsUJxdXVVva6mbomJicqdd955yeOhoaFKcXGxfVMURSkpKbHfHzly5FW/Z0ZGhjJixAjVP7ts7XpTvQDZOuGWlJSkjBs3TgGUoKAgJTY2VnnttddUqUWn0zX7mOrqaqVHjx6N7qcoSpP2a8omgS5bY5s0uQjVpaen880333DdddcBEBgYyMaNG8nLy+PUqVM89thjABiNRsrKyvD19QXghRdeoLq6Gnd3dwBeeeUVlixZAoDBYODNN98kJSWFzMxMPvzwQ0wmEwCjR48mLS2NZ599loyMDFauXHlJTRqNhhdeeIHk5GSysrL45z//iYeHBwaDgeLiYvR6PUePHiUxMbHZn9dkMrFkyRJSU1PJyMjgvffew2AwANC1a1e++eYbCgoKyM3NZceOHQCsX7+eLl26sG3bNoqLi3nqqaea/b6ic1D9fxXZOt9W9ww9JCREOX78uPLyyy8rgLJr1y5l6dKlitFoVAYOHKhkZ2crY8eOVQBl9+7dyj333KMAytatW5XExERl4sSJ9ufuuusuBVCWLFmibNy4UfH29lbc3NyUTZs2Ka+++qoCKKNHj1aqq6uV119/XTEYDIrJZLqkvkceeUQ5deqU0q1bN8XV1VX58ssvldWrV9ufb+qZd0P7ffjhh8rnn3+ueHp6Kh4eHsq3336r/P3vf1cA5e2331beeecdRafTKU5OTsqoUaPsx8kZumxN2FQvQLZOuCUlJSnFxcVKQUGBkpycrCxdulQxmUxKSEiIYrFYFDc3N/u+r776qrJy5UoFUF5++WXl3XffVXQ6nZKRkaE8/fTTymuvvaYYjUalrKxM8fX1VQClpKRE6d69u/01hg8frpw5c0aBmkCvrKxUjEbjZev77rvvlMcff9x+PyoqSqmqqrI3z1xtoOt0OqWyslIJCgqyP3brrbcqcXFxCqC88cYbyvr165Vu3bpd8loS6LI1tukRQiV33XWXvUmhVlBQEPn5+ZSUlNgfS0lJYciQIQDs3r2bt99+mxtuuIHY2Fi2b9/O8uXLGT58OImJieTl5eHv74+rqysHDx60v4ZGo0Gn09nv5+TkUFlZednagoKCSElJqVeDk5MTXbt2JT09/ao/c1BQEAaDgV9//bVebRaLBYB//OMfvPzyy+zcuZPq6mo++OADezOSEI2RQBftSnp6Oj4+Pri5udlDPSwsjLNnzwKwZ88eevXqxd13383u3buJj48nLCyMyZMns3v3bgByc3MpKyujX79+lw1fRVEarSM8PNx+PywsjOrqarKysq7p82VkZFBdXU2PHj3Iz8+/5PmioiLmz5/P/PnzGTBgADt37mT//v3s2bOn0ZqFkIuiol0xm83s2bOH1157DaPRSP/+/Zk1axZr1qwBoLy8nIMHD/Lkk0/aA3zPnj3MnTvXfl9RFKKjo1myZAn+/v5AzZnx7bff3uQ6/v3vf/OnP/2JiIgIXF1defXVV1m3bh1Wq/WaPp/FYmHFihW8++679ou7ISEh3HbbbQD89re/JSIiAoDCwkKsVqv9PbOysujevfs1vb9wfKq3+8jW+ba6F0Uv3oKDg5XNmzcreXl5SmJiojJ37tx6z7/66qtKWVmZYjAYFEB58sknFUVRlC5dutj3MRqNyj/+8Q/l9OnTSmFhoRIXF6c89dRTCtS0oaelpV2xPo1Go7z44otKamqqkp2drXz66aeKl5eX/flruShqMpmUN954Q0lKSlIKCwuV48ePK3/4wx8UQHn22WeV5ORkpbi4WElJSVGeffZZ+3HTpk1TUlNTlYKCAuXJJ59U/d9Qtva3ac7fEEII0cFJk4sQQjgICXQhhHAQEuhCCOEgJNCFEMJBqNYPPTs7u97ADSGEEI0LDw+nS5cuDT6nWqCnpKQwdOhQtd5eCCE6pJiYmMs+J00uQgjhICTQhRDCQUigCyGEg5BAF0IIByGBLoQQDqLRQF++fDlZWVnExsZedp93332XU6dOcfToUQYNGtSiBQohhGiaRgN91apVTJw48bLPT5o0icjISCIjI5kzZw4ffvhhixYohBCiaRrth/7jjz/Wm+j/YlOmTGH16tUA7N+/Hy8vLwICAsjMzGy5KoUQ16R2xSatVmv/2dim0Wjq/ay9Xffxi283daut6XKP1d6+3GMXP97Q/Yb2v9xzTbnflNe8nIv3Kygo4KeffmrSsc1xzQOLgoODSUtLs983m80EBwc3GOizZ89mzpw5APj5+V3rWwvR4Tk5Odk3g8GAXq+336+9rdfrG9x0Oh06na7e7ctton3Ralvn8uU1B3pD/0Ndbqms6OhooqOjgSuPdhKiI9HpdDg7O2Mymeyb0Wi8ZDMYDPafBoMBJyenJp/hXSubzYbVaq3380qboiiX3FYUpd7tix+rfZ/ax2q32sfr/mzoubqPXXz/Ss9d7n5Dt5vzXEP3m/Ial1N3v+rq6iYd01zXHOhms5nQ0FD7/ZCQkGtaRFeI9sJoNOLq6oqrqysuLi72n87Ozvafzs7OODk5XfV7VFdXU1VVRXV1db3NYrFgsVjst+s+ZrVaL7lts9ns9+tutQEuOodrDvRNmzYxb9481q5dy4033khhYaG0n4sOwWg04uHhgYeHB+7u7ri7u+Pm5oa7uzuurq5NDmqLxUJFRYV9q6ysrHe7qqqKyspK++3arbq6WhZ+Fi2q0UD/7LPPuPXWW/Hz8yMtLY2FCxfaf9GXLVvGf//7X37zm9+QmJhIWVkZjzzySKsXLURTaTQaPD098fb2xsvLCy8vLzw8PPD09MRkMl3x2KqqKkpKSigtLaWsrIzS0lLKy8spKyujrKyM8vJyysvLW+3PZyGaq9FAnzFjRqMvMm/evBYpRohrYTKZ8PPzw9fXFx8fH3x9ffH09LzsRcGqqiqKioooLi62/ywuLqakpISSkhIJatHhqDZ9rhDXwsnJiS5duuDv70+XLl3w8/PDzc2twX2Lioo4d+4cBQUFnDt3jsLCQgoLCykvL2/jqoVoXRLookNwcXEhMDCQgIAAAgIC8PHxuaSHSFVVFXl5efYtPz+fgoICLBaLSlUL0bYk0EW7ZDAYCA4OJiQkhMDAQLy8vOo9b7VaycvLIzs7m+zsbHJycigsLFSpWiHaBwl00W74+voSHh5OaGgo/v7+9QZfVFVVkZmZSUZGBllZWeTk5Eh3PCEuIoEuVKPRaAgKCqJbt26Eh4fj6upqf85qtZKeno7ZbObs2bPk5uZKFz8hGiGBLtqURqMhMDCQnj17EhERUa/rYGlpKSkpKaSmppKeni5t30I0kwS6aBPe3t706tWLHj161DsTLygoICkpieTkZHJzc1WsUIiW56l1ppvBjx4GP7rZN19OV+XwVPq6Fn8/CXTRapycnOjZsye9evWiS5cu9seLioo4deoUZ86coaCgQMUKhbh2OrSEG3zoaehCT4M/PYz+NT8N/vjqG+5K66o1tkotEuiixXl7e9O3b1+ioqLso4orKytJTEzk1KlTZGdnq1yhEM3nhI7uRj+iDF2JMnYlytiFKENXuhv9MGgajtIyWxWnq3JIqsolqSqXM3V+tgYJdNFiQkNDGTBgAMHBwfbH0tPTSUhIICkpSXqliA5Bg4YwJ2/6mgLpbQywb90MfjhpGh51bK4uILEym8SqnJqtMpszVblkWNq2K60EurgmGo2GyMhIBgwYgI+PD1Azg+DJkyeJi4uTJhXRrrlqjfQxBtDPFEQ/YyB9jIH0MQU02CRiU2wkVeVyojKLU5VZnKzM5lRVNqcqsylTqlSo/lIS6OKqaLVaoqKiGDRoEO7u7gCUlJRw/Phx4uPjZR4U0e746FzobwqmvymY60zBXGcMopvBF63m0sUmMqoLSajMJL4yk4TKjPMhnk250r5/ryXQRbNoNBqioqIYPHiwfe6UgoICjhw5wunTp+2LFQihJi+tMwOdQxhoqtkGOIcQ4uR9yX5VioUTFRnEVWTwa0U6cZUZxFdmkG8tU6HqayeBLpqsW7duDB061D4MPz8/n0OHDpGUlCSDfoRqjBo915mCuMEUxvXOoQxyDqWb4dIlLkttlfxakc6xirMcr0jneMVZTlZmU43jXNuRQBeN8vf356abbiIgIACAwsJCfvnlF06fPq1yZaIzCtF7MdglnCHO4Qx2DqOfKeiSXibltmqOV5zlaIWZo+VmjlaYOV2Vgw3HPvGQQBeX5ezszLBhw+jVqxcAZWVlHDp0iISEBGlaEW1Ch5brTEEMdQ5nqEsEQ50jCHTyrLePTbERX5HBoYo0DpencqQ8jROVWVjofL+jEujiEhqNhj59+jBs2DAMBgNWq5XY2FgOHz4sFztFqzJq9NzgHMZwl27c6NyNIS7hl/Q4KbCWcbA8hYNlKRwsT+VwRRoltkqVKm5fJNBFPb6+vowaNco+sjMlJYW9e/dSVFSkcmXCEZk0em5wDmeESw9udu3OIFMYRm39WDpTlUtMWTIHypOJKUsmsSoHxcGbTq6WBLoAarohDho0iEGDBqHVaikpKWHPnj0kJyerXZpwIHq0XO8cyijXnoxw6ckQ5/B6AW5TbPxakc6+siT2lZ3hQHky2ZZiFSvuWCTQBT4+PvaFwBVF4fjx48TExEjzimgRUYYu3OIaySjXSG526YGb7kITik2xEVtxlj2lp9lbdob9ZUmcs8nSgFdLAr2TGzBgAEOHDkWn01FUVMSuXbvIzMxUuyzRgXlpnbnFNYrRbpGMdo0i2Kn+alOnKrP4sTSRn8tOs6f0tAR4C5JA76ScnZ0ZM2YMISEhAPz666/s379f5iAXzaZBwwBTMGPdejPGNYobnMPQ1Rl9mWMp5ofSU+wuPcVPpYltPr9JZyKB3gkFBwczduxYnJ2dKS8vZ9euXaSlpaldluhAXLVGbnWNZJxbH8a59aaL3t3+XJViYW/pGXaVnmRXyQniKjPlImYbkUDvZAYNGsSQIUPQaDSYzWZ27txJebn8ySsaF6T35Hb3vkxw68tNLj3qXcw8W13AjpIT7ChJ4Oey05RKN0JVSKB3EgaDgTFjxhAeHo6iKPzyyy8cOnRI7bJEO9fbGMAk9+uY6NaXAc4h9sdtio0DZcl8VxLP9pJ4Eirlukt7IIHeCXh6ejJhwgS8vLyoqKhg586d0sQiLmuQKZTJHv35jft19eZEKbVVsqvkJFtL4vi+JIE8a6mKVYqGSKA7uMDAQMaPH4/JZCIvL49t27ZRXCz9esUFGjQMcQ7jDvcBTPa4juA6sxLmWkrYWhLHN8XH+ak0kUpFLpq3ZxLoDqxXr16MGjUKrVZLcnIy33//vfRiEXY3mMKY4jGQOzz6E1Sna+HZ6nN8U3ycLUWxHChPdvgJrRyJBLqDGjx4MIMHDwbg6NGjHDhwQKa4FfQxBnC3x/VM8bieMIOP/XFzdQGbi47xdVEshyvSpFdKByWB7mA0Gg0jR46kT58+2Gw2fvrpJxISEtQuS6goRO/F3Z6DuMdjEL1NAfbHM6oL2Vx0jI1FRzlUkapihaKlSKA7EJ1Ox9ixY+nWrRsWi4UdO3aQkpKidllCBW5aI3e49+dez8Hc7NrD/ni+pZTNxTUhvq8sSc7EHYwEuoPQ6XRMmDCBkJAQKisr+fbbb8nKylK7LNGGtGgY6dqT+zyHMMn9Opy1TkDNYg9bi39lQ9FhdpWcdKgVekR9EugOQK/XM3HiRIKCgigrK2PLli0UFBSoXZZoI+FOPkz3Gsq9noPrzZvyc+lpPi88yJbiWJkvvJNoUqBPmDCBd999F51OxyeffMIbb7xR73kPDw/+9a9/ERYWhl6vZ/Hixaxatao16hUXcXJyYtKkSQQEBFBaWsrXX39NYaHMleHoTBo9k937c7/XMEbUaVJJrspjfeEvfF54CHO1/Kfe2TQa6FqtlqVLlzJ+/HjMZjMxMTFs2rSJ+Ph4+z5PPvkkcXFx3Hnnnfj5+XHixAnWrFkj06+2Mr1ebw/z4uJivv76a+lj7uCiDF35vfcwpnkOxlvnAkC5rYrNRcdYW/iLtIt3co0G+rBhw0hMTCQpKQmAtWvXMmXKlHqBrigK7u41k/O4ubmRn58v/Z1bmU6nY+LEifYw37x5MyUlJWqXJVqBQaPjDvcBzPQezo0u3eyPHylP47NzB/hP0VGKbRUqVijai0YDPTg4uN4wcbPZzI033lhvn/fff59NmzaRnp6Ou7s79913X4N9nmfPns2cOXMA8PPzu+R50TRarZbx48cTFBREaWkpW7ZskTB3QCFO3jzkNZz7vYbiq3cDoNhawZdFh1hTcIDjlekqVyjam0YDXaPRXPLYxWE9YcIEjhw5wtixY+nRowfbt29n4MCBl/z5Hx0dTXR0NAAxMTHXUnenpdFoGDduHGFhYZSXl7NlyxZZ79OBaNAw2jWSh71vZrxbb7Tn5xWPrTjL6oJ9bCg8TJlSpXKVor1qNNDNZjOhoaH2+yEhIaSn1z8zeOSRR3j99dcBOH36NElJSfTu3VtCuxWMGDGCbt26UVlZyZYtWzh37pzaJYkW4KIxcK/XYGZ530yksSsAlTYLm4sOszJ/rwz8EU3SaKDHxMQQGRlJREQEZ8+eZfr06cyYMaPePqmpqYwbN46ffvqJLl260KtXL86cOdNqRXdWgwYNom/fvlgsFrZu3Up+fr7aJYlrFKL34lGfEczwGoanzhmA9OpzrC7Yx7/O7ZcZDUWzNBroVquVefPmsXXrVnQ6HStWrCAuLo65c+cCsGzZMl555RVWrVrFsWPH0Gg0LFiwgLy8vFYvvjPp1asXQ4cOxWaz8f3338u6nx3cIFMof/C9hd+4X4deowPgQFkS0fk/8W3xr1iwqVyh6Ig0oE4fp5iYGIYOHarGW3c4wcHBTJo0Ca1Wy48//livh5HoODRouN2tD4/7jrb3VqlWrGwuOsbH+T9ytMKscoWiI7hSdspI0XbOy8uL8ePHo9VqOXz4sIR5B2TQ6LjHYxBP+t5KT2MXAAqt5XxasI8VBXtk0WTRYiTQ2zGj0cjEiRMxGAycOXNGLjJ3MC4aAw96D2euzygCnTyBmmlqP877kc8KY2TdTdHiJNDbqdq+5h4eHuTk5LBr1y61SxJN5KV15lGfETzmM9I+mjOuIoMP8naxseiotI+LViOB3k7ddNNN9oFDW7dulZG3HYCvzpU/+NzCw94346YzAnCgLJn38r5nR4nMSS9anwR6OxQZGUm/fv2wWq1s27aNsrIytUsSV+Cvc+NJ31uZ6T0cZ60BgF0lJ3g393v2lSepXJ3oTCTQ2xlfX19GjRoFwE8//UROTo7KFYnLuRDkN9nnHt9aHMc7uTs4UpHWyNFCtDwJ9HbEaDRy++23o9friY+P58SJE2qXJBrgq3PlSd9bedj7ZnuQf1N8nLdzvpP5VYSqJNDbkVtvvRV3d3eys7P5+eef1S5HXMRL68wffG/hMZ+RuGpr2sj/WxTL27nf8WtlhsrVCSGB3m7079+f8PBwKioq2L59Ozab9IRoL1w0Bmb7jORx39H24fnbiuN4M2ebnJGLdkUCvR3w9/e3T0m8e/duSktl/o72wAkdD3rfyHy/cfjra+b7/6H0FG9kb5XJskS7JIGuMoPBwLhx49BqtcTGxpKSkqJ2SZ2eBg1TPAaywH8CEQZfAA6Vp/Jq9jf8XHZa5eqEuDwJdJWNGjXKPnho//79apfT6Y1w6cHfu0xmgHMIACcrs3gt+1u+LflV5cqEaJwEuop69uxJjx49qKqqYseOHdJurqIoQ1de7PobbnPrA0BGdSFv5mxjfeFBrDKyU3QQEugqcXV1ZeTIkQDs3btXVh1SiZ/Ojb/4384DXsPQabQUWyt4P28n0fk/Ua7IIueiY5FAV8mYMWMwGAwkJydLf3MVGDV6HvMZyR99x+KuM2FRrKzM38NbudtlUQnRYUmgq6B///4EBQVRVlbGDz/8oHY5nc4k934s7HIH4ecveH5XEs/LWVs4VZWtcmVCXBsJ9Dbm5eVln5z+hx9+oKKiQuWKOo/exgBe7vpbRrlGApBQkclL2ZvZXXpK5cqEaBkS6G1Io9Fwyy23oNfrOXHiBKmp0pe5LXhqnfl//uN5xPtmdBot+ZZSFuVu418F++WCp3AoEuhtqF+/fgQEBFBaWsrevXvVLsfhadBwn+dgXujyG/z0blgVG8vzf2ZxzjbO2crVLk+IFieB3kbc3d0ZNmwYUDOLYlVVlcoVObbrjEG8FnA3Q1zCAdhbeoYXsv5DfKUsri0clwR6G6ltaklMTJTRoK3IXWviWf/b7c0rmdVFvJz9NV8VHVG7NCFanQR6G+jVqxfBwcGUl5fLLIqt6E73gbwc8Fu66j2wKFY+zvuRN3O3USJrd4pOQgK9lZlMJvvEW3v27KGyUsKlpYU7+fBawN2McesFQExZMgsyN0jziuh0JNBb2Y033ojJZMJsNnP6tEzs1JL0aJnrewvP+I3HWetEgbWM/83+L/8+F4OConZ5QrQ5CfRWFBgYSK9evbBYLPz0009ql+NQrjeFsjhwKv1MQQB8UXiQl7K+llGeolOTQG8lWq3WPlfLkSNHZK6WFuKscWKB/wQe8xmJTqMlpSqPBZkbZHCQEEigt5oBAwbg7e3NuXPnOHJEeli0hJEuPVkcOJVwgy8WxcrSvF28lbNdJtES4jwJ9Fbg6urKoEGDgJo+5zIt7rVx15r4e5fJ/N675uLy8Yp0nsn4nGMVZ1WuTIj2RQK9FQwbNgwnJyfOnDlDerqsOXktxrhG8WbgNIKdvKi0WViS+x1L83ZhkSH7QlxCAr2Fde3alcjISCwWC/v27VO7nA7LXWtiYdc7eMCrZnTtofJU/pT+OSerslSuTIj2SwK9BWk0GkaMGAHA0aNHKSkpUbmijmmUS0+WBN1LsJM3FbZq3szZxrL8H2UiLSEaIYHegnr37o2fnx/FxcVyIfQqOGuceLHLZB7xuRmAw+WpzE9fz0mZp1yIJpFAbyFOTk4MGTIEgH379mG1WlWuqGMZ7BzGe0HT6W7wo0qx8HbOd7yft0vOyoVoBm1TdpowYQIJCQmcOnWKBQsWNLjP6NGjOXz4MMePH2fXrl0tWWOHcP311+Ps7ExGRgZJSUlql9NhOKFjgf8ENoY/QXeDH/EVGfwm6f94N+97CXMhmqnRM3StVsvSpUsZP348ZrOZmJgYNm3aRHx8vH0fT09PPvjgAyZOnEhaWhr+/v6tWnR74+rqSv/+/QHkQmgz9DT4837Q/Qx0DsGm2Fiat4tFOVupUuSvGyGuRqOBPmzYMBITE+1nnWvXrmXKlCn1An3GjBls2LCBtLQ0AHJyclqp3PZp2LBh9qlxO9tnv1ozvYbzUtc7cNYaSKvK5+n0dewrl79shLgWjTa5BAcH24MawGw2ExwcXG+fqKgovL292blzJ7/88gsPPvhgg681e/ZsYmJiiImJwc/P7xpLbx/8/PyIjIzEarUSExOjdjntno/OhZUhD/FG4D04aw18fu4g45LekTAXogU0eoau0WgueUxR6s9kp9frGTx4MOPGjcPZ2Zm9e/eyb98+Tp2qP79GdHQ00dHRAA4TfsOHDwfg+PHjFBcXq1xN+zbKpSfvBU0nwMmDQms5z2ZsYFPxUbXLEsJhNBroZrOZ0NBQ+/2QkJBLRj+azWZyc3MpKyujrKyMH374gYEDB14S6I4mJCSEoKAgKioqOHz4sNrltFt6tCzwn8A8vzEA7Cs7w7yzazlrOadyZUI4lkabXGJiYoiMjCQiIgInJyemT5/Opk2b6u2zceNGRo0ahU6nw9nZmRtvvLFeG7ujql0j9MiRI7JG6GWEOfmwMeIJ5vmNwaJYWZSzjakpyyTMhWgFjZ6hW61W5s2bx9atW9HpdKxYsYK4uDjmzp0LwLJly0hISODbb7/l2LFj2Gw2PvnkE3799ddWL15N3bt3x8/Pj9LSUof/rFfrt+4DWBw4FQ+dM+bqAp44+xkx5bKeqhCtRQPqLO0SExPD0KFD1Xjra6bRaLj33nvx8vLixx9/7BR/jTSHSaPnf7reyUzvmusLW4pieSbjCwpt5SpXJkTHd6XslJGiVyEqKgovLy8KCwtJSEhQu5x2pbvBj4+Df08/UxCVNgsvZW9mVcFetcsSolOQQG8mnU7H4MGDAfjll18u6fHTmd3pPpC3AqfhpjNypiqXueZ/cbxSpg8Woq1IoDdT7969cXNzIzc3VxZ9Ps+g0fFSl9/aJ9X6T+ER/pL5JSW2SpUrE6JzkUBvBp1OZ1+J6ODBgypX0z6EOHkTHfx7rncOpdJm4e9Zm1h9TqY/EEINEujN0KdPH1xcXMjJySElRXprjHXtzfvB0/HWuZBalc+cs//iaIVZ7bKE6LQk0JtIp9Nx/fXXA3J2rkXDM37j+bP/bQBsK47j6fR10otFCJVJoDdR3759cXFxITs7m9TUVLXLUY2X1pkPgmcwxq0XVsXGopyt/F/eLhR1er8KIeqQQG8CvV7PwIEDgc59dt7PGMjykJmEG3zJt5Tyh7Nr+LEsUe2yhBDnSaA3QW3beXZ2dr2ZJzuTqR6DeDNwGs5aJ46Wm3nMvBqzDN8Xol2RQG+ETqfr1GfnOrT8vctk5viOAuCzcwd4PvM/VCoWlSsTQlxMAr0RUVFRuLi4kJub2+nOzn11riwL/j0jXHtQpVh4MVO6JArRnkmgX4FGo7H3bPF/QfoAACAASURBVOls0+P2MwayKvRhQpy8ybYU85h5tUysJUQ7J4F+BT179sTd3Z2CggKSk5PVLqfN3OHen3eD7sNFa+BgeQqPmT8l01KkdllCiEZIoF+GRqOxjwo9cuRIp5izRYOGv/iP509+Nf3L1537hQWZG6S9XIgOQgL9Mrp164aXlxdFRUUkJjp+1zwXjYH3g6czyf06rIqN/8n6muiCn9QuSwjRDBLol1Hbdn706FGHPzsP0XuxKvRh+pmCKLSWM/fsv9hd6tjLBwrhiCTQGxAcHIyfnx9lZWWcPHlS7XJa1RDncFaGPISf3o3EyhweNq/idFWO2mUJIa6CBHoDavudHz9+HKvVqnI1recej0G8HXgvRq2e3SUnmXt2jczHIkQHJoF+EV9fX0JCQqiqqiIuLk7tclrFxRc/V+T/zN+zNmPFpnJlQohrIYF+kdqz84SEBKqqqlSupuWZNHreCbqPKR4DsSo2XszaxMqCPWqXJYRoARLodbi7u9O9e3dsNhuxsbFql9Pi/HRurAp9iMHO4RRZy/nD2TXsLHXsawRCdCYS6HX0798frVbLyZMnKS0tVbucFhVl6MKnoY8SZvAhrSqfB80rOVGZpXZZQogWJIF+ntFopHfv3kBNV0VHMsKlB8tDZuKpc+ZweSoPpa0ix1qidllCiBYmgX5enz590Ov1pKWlUVBQoHY5LeZez8G8FTgNJ42OLUWxPJW+lnKlWu2yhBCtQAId0Gq19OvXD4Bjx46pXE3LecZvPP/PfzwAH+Tt5n+z/ysrCwnhwCTQge7du+Pq6kp+fj5nz55Vu5xrpkfLm4FTme41FKti44WsjfyzYK/aZQkhWpkEOjBgwAAAh+jZ4qY18knwg4x2i6LMVsXcs2v4riRe7bKEEG2g0wd6YGCgfZh/R5+Eq6veg3+FPsp1piByLSX8Pm0FRyvMapclhGgjnT7Q+/fvD0BcXFyHHuYfaejCZ2GzCHHy5nRlDg+kLSelOl/tsoQQbahTB7qHhwfh4eFYLJYOPcx/qHM4/wx9BG+dCwfLU5iZtpJ8a5naZQkh2linDvR+/fqh0Wg4ffo0FRUVapdzVSa49eXD4Adw1jqxtTiOx8+ukW6JQnRSnTbQnZyc6NWrF1Azq2JH9IDXMN4IuAedRsu/CvbzXOZXMsGWEJ1Ypw30yMhIDAYDGRkZ5OXlqV1Os/3JbxzP+k8A4M2cbbyd+53KFQkh1KZtyk4TJkwgISGBU6dOsWDBgsvuN2TIECwWC1OnTm2xAlvLddddB3S8s3MtGv7RdQrP+k/AqthYkLFBwlwIATQh0LVaLUuXLmXSpEn07duX+++/nz59+jS43xtvvMHWrVtbpdCWFBISgpeXFyUlJSQnJ6tdTpMZNDo+CJ7Boz4jqLBVM+fsv1h9bp/aZQkh2olGA33YsGEkJiaSlJREdXU1a9euZcqUKZfs99RTT/Hll1+SnZ3dKoW2pNph/r/++muHWS/URWNgdcgjTPEYSJG1nBlpy/lvccf660II0boaDfTg4GDS0tLs981mM8HBwfX2CQoK4u677+ajjz664mvNnj2bmJgYYmJi8PPzu8qSr42HhwdhYWFYLBYSEhJUqaG5fHQufBE+h9FuUeRYipmasoy9ZWfULksI0c40GugajeaSxy4+q33nnXdYsGABNtuVe1hER0czdOhQhg4dSm5ubjNLbRl9+/a1d1WsrKxUpYbmCNR78lX44wxyDiOlKo87kz/geGW62mUJIdqhRnu5mM1mQkND7fdDQkJIT68fKEOGDGHt2rUA+Pn58Zvf/AaLxcLGjRtbuNxro9PpOlRXxe4GP9aFzSbEyZv4igzuT1tOlqVI7bKEEO1Uo4EeExNDZGQkERERnD17lunTpzNjxox6+3Tv3t1+e+XKlXz99dftLswBevbsidFoJCsrq913VbzOGMS/wx7DT+9GTFkyD6atpNBWrnZZQoh2rNFAt1qtzJs3j61bt6LT6VixYgVxcXHMnTsXgGXLlrV6kS2lb9++AO1+mP8w5wg+DX0ED50zu0pO8Kh5tYz+FEI0SgPqrHgQExPD0KFD2+z9/P39ufvuu6moqGDNmjXtdiKuMa5RLA+ZibPWwOaiYzx59t9U0z5rFUK0vStlZ6cZKVrbVTEhIaHdhvlk9/58EHw/Bo2ez84d4C8ZX2KTFYaEEE3UKQLdaDTSvXt3FEUhPr59LvbwO8/BvB14LzqNlmV5P/BS9tdqlySE6GA6RaD37t0bvV5PamoqxcXFapdziUe8b+bVgLsAWJSzjSUylF8IcRU6RaDXTlXQHi+GPuU7hue7TAJgYdZmPs7/UeWKhBAdlcMHenBwMB4eHhQXF9cb8doePOc/gT/6jcOm2Hg2cwNrzh1QuyQhRAfm8IFe21UxISGhXc3b8nLX3zLbZxQWxcrT6ev4quiI2iUJITo4hw50FxcXwsPDsdlsnDhxQu1yANCgYVHAPfze+0aqFAt/OLuGb4p/VbssIYQDcOhA79WrF1qtlqSkJMrK1F9jU4eWJYH3cq/XYMpt1cwy/5OdpSfVLksI4SAcNtA1Gg29e/cG2sfFUD1algbfz50eAym1VfJg2kqZMVEI0aIcNtBDQkJwd3enqKiIs2fPqlqLQaNjWfADTHS/jiJrOQ+kreCX8hRVaxJCOB6HDfTai6FqDyQyafR8EjKTcW69KbCWcX/qJxytMKtakxDCMTlkoLu6uhIaGorValX1YqizxolVoQ9zi2skeZYSfpcaTVxlhmr1CCEcm0MGelRUFFqtljNnzlBRUaFKDa5aI5+GPMJNrt3JthRzb8rHnKzKUqUWIUTn4JCBXnsxVK3mFnetiTWhjzLUJYKM6kLuTf2Y01U5qtQihOg8HC7Q1b4Y6qE18e+wx7jBOQxzdQHTUpaRUp3f5nUIITofhwv02nlb1Gg799I6sy5sNgOcQ0itymdq6jLM1QVtXocQonNyqEB3dnZWbWSot86FdWGz6W8KJqkql2kpy0i3FLZpDUKI9k1jMqAP8AGthuozLd9BwqECvfZiaHJycpuODPXVubIubDb9TEGcrsxhWuoyMmUxZyE6La2XG06BvugDfdAH+qIP9MUpwBedjzsAlQmp5C76d4u/r0MFeu3F0ISEhDZ7T1+dK5+HzaGPKZDEymympi4j29L+5lwXQrQwnRa9v9eFwK4T3lpnY4OHKNUWLJn5WDJaZ5F6hwn0wMBAPD09KSkpabNpcv10bnwRPodexgBOVmYxLWUZOdaSNnlvIUTbqG0mqQnt2rNuP/RdvNDodQ0eYysppzojD8v5rfa2NbcIWnHWV4cJ9Nqz85MnT7bJNLn+Oje+CJ9LlLErCRWZ3Jv6MbkS5kJ0WFpP1zqh7Ys+qH4zSUMsuYU1oZ2eR3XmhQC3FZe3YeUXOESgGwwGunXrBrRN75Yuene+CJtDpLEr8RUZ3Jv6MXnW0lZ/XyHENdJq0Pl54RTkW6eppOasW+tiavAQpdqCJaug5iw7PQ9L5vmfWfkoVZY2/gBX5hCB3qNHD/R6PWfPnm31NUPrhnlcRQa/kzAXot3RGPToA843jwT51YR2gA/6rt5onBqOvfrNJPlYMvOoTs/DmlvYqs0kLckhAr2tLoZ20bvzZdhcehq7SJgL0Qze3t7Mnz+fiIgINBpNy72wVoNGpwWdrv5Prfbyx9hsKFYbWK3nf9pQrFbwUsAL6NNy5V0tRVFITk7mnXfeoaCg6WNZOnyg+/r64u/vT0VFBcnJya32Pl31HnwRNpeeRn+OV6RzX+rH5FvVXzRDiI5g/vz5/PLLL7z88stYrdbmv4BOi8agR6PX1/x00oOTHs1lg1tBqbaiVFvsG9UWlGprhzjb1ul0TJ48mfnz57Nw4cImH9fhA71Xr14AJCYmXt0vShNImAtxbSIiIpoW5k56NE46NE76OpsONJcJbsV2Ibir6oS3pXWyoK1YrVa2bNnC1KlTm3Vchw50nU5Hz549gda7GHpxmP8u9WMKJMyFaBaNRnMhzDWampCuG9oG/fkugA03xyhW6/kzbEu9M2+strb7EG3MarU2u3mqQwd6REQEJpOJnJwc8vJavqO+hLkQV0/r4WLvRaJ1NaE7f0FSo2u47zYoKJa6TSR1gtvW/ptJ2oMOHei1zS2tcXZe05tFwlyIK9KAztezJrgv6gqodXO+sJvJiNZ0fvSkUie4q6z12rg7iueee45Zs2ZhtVp5+umn2bZt2yX7LFq0iN/+9rdUVVVx+vRpHnnkEQoLC5t8/NXosIHu5uZGcHAwFouFxMTEFn3tC71ZLrSZS5iLTk2vQ9/Vu/7Am/NdAbVGpwYPsZVV2kdJKpNuwJJd4BDt23369GH69On069ePoKAgvvvuO6KiorDZ6jf/bN++nb/+9a9YrVZef/11/vrXv/Lcc881+fir0WEDPSoqCo1GQ3JyMlVVVS32uv46N74Im0NPYxd+lQugopOxD3MPuhDcTkF+6Py9aroDNsB6rpjq9LyaOUrSLwx1t527MHLadut0lPJKAIJXLGiV2s8++sYVn3/ggQd4+umnMRgM7N+/nyeeeILly5czZMgQFEVhxYoVvPPOO42+z5QpU1i7di1VVVUkJyeTmJjIsGHD2LdvX739tm/fbr+9b98+pk2b1qzjr0aHDfTWaG7x17nxZfjceoOGJMyFI9K6u1wY3l5nqLvex6PB/RWbcmG0ZO38JOfDuzao27PevXtz3333MWLECCwWC0uXLuVvf/sbwcHB9O/fHwBPT08A/t//+3888MADl7zGDz/8wB//+EeCg4Prha/ZbCY4OPiK7//oo4+ybt06gKs6vqk6ZKAHBQXh7u5OcXFxi61K5Hd+bhYJc+EwNKDz8bgwUrJOgNdt366rdjbA6tqz7fTcmhDPzG+xppLGzqRbw7hx4xg8eDAxMTFAzdoJ3377Ld27d+e9995jy5Yt9nbsxYsXs3jx4su+VkM9T640f9Tzzz+PxWJhzZo1V3V8czQp0CdMmMC7776LTqfjk08+4Y036v+DzJgxgwULav6MKikp4fHHH+fYsWMtUmBDas/OT5482SKv56tz5YvwOUSdn5tFRoCKDkWnvbR9O6i2fdvQ4CF127ct6bn229acjjPMvTk0Gg3//Oc/ef755+s9/sILLzBhwgSefPJJfve73zFr1qxGz9DNZjOhoaH2x0NCQkhPT2/wfWfOnMkdd9zBuHHj7I815/jmajTQtVotS5cuZfz48ZjNZmJiYti0aVO9BZiTkpIYPXo0586dY+LEiXz88ccMHz68RQq8WEtPxOWrc+WLsLn0MgbYZ02UMBftkcbgZF8wwSnIF32QX014+19+GldrYUnNTIAXNZXUbd/uDHbs2MHGjRtZsmQJOTk5eHt74+7uTkFBARs2bOD06dOsWrUKaPwMfdOmTXz22We8/fbbBAUFERkZyYEDBy7Zb8KECSxYsIDRo0dTXl7e7OOvRqOBPmzYMBITE0lKSgJg7dq1TJkypV6g792713573759hISEtEhxDamdiMtsNlNScm2/lL46V9aHzaG3KYATlZlMS10mYS5Up3U11YR10EU9Svw8G9xfsSlYcs7VD+70PKozclHK2n/7dluIj4/nb3/7G9u2bUOr1VJdXc2f//xnvvrqK7Tnpw/461//2qTXiouLY/369cTFxWGxWHjyySftPVSio6P56KOPOHjwIO+//z5Go9F+cXTfvn08/vjjVzz+WjUa6MHBwfUWjDCbzdx4442X3X/WrFl88803DT43e/Zs5syZA4Cfn19zawWw/2dxrWfnPjoX1ofNoa8p8PziFHJmLtqWztu9pmnkojZunbtLg/srFiuWrILzZ9m59uBuj9O4tkfr169n/fr19R4bPHjwVb3Wq6++yquvvnrJ47Nnz7bfjoyMbPbx16rRQG9OA/6tt97KrFmzGDlyZIPPR0dHEx0dDWC/ONFc27dvJzAwkOzs7Ks6HmoWdK4N81PnVxqSxSlEq9Bq0HfxvjDgxj74xufCQJuL2Coqa6ZvrdOTxJKRhyW7QEZMiitqNNCb2oDfv39/PvnkEyZNmkR+fn7LVnmRjIyrXy3bW+diX9A5sTKbaakfy7Jx4to56XE633+73hl3V+/Lt28Xl9VckLyoqcRaIGvSiqvTaKDHxMQQGRlJREQEZ8+eZfr06cyYMaPePqGhoWzYsIEHH3yQU6dOtVqx18pL68y6sNn0NwWTWJkjCzqLZtO4GHEKrGnftvcqCfJF5+uJRtvwREr2ZcrsvUrOL1NWos4yZcJxNRroVquVefPmsXXrVnQ6HStWrCAuLo65c+cCsGzZMv7+97/j6+vLBx98AIDFYmHo0KGtW3kzedYJ89OVOdwrYS6uQOvldukyZUG+6DzdGtxfsdqozsi/0AWwtqkkMx+lsrqNqxedVZP6oX/zzTeXXOhctmyZ/fbs2bPrXQxobzy0JtaFzWaAcwhnqnKZlrqMTEuR2mUJtTW0vuT521rny7RvV1ZdaN+uc7ZtyS5w6KlcRcfQIUeKNoeH1sTasMcY6BxCUlUu01IkzDuduu3bAXWGu19hfUlrcdklQ9wtGXlY84tArkuKdsqhA91da+LfYY8xyDmMlKo8pqUsI8NSqHZZopVc2r5dMwhH5+d1+fbtvKKG27eLZdoHcXlNmf524cKFzJ49m5ycHKBmCoDalg6ZPreZasJ8FjecD/OpKctIlzB3CFpv9/O9SHya175tH3BTp327ouVm6hSdQ3Omv12yZAlvvfXWVR/fXA4Z6G5aI5+FzmKwczipVflMS1nGWcs5tcsSzaHVoPf3unRGwCa2b9vPuDPysGRJ+3Z7ktFnUau8bmD8s1d8vq2nz22t46/E4QK9NsyHuISTVpXP1NRlmCXM2y2Nwalm/u1An+b135b2bdEMak2fO2/ePGbOnMkvv/zCM888w7lz52T63KZy1RpZEzqLoS4RmKsLuCd1GebqArXLEoDW3bl+F8BG5ieBBvpvZ9S2b0v/7Y6ssTPp1qDG9Lkffvghr7zyCoqi8Morr/DWW28xa9Ys9afP7QhctUY+C32UYbVhnvKRhHlbq7u+5EXTuDZpfpI6oyUtmfkoVdJ/W7QMNabPrTs9SXR0NF9//TWg8vS5HcGFMO/G2eoCpqYsI03CvPXodegDfOqcaZ+/OBngg8ZwmfUlyysvOts+39adc07at0WrU2P63ICAADIzMwG4++67OX78eLOOvxodPtBrmllqw/wcU1OWkVrdunPJdBYaVxNOARfOsmsH3VypG6D1XDGWjPwLXQAzO+f826J9UWP63EWLFnH99dejKArJycn20fWtOX2uBpUuI8XExFzz9AAuGgOfhc3iRnuYf0SKhHnz1C5TVtu+XWcAjs7TtcFDFKsNa865+osmnD/r7gjrS4q2t3r1ambOnKl2GR1OQ9/blbKzw56hu2gMrAl7VMK8qeo2kwT4XLgoGeCD1niZZpKKqpqV3C+6KGnJPtdi60sKIVpOhwz02jAf7tKd9OpzTEtZJmF+ntbNuaZNO6BujxKfKzeT1F2mrHZx4AyZxlWIjqbDBfrFYT41ZRnJ1Xlql9W2NBp0fp7oA2vPuJvQm8Rqw5JZUG+UZPX5Nm5ZpkwIx9DhAv3jkN93mjDXGJ3qNY841d6+wqRS9Va7ybwwalJGSwrh+DpcoL+fu5NwJx8eTFvpMGFeMzeJT014B9TpCujjcdljrAXFddq1L7RzS28SITqvDhfo+8qTuPXM21jpYGebTnr0Xb0vujBZE+KXW1tSqbbUDLrJzK/XVCKTSgkhGtLhAh1o12Gu9XC1h7VTwIVBN1daouzC3CT558M7F0tmPtbcQlkUWIh2qCnT3w4cOJCPPvoIk8mExWLhiSeesE89INPntid6Xc1K7vbg9rnQBdDF1OAhitV2vk07336WXTvwxlZa0cYfQAhxtZo6/e2iRYv4n//5H7799lsmTZrEokWLGDNmjEyfqxath0uddu3zwR3gi87fE8350WUXs5VWXLggWfeipAxxFwKAOXPmtMrrfvzxx1d8vq2nz1UUBQ+Pmutgnp6e9vlaZPrc1lR7tn2+PdupToBf9mzbZjvftn1xcOfLSjdCtENqTJ87f/58tm7dyuLFi9Fqtdx8880AMn1uS9B6nm/brntRMsAHnd8VzrbLKi6MlLQHd37NgsAyUlKIq9LYmXRrUGP63Mcff5w//elPbNiwgXvvvZfly5czfvx4mT63yZz06AO8ay5GBvhc2Lr6oHW5TE8Smw1L1vmLkXXPtjPzsRXJ2bYQjkCN6XMfeugh/vjHPwLw+eef88knnwAyfe4ldN7u9iYSext3gA86H4/L9iSxlZRfCOzMC93/ZF4SIRyfGtPnpqenM3r0aHbv3s3YsWM5depUs46/Gh0u0L1n34HLTf0afE6xWKnOOnchrGvPtrPyZZUbIToxNabPnT17Nu+++y56vZ6Kigr7xWCZPrcO97tG4nrr9fVDu3aTniRCtEsyfe7Vcfjpc4s3/kzxf35SuwwhhGh3Gu7e0Z610NVgIYRwNB0v0IUQHY6iKOh0OrXL6FB0Ol2zuzNKoAshWl1ycjKTJ0+WUG8inU7H5MmTSU5ObtZxHa4NXQjR8bzzzjvMnz+fqVOnNjiwRtRXu7B0U6YiqEsCXQjR6goKCli4cKHaZTg8aXIRQggHIYEuhBAOQgJdCCEchGojRbOzs0lJSVHjrdsFPz8/cnNz1S5DVfIdyHfQ2T8/NP87CA8Pp0uXLpd9XpGt7beYmBjVa1B7k+9AvoPO/vlb+juQJhchhHAQEuhCCOEgdMBLahfRWR06dEjtElQn34F8B53980PLfQeqXRQVQgjRsqTJRQghHIQEuhBCOAgJ9FYWEhLC999/T1xcHMePH+fpp58GwNvbm23btnHy5Em2bduGl5eXypW2Pq1Wy6FDh9i8eTPQ+b4DT09PPv/8c+Lj44mLi2P48OGd7juYP38+x48fJzY2ls8++wyj0ejw38Hy5cvJysoiNjbW/tiVPvNzzz3HqVOnSEhI4Pbbb2/We0mgtzKLxcIzzzxD3759GT58OE8++SR9+vThueeeY8eOHURFRbFjxw6ee+45tUttdX/84x+Jj4+33+9s38G7777Lt99+S58+fRg4cCDx8fGd6jsICgri6aefZsiQIfTv3x+dTsf06dMd/jtYtWoVEydOrPfY5T5znz59mD59Ov369WPixIl88MEH9jVPm0r1jvWdafvPf/6j3HbbbUpCQoISEBCgAEpAQICSkJCgem2tuQUHByvfffedMmbMGGXz5s0K0Km+A3d3d+XMmTOXPN6ZvoOgoCAlNTVV8fb2VnQ6nbJ582Zl/PjxneI7CA8PV2JjYxv9d3/uueeU5557zr7ft99+qwwfPrzJ7yNn6G0oPDycQYMGsX//frp27UpmZiYAmZmZVxzK6wjeeecdnn322Xqrm3em76B79+7k5OSwcuVKDh06RHR0NC4uLp3qO0hPT2fx4sWkpqaSkZFBYWEh27dv71TfQa3Lfebg4GDS0tLs+5nNZoKDg5v8uhLobcTV1ZUvv/yS+fPnU1xcrHY5bWry5MlkZ2d36v7Ger2eG264gQ8//JAbbriB0tJSh2taaIyXlxdTpkyhW7duBAUF4erqygMPPKB2We1KQ4t/NGcZOgn0NqDX6/nyyy9Zs2YNX331FQBZWVkEBAQAEBAQQHZ2tpoltqoRI0Zw5513kpSUxNq1axk7diyffvppp/oOzGYzZrOZAwcOAPDFF19www03dKrv4LbbbiMpKYnc3FwsFgsbNmzg5ptv7lTfQa3LfWaz2UxoaKh9v5CQENLT05v8uhLobWD58uXEx8ezZMkS+2ObNm3ioYceAuChhx5i48aNapXX6p5//nlCQ0Pp1q0b06dP5/vvv+fBBx/sVN9BVlYWaWlpREVFATBu3Dji4uI61XeQmprK8OHDcXZ2Bmq+g/j4+E71HdS63GfetGkT06dPx2AwEBERQWRkpP0koKlUv2DgyNuIESMURVGUo0ePKocPH1YOHz6sTJo0SfHx8VG+++475eTJk8p3332neHt7q15rW2yjR4+2XxTtbN/BwIEDlZiYGOXo0aPKV199pXh5eXW67+Cll15S4uPjldjYWGX16tWKwWBw+O/gs88+U9LT05WqqiolLS1NefTRR6/4mZ9//nklMTFRSUhIUCZOnNis95Kh/0II4SCkyUUIIRyEBLoQQjgICXQhhHAQEuhCCOEgJNCFEMJBSKCLJlEUhcWLF9vvP/PMMyxcuLBFXnvlypVMnTq1RV7rSqZNm0ZcXBzff//9Jc/17duXHTt2cOLECU6ePMnf/vY3+3MLFy7kmWeeafX6GvPQQw8RGBhovx8dHU2fPn1UrEi0NxLookkqKiq455578PX1VbuUepozE92sWbN44oknGDt2bL3HTSYTmzZt4vXXX6dXr14MHDiQm2++mSeeeKKly7Vr7gx6AA8//DBBQUH2+7Nnz643e6UQEuiiSSwWCx9//DF/+tOfLnnu4jPs2rlqRo8eza5du1i3bh0nTpzgtddeY8aMGezfv59jx47RvXt3+zG33XYbP/zwAydOnGDy5MlATegtWrSIAwcOcPToUebMmWN/3e+//541a9bUm2O61vTp0zl27BixsbG8/vrrALz44ouMHDmSjz76iEWLFtXbf8aMGfz8889s374dgPLycubNm1dvrpWBAweyY8cOTp48yWOPPQbUDNnevXs3hw8fJjY2lpEjRwIwfvx49uzZw8GDB1m/fj2urq4AJCUl8eKLL/Ljjz/y7LPPsn//fvvrh4eHc/ToUXutBw4cIDY2lmXLlgEwdepUhgwZwpo1azh8+DAmk4mdO3cyePDgy37m2n+L//3f/+XIkSPs3bvXPgnUtGnTiI2N5ciRI+zevfvSf3DRYak+5vP3oAAABRlJREFUkkq29r8VFxcr7u7uSlJSkuLh4aE888wzysKFCxVAWblypTJ16tR6+0LNqNCCggIlICBAMRgMitlsVl566SUFUJ5++mllyZIl9uO/+eYbRaPRKD179lTS0tIUo9GozJ49W3nhhRcUQDEYDEpMTIwSERGhjB49WikpKVEiIiIuqTMwMFBJSUlR/Pz8FJ1Op+zYsUOZMmWKAig7d+5UBg8efMkxb731lvL0009f8nh+fr7i7u6uLFy4UDly5IhiMpkUX19fJTU1VQkMDFT+/Oc/K88//7wCKFqtVnFzc1N8fX2V3bt3Ky4uLgqgPPvss8qLL76oAEpSUpLyl7/8xf76hw8fVrp162bfr/az1h01uHr1auWOO+5osP7a+1f6zIqi2I9/44037O9x7NgxJSgoSAEUT09P1X+/ZGuZTc7QRZMVFxezevVq+6pLTRETE0NmZiZVVVWcPn2abdu2ARAbG0tERIR9v/Xr16MoComJiZw5c4bevXtz++23M3PmTA4fPsz+/fvx9fUlMjISgAMHDpCcnHzJ+w0dOpRdu3aRm5uL1WplzZo13HLLLVesUaPRXHZGu9rHN27cSEVFBXl5eezcuZNhw4YRExPDI488wsKFC+nfvz8lJSUMHz6cvn378vPPP3P48GEeeughwsPD7a+3bt26ep/5d7/7HQD33Xef/bkxY8awb98+jh07xtixY+nXr98V67/SZ66srOTrr78G4ODBg/bv/Oeff2bVqlU89thj6HS6K76+6Dj0ahcgOpZ33nmHQ4cOsXLlSvtjFoulXpuwwWCw366srLTfttls9vs2mw29/sKv38WBqigKGo2Gp556yv6fQK3Ro0dTWlraYH0NTT/amF9//fWS0O/WrRslJSWUlJRctr4ff/yRW265hcmTJ/Ppp5/y5ptvUlBQwPbt25kxY0aD71W37nXr1vH555+zYcMG+39mRqORDz74gCFDhmA2m1m4cCEmk+mK9V/pM1dXV9tvW61W+3f++OOPM2zYMCZPnsyRI0e4/vrryc/Pv+L7iPZPztBFsxQUFLB+/XpmzZplfyw5OdneljtlypR6gd5U9957LxqNhu7du9O9e3dOnDjB1q1befzxx+0hFBkZiYuLyxVfZ//+/YwePRpfX1+0Wi33339/o23Ea9asYeTIkYwbNw6ouUj63nvv1WtrnzJlCkajER8fH2699VZiYmIICwsjOzubTz75hOXLl3PDDTewb98+RowYQY8ePQBwdna2/1VxsTNnzmC1WnnxxRftZ+e14Z2bm4urqyvTpk2z719cXIy7u3uLfObu3btz4MABFi5cSG5ubr0pW0XHJWfootneeust5s2bZ78fHR3Nxo0b2b9/Pzt27LCf1TbHiRMn2L17N127duUPf/gDlZWVfPLJJ0RERHDo0CE0Gg05OTncddddV3ydzMxM/vrXv7Jz5040Gg3//e9/2bRp0xWPqaioYMqUKfzf//0fS5cuRafT8emnn/L+++/b9zlw4ABbtmwhLCyMV155hYyMDGbOnMlf/vIXqqurKSkpYebMmeTm5vLwww/z73//G6PRCMDf/vY3Tp061eB7r1u3jsWLF9ubQgoLC4mOjiY2Npbk5GRiYmLs+65atYqPPvqI8vJybrrppmv6zG+++SaRkZFoNBp27NhhvyD7/9u5YxoAAACEYf5do4OlFbGDA755WwSIMLkARAg6QISgA0QIOkCEoANECDpAhKADRAwOUUaO7rPkzgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from numpy import array\n",
    "from statsmodels.stats.power import TTestIndPower\n",
    "# from statsmodels.stats.power import TTestIndPower\n",
    "\n",
    "# parameters for power analysis\n",
    "effect_sizes = array([0.2, 0.5, 0.8])\n",
    "sample_sizes = array(range(5, 100))\n",
    "# calculate power curves from multiple power analyses\n",
    "analysis = TTestIndPower() # or TTestIndPower for independent samples\n",
    "analysis.plot_power(dep_var='nobs', nobs=sample_sizes, effect_size=effect_sizes)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**What does this mean?**\n",
    "\n",
    "* Power is the probability of rejecting the null hypothesis when, in fact, it is false.\n",
    "* Power is the probability of making a correct decision (to reject the null hypothesis) when the null hypothesis is false.\n",
    "* Power is the probability that a test of significance will pick up on an effect that is present.\n",
    "* Power is the probability that a test of significance will detect a deviation from the null hypothesis, should such a deviation exist.\n",
    "* Power is the probability of avoiding a Type II error.\n",
    "* Simply put, power is the probability of not making a Type II error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "<img style=\"float: center; width:300px;\" src=\"stuff/higgins-science-creationism.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Hypothesis testing:\n",
    "* parametric\n",
    "* non-parametric\n",
    "* Bayesian\n",
    "* model-based\n",
    "* ...\n",
    "\n",
    "More later!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The research process\n",
    "\n",
    "#### Pearl's research flow\n",
    "<img style=\"float: center; width:700px;\" src=\"stuff/Pearl-flow.png\">\n",
    "Pearl & Mackenzie, “The book of why”, 2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Further readings\n",
    "\n",
    "* [A review of statistical concepts](http://www.sagepub.com/upm-data/49259_ch_1.pdf)\n",
    "* [Bruno Olshausen's probability primer](http://www.rctn.org/bruno/npb163/probability.pdf)\n",
    "* [Handbook of Biological Statistics - Power analysis](http://www.biostathandbook.com/power.html)\n",
    "* [Introduction to Power analysis in Python](https://towardsdatascience.com/introduction-to-power-analysis-in-python-e7b748dfa26)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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
}