Stats with Scipy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python ver: 3.6.7 (default, Feb 28 2019, 07:28:18) [MSC v.1900 64 bit (AMD64)]\n",
      "Python env: dsml\n",
      "Currrent dir: C:\\Users\\catch\\Documents\\GitHub\\Statistics-in-Python\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/markdown": [
       "First code cell from custom template `dsml.ipynb`. \n",
       "_See Tim Paine's [jupyter_lab templates extention](https://github.com/timkpaine/jupyterlab_templates)._\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sys\n",
    "from pathlib import Path, PurePath as PPath\n",
    "\n",
    "print('Python ver: {}\\nPython env: {}'.format(sys.version, Path(sys.prefix).name))\n",
    "print('Currrent dir: {}\\n'.format(Path.cwd()))\n",
    "\n",
    "def add_to_sys_path(this_path, up=False):\n",
    "    \"\"\"\n",
    "    Prepend this_path to sys.path.\n",
    "    If up=True, path refers to parent folder (1 level up).\n",
    "    \"\"\"\n",
    "    if up:\n",
    "        # NB: Path does not have a str method.\n",
    "        newp = str(PPath(this_path).parent)\n",
    "    else:\n",
    "        newp = str(PPath(this_path)) \n",
    "    \n",
    "    if newp not in sys.path:\n",
    "        sys.path.insert(1, newp)\n",
    "        print('Path added to sys.path: {}'.format(newp))\n",
    "\n",
    "# if notebook inside another folder, eg ./notebooks:\n",
    "nb_folder = 'notebooks'\n",
    "add_to_sys_path(Path.cwd(), Path.cwd().name.startswith(nb_folder))\n",
    "\n",
    "\n",
    "def get_project_dirs(which=['data', 'images'], nb_folder='notebooks'):\n",
    "    dir_lst = []\n",
    "    if Path.cwd().name.startswith(nb_folder):\n",
    "        dir_fn = Path.cwd().parent.joinpath\n",
    "    else:\n",
    "        dir_fn = Path.cwd().joinpath\n",
    "        \n",
    "    for d in which:\n",
    "        DIR = dir_fn(d)\n",
    "        if not DIR.exists():\n",
    "            Path.mkdir(DIR)\n",
    "        dir_lst.append(DIR)\n",
    "    return dir_lst\n",
    "\n",
    "DIR_DATA, DIR_IMG = get_project_dirs()\n",
    "    \n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import stats as sps\n",
    "import pandas as pd\n",
    "#pd.set_option(\"display.max_colwidth\", 200)\n",
    "\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "plt.ion()\n",
    "plt.style.use('seaborn-muted')\n",
    "\n",
    "from pprint import pprint as pp\n",
    "\n",
    "# Filtered dir() for method discovery:\n",
    "def filter_dir(obj, start_with_str='_', exclude=True):\n",
    "    return [d for d in dir(obj) if not d.startswith(start_with_str) == exclude]\n",
    "\n",
    "def get_mdl_pkgs(alib):\n",
    "    import inspect\n",
    "    \"Inspect module hierarchy on two levels ony.\"\n",
    "    for name, mdl in inspect.getmembers(alib, inspect.ismodule):\n",
    "        print('\\n{:>13} : {}'.format(mdl.__name__, filter_dir(mdl)))\n",
    "        for mdl_name, mdl_sub in inspect.getmembers(mdl, inspect.ismodule):\n",
    "            if mdl_sub.__doc__:\n",
    "                print('\\n{:>20} : {}'.format(mdl_name, mdl_sub.__doc__.strip()))\n",
    "                \n",
    "\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\"\n",
    "\n",
    "from IPython.display import HTML, Markdown #, IFrame\n",
    "# for presentations:\n",
    "#display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "\n",
    "Markdown( \"\"\"First code cell from custom template `dsml.ipynb`. \n",
    "_See Tim Paine's [jupyter_lab templates extention](https://github.com/timkpaine/jupyterlab_templates)._\n",
    "\"\"\")\n",
    "\n",
    "def new_section(title='New section'):\n",
    "    style = \"text-align:center;background:#c2d3ef;padding:16px;color:#ffffff;font-size:2em;width:98%\"\n",
    "    return HTML('<div style=\"{}\">{}</div>'.format(style, title))\n",
    "\n",
    "\n",
    "def add_div(div_class, div_start, div_text, output_string=True):\n",
    "    from IPython import get_ipython\n",
    "    from IPython.display import HTML, Markdown\n",
    "    \"\"\"\n",
    "    Behaviour with default `output_string=True`:\n",
    "    The cell is overwritten with the output string, but the cell mode is still in 'code' not 'markdown':\n",
    "    ```\n",
    "    [x]\n",
    "    add_div('alert-warning', 'Tip: ', 'some tip here', output_string=True)\n",
    "    [x]\n",
    "    <div class=\"alert alert-warning\"><b>Tip: </b>some tip here</div>\n",
    "    ```\n",
    "    The only thing to do is change the cell mode to Markdown.\n",
    "    If `output_string=False`, the HTML output is displayed in an output cell.\n",
    "    \"\"\"\n",
    "    accepted = ['alert-info', 'alert-warning', 'alert-danger']\n",
    "    if div_class not in accepted:\n",
    "        return HTML(f\"\"\"<div class=\"alert\"><b>Wrong class:</b> `div_start` is one of {accepted}.\n",
    "                    </div>\"\"\")\n",
    "    div = f\"\"\"<div class=\"alert {div_class}\"><b>{div_start}</b>{div_text}</div>\"\"\"\n",
    "    if output_string:\n",
    "        return get_ipython().set_next_input(div, 'markdown')\n",
    "    else:\n",
    "        return Markdown(div) #HTML(div)\n",
    "\n",
    "# autoreload extension\n",
    "from IPython import get_ipython\n",
    "ipython = get_ipython()\n",
    "\n",
    "if 'autoreload' not in ipython.extension_manager.loaded:\n",
    "    %load_ext autoreload\n",
    "\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some nice colors from http://colorbrewer2.org/\n",
    "COLOR1 = '#7fc97f'\n",
    "COLOR2 = '#beaed4'\n",
    "COLOR3 = '#fdc086'\n",
    "COLOR4 = '#ffff99'\n",
    "COLOR5 = '#386cb0'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---\n",
    "# Computational Statistics with Scipy\n",
    "\n",
    "See Scipy 2015 tutorial\n",
    "https://github.com/AllenDowney/CompStats\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical inference\n",
    "\n",
    "## Quantifiers in order of importance:\n",
    "\n",
    " 1. Effect size: usually a single number, ideally comparable acros stidues\n",
    "\n",
    " 2. Confidence interval/SE: quantifies the precision of the estimate\n",
    "\n",
    " 3. p-value: indicates whether the apparent effect might be due to chance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Helper function to obtain the domain and range of a (scipy.stats [continuous | discrete] distribution) random variable (rv) using its pdf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eval_pdf(rv, num=4, points=100):\n",
    "    \"\"\"\n",
    "    Evaluate the normal (Gaussian) probability density function (PDF) within '`num` standard deviations of the mean.\n",
    "    Input: \n",
    "        rv: numpy rv object\n",
    "        num: the number of standard deviations away from the mean to plot\n",
    "        points: number of points to plot\n",
    "    Output: a pair of NumPy arrays, i.e: domain, rv.pdf(domain)\n",
    "    \"\"\"\n",
    "    mean, std = rv.mean(), rv.std()\n",
    "    xs = np.linspace(mean - num*std, mean + num*std, points)\n",
    "    ys = rv.pdf(xs)\n",
    "    return xs, ys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `scipy.stats.norm` to produce normal distributions given a mean and standard deviation.\n",
    "\n",
    "**Note**: other names for mean and standard deviation:  \n",
    "\n",
    "| Mean | Standard deviation |  \n",
    "|----|----|  \n",
    "| mean | std |\n",
    "| $\\mu$ | $\\sigma$ |   \n",
    "| `loc` | `scale` |  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Given (sample) stats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1, sig1 = 178, 7.7\n",
    "male_height = sps.norm(mu1, sig1)\n",
    "\n",
    "mu2, sig2 = 163, 7.3\n",
    "female_height = sps.norm(mu2, sig2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 'Recover' the true population using the pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29eZQkd3Xn+7m51b4v3bV0d1VLrV4k0VraCJCAERgjYYzGniNbPN7DYM/R0SA44OPlgG3GPPvY4zfDm/cMD6MnDbYsLIzB2Fh49Aw2liyzSKi19aKW1Ht3VVfXvi+53vdHZHVlRGRWZ1VlZmRG/j7n5OmMX/wi4mZ25Td+cX/3d6+oKgaDwWDwLwGvDTAYDAZDcTFCbzAYDD7HCL3BYDD4HCP0BoPB4HOM0BsMBoPPCXltQDY6Ozt1YGDAazMMBoOhYnjhhRcmVLUr276yFPqBgQEOHz7stRkGg8FQMYjI+Vz7jOvGYDAYfI4ReoPBYPA5RugNBoPB55Slj95gMBgA4vE4Q0NDrKyseG1K2VBbW0t/fz/hcDjvY4zQGwyGsmVoaIimpiYGBgYQEa/N8RxVZXJykqGhIQYHB/M+zrhuDAZD2bKyskJHR4cR+TQiQkdHx4afcMyI3lC26MIoOvwCLE9ZDYEg0joIfbcgwYi3xhlKhhF5O5v5PozQG8oOjc6hZ56GsePufdPnYOgnMHAH9NxsRMBgyAPjujGUFbo4gb7w51lF/grxRfTkd9HXnkA1VTrjDIYN8vTTT/P+97/fazOM0BvKB12eQY/8FcSX8jtg7FVL8E3xHINhXYzrxlAWaGzREvnYgn1HQxey820QaUSnz1pum1Ribf/Iy2ioFtl9Z2kNNpSUp/5unSe8AnDnz1+fc9+5c+e46667uOOOO3j22Wc5ePAgH/3oR/m93/s9xsbGePzxxwH41Kc+xfLyMnV1dfz5n/85e/futZ1ncXGRT3ziExw9epREIsHnPvc57rnnnqJ+rlXMiN7gOaqKvv4krMzYd/TchNz6K0j3AaR1J4HBdyKH/iNEmuz9Lj6LTp4uncGGquPUqVN88pOf5MiRI7z22mt87Wtf4wc/+AGf//zn+aM/+iP27dvHM888w0svvcTv//7v89u//duuc/zhH/4h73rXu3j++ed56qmn+M3f/E0WFxdLYr8Z0Ru8Z/w1mDplb+s+gOx5LyL2sYjUtcGb7kNf+UuIL19p15Pfhdb/aKJxDEVhcHCQG2+8EYDrr7+ed7/73YgIN954I+fOnWN2dpZf/uVf5uTJk4gI8XjcdY7vfe97PPHEE3z+858HrNDRCxcusH///qLbb0b0Bk/R+DJ66p/sjc19yN73u0R+FWnoRG64F8iIuInOomefKZ6hhqqmpqbmyvtAIHBlOxAIkEgk+OxnP8udd97JsWPH+M53vpM1zl1V+da3vsXLL7/Myy+/XDKRBzOiN3iMnnkK4hmPrxJArnsfEgiue5w096F9h2D4+bXG4cPotuuRpp4iWWvwivV86OXA7OwsfX19ADz66KNZ+7z3ve/li1/8Il/84hcREV566SVuvvnmkthnRvQGz9DFCbj8ir1x59uQhs68jpfBd0BNS+YZrfh7g6HE/NZv/Raf+cxnuP3220kmk1n7fPaznyUej/OmN72JG264gc9+9rMls0/KMTTt0KFDagqP+J/UiSfs8fL1HdbkayD/B02dPI0e+4atTW7+MNLcVygzS0osmmDi0hyxqCUWgaDQ3t1IY0utx5Z5w4kTJ0rm3qgksn0vIvKCqh7K1t+4bgyeoMtTMPaqrU0G3rEhkQeQjmvQlp0we2Ht3Od/iNz4i4Uws2REV+JcPDnJ8NkpUkn74Os0o7Rva2RgbxctHfUeWWioZIzQGzxBL/wYyBC0+k7o3Juz/3rIrtvRI2tCz9RpdH6kYnz1U2MLHHvuIslE7lW+U6MLTI0usGNPB9dcv82kfjBsCOOjN5QcXZmF0WO2Ntn5ts2LV+sucLhq9PyPNmteSZkYmePIjy+sK/KZXDw5yRsvj5jVwIYNYYTeUHJ0+AXIzFFT1wbdm/fDigiy83Z74+Qb6PL0ps9ZCsaGZzn23EU0ZRftSE2IHdd2sOu6Tprb61zHXTo3zWsvDBuxN+SNcd0YSoqmEjB6xNYmO96aM2Y+b9p3Q+N2WLi8dq2RV5Dd/25r5y0SC7MrnDg8jFOrBw90s+PaDoJB6/tQVabHF3n1+SHisbVojssXZ6lvrmHXdV2lNNtQoeT16xKRu0TkdRE5JSKfzrJfROQL6f1HROSWjH3nROSoiLwsIiaUptqZeN22opVQLXQf2PJpRQTpdcQkX34FTWUPdfOSZDLFq4eHSDlG8ntv7mVgb9cVkQfrc7V3N3Lz2weJ1NrHZWdfHWN+ehmD4WpcVehFJAh8CbgbOAB8UEScv8y7gT3p1/3Alx3771TVm3KF/hiqB730sr1h2w1IMP/al+vSfQCCaysYiS/B5MnCnLuAnDk+xuJc1Na29+Zeegfach7T0FzDLW8fJBxZW0imCq8eHsrbv2/YHF/4whfYv38/H/rQh4py/s997nNX0iIUi3xG9G8GTqnqGVWNAV8HnCnX7gEeU4tngVYRqYyQB0PJ0KVJWxgkgPQUbmWgBCOwzb6CUi+9VLDzF4KpsQWGTk/a2rbtaFlX5Fepa4yw9+ZeW9vSQozTxy7nOMJQCP70T/+UJ5988kqWykokHx99H3AxY3sIuC2PPn3ACFYM3fdERIH/V1UfznYREbkf62mAnTt35mW8obLQEcdovrk/71Ww+SI9N6OXXlxrmDmHLk8hde0Fvc5mUFVOHbWLcm19mOsO5j8m6uptpmdXKyPn1zJ9Dp+dpnew3deLqlL/+l+Kev7AOz+Ttf2BBx7gzJkzfOADH+C+++7j9OnTrjTDjz76KN/+9rdJJpMcO3aMX//1XycWi/HVr36VmpoannzySdrb23nkkUd4+OGHicViXHvttXz1q1+lvt6+LuL06dM8+OCDjI+PU19fzyOPPMK+ffu2/vny6JMt5s053b9en9tV9RYs986DIvKObBdR1YdV9ZCqHurqMhNMfkM15Q6pdPrUC4A0dkOTfdSrl4/l6F1aLl+Ycbls9t/aRyi8fl4fJ9e+aTt1DfYsnWeOj27ZPoObhx56iN7eXp566ikWFxdzphk+duwYX/va1/jJT37C7/zO71BfX89LL73EW9/6Vh577DEAfuEXfoHnn3+eV155hf379/OVr3zFdb3777+fL37xi7zwwgt8/vOf52Mf+1hBPkc+I/ohYEfGdj9wKd8+qrr675iI/B2WK8ikGaw2ps/ZK0cFa6Br6yOVbEjPTeh8xp/o+KvowNs9XWSUTKY4e2LM1rZtRwutnQ0bPlcoFOSaG7Zx7Lm1h+jJ0QWmJxZp28T5DPmRK80wwJ133klTUxNNTU20tLTwcz/3cwDceOONHDliRZkdO3aM3/3d32VmZoaFhQXe+9732s6/sLDAj370I+69994rbdGofWCwWfIR+ueBPSIyCAwD9wH/i6PPE8DHReTrWG6dWVUdEZEGIKCq8+n3PwP8fkEsN1QU6kh3QNfeDac7yJuuvXDyu6DpiJvlaSvs0sOVskOnJ4kur1XGkoAwuL970+fr7Gmipb2e2am1m+fpY6Pc+s5Bs2q2SKymGXZWjnruueeumsYY4CMf+Qjf/va3OXjwII8++ihPP/207TypVIrW1lZeftnh4iwAV/2lqWpCRD4OfBcIAn+mqsdF5IH0/oeAJ4H3AaeAJeCj6cO3AX+X/sMLAV9T1X8s+KcwlDWaSsDEG7Y26dp6SGUuJFSLtu+2Rdzo2KuepURIxJNceGPC1ta/u93lftkIIsI1N2zjxWfOXmmbn15m4tI8XX3Nmz5vuZLLh15KtppmeH5+np6eHuLxOI8//viVtMarNDc3Mzg4yDe/+U3uvfdeVJUjR45w8ODBLdueVxy9qj6pqtep6jWq+ofptofSIk862ubB9P4bVfVwuv2Mqh5Mv65fPdZQZUydhmTGI2i4Htp2FfWS4ozNHzvh2UrSS+emScTXQiBD4QC7rtv6JHRLRz2dPfayihdOTuTobdgqW00z/Ad/8AfcdtttvOc978k5wfr444/zla98hYMHD3L99dfz93//94Uw3aQpNhSf1Kt/Z5ULXKX3VgJ7fqao19RkHP3Rn0BqraSbHPwQ0lraiK5USnn2e2/Y3DYD+7sY3Ld5t00mi3Mr/OT79nq5t7xj0DdZLk2a4uxsNE2xyXVjKCqaiMKkvR6sa7RdBCQYhs49dluc8wQlYHx4zibygaDQN1i4UM+G5lratzXa2i6emszR21CtGKE3FJep05BaEzpqWlyZJouF64Yy8boV5lkiVJWLp+yulO07W4nUFHYSese1Hbbt8UtzLC/GCnoNQ2VjhN5QVNSZgqBrX+miQtp2u1MizDkjg4vH7OQS8zP2ItE7runI0XvztHU10NBcY2tzrr6tZMrRvewlm/k+jNAbioamkjBp9x9L53Ulu74EgtBxjd0mR/RPMRk6PWXb7tjeRH1TTY7em0dE2LnHPrk7cm6GRLz8ErptlNraWiYnJ43Yp1FVJicnqa3d2Cpok6bYUDxmzrujbZp7c/cvAtJxnd03P/kGuvvOoj9VxKIJJkbmbG1OF0sh6e5v5vTxUWIrlpssmUwxNjyXVw6dcqa/v5+hoSHGx8e9NqVsqK2tpb+/f0PHGKE3FA2X26Zjz9bzzm+U9t0ggbVCJ8vTsDQJBc6x4+TyhRlbrvn6xgitncWLhAkEAvTsauX862tzAiPnpyte6MPhMIODg16bUfEY142hKKgqTNiFXhxRMKVAQjXQOmBvLLL7RlVtSccAena1Ff0pomenXdTnppZduXUM1YkRekNxWLgMsfm17UDYLbglwnmD0cniCv3c9DJL82sCK2JF2xSbuixPDSPny7ucoqE0GKE3FAXXpGf77sIVGNkoHY4nifkRNDqfvW8BGDlnF9eO7U2u6lDFomeXfVR/+cIMqZQpTFLtGKE3FIcpR7SNU2xLiNQ0uROaTZ0pyrUSiSRjw/ZJWKf4FpOuvmZC4bWfdTyWZHJkoWTXN5QnRugNBUejC7DgyI/uCHMsNdJ+rW1bHTeiQjE5Mm8r7RepDblWrhaTYDBAd3+Lre3yxZkcvQ3VghF6Q+GZdoyWm3qQsMe5V9p327enzxWlcPjo0Kxte/uOFgKB0qYN7nHMB0yNLvgipt6weYzQGwqOa7Tc7u1oHrBcN5k3m2QU5oYLeol4LMHU6KKtzTm6LgVNbXXUNqzNh6RSyvil4s1JGMofI/SGgqKasqpJZSDO0bQHiAi02eOxC+2+Gb80b1vBWd8Y8aSOq4iwrc9+gxkbns3R21ANGKE3FJa5YUhk5HcJ1Xla2SkTcT5ZFHhC1um26e5v8azaU/cOu9BPjy0QiyZy9Db4HSP0hoKiTvFs31361bC5aHessFwcK1iYZXQlzsy43W2zzQO3zSqNzbU0ZOTVUbWyWhqqkzL5BRp8gzOssgzcNqtIuB6aHLl2CjSqH3eEVDa21BYlgdlGcI7qRy8a9021YoTeUDA0tugOqywjoQe3+0adEUKbZGzILvRejuav2ODw089OLhFdjufobfAzRugNhWP6rH27abv3YZVOnO6b6XNbLkYSXYkzO7Vka+vu975Ad11jhKZW+2TwxIiJvqlGjNAbCoY6om1oLcOsg009EMoQv8SK+ylkgzjFs6m1ltr6yJbOWSi6eu03nPER46evRozQGwqCqrrDKtvKT+hFAtC6y97ovEFtkAlHjHpnr/ej+VWctsyMLxKPmcVT1YYRekNhWJ5yZKsMQUtpasNuFGkbsG2r0+W0AeKxJNPj9lwyXT3lI/QNTTW2SWFVmLxs3DfVhhF6Q2FwimXLDiRQpnVtnE8as0NocnOTlJOj864CI876rV7T1dNk2zZhltWHEXpDQXCOisvRbbOK1LVBbUY+GE3C7MVNnWvCIZrl5LZZxWnT1NiCLfGawf8YoTdsGU0lYeaCvbGMhR4Al/vm3IZPkUymmBx1uG3KUOibWmupqcvIfZNUpsZM6uJqwgi9YevMj0AytrYdroeGLu/syQOnn97lesqD6bEFUsk1v01NXcgVzlgOiAhdvXb3jfNJxOBv8hJ6EblLRF4XkVMi8uks+0VEvpDef0REbnHsD4rISyLyD4Uy3FBGOEfDbQOe5XjJG2fkzeIYGlvK3jcHE5fto+LOnuay/dydjgniydEFWwI2g7+5qtCLSBD4EnA3cAD4oIgccHS7G9iTft0PfNmx/5PAiS1bayhLdOa8bVs8qg27ESRcD43b7Y2z57N3zoKquqJXOrY35ejtPS0d9a7KU3NTyx5aZCgl+Yzo3wycUtUzqhoDvg7c4+hzD/CYWjwLtIpID4CI9AM/C/yPAtptKBM0GXfndXe6RcoVx6hep/MX+oWZFWIra9kgg8EAbZ1ltgo4g0BAXJWuJkyYZdWQj9D3AZkhCUPptnz7/N/AbwHrTvOLyP0iclhEDo+Pj+dhlqEsmBu2olZWqW1Far3P85IP4nTfzOQv9E6RbOtuIBAs7ykv5xOHiaevHvL5y8zmdHQ697L2EZH3A2Oq+sLVLqKqD6vqIVU91NVV3hN5hjWcbhuX77ucaemHzBTKy1N5py12imRnGbttVulwjOgX56IsL8Zy9Db4iXyEfgjYkbHdD1zKs8/twAdE5ByWy+ddIvKXm7bWUH440x5UkNBLqMZdFCWPUX10Oc78zIqtrZz986uEIyFaOuzuJTOqrw7yEfrngT0iMigiEeA+4AlHnyeAD6ejb94CzKrqiKp+RlX7VXUgfdy/qOr/WsgPYPAOTUSt0MpMKkjoAbefPg+hd8bON7XVEakt01XADpxPHpOXTTx9NXBVoVfVBPBx4LtYkTPfUNXjIvKAiDyQ7vYkcAY4BTwCfKxI9hrKidmL2Lx49R1ITWPO7uWI6wkkjwnZyRGn26ZyPrPzyWN6YpFEwiQ58zt5DUNU9UksMc9seyjjvQIPXuUcTwNPb9hCQ9lS0f75VZr7QIJrE8rRWXR5Bqlrzdo9lUwx7SgZWAlum1XqmyLUNUSu+OY1pUyPLZblil5D4SjvMAFDeeOKn688oZdg2BL7TGbO5ew/M7lEMrkWQBapDdHYUn6rYXMh4g6znBo17hu/Y4TesCk0vuwu2NG60xtjtoi0Of30F3L0dItix7bGsl0Nmwtn9I1ZJet/jNAbNsfskH27obv8ygbmS4vjSWT2Qk7hc07EOkfHlUBrVwOBwNrNKbocZ2k+6qFFhmJjhN6wKdz++coczQPQ3GMVSlklOg8rM65uy4sxmyCKQFtX5Ql9MBigtavB1ua8gRn8hRF6w+aYtbs3pIKFXgKhLH56t/vG6bZpbq8nHAkW07Si4XTfGD+9vzFCb9gwWf3zLTuyd64QnDeqbPH0zlGvUywrCafLaWZyyYRZ+hgj9IaN46zGVMn++VWcEUMOP30qmWKmgsMqndQ31lDXELmyrSl1fT6DfzBCb9gwrqiUCnbbXKFpfT99trDKcqsNu1Gco3rjp/cvRugNG8dH/vlVLD99v70x44bmh7BKJ9n89CbM0p8YoTdsCD/651dZz0/vrLFaiWGVTlo77WGWK0txlhdMNks/YoTesDFc/vmuyvfPr+J8Mkn76aPLcRbn7HHmlRhW6SQYCtDiKJZiiob7EyP0hg3h9s9XXtqDnOTw0zvFr7mtrmLDKp20dxs/fTVghN6wMZz++ZbK98+vkiue3umf94PbZhWnn35mYpFUct1icIYKxAi9IW80sZIlv40//POrOG9cOnOBqTF72KGfhL6+qYaajFz6qaQyM7nkoUWGYmCE3pA/Tv98fad//POrOPz0qenzJOJrC4lC4SDNbXWltqpomGyW1YERekPe6IxD6H0QVumiudfKT58mEJ+nJrA2wm3rbqj4sEonTj+9mZD1H0boDfkz68w/7z+ht/z0vba21tDElfcd3f5x26zS1m1PcLY4FyW6HPfIGkMxMEJvyAurPqwzft5/Qg+4nlRawpNX3rf5yD+/SjgScrmjzKjeXxihN+TH7BCu+rCRhpzdKxnnhGxryBL6hqYaauvCXphUdFx++jGT98ZPGKE35IU6wip9O5qHdB3ZtZ9GXXCJmsAybT5026zi9NNPj5l0CH7CCL0hP2b8l98mFxIMQ5PdT98SmvBVWKWTprY6QuE1OYjHkizMrHhokaGQGKE3XBVNxmB+xN7o5xE9kKi3L5xqDU/R2uGzUNIMAgGhzVF1yvjp/YMResPVcfrn69qRGv+ObgFmEx227faaSYIhf/9cnK4pI/T+wd9/uYaC4Mv881fh8mwjKV2Ll69hAY3Oe2hR8XH66Wcnl03VKZ9ghN5wdXyc3yYbqVSK6YkY84lW+w7nhLTPqGuI2KtOqTIzbtIh+AEj9IZ1yeqf9/mIfnZymWQy5XLfuJ5sfIg7zNK4b/yAEXrD+swNg2ZkM6xtRWoqt1ZqPqyK20yi076jGoTe+Ol9SV5CLyJ3icjrInJKRD6dZb+IyBfS+4+IyC3p9loR+YmIvCIix0Xkfy/0BzAUF1/nn8/BqrjNxdvRDD89y1No1N/C19pZT2Yqn+WFGMuLpupUpXNVoReRIPAl4G7gAPBBETng6HY3sCf9uh/4cro9CrxLVQ8CNwF3ichbCmS7oRQ4MlaKz9ISO4lFE1fix5OEmE+22Dv43E8fCgdpcYSRTptRfcWTz4j+zcApVT2jqjHg68A9jj73AI+pxbNAq4j0pLdX/0rC6ZdZblchaDIOc5fsjT6fiHWK2lJgm227Gvz0JszSf+Qj9H1A5rBuKN2WVx8RCYrIy8AY8E+q+ly2i4jI/SJyWEQOj4+P52u/oZjMDYNmhNfVtiC1Lbn7+wBXjpcsdWT9jjND5/T4IqmUGZ9VMvkIfbbk287/9Zx9VDWpqjcB/cCbReSGbBdR1YdV9ZCqHurq6srDLEOxqTb/vKq6Rq91Pbux/XkvTaIxf49wG1trbTVxE/EU89PLHlpk2Cr5CP0QkOmY7QcubbSPqs4ATwN3bdhKgzc44+d9Hla5OBcltpK4sh0MBmjpaoWm7faOPnffiIhx3/iMfIT+eWCPiAyKSAS4D3jC0ecJ4MPp6Ju3ALOqOiIiXSLSCiAidcBPA68V0H5DkbD888P2Rp/7551i1tpZTyAYcD3JVIOfvr3b5L3xE6GrdVDVhIh8HPguEAT+TFWPi8gD6f0PAU8C7wNOAUvAR9OH9wB/kY7cCQDfUNV/KPzHMBScbPHzvvfP28VsdfGQtO5ELz67tmPGXmnLjzjj6eemlonHkjaXjqFyuKrQA6jqk1hintn2UMZ7BR7MctwR4OYt2mjwAHWKmc/dNslEitkJ+3L/K2LX3I/lp09PTaXj6f2c2K2mLkxDcw2Lc9ErbdPjC3T3+ftm71fMylhDdlz55/09ETszYY8sqa0PU9do5X2RUA009dgPqILoG9cq2VHjvqlUjNAbXFj5baorft7ltuluRDKXiFajnz5L3htTdaoyMUJvcDObzT/f7J09JcA5WnWKnCviqAr89C0d9QSCaze76HKCpfnoOkcYyhUj9AYXbv+8v902y4sxlhbW8rmI4Kq2REu/rY6s5af3d376YDBAa6f9e5g07puKxAi9wY1D6P3un3e6bZrb6wmF7dElEoy4/fRVMKrvMGmLfYEReoMNTaxkyT/vc6G/itvmCi4/vf+F3lV1amKJZCKVo7ehXDFCb7AzexFbhov6Dl+HEaZSyvS4Pb+NM9fLKtI6YG+YPuf7ycm6xgi19eEr26mUMjOxuM4RhnLECL3Bhk5Xl39+bso+Qg1HgjS21mbv3NIHkuHSic7BykyRLfQWETFVp3yAEXqDHZd/fsAbO0qEc3KxfZsjrDIDCYSsSdlMqtB9YyZkKw8j9IYraHwJFsfsjT4vNOLyz+dw26zinJiuBj99W1eDqTpV4RihN6zhXATUuA0J12fv6wOiy3EWZldsbTknYldxurKmz/veT29VnXKGWfo7tNRvGKE3XKHa4uedo/mmtjoiNVdJ/9TUA4G1yUnii7A0WQTrygtnmOXkZeO+qSSM0BvWmD5n2/R7/LzT1+wUs2xIIOh2Z82cK6BV5Un7dvt3MzOxSDJpwiwrBSP0BgB0ZQ6WpzJaBFr865+3wio3LvTgnqBWxw3SjzQ01VBTlxFmmVRmHNk+DeWLEXqDhXNU2txnZW30KXNTSyTi9rDKpra6/A5uG7Bvz5xH1d+jWxFxr5I1fvqKwQi9AQCdPmvbFqeY+YyNhFW6aOiGzEnqZAzmnNU1/Ydzotr46SsHI/QGK2rE6X5oG/TEllLhnIjN120D1ujWNaqvAveNFWa5djO0ksGZbJaVgBF6AyyOQzzD35otgZePWMkSVukshn013H76s9k7+ohQOEhLpz3c1hQjqQyM0Bvco9HWnVZ0iU+ZvGz3LTe35xFW6cQ5op+/hCb8P7p1PvlMjBg/fSVghN7g9s/7Pe2BQ5w6tzdt+BxS2wJ17WsNmqqK8oLO72pmYolEPOmRNYZ8MUJf5Wgqkc5YmYGP/fPJRMqdrXITQg+4RvXVEGZZ31RDXUPkyraqMjVmslmWO0boq53ZYUjF17YjjVDf4Z09RWZ63F4EvKYuTEPz5sJIxXlDrAI/PbhvjE5XmKH8MEJf5ej0GXtD22D+YYYViFOUOrdvIKzSSetOIOPYpUl0ZXbzxlUInT3OMMt53+f7qXSM0Fc7U3ahl/bdHhlSfFTVJfQdPZt02wASqoXmPntjFYzqWzoaCIbWpCMeSzI3veyhRYarYYS+itHovDstsY8XSi3MrhBdSVzZDgTFVfx6ozhvjOq4cfqRQMC9Sta4b8obI/TVjHP02dTr67TEEw4xau9uJBjc4k/A+QQ0fQ5N+T8KxemnN2GW5Y0R+irGNfps92+0DbjFaNPRNpk0bnekQ4jCfPWlQ1ici5piJGVMXkIvIneJyOsickpEPp1lv4jIF9L7j4jILen2HSLylIicEJHjIvLJQn8Aw+ZQTblG9NLmX//8ylKMhRn7atjNxM87sSKL/wEAACAASURBVNIh2G+Q1eC+idSEaOmwP/2ZUX35clWhF5Eg8CXgbuAA8EEROeDodjewJ/26H/hyuj0B/Lqq7gfeAjyY5ViDF8yPQCJD+EK10NzrnT1FxilCLe31RGo3uBo2B64J7CoQeoCuXqf7Zs4jSwxXI58R/ZuBU6p6RlVjwNeBexx97gEeU4tngVYR6VHVEVV9EUBV54ETgCNMweAJTjFqG0TEv548p9B3biHaxoUznn7hMhrz/yKizp5m2/bMxBLxaCJHb4OX5PPL7gMyl04O4Rbrq/YRkQHgZuC5jRppKDw6ddq27Vr84yPisSQzE3bh7ewtnNBLpMHy1Wfi+H79SF1DxLXYbMKkLi5L8hH6bKtJnKsj1u0jIo3At4BPqWrW5zsRuV9EDovI4fHx8TzMMmwWjS1arptMOq7xxpgSYC3oWduub6qhvrHARVXa7d+f80bqV5yjeuO+KU/yEfohILOmXD/gDCvI2UdEwlgi/7iq/m2ui6jqw6p6SFUPdXV15WO7YbM4RahpOxLZWJreSsLptukqpNsmjXRca2+YOlsVYZZOP/3U2ALJhL+rbVUi+Qj988AeERkUkQhwH/CEo88TwIfT0TdvAWZVdUSsteVfAU6o6n8vqOWGTaOTp2zb0n5tjp6VTzKZcuVM7+xtztF7CzT1uMMsncnifEhjS62rluzUmHHflBtXFXpVTQAfB76LNZn6DVU9LiIPiMgD6W5PAmeAU8AjwMfS7bcD/xvwLhF5Of16X6E/hCF/NJVwL5RyjkZ9xNToAsnk2ggzUhuiqbW24NcREdf36Lyh+hERcU1sj18y7ptyI6/4MlV9EkvMM9seynivwINZjvsB2f33Bq+YuWjVOF0l0uieSPQR48N20enqbS5a0jZpvxa9fGStYeoU8NNFuVY50dXXzPCZqSvbEyPzpJIpAltddWwoGOZ/osrQKccos/0a32arTCZTrrQH3X1FcNus0jYAmSGqy9Po0mTxrlcmtHbU2yp0JRMpk6O+zDBCX0WoKjj98z5220w7JgazreYsJBKqgZad9sZqcd84JmXHh/2frrmSMEJfTSyOw8rM2rYEfZ2tcszptukrnttmFeeNUydPFvV65UJ3X4tte+Ky5b4xlAdG6KuJidft220DSDCSvW+Fk0qmXKlzu4oRbeOkc499e/ZiVaySbemoJxxZKyifiLtLNhq8wwh9FaETb9i2pfM6jywpPlPjiyTiayPKcE2Q1s7ip2CW2lZo3GZvrIJRfSAgrhup84nK4B1G6KsEXZ5xFBkR6NiTs3+lMzZk9xEXM9rGifMGqs4nKZ/S1edeJWvcN+WBEfpqwTGap6XfytHiQ5KJlGs1bFGjbZx07rVvT59HE9HSXd8jWjsbXO6byVGzeKocMEJfJThHlX5220xcnrdH29SGtlwycEPUd0Jd29q2JqsiyVkgIHT32ydlRy+a6JtywAh9FaCxRZgbsjf6WOjHHOKyrb+lpGsFrFWyTvfNGzl6+4ttDqGfvDxPIu7/nD/ljhH6asDpI27cZk0a+pB4LOFyFzjFpxS4npgmT6HJeMntKDXN7XXU1mfkvkkp45dM5SmvMUJfBejYCdu2OH3IPmJ8eM5aGJamvjFCYxFy21yV5j4rvcQqqXhVuG9ExHVjHR2aydHbUCqM0Pscjc7D7AV7Y9d+b4wpAaOOaJvuHaV126wiIq7v2XnD9SvbdtiFfnpskdiKqTzlJUbo/c74a/btxu1Ifbs3thSZ5cUYMxNLtjYv3DarSLfjhjp1qiqibxqaa2lssT9FmVG9txih9zk67nDbOMXHR1y+YBeTpra6wleS2ghNvVCTcaNJJaoi9w24b7Aj52dsLjVDaTFC72N0ZRbmhu2NPnXbqKpL6Ht2eTvhLCLguLHq+KseWVNatu1sIdNjtjgXZX5mxTuDqhwj9H7GMZqnuQ+p9c6VUUxmJhZZWVqLagkExJVoywvEeWOdOoPGl70xpoTU1IZp32YvT3n5/LRH1hiM0PsYHT1u23aJjo8YOW8fzXf1NttWaXpG4zbH4qmUO9zVp/TsbLNtjw7N2qp9GUqHEXqfoguj7tw2Xfs8s6eYJOJJV/m67R67bVax3DcHbG16+ahH1pSWjp5GV0oEZ2oKQ2kwQu9TXGLSNoDUNGXvXOGMDc2SSq5N9NXUhWnrKp88PrLtBnvD3BC6PJW9s48IBAKuUMsR477xBCP0PkRTSRhzuG22v8kja4qLqjJ81i4ePbtay6o8otS1Q3O/rU0vH/PImtLSs8vuvpkeW2Rpwf8hpuWGEXo/Mn0G4hnx5MEa36YknpteZmHWHs3hdbRNNmT7jfaG0WNVEW7Y2FJLU1udre3SWTOqLzVG6H2Iy23TvR8JhrN3rnCGz9hdIB3bm6itL8OqWV37ILBWQJvorHvFsk/pG7Qv0Bs5P2MmZUuMEXqfofEldwHwbTfm6F3ZxKIJVxWjvt1tOXp7i4RqXRlD9fIRj6wpLd39zYTCmZOySVdhGENxMULvNy4ftfKfr1LXZiXY8iEj56fR1Jr7o7YhTHt34zpHeIvrhjv+WlXE1AeDAXoG7O4057yKobgYofcRqoqOvGRrk+0Hy2pislCoqsvX2zfYXt6ftW3AnRJhtDpCLfsG7O6b+ell5qb9f5MrF4zQ+4npc7CcIX4ShO0HPTOnmIxfmnethO3ZWX6TsJmIBJDem2xteumlqpiUrWuMuFbKXjw14ZE11YcReh/hHM3TtReJ1HtjTBFRVS6etItEd38L4ZpQjiPKiO0HQTJ+dstTMHPeO3tKSP9u+6h+fHiO5cWYR9ZUF0bofYJG510FwKXnFo+sKS6zk0uux/6dezo8smZjSKTBVTxcL72Uo7e/aN/WSH3TWjZRVRg6PemhRdVDXkIvIneJyOsickpEPp1lv4jIF9L7j4jILRn7/kxExkSkOlaIeIQ1ms9wAdR3Qkt/zv6VzMVTdnFo39ZIQ7MHVaQ2ifTebG+YfAONzmXv7CNExHVDHjk3QzxmasoWm6sKvYgEgS8BdwMHgA+KyAFHt7uBPenX/cCXM/Y9CtxVCGMN2dFkHIZftLVJ7y3lPTG5SZbmo658KTv3dHpkzSZp2Qn1GYKnKXT4Be/sKSHb+luIZLjYkskUl876Px2E1+Qzon8zcEpVz6hqDPg6cI+jzz3AY2rxLNAqIj0AqvoMYP4ni8noUUhkuDJCteDMr+ITLjh8802ttbR2VtY8hIggfT9lb7z0UlVUnwoEA/RfY/fVD52eIpkwC6iKST5C3wdczNgeSrdttM+6iMj9InJYRA6Pj49v5NCqRjWFDv3E3thzMxLysLJSkVheiLmKi+zY01mZTy7bboBwRmqAZBQuv+KdPSWkd7CdYHBNemLRBMNmVF9U8hH6bL8iZzxYPn3WRVUfVtVDqnqoq6trI4dWN5MnHSGVAaTvVu/sKSLnXh8nMxKxrjFCV2+zdwZtAQmGkV77/5MOPW8lpPM54UjQtYL5wskJM6ovIvkI/RCwI2O7H7i0iT6GAqOq6IVn7Y3d1/syHfHSQtQ1mh/Y10UgUIGj+VV6b3Xkv5lzVwXzKTv2dNpG9fFo0ozqi0g+Qv88sEdEBkUkAtwHPOHo8wTw4XT0zVuAWVUdKbCtBifTZ2Hefj+VHbd5ZExxOfea3Z1X3xhxFaCuNCRSD46slnr+h6j6f2QbqQnR5/DVX3jDjOqLxVWFXlUTwMeB7wIngG+o6nEReUBEHkh3exI4A5wCHgE+tnq8iPwV8GNgr4gMicivFvgzVCWqip77N3tj+zVIg//cXgtzK4xetCfBGtjXXZm+eQfSfxs2z+fylKuWgF/ZeW0HwVDGqD6WdIXOGgpDXksJVfVJLDHPbHso470CD+Y49oNbMdCQg6nT7tH8wNu9saXInD46atuub6qhu78yffNOpK4N3X4jZGSy1PM/tFxw4u/1jOGaEP272zn/xlok1YU3JugZaKWm1p9ptb3C339JPkVV0fOO0XzHHqSpxxuDisjk6DxTYwu2tt0H/DGaX0V2vs2RFmEaRqtjfeGOPZ22urLJZIqzr46tc4RhMxihr0QmXof5y7Ym2XWHR8YUj1RKOeUYzbd21tPZ46/JZqlrA0cKYz33AzSV8Mii0hGOBBnYZ3c3jpyfcVUNM2wNI/QVhqYS6Jmn7I2d1yFN270xqIhcOjfF0rx9EdG1N2731Wh+FdnlGNVHZ8G5PsKn9A62U99orwp28sjlqsjqWSqM0FcaQ4dhJTPMUJCBd3plTdGILsc5c9z+CL99ZytNrXU5jqhspLYVeu1J6PTCj9HoQo4j/EMgIFxzo32gMjOx6JqAN2weI/QVhMYW0Qs/tDf23YI0VFiulzw4eWTEFmoXDAXYfaDbQ4uKj+y6w0pfsUoyhp77V+8MKiEd2xpp626wtZ06eplY1P/uq1JghL6C0DNPQTIjf3eo1pe++fFLc4xfsicu232gm5o6f0diSLgOV+TU5SPo3LAn9pQSEeG6gz22BXDxWJJTRy+vc5QhX4zQVwg6dcZVdk523YGEKyuh19WIxxK88Yp9rV1TWx19jqIVvqXnZntmS0Bff7IqUiPUN9a4JmZHL84yeXk+xxGGfDFCXwFoMoa+8Y/2xvpOl0+30lFVXn/pErGVtcd1Edh3c68vJ2CzIYEgcu177I1LE3DhR94YVGJ27OmkodmekO/Ei8O2vwnDxjFCXwHo2X+1ojAykL3vQwLBHEdUJiPnpl0umx17OmlsqZyiIoVA2gbd4ZYXfoQu+j+rayAg6Rv7Wls8muTEi8MmCmcLGKEvc3TqNAwftjf2/RTSvKEs0GXP4lyUkw5/bFNrLYP7/ZfSIR/kmndDOGNyUlPoib+3isz4nOb2enY5XDhTowum7OAWMEJfxmh0Hn3tH+yNta3I4Du8MahIxGNJjj53gVRybcQWDAY4cKifQKA6/0QlXIfs+Rl74+I4evr73hhUYgb2dtHSYZ9/On1slOnxRY8sqmyq81dUAaim0Ne+A/GljFZB9v4sEozkPK7SSKWUV58fYnkhZmvfc3C7rZB0VdK5F7odVTtHXkLH/J/KWEQ4cKiPUHhNolTh2HMXWV6MrXOkIRtG6MsUPf0vMHPe1iYDb0dad3pjUJE4feyyK5fN9h0tbN/Z6pFF5YOIIHvugjp7kQ59/X+iC6M5jvIPtfUR9t1id1Em4kmO/vgCibj/o5AKiRH6MkQvvQjDz9sbWwdg51s9sadYnHt9nKHT9mITzW11XFdFUTZXQ0I1yP5/D5Ix8Z6Ko0e/iUb9H3bY1dvMoGOh3OJ8lCM/vmBy128AI/Rlhk6eQk9+z94YaUT2/5yv0tYOnZ50ZSmsqQ1xw1t22CoPGUCatiPX/rS9MTaPHvubqigovuu6Tlda6tnJpfS8jhH7fDC/qDJCp86gx/8WW7ndQBi54V4k0uiZXYVm6MwkJ4/YI2yCoQA3vnWnyUOeA+m9BfoO2RsXLqNH/9r3Yi8i7Lu5j+Y2e56j6bFFjj130Yzs88AIfZmgU2fR498CzfQ9CrL/Ht9kplRVzrw6yslX7CIfCApveutO3yYsKxRyzbuh/Vp749wwevQbvhf7YCjAm962y7WmYnJ0gZd/eI54zCyoWg8j9GWAjh1Hj30THPnHZc/PIJ17PLKqsCSTKV5/6RLnX5+wtUtAuOG2nbR2NuQ40rCKSAA5cA809dp3zA2hr3zN9z77cCTIwdt3uaKx5qaWefGZc67ILcMaRug9xKoU9UP0xBOOkTzIte+xHtd9wPJijJeeOcvI+RlbeyAo3HjbDjq2+cctVWwkGEHe9EvgrCa2cBl96S/QBX9XZ4rUhLjpjl2uNAlL81EOP32a8UtzHllW3hih9wiNL6HH/wY994xrn1zz04jTH1uBqCpjQ7McfuoM8zP2ikGhcJCb7higY7u/qkWVAgnVIm+6zy320XlL7Ede9nW6gJraMDe/fdC1oCoRT3HsuYuuFNcGI/SeoFNn0MNfgclT9h0SQPb9HNL/U94YVkCiK3GOPXeR488PuWKeaxvC3PLOQVra/ZV5s5RYYv9BaL/GviOVQN/4/9BXv43G/LuKdNWN093nLhI/dHqK5//ltFlFm4GU453/0KFDevjw4at3rDA0OmctYR9/zb0zVItc/x8qfkFUMpHiwqkJLp6czDqq6tzexL5b+2wFoQ2bRzWFnvpnuPSCe2eoFhl8J/Tc5KvQ3ExUlaHTk5w+Nko2KevqbWb39d3UN/p/lbWIvKCqWV0BRuhLgMaW0Is/hksvuiZcAWjqRfZ/wCoSXaEk4kkunZ3m4qnJrFWBRITBA93s3NNhFkMVAb18xFp/kcqS9Ky+Exm4Azr3+fa7n5lY5NXDw0SX3Z9fxCpDuWNPJw0+TqthhN4jdHEcHX4BRo9l/wEC7HgrMvD2ik05vDC7wsj5aUbOz+T0iza31bHvll4amqsr3XCp0aVJa2J/IUdVpvoOpPdW2HYDEvKf4CXiSc68OsbwmamcfTq2N9Kzq42O7U22alZ+wAh9CdHoHIy/ho69CvMjuTs2bUf23IU4J9TKHFVlcS7KxMg845fmWJhdydk3HAkysK+Lvt3tvh1JlhuaSsLwC+j5f7OXncwkGIGOPUj3AWgbQAKh0hpZZGYmFjl19LIrACCTcE2Qrt5mOnuaaeusJ+CD1dhG6IuIxpasOObZizB9Fq5WHCLSaNV57TlYEX7TVEpZmo8yN73MzMQiM+OLRK9S7ScQFPqv6WDXdZ2EwpX5pFLpaHTOKlgzehzbSmsngTC07kTaBqC5Hxq3VezTZSZWxNccZ06MsrK4fg7/QFBo7Wigtaue5vZ6mlpqK/Lv1gj9FlFNQWwBVuZgeQpdnrYEfWHMVfkpJ5FGZMdt0HMzEiyvZf6qSjyaJLocZ3kpxvJCjKWFGItzKyzOR2154tcjHAnSf007fYPthGv8NUqsVHRpEj3/Axg7wbqCv4oEoaELGruR+g6oa4faVqhtRkKV53pLpZTxS3NcPDmx7gjfSV1jhIamGhqaaqhrjFDXEKG2PkykNlS2NRK2LPQichfwJ0AQ+B+q+seO/ZLe/z5gCfiIqr6Yz7HZ2KjQq6YgsWIlrFYFUun3KfsrlbQWJqWS1qRoKgHJuOU/T8asZeTJKMRXIL5s5YKPLUJ80Tp+M7TsRPputR6VNzFSiseSVky0gpL+VzX98dT2PpX5SqasfxMpkskUyUSKRCJFMp4kHk+RiCWJxRLEo0liK4ktxV23dzeyfVcrnT1NJiFZmaLReXTkJRh5xRq0bIZgBCINVuWrcD2E66zInlAEgjUQDFuvQBgCIeslAQgErRuIBCAQAALWDKmk/8VqL+aNRFWZm17m8vkZxoZnScQ3H2cfjgSJ1IQI1wQJR0KEwgFCkSChUIBgKEgwKARCAYKBAIGgEAgKEggQCAgipP+1Xoi1OlzAei/W+2A4uOE5hC0JvYgEgTeA9wBDwPPAB1X11Yw+7wM+gSX0twF/oqq35XNsNjYs9LND6Mtfzbt/cRFo7kU690L3fqTGHee7EZ77p5MsldnS7kBQaO1soKuniY6eJpOIrIJQTcHsRat4yeQpiJVJ2oTGbQRu/ZWSXCqZTDEzvsj4yDyTl+fLsvD4m999zYaDF9YT+nyer98MnFLVM+mTfR24B8gU63uAx9S6azwrIq0i0gMM5HHs1vF0ok+gcRu09CHNO6zJrXABk3OVwSRmpCZEU1stzW11tHY20NxW54vJq2pEJACtu5DWXdZT3NIETJ9FZ4dhbmjzo/2tW1ayKwWDATq2N9GxvQlVZWkhxsz4IrNTS8xPL5fJwKqw30c+Qt8HXMzYHsIatV+tT1+exwIgIvcD9wPs3LnBRUOlmNQM1UFts+WvrGtH6tuhoRsaOosatVAqnQ+FA0Rqw9TVh6ltiFDfGKG+qYaG5hoiNSETNeNDRMTyxzd0If1Wm8YWrbmnpXF0eQqWpq15qJU5Vz6mAhtTvHOve1m54ovv290OWGGaq3NUSwsxVhZjLC/GiC4nsq4RKY5dhT1fPgqV7ZJOf0+uPvkcazWqPgw8DJbrJg+7Mq4egFBtWvDF4f+TtI8wYPkJV/2FgRAEQ5Y/Me1XlFANhGqsc4XqLB9kpBEiDZ5NoIbCQWsVqYAga3681X/Tfj8JCIH0diAoBAJCMBggELT8hKv+w1A4QCgctM5bY/kaIzUhgiEzQjeARBqgfRDaB20/XlWFxLI1ZxVbtN7HlyERRZMxSETTc13x9PxXPD0XljEvdmW+zDGPhlq/uTIhFA7S3Fbnyn8P1uRuPGoJfjyaJB5PkognScSSV+bBUskUyWR6niyppFLWfNnqPJpemWdLz7ml1BLFjPk3KXCMfz5CPwTsyNjuBy7l2SeSx7FbRhq3Ibf/WqFPWxbc8o5Br00wGKzRf7jeejV02fd5ZJMXBAJCTV2YmrrKmpfKZxj3PLBHRAZFJALcBzzh6PME8GGxeAswq6ojeR5rMBgMhiJy1RG9qiZE5OPAd7FCJP9MVY+LyAPp/Q8BT2JF3JzCCq/86HrHFuWTGAwGgyErZsGUwWAw+ID1wivNDJzBYDD4HCP0BoPB4HOM0BsMBoPPMUJvMBgMPqcsJ2NFZBw477Ud69AJTHhtxFUodxvL3T4ofxvL3T4ofxvL3T7I38ZdqtqVbUdZCn25IyKHc81ulwvlbmO52wflb2O52wflb2O52weFsdG4bgwGg8HnGKE3GAwGn2OEfnM87LUBeVDuNpa7fVD+Npa7fVD+Npa7fVAAG42P3mAwGHyOGdEbDAaDzzFCbzAYDD7HCH0WROTPRGRMRI5l2fcbIqIi0pnR9hkROSUir4vIe72wT0Q+JyLDIvJy+vU+r+zLZWO6/RNpO46LyH/1ysYc3+FfZ3x/50TkZa/sW8fGm0Tk2bSNh0XkzV7ZmMO+gyLyYxE5KiLfEZHmjH1efIc7ROQpETmR/pv7ZLq9XUT+SUROpv9t88LOdey7N72dEpFDjmM2bt9qpRPzWnsB7wBuAY452ndgpVw+D3Sm2w4ArwA1wCBwGgiW2j7gc8BvZOlbcvvWsfFO4J+BmvR2dzl9h479/yfwn8vwO/wecHf6/fuAp8vpO8SqQfHO9PtfAf7A4++wB7gl/b4JeCNty38FPp1u/zTwf3hh5zr27Qf2Ak8DhzL6b8o+M6LPgqo+A0xl2fV/Ab+FvRziPcDXVTWqqmexcvK/OcuxpbAvGyW3D3La+J+AP1bVaLrPmFc2rvcdilUg9xeBv/LKvnVsVGB1lNzCWsW2cvkO9wLPpN//E/AfvLIvbeOIqr6Yfj8PnMCqZX0P8Bfpbn8B/Hsv7Mxln6qeUNXXsxyyKfuM0OeJiHwAGFbVVxy7chVG94KPi8iR9CP16qNoOdl3HfB2EXlORP5VRH4q3V5ONgK8HRhV1ZPp7XKy71PAfxORi8Dngc+k28vFxmPAB9Lv72WtlKjn9onIAHAz8BywTa0qeKT/7U5388xOh3252JR9RujzQETqgd8B/nO23VnavIhZ/TJwDXATMILleoDysQ+simZtwFuA3wS+kR49l5ONAB9kbTQP5WXffwJ+TVV3AL8GfCXdXi42/grwoIi8gOWKiKXbPbVPRBqBbwGfUtW59bpmaSu6ncW2zwh9flyD5Q97RUTOYRU5f1FEtpNf8fSio6qjqppU1RTwCGuPc2VhX4Ytf6sWPwFSWAmbysZGEQkBvwD8dUZz2dgH/DLwt+n336TM/p9V9TVV/RlVvRXrZnnaa/tEJIwloo+r6up3NyoiPen9PcCqG7HkduawLxebss8IfR6o6lFV7VbVAVUdwPqyb1HVy1jFzu8TkRoRGQT2AD8ptY2rf7Rpfh7rEZpysS/Nt4F3AYjIdUAEKytfOdn408BrqjqU0VZO9l0C3pl+/y5g1b1UFjaKSHf63wDwu8BDXtqXfmL8CnBCVf97xq4nsG6apP/9ey/sXMe+XGzOvmLPelfiC2skMgLEsUT9Vx37z5GOuklv/w7WyOV10hERpbYP+CpwFDiS/mPo8cq+dWyMAH+JdRN6EXhXOX2H6fZHgQey9C+X7/AO4AWsyIvngFvL6TsEPokVOfIG8MekV997+B3egeXaOAK8nH69D+gAvo91o/w+0O6FnevY9/Pp7zQKjALf3Yp9JgWCwWAw+BzjujEYDAafY4TeYDAYfI4ReoPBYPA5RugNBoPB5xihNxgMBp9jhN5QsYjIgDM7Zh7HPCAiH75Kn4+IyP+TY99vr3OciMi/ZGZs3Cwi8s+ZGRUNhq1ghN5QVajqQ6r62BZOkVPoseKfX9H1l7Dny1eBjxXgPAaDEXpDxRMUkUfSubu/JyJ1ACJyjYj8o4i8ICL/JiL70u2fE5HfSL//qXQSuB+LyH9zPB30po8/Kem8+SLyx0CdWLngH89iy4dYW2GJiHw4ff5XROSr6bZHReTL6RzkZ0TknekkdCdE5NGMcz2BlXPHYNgyRugNlc4e4Euqej0ww1pa3IeBT6iVc+U3gD/NcuyfY62CfSuQdOy7Cfgl4Ebgl0Rkh6p+GlhW1ZtU9UNZznc71qpVROR6rBWM71LVg1grRldpw0pf8GvAd7DSX18P3CgiNwGo6jRQIyId+X8VBkN2Ql4bYDBskbOquloJ6gVgIJ0J8G3AN61UIoBVqOEKItIKNKnqj9JNXwPen9Hl+6o6m+77KrALe3rYbLSrlVMcLCH/G1WdAFDVzLzt31FVFZGjWOmQj6avcxwYwFoGD1airV5g8irXNRjWxQi9odKJZrxPAnVYT6ozqnrTOsdlS/e63nnz+a0kRCSgVgZRIXf62NVzpxzXSTmuUwss53Fdg2FdjOvG4DvSk6FnReReuBINc9DRZxqYF5G3pJvuy/P08XRa2Wy8DuxOv/8+8IurrhcRad/IZ0hnNdyOlUDPYNgSRugNfuVDwK+KyCvAcawSw4fPvQAAALJJREFUbE5+FXhYRH6MNQKfzeO8DwNHckzG/k/g3wGo6nHgD4F/TduQTwraTG4FnlXVxAaPMxhcmOyVhqpFRBpVdSH9/tNYqZ0/eZXD1jtfD/CYqr6nALb9CfCEqn5/q+cyGIyP3lDN/KyIfAbrd3Ae+MhWTqaqI+lQz+YCxNIfMyJvKBRmRG8wGAw+x/joDQaDwecYoTcYDAafY4TeYDAYfI4ReoPBYPA5RugNBoPB5/z/6+PdRz7JcWkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs1, ys1 = eval_pdf(male_height)\n",
    "plt.plot(xs1, ys1, label='male', linewidth=4, color=COLOR2)\n",
    "\n",
    "xs2, ys2 = eval_pdf(female_height)\n",
    "plt.plot(xs2, ys2, label='female', linewidth=4, color=COLOR3)\n",
    "\n",
    "plt.xlabel('height (cm)')\n",
    "plt.legend()\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sample the population using the `.rvs()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "\n",
    "male_sample = male_height.rvs(N)\n",
    "\n",
    "female_sample = female_height.rvs(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stats from sampling the population:  \n",
    "The sample mean is close to the population mean, but not exact, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(178, 177.98701656715224, 7.723952507001857)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "(163, 163.3122057227971, 7.364243464560319)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean1, std1 = male_sample.mean(), male_sample.std()\n",
    "mu1, mean1, std1\n",
    "\n",
    "mean2, std2 = female_sample.mean(), female_sample.std()\n",
    "mu2, mean2, std2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Percentage differences\n",
    " \n",
    "## Compared to what? \n",
    "\n",
    "#### The answer to that question is the name of the series which is subtracted from another to obtain a difference AND which is used to normalize this difference, hence this series has to be the same. Even Allen Downey got this wrong...\n",
    "#### Another notable fact: a percentage differece is a _signed_ quantity!  \n",
    "\n",
    "In general a pct difference formula explicitely tells you which quantity is the reference from the difference formula (but only if you can see the formula). \n",
    "\n",
    "In the `pct_means_diff1` formula below, the reference is the mean of female_sample because that is the quantity subtracted from the other (male_sample) mean AND it is used to normalize the diff viz the female_sample data.  Therefore, it is a consistent formula:\n",
    "```\n",
    "# compared to the female variable:\n",
    "pct_means_diff1 = (male_sample.mean() - female_sample.mean()) / female_sample.mean()\n",
    "\n",
    "# compared to the male variable:\n",
    "pct_means_diff2 = (female_sample.mean() - male_sample.mean()) / male_sample.mean()\n",
    "```\n",
    "Here are a few other examples:\n",
    "```\n",
    "# %change in stock market values (compared to previous day):\n",
    "(closing_prices_today - closing_prices_yesterday) / closing_prices_yesterday\n",
    "\n",
    "# guess what this one is:\n",
    "(Jetter - DiMaggio) / DiMaggio\n",
    "```\n",
    "## The *correct*  pct difference formula, will yield the *correct* narrative, e.g.:\n",
    "\n",
    "```\n",
    "On average, the male height (cm) is 9.47% higher than that of the female.\n",
    "\n",
    "# see code cell below\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_pct_change_story(series, series_names, measure, units):\n",
    "    \"\"\"\n",
    "    Calculate the % difference between the means of series[0] and series[1]\n",
    "    using the latter as reference and printout a meaningful narrative.\n",
    "    Input:\n",
    "      series (list[np arrays]): series[1] is the reference\n",
    "      series_names (list): [series1_name, series2_name]\n",
    "      measure (string): what is measured, e.g. height\n",
    "      units (string): series measurement units (common to both)\n",
    "    Output:\n",
    "      (string): Formatted narrative\n",
    "    TODO:\n",
    "      Complete input validation\n",
    "    \"\"\"\n",
    "    if len(series) != 2:\n",
    "        return 'Expecting a list of 2 np arrays.'\n",
    "    if len(series_names) != 2:\n",
    "        return 'Provide a name for each series in this order: [series1_name, series2_name].'\n",
    "    \n",
    "    mu1, mu2 = series[0].mean(), series[1].mean()\n",
    "    pct_diff = (mu1 - mu2) / mu2\n",
    "    comparator = 'higher' if pct_diff > 0 else 'lower'\n",
    "\n",
    "    s = 'On average, the {} {} ({}) is {:.2%} {} than that of the {}. (Value: {:.2f}.)'\n",
    "    print(s.format(series_names[0], measure, units, abs(pct_diff), comparator, series_names[1], pct_diff))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "On average, the female height (cm) is 8.24% lower than that of the male. (Value: -0.08.)\n",
      "On average, the male height (cm) is 8.99% higher than that of the female. (Value: 0.09.)\n"
     ]
    }
   ],
   "source": [
    "units = 'cm'\n",
    "what = 'height'\n",
    "\n",
    "which1 = ['female', 'male']\n",
    "case1_series = [female_sample, male_sample]\n",
    "\n",
    "mean_pct_change_story(case1_series, which1, what, units)\n",
    "\n",
    "which2 = ['male', 'female']\n",
    "case2_series = [male_sample, female_sample]\n",
    "\n",
    "mean_pct_change_story(case2_series, which2, what, units)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part II\n",
    "\n",
    "An alternative way to express the difference between distributions is to see how much they **overlap**.\n",
    "\n",
    "### The \"overlap\" is the area under the curves that ends up on the wrong side of a threshold.  \n",
    "\n",
    "There are different way to compute a threshold:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Simple threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "170.59325465118894"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Simple:\n",
    "simple_thresh = (mean1 + mean2) / 2\n",
    "simple_thresh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* pdf crossing point threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "170.30008218110117"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pdf_cross:\n",
    "pdf_thresh = (std1 * mean2 + std2 * mean1) / (std1 + std2)\n",
    "pdf_thresh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simple threshold:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "182"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "159"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pdf cross threshold:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "174"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "165"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# How many in each group are below/above the threshold:\n",
    "print('Simple threshold:')\n",
    "male_below_simple_thresh = sum(male_sample < simple_thresh)\n",
    "male_below_simple_thresh\n",
    "\n",
    "female_above_simple_thresh = sum(female_sample > simple_thresh)\n",
    "female_above_simple_thresh\n",
    "\n",
    "print('Pdf cross threshold:')\n",
    "male_below_pdf_thresh = sum(male_sample < pdf_thresh)\n",
    "male_below_pdf_thresh\n",
    "\n",
    "female_above_pdf_thresh = sum(female_sample > pdf_thresh)\n",
    "female_above_pdf_thresh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overlap:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.174"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.165"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "male_overlap = male_below_pdf_thresh / N\n",
    "male_overlap\n",
    "\n",
    "female_overlap = female_above_pdf_thresh /N\n",
    "female_overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Misclassification rate:  \n",
    "The fraction of people who would be misclassified if you tried to use height to guess sex, is the average of the male and female overlap rates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.16949999999999998"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "misclassification_rate = (male_overlap + female_overlap) / 2\n",
    "misclassification_rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability of superiority (for a statistics):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In 92.20% of randomly chosen pairs (1 man, 1 woman), the man will be taller.\n"
     ]
    }
   ],
   "source": [
    "prob_male_taller = (male_sample > female_sample).mean()\n",
    "\n",
    "print(f'In {prob_male_taller:.2%} of randomly chosen pairs (1 man, 1 woman), the man will be taller.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overlap, misclassification rate and \"probability of superiority\" have two good properties:\n",
    "\n",
    "* As probabilities, they don't depend on units of measure, so they are comparable between studies.\n",
    "\n",
    "* They are expressed in operational terms, so a reader has a sense of what practical effect the difference makes.\n",
    "\n",
    "### There is one other common way to express the difference between distributions: Cohen's effect size.\n",
    "\n",
    "* Cohen's $d$ is the difference in means, standardized by dividing by the standard deviation. Mathematically:  \n",
    "\n",
    "\n",
    "$ d = \\frac{\\bar{x}_1 - \\bar{x}_2} s $ where $s$ is the pooled standard deviation:  $s = \\sqrt{\\frac{n_1 s^2_1 + n_2 s^2_2}{n_1+n_2}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation of Cohen's d value:\n",
    "\n",
    "I dug a bit deeper into the effect size calculations and visited these resources:  \n",
    "* https://www.statisticshowto.datasciencecentral.com/pooled-standard-deviation/\n",
    "> You can only use the above formulas if the standard deviations for the two groups are the same \n",
    "(this is because it would otherwise be violating the assumption of homogeneity of variances. \n",
    " If the standard deviations are different, run Hedge’s g or Glass’s Delta instead.\n",
    "* https://stackoverflow.com/questions/21532471/how-to-calculate-cohens-d-in-python\n",
    "* https://en.wikipedia.org/wiki/Effect_size#Cohen's_d \n",
    "* https://en.wikipedia.org/wiki/Pooled_variance:\n",
    " > Definition: The pooled variance is an estimate of the fixed common variance {\\displaystyle \\sigma ^{2}}\\sigma ^{2} underlying various populations that have different means.\n",
    "* https://github.com/statsmodels/statsmodels/issues/5896 \n",
    "\n",
    "Thanks to the last one (Statsmodels repo), I found out that a statistical analysis package already has an effect size function implemented for multiple cases.  \n",
    "This package is mysteriously call [Pengouin](https://pingouin-stats.org/index.html), but does wonder in whetting one's appetite:  \n",
    ">Pingouin is designed for users who want simple yet exhaustive stats functions.\n",
    "For example, the ttest_ind function of SciPy returns only the T-value and the p-value. By contrast, the ttest function of Pingouin returns the T-value,  \n",
    "the p-value, the degrees of freedom, the effect size (Cohen’s d), the 95% confidence intervals of the difference in means, the statistical power and the Bayes Factor (BF10) of the test.\n",
    "```\n",
    "> conda install -c conda-forge pingouin\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CohenEffectSize(group1, group2):\n",
    "    \"\"\"Compute Cohen's d.\n",
    "\n",
    "    group1: pd Series or NumPy array\n",
    "    group2: pd Series or NumPy array\n",
    "\n",
    "    returns: float\n",
    "    \"\"\"\n",
    "    n1, n2 = len(group1), len(group2)\n",
    "    \n",
    "    diff = group1.mean() - group2.mean()\n",
    "    var1, var2 = group1.var(), group2.var()\n",
    "\n",
    "    pooled_var = (n1 * var1 + n2 * var2) / (n1 + n2)\n",
    "    d = diff / np.sqrt(pooled_var)\n",
    "    return d\n",
    "\n",
    "\n",
    "def overlap_superiority(control, treatment, n=1000, use_pdf_threshold=True):\n",
    "    \"\"\"Estimates overlap and superiority based on a sample.\n",
    "    \n",
    "    control: scipy.stats rv object\n",
    "    treatment: scipy.stats rv object\n",
    "    n: sample size\n",
    "    use_pdf_threshold: bool, \n",
    "       if True: pdf_cross threshold is used else, simple threshold.\n",
    "       Implemented to plot the effect of the thresholding method.\n",
    "    \"\"\"\n",
    "    control_sample = control.rvs(n)\n",
    "    treatment_sample = treatment.rvs(n)\n",
    "    \n",
    "    mu1, mu2 = control.mean(), treatment.mean()\n",
    "    \n",
    "    if use_pdf_threshold:\n",
    "        std1, std2 = control.std(), treatment.std()\n",
    "        # pdf_cross threshold:\n",
    "        thresh = (std1 * mu2 + std2 * mu1) / (std1 + std2)\n",
    "    else:\n",
    "        # simple threshold:\n",
    "        thresh = (mu1 + mu2) / 2\n",
    "    \n",
    "    control_above = sum(control_sample > thresh)\n",
    "    treatment_below = sum(treatment_sample < thresh)\n",
    "    \n",
    "    overlap = (control_above + treatment_below) / n\n",
    "    \n",
    "    superiority = (treatment_sample > control_sample).mean()\n",
    "    return overlap, superiority\n",
    "\n",
    "\n",
    "def plot_pdfs(cohen_d=2, pdf_threshold=True):\n",
    "    \"\"\"Plot PDFs for distributions that differ by some number of stds.\n",
    "    \n",
    "    cohen_d: number of standard deviations between the means\n",
    "    pdf_threshold: bool, use pdf_cross threshold (to pass to overlap_superiority().)\n",
    "    \"\"\"\n",
    "    control = sps.norm(0, 1)\n",
    "    treatment = sps.norm(cohen_d, 1)\n",
    "    \n",
    "    xs, ys = eval_pdf(control)\n",
    "    plt.fill_between(xs, ys, label='control', color=COLOR3, alpha=0.7)\n",
    "\n",
    "    xs, ys = eval_pdf(treatment)\n",
    "    plt.fill_between(xs, ys, label='treatment', color=COLOR2, alpha=0.7)\n",
    "    \n",
    "    o, s = overlap_superiority(control, treatment, use_pdf_threshold=pdf_threshold)\n",
    "    plt.text(0, 0.05, 'overlap ' + str(o))\n",
    "    plt.text(0, 0.15, 'superiority ' + str(s))\n",
    "    \n",
    "    plt.text(-4, 0.35, f\"Cohen's d: {cohen_d:.2f}\")\n",
    "    if pdf_threshold:\n",
    "        tresh = 'Pdf cross'\n",
    "    else:\n",
    "        tresh = 'Simple'\n",
    "    plt.text(-4, 0.39, f\"Threshold: {tresh}\")\n",
    "    \n",
    "    plt.legend()\n",
    "    plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZgU1bn/P29V9yzAgMgiIhpwRdZBFtewGNeoUeK+E1yTYGK8MXI1vxtj1JjovVe9wRhNjFFj1GgwRNEo7gsqoLgAgggDDAyzb909Pb29vz+qZxiGnpke6Jnurj6f55lnuqrOqTpVXfWt0+95z/uKqmIwGAwG92KluwEGg8Fg6FmM0BsMBoPLMUJvMBgMLscIvcFgMLgcI/QGg8HgcjzpbkAiBg8erCNHjkx3MwwGgyFrWLFiRZWqDkm0LSOFfuTIkSxfvjzdzTAYDIasQUQ2dbTNmG4MBoPB5RihNxgMBpdjhN5gMBhcTkba6A0GgzsIh8OUlpYSDAbT3RTXUFBQwIgRI/B6vUnXSUroReQU4D7ABv6oqnd1UG4q8AFwvqo+2526BoPBfZSWllJUVMTIkSMRkXQ3J+tRVaqrqyktLWXUqFFJ1+vSdCMiNrAAOBUYA1woImM6KPcb4N/drZtKqqurKS4upri4mGHDhrHffvtRXFzMXnvtxZgxqT/0m2++yemnn96tOjNnzkzoVfToo48yb968TuuWlJRQWFhIcXExY8aM4dprryUWi+1Sbs6cOTz77LMAvPPOO4wdO5bi4mKampq61VaDYU8IBoMMGjTIiHyKEBEGDRrU7V9IydjopwHrVXWDqoaAp4AzE5S7DngOqNiNuilj0KBBrFy5kpUrV3Lttdfyk5/8pHXZsro+3Ugk0pPNSwkHHXQQK1eu5LPPPmP16tU8//zznZb/61//yk9/+lNWrlxJYWFht46VDdfDkNkYkU8tu3M9kxH6/YAtbZZL4+vaHng/YDbwYHfr9ibRaJSrrrqKsWPHctJJJ7X2bmfOnMnNN9/MjBkzuO+++6isrOTss89m6tSpTJ06lffeew+At956q/XXwqRJk2hsbATA5/NxzjnnMHr0aC6++GJaQj+/9tprTJo0ifHjxzN37lyam5t3adOf//xnDj30UGbMmNF6nGTxeDwcc8wxrF+/HlVl3rx5jBkzhtNOO42KCud9+8c//pFnnnmG2267jYsvvniXfTz22GNMmDCBiRMncumllwLOr4EbbriBWbNmcdNNN1FTU8NZZ53FhAkTOOqoo/jss886vB5lZWVMnz6d4uJixo0bxzvvvNOtczIYDKknGRt9otdH+yD29wI3qWq03dsmmbpOQZGrgasBDjjggCSa1X2++uor/va3v/Hwww9z3nnn8dxzz3HJJZcAUFdXx1tvvQXARRddxE9+8hOOO+44Nm/ezMknn8yaNWu45557WLBgAcceeyw+n4+CggIAPvnkE1atWsXw4cM59thjee+995gyZQpz5szhtdde49BDD+Wyyy7j97//Pddff31re8rKyvjFL37BihUrGDBgALNmzWLSpEkALFq0iOXLl3Pbbbd1eD6BQIDXXnuN2267jYULF7J27Vo+//xzysvLGTNmDHPnzuXKK6/k3Xff5fTTT+ecc87Zqf6qVau44447eO+99xg8eDA1NTWt29atW8eSJUuwbZvrrruOSZMm8fzzz/P6669z2WWXsXLlyoTX46GHHuLkk0/mlltuIRqNEggEUvPlxdFADbr9UwjF91u4FzJsApJflNLjZALBQIiyTXUEA2EA8gs97PuNgRT2zUtzy3af2BfPQbA+dTssGIA17uzU7a8DSkpKeP/997nooou6Xe/000/niy++6KGWJUcyQl8K7N9meQSwrV2ZKcBTcZEfDHxbRCJJ1gVAVR8CHgKYMmVKj2RDGTVqFMXFxQBMnjyZkpKS1m3nn39+6+clS5awevXq1uWGhgYaGxs59thjueGGG7j44ov57ne/y4gRIwCYNm1a6+fi4mJKSkooKipi1KhRHHrooQBcfvnlLFiwYCeh//DDD5k5cyZDhgxpbcO6desA+M53vsN3vvOdhOfx9ddfU1xcjIhw5plncuqpp3L99ddz4YUXYts2w4cP5/jjj+/yerz++uucc845DB48GIC99967ddu5556LbdsAvPvuuzz33HMAHH/88VRXV1NfX5/wekydOpW5c+cSDoc566yzWq/3nqLhALr5Q6grARGQ+K3rr0Srv0IHH4bsNwWxk/dEyFQikSilX1VTtd35xWjZTn/J39BMTbmPvffpxwGHDMbjtdPZzN0jWA+pfCmn8qXRCSUlJTz55JMJhT4SieDxZLYDYzKtWwYcIiKjgK3ABcBOZ6uqrcO/IvIo8IKqPi8inq7q9ib5+fmtn23b3mlgsm/fvq2fY7EYS5cu3cWePX/+fE477TQWL17MUUcdxZIlSxLuNxKJkGzmrt2xt7XY6Pd0X6raYZ221yPRuYhIwusxffp03n77bV588UUuvfRSbrzxRi677LJutWuXdoab0HUvQbABvH0doW8lHzQGlWvQYD0cfAJiZfZD1xnRSIyvPi3DV99MXoG98/fjjXtdbG+kyR/isOLh2Sn2aeCxxx7jnnvuQUSYMGECt99+O3PnzqWyspIhQ4bw5z//mQMOOIA5c+bQv39/li9fzvbt2/ntb3/LOeecw/z581mzZg3FxcVcfvnlDBw4kBdffJFgMIjf7+e1117jZz/7GS+99BIiws9//vOdOo/ppksbvapGgHk43jRrgGdUdZWIXCsi1+5O3T1vds9y0kkn8bvf/a51uUVUv/76a8aPH89NN93ElClT+PLLLzvcx+jRoykpKWH9+vUAPP7448yYMWOnMkceeSRvvvkm1dXVhMNh/v73v+92m6dPn85TTz1FNBqlrKyMN954o8s63/rWt3jmmWeorq4G2Ml0037ff/3rXwHHy2jw4MH0798/4fXYtGkTQ4cO5aqrruKKK67g448/3u1zAtBIEP3qZWhuhLx+7UQ+jljg7QeNZejXr6Ox6B4dM11EozG++qwMf0MCkY8jIuQVeGjyhVj3aRmRSHaea2/SYqJ8/fXX+fTTT7nvvvuYN28el112GZ999hkXX3wxP/rRj1rLl5WV8e677/LCCy8wf/58AO666y6++c1vsnLlSn7yk58AsHTpUv7yl7/w+uuv849//IOVK1fy6aefsmTJEm688UbKysrScr6JSGpmrKouVtVDVfUgVb0jvu5BVW0/+Iqqzmnxoe+obqZz//33s3z5ciZMmMCYMWN48EHnNO+9917GjRvHxIkTKSws5NRTT+1wHwUFBfz5z3/m3HPPZfz48ViWxbXX7vxe3Hfffbn11ls5+uijOeGEEzjiiCNaty1atIj/+q//SrrNs2fP5pBDDmH8+PF8//vf3+WlkoixY8dyyy23MGPGDCZOnMgNN9yQsNytt97aej3mz5/PX/7yFyDx9XjzzTdbB2efe+45fvzjHyd9Du1RjaHrl0BTndOT7wwRp0xDKbrp/d0+ZrpQVTaurqCxPog3P7HItyAiePNtAo3NfP1FedK/HnOVRCbKpUuXtpphLr30Ut59993W8meddRaWZTFmzBjKy8s73O+JJ57Yau589913W02n++yzDzNmzGDZsmU9eFbdQzLxJpkyZYqa6JWGWPkqKP0ogbmmE1Qh7IODTsTaa/+uy2cI1eWNbFhV0WFPPhGqSqg5yjcOG8LQ4f17uIW7x5o1azj88MNbl2PLH0mtjb65EWvK3E6L3H///VRUVHD77be3rhs8eDBlZWV4vV7C4TDDhw+nsrKSOXPm7OS40K9fP3w+H2+++Sb33HMPL7zwAuDMeVm+fHnrL//rr7+eCRMmMHeu05ZLL72Uc889lwkTJvTIYGz76wogIitUdUqi8ibWjSEj0WA9bF0OnsLkRR6csnY+bHoXjezqzpqJhEMRNq2twuO1ujXOIiJ4vBZbvqoiFDTzHToikYnymGOO4amnngKceSbHHXdcp/soKipqdadOxPTp03n66aeJRqNUVlby9ttvM23atNSdxB6SvaNWBteiGkNL3gEUrN0YbLTzIORDSz9CRn4z1c1LKarKpnVVxKIx8gq6/zjatkU0HKFkbQWHTNg38ycnFQxIuXtlV7Q1Udq2zaRJk7j//vuZO3cud999d+tgbGdMmDABj8fDxIkTmTNnDgMHDtxp++zZs1m6dCkTJ05ERPjtb3/LsGHDdvLsSyfGdGPIOGLVG6Dkre6ZbNqjCmE/Mvp0pG/CpDsZQUNtE+tWbuvSLt8ZLSacg8cPY+DgLsYyeplEJgbDnmNMN4asRjUG25Y7vfI96Z2KgFjo1hWpa1yKUVVKv65GLNmjnriIYNvC1q9rzMCsISFG6A0ZhdZshJDfEfo9xVPouFz6KroumwYaa5sINDbj8e75Y2h7LIKBEHVV/hS0zOA2jNAbMgaNRZ0B2FSIPOzo1W9bkXE9XVVly9c1e9ybb0FEsGyh1PTqDQkwQm/IGLR2I4QDqRN6iPfqt4O/MnX7TAH1NQECvtT05luwPRbBpjC1laZXb9gZI/SGjEBVYdvHqRV52NGrL9s1ZEQ6KSupxUpRb76FFlv9to21pldv2Akj9IbMoLHMiUiZaqEHp1ffsA1t7tgPujdp8oXwN6S2N99Ci60+0JgdcwgMvYPxozdkBFqxes+8bDpDBFC0ai2yX0Lvs16lYqvjR94TPu/OPoXy0noOHFOQ8v3vKZ9/sLk17HIqKOjjZfxRHYc1r6ur48knn+QHP/hBSo537733cvXVV9OnT5+U7K+3jmN69Ia0oyEfNJQ6Pe+ewlMAlV+isfTOII1EolRtb8ST13NRJ715FjUVfsKhzJstGwyEyS/0pOyvq5dGXV0dDzzwwC7ro9HdCwZ37733pjzHQm8cx5VCv337di644AIOOuggxowZw7e//e3WOO+J2J28r53ta86cObtV99Zbb+Wee+7ptEx1dTWzZs2iX79+neaX/fTTTzn66KMZP348Z5xxBg0NDQB89NFHrVmhJk6cyMKFC3erralEK9c5E5x6clan5YFoBK3d3HPHSILaCj8aUyyr585VLIF4OONcZ/78+a35G6ZOncqsWbO46KKLGD9+PABPPPEE06ZNo7i4mGuuuab1BfD973+fKVOmMHbsWH7xi18ATsycbdu2MWvWLGbNmgU4sXBuuukmJk+ezAknnMBHH33EzJkzOfDAA1m0aBHgvFRuvPFGpk6dyoQJE/jDH/4AOFoxc+bMXbLTJTrOnuI6oVdVZs+ezcyZM/n6669ZvXo1d955Z6dR6LKJgoICfvWrX3X5Qrjyyiu56667+Pzzz5k9ezZ33303AOPGjWP58uWsXLmSl19+mWuuuSateWE1FoXKNU6Pu6exbKj4PG0DlapK2aY6LE/PP3a2x2L7lvqcH5S96667WvM33H333Xz00UfccccdrF69mjVr1vD000/z3nvvsXLlSmzbbg3Hfccdd7B8+XI+++wz3nrrLT777DN+9KMfMXz4cN54443WMOB+v5+ZM2eyYsUKioqK+PnPf86rr77KwoULW6PP/ulPf2LAgAEsW7aMZcuW8fDDD7Nx40bAyU537733snr1ajZs2MB7772X8Dh7iuuE/o033sDr9e4UEri4uJhvfvObqCo33ngj48aNY/z48Tz99NOtZTrK+7pixQpmzJjB5MmTOfnkk1tjTM+cOZObbrqJadOmceihh7bmRs3Ly2PAACf+Rkc5Zttyxx13cNhhh3HCCSewdu3aLs+vb9++HHfcca1pDDti7dq1TJ8+HXDCqbZkiOrTp09rNpxgMJj+2CgNWyEWdnrcPY2dD4EaCNb2/LES4G9sJhQMY9s9f81tj0UkFKWxLtjjx8ompk2bxqhRTp6k1157jRUrVjB16lSKi4t57bXX2LBhAwDPPPMMRxxxBJMmTWLVqlU7ZZxrS15eHqeccgoA48ePZ8aMGXi9XsaPH98a5+aVV17hscceo7i4mCOPPJLq6mq++uqr1vaMGDECy7Jas9P1BK4bjP3iiy+YPHlywm1tkwNUVVUxderUVjFMlPf1yCOP5LrrruOf//wnQ4YM4emnn+aWW27hkUceAZwUYh999BGLFy/ml7/8JUuWLOGYY47hmGOOAegwx2wLK1as4KmnnuKTTz4hEolwxBFHtLa9JQZ++xj2yTJu3DgWLVrEmWeeyd///ne2bNmRo/3DDz9k7ty5bNq0iccffzytadC0ai2JUwv3APGBSq3ZgOy3d5fFU011WSNKzwzCJkSgansj/Qf24NhHltE+c9rll1/Or3/9653KbNy4kXvuuYdly5YxcOBA5syZQzCY+IXp9Xpbv0/LslqzzVmW1fpLWVX5v//7P04++eSd6r755psJs9P1BEn16EXkFBFZKyLrRWR+gu1nishnIrJSRJaLyHFttpWIyOct21LZ+O7SWXKARG/WtWvX8sUXX3DiiSdSXFzM7bffTmlpaev+vvvd7wK75p9toSWn6v33309dXd0ugvrOO+8we/Zs+vTpQ//+/XfKEXvttdfutsgDPPLIIyxYsIDJkyfT2NhIXt4Ot8UjjzySVatWsWzZMn796193eBP3NBoJOj363jDbtODJh6qvnJg6vUgsFqO63Ie3Bwdh2+P12tSW+4hGevdcM4nOwgt/61vf4tlnn6WiwgmRUVNTw6ZNm2hoaKBv374MGDCA8vJyXnrppaT21xEnn3wyv//97wmHnYHjdevW4fd3Pqltd47TGV125UTEBhYAJ+Ik+14mIotUte1vmdeARaqqIjIBeAYY3Wb7LFWtSlmrO2Hs2LE8++yzCbd1Zq/sKO/r2LFjWbp0aad1OnoTJ8qpOnr06J3K9FTvbvTo0bzyyiuAc2O9+OKLu5Q5/PDD6du3L1988QVTpvS+26HWbY4PwvaiBdHyQNgPvgooGtZrh22sDRKLao/4zneEWIKqUl/jZ++hKUz2sQcU9PGm3L2yMwYNGsSxxx7LuHHjKCwsZJ999mndNmbMGG6//XZOOukkYrEYXq+XBQsWcNRRRzFp0iTGjh3LgQceyLHHHtta5+qrr+bUU09l3333Tdp+fuWVV1JSUsIRRxyBqjJkyBCef/75TuvsznE6I5nf7NOA9aq6AUBEngLOBFqFXlV9bcr3BdI2AnT88cdz88038/DDD3PVVVcBsGzZMgKBANOnT+cPf/gDl19+OTU1Nbz99tvcfffdHeZ+Peyww6isrGTp0qUcffTRhMNh1q1bx9ixY5NqS0tO1fHjx7N06VK+/PLLnYR++vTpzJkzh/nz5xOJRPjXv/7FNddcs+cXAaioqGDo0KHEYjFuv/321l8HGzduZP/998fj8bBp0ybWrl3LyJEjU3LMblO5tnds8wnQmvVILwp9ZVkDIr3/WIglVG5rzBih78znvad48sknO9x2/vnnJ0zi/eijjyYsf91113Hddde1Lvt8O6Tv1ltv3alsyzbLsrjzzju58847d9o+c+ZMZs6c2brcNk91++PsKcl0L/YDtrRZLo2v2wkRmS0iXwIvAm1zeynwioisEJGrOzqIiFwdN/ssr6zc/bgkIsLChQt59dVXOeiggxg7diy33norw4cPZ/bs2UyYMIGJEydy/PHHtyYH6Ii8vDyeffZZbrrpJiZOnEhxcTHvv598PtKucsweccQRnH/++RQXF3P22WfzzW/uSJLx4IMPttrp2zNy5EhuuOEGHn30UUaMGNE6UHTllVfSEsf/b3/7G4ceeiijR49m+PDhfO973wMc81XLucyePZsHHnigNZdmb6LNjRCocgZIexu7AGo29JpPfSQcpb4q0KO+8x3h8Vo01gYz0qfe0Ht0mXhERM4FTlbVK+PLlwLTVDXh60ZEpgP/paonxJeHq+o2ERkKvApcp6pvd3ZMk3jE/cS2fwZbV0Bev/Q0IORHDpqF7PWNHj9U1fZGStZU7FYGqVQQCkY44NDBDN2v62xMqcYkHukZeiLxSCnQNsvyCGBbR4XjIn6QiAyOL2+L/68AFuKYggw5jKpC1bqeiWuTLCJo9de9cqjqskZnElOasGyhalv6Jk/lui9/qtmd65mM0C8DDhGRUSKSB1wALGpbQEQOlviooogcAeQB1SLSV0SK4uv7AicBqU2Hbsg+mhudP6vzgbQexVMA9Vt63HwTCUdprG/q1UHY9tgei4CvmVBz75tvCgoKqK6uNmKfIlSV6urqLufRtKfL35KqGhGRecC/ARt4RFVXici18e0PAmcDl4lIGGgCzo974OwDLIy/AzzAk6r6crdaaHAdWh8PQ5DOyVpiAepEzRywf5fFd5f6mgCQ2nDE3aXl2PXVAYYM79+rxx4xYgSlpaXsybibYWcKCgoYMWJEt+okZTRU1cXA4nbrHmzz+TfAbxLU2wBM7FaLDO6nen16e/MtKM7kqR4U+urtvrS+z1qwbKF6e2OvC73X622diWpIH64LgWDIbDTkg6ba9NrnW/AUQN1mJ95ODxCNxGioTa/ZpgXbY+GrDxIO9cy5GjKb9N+BhpxC6+MzizOjmwsaBV/PBLxrqGsCejHkQSe0mm9qej7EriHzMEJv6F2q16dtklRCNJ6rtgeo2e4jjXMHd0Esoabc13VBg+swQm/oNTTc5CTpTsckqY7wFEDtxpTHvolFY9RV+/F6e3+SVEd4vBYNtU05HfsmVzFCb+g1tGGbE6gyA0wZrVg2RCMQqE7pbn0NzWhM0+o/354W801j3KRkyB2M0Bt6j7oSMvOW0x1jBymirsqfQUabtii1lZ1HTjS4j0x86gwuRGPR3g9JnCx2HtRuSNnuVJXaCn9GeNu0x+O1nZeQmcCUU2TenWhwJ/4K0FjvhiROFssLwQbH9TMFBANhwqFIj+aF3V0sS4hGYgQaQ+luiqEXycCnzuBGtG6LE3s+E4nbrrWhwxBO3aLFhTET3Co7or7GmG9yCSP0hh5HVR37fCZ527RHbEiRm2VthQ+rF/LC7i6Wx6Kmwgh9LmGE3tDzNDdAOJBZ/vPt8eRDY9keBzmLhKP4G5qxPZn7aNm2EPSH0hLkzJAeMvduNLgGbdgaTxmYub3c1rGDPZwl25hBs2E7oqVtDWaWbM5ghN7Q89SWZHZvvgVVtH5L1+U6obbS78wVyHDEEuNmmUMYoTf0KBoNOx43mRDErCvsfKjdtNuuh6rqpAzMoNmwHeHxOLNkY7EMHSA3pBQj9IaexV/h/M9Et8r2WB4IN0Fo97IxNflDRKOakW6V7RFLUFUCjc3pboqhF8iCp8+QzWj91sx1q2yPCKBo4/bdqt5Q25RVE5E0ZqJZ5gpJCb2InCIia0VkvYjMT7D9TBH5TERWishyETku2boGl5PpbpXtEdsZU9gNaiv82J7M7823YHst6oydPifoUuhFxAYWAKcCY4ALRWRMu2KvARNVtRiYC/yxG3UNLkWbfRDyZ8dAbAuefPBt73Yykmgkhr8hmNFule2xbaHJHzLJSHKAZO7KacB6Vd2gqiHgKeDMtgVU1ac7frP2ZUcQ7i7rGtyLNpY5HzLY1XAXxHJsGoGqblXz1QdBMtutsj0tbfXVB9PcEkNPk4zQ7we09Tkrja/bCRGZLSJfAi/i9OqTrhuvf3Xc7LPcJBJ2CfWbHFNItqHa7XAI9TWBTMoxkjziRNo0uJtkhD5RF2WXW1pVF6rqaOAs4FfdqRuv/5CqTlHVKUOGDEmiWYZMRjUGDdscU0i2Yed1KxyCqhP6187AaJVd4fFYJpplDpDMnVkK7N9meQTQYXdHVd8GDhKRwd2ta3ARgerMjVbZFZYXgvVORqwkCAUjhJszM1plV1i2RTQaIxgIp7sphh4kmadwGXCIiIwSkTzgAmBR2wIicrDEDX4icgSQB1QnU9fgTrSxLHvcKtsj4vz5KpIqng1hDzpFTdYpt9OlO4SqRkRkHvBvwAYeUdVVInJtfPuDwNnAZSISBpqA8+ODswnr9tC5GDKJuk1OzzhbUdCGUmTgN7osWl/dlFXjze0RS6ir9DN0vwHpboqhh0jK701VFwOL2617sM3n3wC/Sbauwd1oNAT+avD2SXdTdh9PPtRvQVU77amrKvU1AewsCHvQER6PRWNdkFgshmVloanN0CXmWzWkHn9l5iUB7y5iQzjYZTiEJn+IWCw7wh50REsCc5N1yr0YoTekHG3Ylr32+RZEQBRt7DxscWNd0BUeK7GYGju9izFCb0g9dZuzK+xBh9hQv7nTErWVfuwMziaVLLbHhC12M0boDSlFwwEno1Q2hT3oCDsPGrY5cwISEIvG8NVnV9iDjrA9FoHGZiIREw7BjWT/HWrILHwVO9wTsx3LBo1CU23Czf54iN+sdatsg4iAgL/ehC12I0boDSlF60uzMxRAR2jHYYuzLSxxV6iasMVuxQi9IWWoKjRsyc6wBx1heTq009dVucM+34LHY1FfbYTejRihN6SOkA/CzdkZyKwj7DzwVewStjgSidLkC7nCPt+CZQvNTWHCoUi6m2JIMe65Sw1pR33lgLrDPt9CS6yeQPVOq/0NzVkXlrgrdoQtNnZ6t2GE3pA66rdkZxCzrtDYLnb6hppA1k8VSIhgzDcuxIVPpSEdOPb5rS7xn2+H5XVi67ehriqAx0VmmxbsuJ3eTYPMBiP0hlQRrINYxHFJdBt2HgSq0agTyjccitLcFMZy0UBsC5YlhEMRQs3GTu8mjNAbUoL6yrM/7EFHtNjh/U7mM39DML7afULfaqevM+kF3YQRekNqqNvsLm+b9qi25sB1vQ1bjD+92zBCb9hjVGPg2+4u//n22HnOywyoqw5kZdrAZPEYO73rcO/daug9AjXZmzYwWSwvBOto9vuzNm1gsli2RTSiNDeZ9IJuIaknU0ROEZG1IrJeROYn2H6xiHwW/3tfRCa22VYiIp+LyEoRWZ7KxhsyA/W72D7fggiIha+iMr7oXqF3UBqNnd41dBliUERsYAFwIk6y72UiskhVV7cpthGYoaq1InIq8BBwZJvts1S1KoXtNmQSdZvdEa2yKzRGfUUNiPtT7okl1FcHGDK8f7qbYkgByfTop8a6p14AACAASURBVAHrVXWDqoaAp4Az2xZQ1fdVtSXE3wfAiNQ205CpaCzqRKy089LdlB5HrXwa6sKu9J9vj+2xXBe0LZdJ5o7dD9jSZrk0vq4jrgBearOswCsiskJEru6okohcLSLLRWR5ZWVlEs0yZASBapywB+4Xv5AWEImAuCo8Z2IsS4jFlKDf2OndQDK/txMZIxPe6SIyC0foj2uz+lhV3SYiQ4FXReRLVX17lx2qPoRj8mHKlCnuf5Jcgvq2u98+H8cXKkABiYXALkx3c3ocVaWxvonCfu7/teZ2kumGlQL7t1keAWxrX0hEJgB/BM5U1dYIUKq6Lf6/AliIYwoyuIX6LY5HSg5QH+yDKGgkNwYprbid3pD9JCP0y4BDRGSUiOQBFwCL2hYQkQOAfwCXquq6Nuv7ikhRy2fgJOCLVDXekF40FnFmi9ruF3pVaGwqxLZjEM4NoW+x08diufGLzc10abpR1YiIzAP+DdjAI6q6SkSujW9/EPgvYBDwQNztLKKqU4B9gIXxdR7gSVV9uUfOxND7+KtwYvW63z7fHPESVQuvrU5Mn1jUnXF92mBZQkSVoD9EnyIXT4bLAZLyiVPVxcDidusebPP5SuDKBPU2ABPbrze4A8c+nzhxttvwNefvGJgSgWgzWH3S2aReQRUa65uM0Gc57u+KGXqOus25ZZ9vWVDQHDHfWJZQV2Xs9NmOEXrDbqHRkBP6IBf85xUag4XYVjydoGVBJOiuJOgd4PFY+OqCxk6f5RihN+we/irH8db1oQAgGPESU2FHeBtxTFbq/pjtYjmzBgI+k14wmzFCb9gttDGH/OeDBWj76SSqaCSUngb1NqomPn2WY4TesHvUb84Jsw202OfbDTqLOOabHMCyhboqf7qbYdgDjNAbuo1GmqGpNicGYlWhsbkQj9Ve6C2INOWEnd72WPgamolFc8PDyo0YoTd0H39lPGyv++3zTeE8VBOdqjhvgVgO2OnjJ+9vNHb6bMUIvaHbaMO23LHPNxeg2sELTRWN5ob4aUxprGtKdzMMu4kRekP3qd+SM/b5uqY+WNKByUIkh8IhGH/6bMYIvaFbaLgJmhtywj4fU6dHb7e3z7cgueNPb3ssAo3NRCPGTp+NGKE3dI8css8HQs60/45PVQCFmPtjtkv8Ozd2+uzECL2hW2hDaU70YKEL+3wLqo4XUi4QUxpqjPkmGzFCb+ge9aU5ZJ/vi9WR2aYFEcfNMgewvRZ1Jj59VmKE3pA0GvJDyJ8TicCjMSEQysfuaCC2BbEg0pwTv3IsWwj6Q0TC0XQ3xdBNjNAbksdXnjPxbfyhfECTONV4gaj7wyG02Ol99bnhaeQmjNAbkkbrS8mVW6YxUXybjsghf3pUaajNDVOVm0jqqRWRU0RkrYisF5H5CbZfLCKfxf/eF5GJydY1ZAeqCg25Y5+vD/bt2mzTglgQzg3xsz2WiXuThXQp9CJiAwuAU4ExwIUiMqZdsY3ADFWdAPwKeKgbdQ3ZQHODY57IAft8JGbRFM7reKJUe0Sca5MDs4UtWwgFI4Sa3R/6wU0k06OfBqxX1Q2qGgKeAs5sW0BV31fV2vjiB8CIZOsasgP1leeEkAH4m/ORpOzzLRg7vSGzSUbo9wO2tFkuja/riCuAl7pbV0SuFpHlIrK8srIyiWYZepX6zSDuTobdQkOwsGv/+fbkkj89Sr1xs8wqkhH6RHd8wq6diMzCEfqbultXVR9S1SmqOmXIkCFJNMvQW6jGoKEsZ+zzdU19d6QNTBaxIJwb4ufxWNRXB5xxG0NWkIzQlwL7t1keAWxrX0hEJgB/BM5U1eru1DVkOIEa0ChY7u/RhyI2oagHS7opYmJBNOykGHQ5lm0RCcdobnJ/6Ae3kIzQLwMOEZFRIpIHXAAsaltARA4A/gFcqqrrulPXkPmoL3fSBjY2F3TTPt8GEWfyVE6gNJr0gllDly4UqhoRkXnAvwEbeERVV4nItfHtDwL/BQwCHognKYjEzTAJ6/bQuRh6irpNOeFtA07awN1GQSNBxFuYugZlKGI56QWHDO+f7qYYkiCpp1dVFwOL2617sM3nK4Erk61ryB40GnYiVnr2QACzBFVoaOrTcVjirhBx/OkLBiYenXIRHo9FQ20TsZhiWS4/WReQG9McDbuPvxLIjbDEwYiXmEr37fMtiAWxqDOe4XLEElShyZcrpqrsxgi9oVO0cVtODDBCN8MedIjkjpulCYeQNRihN3RO3eaccqu02MOXmpAzYYstW6itNOEQsgEj9IYO0XAAgvU5kjZQOk8bmCxiOXlkc8BJyaQXzB6M0Bs6xleRM2kD/c1dpQ1MllxLL4gJh5AFGKE3dIjWbcb17iNxGoKFqeuEq6KRHBE/xWSdygKM0BsSoqpOfBs7P91N6RXqmroRlrgrcigcgu11whabcAiZjRF6Q2KaaiEWyYmwB+GoTTDi3X23yva0hkNwv/hZlhBujhAKmrDFmYwRekNCtLEsJ4QKnLAHyaUN7CY54GYZnwlv3CwzHCP0hsTUleRM2IO6pj6pH4lQ0Bxxs2wJh2DIXIzQG3ahNexBDtjnd4Q9SPFsViueXjAHfhS1DYdgyEyM0Bt2xV/h/M8Bt8qmcF487EHXZX/4/26gwdeY5J4FYlF++B+X09BYT0NjPU8/93i327do8bOccd5MzjhvJosWP5uwzLayUq667iLOufQUrvjh+ZRXlAHw5bpVXHrVbGZffCLnXHoKLy/5V7ePnwxiCSgEGt1vqspWJBNHy6dMmaLLly9PdzNyltiWD6ByDXj7pbspPc72hv5srRtEnqfjwURVRVWxrOT7RaqKRsNYfQcj+X3ZWraF6356Bf/46ytJ76O+oY4L557B3x75F4JwwdzTeeqRF+jff8BO5X56yw+YfuzxfOfb5/Dh8vf554t/585f/C8lmzcgInxj/1FUVJZz4dzTWfjkEvoXDejgiLtPczDC8JED2W/U3inftyE5RGSFqk5JtM306A07oapQWwJ2Qbqb0in+QJDTf/hrJp1zIxNm/wdPv/w+AAee8kOqahsAWL7qa46feysAv3zgGS67+f844YpfctjpP+LhZ5cAUBvoxxP/+AsXXTeXc6+9hAcefxiArdvLmH3VBdzxu7u5YN7lbK8s59TLZlNbXwfA48/9jbOvuZizr7mYJxY+lbhOVRWnXnACtXU13PfAbyjduonzLj+V//ndndz8y5/wxts7RP8/b/0xb77z6k7n+P4Hb3HU1OMY0H8v+vcfwFFTj+O9D97c5Vp8XfIVR045FoBpk49u3c/IAw7kG/uPAmDokH3Ye+Agautq9vjaJ8Ljsaip8PXIvg17Tm6MthmSp7nemcLvzeywxC+/t5LhQwbywoL/BKC+sWu/9c/Xbeb9J+7A3xRk8nk3cfKxU3jjo3JKt23hr/f/CVXlx7feyIrPP2HYkGGUlG7mlzf8nFvm3bjTflZ/9SX/fPUFnrjvj6gql/z4SqaMn0RRv/7t6qgzCKDKj39wE+s3rOOZvzjplJd/8gFPPPUnZk0/iUZfAys/X8Gvfv7fOx2noqqcYUOHty7vM3RfKqrKdzmvww4+nCVvvMTF58/ltbf+jT/go66+lr0GDNxx7qtXEg6H2X+/byR9jbuDZQvNgTCh5gh5+UZWMo2kevQicoqIrBWR9SIyP8H20SKyVESaReSn7baViMjnIrJSRIw9JsPRhm2AZrx9fvwhB/DaB58z/3+f4J0VaxhQ1PWL6TuzplBYkMfggf2ZOXUs76zcxAcff8DSjz/k/B9ezgXz5lCyZRObtzr57PcdOowJh4/bZT+frPqU44+ZQWFBIX0K+/CtY2fw8RefJqgTv4bR0C77mDLpKDaXbqK6poqXXl3ECTNPxePZWSATm1V3/V5umHcLy1d+yHmXf5sVn3zA0CHDsO0d8x8qqyq45bYbuO2Wu7tlfuoOrW6WNbkxUSzb6PLVKyI2sAA4EScH7DIRWaSqq9sUqwF+BJzVwW5mqWrVnjbW0AvUlmSFW+WhI4ez7OnfsPidj7nl/ic58eiJ/L9rz8FjW63eH8HmnePNSDuR9IcKUJQrzr+Mc06bvdO2rdvLKCxInCmqs3GtRHWccAi7XtPTT5nN4lee5+Ul/+KXN9+9y/Z9hgxj2ScftC6XV5QxddJRu5QbOmQf/vfXfwAgEPCz5M2XKernZH7y+RuZ99PvMe/q/2DCuCM6bHcqEMuJZjl4X5N1KtNI5vU+DVivqhtUNQQ8BZzZtoCqVqjqMsD9kZxcjEZDTiCzLHCr3FZRQ5+CPC45fTo3XH4GH6/ZAMDI4UNZsdr5/I8lH+xUZ9Gbywg2h6iua+St5asYOXISx06exvOvvECgyemJlldVUNOFHXvyuEm88f7bNAWDNAWbeP39tzli3MSOK4Sb6FvYj0BgZ1/zM087h78+/WcADj7w0F2qHXPUDJZ+9A4NDfU0NNSz9KN3OOaoGbuUq62rIRZzwjf86bEHOOv085zDhkP8ZP41nHHqdznp+NM6PadUsMPN0kSzzDSS6brtB2xps1wKHNmNYyjwiogo8AdVfagbdQ29ia/csQxkuNkG4POvNnPT/zyBZQlej4cFP3cyWf6/a8/hqlsf5K4/LmTa+IN3qjN13MGc8cO72Ly9ihuvOJ/Bew9l3yGDKCkt4bKfXAVAn4I+3PGzX2B1Evrh8EMO4zsnnsYlP74CgNmnnMHogw9j6/ayxBU0xl5FRRRPmMx3Lz6J446eyQ3zbmbQ3kMYNfIgZk0/KWG1Af334urv/YiLrvgOANd870cM6L8XAAse/h/Gjh7PzG+eyPKPP+D+B38LIkwunsbN/3EbAP9+7UU+XvkR9Q21ra6Zt91yD6MPHZvMJe42TtYpxV/fTNFA9+fNzSa6dK8UkXOBk+N5YRGRS4FpqnpdgrK3Aj5VvafNuuGquk1EhgKvAtep6tsJ6l4NXA1wwAEHTN60adPun5Vht4hteg+qvoK8vuluSsr55QPP0K9PAf8xxxHNsvq92FY/sFO3ypQRi0DhQCR/Z3fVpmAT51xyMk89+kKrqSXbaQ5GGHbAXux/0KB0NyXn2FP3ylJg/zbLI4BtyR5cVbfF/1cAC3FMQYnKPaSqU1R1ypAhQ5LdvSFFqKoT9sCT+WabVFAT6LvnSUaSJUE0yw+WvctZF36LC8+93DUiD475prbcZ6JZZhjJmG6WAYeIyChgK3ABcFEyOxeRvoClqo3xzycBt+1uYw09SFMNRMKu7M0D/OIH57V+DkVsgpE8vFYvRVwUCyIhJ/euOH2ro6Yex78Xvt87x+9FLFsINUdobgpT0Cc3UlBmA10KvapGRGQe8G/ABh5R1VUicm18+4MiMgxYDvQHYiJyPTAGGAwsjLteeYAnVfXlnjkVw56g9VvIicAsOElGpCeiVXaGgIaDSF5mz0/YU1rcLOuqAwwzQp8xJOVHp6qLgcXt1j3Y5vN2HJNOexqATtwRDBlDzdc5kwS8JtAPScdLLdwELhd6cHr1NeU+hu2/V7qbYohjQiAY0OZGCDbkRBLwaCxFScC7i1gQacqJGP8tScPDIZOMJFMwQm9AG7Y6H7LArXJPaWx23P56/1Q7niXrNnbMks2NePzZgBF6A9RsyIrZsKmgJpDOwWZFcySXrFhigpxlEEbocxyNNDvx57NgNuyeogr1PZFkJFnEdtws3W+9cWbJ1jQRjZpZspmAEfocxwliRk6YbXyhAjTJJCM9gzhvm5j7I4WI5Qx3++qC6W6KASP0htoNJIqI6EbqAn1QTfO5au6Yb0CN+SZDMEKfw2g0DPWl4MnsJCOpQBVq/P3w2Gky27QgFoRyxHzjtamt8JlcshmAEfpcprEl9rz7bwN/KJ+oWliSZtERCzSaE+YbyxJiMcVXZ7xv0o37n3BDh2jN1+SK2aa2qS+aKeeqioZzR/yqjfkm7Rihz1E0Goa6LblltpE0m21aEAvC/tww3+TZ1Fb4jfkmzRihz1Uat5MrZptAKJ9ozMKyMkRsxIJo1Alf7HJazDf+euN9k07c/5QbEpJLZpu6pj4Z2HlWNJIj5hs13jfpxgh9DuJ422zOGbNNtb8Ij2TYxB2xIOTLGfNNdbnxvkknRuhzEG3YulNsdDfjD+UTySSzTQtiQSyHvG+iMRqN903acP+TbtiVqq/Ila++xt8v/ZOkOiKXJk8JVJU1prsVOUtuPO2GVjQShMatOZEyMKZQHeiHx87QQU/LhlBueN94vTZ1lX6ikQwzoeUISQm9iJwiImtFZL2IzE+wfbSILBWRZhH5aXfqGnoXrdviGK5zwGzTGCwkltbYNl0hjgkt2pzuhvQ4YgmqSn1NjvyCyTC6fNpFxAYWAKfipAe8UETGtCtWA/wIuGc36hp6k6q1OROSuNrfLwv8ihQN5Yb4iSXGfJMmkunWTQPWq+oGVQ0BTwFnti2gqhWqugxoP7LUZV1D76EhH/grcyIkcTQm1DX1xZOukMTJ0hq62P32G4/XoqEmQDiU4d+JC0lG6PcDtrRZLo2vS4ak64rI1SKyXESWV1ZWJrl7Q3fQmo3OhxwISVzX1AeQLDhVARQi7p9Q1JJ5qrbS+NT3NskIfaJHJdnuR9J1VfUhVZ2iqlOGDBmS5O4NyaKqULkmJ3rzAJW+AQjZM/CnodwQP8sWKkobnPvR0GskI/SlwP5tlkcA25Lc/57UNaQSf4UTHtd2fwLwYNiDP5Tf+wnAdxexIRx0/Opdju2xCAZCNPncnzs3k0hG6JcBh4jIKBHJAy4AFiW5/z2pa0ghWrk23U3oNar9RU4Yn4w32+xMLvjUt5hvzKBs79Kl+4WqRkRkHvBvwAYeUdVVInJtfPuDIjIMWA70B2Iicj0wRlUbEtXtqZMxJEYjzVC7EbyF6W5Kj6MKVb6i9CcY6S5iQXMj5BW5PgSRJ8+mqqyB/Q7aG9t2v5tvJpCUn52qLgYWt1v3YJvP23HMMknVNfQuWrc5Z0IeNAQLiaqN18rQSVId0RISIdrs+slsliVEYkp9tZ+9hxaluzk5gfuf/BxHVaFiFVjut80DVPr6k81TTTXkT3cTegWxnEFZQ+9ghN7tBKqhqRbsvHS3pMdpjnioD/bJfN/5jrAsx6c+liWDyHuAx2vRWB+kyW8GZXsDI/QuRytWO6OS2TYyuRtU+bJzEHYHEg905v5evYggQMVW06vvDYzQuxgNN0HtBvC4fxA2FhMqff2zbxC2PS2DstlrfUqalkHZSCTLv7MswAi9i9Hq9TkTwKy2qQ9RtbAkyxUyPiibC9mnLEvQmFKzPTcmi6UT9ytAjqIag/IvcmImrCpsb9gLK9OySO0uIk6vPgewPRbbt9SZmbI9jBF6l6J1m534KTkwE9YfyicYycN2jdDbEGmGaA5kn7KFUDBCQ637f8GkEyP0LkRVoeyTnAlHXNawF4Jm8SBsAkTQZvcPVIoIYgnbSmpNr74HMULvRhrLoKk+J8w2TSEvDdnsUtkRYkGoCWJZNvFrN/B4Lfz1QQKN7k/Aki6M0LsQLVvpCIWruriJ2d64F5LVLpUdISCK5oCtXkRAYFtJXbqb4lqM0LsM9VeBrwI8BeluSo/THPFQm8k5YfcUieeUzYGolt48m/pqv5lA1UMYoXcZun2lExTLfV3cXShvGIC6zTa/E/EJVDkQq74lqmXZZtOr7wmM0LsIDVRD3Rbw9El3U3qc5oiHKn8R3h6wzZ962Wxq6/dccEKhED+78+ec8b1zuOTHV7B1e1nCcj+45XrO+/6lfPfqi7j9/t8QjbY5J8vm1SWLmHjMSFat+WyP25TJePNtarY3ml59D2CE3iWoKlq6LGds82UNewGpPVVVJZbCODML//0v+vcr4l9/fpZLZl/AfY8sSFjutzffwTO/f5zn/vBXauvrePWd11u3+QMB/vbPZxl/+PiUtStTaenVb91Yk+aWuA8j9G7BX+F42+RAuINg2EuNvx9/e/4Jzr7mYs6+5mKeWPgUAPf+aQFP/+u51rK/f/yPPPbckwA8+vcnuOi6uZx77SU88PjDAGzdXsbsqy7gjt/dzQXzLmd7ZflOx7r+lzdx4bw5fPfqi3h28fOt648+63j++6H7ueCHl3P1/HnU1NXu0s43l77DGSd8G4ATvjmLj1YuT+hC2K9vXwAi0SjhSLhV8AAWPPYQc867hDyPDZoDtvp8m7pKv/HASTFG6F2A05v/yBm8y4He/Nb6gaz56ksWvfoCT9z3Rx6/92H+8dIivly/llNmnMArby9pLfvqO69x4jeP5/0VH7J5Wyl/vf9PPP3AY6z56ktWfP4JACWlmznjW6fy9ILHGL7Pvjsd65c33MzffvcoT/7fI/ztn3+nrqEegKZgE6MPPoynFvyFyeMn8Ye//mmXdlZUVzJsyD4AeGwP/fr2a63fnu/ffD3HX/Bt+hT24YTjZgHw5fq1lFdWMP3I4wDQ5twIdobAlq+r090UV5GU0IvIKSKyVkTWi8j8BNtFRO6Pb/9MRI5os61ERD4XkZUisjyVjTc4aP0W8FflhKeNP5RHXVNfPl/zMccfM4PCgkL6FPbhW8fO4OMvPmX0wYdRU1dLRXUlazd8RVG/IvYdOowPPv6QpSs+5PwfXs4F8+ZQsmUTm7duAWDfocOYcPi4hMd78vm/c973L+Wy66+ivLK8tY5lWZw841sAnHb8KXyyalf7eaLeu3TwIv79nfey5Ml/EQ6H+ejTFcRiMe7+w33ccNWPWmpCuCknZst682waaptoNLNlU0aXUydFxAYWACfiJPteJiKLVHV1m2KnAofE/44Efh//38IsVa1KWasNrWgsAps/cBKLuLw3rwqbawYjKNpJeMcTjpvFknfeoKq2mlNmnNha94rzL+Oc02bvVHbr9jIKCxKbu5Z9+jEffrKMv/zvwxQWFHDFjT+gOZR4oDDRpd9n8FC2V5azz5ChRKIRfH4fA4r6d9ju/Lx8Zhx1HG8ufZtxh43h600buPJnPwCguraG62+bz723/TfjJs10dbpBEcG2hZK1lYydtj+W5eKT7SWS6dFPA9ar6gZVDQFPAWe2K3Mm8Jg6fADsJSL7tt+RIfVoxWoI+12ffg6gxt+PQDgfjxVl8rhJvPH+2zQFgzQFm3j9/bc5YtxEAE6ecSIvv/UqS959o9UMcvTkI3n+lRcINDkJuMurKqip63zQzxfw0b+oiMKCAjZuKeHzL3ekO47FYix55w0AXnrjFSaNnbhL/RlHHce/ljhZNJe88wZTJ07epUcfaApQWe30gSLRCO8uW8qo/b9BUd9+vPnMy7z02EJeemwh40eP5d5b72bsgQflRGRLj9emuSlM5Tb3h4HoDZIJhrIfsKXNcik799Y7KrMfUIYTWfsVEVHgD6r6UKKDiMjVwNUABxxwQFKNz3U05HNi2uTAAGwkZrGlbhC2RBGBww85jO+ceBqX/PgKAGafcgajDz4MgINHHkigKcDQQUMYMmgwAMdMPpKNW0q47CdXAdCnoA93/OwXWJbd4TGPnXwUz764kHOvvYRvjDiA8aPHtm4rLCjk600buHDeHPr17cdv/vNXu9SffcoZ3PLbX3LG986hf1H/ncqc94PLeOaBx2gKBvnxrT8jHA4RjcWYVjx5l18dOyGWkzHMU+D6X3Aer8XWr6vZe2hfvHm5Ebepp5CuAgmJyLnAyap6ZXz5UmCaql7XpsyLwK9V9d348mvAz1R1hYgMV9VtIjIUeBW4TlXf7uyYU6ZM0eXLjTm/M1QV3fAG1G2GvL7pbk6Ps6V2byp8A8jLkFmwR591PEuff73rgj1BLAIF/ZGCAek5fi8SCkbYe1gRBx4+NN1NyXhEZIWqTkm0LRnTTSmwf5vlEcC2ZMuoasv/CmAhjinIsIdo3Sao2wRe90+O8jXnU+EbgNfKDJFPO5YNwYbcGJjNt6kua6S+2v0eRz1JMkK/DDhEREaJSB5wAbCoXZlFwGVx75ujgHpVLRORviJSBCAifYGTgC9S2P6cRMNNsOk9Jzqly3++R2PCxuqhWBLLqFNNW28ecAKeiTMT2uWhfUUE22uxcU0lkbD75xH0FF0avlQ1IiLzgH8DNvCIqq4SkWvj2x8EFgPfBtYDAeB78er7AAvjA1Ae4ElVfTnlZ5FDqCq6+X2nN5cDJptt9QMJRTzkeUxvfifEhmgYbW5wvQnH47EIBSNsWV/NKGPC2S2SGuFQ1cU4Yt523YNtPivwwwT1NgC7uiMYdhutWR832fRLd1N6nPqmQioaB+DNELt8xhE34ainAHG515U336aqrJEBg/qw91D33/upxsyMzSI0UAOb3ne8bDLJjtEDNEc8bKweim1F3X6qe4A4XjiBKteHMhYRPF6LjWsqCAbcPzaRaozQZwkaaUa/XgKI61MExhQ2VA0lpoJtudsGvceI5YQyDlTTyRwyV2B7LFBY/3kZ0ahL8gP3EkboswDVGFryjpOEwutun3lV2FI7iEAo333pAXuKeDJxDda5Xuy9+TZBf5hNX1aaHLPdwAh9hqOq6JYPoX4zeN0/+FreMIAqX3+8dsSYbLqDZUNzY04kKfEW2FSX+9i6wYQzThYj9BmOln8BFWsckXe58lX7+7K1fm88RuR3A3HEvqnWcb91MSJCXr5N2aY6KraaEAnJYIQ+g4lVfglblzmTosTdX1VtoA8lNUPx2FFMDKvdpWVwthoNB9PdmB5FLMGbb7NpbSVV292fQH1Pcbd6ZDGx8tWw+X0nLWAn8VjcQI2/Lxuq98GWKJYYu+seIRYgEKhyfc/esgRvnuOJY4KfdY4R+gxDVYmVfQqlH+SEyFf5+rGxZigeiRoPm1TRIvb+KjQUSHdrehTLtvB6LUq+rKS8tM4M0HaAu/30sgyNRdBN70PNeteLvCqU1u9NeeMAvJbpyaecltzBgWo0FkbyB7g2hr1lW3jzhc3rqmnyhzngkMEmhn07S96j9gAACQ5JREFUjNBnCNrsQze+Af5KZ9ari0cjw1GLkuohNDT3Ic8yA689h+yYPRsNI332du1Yj2UJeQU2lVsbCAbCHDhmKHn5Rt5acOe3nkWoKrGajejqhRCodr3I1zcVsnr7CBqaC/Eake8F4hPswk1o43Y04t6k2yKO2Pvqmlj10RZqK93vapos5pWXRjTkQ7d8BHUlYBe4OudrOGqxrX4gVf7+2BIjzzaToXoVywMaA18Fml+EFPR3Ze/eEXsP0UiM9Z+XM2iYnxEHDcr53n1un32a0GgYrVwNZSud+f4u7sXHFKp9RZTW700sZpmJUOmkxW7f3IiG/FC4F+Lt48p7z/ZYWLZQU+6jttLPfgfuzZDh/bFt973cksEIfS+i0RBauQ62r3TCDHsKwePOAddYTKgO9KOsfiDhmI1HonhMqOEMQHb07gM1qFUPBQNcKfgtvftYNEbp+mrKSmoZPmpvBg8rcuLm5BBG6HsYVXVmK1Z9CdXrnSiDngLIc2dY2WDYS5W/iCpfEVG1sCWaMen/DG0QKx4QLQaBWlTqIL8fktfXdUHzLNsiz7aIRmNs+aqK0q9rGDK8iCHD+1PQx7tLwnY34q5vNENoFff6LVC9Dprjg0Iu7MGrQnPES11TIdX+IoKRPEDxWDHyTOq/zEesuNulOt45wQawvZDXD/EUOKLvEh20bQvbtojFlIrSespL6yno42XwsCL2GtzX1aKflNCLyCnAfTgZpv6oqne12y7x7d/GyTA1R1U/TqauG9BYFIJ1TqjYhm3QUOokcFYFO89VcWpUIRjxEmjOp7G5gPpgH6IxCwVsUeNJk7W0CX8dizodFYmHVPAUOolN7DxXCL/jiulBVYmEopRuqGHrxhps22LAoD4UDSykb1E+BX3yXOOP36XQi4gNLABOxEkCvkxEFqnq6jbFTgUOif8dCfweODLJulmBagzCTRAOQMiPNjdCU42T9CHY6Nz8qs6DYeeDN3tNM6oQjtmEozbhqIfmsIemcB6BcD7BsHdHOcBjxfAaDxp30WLWAUCd+z3k3yHwlscRfTsPsey4+FtOuOQs0kURwfZIq70+Fo1RU+GjptzX2jErKPTSpyiPPv3yyS/0kpdv48334PHaWfUSSKZHPw1YH08LiIg8BZwJtBXrM4HH4ikFPxCRvURkX2BkEnVThqpCsG5Hbxp1bJAa/9+6HINYDDTq/MWcP42FIRyEaAgiQYg07/ivUad+ywOgMefGtryJ3SJ7IeNPMOJ1etPq3HAKqAoaf9pUxTl1hJgKqkJMLWLqLEdiFpGYTSRmE43/j8RsYirOzlp+0qsgolgSw7bi16HtqZpJre6mvRtmNAoRP+BPUFZ2vCjEBssCLLAsJJ7UfMfbQNr80pX46rbiGf9sWb0zbiCCx7uzaTXUHKYpEKKqrBFEENTpz4kgAh6vjSfv/7d3dyFS1WEcx7+/2dUiQ7owKHwpg6CshGIKQ6jACCvxthd6wQgRShSyMCXoPqiELkKsCPLOiiSsNOiiG6XVshAzJXpZM9qIXi4Cmfbp4pxZR3dGD47nnN3//D6w7Mw553/2eYadZ/5z5sx5hrLfMxoMDzeyF4Gh7KyfRkNI+e/8tpSPb2giZSnLf2ioUcohpCKP3lzg5477o2Sz9nNtM7fgWAAkrQZWAyxYsKBAWF20/iWO7mZy94XIF+XLJ14E2svi1LKuwZHN0ntNV6L6Ge14iKNjV2ZF/YywoyPObJVOpdte37M4B8Mahx7XnckuJTJ9ZjJWBvU+Bz/a/2jj+U/HqtPG99hvr783a049hz/zWf8kAePj47Ra/9FqTX7+qyPHdhFXe02X17KJmxI33Dpv0gtOv4oU+m6P7plVoNc2RcZmCyO2AlsBms3mec0RNeMStPiB8xk67TRw13UzK6ZIoR8F5nfcnwf8UnCbmQXGmplZiYp8a+AL4FpJCyXNBB4Edp6xzU7gMWWWAH9FxImCY83MrETnnNFHREvS08AnZKdIvhkRhyStyde/DuwiO7XyGNnplavONraUTMzMrCtNxQv1N5vNGBkZqTsMM7NpQ9L+iGh2WzdYF3wwMxtALvRmZolzoTczS5wLvZlZ4qbkh7GSxoAfz3P4HOD3CxjOdOCcB8cg5u2ci7kqIi7vtmJKFvp+SBrp9clzqpzz4BjEvJ1z/3zoxswscS70ZmaJS7HQb607gBo458ExiHk75z4ld4zezMxOl+KM3szMOrjQm5klLulCL2mDpJA0p+5YyibpJUnfSvpa0vuSLqs7prJIWi7piKRjkjbWHU/ZJM2X9Jmkw5IOSVpXd0xVkTQk6UtJH9YdS1XyVqw78ufzYUm397vPZAu9pPlkTcl/qjuWiuwBboyIxcB3wPM1x1OKjobz9wKLgIckLao3qtK1gGci4npgCfDUAOTctg44XHcQFdsCfBwR15E1kus7/2QLPfAK8Bxn6Y6akojYHRGt/O5esm5eKZpoVh8RJ4F2w/lkRcSJiDiQ3/6H7Ik/t96oyidpHnA/sK3uWKoiaTZwB/AGQEScjIg/+91vkoVe0krgeEQcrDuWmjwBfFR3ECXp1Yh+IEi6GrgZ2FdvJJV4lWyyNn6uDRNyDTAGvJUfstomaVa/Oy3SM3ZKkvQpcEWXVZuBTcA91UZUvrPlHBEf5NtsJnurv73K2CpUuOF8aiRdCrwLrI+Iv+uOp0ySVgC/RcR+SXfVHU+FhoFbgLURsU/SFmAj8EK/O52WIuLubssl3QQsBA5KguwQxgFJt0XErxWGeMH1yrlN0uPACmBZpPsFiSLN6pMjaQZZkd8eEe/VHU8FlgIrJd0HXAzMlvRORDxSc1xlGwVGI6L9jm0HWaHvS/JfmJL0A9CMiKSvfidpOfAycGdEjNUdT1kkDZN92LwMOE7WgP7hlHsRK5uxvA38ERHr646navmMfkNErKg7lipI+hx4MiKOSHoRmBURz/azz2k7o7dJXgMuAvbk72T2RsSaekO68Aa04fxS4FHgG0lf5cs2RcSuGmOy8qwFtkuaCXwPrOp3h8nP6M3MBl2SZ92YmdkpLvRmZolzoTczS5wLvZlZ4lzozcwS50JvZpY4F3ozs8T9D5+3BUeMwRO1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZgU1fW/31PVPRsMO4iIBnfCvrtABMQFXGKIRo0rMW5J1CQmRqJZTKLGRH/5qgmGaGJMYoxbgiEGjaLiCgooooAIwgADA7Mv3dM9vZ3fH9UzzD49MDPdXX3f55lnuqrurbrVXfWpW+eee46oKgaDwWBwL1ayG2AwGAyG7sUIvcFgMLgcI/QGg8HgcozQGwwGg8sxQm8wGAwux5PsBrTGoEGDdMSIEcluhsFgMKQNa9euLVXVwa1tS0mhHzFiBGvWrEl2MwwGgyFtEJEdbW0zphuDwWBwOUboDQaDweUYoTcYDAaXk5I2eoPB4A7C4TCFhYUEg8FkN8U15OTkMHz4cLxeb8J1EhJ6EZkLPADYwB9V9Z42yk0FVgEXqeqznalrMBjcR2FhIfn5+YwYMQIRSXZz0h5VpaysjMLCQo488siE63VouhERG1gEzANGAV8VkVFtlPsV8L/O1u1KysrKmDBhAhMmTGDo0KEcdthhTJgwgX79+jFqVNcfesWKFZxzzjmdqjNr1qxWvYoee+wxbrjhhnbrxmIxbrrpJsaMGcPYsWOZOnUq27dvB+Css86isrKyU21pjYKCAsaMGXPQ+zEYgsEgAwcONCLfRYgIAwcO7PQbUiI9+mnAVlXdFj/Qk8B5wMZm5W4E/glMPYC6XcbAgQNZt24dAHfccQe9e/fm+9//PgUFBQkJciQSweNJXYvWU089xZ49e1i/fj2WZVFYWEivXr0AWLZsWZJbZzC0xIh813Ig32cig7GHAbsaLRfG1zU+8GHAfGBxZ+v2JNFolGuuuYbRo0dzxhlnEAgEAKeHfdtttzFz5kweeOABSkpKOP/885k6dSpTp07l7bffBuD1119veFuYOHEiNTU1APh8Pi644AJGjhzJpZdeSn3o51deeYWJEycyduxYrrrqKurq6lq06c9//jPHHXccM2fObDhOexQVFXHooYdiWc5PN3z4cPr37w848w9KS0spKChg5MiRXH311YwZM4ZLL72U5cuXM336dI499ljee+89wHkQXn755Zx66qkce+yxPPLII61+Z7fccgtTp05l3Lhx/OEPf+js124wGJJMIl3X1h4fzYPY3w/cqqrRZk+bROo6BUWuBa4FOOKIIxJoVufZsmUL//jHP3jkkUe48MIL+ec//8lll10GQGVlJa+//joAl1xyCd/97neZMWMGO3fu5Mwzz2TTpk3cd999LFq0iOnTp+Pz+cjJyQHggw8+YMOGDQwbNozp06fz9ttvM2XKFBYsWMArr7zCcccdxxVXXMHvf/97vvOd7zS0p6ioiJ/+9KesXbuWvn37Mnv2bCZOnAjA0qVLWbNmDT//+c+bnMOFF17IjBkzePPNN5kzZw6XXXZZQ53GbN26lWeeeYaHH36YqVOn8sQTT/DWW2+xdOlS7r77bp577jkA1q9fz6pVq/D7/UycOJGzzz67yX7+9Kc/0bdvX1avXk1dXR3Tp0/njDPO6JR9sCvR2nJ074cQqnVW5PZDho5DsvOT0p7uJFgbomhHJcHaMADZuR4O/Vx/cntlJbllB07s439CsKrrdpjTF2vM+V23vzYoKCjgnXfe4ZJLLul0vXPOOYePP/64m1qWGIkIfSFweKPl4cCeZmWmAE/GRX4QcJaIRBKsC4CqPgw8DDBlypRuyYZy5JFHMmHCBAAmT55MQUFBw7aLLrqo4fPy5cvZuHG/dam6upqamhqmT5/OzTffzKWXXsqXv/xlhg8fDsC0adMaPk+YMIGCggLy8/M58sgjOe644wC48sorWbRoUROhf/fdd5k1axaDBw9uaMOnn34KwBe/+EW++MUvtjiH4cOHs3nzZl599VVeffVV5syZwzPPPMOcOXNanOvYsWMBGD16NHPmzEFEGDt2bJPzPu+888jNzSU3N5fZs2fz3nvvNXxHAC+99BLr16/n2WefBaCqqootW7b0uNBruBbd+S5UFoAISPzS9ZegZVvQQccjh01B7MQ9EVKVSCRK4ZYySvc6b4yW7fSX/NV1lO/zMeCQ3hxx7CA8XjuZzTwwglXQlQ/lrnxotENBQQFPPPFEq0Kf6uZeSEzoVwPHisiRwG7gYqDJ2apqw10vIo8Bz6vqcyLi6ahuT5Kdnd3w2bbtBtMN0GDnBmfAc+XKleTm5japv3DhQs4++2yWLVvGiSeeyPLly1vdbyQSIdHMXQdib8vOzmbevHnMmzePQw45hOeee66F0Dduk2VZDcuWZRGJRNo8fvNlVeW3v/0tZ555Zqfb2VVoOIB++gIEq8HbyxH6BrJBY1CyCQ1WwTGnIVZq33TtEY3E2PJhEb6qOrJy7Ka/hzfudbG3hoA/xPEThqWn2CeBv/71r9x3332ICOPGjePOO+/kqquuoqSkhMGDB/PnP/+ZI444ggULFtCnTx/WrFnD3r17+fWvf80FF1zAwoUL2bRpExMmTODKK6+kf//+/Pe//yUYDOL3+3nllVf4wQ9+wAsvvICI8KMf/ahJ5zHZdGijV9UIcAOON80m4GlV3SAi14vI9QdS9+Cb3b2cccYZ/O53v2tYrh/c/eyzzxg7diy33norU6ZM4ZNPPmlzHyNHjqSgoICtW7cC8Le//Y2ZM2c2KXPCCSewYsUKysrKCIfDPPPMMx227f3332fPHuelKBaLsX79ej73uc91+hzr+fe//00wGKSsrIwVK1YwderUJtvPPPNMfv/73xMOO+aDTz/9FL/ff8DH6ywaCaJbXoS6Gsjq3Uzk44gF3t5QU4R+9ioai/ZY+7qSaDTGlvVF+KtbEfk4IkJWjoeAL8SnHxYRiaTnufYkGzZs4K677uLVV1/lww8/5IEHHuCGG27giiuuYP369Vx66aXcdNNNDeWLiop46623eP7551m4cCEA99xzD1/4whdYt24d3/3udwFYuXIlf/nLX3j11Vf517/+xbp16/jwww9Zvnw5t9xyC0VFRUk539ZIaGasqi5T1eNU9WhVvSu+brGqNh98RVUX1PvQt1U31XnwwQdZs2YN48aNY9SoUSxe7Jzm/fffz5gxYxg/fjy5ubnMmzevzX3k5OTw5z//ma985SuMHTsWy7K4/vqmz8VDDz2UO+64g5NOOonTTjuNSZMmNWxbunQpP/nJT1rst7i4mHPPPZcxY8Ywbtw4PB5Phy6Z7TFt2jTOPvtsTjzxRH784x8zbNiwJtuvvvpqRo0axaRJkxgzZgzXXXddkzeC7kQ1hm5dDoFKpyffHiJOmepCdMc7PdK+rkRV2b6xmJqqIN7s1kW+HhHBm21TW1PHZx/vS/jtMVN59dVXueCCCxg0aBAAAwYMYOXKlQ1mmMsvv5y33nqrofyXvvQlLMti1KhR7Nu3r839nn766QwYMACAt956i69+9avYts0hhxzCzJkzWb16dTeeVedI33fcBLjjjjsaPo8YMaLJgMj3v//9hs8rVqxoUm/QoEE89dRTLfb329/+tsW6WbNmMWvWrIblxm8Cc+bM4YMPPmhRp/Hxvva1r/G1r32tRZm2bPRz585l7ty5LdYDDbb3QYMGNTnXxx57rOFz8+/huOOO4+GHH26yn8ZlLMvi7rvv5u677271mN2JFm8Cf0nHIl+PiNOzL99CrP8IrH6Hd1wnRSgv9lFR4m+zJ9+cerGvrghQUlTDkGF9eqCV6YmqdvidNt7e2OzZ3kO0sbk31R+2JtaNISXRYBXsXgOe3NbNNW0hAnY27HgLjbR0Z01FwqEIOzaX4vFanRqzERE8XotdW0oJBXvmLSsdmTNnDk8//TRlZWUAlJeXc/LJJ/Pkk08C8Pe//50ZM2a0u4/8/PwGd+rWOOWUU3jqqaeIRqOUlJTwxhtvMG3atK47iYPE1T16Q/s0fuNJJVRjaMGbgIJ1AIONdhaEfGjhe8iIL3R187oUVWXHp6XEojGycjp/O9q2RTQcoWBzMceOOzT1Jyfl9O1y98qOGD16NLfffjszZ87Etm0mTpzIgw8+yFVXXcW9997bMBjbHvVm0vHjx7NgwYKGuSv1zJ8/n5UrVzJ+/HhEhF//+tcMHTq0iYdbMpFUfOWYMmWKmsQjmUusbBsUvN6Kh00nUIWwHxl5DtKr1aQ7KUF1RYBP1+3p0C7fHqpKqC7KMWOH0n9QgmauHmLTpk18/vOfT3YzXEdr36uIrFXVKa2VN6YbQ0qhGoM9a5xe+cH0TkVALHT32q5rXBejqhR+VoZYclA9cRHBtoXdn5WnvK3YkByM0BtSCi3fDiG/I/QHiyfXcbn0FR/8vrqBmooAtTV1eLwHfxvaHotgbYjK0p5zfTWkD0boDSmDxqLOAGxXiDzs79XvWZtyPV1VZddn5Qfdm69HRLBsodD06g2tYITekDJoxXYI13ad0EO8V7/XcdNMIarKa6n1dU1vvh7bYxEMhKkoMb16Q1OM0BtSAlWFPe93rcjD/l590bqu3e9BUlRQgdVFvfl66m31e7ZXmF69oQlG6A2pQU2RE5Gyq4UenF599R60rm0/6J4k4Avhr+7a3nw99bb62pr0mENg6BmMH70hJdDijQfnZdMeIoCipZuRw1r1PutRinc7fuTd4fPu7FPYV1jFUaNyunz/B8tHq3Y2hF3uCnLyvIw9se2w5pWVlTzxxBN885vf7JLj3X///Vx77bXk5eV1yf566jimR29IOhryQXWh0/PuLjw5UPIJGkvuDNJIJErp3ho8Wd0XddKbZVFe7CccSr3ZssHaMNm5ni776+ihUVlZyUMPPdRifTR6YMHg7r//fmpraw+objKP40qh37t3LxdffDFHH300o0aN4qyzzmqI894aB5L3tb19LViw4IDq3nHHHdx3333tlikrK2P27Nn07t273WBmH374ISeddBJjx47l3HPPpbq6GqAh3vyECRMYP348S5YsOaC2diVa8qkzwak7Z3VaHohG0Iqd3XeMBKgo9qMxxbK671zFEoiHM850Fi5cyGeffcaECROYOnUqs2fP5pJLLmnI1fD4448zbdo0JkyYwHXXXdfwAPjGN77BlClTGD16ND/96U8BJ9jhnj17mD17NrNnzwagd+/e3HrrrUyePJnTTjuN9957j1mzZnHUUUexdOlSoO0sbStWrGDWrFktstO1dpyDxXVCr6rMnz+fWbNm8dlnn7Fx40buvvvudqPQpRM5OTn84he/6PCBcPXVV3PPPffw0UcfMX/+fO69914AxowZw5o1a1i3bh0vvvhij0ajbA2NRaFkk9Pj7m4sG4o/StpApapStKMSy9P9t53tsdi7qyrjB2Xvuecejj76aNatW8e9997Le++9x1133cXGjRvZtGkTTz31FG+//Tbr1q3Dtm3+/ve/A3DXXXexZs0a1q9fz+uvv8769eu56aabGDZsGK+99hqvvfYaAH6/n1mzZrF27Vry8/P50Y9+xMsvv8ySJUsaos82ztK2evVqHnnkEbZv3w442enuv/9+Nm7cyLZt23j77bdbPc7B4jqhf+211/B6vU1CAk+YMIEvfOELqCq33HILY8aMYezYsU0iVLaV93Xt2rXMnDmTyZMnc+aZZzbEmJ41axa33nor06ZN47jjjuPNN98EICsri759nfgbbeWYbcxdd93F8ccfz2mnncbmzZs7PL9evXoxY8aMhjSGbbF582ZOOeUUwAmn+s9//hOAvLy8hmw4wWAw+bFRqndDLOz0uLsbOxtqyyFY0f3HagV/TR2hYBjb7v7v3PZYREJRaiqD3X6sdGLatGkN2dFeeeUV1q5dy9SpU5kwYQKvvPIK27ZtA+Dpp59m0qRJTJw4kQ0bNjTJONeYrKyshmiyY8eOZebMmXi93iaZ3F566SX++te/MmHCBE444QTKysrYsmVLQ3uGDx+OZVkN2em6A9cNxn788cdMnjy51W2NkwOUlpYyderUBjFsLe/rCSecwI033si///1vBg8ezFNPPcXtt9/Oo48+CjgpxN577z2WLVvGz372M5YvX87JJ5/MySefDNBmjtl61q5dy5NPPskHH3xAJBJh0qRJDW2vj4HfPIZ9oowZM4alS5dy3nnn8cwzz7Br1/4c7e+++y5XXXUVO3bs4G9/+1tS06Bp6WZaTy3cDcQHKrV8G3LYgJ45ZiPKimpQumcQtlUESvfW0Kd/N459pBnNQwtfeeWV/PKXv2xSZvv27dx3332sXr2a/v37s2DBAoLB1h+YXq+34fdsK5NbW1naVqxY0Wp2uu4goR69iMwVkc0islVEFray/TwRWS8i60RkjYjMaLStQEQ+qt/WlY3vLO0lB2jtybp582Y+/vhjTj/9dCZMmMCdd95JYWFhw/6+/OUvAy3zz9ZTn2P2wQcfpLKysoWgvvnmm8yfP5+8vDz69OnTJP789ddff8AiD/Doo4+yaNEiJk+eTE1NDVlZ+90WTzjhBDZs2MDq1av55S9/2eZF3N1oJOj06HvCbFOPJxtKtzgxdXqQWCxG2T4f3m4chG2O12tTsc9HNNKz55pKtBdeeM6cOTz77LMUFzshMsrLy9mxYwfV1dX06tWLvn37sm/fPl544YWE9tcWB5Kl7UCO0x4dduVExAYWAafjJPteLSJLVbXxu8wrwFJVVREZBzwNjGy0fbaqlnZZq9th9OjRDYmsm9OevbKtvK+jR49m5cqV7dZp60ncWo7ZkSNHNinTXb27kSNH8tJLLwHOhfXf//63RZnPf/7z9OrVi48//pgpU3re7VArd8YHYXvQgmh5IOwHXzHkD+2xw9ZUBIlFtVt859tCLEFVqSr3M2BIFybkPghy8rxd7l7ZHgMHDmT69OmMGTOG3NxcDjnkkIZto0aN4s477+SMM84gFovh9XpZtGgRJ554IhMnTmT06NEcddRRTJ8+vaHOtddey7x58zj00EMTtp9fffXVFBQUMGnSJFSVwYMH89xzz7Vb50CO0x6JvLNPA7aq6jYAEXkSOA9oEHpV9TUq3wtI2gjQqaeeym233cYjjzzCNddcA8Dq1aupra3llFNO4Q9/+ANXXnkl5eXlvPHGG9x7771t5n49/vjjKSkpYeXKlZx00kmEw2E+/fRTRo8enVBb6nPMjh07lpUrV/LJJ580EfpTTjmFBQsWsHDhQiKRCP/5z3+47rrrDv5LwEk5OGTIEGKxGHfeeWfD28H27ds5/PDD8Xg87Nixg82bNzNixIguOWanKdncM7b5VtDyrUgPCn1JUTUiPX9biCWU7KlJGaFvz+e9u3jiiSfa3HbRRRe1msS7cVa2xtx4443ceOONDcs+337pa57foX5bW1na2stO1/w4B0si3YvDgF2Nlgvj65ogIvNF5BPgv8BVjTYp8JKIrBWRa9s6iIhcGzf7rCkpOfC4JCLCkiVLePnllzn66KMZPXo0d9xxB8OGDWP+/PmMGzeO8ePHc+qppzYkB2iLrKwsnn32WW699VbGjx/PhAkTeOedxPORdpRjdtKkSVx00UVMmDCB888/ny98YX+SjMWLFzfY6ZszYsQIbr75Zh577DGGDx/eMFB09dVXUx/H/x//+AfHHXccI0eOZNiwYQ3pCt96662Gc5k/fz4PPfRQQy7NnkTraqC21Bkg7WnsHCjf1mM+9ZFwlKrS2m71nW8Lj9eipiKYkj71hp6jw8QjIvIV4ExVvTq+fDkwTVVbfdyIyCnAT1T1tPjyMFXdIyJDgJeBG1X1jfaOaRKPuJ/Y3vWwey1k9U5OA0J+5OjZSL/PdfuhSvfWULCp+IAySHUFoWCEI44bxJDDOs7G1NWYxCPdQ3ckHikEGmdZHg7saatwXMSPFpFB8eU98f/FwBIcU5Ahg1FVKP20e+LaJIoIWvZZjxyqrKjGmcSUJCxbKN2TvMlTme7L39UcyPeZiNCvBo4VkSNFJAu4GFjauICIHCPxUUURmQRkAWUi0ktE8uPrewFnAB93upUGd1FX4/xZ7Q+kdSueHKja1e3mm0g4Sk1VoEcHYZtjeyxqfXWE6nrefJOTk0NZWZkR+y5CVSkrK+twHk1zOnyXVNWIiNwA/A+wgUdVdYOIXB/fvhg4H7hCRMJAALgo7oFzCLAk/gzwAE+o6oudaqHBdWhVPAxBMidriQWoEzWz7+EdFj9Qqsprga4NR9xZ6o9dVVbL4GF9evTYw4cPp7CwkIMZdzM0JScnh+HDh3eqTkJGQ1VdBixrtm5xo8+/An7VSr1twPhOtcjgfsq2Jrc3X4/iTJ7qRqEv2+tL6vOsHssWyvbW9LjQe73ehpmohuThuhAIhtRGQz4IVCTXPl+PJwcqdzrxdrqBaCRGdUVyzTb12B4LX1WQcKh7ztWQ2iT/CjRkFFoVn1mcGt1c0Cj4uifgXXVlAOjBkAft0GC+Ke/+ELuG1MMIvaFnKduatElSraLxXLXdQPleH0mcO9gCsYTyfb6OCxpchxF6Q4+h4YCTpDsZk6TawpMDFdu7PPZNLBqjssyP19vzk6TawuO1qK4IZHTsm0zFCL2hx9DqPU6gyhQwZTRg2RCNQG1Zl+7WV12HxjSp/vPNqTff1MRNSobMwQi9oeeoLCA1LzndP3bQRVSW+lPIaNMYpaKk/ciJBveRinedwYVoLNrzIYkTxc6Cim1dtjtVpaLYnxLeNs3xeG3nIWQmMGUUqXclGtyJvxg01rMhiRPF8kKw2nH97AKCtWHCoUi35oU9UCxLiEZi1NaEkt0UQw+SgnedwY1o5S4n9nwqErdda3WbIZw6Rb0LYyq4VbZFVbkx32QSRugN3Y6qOvb5VPK2aY7Y0EVulhXFPqweyAt7oFgei/JiI/SZhBF6Q/dTVw3h2tTyn2+OJxtqig46yFkkHMVfXYftSd1by7aFoD+UlCBnhuSQulejwTVo9e54ysDU7eU2jB0c5CzZmhSaDdsW9W2rNrNkMwYj9Ibup6IgtXvz9aiiVbs6LtcOFSV+Z65AiiOWGDfLDMIIvaFb0WjY8bhJhSBmHWFnQ8WOA3Y9VFUnZWAKzYZtC4/HmSUbi6XoALmhSzFCb+he/MXO/1R0q2yO5YFwAEIHlo0p4A8RjWpKulU2RyxBVamtqUt2Uww9QBrcfYZ0Rqt2p65bZXNEAEVr9h5Q9eqKQFpNRNKYiWaZKSQk9CIyV0Q2i8hWEVnYyvbzRGS9iKwTkTUiMiPRugaXk+pulc0R2xlTOAAqiv3YntTvzddjey0qjZ0+I+hQ6EXEBhYB84BRwFdFZFSzYq8A41V1AnAV8MdO1DW4FK3zQcifHgOx9Xiywbe308lIopEY/upgSrtVNse2hYA/ZJKRZACJXJXTgK2quk1VQ8CTwHmNC6iqT/e/s/ZifxDuDusa3IvWFDkfUtjVsAViOTaN2tJOVfNVBUFS262yOfVt9VUFk9wSQ3eTiNAfBjT2OSuMr2uCiMwXkU+A/+L06hOuG69/bdzss8YkEnYJVTscU0i6odrpcAhV5bWplGMkccSJtGlwN4kIfWtdlBaXtKouUdWRwJeAX3Smbrz+w6o6RVWnDB48OIFmGVIZ1RhU73FMIemGndWpcAiqTuhfOwWjVXaEx2OZaJYZQCJXZiFweKPl4UCb3R1VfQM4WkQGdbauwUXUlqVutMqOsLwQrHIyYiVAKBghXJea0So7wrItotEYwdpwspti6EYSuQtXA8eKyJEikgVcDCxtXEBEjpG4wU9EJgFZQFkidQ3uRGuK0setsjkizp+vOKHi6RD2oF3UZJ1yOx26Q6hqRERuAP4H2MCjqrpBRK6Pb18MnA9cISJhIABcFB+cbbVuN52LIZWo3OH0jNMVBa0uRPp/rsOiVWWBtBpvbo5YQmWJnyGH9U12UwzdREJ+b6q6DFjWbN3iRp9/Bfwq0boGd6PREPjLwJuX7KYcOJ5sqNqFqrbbU1dVqsprsdMg7EFbeDwWNZVBYrEYlpWGpjZDh5hf1dD1+EtSLwl4ZxEbwsEOwyEE/CFisfQIe9AW9QnMTdYp92KE3tDlaPWe9LXP1yMComhN+2GLayqDrvBYicXU2OldjBF6Q9dTuTO9wh60iQ1VO9stUVHix07hbFKJYntM2GI3Y4Te0KVouNbJKJVOYQ/aws6C6j3OnIBWiEVj+KrSK+xBW9gei9qaOiIREw7BjaT/FWpILXzF+90T0x3LBo1CoKLVzf54iN+0datshIiAgL/KhC12I0boDV2KVhWmZyiAttC2wxanW1jijlA1YYvdihF6Q5ehqlC9Kz3DHrSF5WnTTl9Z6g77fD0ej0VVmRF6N2KE3tB1hHwQrkvPQGZtYWeBr7hF2OJIJErAF3KFfb4eyxbqAmHCoUiym2LoYtxzlRqSjvr2AeoO+3w99bF6asuarPZX16VdWOKO2B+22Njp3YYRekPXUbUrPYOYdYTGWtjpq8tr036qQKsIxnzjQlx4VxqSgWOf3+0S//lmWF4ntn4jKktr8bjIbFOPHbfTu2mQ2WCE3tBVBCshFnFcEt2GnQW1ZWjUCeUbDkWpC4SxXDQQW49lCeFQhFCdsdO7CSP0hi5BffvSP+xBW9Tb4f1O5jN/dTC+2n1C32CnrzTpBd2EEXpD11C5013eNs1RbciB63obthh/erdhhN5w0KjGwLfXXf7zzbGznIcZUFlWm5ZpAxPFY+z0rsO9V6uh56gtT9+0gYlieSFYSZ3fn7ZpAxPFsi2iEaUuYNILuoWE7kwRmSsim0Vkq4gsbGX7pSKyPv73joiMb7StQEQ+EpF1IrKmKxtvSA3U72L7fD0iIBa+4pL4onuF3kGpMXZ619BhiEERsYFFwOk4yb5Xi8hSVd3YqNh2YKaqVojIPOBh4IRG22eramkXttuQSlTudEe0yo7QGFXF5SDuT7knllBVVsvgYX2S3RRDF5BIj34asFVVt6lqCHgSOK9xAVV9R1XrQ/ytAoZ3bTMNqYrGok7ESjsr2U3pdtTKproy7Er/+ebYHst1QdsymUSu2MOAXY2WC7+wToUAACAASURBVOPr2uLrwAuNlhV4SUTWisi1bVUSkWtFZI2IrCkpKUmgWYaUoLYMJ+yB+8UvpDlEIiCuCs/ZOpYlxGJK0G/s9G4gkfft1oyRrV7pIjIbR+hnNFo9XVX3iMgQ4GUR+URV32ixQ9WHcUw+TJkyxf13kktQ31732+fj+EI5KCCxENi5yW5Ot6Oq1FQFyO3t/rc1t5NIN6wQOLzR8nBgT/NCIjIO+CNwnqo2RIBS1T3x/8XAEhxTkMEtVO1yPFIygKpgHqKgkcwYpLTidnpD+pOI0K8GjhWRI0UkC7gYWNq4gIgcAfwLuFxVP220vpeI5Nd/Bs4APu6qxhuSi8YizmxR2/1Crwo1gVxsOwbhzBD6ejt9LJYZb2xupkPTjapGROQG4H+ADTyqqhtE5Pr49sXAT4CBwENxt7OIqk4BDgGWxNd5gCdU9cVuORNDz+MvxYnV6377fF3ES1QtvLY6MX1iUXfG9WmEZQkRVYL+EHn5Lp4MlwEk5BOnqsuAZc3WLW70+Wrg6lbqbQPGN19vcAeOfb71xNluw1eXvX9gSgSidWDlJbNJPYIq1FQFjNCnOe7vihm6j8qdmWWfr19Q0Awx31iWUFlq7PTpjhF6wwGh0ZAT+iAT/OcVaoK52FY8naBlQSToriTobeDxWPgqg8ZOn+YYoTccGP5Sx/HW9aEAIBjxElNhf3gbcUxW6v6Y7WI5swZqfSa9YDpjhN5wQGhNBvnPB3PQ5tNJVNFIKDkN6mlUTXz6NMcIveHAqNqZEWYbqLfPNxt0FnHMNxmAZQuVpf5kN8NwEBihN3QajdRBoCIjBmJVoaYuF4/VXOgtiAQywk5veyx81XXEopnhYeVGjNAbOo+/JB621/32+UA4C9XWTlWcp0AsA+z08ZP31xg7fbpihN7QabR6T+bY5+tyUG3jgaaKRjND/DSm1FQGkt0MwwFihN7Qeap2ZYx9vjKQhyVtmCxEMigcgvGnT2eM0Bs6hYYDUFedEfb5mDo9eru5fb4eyRx/ettjUVtTRzRi7PTpiBF6Q+fIIPt8bciZ9t/2qQqgEHN/zHaJ/+bGTp+eGKE3dAqtLsyIHix0YJ+vR9XxQsoEYkp1uTHfpCNG6A2do6owg+zzvbDaMtvUI+K4WWYAttei0sSnT0uM0BsSRkN+CPkzIhF4NCbUhrKx2xqIrUcsiNRlxFuOZQtBf4hIOJrsphg6iRF6Q+L49mVMfBt/KBvQBE41XiDq/nAI9XZ6X1VmeBq5CSP0hoTRqkIy5ZKpaS2+TVtkkD89qlRXZIapyk0kdNeKyFwR2SwiW0VkYSvbLxWR9fG/d0RkfKJ1DemBqkJ15tjnq4K9Ojbb1CMWhDND/GyPZeLepCEdCr2I2MAiYB4wCviqiIxqVmw7MFNVxwG/AB7uRF1DOlBX7ZgnMsA+H4lZBMJZbU+Uao6I891kwGxhyxZCwQihOveHfnATifTopwFbVXWbqoaAJ4HzGhdQ1XdUtSK+uAoYnmhdQ3qgvn0ZIWQA/rpsJCH7fD3GTm9IbRIR+sOAXY2WC+Pr2uLrwAudrSsi14rIGhFZU1JSkkCzDD1K1U4QdyfDrqc6mNux/3xzMsmfHqXKuFmmFYkIfWtXfKtdOxGZjSP0t3a2rqo+rKpTVHXK4MGDE2iWoadQjUF1UcbY5ysDvfanDUwUsSCcGeLn8VhUldU64zaGtCARoS8EDm+0PBzY07yQiIwD/gicp6plnalrSHFqy0GjYLm/Rx+K2ISiHizppIiJBdGwk2LQ5Vi2RSQcoy7g/tAPbiERoV8NHCsiR4pIFnAxsLRxARE5AvgXcLmqftqZuobUR32Zkzawpi6nk/b5Rog4k6cyAqXGpBdMGzp0oVDViIjcAPwPsIFHVXWDiFwf374Y+AkwEHgonqQgEjfDtFq3m87F0F1U7sgIbxtw0gYeMAoaCSLe3K5rUIoilpNecPCwPsluiiEBErp7VXUZsKzZusWNPl8NXJ1oXUP6oNGwE7HScxACmCaoQnUgr+2wxB0h4vjT5/RvfXTKRXg8FtUVAWIxxbJcfrIuIDOmORoOHH8JkBlhiYMRLzGVztvn6xELYlFnPMPliCWoQsCXKaaq9MYIvaFdtGZPRgwwQifDHrSJZI6bpQmHkDYYoTe0T+XOjHKrtDjIh5qQMWGLLVuoKDHhENIBI/SGNtFwLQSrMiRtoLSfNjBRxHLyyGaAk5JJL5g+GKE3tI2vOGPSBvrrOkobmCiZll4QEw4hDTBCb2gTrdyJ691H4lQHc7uuE66KRjJE/BSTdSoNMEJvaBVVdeLb2NnJbkqPUBnoRFjijsigcAi21wlbbMIhpDZG6A2tE6iAWCQjwh6EozbBiPfA3Sqb0xAOwf3iZ1lCuC5CKGjCFqcyRugNraI1RRkhVOCEPUgsbWAnyQA3y/hMeONmmeIYoTe0TmVBxoQ9qAzkdf1IhIJmiJtlfTgEQ+pihN7QgoawBxlgn98f9qCLZ7Na8fSCGfBS1DgcgiE1MUJvaIm/2PmfAW6VgXBWPOxBx2W/9eObqfbVJLhngViUb33vSqprqqiuqeKpf/6t0+1buuxZzr1wFudeOIuly55ttcyeokKuufESLrh8Ll//1kXsKy4C4JNPN3D5NfOZf+npXHD5XF5c/p8WdX/5m59y4pyDy+4ploBCbY37TVXpihF6Qwu0ahcZ0RUFqoM5HWaTUlVisRiLfvEb+vTOT2i/9XV+98uH6JPflxpfNU/9q3NCX1VdyeJHH+DxPz7H3//4bxY/+gDV1VUtyv3md3dz7rwv8+zfXuTar32bB37/awBycnK58ye/YcnfX+ah3/yFex/4OdU1++tv2LSemprqTrWpLWKqVJVnhqdROmKE3tAEVYWKArBzkt2UdvHXBjnnW79k4gW3MG7+93jqxXcAOGrutyitcMRrzYbPOPWqOwD42UNPc8Vtv+W0r/+M48+5iUeeXQ5ARW1vHv/XX7jkxqv4yvWX8dDfHgFg994i5l9zMXf97l4uvuFK9pbsY94V86moqgTgb//8B+dfdynnX3cpjy95svU6paXMu/g0KirLeeChX1G4ewcXXjmP3/zubm772Xd57Y2XGs7nh3d8mxVvvtzkHN9Z9TonTp1B3z796NOnLydOncHbq1a0+C4+K9jCCVOmAzBt8kkN+xlxxFF87vAjARgy+BAG9B9IRWU5ANFolN8supvvfuuHB/dDxPF4LMqLfV2yL0PXkxmjbYbEqatypvB7Uzss8Ytvr2PY4P48v8gRqqqajnuTH326k3cevwt/IMjkC2/lzOlTeO29fRTu2cXfH/wTqsq377iFtR99wNDBQyko3MnPbv4Rt99wS5P9bNzyCf9++Xkef+CPqCqXfftqpoydSH7vPs3qqDMIoMq3v3krW7d9ytN/cdIpr/lgFY8/+Sdmn3IGNb5q1n20ll/86P81OU5x6T6GDhnWsHzIkEMpLt3X4ryOP+bzLH/tBS696Cpeef1/+Gt9VFZV0K9v//3nvnEd4XCYww/7HABPPvsXZs04jcGDhiT2hXeAZQt1tWFCdRGyso2spBoJ9ehFZK6IbBaRrSKysJXtI0VkpYjUicj3m20rEJGPRGSdiKzpqoYbuget3gNoytvnxx57BK+s+oiF//c4b67dRN/8jh9MX5w9hdycLAb178OsqaN5c90OVr2/ipXvv8tF37qSi29YQMGuHezc7eSzP3TIUMZ9fkyL/Xyw4UNOPXkmuTm55OXmMWf6TN7/+MNW6sS/w2ioxT6mTDyRnYU7KCsv5YWXl3LarHl4PE0FsvVJSC1/l5tvuJ01697lwivPYu0HqxgyeCi2vX/+Q0lpMbf//GZ+fvu9WJZFcck+XnptGV+9YEGH31miNLhZGvNNStLho1dEbGARcDpODtjVIrJUVTc2KlYO3AR8qY3dzFbV0oNtrKEHqChIC7fK40YMY/VTv2LZm+9z+4NPcPpJ4/nx9Rfgsa0G749gXdN4M9JMJP2hHBTl6xddwQVnz2+ybffeInJzWs8U1d4s0NbqOOEQWn6n58ydz7KXnuPF5f/hZ7fd22L7IYOHsvqDVQ3L+4qLmDrxxBblhgw+hP/75R8AqK31s3zFi+T3djI/+fw13PD9r3HDtd9j3JhJgDNIu6uwgHMvnAlAMBjgnK/M5PlnXm/zvBJBLCea5aBDTdapVCORHv00YKuqblPVEPAkcF7jAqparKqrAfdHcnIxGg05gczSwK1yT3E5eTlZXHbOKdx85bm8v2kbACOGDWHtRufzv5avalJn6YrVBOtClFXW8PqaDYwYMZHpk6fx3EvPUxtweqL7Sospj9ux22LymIm89s4bBIJBAsEAr77zBpPGjG+7QjhAr9ze1NY29TU/7+wL+PtTfwbgmKOOa1Ht5BNnsvK9N6murqK6uoqV773JySfObFGuorKcWMwJ3/Cnvz7El8650DlsOMR3F17HufO+zBmnnt1Q/pTpp/Lq82t44V9v88K/3iYnJ/egRR4au1maaJapRiJdt8OAXY2WC4ETOnEMBV4SEQX+oKoPd6KuoSfx7XMsAylutgH4aMtObv3N41iW4PV4WPQjJ5Plj6+/gGvuWMw9f1zCtLHHNKkzdcwxnPute9i5t5Rbvn4RgwYM4dDBAykoLOCK714DQF5OHnf94KdY7YR++Pyxx/PF08/msm9/HYD5c89l5DHHs3tvUesVNEa//HwmjJvMly89gxknzeLmG25j4IDBHDniaGafckar1fr26ce1X7uJS77+RQCu+9pN9O3TD4BFj/yG0SPHMusLp7Pm/VU8uPjXIMLkCdO47Xs/B+B/r/yX99e9R1V1RYNr5s9vv4+Rx41O5CvuNE7WKcVfVUd+f/fnzU0npKNgRCLyFeDMeF5YRORyYJqq3thK2TsAn6re12jdMFXdIyJDgJeBG1X1jVbqXgtcC3DEEUdM3rFjx4GfleGAiO14G0q3QFavZDely/nZQ0/TOy+H7y1wRLOoqh97qvqT5emBGC2xCOT2R7J7N1kdCAa44LIzefKx5xtMLelOXTDC0CP6cfjRA5PdlIxDRNaq6pTWtiViuikEDm+0PBzYk+jBVXVP/H8xsATHFNRauYdVdYqqThk8eHCiuzd0EarqhD3wpL7Zpisor+118ElGEqWVaJarVr/Fl746h69+5UrXiDw45puKfT4TzTLFSMR0sxo4VkSOBHYDFwOXJLJzEekFWKpaE/98BvDzA22soRsJlEMk7MrePMBPv3lhw+dQxCYYycJr9VDERbEgEnJy74rTtzpx6gz+t+Sdnjl+D2LZQqguQl0gTE5eZqSgTAc6FHpVjYjIDcD/ABt4VFU3iMj18e2LRWQosAboA8RE5DvAKGAQsCTueuUBnlDVF7vnVAwHQ2bNhs1FuiNaZXsIaDiIZKX2/ISDpd7NsrKslqFG6FOGhPzoVHUZsKzZusWNPu/FMek0pxpoxx3BkDKUf5YxScDLa3sjyXiohQPgcqEHp1dfvs/H0MP7JbsphjgmBIIBrauBYHVGJAGPxrooCXhnEQsigYyI8V+fNDwcMslIUgUj9Aa0erfzIQ3cKg+WmjrH7a/nT7XtWbJuY/8s2cyIx58OGKE3QPm2tJgN2xWU1yZzsFnRDMklK5aYIGcphBH6DEcjdU78+TSYDXuwqEJVdyQZSRSxHTdL91tvnFmy5QGiUTNLNhUwQp/hOEHMyAizjS/kxJ5PJMlI9yDO0ybm/kghYjnD3b7KYLKbYsAIvaFiG61FRHQjlbV5HSYZ6XY0c8w3oMZ8kyIYoc9gNBqGqkLwpHaSka5AFcr9vfHYSTLb1CMWhDLEfOO1qSj2mVyyKYAR+kympj72vPsvA38om6haWJJk0RELNJoR5hvLEmIxxVdpvG+SjfvvcEObaPlnZIrZpiLQC02Vc1VFw5kjfmXGfJN0jNBnKBoNQ+WuzDLbSJLNNvWIBWF/Zphvsmwqiv3GfJNkjNBnKjV7yRSzTW0om2jMwrJSRGzEgmjUCV/scurNN/4q432TTNx/lxtaJZPMNpWBvBTsPCsayRDzjRrvm2RjhD4DcbxtdmaM2abMn49HUmzijlgQ8mWM+aZsn/G+SSZG6DMQrd7dJDa6m/GHsomkktmmHrEglkHeN9EYNcb7Jmm4/043tKR0C5ny05f7eyd/klRbZNLkKYHSoppktyJjyYy73dCARoJQszsjUgbGFMpqe+OxU3TQ07IhlBneN16vTWWJn2gkxUxoGUJCQi8ic0Vks4hsFZGFrWwfKSIrRaRORL7fmbqGnkUrdzmG6www29QEc4klNbZNR4hjQovWJbsh3Y5YgqpSVZ4hbzApRod3u4jYwCJgHk56wK+KyKhmxcqBm4D7DqCuoScp3ZwxIYnL/L3TwK9I0VBmiJ9YYsw3SSKRbt00YKuqblPVEPAkcF7jAqparKqrgeYjSx3WNfQcGvKBvyQjQhJHY0JloBeeZIUkTpSG0MXut994vBbV5bWEQyn+m7iQRIT+MGBXo+XC+LpESLiuiFwrImtEZE1JSUmCuzd0Bi3f7nzIgJDElYE8QNLgVAVQiLh/QlF95qmKEuNT39MkIvSt3SqJdj8SrquqD6vqFFWdMnjw4AR3b0gUVYWSTRnRmwco8fVFSJ+BPw1lhvhZtlBcWO1cj4YeIxGhLwQOb7Q8HNiT4P4Ppq6hK/EXO+FxbfcnAA+GPfhD2T2fAPxAERvCQcev3uXYHotgbYiAz/25c1OJRIR+NXCsiBwpIlnAxcDSBPd/MHUNXYiWbE52E3qMMn++E8Yn5c02TckEn/p6840ZlO1ZOnS/UNWIiNwA/A+wgUdVdYOIXB/fvlhEhgJrgD5ATES+A4xS1erW6nbXyRhaRyN1ULEdvLnJbkq3owqlvvzkJxjpLGJBXQ1k5bs+BJEny6a0qJrDjh6AbbvfzTcVSMjPTlWXAcuarVvc6PNeHLNMQnUNPYtW7syYkAfVwVyiauO1UnSSVFvUh0SI1rl+MptlCZGYUlXmZ8CQ/GQ3JyNw/52f4agqFG8Ay/22eYASXx/SeaqphvzJbkKPIJYzKGvoGYzQu53aMghUgJ2V7JZ0O3URD1XBvNT3nW8Ly3J86mNpMoh8EHi8FjVVQQJ+MyjbExihdzlavNEZlUy3kckDoNSXnoOw+5F4oDP39+pFBAGKd5tefU9ghN7FaDgAFdvA4/5B2FhMKPH1Sb9B2ObUD8qmr/UpYeoHZSORNP/N0gAj9C5Gy7ZmTACzikAeUbWwJM0VMj4omwnZpyxL0JhSvjczJoslE/crQIaiGoN9H2fETFhV2FvdDyvVskgdKCJOrz4DsD0We3dVmpmy3YwRepeilTud+CkZMBPWH8omGMnCdo3Q2xCpg2gGZJ+yhVAwQnWF+99gkokReheiqlD0QcaEIy6q7oegaTwI2woiaJ37BypFBLGEPQUVplffjRihdyM1RRCoygizTSDkpTqdXSrbQiwIBSCWZhO/DgCP18JfFaS2xv0JWJKFEXoXokXrHKFwVRe3dfbW9EPS2qWyLQRE0Qyw1YsICOwpqEx2U1yLEXqXof5S8BWDJyfZTel26iIeKlI5J+zBIvGcshkQ1dKbZVNV5jcTqLoJI/QuQ/euc4Jiua+L24J91X1Rt9nmmxCfQJUBserro1oW7TS9+u4gM0brMgStLYPKXeDtleymdDt1EQ+l/ny83WCbn3fFfJ747Z/p37ffQe0nFArxo/t+zqYtn9C3T19+9cM7OWzooS3KffP271BaXkYkGmXSmPH88Fvfx7Zt/v3Sf7n/T79j8MBBiOXl4guu5MtfvPig2pTKeLNtyvfWcOgR/cjt5f6QHT2JEXqXoKpo4eqMsc0XVTsi3JWnqqpd6vmx5H//oU/vfP7z52d5ccXLPPDoIn59250tyv36trvo3asXqsr377yNl998lbmzTgfgjFPm8MNvfAey85Hcg3vwpDr1vfrd28s5ZszQJLfGXRihdwv+YsfbJgN688Gwl3J/b/7x3OP8+6XnAZg/91wum38x9/9pEYcOGcpF554PwO//9kd65eVxxfmX8Ngzj/PSG68SDoeYPX0m37z8GnbvLeKGH3+XKeMns37TR/zfT37V5Fjf+dmt7CvZR10oxCVfupALzvoSACd96VQuOOtLrP5wLX3y87ln4S8Y0K9/k7orVr7J9Zd9HYDTvjCbex76f6hqg6DV07uX85tFolHCkXCL7Vg2hHyQne98djHebJvKEj+1NXXk5bvfa6ynMDZ6F+D05t9zBu8yoDe/u6o/m7Z8wtKXn+fxB/7I3+5/hH+9sJRPtm5m7szTeOmN5Q1lX37zFU7/wqm8s/Zddu4p5O8P/omnHvorm7Z8wtqPPgCgoHAn586Zx1OL/sqwQ5qaVn52823843eP8cRvH+Uf/36GyuoqAALBACOPOZ4nF/2FyWMn8oe//6lFO4vLShg6+BAAPLaH3r16N9Rvzjdu+w6nXnwWebl5nDZjdsP6V95awVeuv5zv33k7RTvdnyWs3gNn12dlyW6Kq0hI6EVkrohsFpGtIrKwle0iIg/Gt68XkUmNthWIyEcisk5E1nRl4w0OWrUL/KUZ4WnjD2VRGejFR5ve59STZ5Kbk0tebh5zps/k/Y8/ZOQxx1NeWUFxWQmbt20hv3c+hw4Zyqr332Xl2ne56FtXcvENCyjYtYOdu3cBcOiQoYz7/JhWj/fEc89w4Tcu54rvXMO+kn0NdSzL4syZcwA4+9S5fLBhfYu6rZmBWvTW4/z+7vtZ/sR/CIfDvPfhWgBmnjiDZX/5F88sfpwTJk3jx/f8KCNmy3qzbKorAtSY2bJdRoemGxGxgUXA6TjJvleLyFJV3dio2Dzg2PjfCcDv4//rma2qpV3WakMDGovAzlVOYhGX9+ZVYWf5IARF2wnveNqM2Sx/8zVKK8qYO/P0hrpfv+gKLjh7fpOyu/cWkZvTenTP1R++z7sfrOYv//cIuTk5fP2Wb1IXat39r7Wv/pBBQ9hbso9DBg8hEo3g8/vom9+nzXZnZ2Uz88QZrFj5BidNmka/Pn0btn157nk88KeH0EAF0muIq9MNigi2LRRsLmH0tMOxLBefbA+RSI9+GrBVVbepagh4EjivWZnzgL+qwyqgn4i0dC8wdDlavBHCftennwMo9/emNpyNx4oyecxEXnvnDQLBIIFggFffeYNJY8YDcObM03nx9ZdZ/tZrDWaQkyafwHMvPU9twEnAva+0mPLK8naP56v10Sc/n9ycHLbvKuCjT/anO47FYix/8zUAXnjtJSaOHt+i/swTZ/Cf5U4WzeVvvsbU8ZNb9OhrA7WUlDl9oEg0wlurV3Lk4Z8DaFgP8PqqNznyiBEQqcuIyJYer01dIEzJHveHgegJEhmMPQzY1Wi5kKa99bbKHAYU4UTWfklEFPiDqj7c2kFE5FrgWoAjjjgiocZnOhryOTFtMiDefCRmsatyILZEEYHPH3s8Xzz9bC77tjPYOX/uuYw85ngAjhlxFLWBWoYMHMzggYMAOHnyCWzfVcAV370GgLycPO76wU+x2hncnD75RJ797xK+cv1lfG74EYwdObphW25OLp/t2MZXb1hA7169+dUPf9Gi/vy553L7r3/GuV+7gD75fZqUufCbV/D0Q38lEAzy7Tt+QDgcIhqLMW3C5Ia3jn/8+2lWrHoLj23TJ78PP//ejxyvqkCFY6Zz+Rucx2ux+7MyBgzphTfL+I0cDNKRO5mIfAU4U1Wvji9fDkxT1Rsblfkv8EtVfSu+/ArwA1VdKyLDVHWPiAwBXgZuVNU32jvmlClTdM0aY85vD1VFt70GlTshy/2eNrsqBlDs60tWisyCPelLp7LyuVeTc/BYBHL6IDl9Oy6b5oSCEQYMzeeozw9JdlNSHhFZq6pTWtuWiOmmEDi80fJwYE+iZVS1/n8xsATHFGQ4SLRyB1TuAG9espvS7fjqsin29cVrpYbIJx3LhmB1ZgzMZtuUFdVQVeb+9IrdSSJCvxo4VkSOFJEs4GJgabMyS4Er4t43JwJVqlokIr1EJB9ARHoBZwAfd2H7MxINB2DH2050Spe/vkdjwvayIVgSS6lTTVpvHnACnokzE9rloX1FBNtrsX1TCZGw+2P+dBcdGr5UNSIiNwD/A2zgUVXdICLXx7cvBpYBZwFbgVrga/HqhwBL4gNQHuAJVX2xy88ig1BVdOc7Tm8uA0w2e6r6E4p4yPKY3nwTxIZoGK2rdr0Jx+OxCAUj7NpaxpHGhHNAJDTCoarLcMS88brFjT4r8K1W6m0DWrojGA4YLd8aN9n0TnZTup2qQC7FNX3xpohdPuWIm3DUk4O43OvKm21TWlRD34F5DBji/mu/qzEzY9MIrS2HHe84XjapZMfoBuoiHraXDcG2om4/1YNAHC+c2lLXhzIWETxei+2bignWun9soqsxQp8maKQO/Ww5IK5PERhT2FY6hJgKtuVuG/RBI5YTyri2jHbmkLkC22OBwtaPiohGXZIfuIcwQp8GqMbQgjedJBRed/vMq8KuioHUhrLdlx6wu4gnE9dgpevF3pttE/SH2fFJickx2wmM0Kc4qorueheqdmZEZMp91X0p9fXBa0eMyaYzWDbU1WREkhJvjk3ZPh+7t7U/s9mwHyP0KY7u+xiKNzki73LlK/P3YnfVADxG5A8AccQ+UOG437oYESEr26ZoRyXFu02IhEQwQp/CxEo+gd2rnUlR4u6fqqI2j4LyIXjsKCaG1YFSPzhbhoaDyW5MtyKW4M222bG5hNK97k+gfrC4Wz3SmNi+jbDzHfDkuT7ZRLm/F9vKDsGWKJYYu+tBIRYgUFvq+p69ZQneLMcTxwQ/ax8j9CmGqhIr+hAKV2WEyJf6erO9fAgeiRoPm66iXuz9pWioNtmt6VYs28LrtSj4pIR9hZVmgLYN3O2nl2ZoLILueAfKt7pezLtLqQAACUJJREFU5FWhsGoA+2r64rVMT77Lqc8dXFuGxsJIdl/XxrC3bAtvtrDz0zIC/jBHHDvIxLBvhhH6FEHrfOj218Bf4sx6dfFoZDhqUVA2mOq6PLIsM/Dafcj+2bPRMJI3wLVjPZYlZOXYlOyuJlgb5qhRQ8jKNvJWjzt/9TRCVYmVb0c3LoHaMteLfFUgl417h1Ndl4vXiHwPEJ9gFw6gNXvRSF2yG9RtiDhi76sMsOG9XVSUuN/VNFHMIy+JaMiH7noPKgvAznF1ztdw1GJPVX9K/X2wJUaWbSZD9SiWBzQGvmI0Ox/J6ePK3r0j9h6ikRhbP9rHwKF+hh89MON795l99klCo2G0ZCMUrXPm+7u4Fx9TKPPlU1g1gFjMMhOhkkm93b6uBg35Ibcf4s1z5bVneywsWyjf56OixM9hRw1g8LA+2Lb7Hm6JYIS+B9FoCC35FPauc8IMe3LB484B11hMKKvtTVFVf8IxG49E8ZhQwymA7O/d15ajVhXk9HWl4Nf37mPRGIVbyygqqGDYkQMYNDTfiZuTQRih72ZU1ZmtWPoJlG11ogx6ciDLnWFlg2Evpf58Sn35RNXClmjKpP8zNEKseEC0GNRWoFIJ2b2RrF6uC5pn2RZZtkU0GmPXllIKPytn8LB8Bg/rQ06et0XCdjfirl80RWgQ96pdUPYp1MUHhVzYg1eFuoiXykAuZf58gpEsQPFYMbJM6r/UR6y426U63jnBarC9kNUb8eQ4ou8SHbRtC9u2iMWU4sIq9hVWkZPnZdDQfPoN6uVq0U9I6EVkLvAAToapP6rqPc22S3z7WTgZphao6vuJ1HUDGotCsNIJFVu9B6oLnQTOqmBnuSpOjSoEI15q67KpqcuhKphHNGahgC1qPGnSlkbhr2NRp6Mi8ZAKnlwnsYmd5Qrhd1wxPagqkVCUwm3l7N5ejm1b9B2YR37/XHrlZ5OTl+Uaf/wOhV5EbGARcDpOEvDVIrJUVTc2KjYPODb+dwLwe+CEBOumBaoxCAcgXAshP1pXA4FyJ+lDsMa5+FWdG8POBm/6mmZUIRyzCUdtwlEPdWEPgXAWteFsgmHv/nKAx4rhNR407qLerAOAOtd7yL9f4C2PI/p2FmLZcfG3nHDJaaSLIoLtkQZ7fSwao7zYR/k+X0PHLCfXS15+Fnm9s8nO9ZKVbePN9uDx2mn1EEikRz8N2BpPC4iIPAmcBzQW6/OAv8ZTCq4SkX4icigwIoG6XYaqQrByf28adWyQGv/fsByDWAw06vzFnD+NhSEchGgIIkGI1O3/r1Gnfv0NoDHnwra8rbtF9kDGn2DE6/Sm1bngFFAVNH63qYpz6ggxFVSFmFrE1FmOxCwiMZtIzCYa/x+J2cRUnJ3Vv9KrIKJYEsO24t9D41M1k1rdTXM3zGgUIn7A30pZ2f+gEBssC7DAspB4UvP9TwNp9KYr8dX/v727C5GqDuM4/v3Nzm6RIV0YFL6UQVBWQjGFIVRghJV42wu9YIQIJQpZmBJ0H1RCFyFWBHlnRRJWGnTRjdJqWYiZEr2sGa1EJSRt0z5dnDPr6M7o4Hjm7P7n94FlZ845/7PPM+w8858zZ87TXDzz25VKbz43kKgOnn5odeyffzn59xjHj50ACRHZfE5CgurgANWhgez3YIVqtZK9CAxkZ/1UKkLKf+e3pXx8RRMpS1n+AwOVQg4hdfLozQZ+bro/QjZrP9c2szscC4CklcBKgHnz5nUQVgv1k8ThnUzuvhD5onz5xItAY1mcWtYyOLJZervpSvR+Rjse4vDolVlRPyPsaIozW6VT6TbWty3OQVXj0Oa6M9mlRKbPTMaKoPbn4EfjH208/2laddr4Nvtt9/dmzCrn8Gc+658kYHx8nHr9P+r1yc9/NeXYKOJqrGnxWjZxU+KGW+dMesHpVieFvtWje2YVaLdNJ2OzhRGbgc0AtVrtvOaIGrwELXzgfIZOOxXcdd3MOtNJoR8B5jbdnwP80uE2Qx2MNTOzAnXyrYEvgGslzZc0BDwIbD9jm+3AY8osAv6MiGMdjjUzswKdc0YfEXVJTwOfkJ0i+WZEHJC0Kl//OrCD7NTKI2SnV64429hCMjEzs5Y0FS/UX6vVYnh4uOwwzMymDUl7I6LWal1/XfDBzKwPudCbmSXOhd7MLHEu9GZmiZuSH8ZKGgV+PM/hs4DjFzCc6cA5949+zNs5d+aqiLi81YopWei7IWm43SfPqXLO/aMf83bO3fOhGzOzxLnQm5klLsVCv7nsAErgnPtHP+btnLuU3DF6MzM7XYozejMza+JCb2aWuKQLvaR1kkLSrLJjKZqklyR9K+lrSe9LuqzsmIoiaamkQ5KOSFpfdjxFkzRX0meSDko6IGlN2TH1iqQBSV9K+rDsWHolb8W6LX8+H5R0e7f7TLbQS5pL1pT8p7Jj6ZFdwI0RsRD4Dni+5HgK0dRw/l5gAfCQpAXlRlW4OvBMRFwPLAKe6oOcG9YAB8sOosc2AR9HxHVkjeS6zj/ZQg+8AjzHWbqjpiQidkZEPb+7m6ybV4ommtVHxBjQaDifrIg4FhH78tsnyJ74s8uNqniS5gD3A1vKjqVXJM0E7gDeAIiIsYj4o9v9JlnoJS0HjkbE/rJjKckTwEdlB1GQdo3o+4Kkq4GbgT3lRtITr5JN1sbPtWFCrgFGgbfyQ1ZbJM3odqed9IydkiR9ClzRYtVGYANwT28jKt7Zco6ID/JtNpK91d/ay9h6qOOG86mRdCnwLrA2Iv4qO54iSVoG/BYReyXdVXY8PVQFbgFWR8QeSZuA9cAL3e50WoqIu1stl3QTMB/YLwmyQxj7JN0WEb/2MMQLrl3ODZIeB5YBSyLdL0h00qw+OZIGyYr81oh4r+x4emAxsFzSfcDFwExJ70TEIyXHVbQRYCQiGu/YtpEV+q4k/4UpST8AtYhI+up3kpYCLwN3RsRo2fEURVKV7MPmJcBRsgb0D6fci1jZjOVt4PeIWFt2PL2Wz+jXRcSysmPpBUmfA09GxCFJLwIzIuLZbvY5bWf0NslrwEXArvydzO6IWFVuSBdenzacXww8Cnwj6at82YaI2FFiTFac1cBWSUPA98CKbneY/IzezKzfJXnWjZmZneJCb2aWOBd6M7PEudCbmSXOhd7MLHEu9GZmiXOhNzNL3P/pH8V+QBkn2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cohen_d = CohenEffectSize(male_sample, female_sample)\n",
    "\n",
    "plot_pdfs(cohen_d, pdf_threshold=True)\n",
    "plot_pdfs(cohen_d, pdf_threshold=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In summary, the best way to report effect size depends on the audience and your goals.   \n",
    "### There is often a tradeoff between summary statistics that have good technical properties and statistics that are meaningful to a general audience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}