{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Association Between Categorical Variables\n",
    "\n",
    "A standard approach for testing the independence/dependence of a pair of categorical variables is to use a $\\chi^2$ (Chi-square) test of independence. \n",
    "\n",
    "The null and alternative hypotheses for the $\\chi^2$ test are as follows:\n",
    "\n",
    "* $H_0$: the two categorical variables are independent\n",
    "* $H_A$: the two categorical variables are dependent\n",
    "\n",
    "## Contingency tables of observed counts\n",
    "\n",
    "We typically depict the relationship between categorical variables using a \"contingency table\".\n",
    "\n",
    "|       | B1                  | B2                  | Total                               |\n",
    "|-------|---------------------|---------------------|-------------------------------------|\n",
    "| A1    | $O_{11}$            | $O_{12}$            | $O_{11}$+$O_{12}$                   |\n",
    "| A2    | $O_{21}$            | $O_{22}$            | $O_{12}$+$O_{22}$                   |\n",
    "| Total | $O_{11}$+$O_{21}$   | $O_{12}$+$O_{22}$   | $O_{11}$+$O_{12}$+$O_{12}$+$O_{22}$ |\n",
    "\n",
    "The rows and columns indicate the different categories for variables A and B respectively, and the cells, $O_{ij}$, give the counts of the number of observations for the corresponding combination of A and B.  For example, the cell $O_{11}$ gives the number of observations that that belong to both the category A1 and B1, while $O_{12}$ gives the number that are both A1 and B2, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sbn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load the bumpus dat set\n",
    "dataurl = \"https://github.com/Bio204-class/bio204-datasets/raw/master/bumpus-data.txt\"\n",
    "bumpus = pd.read_table(dataurl)\n",
    "\n",
    "# for clarity of presentation let's replace the sex designation \"m\",\"f\" with \"male\",\"female\"\n",
    "bumpus.sex.replace(dict(m=\"male\", f=\"female\"), inplace=True)\n",
    "\n",
    "# similarly replace survived variables\"F\",\"T\" with \"False\" and \"True\"\n",
    "bumpus.survived.replace(dict(F=False, T=True), inplace=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Table of observed counts\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>sex</th>\n",
       "      <th>female</th>\n",
       "      <th>male</th>\n",
       "      <th>All</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>survived</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>False</th>\n",
       "      <td>28</td>\n",
       "      <td>36</td>\n",
       "      <td>64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>True</th>\n",
       "      <td>21</td>\n",
       "      <td>51</td>\n",
       "      <td>72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>All</th>\n",
       "      <td>49</td>\n",
       "      <td>87</td>\n",
       "      <td>136</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "sex       female  male  All\n",
       "survived                   \n",
       "False         28    36   64\n",
       "True          21    51   72\n",
       "All           49    87  136"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Construct a contingency table for the sex and survival categorical variables.\n",
    "# The pandas.crosstab function will create the contingency table\n",
    "# if margins = True, will including row, column, and overall sums in table\n",
    "# if margins = False only includes the the joint counts\n",
    "observed = pd.crosstab(bumpus.survived, bumpus.sex, margins=True)\n",
    "\n",
    "print(\"Table of observed counts\")\n",
    "observed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(49, 87)"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(bumpus.sex == \"female\"), np.sum(bumpus.sex == \"male\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Expected counts\n",
    "\n",
    "If the two categorical variables were independent than we would expect that the count in each cell of the contigency table to equal the product of the marginal probalities times the total number of observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Table of Expected Counts\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>sex</th>\n",
       "      <th>female</th>\n",
       "      <th>male</th>\n",
       "      <th>All</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>survived</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>False</th>\n",
       "      <td>23.058824</td>\n",
       "      <td>40.941176</td>\n",
       "      <td>64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>True</th>\n",
       "      <td>25.941176</td>\n",
       "      <td>46.058824</td>\n",
       "      <td>72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>All</th>\n",
       "      <td>49.000000</td>\n",
       "      <td>87.000000</td>\n",
       "      <td>136</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "sex          female       male  All\n",
       "survived                           \n",
       "False     23.058824  40.941176   64\n",
       "True      25.941176  46.058824   72\n",
       "All       49.000000  87.000000  136"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nobs = observed.ix['All','All']\n",
    "\n",
    "prob_female = observed.ix['All','female']/nobs\n",
    "prob_male = observed.ix['All', 'male']/nobs\n",
    "\n",
    "prob_surv = observed.ix[True, 'All']/nobs\n",
    "prob_died = observed.ix[False, 'All']/nobs\n",
    "\n",
    "expected_counts = []\n",
    "for i in (prob_died, prob_surv):\n",
    "    row = []\n",
    "    for j in (prob_female, prob_male):\n",
    "        row.append(i * j * nobs)\n",
    "    expected_counts.append(row + [np.sum(row)])\n",
    "expected_counts.append(np.sum(expected_counts,axis=0).tolist())\n",
    "\n",
    "expected = pd.DataFrame(expected_counts, index=observed.index, columns=observed.columns)\n",
    "\n",
    "print(\"Table of Expected Counts\")\n",
    "expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Difference between Observed and Expected counts\n",
    "Since the  margins of the table (the row and column labeled \"All\") are the same for both the observed and expected counts, let's work with just the inner cells of the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Observed = observed.ix[:-1, :-1]  # take all but the last row and colum \n",
    "Expected = expected.ix[:-1, :-1]  # as above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Observed Counts</h3>\n",
       "\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>sex</th>\n",
       "      <th>female</th>\n",
       "      <th>male</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>survived</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>False</th>\n",
       "      <td>28</td>\n",
       "      <td>36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>True</th>\n",
       "      <td>21</td>\n",
       "      <td>51</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "\n",
       "<h3>Expected Counts</h3>\n",
       "\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>sex</th>\n",
       "      <th>female</th>\n",
       "      <th>male</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>survived</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>False</th>\n",
       "      <td>23.06</td>\n",
       "      <td>40.94</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>True</th>\n",
       "      <td>25.94</td>\n",
       "      <td>46.06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This code cell just shows how you can do some fancy HTML output\n",
    "# from Jupyter code cells.  Here I wanted to output both observed and expected\n",
    "# count tables in the same cell with some nicely formatted headers.\n",
    "# This is purely for presentation, we're not doing an analysis here.\n",
    "\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "def fmt(x):\n",
    "    return \"{:04.2f}\".format(x)\n",
    "\n",
    "Obs_as_html = Observed.to_html(float_format=fmt)\n",
    "Exp_as_html = Expected.to_html(float_format=fmt)\n",
    "\n",
    "HTML(\"<h3>Observed Counts</h3>\\n\\n\" + Obs_as_html + \"\\n\\n\" +\n",
    "     \"<h3>Expected Counts</h3>\\n\\n\" + Exp_as_html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can quantify the difference between the observed counts and expected counts on a per cell basis using something which is essentially a squared deviation (but for counts rather than values of a random variate) normalized by the expectation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>sex</th>\n",
       "      <th>female</th>\n",
       "      <th>male</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>survived</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>False</th>\n",
       "      <td>1.058824</td>\n",
       "      <td>0.596349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>True</th>\n",
       "      <td>0.941176</td>\n",
       "      <td>0.530088</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "sex         female      male\n",
       "survived                    \n",
       "False     1.058824  0.596349\n",
       "True      0.941176  0.530088"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff_squared = (Observed-Expected)**2/Expected\n",
    "diff_squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The $\\chi^2$-statistic\n",
    "\n",
    "We can calculate a new test statistic, we'll call $\\chi^2$ (chi-squared), as the sum of\n",
    "\n",
    "\\begin{align}\n",
    "\\chi^2 &= \\frac{(O_{11} - E_{11})^2}{E_{11}} + \\frac{(O_{12} - E_{12})^2}{E_{12}} + \\cdots + \\frac{(O_{mn} - E_{mn})^2}{E_{mn}} \\\\\n",
    "\\\\\n",
    "&= \\sum_{i=1}^{m}\\sum_{j=1}^{n} \\frac{(O_{ij} - E_{ij})^2}{E_{ij}}\n",
    "\\end{align}\n",
    "\n",
    "where $m$ and $n$ are the number of categories of the two variables under consideration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi-squared 3.12643678161\n"
     ]
    }
   ],
   "source": [
    "# cast to matrix so we can calculate sum across all cells\n",
    "chi2 = np.sum(diff_squared.as_matrix())  \n",
    "\n",
    "print(\"Chi-squared\", chi2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can think of the $\\chi^2$-statistic as the sum of the normalized squared differences between observed and expected counts (where each cell is normalized by the expected counts)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## The $\\chi^2$-distribution\n",
    "\n",
    "Under the null hypothesis that the categorical variables are independent the $\\chi^2$-statistic has a sampling distribution called the $\\chi^2$-distribution.  The $\\chi^2$-distribution takes a single parameter, degree of freedom (df), which for the $\\chi^2$ test of independence is equal to $(m-1)(n-1)$ where $m$ and $n$ are the number of categories in the two variables.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFrCAYAAADFOmBlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8HNW5//HPzPamLrnIklwkWZZ7oWMwBgcHHAIBEwcC\naRcIFy4kkFxCLj8CAWJSSCAhJCGFhITgQBJKTAebasDGXZaLbLlILupte5vfHyvLNrZVdyXt7PN+\nvWSVXc0+x1rpu+fMmXMUTdM0hBBCCJG01KEuQAghhBADI2EuhBBCJDkJcyGEECLJSZgLIYQQSU7C\nXAghhEhyEuZCCCFEkkt4mG/cuJFrrrnmuK+vWLGCK664giVLlvDss88mugwhhBBCt4yJPPgf/vAH\nXnjhBRwOxzFfD4fDPPjgg/z73//GYrHwpS99ifPPP5+srKxEliOEEELoUkJ75kVFRfz6178+7uu7\ndu2iqKgIp9OJyWRi9uzZrFmzJpGlCCGEELqV0DBfsGABBoPhuK+73W5cLlfX5w6Hg46OjkSWIoQQ\nQujWkEyAczqduN3urs89Hg9paWk9fl93K89e/6M3+cq9r8WlPiGEECKZJPSc+WGfDuEJEyawd+9e\n2tvbsVqtrFmzhm984xs9HkdRFBoaTtyDNxoUWt2hk96eDHJzXUldf0+kfclLz20DaV+yS4X29WRQ\nwlxRFACWL1+Oz+dj8eLF3HnnnXz9619H0zQWL15MXl7egB7DYjbgD0bQNK3r8YQQQohUkPAwz8/P\nZ9myZQAsWrSo6+vz5s1j3rx5cXscq8mApkEwHMViOv48vRBCCKFXulk0xmqOBXggGBniSoQQQojB\npZswt3SGuT8kYS6EECK16CbMrabYGQN/IDzElQghhBCDSz9hbukcZpeeuRBCiBSjmzA/POlNzpkL\nIYRINfoJ88PnzCXMhRBCpBjdhLlVwlwIIUQCBINBFi++BIDKygqWLLmM3/3u+H1HeuOdd1Zy7713\nxbM8YJAWjRkMVnOsKXLOXAghRDzFVjGNLUa2evVHLF78JS6//Mo+H+eRRx5izZqPKC4ujXOFOgrz\nw+fM/UGZzS6EEHryzIqdrNlWf9LbDQaFSOTke3ecyClleVw5v/ikt/t8Pn74w7vo6OggP38MANu2\nVfLSSy9iMpnIy8tj7tx5AOzfX8uDD953zOqjCxYs5HOfu/SYY06dOp1zzpnHCy/8u0+19oZuwlyG\n2YUQQsTL88//i/Hji7nuuhuprKxg3bq1lJWV89nPLiI7O6cryAHy88fwq1/9rsdjzp9/AevXr01I\nvboLc5nNLoQQ+nLl/OJue9GJ2GilpmYvZ545F4Dy8ikYjSdfJvzonvnh/UFO1DNPJN2EuawAJ4QQ\nIl7Gjh1PRcUmzj77HHbs2EY4fPJs6W3PPJF0NJu9cwU46ZkLIYQYoEsvvZwDB/Zz003X8dxz/8Ri\nMQ91Sd3STc/cKovGCCGEiBOz2cwPf7j0uK9//evXD+i4M2fOZubM2QM6xonopmduNqkoQEBmswsh\nhEgxuglzRVGwmA1yzlwIIUTK0U2YQ2wSnJwzF0IIkWp0FeZWk0HOmQshhEg5+gpzs1GG2YUQQqQc\nXYW5xRzrmUe1vi3rJ4QQQiQzXYX54VXggtI7F0IIESfx2DXN43Fzxx3f5uabr+fGG79ORcXmuNao\nm+vM4UiY+wKRrkVkhBBCiIGIx65py5Y9xZw5p7F48RL27dvLPff8H3/609/iVqOuEs9uiTXHFwiT\n6bIMcTVCCCHi4d87l7O+/uQ9WYOqEIn27fTqzLypfKF40Ulvj/euaUuWXI3JFFtFLhwOY7HEN6N0\nFeY2a6w53oAsHCOEEKL/4r1rmsPhBKCpqZH777+bW2/9blzr1VWYH90zF0IIoQ9fKF7UbS86WXZN\n27VrJ/fe+3/cfPO3mT59Rlzr1WWYe/0S5kIIIfov3rum7d5dzd13f48f/vBBJkw4+Xau/aWrMJdh\ndiGEEPFw6aWXc//9P+Cmm66jsLBowLumPf74rwkGQzzyyM/QNA2n08XSpT+LU7U6C3MZZhdCCBEP\n8d41benShwZaUrd0dZ253WICZJhdCCFEatFVmMswuxBCiFSkqzCXYXYhhBCpSJdhLsPsQgghUomu\nwtxsUjGoCt5AaKhLEUIIIQaNrsJcURRsFiO+gGy0IoQQInXoKswhNtTu9UvPXAghRHzEY9c0v9/P\nnXfezs03X8+3v30TjY2Nca1Rd2FusxplNrsQQoi4OdGuaTfccFOfjvHii88xceIkHn30cT7zmc/y\n1FN/iWuNulo0BmI982AoSjgSxWjQ3WsVIYRIOQ3PLqPjkzUnvX2vQSUSifbpmK45p5C7eMlJb4/3\nrmlXXvmlzhcFUFd3CJfL1ad6e6LLMIfY5Wku+8CW3xNCCJGa4r1rGsTmdd16641UV+/iF7/o2zB9\nT3QX5kcvHCNhLoQQyS938ZJue9HJsmsawCOP/IZ9+/bw3e9+i3/84/m41au7MJeFY4QQQgxUvHdN\n++tf/0xeXh4XXngRVqsNg+HkLw76Q7dhLgvHCCGE6K9475q2aNEl3H//PSxf/gKapvH97/8gLnUe\nprsw7xpmlzAXQgjRT/HeNS0zM4uHHvrlQMs6Kd1N95ZhdiGEEKlGt2Eu15oLIYRIFfoLcxlmF0II\nkWJ0F+Y26ZkLIYRIMboLc6fNBIDHJ+uzCyGESA26C3NHZ5i7JcyFEEKkCN2FucVkwGRUJcyFEELE\nRTx2TTts7949LFw4j1Aovhmlu+vMITbULmEuhBAiHk60a9rll1/Z5+N4vR5+/euHMZstca5Qx2He\n0Oob6jKEEELEwaoVu6jeVn/S21WDSrSPu6aNL8vjzPkTTnp7vHdNA/jxjx/ghhtu5nvfu71PtfaG\nbsO8pt4t26AKIYTol3jvmvanPz3OWWfNZcKEYkCLe726DXOITYLLcMZ/OEMIIcTgOXP+hG570cmw\na9rrr79CXt4I/vOf52lqauLb376JRx99PG71SpgLIYQQnxLvXdOWLXuu6+PFiy/h4Ycfi1utkODZ\n7Jqm8YMf/IAlS5Zw7bXXUlNTc8ztL774Il/4whdYvHgxTz/9dNwe1yHXmgshhBiASy+9nAMH9nPT\nTdfx3HP/HPCuacdSOifVxU9Ce+ZvvvkmwWCQZcuWsXHjRpYuXcpjjx15NfKTn/yEV155BavVysUX\nX8yiRYtwuVwDflyXXGsuhBBiAOK9a9rRnn32hQEf49MSGuZr165l7tzYOYfp06dTUVFxzO1lZWW0\ntbV1zQA8eibgQBweZu+QMBdCCJECEhrmbrf7mJ620WgkGo2iqrHR/ZKSEi6//HLsdjsLFizA6XTG\n5XFlmF0IIUQqSWiYO51OPB5P1+dHB/n27dt5++23WbFiBXa7ne985zu89tprXHjhhd0eMze352H4\nAl9sk5UISq/uP5wkW719Je1LXnpuG0j7kp3e29eThIb5rFmzWLlyJQsXLmTDhg2UlpZ23eZyubDZ\nbJjNZhRFISsri/b29h6P2ZvLD0KBWI+8odkT98sVEikRl1cMJ9K+5KXntoG0L9mlQvt6ktAwX7Bg\nAR988AFLliwBYOnSpSxfvhyfz8fixYu58sorueqqqzCbzRQWFnLZZZfF5XGd1s4JcF4ZZhdCCKF/\nCQ1zRVG49957j/nauHHjuj5esmRJV9DHk81iwKAquP0S5kIIIfRPl2udKoqCw2bC3XnuXAghhOiv\neO2adtllF3HLLd/kllu+2e9d105GlyvAQezytHZPcKjLEEIIkeTisWva/v21TJxYxoMP/jwBFeo5\nzK1GDjZ6iEY1VDU+168LIYQYfC3738DbWnnS2w+pKpFo33ZNs2eUk5m/4KS3x3vXtG3btlJfX88t\nt3wTq9XKzTd/m8LCoj7V3B39hrndjAZ4/CFc9nguwyeEEELv4r1rWk5OLtde+zXmzTufTZs2cN99\n/4/f//7JuNWr2zBPs8dmtLd7ghLmQgiRxDLzF3Tbi06GXdPKysowGGKRO23aDJqamuJar37D3BEL\n8HZPkPzcIS5GCCFEUon3rml/+tPvSU9P56qrrqWqagd5eSPiWq/uw7zNK5PghBBC9M2ll17O/ff/\ngJtuuo7CwqIB75r25S9/lfvu+3+sWvU+RqOR73//B3GqNEa/YW4/3DOXa82FEEL0Tbx3TXO5XPzk\nJw8PtKyT0uV15gDpziPD7EIIIYSe6TbMjz5nLoQQQuiZfsP88DC7nDMXQgihc7oNc6vZgNmo0iY9\ncyGEEDqn2zBXFIU0h1mG2YUQQuiebsMcYufNO7zBznV1hRBCCH3Sd5jbzYQjGt6A7J4mhBBCv/Qd\n5jKjXQghRAqQMBdCCCGSnK7DPP3wkq4S5kIIIXRM12GeJmEuhBAiBeg7zI/aBlUIIYTQK12HeYbT\nAkCbW8JcCCGEfuk7zF2xMG9xB4a4EiGEECJxdB3mFpMBu8VIa4eEuRBCCP3SdZgDZLostEiYCyGE\n0DHdh3mGy4I3ECYQigx1KUIIIURC6D7MMzsnwclQuxBCCL3SfZh3TYKTMBdCCKFTug/zTJnRLoQQ\nQuf0H+YyzC6EEELn9B/mMswuhBBC53Qf5rJwjBBCCL3TfZi77CYMqiI9cyGEELql+zBXFYUMp1nC\nXAghhG7pPswhNtTe5g4SjWpDXYoQQggRdykR5lkuK1FNk33NhRBC6FJKhHl2uhWApjb/EFcihBBC\nxF9qhHlaLMwb231DXIkQQggRf6kR5tIzF0IIoWMpEeY5nT3zpnaZ0S6EEEJ/UiLMpWcuhBBCz1Ii\nzG0WI3aLkcY2OWcuhBBCf1IizCHWO29q96Npcq25EEIIfUmdME+zEgxFcftCQ12KEEIIEVcpE+Y5\nh8+bt8t5cyGEEPqSMmEuk+CEEELoVeqE+eGFYyTMhRBC6EzKhHlOhoS5EEIIfUqZMM/NsAHQ0CqX\npwkhhNCXlAlzh9WEw2qkrkXCXAghhL6kTJgD5GXaaWz1yb7mQgghdCXFwtxGJKrR3CHnzYUQQuhH\naoV553nzehlqF0IIoSOpFeaZnWEuk+CEEELoiDGRB9c0jXvuuYft27djNpt54IEHKCgo6Lp906ZN\n/PjHPwYgJyeHn/70p5jN5oTV0xXm0jMXQgihIwntmb/55psEg0GWLVvG7bffztKlS4+5/e677+bB\nBx/kqaeeYu7cuRw4cCCR5XQNszdImAshhNCRhPbM165dy9y5cwGYPn06FRUVXbft3r2bjIwMnnji\nCaqqqpg3bx5jx45NZDmkOcxYTAYZZhdCCKErCe2Zu91uXC5X1+dGo5FoNApAS0sLGzZs4JprruGJ\nJ55g1apVfPzxx4ksB0VRyM2wUd/ik61QhRBC6EZCe+ZOpxOPx9P1eTQaRVVjrx8yMjIoLCxk3Lhx\nAMydO5eKigpOO+20bo+Zm+vq9vaeFI5yUdvgxmAxkZ1uG9CxEmGg7RvupH3JS89tA2lfstN7+3qS\n0DCfNWsWK1euZOHChWzYsIHS0tKu2woKCvB6vdTU1FBQUMDatWu54oorejxmQ0PHgGrKcsYm2G2p\namBSUeaAjhVvubmuAbdvOJP2JS89tw2kfckuFdrXk4SG+YIFC/jggw9YsmQJAEuXLmX58uX4fD4W\nL17MAw88wG233QbAzJkzOffccxNZDgCjshwAHGryDLswF0IIIfojoWGuKAr33nvvMV87PKwOcNpp\np/Hss88msoTjjMy2A3CwyTuojyuEEEIkSkotGgMwMqszzJslzIUQQuhDyoW5zWIk02XhUJOn5zsL\nIYQQSSDlwhxivfOm9gCBYGSoSxFCCCEGLCXDfFTnefNDMtQuhBBCB1I0zGMz2g82y1C7EEKI5JeS\nYd41o71ReuZCCCGSX0qG+ejOnvmBRumZCyGESH4pGeYZTjMOq5HaBvdQlyKEEEIMWEqGuaIo5Oc6\nqW/xEQzJjHYhhBDJLSXDHGBMrgMNOCDXmwshhEhyKRvm+blOAPY3SJgLIYRIbr0K8z/84Q80NDQk\nupZBNSY3NglOzpsLIYRIdr0Kc7/fz5e//GWuv/56XnnlFUKhUKLrSrj8nFiYS89cCCFEsutVmN98\n88289tprXH/99Xz88cd8/vOf54c//CFbt25NdH0JY7eayEqzSM9cCCFE0uv1OXOfz0dtbS01NTWo\nqkpaWhr3338/Dz30UCLrS6gxuU5a3UE6vMGhLkUIIYTot17tZ3777bfz8ccfc84553DjjTcyZ84c\nAILBIGeffTa33357QotMlMIRTjbtamJfvZvJY7OGuhwhhBCiX3oV5meccQb33Xcfdru962vBYBCz\n2cxLL72UsOISrTDPBcC+Qx0S5kIIIZJWr4bZn3322WOCPBqNcvnllwOQm5ubmMoGQdHIWJjvresY\n4kqEEEKI/uu2Z37ttdeyevVqAMrKyo58k9HI/PnzE1vZIMhJt2K3GNlXJ5PghBBCJK9uw/zJJ58E\n4P777+euu+4alIIGk6IoFI5wsn1fK/5gGKu5V2cdhBBCiGGl2/RauXIl5513HpMnT+b5558/7vZL\nL700YYUNlsIRLrbta6Wm3k3JmIyhLkcIIYTos27DfPPmzZx33nldQ+2fpocwLxrRed78UIeEuRBC\niKTUbZjfcsstACxdurTra263m4MHD1JSUpLYygbJ4Ulwew7JJDghhBDJqdez2e+8806am5u56KKL\nuOWWW/jFL36R6NoGxchsO1azgd0H24e6FCGEEKJfehXmTz/9NHfccQfLly/n/PPP5z//+Q/vvfde\nomsbFKqiMG5UGgebvHj94aEuRwghhOizXi/nmpGRwTvvvMO8efMwGo0EAoFE1nVC/nBiHnPcqDQA\n9hyS3rkQQojk06swLy4u5oYbbqC2tpYzzjiDW2+9lSlTpiS6tuNc/8Id1HTsj/txD4e5DLULIYRI\nRr26sPpHP/oR69evp6SkBLPZzOc//3nOPffcRNd2HH84wB8r/sYdp9yKzWiN23HHj46FefUBCXMh\nhBDJp1dh7vV62bFjB6tXr0bTNAAqKyu5+eabE1rcp11S9hle3PY6T2/7F1+bfBWKosTluJkuC5ku\nC9XSMxdCCJGEejXMfuutt/Lxxx8TjUYTXU+3lky9hPHpRayt38gHBz6O67HHj0qjzR2kqc0f1+MK\nIYQQidarnnljYyNPPPFEomvpkVE18LXJV/Hg6kd4tupFxqYVMsY1Oi7HnpCfztodDVTVtpKdPjIu\nxxRCCCEGQ6965pMmTWLbtm2JrqVXsqyZXFN+JeFomD9u+Rv+cHx60iUF6QBU7W+Ly/GEEEKIwdKr\nnnlVVRWXXXYZ2dnZWCwWNE1DURTeeuutRNd3QlNzyjm/4BzeqnmXZduf4yvlSwZ8/rxohAuTUWVn\nrYS5EEKI5NKrMH/00UcTXUeffX7CZ9nVtoc1despzZzAmaNPHdDxjAaVcaPSqKptxRcIY7PIDmpC\nCCGSQ6+G2fPz81m3bh3PPPMMWVlZrFmzhvz8/ETX1i2DauDrk6/CZrTxzI7nqek4MOBjloxJR9Ng\n1wHpnQshhEgevQrzn/3sZ7zzzju8/vrrRCIR/vWvf/Hggw8murYeZduy+Er5FwlFw/x+85N4Qt4B\nHa84P3befEdNazzKE0IIIQZFr8L8/fff56c//SkWiwWn08kTTzzBu+++m+jaemVqTjmfHXs+Tf5m\n/rzlaaJa/y+fKy3IQFUUtu2VMBdCCJE8ehXmqnrs3YLB4HFfG0oXjVtAedZEKpu38/LuN/p9HJvF\nyNhRLnYfbMcXkE1XhBBCJIdeJfLChQv51re+RXt7O3/+85+5+uqrWbRoUaJr6zVVUfnq5C+Rbc3i\nlT1vsalhS7+PNakok0hUk6F2IYQQSaNXYT5v3jzOO+88MjIyWLt2Lbfeeivf/OY3E11bnzhMdq6f\nei0m1cRfKv9BnbehX8cpL8oEYOvelniWJ4QQQiRMt2He1NTE1VdfzTXXXMPf//53jEYjH330EU89\n9RTt7cNvHfMxrtFcVXY5/oif329+sl9bphaPScdkVKncI2EuhBAiOXQb5vfddx+zZ8/m/fff55ln\nnuGZZ55h1apVlJWV8aMf/WiwauyTU0fOYt6YszjoqeOpbc92bQzTWyajgeL8dGob3LR7ggmqUggh\nhIifbsN8+/bt3HbbbZhMpq6vmUwmbrvtNiorKxNeXH99oXgRE9LHsq5+E6/tXdHn7y8fGxtq37ZP\neudCCCGGv27D3GKxnPDriqIMq9nsn2ZQDVw39VoyLRn8p/o1NtRv7tP3l3WeN5ehdiGEEMmg20Tu\nbr3zeO0lnigus5NvTvsqZoOZv1Quo6Zjf6+/d+xIFzaLga17mxNYoRBCCBEf3S5AXlVVxfnnn3/c\n1zVNo6Ghf7PFB9MY12i+Wr6Exzc/yW83/Zn/nXML6RZXj99nUFUmFmSyYWcjja0+cjJsg1CtEEII\n0T/dhvlrr702WHUkzPTcKVwyfiEvVr/K45v/wrdm3oDJYOrx+yaNjYV55d4WzpEwF0IIMYx1G+ZD\nvZlKvHym6DwOeupZU7eOp7b9s1dbpk4ZlwXA5l1NnDN99GCUKYQQQvTL8J3FFkeKonB12eWMSytk\nTd16Xt+7ssfvGZllJy/TRsXuZkLhyCBUKYQQQvRPSoQ5gMlg4rqpXyHDks6L1a+ytm5jt/dXFIUZ\nxTkEQhG27ZOlXYUQQgxfKRPmAOkWFzdO+xpWg4UnK5exs3V3t/efWZIDwIaqxsEoTwghhOiXlApz\niM1w/6+p1xBF43eb/swhT/1J71s8Jh2H1ciGnY19XklOCCGEGCwpF+YAk7JKuarsCrxhH49t/CPt\nwY4T3s+gqkydkE1LR4B9de5BrlIIIYTonYSGuaZp/OAHP2DJkiVce+211NTUnPB+d999Nz//+c8T\nWcpxzhg1h4vGLaDJ38JvNj5BIHLiddhnFHcOte+UoXYhhBDDU0LD/M033yQYDLJs2TJuv/12li5d\netx9li1bxo4dOxJZxkldNPYCTh81h30dtfyp4iki0eNnrU8Zl41BVeS8uRBCiGEroWG+du1a5s6d\nC8D06dOpqKg45vb169ezefNmlixZksgyTkpRFK6aeDmTskqpaNrKM1UvHHdu3G41MrEwg711HTS3\n+4ekTiGEEKI7CQ1zt9uNy3Vk+VSj0Ug0GgWgoaGBRx99lLvvvntIJ5cZVAPfmPJl8p2jeH//R7y0\n+43j7jO9c6h9owy1CyGEGIa6XQFuoJxOJx6Pp+vzaDTatdvaq6++SmtrK9dddx0NDQ0EAgHGjx/P\npZde2u0xc3N7Xlu971zcPf9W7n7rZ7yy501GZGZyUen8rlsXnD6Op9+sYv2uJq68cFICHv+IxLRv\n+JD2JS89tw2kfclO7+3rSULDfNasWaxcuZKFCxeyYcMGSktLu2675ppruOaaawB47rnn2L17d49B\nDtDQcOKZ5wOn8t/TvsFDax/jz+ufJepXOW3U7K5bS8aks2VXEzuqG8l0nXhr2IHKzXUlsH1DT9qX\nvPTcNpD2JbtUaF9PEjrMvmDBAsxmM0uWLOHBBx/kzjvvZPny5Tz77LOJfNh+y7Flc/OM/8JmtPG3\nbc+yubGy67bTykegAau31g1dgUIIIcQJKFqSrYYyGK++qtv28Mv1vwc0bpr+X5RkjqfdG+S2X31A\n4Qgnd3/1lIQ8biq8upT2JSc9tw2kfckuFdrXk5RcNKYn49PHct3Ua4lqGr/d9GdqOvaTZjdTPjaT\nPYc6qGv2DnWJQgghRBcJ85OYnD2Rr5R/kUAkwKMb/sBBTx2nlY8A4GMZahdCCDGMSJh3Y/aIGSyZ\neBnukIdfrn+cggIFo0Hl48o6WatdCCHEsCFh3oOz809nccnnaQ928Lstf2RSiZmDTV5qGzw9f7MQ\nQggxCCTMe2FewVlcVnwxrYE2Dma+hWL2sqri4FCXJYQQQgAS5r12QeG5XDJ+IZ5IB9byT/hg+27C\nkehQlyWEEEJImPfFhWPnc9G4BWD2EipaxQfb9gx1SUIIIYSEeV9dNPYCzsw7G9Xq5bn9f6c10DbU\nJQkhhEhxEuZ9pCgKV03+HE53GSFjOw998hua/S1DXZYQQogUJmHeD4qisKhoIaH942kONPPwut/S\n6Gse6rKEEEKkKAnzfjp10ghMjZMw1pfR5G/hF+t+Q723YajLEkIIkYIkzPvJYjZwWvlIOvaM5dT0\nc2kNtPGLdb/loEdWhxNCCDG4JMwHYN6M0QA0VI3mipJLaA928PC637LfLdegCyGEGDwS5gNQOMLF\nxIIMKve0UGKZwZKJX8Ad8vDwut9S3bZ3qMsTQgiRIiTMB+gzpxYA8PonNczNP51rJ30RfyTAr9Y/\nTmXT9iGuTgghRCqQMB+g6cU55GXa+GhLHe2eIKeNms11U64hSmz71E/qNgx1iUIIIXROwnyAVEVh\nwZwCwpEoK9fvB2Ba7mRunv4NTKqJP295mndrVw1xlUIIIfRMwjwOzpo6ErvFyMp1tYTCEQBKMifw\nrVk34DQ5+MeO53l59xuybaoQQoiEkDCPA6vZyDkzRtPuDfHhliOXphW48rlt9o1kWzN5afcbPFv1\nAlFNNmcRQggRXxLmcbJgTgFGg8LLH+4lEj0S2Hn2XG6b/d+MdozkndpV/GHzXwlGgkNYqRBCCL2R\nMI+TTJeFs6eNpr7Vx+rK+mNuy7Ck8+1ZN1KaMYGNjVt4eP3v6Ai6h6hSIYQQeiNhHkcXnV6IQVX4\nz6o9RKPHnh+3m2zcNOMbnDpyFnvba/jZJ49S56k/yZGEEEKI3pMwj6OcdBtnThnJoWYva7YdH9RG\n1ci1k77IZ8eeT6O/mYfWPsbO1t1DUKkQQgg9kTCPs4vPKEJVFJav2kP0BLPXFUVh0fgLubpsMb6I\nn1+tf1yuRRdCCDEgEuZxlpdp54zJI9jf6OGTE/TODztz9Cn897SvY1SNPLHl7yyvfl1mugshhOgX\nCfMEWHTWWAyqwnPvVhOOnDygJ2WXcvvsm8i2ZvLKnjf5xao/EJCZ7kIIIfpIwjwBRmTaOWfGaOpa\nfLy78UCqof5EAAAgAElEQVS39x3tHMl35/wPxRnj+Lh2PT9f+xgt/tZBqlQIIYQeSJgnyCVnjcNi\nNvDi+7vxBcLd3tdldvI/M65j/vizqHUf4Mef/FJ2XRNCCNFrEuYJku4ws/DUQtq9IV5bva/H+xtV\nIzfMuZorSi7BHfTwyLrf8tHBTwahUiGEEMlOwjyBLjy1gDSHmddW19DmDvR4f0VROK/gbG6a8Q1M\nBjN/3foM/9j+HOFo9z17IYQQqU3CPIGsZiOfP2ssgVCEf71b3evvm5RVyv/O+R9GO0by7v4PeXjd\n72gNtCWwUiGEEMlMwjzBzpkxmjG5Tt7fdJBdB3ofyHn2HL4z52bmjJjB7va9PLj6EapadiWwUiGE\nEMlKwjzBDKrKlz9TCsDfXt9x3DKv3bEYzHy1/EtcUXIJnrCXX274PW/te1e2UhVCCHEMCfNBUFqQ\nwRmTR7D3UEePl6p92uHz6LfOjO2N/u+dy/njlqfwhX0JqlYIIUSykTAfJIvPK8ZqNvCvd3bR4e37\nwjDFGeP43im3MiF9LOvrN/Hg6kfY216TgEqFEEIkGwnzQZLhtPD5s8fh8Yf5x4qd/TpGuiWNW2fe\nwIVF82nyt/DQ2sdYWfO+DLsLIUSKkzAfRBfMGUPRSBerKg6xubqpX8cwqAYumbCQm6Z/A5vRyj+r\nXuTxzU/iCXnjXK0QQohkIWE+iAyqytc+W4ZBVXjy1W09rgzXnUnZpXz/1G9TmjGBTY1bWLr6YVk1\nTgghUpSE+SArHOHiotOLaGoP8M93BnapWboljf+ZeR0Xj1tAa6CNX6z7DS9Vv04kGolTtUIIIZKB\nhPkQWHTmWEbnOFi5bj/b9rYM6FiqonLRuAXcOvMG0s1pvLznTX6+7jfUexvjVK0QQojhTsJ8CJiM\nKl+7qAxVUfjDS5V4/KEBH7MkczzfP/XbnDJiJnva97F0zcOsOrBaJscJIUQKkDAfIhNGp3PJWWNp\nbg/w5Kvb4xK6dpONr07+El8r/xIGReWpbf/k95ufxB30xKFiIYQQw5WE+RC6+MwiivPTWbOtnlUV\nh+J23DkjZ/L9U79NScZ4NjZu4f6PH2JD/ea4HV8IIcTwImE+hAyqynWfK8dqNvDUGzs42Bi/HnSW\nNZNbZl7PZcUX44v4+X3FX/lTxVN0BN1xewwhhBDDg3GoC9CzaCBAqL6OUFMTocZGQk2NhJubiHp9\nRP0+oj4fEb+PWwIh/IEwld98GrvFhKIqqBYLqtWKarWhWKwYnA6M6ekY0tIxpqVjzMjAmJODKTML\nxXjiH6OqqFxQeC5Tsifxt63PsrZ+I9tbdvLFiZcxK2/aIP9vCCGESBQJ8zjRolECtTX4tm/Dv7sa\n/769hOrq4CTnwhWjEdVmQ7XasNrt+DuCuP0hooqRdJuJaMBPpKODaCBw0mMAoKqYsrIx5eZhys3F\nlJuHOX80lvwxGLOyURSFkY48bpt9Iytr3uc/1a/yx4q/sS5vGl8svRSX2Zmg/xEhhBCDRcJ8ACJe\nL57NG3GvW4t3ayVR75FV2FSbDVtJKebR+ZiyczDl5GDMzsGUlYXqcKCaTMccKxSO8OOnN1C9v42v\nXVTG3GmjgdiLBC0YJOJxE25rI9LeTritjXBrC6HGBkINDYQa6vFu3QJbj61PtVox54/Bkp+PeUwB\nZxaNZfKMm3lq1/Osr99EVcsuriz9PLPypqMoSsL/v4QQQiSGhHkfaZEInorNtL33Dp7NmyASW6DF\nmJODc9Zs7BPLsBaXYMrJ7VNAmowG7vzKKdz60Nv87fUdFOa5KBrpQlFVFKsV1WrFlJ1z0u+PBgKE\nGhsJ1h0ieGA/wf21BPbXxkYJdh21FrzBwOWj82nJzWS1qY4X6//C6nGTuXLSF8i2ZfX7/0UIIcTQ\nkTDvpYjXQ+uKt2h7ZyXhlthCL5aCApyz5uCcNQfz6NED7t2OzHbwX58r51f/3MQv/7WJ//eVOWQ4\nLb36XtViwZKfjyU/H2bN7vp6NBQiVHeIQM0+/Hv24N9dTWDfXhw1Yc7rvE/Y8D6bs1djL5lI2azz\ncZSUYrDbB9QWIYQQg0fCvAcRt5uWN16jdcWbRH0+VKuV9HPPI/2cc7EWjY37480ozuEL547nX+9U\n8+i/N3PHVTMxGQ39Pp5qMmEZU4BlTAFpZ5wFgBYOEziwH//u3fj3VNNSVcmoQ40o9Zs59MFmUBQs\nYwqwlZR2vRkzMuLVRCGEEHEmYX4SWiRC69sraHrhOaJeLwaXi5zLF5M+bz4Gmy2hj33R6UUcaPTy\n4ZZDPPHyNq77XHlcz2krRiPWwiKshUVw7jxGAu72Jt7+4Blat20mvz7EqAO1BGr20briTQDMo0Zj\nnzQJ+6RybBPLMNgdcatHCCHEwEiYn4Cvqoq6v/2F4P5aVJuNnMVfJGPefFRL74a8B0pRFL762YnU\nt3j5qLKOvEwbl84dn9DHdKZls+izN1J95h6e3vZv6toPUtRuZH54LLkHO/Dt2EHrirdoXfFWrOde\nNBZ7WWe4F5cM2v+NEEKI4ylaki3e3dDQkbBja+EwTS8+T/MrLwGQdvZcci67AmNaWsIe82i5ua5j\n2tfmCfKjv35CQ6ufay+cyLyZ+YNSRyQaYUXNe7y0+w1C0RBFaQVcMe5iRjaF8G6txLdtK77qXV2T\n/xSjEWtxCY7JU3BMmYp5TMEJRxI+3T690XP79Nw2kPYlu1RoX08kzDsFDx3i4O9/S2DvHky5uYz8\nxvXYiksS8lgnc6InZF2zlwf+uhaPP8TNl01lZmnuoNXT7G/h+Z0vs7Z+IwCnjpzF5yd8lgxLOtFA\nAF/Vdrxbt+LdWkmgZl/X9fCG9IyuYLeXT8bgdJ60fXqi5/bpuW0g7Ut2qdC+niQ0zDVN45577mH7\n9u2YzWYeeOABCgoKum5fvnw5Tz75JEajkdLSUu65554ej5mIH5incgsHf/MoUZ+PtDPPJu+qq1Gt\niT0vfiIne0LuPtjOj/++Dk2D266czsTCzEGta2frbv654wVq3AcwG8wsLJrP/IK5mAxHrpUPd7Tj\nrdyCp2Iz3orNRDo626EoWMeNxzFlKvlnn4YvYwSKqs9VhPX8B0XPbQNpX7JLhfb1JKFh/sYbb7Bi\nxQqWLl3Kxo0b+d3vfsdjjz0GQCAQ4HOf+xzLly/HbDZz++23s2jRIs4777xujxnvH1jrO29T/9ST\nKKrKiGu/RtqZZ8X1+H3R3RNyc3UTv/znJowGldu+OJ2SMYM7uzyqRfnwwBperH4Vd8hDjjWLL5Qs\nYlrO5OOG1LVolMC+fXgqNuHdUoFv106IRgFQHQ4ck6dgnzwVx9Rpg3YKYzDo+Q+KntsG0r5klwrt\n60lCJ8CtXbuWuXPnAjB9+nQqKiq6bjObzSxbtgyz2QxAOBzGMoiTqDRNo+m5f9H88nJUp5P8m27B\nVlI6aI/fV1PHZ3PjpVP4zfMV/PyZjXznizOYkJ8+aI+vKipn5Z/GzLxpvLLnTd6u/YDHNz/JhPSx\nXFp8MePTi7ruq6gq1rFjsY4dS/aiS4h4PbEV8nZup+mTdXSs/piO1R8f6bVPm45j2nQsBYWyEp0Q\nQvRDQsPc7Xbjch15RWE0GolGo6iqiqIoZGXFVhz761//is/n48wzz0xkOV00TaPp+X/T/PJyTCNG\nkH/r7Zjz8gblsQdiVmkuN1wymd++sIWfP7OB7yyZybhRg9uztZtsXF7yOc4afRrP73qZzY2VPLT2\n18zIncIl4xcywnH8/6PB7sA1+xRyF84nrb6d4IH9eDZvwrNpI76dVfird9H0/L8xZmbimBoLdvuk\ncpkhL4QQvZTQYfYHH3yQGTNmsHDhQgDmzZvH22+/3XW7pmn85Cc/Ye/evTz88MNdvfRE27fsGWqe\n/gfWUSOZ8sB9WLKTaxnTd9fX8tBTa7FZTdz/zTMpHuQh96Nta9jJ3zY+x46malRF5fzxZ7F48sVk\n2Ho3ahB2u2lZt4GWT9bSsm4d4Y7YFq2q2Uz61MlkzplD5pxZWJPgxZYQQgyVhIb566+/zsqVK1m6\ndCkbNmzgscce4/HHH++6/a677sJqtXLXXXf1+pgDPS/S/NorND77D0w5uYz53+9hysoe0PHiqS/n\nfVZVHOSPy7ditRi59YpplBYMXaBrmsbGxi28uOsV6rwNmFUT5xeew/mF52IzWrvu11P7tEgEf3U1\n7k0b8GzaSHB/bddt5vwxOKZNxzltBtYJE4blJDo9n7fTc9tA2pfsUqF9PRm02ewAS5cuZcuWLfh8\nPiZPnswVV1zB7NmxdcQVReHaa6/lggsu6PaYA/mBuTdu4MCjj2DMyKTgjjsx5QzeZV690dcn5Idb\nDvGnl7ZiUBX++7IpTJtw8o1YBkMkGmHVwTW8vPsN2oMdOIx2Lig6l3Pyz8RqtPS5faGmRjybNuLZ\ntBHv1kq0cBjonEQ3ZRqO6dNxTJ6KwTE8VqPT8x8UPbcNpH3JLhXa15OUuc48ePAA+x74IVokQsH3\n/i8h66oPVH+ekBt3NvLY8xVEoxpfv3gSZ0wemaDqeqZpGtGohsfv472aj3i39iP8YT8Oo51zxpzB\nwknn4m4LxO6sKBiNKgaDisGooPbQ044GAni3VuLZHAv3w5vdoKrYikti59qnT8c8auAb3vSXnv+g\n6LltIO1LdqnQvp6kRJhHvB72PXAfobpDjLzum6SddnoCKhu4/j4hd9S08sg/N+ELhLnqghIumFPQ\n8zd1Q9M0goEIXncAjzuIpyOAxx3A7wsT8IcI+MOdb7GPw6EI4XCUSDhKf59NigKGznA3GlVMZgNm\ni7Hz7eiPjVgsBoy+djiwh2j1NrR9OzGF/ahosa1oO2fH2yaWoZoGZx4G6PsPip7bBtK+ZJcK7euJ\n7tdm1zSNQ3/6A6G6Q2QuvGjYBvlAlBZk8L2rZ/Hzf2zg729W0dTuZ/G8YlT15D3UaDSKuz1Aa7OP\n9hYfrS1e2lt8tLX4cHcECIeiPT6uyWzAYjVis5sxGFWMJhWj0RD72KiidD5+JBrmkLeBem8DES2C\nSTUxyj6CTHMmWhQinS8EIpFo14sCvz9Me5ufaKS7VwcOUGfD2NipGpMSxhT0Yt7uxVy5Dkt0Fa6c\nNDLGjiarbALpo3NwOC0YjMPvfLsQQgyE7sO87b138GxYj61sEjlfuGKoy0mYgjwnd14zm0ee3chr\nq2s42OTlhksmY7MYCfhDNNV7aKx309T51tzgIXKCoLRYjWRk2XE4LThcZuxOCw6nGYfTgtVuwmI1\ndr31NDR+3LHTFJ5Z/zJv17zP7mgIl8nJ/MK5zM8/45iJckeLhKMEAmGCXW8RAv4QPm/sze8Ndn3s\n8wbxeWy0+o66XC8I7AB27AJ2AWC1qDjT7TjTLaSl23BlWElLt3a+t2Ey93/LWSGEGAq6HmYPNTSw\n5567UFSVonvvH1Yz108kHkNFHl+Q3/1zM4f2tzHCYiLXYsTTHjjmPgaDQmaOg6wcB+mZNtIybWRk\n2UjPtGGxmk5y5IE73L6OoJuVNe/zTu0q/BE/NqOVc/PPZF7B2bjMzgE/TjSq4fPGTg+01dTTsnMP\nbbUNuNt9BFQbfqOdgMlBVDnxa1mb3dQV7K4MK2kZVjIy7aRn2bA7zCc9J6/noT49tw2kfckuFdrX\nE932zDVNo+5vf0ELBBjxjeuGfZD3l6ZpNDd4qNndwoGaVg7VtmHzhxmHCoEI7YEIuSNdFBZlkJ3n\nJCfPSUa2rc+96nhymZ1cMmEhC4rO5d3aD1lR8x6v7l3BWzXvcdboU7mg8Fwyrf2/1E5VldjIgtNC\n3qg0OLUYgKjfh6eyEs+mjbg3r8Df4cdvcuKzpBMeOZZQ1mh85jTcPo3GQ27qDxz/x8FkNpCeaSMj\nKxbuGVn2zhdC9n7XK4QQA6XbMO9Y/THeLRXYJ0/BdfrgrCw3WIKBMDW7m9lX3UxNdTMed7DrNle6\nlaIJ2Ywck05Nu49/frgPpb6dK8pzOb08b1gtl2oz2rhw7HzOKzibVQfX8Obed3i79gPe3f8hc0bM\nYH7BXApc8dv2VbXacM2ajWvW7CPrx2/agHvTRgI7VnTdz5Q3AvvUaVAylXB2Ph0dQVpbfLQ1e2lt\n9tHS5KWxzn3c8R1OM66M2ChHZraDrFw7mdkOnGmWYfX/LoTQH10Os0cDAXb/3x1EPR6KfvgA5tzk\nWD2su6EinzfInqomqnc0ULunpWtimNVmomBcJgXjs8gvysTpOnYJ1O37WvjtC1to8wSZVZrL1y8q\nw57AofTu9DQUFolGWFO3njf2vs0hbz0ApRkTmF84l8nZZahK4kYTwm2teDZvjl36VlGBFvADoFis\nOCZPjq0fP3UaxvQMNE3rnDzopa3ZR2uzl9YWH+42P63N3uNm9JvMBjKz7WTmOMjMsZOV4yAz244r\n3Zo0IZ8Kw5jSvuSVCu3riS7DvOnF52l68XmyFn2OnEsvH4Sq4uPTT8hQMEz19kZ2bKlj/96WrpDI\nznMwrjSXoglZ5I509RgIbe4Av3txC9v2tZKbYeW/L51K0cienxzx1ttfuKgWZWvzDlbse49tLVUA\n5NlzOG/MXE4bNRuLIbGXm2nhML6qHbg3xlaiC9XXdd1mGTsOx5SpOKZMxTpuPIrhyGS53FwXhw62\n0dbio6XJQ0ujl+ZGDy1NXlqbvESjx/6qGU0qmdlHBXzn++EY8qnwx1Lal7xSoX090V2Yh1tb2f39\n/0W12Rj3wI9RrSeeJT0c5ea6qK9vp3ZPCzsq6qje0dB1iVjeaBfjJ+YyvjSX9My+77UeiUZ5/r3d\nvPThXowGlSvOHc8FpxSgDmJo9OcXbr/7ICtq3uOTQ+sJaxHsRhtn55/O2aNPJ9s2OPu6Bw8dip1n\n37QBX9UOiEQAUO127JPKcUyZin3yVEZPLDpp+yKRKO2tPloavbQ0emjufN/S7D3u8juT2UBWroPs\nXAfZuU6y8mIfJ3JyYk9S4Y+ltC95pUL7eqK7MK9f9hStb75B3jVfJePceYNTVBz4fSFqdjWz+v3d\ntLfGhnjTMqyUThlJ6eQR/QrwE9m0q4k/vlRJhzfEpKJMvnHxJLLSBucFz0B+4doCHby3/0Pe2/8h\n7pAHBYUpOWWck38mZVklJx2C16JhImEv0bCXaMRHNOInGgmgRYPHvI9GgmjRzvdaCLQoWjSCpoXR\ntAhoETQtEvtaNATRKJ/+1VGUzn9iH4CioCgGFMUAne8V9fDnauw9KpGoSjikEAgo+P0KPm8Uj1sj\nHFaJRAyEIwYiERWT2YrD5cSZnkF6djpZOdmkZ7swGhN/KV0q/LGU9iWvVGhfT3QV5uHWVnbf+V0M\n6emMu/9BFOPwn9/XcKiDzZ/UsnNrPZGIhsGoUjIpj7LpoxiZn5aQ4dZ2T5A/v7KNDTsbsVuMfPnC\nUk4vT/wysPH4hQtGQqyr28Dq/ato9R7EqSqMtLiYlDaa0ZY0lKifaNjbFeBaNNjzQY8TC2FUY1cY\nxz4/6mNFBQ20UIiI10PU6yXq89F1LkRRUO222JvVDAb1uBcFhz+G/v8KRqMKkYgJDTOqwYrRbMds\ntWOyODAYragGO6rRjsHkwGC0oxodqEY7qtq3Xn4q/LGU9iWvVGhfT4Z/2vVB86svo4VCZF/0uWEd\n5JqmUbunhfUf7WP/3lYA0jNtnDZ3HGPGZyZ8ODXNYeZ/Lp/Ke5sO8vSbVTz+YiXrtjdw9YJS0p1D\nv4e4pmlEw17CwRbCgdbY+2Ar4UDs/ehgO5cao5B2+HKwMPj2EfB1fr9iwGh0YLRkYTDaugJNNdpQ\nDTZUgxlVtaAYzKiqGcVgOepzC4rav55uVpqZmg/X4anYjLdiM/6DO7tuM+XmYZ8yFcfkKdjLJh1z\n+kfTorERgmgILRI88nE0iBYNdY4axN6Cfi8+j5uAz0Mo6EWL+IEARmMAg8ED4UaCbggeP9n+GIpq\njoV8Z7gbjJ1hb3JiMLkwmJwYjLGP1QTPURBCDJxueuYRr4fq796GwW5n3NKfDsswj0Y1qrc3sP6j\nfV2XNo0Zm8n0UwsoGJdJXl7aoL+6rG/x8seXtlJV24bdYuTK+cXMnTYqISMCn371rGlRwoEWQv5G\nQv4GQv5GwoFGQv7Gk/aoDSYXBnN6Z+C4MJpchFULW9tr+ahhKzXeFkLASMcIzhg1h9NGzo7LQjT9\naV+oqQnPlliwe7dWxnruAAYDtvETsJdPxj6pHOvYcQN6vkajGu2tPprq3DQ1ttHe3EpHaxsBnwez\nKYTZEsJsCmGxhHA6o9jtESyWMAZDAEXzAd0v3auoZszWdFDtsf93ozMW9odDv/Nn0tfe/nCSCj07\naV/ySqlh9sP7lOdcvpisz148yFV1T9M09lQ18vG7u2lp9AIwoSyXmacXknvUrPKhekJGNY131u/n\n2bd34Q9GKCvM4CsLyxiRFb+FUCJhL3ZzGw0HdxP0HSLkqycUaAItcuwdFQMmSxZGSxZGcyZGS0bX\ne4M5o9vAiGpRtrfs5MMDa9jYUEFYi6AqKtNyyjlj1ClMyirF0M9ed2909/PTwmF81bu6gt2/Z3fX\nkLxqtWKbWIZ90mTs5eVx2/ktFAzT1OChqd4TW8a3wU1TvYdQ8Oj/cw2HUyF3hJGcXAMZmRpOZwSL\nJYAW8RIJdRAJudEibsIhT7ePpxpsGMxpGExpGM3pne+P/tyFog6/F9mQGmEg7UteKRPmWiTC7jv/\nl4i7g/E//cWw2d/68HD66nd3U3+wA0WB0ikjmXVGIRknCMqhfkI2t/v52+s72LCzEaNBZeFphVx8\nRhEWU98CMBL2EvTUEvQeJOg7RNB7iEio7Zj7KKoZkzUHkzUHoyUHkzW38+NMlDhcT+4OeVhzaD0f\nHlzDfvdBANLNaZw+ag5njDqFXHv8VwTsy88v4vHg3b4N79YteCu3EKo7cvmbIT0De3k5jkmTsU0q\nx5QZv1n7mqbR0eanqd5N4+GQr3d3Tbo8TFEgI9tOdp6T7FwH40tyMZrBYgkRDbuJhNxEwrGgj4Q6\niATbCIfaiQTbu52noBodXUFvMKdhNKVhMKdjNKdjNGegGh1DclneUP/uJZq0L7mlTJi7N6znwKOP\nkD5vPiO+fO0QVHW8pgY3H7y5s+uc+ISyXE6ZO5bM7JO/0BgOT0hN0/hkewNPv7mDVneQrDQLV55X\nzCllJ149TtOihPwNBDw1BD21BDy1hAPNx9xHNTox20eSkV1IWMvCZB+J0Zw5KH+0NU1jX0ctqw6u\n4ZNDG/BHYqE1Pr2IU0bMZFbedJzm+Lz4G8jPL9TUiHdrZeytspJIR3vXbebRo2O99knl2EonYrDH\nf+nY43rx9W6aGj7diwerzdgZ8E6y8xxk5znJzLF3zajXNA0tEiAcaiMSbO8M+DYioXbCwfau0D9u\nRKaTohgxmDNi4W7JwGjOOPJ5AsN+OPzuJZK0L7mlTJjv//Uv8axfR+Hd92ItLBqCqo7w+0KseW83\nW9YfQNOgYHwWp50z7pjh9JMZTk9IfzDM8lV7eX3NPsIRjbLCDL50QSljcm0EPfvxu/cQcO8l4Nl/\nTE9MMViw2MdgcYzBbB+N2T4Kgyl2znqo2xeMBNnQUMGHBz+hqmUXGhqqolKeNZFTRs5kWk455gFM\n9opX+zRNI7i/Fm9lZaznvmM7WqBzsxxFwVI0FnvpRGwTy7CVlCYk3A/XcbgX73OH2LenucdefE6e\nk6xcBzl5TuzOE29Kc3iCY1fAh9piExyDbUQCrYSDrUQjvhPW1BX2nUFvNKd3hn1GZ9jb+xX2Q/3c\nTDRpX3JLiTAPt7dT/d1vYxmdT9EPfjhEVcUmIVVuOMDqd3cT8IdJz7Jx1vnFFE3o/XDucHxCHmpy\n8+oHawl79zIuq41x2e0Y1SO9KqMlB4tjTOdbAUZrTlLsKtYaaOOTug2sObSeWvcBACwGM9Nzp3Dq\niFmUZk7o8/n1RLWv63x75RZ827fhq97VtXANioKlsAj7xDJspROxlZZisMf/NNPRbQsGwjQ3Hu7B\ne7rpxZu6eu/ZubH3WTmOXu0nH40ECAfbCAdbYr35zpAPB9uI9BD2RktmZ+BnYjJnYrRkYjRnYrCc\nfM7FcHpuJoK0L7mlRJi3vP4qDc8sI3fJ1WResGBIamqqd/P2K9upP9iB2WJg9pljmTonH4Ohb+d+\nh8sTMhJy42uvwte+k0DHbqKRIz2xBreNPS2ZpGVO4LSZs0hzpvf6uMOlfZ920FPHmkPr+aRuPU3+\nFiC2s9vM3KnMzJtGcca4Xq0LP1jtiwYC+Kt34d2+Dd+O7fird6GFw7EbFQVLQWEs3CeWxS3ce2rb\n4V58Y527c7i++3PxOXnOWMjnxVa5O1kv/mRiYR8L+GPDPnYpoxYJnPD7DCZXrBdvyep8Hwv6EaML\naGnTht0yuvEyXH/34iUV2teTpA/zvfffS6BmHxN+9jAG1+CuNx4JR1m7ai/rP9pHNKpRUp7HmecX\nY3f0b6h2qJ6QmqYR9B7A116Fv62KoO9g120GczpW5zisrnGYHUWs3uHh+feqaWoPYLMYWTBnDJ85\npaBXm7cM9184TdOobtvLmrr1rKvfiCcUu/LAZXIyPXcyM/OmUZIx/qQ99iG7GiEYPBLu27edMNxt\nE8uwl5ZiLS7B6Err82P0t23BQJjmBg9NDbEJd8296cV39uR724s/kWjYRyjYElubIHB4nYLmrvA/\n0UI9imrqvHIi85igP/z5cJ2J3xvD/XdvoFKhfT1J6jAPNTWx+47bsZdPZsxt3x3UOg7tb2Ply9tp\nbfLiTLNwzmdKKSoe2AzpwXxCRqMh/O3V+Nq24WvfSTR8+LIjFYuzEFt6Cba0EoyW7ON6K6FwhLfW\n7uflj/bi9oWwWQxcMLuABacU4LSdPNST6RcuEo1Q1VrN+obNbKyvoCMUWxfAYbIzPWcyM/KmUZZZ\nfEywD5f2dRvugGnESGwlJdiKS7AVl2IaMaLHHmk826ZpGu2t/mMm2p2oF6+qSuxcfK7jSC8+z4nd\n0XCytu4AACAASURBVLde/PGPH4kN4QdaOhcmasFAB56ORsKB5m7WOEg7KuAzjgr6zH6fqx8sw+W5\nmSip0L6eJHWYt7z5Og3L/k7eNV8h49zzBuXxo9Eon3ywl3Wr9qJpMGXWaE47dzxmy8BftSf6CRkL\n8J14WyrxtVd1/dFSjQ5saSXY0kuwusajGnq3Cpw/GGbl+v28+vE+OrwhrGYD588ew4WnFp4w1JP1\nFy6qRdnZupv19ZvZ0LCZ9mCsDTajjcnZE5mWU0559kQKR+UNy/YdDnffzip8VTvw79pJ1H8kOA0u\nF9biznAvKcVaWHTcIjaD8bM73ItvrD92qP7wZkOHWe0msjsn2WXlOcnJc5CZ3f9ePBxpn6ZpRCO+\nY4L+6I8jofYTfr+imo/tyR/Tq0/vXId/6CTr715vpUL7epLUYV7zk6X4qnYw/me/wJiekfDHbmvx\n8dZ/tlJ3oB1nmoXzF01idGH8HjcRT8hoJIivvQpvayX+9p2xTUIAozkTe8YkbBmTMNsHtkhJIBjh\n7Q37eeXjfbR7glhMBs6eNooFc8aQl3lkprUefuGiWpTqtr1sqN/MhoYKWgKxSw9VRWVyXgkT0yYy\nNaecHFvWEFd6clo0SnB/Lb6qHZ0BX0W45cjlhIrJhHXceGwlpdiKS7COG8/IcaOG7BTQMb34+tiQ\n/Ul78Z3n4LM7Q97Wy158b5+bWjTcNTHv00EfDrZ0/X4dS+m8lj4T06eD3pKJakj8Rkd6+N3rTiq0\nrydJG+bhjnaqb7sVW3EJBXd8P6GPqWkaOyrqeO+NKkLBCMXleZzzmZK4r6Eev0ubovg7qvE0b8LX\nug1Niw2xGi3Z2DMmYc8ox2TreWi1rwKhCO9sOMBrq/fR0hFAAWaV5vKZUwsozk8fkuVqE0nTNGrd\nB9ncuIXNjVvZ11Hbddtox0j+f3tnHiRHdef5b96ZlVl3Vd+X1JIQSCCb0xdgwFwejxeM2BmwHTiW\nGY89diwxMevAMMQwNmEbO8YzG15DGMLr8RpHDOEDBnaYsWysNbeNMJeEaCS11K1W33Xfee8fmV2H\nuvqSqlRdrfeJyHgvX2Zlv+yqzO/7veP3Oz9yHs6PnItBX/+qJtC1Ej0eR/FIRdy1yROVoDEAxJ4e\ncINDkDYPQ9w8DKG3r6UukzXVqLHeF/L1rPhI1Wz68BJWfCOePWfJXX6RwC+IvmnUd5hPM9KSVj3D\nNSbY0tkgdhv9/laibcU88+rvMfPoD5ruvtXQTbzw68MY2T8DjmdwxfXbsG1HZ1P+1un8IG3bhl6c\nQT7xNvLJA+UxcFYIwRPc4Qi4WN/xS6MxTAuvvTeHPa9OYHzGuZ/NPT7svnobtnQrYNc4y79dYGQT\nvzu0DwdiBzGSPALDchpRMufB9uBWnBc+B+eGtsEvrH0C2pnGLOTdrvkjKB07CnXsKMx8oXyc4nmI\ng0MQN2+GuMkReC7U2t6Iais+NpdDYs7pss+mF1vx/pCEUERGMCIjFPFgeGsHTFhrXoGyFixLX7L7\n3tBS9R3pUMziyXhV+dX6wz8bxG6j399KtK2Yz/yfHyHzwvMYuO9+iEObmvK3Mqki9jz5DmKzOUS7\nFFx30w74Ao2JK16PU/lBmnrWEfDE29BL8wCclr4nuANy6ALwnt6WTcyxbRuHJlL49b4JvHk4BhtO\nxLbLL+jGlbt6EGni/7IVVH9/qqlhJHEYB2IHcTBxCCm14s62V+nGeaFzcF74HGz2D4Jtg1nSkbCM\nyf2HUTo2itJRZ1NP1FrvbDAIcfMwxE2bIQ5tgjA4BEZq/Xe8yIqfyyMRWzyjvkbkwx6Eoo4V7w9J\nTRV5wOlNM/XskmK/1Lp6hlXqjNG7Ql/lLe9sELuNfn8r0bZifuyer8DM5zH8P78Pim78gzZxLIHf\nPHUQasnA9gu6cPl1W8suK5vFqsft3G70XOx1FNPvAbABioHk2wo5dAEk39ZTDuPZLGYTBfx+ZB7P\nvjqOfMkABWDn5jA++v4eXDAcBtOE7/BMs9T3Z9s2pvOzeDdxCAfj7+FI6igM1woTGB7bgsM4J7gV\n24LD6JY712WXfL17s0ollMbHHAveFXgzXeuDn+vqciz4wU0QhoYgDgzWhH9tFbZtI59VkYgVkIzl\nUczpmDqRQjKeh6bWEfmghGBERjDicS16DwIhT9NFfgHLKC3qutfd/GqW2nkDndBNuSL2fGDdvSNO\nByLmbSrmemwex776FcjvvxC9X/rvDb2+bdt44/fH8YfnjoFmKFx+7Vac976ehv6NpVjpB2loGeQT\nbyIXf8N9gAFO6oISvhBycAdotvVW0HJEo15MTqWwb2QOv3tzEqOTzszgoFfAB3Z04kM7u9EbWR9B\nck6F1b5QVFPD4eSoI+6J9zBXiJWPKZyMrcFhnBMcxrbgFnRIS3vUO5Os5t5s24aRSDiW+/gYSmPH\noI6PVUK/AgBFge/qdoR9cBPEwSEIAwOghdWtoGgW1bPZ81kVyXgBifm8k8bySMYWizxFAf6QB6GI\nB8FwReADIQ84/swJZb2ldtX5+kvtKCfQDb94Qh7HB9f9u+RkiJi3qZinX3wesz/+EaK3fRrBaxrn\n9c00LTz/q0MY2T8D2Svghk/tQEf3mRvfrPeDtG0bpcwR5OJ/RDF9GIANiubgCZ4PJXIheKk5sceb\nwcn3d3w2i+fenMLvD86iqDrjy0NdXnxoZxcuO68TXs+p+0lvBaf6QokXkziUGsWh5BEcSo7WdMn7\neR+2Bbe44j6McItmyZ/qvdmWBX1+DqXxMahjjsCXxsdhq1Xj2BQFvqcX4sAghP4BCAMDEPr6wShn\nJg49sDoPd/mchmQs74p7oZw/WeQBQPEJZWEPhj0IhCUEQh7IXuGMPq8LfvAVj4rY7OQisTf1+vdM\nM+ISy+wCYHhfy5fanQwR8zYV8+n//Siyr7yMwa99A0Jvb0Ouq5Z07HnyHUyOpxDt8uLG3TshK2fW\nWqj+QVqminziLWTnXy1HIeOlbiiRC+EJ7lz1WvD1xFIPnKabePNIDC8fmMGBowlYtg2GpnD+5jAu\n2d6B922NQGrAOv5m06gZ0XPFGA4lK+Keq4ojHhQCGA4MYdi/CcOBoTPWLd9QpzGWBX1u1hH2sTGo\n42MoHR+vBJNxYUMhR9z7ByD090PoHwQXiTRlWO2UGyu2jUJOQyKWRypeQCpRQNJN89nFFjHHM06X\nfdi14sOO2PuDEtg1hhpeC0vdn2XpMLVUWeD1k7ry60e3o1wHOs5yO8aNaLewOWJ/ZoeKiJi3qZiP\n/f29MBIJDH/v4YY82JlUEf/x8/1IxgvYtDWCa/703DPaTbZANOrF1MQYsrF9yMffdLrHKAZy8Hx4\noxeD95yZ7v5msZoHLp1T8YeDs3jpwAwm5pylPCxDYecmR9h3bYnAI65PYW+KnwDbwnR+1hX3UYym\nj5XdzAKO45rN/kEM+4cwHNiEQW8fOKaxSyaB5r8sbcuCPjsDdWICpYnjUCeOQ52YgJlO1ZxHiyL4\nvn4I/QMQXZHne/tA86fXi9OM+9NUA+lk0RF3V+BT8QJSySJMw1p0vtcvIhCSEFgQ+pAEX0CC4hNB\n06dnzZ/K/dm27UzKq7Hm0zA1N13CgU6lC79a4CthbZ3ldo0VeyLmbSjmsydiOPLlLzRsffn8TBbP\n/OxtFAs6dl3Shw9cNXzaD85asW0bpewo1PTryMRGADgBIZTIxVDCF4Lh2nccuZq1PnDT8Tz2jcxh\n38gcJufdpXausF+8PYpdWyKQG7zW/3Q4Ey8U27YxW5jHaPoYRlNjGE2PIVaMl4+zFIMBXx82+4cw\n5BvAkK8fAcF/2l27rXpZGuk01BMTrrg7mzYzA1hVYkhR4Du7wPf2gu/phdDbB76nF3xnJyhmdY3y\nM3l/C0FpKlZ8sSz2hdxia55mKPgDEnxBCYGgm4Yk+INOt/1q3lfNuD/bMt249Sl3zH4h8E1qFWLv\nd2PUBythbBfi13PeNYs9EfM2FPPjf3gTE998AIGrP4aO2z9zWteankjhP36xH5pq4vJrt2LnRY3p\nsl8ttm2ikHwHmdmXoZfmAAC83Adv9DJ4AtvX3bjU6XI6D9yCsL82MocTrrDTFIUtfX7s2hLGruEI\nusOt9Y/dqhdKWs3iaHqsLPAnclOw7IrY+XkvBl1hH/T1Y9DXB2mNE5zW08vS0jRoU5NVAj8B9cRE\n7UQ7ABTLguvscsS9txdCTy/43r66XfXr5f7UkoFUooB0wrHgM8ki0skiUokiNNVYdP7JQu8PSfAH\nFwt9K+7PtgwYeqYco746hK2hpZYcrwdo17L3g+H8YBfyvB8s5wfD+xYNM66X769ZbEgxP/zzpzD3\n05+g83P/Df6PXHHK1zl+NIE9TxyAZdm4+hPbsfW85jiCqYdlasjFX0d27vdu65WCJ7gDA9uuRkFr\nvlvaVtGoB246nsdr783j7SMxHJ3KlBflRAMidg1HsGtLBNv6A+BOw1f3qbBeXiglQ8Xx7ATGMhMY\nzzhp9aQ6ChQ6PVEM+vox5BvAgK8XPXI3+GW659fLvS2FbdswUilokyegTp6ANjXlpNNTi8biKZ4H\n393jinsv+J4edJ23FVlaasp4fCOwbbss9JlkESlX5DPLCD3DUPAFHHHv7PGB4xl4AyJ8AQlevwiu\niWP0q6HiGtcNY6tVRN/U0suIPUAxYlnYWd4Pf7ADJU2sagSs3bpfz2xIMT/w3f+F9PO/w8Dffw3i\nwOApXWN0ZB7PPn0QFE3h+pt2nHa0s9Vi6nlk519FLrYPllkCRXOQw++Hr+MDYPnAun9hni7NuL9M\nXsP+o3G8dSSGA8cSKLmOQHiOxjn9QZw3FMR5QyH0ReWmW+3r+ftLqemKuKePYzw7AdWsdOnSFI0u\nTwf6vD3oV3rQ7+1Fn7enbMGv53tbDtuyYMTjUKcmHaGfmoQ2OemIvFErgBTLguvoBN/d7XTbd3eD\n7+oG19W9LpzfLIVt2ygVdaRdgS9viSLSyULd2fYA4JF5+AKiI/B+CT5X6H0BER5ldd33zcS2TJh6\nxrHutbRr1Wfcrn1nf6kId84kPW+VRe8Dw7vizzllNCO21UqglWg7MX/trv8BdeI4tj70yCn5hn5v\n/wz+33+MgOUY3HjLTvQOBptQy1p0NYHs3O+dSW22AZr1wBu5BEr0EjDsxgpEshzNvj/DtHBoIoW3\njsRx4Fgc0/HKRDGfh8N5QyGcOxTEeYMhhP2Nd1zSTt+fZVuYLcxjLH0cE7lJTGSncCI3Bc2sfTlG\nxBD6vL3Y3rUJITqCPqUXfmHlF8t6xzZN6PPzjvU+Mw06GUNm7Dj0memaiHILMP7AIpHnu7vBBkPr\n1poHKkLPUAyOj8WRSRWRSZXKaS5TQj0FoGkKXr/oir0r9P6K2PMC23IhtG0btqnC0NNQJB3x+RmY\n+oLop2FoGbfns77EURQLhveB4bzO7HzO6+67ZbwPDCuvCwt/w4m5bZp45c8+Db6rG4P3f33Nnz/0\nzix++3/fhSCy+JP/egE6e5q7hlwvxZCeeQGF5AEANhg+AF/HByGH31fXp3I7icGpcKbvL5lVcXAs\ngYNjSRwcTyBdNbko4hextS+Abf1+bO0LNGS8vd2/P8u2MF+IYSI3hRPZKUxkJzGRm6yZPQ8AXk5B\nt9KFHrkTPXIXepQudMudENnWe3Y7VaqdxpjpNLSZaWjT09Bm3XRmGkY8vuhzjjXf4Vj0USflOjrA\nd3SCDYVWPQGv2Sy9NM1CLqM6Ap92BD5bJfalYr0ocM4SO8UnwOsTofhFeH0CFJ8Ir5s/05b90t4X\nHTe5tRa9a+nrWZh6phzHoj6UK/besvCz1WLv7lNNdsm84cS8OD2N17/wZXg/+CF03/n5NX12dGQO\nv3nqIDiexSdv24VoV/OsC700XyXiACd2wtf1EXgC5y7bymt3MViJVt6fbduYiuVxcCyJd8eTOHwi\nhXyp0s2qSBy29vldgQ9goHPtAWE24vdn2zZSahoZOoEDk6M4kZ3CVG4asVJi0bkhMegIvNKNblfo\nO+UOcG3ge341352lqtDnZl2Rn4E2PQVtbg767MyiyXcAAIYBF4mC76gVea6jA1w4ckajzp3qb1NT\nDWTTpSqLvoRs2rHos5nSkl34FAUoXsEVehGK3xV+nwiv3xH+Ro7Zn1aQKst0BF/PlAXe1DLOsjw9\n4zYGsgAWLydcgGakxVY+5wXDKWA4BTTnPS0rfzVivv6fsir0tLPUYa2xy48dmsezT78LlmPwiT+7\noGlCrhfnkZ55HoXUOwAcV6v+risg+c9peZfU2Q5FUeiNKuiNKrj2kn5Yrrgfnkjh8Ik0Dp1I4Y3D\nMbxx2HGtynM0hjq9GOr2YXOPD5u6fYj422eMrVFQFIWgGMC2aD8G+c3l8pKhYrYwh8ncDKbzM5jK\nzWAqP4MD8REciI+Uz6MpGhExhE45ik5PBzo9HeiSo+jwRKG02ZJLWhDKTmyqsW0bVj4PbW4W+tws\n9Lk5Nz8HfW4O+f1vL74YRYENhcBFouDCEXBRJ2UjEXCRCNhAcF103/MC64aOre+NTy0ZZWHPZVRX\n6FV3v4TpiTSmka77WVHioPgEKD4BsleA4nVSWXHLFOGM+PugaMZZFicEsJQrroXwtmXBd8Xe1DMw\ntGw5SI5enF3uL4Fm5bLAl8WerRV+hlNOaSVTW4m5kXVaXoy8ejeP40fi+PW/HQTNUE3rWteKc8jM\nvHCSiF8Jyb/trHv5tws0RaEvqqAvquCqC/sAALF0EYdPpB2Bn0y7Il95ESkSh03dPmzq9rqpDz65\nvVzONgqRFdxlbv015Tktj+n8DCbzM5jOzWAqP4vZ/Bz2x97Ffrxbc67Cyej0uCIvR9Hl6UCHJ4qw\nGATTRkFAKIoCoyiQFAXS5uFFx81CHvrcvGPVV4m9EY+heOg9FO2RxRdlGHChMLhIFGwkXCP4XCQK\nxudbF2IviCwEcWmxN00L+WytyFfnk/ECYrP147wDTmNiQexlhXcE3xV6xeuIPn8GvENSFFUW2uWe\neMtU3S79DEw9B8vIwdRzrvBnYeo5GKUY9OLMsn+PZj2uyCtgeB+i0dtXrGN7irmyuhb91PEU9jx5\nADRN4U9uvQDdff6G1scR8edRSB0EAHBSN/zdV0DyERFvRyJ+CRG/hA/u6AIAlDQD4zNZHJvO4th0\nBsemM9h/NI79RytjpwGFR3+HFwOdCnYMR+GXGHQGPS2fCdwqFF7GVn4YW4O1opbT8pgpzGGuMI+Z\nwhxm8/OYLczhaHoco+mxmnNpikZIDCIqhRGVwohUpREpvOwSuvUI45HBDMkQh4YWHbN0HUYiDj0W\nc7d5GHEn1WMxFN59p+41KY4DGwqDC4XABkNgQ87m7IfBhkLrYgY+w9DupLn6dVlYcpfPqshl1Zp0\nIZ/LlJCYX3pcm2WdvyF6OMheHh5ZgEfm4VF4yIqTemQBvMA0/b1MMwJoKQpOii55jm3bsC21SuSr\nUqOSN7Qk9NKCpb/BxFzPLIj5yt3k8bkc/vOX+2HbwI27d6JnoHHrt7XiLNIzz6OYciwNXuqGv/tK\niL6tRMQ3ECLP4pyBIM4ZqKx4yBQ0jLniPjadwfG5XFngn3llHIDTRd8XVdDf4Wy9ERndYRleD3fW\n/j4UXsYWfhO2BDbVlOuWgVgxjtn8HGYKjsDHinHMF+J4N3HoJFveISD464h8CCExCIVr/hLERkJz\nnDNDvrOr7nFLVaHH44tEXo/FYMTjKMwubeHRkgQ2FAYbDCHd2wlD8oINuoLvNgJO1w3u6UJRFESJ\ngyhxS1r3AKBrBnJZra7Y57MqinnHP/5ysCztCntF4GsFn4dHESA1+TmlKAoUI4JmRHBiZNlzLVOD\nZS5eXVGPthLzimW+fDd7LlPCMz9/G5pq4po/PRf9mxoTaUorzDginna6xXhPD/xdVxARP4vweXhc\nMBzGBcMV3wS5oo6J2SwSBQMjR2M4PpfD+EwWR6dq3VkqEoeesAfdERk9YRk9ERndYQ+CZziS1nqC\no1l0y53olhc7bSoZJcwXE4gV447AF2OYLyYwX4jhSOoYDqeOLvoMT3MIiUF3CyAshhASAwhJTpmP\n967LePFLQQsChJ4eCD314zJYmgYjmYCRSEBPxGEkEjCSCegJp8xIxKFNnkDhwBLXl2WwgSDYQACs\nPwA2EABTlWcDQbB+/xmdrFcPjmcRDLMIhj11j0ejXsxMp1HIa8jnVBRyGgp5rSpVkXfzs1OZusvx\nFqAo1Ai+JHOQZB6Sh4Pk4eGRnVSSeYgSC7qJwx00w4NmVtfgaisx110xp5cZMy8Vdfz7z95GPqvh\ng1dtxrYdp+/Zrb6IXwnRt+WsfQkTKigSh3OHQs6M2vM6ADhr3qdieZyYz2EqVsB0PO9MuJusHYcH\nAJFn0BOR0RXyoCMooTNYSddrUJkzgciK6Pf2oN+7WMh0y0DcFfr5YhzxYgKJUhKJUhLxUhIzhbm6\n12QpBgExgHCV4AfFIIbMLlAlDn7BD6mNltjRPL+sZQ8AZrEIH0qYG52oEnlH9I1Uqiz4y8F4vY7Q\nl0U+ANYfrJT5/WB8XtBc6yx9hqWdpXEr+JCwLGftfVnga4TfaQjkcxoSsQLmZ5Yez19AlDhH8D0V\nwa/Zr2oINLOrv63eFEbGsXQYb30xNwwT//nLA0jGCjj/4l7surS/7nmrRStMuyL+HgCA9/Q63ene\nYSLihGVhGRoDnV4MdNYOCemGiZlEsSzuU7E8puKFupY84DQUOoMSOoIeJw1J6Ah4EAmI8Epnb7c9\nR7PokjvQJXfUPV4ySkiUUoiXEkiUUmWRXxD895JHaj9Q1Z8vMDwCQgABwYeA4Hc3H/yCH0HBD7/g\nh5eX28bCZyQJnmgHZGlpB1mWqjrCnk7BTKXcfNJJ3XJ9fh7qxMSyf4uWJDBeHxifD6zPV5V3xN5J\nnXJaklry+6VpyrG8ZR6RZWw927ahqQaKBR3FvOakBQ3FvJu65YWCjkJOQzJWWPpiC3+bocoiL0oc\nRA8HSarKL5S7+6LEgVnlEtm2EnM967SSGM/iCXC2bWPvv49g5kQaw9uj+PA1p241a4UpV8QPAXCC\nn/i7roTo3XzWvjwJjYFjmfJYejWGaSGeKWE2UcRcsoDZZBFzySJmkwUcm85itI7Q8xyNsE9ExC8h\n7BcRcbewX0TEJ8In82ft71VkRfQojkObemimjmQpiUQphaSags6WMJmYQ0rNIKWmkVLTmF3Cugec\nSXp+3ge/4IOP98LHK/DxXnh5L3yC1y1zyvlVdpO2EloQwHd2gu9cvifTKhVhpNIwUkkYaSc1UykY\nmQzMbMZJM2no83NYti8bjsMdxucD4/PXCr/XC0bxglYUMIoCRvGCUZQzLv4URUEQOQgih0Cofvd+\nNaZpoVTUF4n9okZAXkMqUYChL71uvRpBZHH3N25c8by2EnMjkwHt8dT1qvTHl8cxOjKPrj4/rv7E\n9lP60tX8JNIzz6OUOQwAEOR++LquICJOaDosQ6Mz6EFn0AOgNlbAgtDPJYuYTRQwlyoini45W6ZU\n47a2Go51xD7sExDwCgh6BQSV2rxX5kGfhb9tnuHQKXeg07Xs6zkd0U0daS2DZCmNtJpGUk0jXSX2\nKTWDiewkTLu+45QFREY4SeiVsth7eQUKJ0PhFCi8B+I69xdOixL4Lgl819Ld+oDjE9/M52BmMjAz\nmbLYO/l0TZk2eQLqWH1PczUwDBhZLot7tdBrXWGUwLe0AcAwNGTFWTa3GgzdRKmol7diwc0Xast0\nffnf1wJtJeZ6Nlt3JvvR9+ax74UxKD4B19+8Ayy7tjWqau440jMvoJQdBeCIuL/7SgjKpnX9YBHO\nDqqF/vzNi4MCFUoG4hlH3GPpIuKZEmJpZ4unS5hJLN39x9AU/ApfEXnFEfqAV4Bf5uFzt7DVNo4i\nGwbHcOXlcEth2zYKRhEZLYuMmnVSLYuslivnF47NF+Owl/ATvgBDMVA4D2ROhsIrUDiPI/ScBzIv\nw8vJzjFOhsI7KbsOPexRNA3W6wPr9QErRJa2bRtWqVQWfjOXgZnLuVsWZi7vps6+kU5Bm5qsucZi\nf4QuDANG8oCWPWA8MmiPB4wsg/bIYDwe0LKbemS3vHKcFpvbsGI5BgrHQPE1Zo7G+vsVLIORzUEY\nqJ2ZHp/L4bf//i5YjsaNt+yEZ5VOPGzbhpobQ3rmBai5MQCAoAzC33UFBGWIiDihbfCILDzi4q77\nBVTdRCqrIplVkcypi/M5FWMzWZh1uvIXoClA8fDweXj4ZQ4+mYfXw9cIvs/Du+Xcml3htisURUHm\nPJA5T90Z+dWYlomcXqgR+ayWRU7PI68VkNNzyOkF5PQ8EqUUpvLLOxZZQGREyJwHHk6CzHogcRI8\nrLPJnAceVkJXKQS9gPK+h5PWTS8ARVFgJMlZF79CN/8CtmnCLORhuaIv0waSk/NOb0AuCzPrpvk8\nrHweZiEPIx5fFClvWWjaEfeqRoAj/B7Qkge0JDlzBCQJtCiV951jIhjJc0ZXAbSVmNuGUeP9rVjQ\n8J+/2A9Dt3DdTTsQ6Vx5/blt2yhlR5GZeQFq3pnMIXo3O93pysAKnyYQ2g+BY9AZ8qBzmXE/y7aR\nzWtI5hyhT2VVpPMasgUdmbyGgmYini4ininixPzKM3wlgYUisVAkDrLIOankpE6ehSLWlol88516\ntBKGZuAXvKuOOmdYBvKuuOe0vJO6W/6ksrxewGx+Dpq1iu5qF5qiIbGiK/xOY8AReg8kVoTEiBBZ\nASIrQmJFiIwAkZUgsQJExilrlac+imEcq9/rePQMR72wVvDNbts2bE2DWSjAKuQdoS8U3DTvlLvC\nXykvOA2BxBobAgv1ZNka4S83AERpURktSaBF0WkYiCJoQXTSVTr/aSsxByoOY0zTwq+ffAfZjIqL\nPzKE4e1Le9wBANs2UUi+i+zcK9CK0wAA0bfVscTlFfqBCIQNDk1R8CsC/IqAoTrDodVjyppu2esT\nmwAADblJREFUIlPQkMnrblq1ufu5ooF8SceJ+Tx0Y3UTfRiagixxkEUWssTBI7DwCCwkgYVHdNNl\n9nmW3lCNAZZm4RecSXarRbcMFPQiikYBBaOIvF5AQS+CFi3MpVIo6E55QS+6qbOfVNMwrLWLFUez\njtgzIkRWdPNCJV9uBDgpz/AQGB4CI1RS1klZqrmNOYqiQAkCaEEAgmsPfW1pWkX0i0VnKxVhFZzU\nXChbVF6AVSzBSCVha0vFX1+24uj8t1+seFrbiTntOoz5w3PHMDWRxuZzIrj4w4NLnm+ZKnLx15Gd\n+4Mb2xaQAufC3/kR8J7uM1JnAmEjwXNM2fXtalB1E/mijpy75UtGJV9dXtSRKxnIFnTMJAorTYZe\nBENTFYEXnVTkGYj8QspA4BmIHANRYCFwTLk8rZoo5lXnHM45rx0nBnI0W9f6X01UMc3UUTAc8S+Z\nJRSNEkqGm5pqJW+oKJrVx5w0qaahr6Fn4GRoiq4S+Tqi76YnNwh4hkdE86GYM8AzHDiaB0+z4Bke\nHMOBpzmwNHvaSwlpngfN82ADa28ILGCbZkXsi1UNgLL4l2CpbloswSoVYZsbcAIcADCyjGOHYnjr\n1Qn4QxKu+nj9meuGlkJ2fh9ysddhWyoomoMSuQTejsvACY3xCEcgEFZG4ByBDK1hoo9t2yhpJoqq\ngYJqOGnJqN1XDRRLTnryOam4Cm2VS3+Wq3e5AeA2AvjyRoNnGfAsXbMvcHTlHJauOr74HIam1lVP\nAs9w4BlnXf2pYlqmK/RquTGwIPaqqUE1VWim5uadfdXUoBoqNEsr5wt6EclSak3DBivB0Sx4uiLw\nldQRfyetlJ98jKVZcDRbSSkWHOOkThkHlmZqUpqq7S2iGKY8E7/RNFXMbdvGP/zDP+C9994Dz/P4\nxje+gf7+iiOXvXv34uGHHwbLsrjllltw6623rnjNIufF3mdGwLA0rr9pR03EHNsyUcwcQi72enlm\nOs3K8HV+CErkYjBs6wMPEAiElaEox8qWBBan2vQ2TAslzYSqmShpBkq6WbPvpM5GsTSSqSJUvVJW\nfU46r0HVVmchrRaaosBxNASWBsfSYFkGHEODYylwDA2WpStpdZ6pv8+elHIMBY5lwLE0CqaNbKYI\nlqbBMBRYhgZDU2AZCgxDN6wXgqEZKLTcsPC2lm1BM/Va4a9qEJRMFbqpgZcYJDNZaJYO3dShWRo0\n04BuadBMHbqlV6XO57J6Dpqpr7i08HSgQIGtbgBUpfUbBRw4mgHrNghYmoXA8Phs9KYV/1ZTxfzZ\nZ5+Fpml4/PHH8dZbb+Fb3/oWHn74YQCAYRh48MEH8cQTT0AQBNx222245pprEAot/ehaoPHSKANN\nNXDVx89BuENxvPQUp1FIvoN84m1YhuNsn5f7oIQvhBzcCWodLt0gEAjNhWVoKBINRVo5ytpquqFt\n24ZuWNAMC5puQtVNZ1+3oBomNN2EpjvHNMOCZlTtn3SObphQy9exYBgWiqqBrGFBN61VzzNoFBTl\n/L9YhgKzIPhLCD9Ln5QyFJhyvnLuQp6hKdC0c126Zp8CTZ20f1Keoar3eTC0AJH2QaYp0CwFmneO\nRaNeJKV83Wuu1PthWqYj8uWGgCP4umW4qdMQMCwDumXAsA0YpgHdNmBYRqXcMqBbOgzLhGHpVWW1\n5xX0Igy7Ur4aPntJi8X8j3/8Iy6//HIAwK5du3DgQMXb/+joKAYHB6G43Q0XXXQR9u3bh+uvv37J\n6x2OXIJExsK5FwQxMJBB8sTbKKRHYGqOr2uakeCNXgY5/H7wUn03jwQCgXAqUBRV7jLHKhoIp4Nt\n2zAtp/Ggm47YL4i8blgwzNpUL++7nzHMct4wLTAcg2xOhWnaMC3nPNO0YFi1qWnaMCwbhptXdROF\nklH1GRvWWicztBiKwqIGBEVRoCmAcvdpCm4Z5ZahznkATXGgaB40HLewi69DgXI/y9IU+CWv4/wd\nCgBoC6CczabM2hQWOHYduHPN5XLweisTMViWhWVZoGl60TFZlpHNLt8y7r7GwDZ5HySpiNgxp4yi\nBXiC58MTOBeSbwuxwgkEQttDURXrthGDg6vpeVgtllXVIFhoDJg2DMttDJiWW17JW27jxHI3szq1\nT9p3r7/c8ZpzbRscx6JY1GrK6v2thevZtnMftu3sW5azPNMwLeeY+/lyvuYzOOMNmrv+y8rnNFX5\nFEVBPl+JMbsg5AvHcrnKetV8Pg+fb/klGDfefm9zKrqOiEZXtwa1XSH3175s5HsDyP0R2pumumm6\n8MIL8dxzzwEA3nzzTWzbtq18bHh4GOPj48hkMtA0Dfv27cP73ve+ZlaHQCAQCIQNCWXbzesvqJ7N\nDgDf+ta38M4776BYLOLWW2/F7373O3z/+9+HbdvYvXs3brvttmZVhUAgEAiEDUtTxZxAIBAIBELz\nOTuiIRAIBAKBsIEhYk4gEAgEQptDxJxAIBAIhDaHiDmBQCAQCG1OW3hYWcnH+0bgrbfewj/+4z/i\nsccea3VVGophGLj33nsxOTkJXdfxhS98AVdffXWrq9UwLMvCfffdh2PHjoGmaXzta1/Dli1bWl2t\nhhOPx3HLLbfgX/7lX7Bp06ZWV6ehfOpTnyp7ouzr68M3v/nNFteosTz66KPYu3cvdF3H7bffjltu\nuaXVVWoYTz75JJ544glQFAVVVTEyMoKXXnqp/H22M4Zh4O6778bk5CRYlsUDDzyw7LPXFmK+nI/3\njcAPf/hDPPXUU5DlxgQnWE88/fTTCAaD+M53voN0Oo2bbrppQ4n53r17QVEU/vVf/xWvvvoq/umf\n/mlD/TYB56Vy//33QxRXH/WsXdDc+NI/+clPWlyT5vDqq6/ijTfewOOPP45CoYAf/ehHra5SQ7n5\n5ptx8803AwC+/vWvY/fu3RtCyAHgueeeg2VZePzxx/Hyyy/jn//5n/G9731vyfPbopt9OR/vG4HB\nwUE89NBDra5GU7jxxhtx1113AXCsWJZti/bjqvnYxz6GBx54AAAwOTkJv//Uw0euV7797W/jtttu\nQ0fHxot3MDIygkKhgDvvvBOf+9zn8NZbb7W6Sg3lxRdfxLZt2/DXf/3X+OIXv4irrrqq1VVqCvv3\n78eRI0dWFXmzXRgaGoJpmrBtG9lsFhy3fDyAtnizLufjfSNw7bXXYnJystXVaAqS5HiWzuVyuOuu\nu/A3f/M3La5R46FpGl/96lfx7LPPLttybkeeeOIJhMNhfPjDH8YPfvCDVlen4YiiiDvvvBO33nor\nxsbG8Jd/+ZfYs2fPhnm3JJNJTE1N4ZFHHsHExAS++MUv4le/+lWrq9VwHn30UXz5y19udTUaiizL\nOHHiBG644QakUik88sgjy57fFr/Y5Xy8E9Y/09PTuOOOO3DzzTfj4x//eKur0xQefPBB7NmzB/fd\ndx9KpVKrq9MwnnjiCbz00kv47Gc/i5GREdx9992Ix+OtrlbDGBoawic/+clyPhAIYH5+vsW1ahyB\nQACXX345WJbFpk2bIAgCEolEq6vVULLZLMbGxnDppZe2uioN5cc//jEuv/xy7NmzB08//TTuvvvu\n8rBQPdpCEZfz8b6R2IjO+GKxGO6880585StfKY9tbSSeeuopPProowAAQRBA0/SGamj+9Kc/xWOP\nPYbHHnsM27dvx7e//W2Ew+FWV6th/PKXv8SDDz4IAJidnUU+n0c0Gm1xrRrHRRddhBdeeAGAc3+l\nUgnBYLDFtWos+/btwwc+8IFWV6Ph+P3+8vi/1+uFYRiwrKXj3LdFN/u1116Ll156CX/+538OwPHx\nvhGhKKrVVWg4jzzyCDKZDB5++GE89NBDoCgKP/zhD8HzfKur1hCuu+463HPPPfjMZz4DwzDwd3/3\ndxvm3k5mI/4+d+/ejXvuuQe33347aJrGN7/5zQ3VGPvoRz+K1157Dbt374Zt27j//vs33Pd47Nix\nDbe6CQDuuOMO3Hvvvfj0pz8NwzDwt3/7t8tOQiW+2QkEAoFAaHM2ThOUQCAQCISzFCLmBAKBQCC0\nOUTMCQQCgUBoc4iYEwgEAoHQ5hAxJxAIBAKhzSFiTiAQCARCm0PEnEAgEAiENqctnMYQCIT1wejo\nKI4ePYoDBw7g85///IaM9EcgtCPEMicQCKvmN7/5DbZv345du3aVXSwTCITWQzzAEQiENfOzn/0M\nW7ZswYUXXtjqqhAIBBDLnEAgrJFSqYR8Pk+EnEBYRxAxJxAIa+LFF1/EHXfcgYmJiVZXhUAguJBu\ndgKBUMP4+Dj27duHV155Bd/97ncBOBPfcrkcRkdH8fOf/xySJOGv/uqvcNlll7W4tgQCASBiTiAQ\nluBLX/oS/uIv/gK9vb149913ceWVV7a6SgQCYQmImBMIhLrs3bsXzz77LK644grccMMNra4OgUBY\nBiLmBAKhLtlsFjfeeCNefPHFVleFQCCsAJkARyAQ6vLb3/4WO3fuxKFDh1pdFQKBsAJEzAkEwiKe\neeYZXHfddfjEJz6BX/ziF62uDoFAWAEi5gQCoYYXX3wRl156KTweD6677jpMTk7iqaeeQiKRaHXV\nCATCEpAxcwKBQCAQ2hximRMIBAKB0OYQMScQCAQCoc0hYk4gEAgEQptDxJxAIBAIhDaHiDmBQCAQ\nCG0OEXMCgUAgENocIuYEAoFAILQ5RMwJBAKBQGhziJgTCAQCgdDm/H/8NQO+F2UW8QAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11a95f940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0,8,300)\n",
    "\n",
    "chi2_df1_pdf = stats.chi2.pdf(x, df=1) \n",
    "chi2_df2_pdf = stats.chi2.pdf(x, df=2)\n",
    "chi2_df3_pdf = stats.chi2.pdf(x, df=3)\n",
    "chi2_df4_pdf = stats.chi2.pdf(x, df=4)\n",
    "chi2_df5_pdf = stats.chi2.pdf(x, df=5)\n",
    "\n",
    "plt.plot(x, chi2_df1_pdf, label=\"df = 1\")\n",
    "plt.plot(x, chi2_df2_pdf, label=\"df = 2\")\n",
    "plt.plot(x, chi2_df3_pdf, label=\"df = 3\")\n",
    "plt.plot(x, chi2_df4_pdf, label=\"df = 4\")\n",
    "plt.plot(x, chi2_df5_pdf, label=\"df = 5\")\n",
    "\n",
    "plt.xlabel('$\\chi^2$')\n",
    "plt.ylabel('Density')\n",
    "plt.ylim(0,1)  # note this will cut off the curve for df=1\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFrCAYAAADFOmBlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VOW9//HP3nNNMpMrCUrEcC+iLSpWpRUIUCxoqlCI\nB6xAj7Reqj+pus7yoJ5asUcore065xTWgcOp/LTnJ/VCheqy9ChKBYuk3Cr3i4AhIAQCSSa3yVx+\nf0wySbgkATIMe+b9WitrZvbes+f7GMxnnmfv/WwjHA6HBQAALMuMdwEAAODiEOYAAFgcYQ4AgMUR\n5gAAWBxhDgCAxRHmAABYnD2WOw8EAnr66adVVlamxsZGPfTQQxo1alR0/ZIlS/Tmm28qOztbkjR7\n9mz16tUrliUBAJBwYhrmK1asUFZWlubNm6fKykqNHz++TZhv27ZN8+bN06BBg2JZBgAACS2mYT5u\n3DiNHTtWkhQKhWS3t/24bdu2aeHChSovL1dhYaEeeOCBWJYDAEBCimmYp6SkSJJ8Pp9mzpypxx9/\nvM36O++8U9/73vfk8Xj0yCOPaPXq1RoxYkQsSwIAIOHE/AS4I0eOaPr06ZowYYLuuOOONuumT5+u\nzMxM2e12jRgxQtu3b293X8w8CwDAmWLaMz9+/LhmzJihn/zkJ7r11lvbrPP5fCoqKtJ7770nt9ut\ndevWadKkSe3uzzAMlZdXx7LkuMrN9dI+C0vk9iVy2yTaZ3XJ0L6OxDTMFy5cqKqqKi1YsEDz58+X\nYRi65557VFdXp+LiYj3xxBOaOnWqXC6Xhg4dquHDh8eyHAAAEpJhtbumJfq3L9pnXYncvkRum0T7\nrC4Z2tcRJo0BAMDiCHMAACyOMAcAwOIIcwAALI4wBwDA4ghzAAAsznJh3hgIxrsEAEAS8fv9Ki6+\nS5K0fftWTZ48QQsXzr+gfa1e/aGef/7ZrixPksXCvLa+UT/+jzV6968H4l0KACBJRKZjMSRJ69ev\nU3HxFD344CPnvZ9/+7eX9F//tSAmU5PHdAa4rlZV41ddQ1BfnqiNdykAgEvk9VV7VbLz2DnX22yG\ngsHzC8ivD8zTPaP6nXN9XV2dZs9+VtXV1crPv0qStHPndr377go5HA7l5eVp2LBCSVJZ2SHNnfuC\nDMOIvn/MmLH6znfGt9nnV786WMOHF2r58mXnVWtnWCrM7bbIQELQWpPWAQAs5u2331KfPv30wx8+\nrO3bt2rjxg0aOHCQxo0rUk5Ot2iQS1J+/lX6j/9Y2OE+R436ljZt2hCTei0V5jYz8q0nFCLMASBZ\n3DOqX7u96FhM51paelDf+MYwSdKgQdfJbredc9vWPfNwOCzDMM7aM48lS4W52RTm5zucAgDA+ejV\nq4+2bv27brttuHbv3qlAOydfd7ZnHkuWOgEuOsxOzxwAEEPjx0/U4cNleuSRH+oPf3hTLpcz3iW1\ny1I98+gwO8fMAQAx5HQ6NXv2nDOW33//Axe13xtuGKIbbhhyUfs4G0v1zG225mH2UJwrAQDg8mGp\nMDdNhtkBADidpcLcZhoyxNnsAAC0ZqkwlyJntHOdOQAALSwX5jbz/Gf6AQAgkVkvzG0Gw+wAALRi\nuTA3DYbZAQBozXJhzjA7AABtWS7MTZNhdgAAWrNcmNtMk+vMAQBoxYJhbjCdKwAArVguzE3TYDpX\nAABasVyY22wGw+wAALRivTA3GGYHAKA1y4W5yaVpAAC0Ybkwt5kMswMA0Jolw5zrzAEAaGG5MDdN\nQ2FxG1QAAJpZLsxtpiFJDLUDANDEemFui5RMzxwAgAjLhblp0DMHAKA1y4V5yzA7s8ABACBZMcxt\nkTBnmB0AgAjLhbnJCXAAALRhuTC3ccwcAIA2LBfmzT1zhtkBAIiwXJg3X5pGzxwAgAjrhblBzxwA\ngNYsF+acAAcAQFuWC/PmS9MIcwAAIqwX5pwABwBAG5YL85bpXJkBDgAAyYJhzl3TAABoy3phznSu\nAAC0Ybkw52x2AADaslyYM50rAABtWS/Mm2aAY5gdAIAIy4U5w+wAALRluTBvOZudS9MAAJAsGOYm\nx8wBAGjDHsudBwIBPf300yorK1NjY6MeeughjRo1Krp+1apVWrBggex2uyZOnKji4uIO98mlaQAA\ntBXTMF+xYoWysrI0b948VVZWavz48dEwDwQCmjt3rpYtWyaXy6UpU6Zo9OjRys7ObnefTOcKAEBb\nMR1mHzdunGbOnClJCoVCsttbvjvs27dPBQUF8ng8cjgcGjJkiEpKSjrcZ/Mwe4AwBwBAUox75ikp\nKZIkn8+nmTNn6vHHH4+u8/l88nq90ddpaWmqrq7ucJ8MswMA0FZMw1ySjhw5okcffVT33Xef7rjj\njuhyj8cjn88XfV1TU6P09PQO95edlSZJSklxKjfX28HW1pOIbWqN9llXIrdNon1Wl+jt60hMw/z4\n8eOaMWOGfvKTn+jWW29ts65v3746ePCgqqqq5Ha7VVJSohkzZnS4z+rqOklSZXW9yss77slbSW6u\nN+Ha1Brts65EbptE+6wuGdrXkZiG+cKFC1VVVaUFCxZo/vz5MgxD99xzj+rq6lRcXKxZs2bp/vvv\nVzgcVnFxsfLy8jrcZ/N0rgyzAwAQEdMwf+aZZ/TMM8+cc31hYaEKCwvPa5/N07lynTkAABHWmzSG\nS9MAAGjDcmHOdK4AALRl4TCnZw4AgGTBMGeYHQCAtiwX5vTMAQBoy3Jhzv3MAQBoy3Jhzo1WAABo\ny4JhHimZMAcAIMJyYd48zM5d0wAAiLBcmDPMDgBAW4Q5AAAWZ7kwbxlmZwY4AAAkC4Y5PXMAANqy\nXJgzAxwAAG1ZL8wNQ4bB2ewAADSzXJhLkaF2euYAAERYNMxNpnMFAKCJJcPcpGcOAECUJcPcZhr0\nzAEAaGLJMDcJcwAAoiwZ5pET4Jg0BgAAycJhTs8cAIAIS4Y5w+wAALSwZJhznTkAAC0IcwAALM6S\nYW6aBtO5AgDQxJJhbjNNeuYAADSxaJgzzA4AQDNLhnnz2ezhMIEOAIAlw9zWfE9zwhwAAIuHOUPt\nAABYM8zNpjBn4hgAACwa5jbCHACAKMIcAACLs2SYmxwzBwAgypJhbrdFym4McBtUAAAsGeZOh02S\n5CfMAQCwaJjbm3vmwThXAgBA/FkzzB2Rsv2N9MwBALBkmDvszcPs9MwBALBkmLvs9MwBAGhmyTB3\nOOiZAwDQzJJh7qRnDgBAlDXDvKlnznXmAABYNcyjPXOG2QEAsHaY0zMHAMCiYd58Ahw9cwAArB7m\n9MwBALBmmEeH2emZAwBgzTDnRisAAERZMswdnM0OAECUJcPc5eBsdgAAmlkyzO02U4akRnrmAADE\nPsy3bNmiqVOnnrF8yZIlKioq0rRp0zRt2jQdOHCg0/s0DEMOh6kGeuYAAMgey50vXrxYy5cvV1pa\n2hnrtm3bpnnz5mnQoEEXtG+n3cZ0rgAAKMY984KCAs2fP/+s67Zt26aFCxfq3nvv1aJFi857306H\nyQlwAAAoxmE+ZswY2Wy2s66788479fzzz+uVV17Rhg0btHr16vPat8Nu4wQ4AAAU42H29kyfPl0e\nj0eSNGLECG3fvl0jRozo8H25uV5JUprboaoaf/R1oki09pyO9llXIrdNon1Wl+jt68glCfNwONzm\ntc/nU1FRkd577z253W6tW7dOkyZN6tS+ysurJUmGEZa/MRh9nQhyc70J1Z7T0T7rSuS2SbTP6pKh\nfR25JGFuGIYk6Z133lFdXZ2Ki4v1xBNPaOrUqXK5XBo6dKiGDx9+Xvt02m0KhsIKBEOy2yx5hR0A\nAF0i5mGen5+vpUuXSpKKioqiy++66y7dddddF7zflnuaE+YAgORm2RRsnp+9kZutAACSnHXDvKln\nzsQxAIBkZ90wb+6Zc605ACDJWTjMudkKAACShcPcYW+6pzk9cwBAkrNsmHMbVAAAIiwb5i09c8Ic\nAJDcLBvmLcfMGWYHACQ364Z506Vp3AYVAJDsLBzmkWH2Bk6AAwAkOeuGefMwO2EOAEhynQrzxYsX\nq7y8PNa1nJfmnjnD7ACAZNepMK+vr9d9992nBx54QO+9954aGxtjXVeHHI6WG60AAJDMOhXmjz76\nqFauXKkHHnhAn376qe6++27Nnj1bO3bsiHV95+RqvjSNs9kBAEmu08fM6+rqdOjQIZWWlso0TaWn\np+tnP/uZXnrppVjWd04uZyTM6/2EOQAguXXqfuZPPvmkPv30Uw0fPlwPP/ywbrrpJkmS3+/Xbbfd\npieffDKmRZ5NiitSel1D4JJ/NgAAl5NOhfnQoUP1wgsvKDU1NbrM7/fL6XTq3XffjVlx7UlxRXrm\ntfWEOQAguXVqmP2NN95oE+ShUEgTJ06UJOXm5samsg7YTFNup42eOQAg6bXbM582bZrWr18vSRo4\ncGDLm+x2jRo1KraVdUKq265awhwAkOTaDfNXXnlFkvSzn/1Mzz777CUp6HykuuyqqGqIdxkAAMRV\nu2H+4YcfauTIkbr22mv19ttvn7F+/PjxMSusM1JddpU11CgUDss0jLjWAgBAvLQb5p999plGjhwZ\nHWo/XbzDPMVlV1hSfUNQqe5OncsHAEDCaTcBH3vsMUnSnDlzost8Pp+OHDmi/v37x7ayTmgO8LqG\nAGEOAEhanT6bfdasWaqoqNAdd9yhxx57TL/+9a9jXVuHUl0OSeIkOABAUutUmL/22mt66qmn9M47\n72j06NH64x//qI8//jjWtXUopak3Xlsf/7niAQCIl05P55qZmanVq1ersLBQdrtdDQ3xP4s8tWkW\nOHrmAIBk1qkw79evnx588EEdOnRIQ4cO1cyZM3XdddfFurYOpUZ75oQ5ACB5deqssRdffFGbNm1S\n//795XQ6dffdd2vEiBGxrq1DqczPDgBA58K8trZWu3fv1vr16xUOhyVJ27dv16OPPhrT4joSPWZO\nmAMAklinwnzmzJnyer3q37+/jMtocpboMXOG2QEASaxTYX78+HG9/PLLsa7lvKXSMwcAoHMnwF1z\nzTXauXNnrGs5bxwzBwCgkz3zPXv2aMKECcrJyZHL5VI4HJZhGPrggw9iXV+7UhhmBwCgc2H+m9/8\nJtZ1XBC7zZTLYWOYHQCQ1Do1zJ6fn6+NGzfq9ddfV3Z2tkpKSpSfnx/r2jolxWVTHT1zAEAS61SY\n//KXv9Tq1av15z//WcFgUG+99Zbmzp0b69o6JdXtoGcOAEhqnQrzNWvW6Be/+IVcLpc8Ho9efvll\n/eUvf4l1bZ2S6rKrriEQvf4dAIBk06kwN822m/n9/jOWxUuq265gKCx/YyjepQAAEBedOgFu7Nix\n+vGPf6yqqiotWbJEy5cvV1FRUaxr65Tma819dY1yOW1xrgYAgEuvU2FeWFiovLw8lZaWasOGDZo5\nc6YKCwtjXFrnpKc6JUlVtX7lZLjjXA0AAJdeu2F+4sQJPfbYY9q7d68KCgpkt9u1bt061dfX68Yb\nb1R6evqlqvOc0tOawrzGH+dKAACIj3YPfL/wwgsaMmSI1qxZo9dff12vv/66PvnkEw0cOFAvvvji\npaqxXdGeOWEOAEhS7Yb5rl279MQTT8jhcESXORwOPfHEE9q+fXvMi+uMaM+8ljAHACSndsPc5XKd\ndblhGJfN2ewZTWFeSc8cAJCk2k3k9m53erncCpVj5gCAZNfuCXB79uzR6NGjz1geDodVXl4es6LO\nhzc1cgiAMAcAJKt2w3zlypWXqo4LZreZSnPbVVXbGO9SAACIi3bD/HK5mUpH0tOc9MwBAEnr8jiL\n7SKlpzrlq2tUIMiUrgCA5JMYYd50Elw1Q+0AgCSUUGHOUDsAIBklVpgzcQwAIAklRJhn0DMHACSx\nmIf5li1bNHXq1DOWr1q1SpMmTdLkyZP1xhtvXNRnMD87ACCZdeoWqBdq8eLFWr58udLS0tosDwQC\nmjt3rpYtWyaXy6UpU6Zo9OjRys7OvqDPSWdKVwBAEotpz7ygoEDz588/Y/m+fftUUFAgj8cjh8Oh\nIUOGqKSk5II/J9MTCfNTvoYL3gcAAFYV0zAfM2aMbDbbGct9Pp+8Xm/0dVpamqqrqy/4czI9LtlM\nQycq6y94HwAAWFVMh9nPxePxyOfzRV/X1NQoPT29U+/NzfWedXm3zBSd9DWcc71VWL3+jtA+60rk\ntkm0z+oSvX0duSRhHg6H27zu27evDh48qKqqKrndbpWUlGjGjBmd2ld5+dl78JlpTu0qPaXDRyrl\nsFvzJP3cXO8525cIaJ91JXLbJNpndcnQvo5ckjBvvl3qO++8o7q6OhUXF2vWrFm6//77FQ6HVVxc\nrLy8vIv6jJwMt1QqnayuV15WaleUDQCAJcQ8zPPz87V06VJJUlFRUXR5YWGhCgsLu+xzctLdkqQT\nlYQ5ACC5WHM8+ixyMiJhfryKk+AAAMklccK8qWdeUcXlaQCA5JI4YZ7RMswOAEAySZgwz/a6JEkn\nGGYHACSZhAlzp8Om9FQHPXMAQNJJmDCXIkPtFdX1Cp12XTsAAIkswcI8RYFgWKeqOQkOAJA8EirM\nr8hOkSQdraiNcyUAAFw6CRbmkcliviTMAQBJJKHC/MqcyH3TjxDmAIAkklBhHu2ZnyDMAQDJI6HC\nPMVlV0aak2F2AEBSSagwlyK98xOV9fI3BuNdCgAAl0TChfmVOakKSzp2si7epQAAcEkkXJhzRjsA\nINkkXpjnRML8yImaOFcCAMClkXBh3qPp8rSy44Q5ACA5JFyY52S4leKyq/SYL96lAABwSSRcmBuG\noavzPPryRK0a/JzRDgBIfAkX5pLUs7tHYUmHjtM7BwAkvoQM86vzvJKkL44S5gCAxJeYYd7dI0kq\nPVod50oAAIi9hAzzHt3SZDMNfcFJcACAJJCQYW63merRLU2HjvkUCoXjXQ4AADGVkGEuSQXdvfIH\nQlxvDgBIeAkb5v2uypAk7S2rjHMlAADEVsKGed/8pjA/RJgDABJbwob5lTmpSnXZte8wYQ4ASGwJ\nG+amYahPfrqOnaxTVY0/3uUAABAzCRvmktSvaah9H8fNAQAJLKHDvPm4+R7CHACQwBI6zPv1yJDN\nNLTj4Ml4lwIAQMwkdJi7nDb1zc/QF19Wy1fXGO9yAACIiYQOc0ka1CtLYUk76Z0DABJUEoR5tiRp\nO2EOAEhQCR/mva/0yu20aceBiniXAgBATCR8mNtMUwOvztLRk3U6dqou3uUAANDlEj7MJWlwvxxJ\n0pY9x+NcCQAAXS9JwrybJGnzXsIcAJB4kiLMMz0u9b4yXbu+OKWaei5RAwAklqQIc0m6vl+OQuGw\nPvv8RLxLAQCgSyVNmN/QP1eStGFXeZwrAQCgayVNmOfnpunKnFRt2XtCdQ2BeJcDAECXSZowNwxD\ntwzqrkAwpI276Z0DABJH0oS5JN1yTXdJ0qc7jsa5EgAAuk5ShXn37FQVXOHV9v0nVelriHc5AAB0\niaQKc0m67atXKhQO65OtX8a7FAAAukTShfmt13aX3WbqL38/onA4HO9yAAC4aEkX5mluh24amKuj\nFbXac6gy3uUAAHDRki7MJWn413pIklZtPBTnSgAAuHhJGeZfuTpT+blp2rCrXCerOREOAGBtSRnm\nhmFozE09FQyF6Z0DACwvKcNckm4d1F2eFIc+2lSmBn8w3uUAAHDBkjbMnQ6bRt6Qr5r6gFZvORzv\ncgAAuGAxDfNwOKznnntOkydP1rRp01RaWtpm/ZIlS1RUVKRp06Zp2rRpOnDgQCzLOcOYr/eUy2HT\nnz49qMZA6JJ+NgAAXcUey52///778vv9Wrp0qbZs2aI5c+ZowYIF0fXbtm3TvHnzNGjQoFiWcU6e\nFIcKb+ihletLteazIxp5Q35c6gAA4GLEtGe+YcMGDRs2TJI0ePBgbd26tc36bdu2aeHChbr33nu1\naNGiWJZyTmNvvlpOu6k/rt0vfyPHzgEA1hPTMPf5fPJ6vdHXdrtdoVDLcPadd96p559/Xq+88oo2\nbNig1atXx7Kcs8rwuPStm3rqlM+vDzizHQBgQTEdZvd4PKqpqYm+DoVCMs2W7w/Tp0+Xx+ORJI0Y\nMULbt2/XiBEj2t1nbq633fUXYuqdg7R6y2G9t+4LjR85QOlpzi7/jM6KRfsuJ7TPuhK5bRLts7pE\nb19HYhrmN954oz788EONHTtWmzdv1oABA6LrfD6fioqK9N5778ntdmvdunWaNGlSh/ssL6+OSa1F\nQwv0+1V7tWjZFk0fOzAmn9GR3FxvzNp3OaB91pXIbZNon9UlQ/s6EtMwHzNmjNauXavJkydLkubM\nmaN33nlHdXV1Ki4u1hNPPKGpU6fK5XJp6NChGj58eCzLadfoIVdpzd+P6C+bD2vY13qoT4/0uNUC\nAMD5MMIWu3VYLL997fripH7+/zap4Aqv/mXaTTJNI2afdTbJ8O2S9llTIrdNon1Wlwzt60jSThpz\nNl+5Oku3XttdB7+sZiIZAIBlEOanuWdkP7mdNr310T5VVNXHuxwAADpEmJ8m0+NS8ch+qm0I6L/f\n3aGQtY5CAACSEGF+FoXX99DgvjnacfCk3i8p7fgNAADEEWF+FoZh6B/vuEbpqQ69uXqfSo/54l0S\nAADnRJifQ3qaU/94xzUKBMNatGKb6v2BeJcEAMBZEebtGNyvm0YPuUplx2v03+/ukMWu4gMAJAnC\nvAP/MKqfBvTM1IZd5XrnkwPxLgcAgDMQ5h2w20z9aMJ1ykl36Q8f79em3eXxLgkAgDYI805IT3Xq\n/0z8mpx2U4v+uF2fH66Kd0kAAEQR5p10dXevHrjrWvkDQf369c0qO17T8ZsAALgECPPzcOOAXH1/\n7EDV1Af0q99v1vHKuniXBAAAYX6+hg3uoXtG9tPJ6gb98rXNOlHJlK8AgPgizC/A2Fuu1ne+0UvH\nTtVp7v9s0LGTtfEuCQCQxAjzCzRheB9NGN5HJ6oaNOd/Nuowx9ABAHFCmF+E73yjlyaP7q9Kn19z\n/2ej9hw6Fe+SAABJiDC/SLd/vae+P26gausD+sVrm/Xp9qPxLgkAkGQI8y4wfHAPPX7PYDnshhau\n2KY/rt3PrVMBAJcMYd5Fru2drVn3DYnOFPebtz5TbX1jvMsCACQBwrwLXZXr0b98/+u6piBLm/ce\n1/NLSvTF0ep4lwUASHCEeRdLT3XqyX+4XkXfKFD5qXr97JW/aeX6Lxh2BwDEDGEeA6Zp6LvD++rH\nxV9Tqsuu36/aq5eWblZFFRPMAAC6HmEeQ1/r202zZ9yi6/t1046DJ/Uv//2pVm08pFCIXjoAoOsQ\n5jGWnubU/5n4VX1/3EBJhn7359168XcbVHrMF+/SAAAJgjC/BAzD0PDBPfTiD2/Rzdfk6fPDVZq9\npERvfLRX9f5AvMsDAFgcYX4JZXhceuju6/Tj4sHK8rr03rovNGvhOq3eXKZgKBTv8gAAFkWYx8HX\n+ubohR/coru+2Ut1/oD+75926ae/LdHfdhxVmLPeAQDnyR7vApKVy2HT+GF9NOL6fL398eda89kR\nPb94nb7SM1N33dZbA6/OlGEY8S4TAGABhHmcZXld+sc7rtGYm3pqxV8P6m87juoXr21Sv6sydNc3\ne+naXtmEOgCgXYT5ZeKqPI+e+8Gt+nRLmf64dr+27DuhX/1+i3pfma5xt1ytGwZ0k83kqAgA4EyE\n+WWmT490zSwerINfVuudTw5ow+5yLXh7q3LS3Ro95CoNH3ylUt2OeJcJALiMEOaXqYIrvHrku1/V\nkRM1en/DIa397Ihe/3Cvlq/Zr29cd4WGD+6hgiu88S4TAHAZIMwvc1fmpGnq7V/RhGF99PGWw3p/\nwyF9uKlMH24qU8EVXo0Y3EO3DOquFBe/SgBIViSARXhSHBp3a4Fuv7mnPttXob9sOawt+47rlZW7\ntHTVHt3QP1e3DOqu63pny27j2DoAJBPC3GJspqnr+3fT9f276WR1g9Z8dkRr/35En24/qk+3H1Wa\n266vD8zTLYO6q3/PTJmcCQ8ACY8wt7Asr0vf+UYvFQ0t0P4j1Vq3/Uut33FMH20+rI82H1aW16Ub\n++fq+gHd9JWemfTYASBBEeYJwDAM9emRrj490jV5VH/t+OKkPt1+VBt3leuDjYf0wcZDSnHZ9bW+\nObqhfzd9tU8Ox9gBIIHwFz3BmKaha3tl69pe2Zr27a9od+kpbdpzXJv3lEeH4m2moX75GbquT7au\n7Z2tq7t7GY4HAAsjzBOY3WZqUK9sDeqVrXu/1V+lx3zatOe4tuw9rt2lp7Sr9JTeWv25PCkOXds7\nW4N6ZekrV2cpN8PNrHMAYCGEeZIwDENXd/fq6u5e3X1bb1XX+rX9wElt21+hbQcqor12KXIsfkDP\nzOhPj5xUwh0ALmOEeZLypjp1y6DuumVQd4XDYR0+XqPtB09qT+kp7S491SbcPSkODeiZqX75Gep9\npVcFV3jldvJPBwAuF/xFhgzDUH6uR/m5Ho25qafC4bC+rKjV7qZg31V6Sht3l2vj7vKm7aX8bmnq\nfWW6evdIV58r05Wfm8bc8QAQJ4Q5zmAYhq7MSdOVOWkacX2+JOl4ZZ0+P1yl/UeqtP9wlQ4crdah\n8hp9/PcjkiSH3VSPbmnqmedRzzyPrm56ZB55AIg9whyd0i0jRd0yUnTzNd0lScFQSIeP12r/kSp9\nfrhKB76sUlm5Twe/rG7zvpx0l3rmeXVVnkfX9uumNIep7lmpctjpxQNAVyHMcUFsphnthQ8f3EOS\nFAiG9GVFrUqP+dr8bN57XJv3Htc7nxyQJJmGodxMd1PvP1VX5KSqR9NzevIAcP4Ic3QZu83UVbke\nXZXr0dBrW5ZX1vhVeqxalXUB7TlYoSMnanXkRG1TyLfdR3qaU3mZKcrNdCs3M0V5WSnKy0xVblaK\n0lMdnFUPAGdBmCPmMtKcyuido9xcr8rLW4bhq2v9TcFeEw34Lytq9PnhKu0tqzxjPy6HrU3I56S7\nlZ3uVk66W1npLnlTCHsAyYkwR9x4U53ypjo1oGdmm+WBYEgV1Q0qP1mnY6fqoo/HTtap/FSdDpXX\nnHV/DruCdmsnAAALs0lEQVSpbK9L2eluZae7lO11KyfDrWyvS1npbmV5nEpx2Ql8AAmHMMdlx24z\nlZeZorzMFF172rpwOKzq2kYdO1Wniqp6VVQ16ERVfeR5dYMqqup19ODJc+7bYTcjIwUepzLTXMrw\nOJXhcSmzaVlGmkuZnsiXDNMk9AFYA2EOSzEMQ+lpTqWnOaX8jLNu428M6mRTsJ+oalBFdST0K30N\nOlXjV6WvQfsPVysUrmrncyKT5XhTnfKmOORJbXnuTT39tVPeVE7cAxA/hDkSjtNhU/fsVHXPTj3n\nNqFwWL7aRp3yNaiyxq9Kn1+VNQ065YuEfWWNX9W1jar0Nejw8bMP658u1W1XqsuuNLdDqW670tx2\npaU0Pz/90a5Ut0NpbrtSXHZudAPgohDmSEpm6x5+B4KhkHx1AflqIwFfXdeo6qbnvtpGVddFntf5\ng6r0NehIRY38jaFO12JISnHZo2Gf4rLJ7bTL7bIppenR7bQrxWlTissut9Mmt8uuFKe9ZVunTS6n\njS8FQJIizIEO2Mym4+wdBH/rs/UbAyHV1jeqpj6g2vqAauobT3sMtFrfqJqGyHbn+0WgNUOSq1Xg\nOx02ueymnE6bXI7Ij9Nhk8thtnpuk/O0166mbVq/DoXCF1QTgEsjpmEeDof105/+VLt27ZLT6dS/\n/uu/qmfPntH1q1at0oIFC2S32zVx4kQVFxfHshzgknHYTWV4XMrwuM77vaFQWPX+oOr9AdU1BFTv\nD6rOH1B9w2mP/qDqGwKqO+2x3h+Ur65R/qoG+RuD6qoYdthNOe2m7HZTDpsph91sWmaLPm9ebm/a\ntu3yVtudto/mH7tpymYzZLeZstuanpum7E3LOCkROLuYhvn7778vv9+vpUuXasuWLZozZ44WLFgg\nSQoEApo7d66WLVsml8ulKVOmaPTo0crOzo5lScBlzzSNyPF398X/7xkOh+UPhORvDKqhMaiGxsjz\n0183NP/4g/IHQmpoDMrvD6qh6b0hSTW1fjUGQpGfYEi+usbo6+Al6rkbhuSwmbLZWgLeZjaHv9Gy\n3Gz9uul585eCVl8abGbkx+t1q77OL5sZ+cLQvNw87bHt885ua55zPZdJoqvENMw3bNigYcOGSZIG\nDx6srVu3Rtft27dPBQUF8ng8kqQhQ4aopKRE3/72t2NZEpBUDMOIDpV7L2I/p0/4c7pgKKRAIKzG\nYCTc/YFgS/A3hX9jY9Pj6eua1wdCCgbDCoRCCgQjzxubHgPByLJAKKxgMKRA07Lm7WsbAm2XW+Sw\ngGFEzt8wDEOmGXluGpHAb15nmoZMQ03bNG1rqGm50fT+1q/VZjvDNOR2ORRoDJ6xL7P588943lxH\n07ZNn2EYkqGW2tS8r6bvJM1taVrV8vz0x9brWu3PMCS1+m9iNP03On1bo3U9hqGsU/WqrKyLtq2l\nltb7O8v7m34JzZ+j5roUeV+b9dFlRpttm3+Prbdtfu8Z69X2v2PLdi3rWrZvWd8ZMQ1zn88nr7fl\nT4jdblcoFJJpmmesS0tLU3X1uf9YALh82UxTNqfkki3epUiKjEgEQ+HoqEHzl4GWLwaRLwGhUFih\nUFger1sVJ2sVDIUVDIYVCocVbFofDIYVDEe2a7u+eVmo1fNwm+dnLAuGovsKhcIKhMIKh8MKhSJX\nWIRDkX2Hw5HXoabXoZCi7QiFI4diwuGWda2fh8LW+CKDzrGZht7+xV0dbhfTMPd4PKqpabmspznI\nm9f5fL7oupqaGqWnp3e4z9zci+lfXP5on7UlcvsSuW2A1cX0PpQ33nijVq9eLUnavHmzBgwYEF3X\nt29fHTx4UFVVVfL7/SopKdH1118fy3IAAEhIRjgcuzGZ1mezS9KcOXO0bds21dXVqbi4WB999JF+\n85vfKBwOa9KkSZoyZUqsSgEAIGHFNMwBAEDsxXSYHQAAxB5hDgCAxRHmAABYHGEOAIDFWeJGKx3N\n8Z4ItmzZol/+8pd69dVX411KlwoEAnr66adVVlamxsZGPfTQQxo1alS8y+oyoVBIzz77rPbv3y/T\nNPX888+rX79+8S6ry504cUITJ07Uyy+/rN69e8e7nC713e9+NzoT5VVXXaUXX3wxzhV1rUWLFmnV\nqlVqbGzUvffeq4kTJ8a7pC7zhz/8QcuWLZNhGGpoaNDOnTu1du3a6O/TygKBgJ566imVlZXJbrfr\nhRdeaPf/PUuEeXtzvCeCxYsXa/ny5UpLS4t3KV1uxYoVysrK0rx581RZWanx48cnVJivWrVKhmHo\ntdde0/r16/WrX/0qof5tSpE/Ks8995zcbne8S+lyfr9fkvTKK6/EuZLYWL9+vTZt2qSlS5eqtrZW\nv/3tb+NdUpeaMGGCJkyYIEmaPXu2Jk2alBBBLkmrV69WKBTS0qVL9cknn+jXv/61/v3f//2c21ti\nmL29Od4TQUFBgebPnx/vMmJi3LhxmjlzpqRIL9Zut8T3x0771re+pRdeeEGSVFZWpoyMjDhX1PV+\n/vOfa8qUKcrLy4t3KV1u586dqq2t1YwZM/T9739fW7ZsiXdJXWrNmjUaMGCAfvSjH+nhhx/WyJEj\n411STHz22Wfau3dvQt15s1evXgoGgwqHw6qurpbD4Wh3e0v8ZW1vjvdEMGbMGJWVlcW7jJhISUmR\nFPkdzpw5U48//nicK+p6pmnqn//5n/X++++3+83ZipYtW6acnBx985vf1H/+53/Gu5wu53a7NWPG\nDBUXF+vAgQP64Q9/qJUrVybM35aTJ0/q8OHDWrhwoUpLS/Xwww/rT3/6U7zL6nKLFi3So48+Gu8y\nulRaWpoOHTqksWPH6tSpU1q4cGG721viX2x7c7zj8nfkyBFNnz5dEyZM0B133BHvcmJi7ty5Wrly\npZ599lnV19fHu5wus2zZMq1du1ZTp07Vzp079dRTT+nEiRPxLqvL9OrVS3fddVf0eWZmpsrLy+Nc\nVdfJzMzUsGHDZLfb1bt3b7lcLlVUVMS7rC5VXV2tAwcO6Oabb453KV1qyZIlGjZsmFauXKkVK1bo\nqaeeih4WOhtLJGJ7c7wnkkScjO/48eOaMWOG/umf/il6bCuRLF++XIsWLZIkuVwumaaZUF80f/e7\n3+nVV1/Vq6++qoEDB+rnP/+5cnJy4l1Wl3nrrbc0d+5cSdLRo0dVU1Oj3NzcOFfVdYYMGaKPP/5Y\nUqR99fX1ysrKinNVXaukpES33nprvMvochkZGdHj/16vV4FAQKFQ6JzbW2KYfcyYMVq7dq0mT54s\nKTLHeyJqvvdtIlm4cKGqqqq0YMECzZ8/X4ZhaPHixXI6nfEurUvcfvvtmjVrlu677z4FAgE988wz\nCdO20yXiv89JkyZp1qxZuvfee2Wapl588cWE+jJWWFiov/3tb5o0aZLC4bCee+65hPs97t+/P+Gu\nbpKk6dOn6+mnn9b3vvc9BQIBPfnkk+2ehMrc7AAAWFzifAUFACBJEeYAAFgcYQ4AgMUR5gAAWBxh\nDgCAxRHmAABYHGEOAIDFWWLSGACXh3379unzzz/X1q1b9cADDyTknf4AK6JnDqDT/vd//1cDBw7U\n4MGDo1MsA4g/ZoADcN5ef/119evXTzfeeGO8SwEgeuYAzlN9fb1qamoIcuAyQpgDOC9r1qzR9OnT\nVVpaGu9SADRhmB1AGwcPHlRJSYn++te/6qWXXpIUOfHN5/Np3759euONN5SSkqIHH3xQt9xyS5yr\nBSAR5gDO4ZFHHtEPfvAD5efna8eOHRoxYkS8SwJwDoQ5gLNatWqV3n//fQ0fPlxjx46NdzkA2kGY\nAzir6upqjRs3TmvWrIl3KQA6wAlwAM7qgw8+0HXXXafdu3fHuxQAHSDMAZzh3Xff1e23366ioiK9\n+eab8S4HQAcIcwBtrFmzRjfffLNSU1N1++23q6ysTMuXL1dFRUW8SwNwDhwzBwDA4uiZAwBgcYQ5\nAAAWR5gDAGBxhDkAABZHmAMAYHGEOQAAFkeYAwBgcYQ5AAAWR5gDAGBx/x+6K9mAPmFuaQAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11add6e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0,8,300)\n",
    "\n",
    "chi2_df1_pdf = stats.chi2.pdf(x, df=1) \n",
    "chi2_df2_pdf = stats.chi2.pdf(x, df=2)\n",
    "chi2_df3_pdf = stats.chi2.pdf(x, df=3)\n",
    "chi2_df4_pdf = stats.chi2.pdf(x, df=4)\n",
    "chi2_df5_pdf = stats.chi2.pdf(x, df=5)\n",
    "\n",
    "plt.plot(x, chi2_df1_pdf, label=\"df = 1\")\n",
    "\n",
    "\n",
    "plt.xlabel('$\\chi^2$')\n",
    "plt.ylabel('Density')\n",
    "plt.legend()\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our bumpus example, the degrees of freedom is df=1. We can calculate the corresponding p-value for our sample chi2 value of ~3.13."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.077031938117561066"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.chi2.sf(chi2, df=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $X^2$-square test of independence using scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Same analysis with scipy.stats fxns\n",
    "observed_nomargin = pd.crosstab(bumpus.survived, bumpus.sex, margins=False)\n",
    "Chi2, Pval, Dof, Expected = stats.chi2_contingency(observed_nomargin.values, \n",
    "                                                   correction=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.1264367816091938, 0.077031938117561066, 1)"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Chi2, Pval, Dof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}