notebooks/paired t-test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Paired t-test   \n",
    "```bash\n",
    "by Jeffrey Chijioke-Uche, MSIS, MSIT, CPDS (IBM Sr. Solution Architect, Hybrid Cloud & Multicloud)\n",
    "IBM Information Technology PhD Scholar at Harvard University & Walden University \n",
    "``` \n",
    "\n",
    "\n",
    "## Contents:    \n",
    "1. [Theory](#theory)   \n",
    "2. [Numerical example](#example)\n",
    "3. [Scipy's example](#scipy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='theory'></a>\n",
    "\n",
    ">Paired t-test is a statistical procedure used to determine whether the mean difference between two sets of observations is zero. In a paired sample t-test, each subject or entity is measured twice, resulting in pairs of observations. Source: [statisticssolutions](http://www.statisticssolutions.com/manova-analysis-paired-sample-t-test/)\n",
    "\n",
    "Paired t-test, also called dependent t-test, is used to \n",
    "* Determine if a given medicine has an effect on patients, in which samples are collected before and after medication for the same people, and we can evaluate the statistical difference between samples    \n",
    "\n",
    "\n",
    "### Assumptions:\n",
    "1. The dependent variable must be continuous (interval/ratio)  \n",
    "2. The observations are independent of one another.   \n",
    "3. The dependent variable should be normally or approximately normally distributed.   \n",
    "4. The dependent variable should not contain outliers.   \n",
    "\n",
    "### Hypothesis\n",
    "|Set|\tNull hypothesis\t| Alternative hypothesis |No. of tails |\n",
    "|---|:-----------------:|:----------------------:|:-----------:|\n",
    "|1\t| $\\mu_1 - \\mu_2 = 0$ | $\\mu_1 - \\mu_2 \\neq 0$ |\t2|\n",
    "|2\t| $\\mu_1 - \\mu_2 > 0$ | $\\mu_1 - \\mu_2 < 0$ | 1|\n",
    "|3\t| $\\mu_1 - \\mu_2 < 0$ | $\\mu_1 - \\mu_2 > 0$ | 1|\n",
    "\n",
    "\n",
    "### Test method\n",
    "\n",
    "We will use the ** paired t-test ** method.   \n",
    "\n",
    "#### Steps:   \n",
    "\n",
    "We will use the following nomenclature:   \n",
    "\n",
    "$D =$   differences between two paired samples   \n",
    "$d_i =$ $i^{\\textrm{th}}$ observation in $D$     \n",
    "$n =$   sample size    \n",
    "$\\overline{d} =$ sample mean of the differences    \n",
    "$s_d =$ sample standard deviation of the differences   \n",
    "\n",
    "\n",
    "#### Standard error (SE):   \n",
    "SE$ =  \\frac{\\displaystyle s_d}{ \\sqrt{n}} $\n",
    "\n",
    "\n",
    "#### Degrees of freedom (DF):   \n",
    "DF $ = n-1$\n",
    "\n",
    "\n",
    "#### Test statistic:   \n",
    "$t = \\frac{\\displaystyle \\overline{d} }{\\displaystyle \\textrm{SE}}$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='example'></a>\n",
    "## Numerical example    \n",
    "In this test, let's use a hypothethical example of blood glucose (BG) readings of persons with diabetes.    \n",
    "The test is done twice, say before and after a given clinical treatment and we want to know if there is a significant difference after the treatment.     \n",
    "\n",
    "The database is not real but just to be used as a toy example.   \n",
    "This contains BG readings (in mg/dL) for 40 persons with diabetes at different times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BG 1</th>\n",
       "      <th>BG 2</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Person ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>89.162573</td>\n",
       "      <td>94.023517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>90.857629</td>\n",
       "      <td>95.273755</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>94.912999</td>\n",
       "      <td>96.612870</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>95.254064</td>\n",
       "      <td>96.818673</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>97.136291</td>\n",
       "      <td>97.760342</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                BG 1       BG 2\n",
       "Person ID                      \n",
       "0          89.162573  94.023517\n",
       "1          90.857629  95.273755\n",
       "2          94.912999  96.612870\n",
       "3          95.254064  96.818673\n",
       "4          97.136291  97.760342"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "# This is to print in markdown style\n",
    "from IPython.display import Markdown, display\n",
    "def printmd(string):\n",
    "    display(Markdown(string))\n",
    "\n",
    "# Read dataset locally:\n",
    "#df = pd.read_csv('../data/BG-db.csv',index_col=0)\n",
    "\n",
    "\n",
    "# Read dataset from url:\n",
    "import io\n",
    "import requests\n",
    "url=\"https://github.com/schijioke-uche/data-analysis-with-python-an-spss-model/blob/3b9249b72e078eca78a0c38c5aa215a67f21656f/data/BG-db.csv\"\n",
    "s=requests.get(url).content\n",
    "df=pd.read_csv(io.StringIO(s.decode('utf-8')),index_col=0)\n",
    "\n",
    "\n",
    "del df['BG 3']\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAFjCAYAAACgxwiQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VOXZ8PHfmS3rZN83IAk7ylpABRQUtdr6KKgJILbqq6++hdYdrYo8LoC1PlbBpbb1UXEp7lVbqwJWrAsFZEtCEhKSkH3fJsus5/0jMBBJmBAyM0nm+n4+fEzmnDnnGjCZa+77uu9LUVVVRQghhBA+S+PtAIQQQgjhXZIMCCGEED5OkgEhhBDCx0kyIIQQQvg4SQaEEEIIHyfJgBBCCOHjJBkQQgghfJwkA0IIIYSPk2RACCGE8HGSDAghhBA+TpIBIYQQwsfpvB2AO3R2dpKVlUV0dDRardbb4QghhBBuZbfbqa2tZdKkSfj7+5/284dlMpCVlcWyZcu8HYYQQgjhUW+88QYzZsw47ecNy2QgOjoa6PpLiYuL83I0QgghhHtVVVWxbNky5/vf6RqWycCxqYG4uDiSkpK8HI0QQgjhGf2dGpcCQiGEEMLHSTIghBBC+DhJBoQQQggfJ8mAEEII4eMkGRBCCCF8nCQDQgghhI+TZEAIIYTwcZIMCCGEED7ObZsOORwO1qxZQ15eHgaDgccee4wRI0Y4j2/bto3nnnsOnU7H4sWLufbaawG46qqrCA4OBiApKYl169a5K0QhhBBC4MZkYMuWLVgsFjZv3szevXtZv349L7zwAgBWq5V169bx7rvvEhAQwJIlS1iwYAFGoxFVVdm0aZO7whJCCCHEj7htmmD37t3MnTsXgClTppCVleU8VlhYSEpKCqGhoRgMBqZPn87OnTvJzc2lo6ODG2+8keuvv569e/e6KzwhhBBCHOW2kQGTyeQc7oeu/ZJtNhs6nQ6TyYTRaHQeCwoKwmQy4e/vz0033cQ111xDcXExN998M//85z/R6XoPc8OGDWzcuNFdL0MIn7Gl8GuP3euitLkeu5cQwjW3JQPBwcG0tbU5v3c4HM439R8fa2trw2g0MmrUKEaMGIGiKIwaNYqwsDBqa2uJj4/v9T4rV65k5cqV3R4rKyvjwgsvHOBXJIQQQgxPbpsmmDZtGtu3bwdg7969jBkzxnksLS2NkpISmpqasFgs7Nq1i6lTp/Luu++yfv16AKqrqzGZTP1uxyiEEEKIvnHbyMDChQv55ptvyMzMRFVV1q5dy8cff0x7ezsZGRncd9993HTTTaiqyuLFi4mNjeXqq6/m/vvvZ8mSJSiKwtq1a085RSCEEEKIM+e2d1qNRsMjjzzS7bG0tDTn1wsWLGDBggXdjhsMBp566il3hSSEEEKIHsjHbiEEAAeyOwf0emdN9D/l8R07dnD77beTnp6OqqpYLBbWrFnDhAkTen3OnXfeSUlJCb/73e+6fbjor82bN7No0SL0en2/r5GXl0dLSws/+clPzjieE5nNZj766COuueaaAb2uED2RHQiFEF4ze/ZsNm3axOuvv86vf/1rnnnmmVOe/+233/Lee+8NSCIA8Mc//hGHw3FG1/j8888pKCgYkHhOVFtbyzvvvDPg1xWiJzIyIIQYFFpaWoiIiAC6Pm0/9thjAISFhbF27VqeeuopTCYTt912G88++yz3338/ZWVl2O12brjhBi677DKWL19OREQEzc3NvPTSS6xZs4aSkhIcDge33347s2bNct7vnXfeoba2ljvuuINf/OIX/P73v0ev13PttdeSkJDA008/jVarJTk5mUceeQSz2cwDDzxAa2srNTU1LF26lAsvvJAPPvgAvV7PxIkTeeCBB5gxYwZ5eXmkpqYSGRnJrl27MBgMvPTSS3R2dvLAAw/Q2NgIwIMPPsjYsWO5+OKLmTZtGkVFRURGRrJhwwZefPFFCgoK2LhxIytWrPD8P4jwKZIMCCG85vvvv2f58uVYLBZyc3N57rnnAHjooYdYu3Yt6enpvPPOO/z5z39mzZo1fPHFF7zwwgu8/vrrRERE8Pvf/x6TycSiRYuYPXs2AD/72c9YuHAhb775JuHh4axdu5bGxkauu+46/v73vzvvfc011/DCCy/w9NNPs3fvXsxmM++88w6qqnLppZfy5ptvEhkZyR/+8Ac++OADJk6cyOWXX87FF19MdXU1y5cvZ+nSpVx11VVERUVx9tln09bWxs9+9jMefvhhLr30Uu6//37uuOMOrrvuOgoKCvjkk0+YPXs2S5cupbi4mPvvv5+33nqL0tJSXn31VeLj48nMzOTAgQPceuut5OfnSyIgPEKSASGE18yePZunn34agMOHD5OZmcn27dspLCzkv//7v4Gu7ctHjhzZ7XmFhYWce+65QNe+JWlpaZSWlgIwatQoAPLz89m9ezf79+8HwGaz0dDQ4Bx9+LFjz2toaKCmpobbb78dgM7OTs4991zOP/98Xn31VT7//HOCg4Ox2Ww9XmfixIkAhISEOKczQkJCMJvN5Ofn8/333/Ppp58C0NzcDEB4eLhzP5X4+HjMZnOf/w6FGAiSDAghBoWoqCjn16NGjeKJJ54gISGB3bt3U1tb2+3ctLQ0du3axcKFCzGZTOTn55OUlASAoigApKamEhcXx6233kpnZycvvPACYWFh3a6jKIqzZkCj6SqhCg8PJy4ujueffx6j0cjWrVsJDAzk5ZdfZsqUKSxdupTvv/+er7766qRrnHj/nqSmpnLFFVfw85//nPr6emdNQE/P0Wg0Z1zPIERfSTIghPCaY9MEGo2GtrY27rvvPvz9/VmzZg2rVq3CZrOhKAqPP/54t+dde+21PPTQQyxZsgSz2cyKFSuIjIzsdk5mZiYPPvgg1113HSaTiaVLlzrf8I+ZMWMGt9xyC7/61a+cj2k0Gh544AFuueUWVFUlKCiI3/3udyiKwmOPPcY//vEPjEYjWq0Wi8XCpEmT+ry64dZbb+WBBx7g7bffxmQynXIKIDIyEqvVypNPPsk999zTl79OIfpNUVVV9XYQA+3YdsRbt251floQQpya9CYQYug60/c9WVoohBBC+DhJBoQQQggfJ8mAEEII4eMkGRBCCCF8nCQDQgghhI+TZEAIIYTwcbLPgBACgIO1A9tsZ3x0ep/P/dOf/sSrr77K1q1b8fPzY9++fdx9991ceumlzJs3D6PRyLhx407r/gPV9e+LL77g7LPPJjY29oyu82MVFRXk5uae1MpdCG+QkQEhhNd99NFHXHbZZc7eAV9//TXXX389d911F++99x41NTWnfc2B6vr32muvYTKZzvg6P/b999/zww8/DPh1hegPGRkQQnjVjh07SElJITMzk3vuuYf09HTef/999Ho9VquVr7/+muzsbNLT09m3bx+vvPIKGo2G6dOnc/fdd7Nhwwb27NlDe3s7jz/+uHMnwBO7/qmq2u2cb7/9lk8++QRFUbjsssu4/vrryc/PZ/369djtdhobG1mzZg0tLS0cPHiQVatW8eSTT3LvvfcSHx9PWVkZl19+OYcOHSInJ4cLLriAO++8s8duizk5OfzpT39Cr9dTVlbGZZddxi233OLsYjh16lQuvPBCb/4TCCHJgBDCu9555x2uueYaUlNTMRgMqKrq7AS4ZMkS8vPzueyyywgMDGTDhg289957BAQEcM899/DNN98AXXv+P/jgg92ue2LXvw0bNjjPKSgo4B//+AdvvvkmADfccANz5syhoKCAVatWMXbsWD7++GPef/99HnvsMcaPH8+aNWvQ6/WUlpby8ssv09nZyYUXXsj27dsJCAhg/vz53HnnnT12Wzz33HOpqKjgo48+wmKxMHfuXG677TZuueUWDh8+LImAGBQkGRBCeE1zczPbt2+noaGBTZs2YTKZeP3110lJSTnp3CNHjtDQ0MAtt9wCQFtbG0eOHAGOdxw8lRO7GVZUVPDLX/7SGUNJSQkxMTE8//zz+Pv709bWRnBw8EnXSE5Oxmg0YjAYiIqKcjY+OtZoqLdui2PGjEGn06HT6fD39z+NvyEhPEOSASGEx/3zu2IAtn/xN6afs5D/yrwZAIu5k/++5xdMmTGX+CSFf35XTEVdGztzqhgxKoSgkEgyb3kYrU7Hjq8/p02JoaC0AGOL4rzmMQ31NTS2dDi/P9akKDU1lfT0dP785z+jKAqvvPIKY8eO5Ve/+hW///3vSUtL49lnn6W8vBzoeqM/1sLlVB0Jofdui9KVUAx2UkAohPCa77Z/yozzjg+TG/z8mTxjDt999anzsRGp4/j4nZcxmVq44JJFPLv+bv7nkV9z8MBOYuISe7220RiG3WbjySef7Pb4uHHjOOecc1iyZAmLFi2iuLiY2NhYrrjiCn7zm9+wdOlSiouLnUWLU6dO5d5776W5udnl6znWbXHJkiU89dRTjB07ttdzx4wZw9atW51Fk0J4k3QtFEIAnu1aaKtJ9ti9Lj1npMfuJYS3nOn7nkwTCCH6zWa3UdZShdluQavRkBySgJ/O4O2whBCnSZIBIcRps9qt7CzfT15dIVaHzfm4VtGQHjmKc5KnYdDqvRihEOJ0SDIghDgtNW31bD38Da1mE0GGQCZGjCHU30i7tZP8usPk1RVSY6rjkvR5hPgbvR2uEKIPJBkQQvRZeUsVnxVsx+6wMyVuAtMTzkKr0TqPT44dz/dlP5BVk88n+dtYNP4S/PWylE6IwU5WEwgh+qS8pYpPD/0Lh+pgYdpcZiZN6ZYIQNdyuXNTZjA94SxMlja2HP4GhyrL54QY7CQZEEK41NDRxOdHVxtcmn4+I8NPXa08LX4SI8OSqGitZm9ltidCFEKcAUkGhBCn1Gkz889DX2G1W7lg5GySQuNdPkdRFC4YOZsAvT97KnNoNQ98ox8hxMCRZEAI0StVVfmy6DtMljamxU8iPXJkn59r0BmYnTQVu2rnu9I97gtSCHHGJBkQQvRqX9VBSpsrSAqJZ3rCWaf9/PSIkcQFR1PcVEpVa60bIhRCDARJBoQQPapra2BnxT4C9QHMH3WOy335e6IoCjMTJwOwt0pqB4QYrCQZEEKcxOawsa3oW1RVdc7991ecMYa44GiONFdQ3944gFEKIQaKJANCiJPsKt9PU2cLE2PG9Klg0JUpcRMA2FuVc8bXEkIMPEkGhBDd1LTVc6A6jxC/YGYlThmQayaHJhAREMbhxiO0WztcP0EI4VGSDAghnOwOO9uLd6CiMm/ELHTagdmkVFEUxkeno6oq+XWHB+SaQoiBI8mAEMJpX9VBGjqaGBeVRkJI7IBeOz1iJFqNlty6QoZh53QhhjRJBoQQADR2NPNDZRaB+gBmJU0d8Ov76QykhafQYjZRaS4Z8OsLIfpPkgEhBA7VwfaS/+BQHcxJ+Ql+OoNb7jMuKh2AfNMBt1xfCNE/kgwIIdh2+FuqTbWMCk922XfgTMQGRxHiF8yRjkPYHFa33UcIcXokGRDCxzV3tvDG/g/Qa3ScmzzdrfdSFIW0iBHYVCtHOgrcei8hRN9JMiCEj9u0933aLO38JHEyQYZAt98vLWIEAIfbD7r9XkKIvpFkQAgfllWdx/aSHaSGpzAhZrRH7hkREEa4PoqyjsOYHZ0euacQ4tQkGRDCR1ntVv60+00UReGWGUvRKJ77dTAqcDwO7Bxpl6kCIQYDSQaE8FF/y/2CytYaLkk/n9SjQ/eeMjJwLABHOg559L5CiJ5JMiCED6psreGDnE8J9w8l86wrPH7/MH0koboIyjuLZFWBEIOAJANC+BhVVfnL7r9iddj45bRrCNQHeCWOlMDR2FQrFZ2yAZEQ3ibJgBA+5psju9hffZCp8ROZnTTNa3GMCOgqWJSpAiG8T5IBIXxIm6WdV/e+i16r58ZpGSiK4rVYog0JBGiCONJRgEN1eC0OIYQbkwGHw8Hq1avJyMhg+fLllJR0Hwrctm0bixcvJiMjg7fffrvbsfr6es4//3wKCwvdFZ4QPumtA3+jubOFxRN+SmxwtFdjURSF5IA0Oh3t1FuqvBqLEL7ObcnAli1bsFgsbN68mbvuuov169c7j1mtVtatW8fLL7/Mpk2b2Lx5M3V1dc5jq1evxt/f312hCeGTCuqL+aLgaxJD4rhi7EJvhwNAUkAqAGWdRV6ORAjf5rZkYPfu3cydOxeAKVOmkJWV5TxWWFhISkoKoaGhGAwGpk+fzs6dOwF44oknyMzMJCYmxl2hCeFz7A47L+16AxWVm6cvRafVeTskAOL9R6CgUN4hyYAQ3uS23wgmk4ng4GDn91qtFpvNhk6nw2QyYTQanceCgoIwmUy8//77REREMHfuXF566aU+3WfDhg1s3LhxwOMXYjj556F/UdxUxgWjzvHYToN94afxJ9ovgVpzBWZHJ34aGREUwhvcNjIQHBxMW1ub83uHw4FOp+vxWFtbG0ajkffee49vv/2W5cuXc/DgQVatWkVtbe0p77Ny5Ury8vK6/dm6dat7XpQQQ1BTRzNvZ39CkCGQ6yYv8nY4J0nyT0VFpaKz2NuhCOGz3JYMTJs2je3btwOwd+9exowZ4zyWlpZGSUkJTU1NWCwWdu3axdSpU3njjTd4/fXX2bRpE+PHj+eJJ54gOtq7RU5CDHVv7P+QDmsnS866ghC/YNdP8LBE/1EAMlUghBe5bZpg4cKFfPPNN2RmZqKqKmvXruXjjz+mvb2djIwM7rvvPm666SZUVWXx4sXExsa6KxQhfFZeXSFfFX/PyLAkLkqd6+1wehRliMNfE0B5ZxGqqnp1uaMQvsptyYBGo+GRRx7p9lhaWprz6wULFrBgwYJen79p0yZ3hSaET3A4HLz8w2YAbpyWiUYzOLcVURSFBP+RHG4/SJO1jnCDjAYK4WmD87eDEOKMbT38DUWNpcwbMYtx0Wmun+BFif5dSwzLZYmhEF4hyYAQw1Cr2cRbB/5GgM6fZZOv8nY4LiUGjARkvwEhvEWSASGGob8e+AiTpY1rJl1OeECot8NxKVAbTIQ+hurOUuliKIQXSDIgxDBzuOEIWwr/TWJIHJeOnu/tcPosMWAUduxUmo94OxQhfI4kA0IMI6qq8sqet1FRuXFaBjqN1tsh9VnS0SWGst+AEJ4nyYAQw8iuiv3k1hUyI+Fszood5+1wTkuMXyJaRUdlp4wMCOFpkgwIMUzYHXbe3PchGkUzJIoGf0yr6IgxJNJgraHT3u7tcITwKZIMCDFMfFn0LeWtVSwYdS6JIXHeDqdfEvxHAEjdgBAeJsmAEMNAp83M21mf4Kc1cM2kn3k7nH6L908BkKkCITxscPQxFUL06A9ffNin80o7CmnqbCHJP5VX//1lv+41KTWyX88bSFGGOHSKgcrOEm+HIoRPkZEBIYY4q8NCeUcROsVAYsAob4dzRjSKlji/JJptDbTbWr0djhA+Q5IBIYa4is5iHNhJDkhFqwz9wb54qRsQwuMkGRBiCLM6LFR2HkGv+BHrl+TtcAbEsbqBCpkqEMJjJBkQYgg7NiqQGDAKjTJ0Nhg6lQh9DAaNvxQRCuFBkgwIMUQdHxUwDJtRAQCNoiHeLwWTvZlWW5O3wxHCJ0gyIMQQVdlZ0jUq4D8K7TAZFThGlhgK4Vl9SgZqamoA2LVrF2+88Qbt7bI7mBDeZHNYqTQfHRXwT/Z2OAMu3u9YMiB1A0J4gstk4OGHH+aFF16goKCAu+66i+zsbFatWuWJ2IQQvag0l2BXbST4jxx2owIAYfoo/DWBVJqPoKqqt8MRYthzmQwcOHCA1atX8+mnn3L11Vezdu1aKioqPBGbEKIHDtVOZWcpOkVH3DAcFQBQFIV4/xTa7SZabA3eDkeIYc9lMmC323E4HGzdupV58+bR0dFBR0eHJ2ITQvSg1lKJTbUQ65c8LPYV6M2x/QYqpG5ACLdzmQxceeWVzJkzh8TERCZPnsyiRYvIyMjwRGxCiB9RVfXoPLpC3NEiu+Eqwe/Y5kNSNyCEu7n8WHHDDTdw/fXXo9V2zUu+8cYbREREuD0wIcTJmm0NtNtNRBni8NP4ezsctzLqwgjSGqnq7KobUBTF2yEJMWz1mgwsX778lD98r732mlsCEkL0rqKzGDg+hD6cddUNjKCgLYtGay0RhhhvhyTEsNVrMrBy5UpPxiGEcKHD3kaTtQ6jLgyjLszb4XhEvF8KBW1ZVHaWSDIghBv1WjMwc+ZM5x+tVkthYSFTpkxBURRmzpzpyRiFEBzfq98XRgWOcRYRStMiIdzKZQHhq6++yh/+8AdeeeUV2traWL16NX/5y188EZsQ4iibw0qtuQI/jT+Ret/5hBysCyFEF05V5xEcqsPb4QgxbLlMBj744AP+8pe/EBAQQHh4OO+++y7vvfeeJ2ITQhxVa6nEgZ1Yv2QUxbd2EY/3T8GqWqi3VHs7FCGGLZe/VTQaDQaDwfm9n5+fc2WBEMIzasxlgEKMX6K3Q/E459bEssRQCLdxmQzMnDmTJ554go6ODrZs2cJtt93G7NmzPRGbEAIw2Vpos7cSoY/GoPHzdjged6xuQPoUCOE+LpOBe++9lxEjRjB27Fg+/PBDzj//fOlNIIQHVZvLAIgZJm2KbTaVJpOZhpZO6po6aGmznLL/QIA2iDB9FFXmMuyq3YORCuE7XG46tG7dOq644goyMzM9EY8Q4gR21UadpRKD4ke4PtLb4Zw2u12lqtpGeaWNykob9Y12WlsdQGO383RahTCjPyPijIwbEUGYsfsISILfCHKsddSaK4ZtPwYhvMllMjBy5EjWrl1Lc3MzP/vZz7jiiitIShoen1CEGOzqLdXYVRvx/ilDpnDQalUpKbVSeNjC4WIrFsvxT/1BQQqJCTqMeiM6rQZFgU6LncbWThqPjhTszq0hPjKI8yYnEBsRCHQVEeaYdlNpLpFkQAg3cJkMLFu2jGXLllFRUcGnn37Kr371KwIDA3nrrbc8EZ8QPu34FMHgLhw0mx0UFVspOGyhpNSKzdb1uDFYw/ixBhIT9CTE6wgK7EpobDUnv6Hb7A4OlzeTW9xAaY2Jd7cd4qy0SGZNinf2YajsPMLUUI+9LCF8Rp9anrW2tvLtt9/yzTffYLfbmTNnjrvjEmJQ21L4tdvv0dTRTKutiVBdJP7aQLff73S1tTs4XGSh8LCV0nIrjqPbAISHaUhPNZCWaiAmWtvjtua5pr09XzQCxkRAdAPkZsOBwnoKKuuYMgOCtCFUm8vIbt2NVun7iqZLGdmPVyeEb3GZDNx6663k5ORw8cUX85vf/IbJkyd7Ii4hfF5efREAsYNoVKCl1U7h4a4RgIpKm/PxmGgtaaMMpKcaiIg486XH4REw6zwoPKRypEhh1/cqcVMjaKOFVlsTYUOwfkKIwcxlMnDttdcyb948dLrh2zddiMFGVVUKG4rRoiXcy3vyt7U5yC+wkHfITHXN8Wr+hHgdaakG0kfpCQkZ+L1HNBoYPRYCAlTycqAyPwL96GKarfWSDAgxwFy+w0+bNo01a9Zw5MgRnnnmGX73u99x3333ERoqE3dCuEu1qQ6TpZ1oQ8JpDYkPFLPZQcFhK3mHzJSV21BVUBRISdKRnmYgdZTBOf/vbkkpYDBA1oFwdKpCo7mBEYNv1kSIIc1lMrB69WrOO+889u/fT1BQEDExMdxzzz289NJLnohPCJ9U0FAMQJQhzqP3ra2zsT/LTG6+2VkEGB+rY+wYA+lpnksAfiwmDsZYdBSZQmkPbqLTYsXfoPdKLEIMRy6TgbKyMjIyMnjrrbcwGAzccccdXHHFFZ6ITQif5FAdFDWW4q/zI9RDw+HNTfDhx62UlFoBCDFqmDjej7FjDIS6YQqgP5JSoPZIBB1KE9kFjUwbH0MPtYlCiH5wmQxotVpaW1udFcHFxcVoNENjvbMQQ1FFSzUdtk4mRI9GY3fvz1pHO+TnQl2NAlhJStQxdbI/I1P0aDSD7512VFwEOabDtNFA2ZEYkn2nm7MQbuUyGVi5ciXLly+nsrKS//f//h979+5l7dq1nohNCJ9U0NC1B396xAjqat1zD4cDSkvg8CFwOBRCw1QuOj+EpMTBPfQeog9DQYM2pIGCgxAZBYFB3o5KiKHPZTIwb948Jk2axP79+7Hb7TzyyCNERUV5IjYhfI7NYaeoqZQgQyCxwdHU1TYM+D06O+DAXmhpVtAbVMaPU4mNZ9AnAgAaRUuILozmwAYcipWcLD3TZyLTBUKcoV6TgY0bN/b4+MGDBwFYsWKFeyISwoeVNldgtVsZH5Xe42Y9Z6qxHg7sA6tFITZeZex40BtcP28wCdVH0GxrIDy5nsaieEpLVFJGejsqIYY2mfwXYhApPDZFEDnwk+EVZbBnF9isMGa8ysSzh14iABCq6yqqNMY3oNOrFBWAxeLloIQY4nodGZBP/kJ4ls1uo6S5nFA/I5EB4QN67bIjkJejoNernD0Nwgb28h4VpAtBg5ZWRz2paZCfq1BUoDJ2grcjE2LokpEBIQaJspZK7A47o8KTB3SKoLT4aCJgUJk2c2gnAgAaRUOoPoJORztRSR0EBKqUl0KbyduRCTF0STIgxCBR3NTVoXBk2MC1CK+q6PrkbPBTmT4Tgo0DdmmvOrYdcYu9ntFjQVUVCvK8HJQQQ5gkA0IMAg7VQUlzBYH6AKKDBmajoaZGyDkAWp3K1J9AUPCAXHZQOLYZU5O1nqgYCAtXqatVaBz4xRdC+IRekwGr1cqTTz7JO++8A8CcOXMYN24cEydOJDs722MBCuELqky1mG1mRoQlDsgUQUc77P+h6+uzpkDwMEoEAAI0QRg0/jRb6wGV9LFdjxcXejUsIYasXpOB//mf/6G6upqLLroIgKioKHJzc3nmmWf405/+5PLCDoeD1atXk5GRwfLlyykpKel2fNu2bSxevJiMjAzefvttAOx2O/fffz+ZmZksWbKE/Pz8M3ltQgwZxY0DN0XgcMD+PWC1KowZ37Uxz3CjKAphukhsqpU2ewuhYRARqdJQr9Dc5O3ohBh6ek0Gtm3bxvr16wkP715tdNFFF3Ho0CGXF96yZQsWi4XNmzdz1113sX79eucxq9XKunXrePnll9m0aRObN29RfaLqAAAgAElEQVSmrq6OL7/8EoC//vWv3H777Tz99NP9fV1CDBmqqlLSVIZeqyfBGHvG1yvMB1OrQkKSSlLKAAQ4SIWdMFUAMDKt6/EiGR0Q4rT1urRQr9ej0x0//PDDDzu/NhhcL07evXs3c+fOBWDKlClkZWU5jxUWFpKSkuJsgzx9+nR27tzJT3/6Uy644AIAKioqCAkJOb1XI8QQ1NDRRKuljbSIEWg1Z9YUqL4OjhQrBAaqjBk3QAEOUsfrBupICkglPKKrdqC+VqG1RcUovz6E6LNekwGtVktdXZ1z6+GpU6cCUF1djVbr+heWyWQi+ISJSq1Wi81mQ6fTYTKZMBqPlzUHBQVhMnWtC9LpdKxatYovvviCZ5991uV9NmzY0OtuiUIMBcdWEYwISzyj61gsXQWDiqIycTJoXW42PrTpNQaCtCG02pqwqza0io6RabB3V9fowNlTvR2hEENHr9MEGRkZrFixgsOHDzsfKykp4Y477mDZsmUuLxwcHExbW5vze4fD4Rxp+PGxtra2bsnBE088wWeffcZDDz1Ee3v7Ke+zcuVK8vLyuv3ZunWry/iEGCxKmsrQKBpSQhLO6DoFuWAxK6SOhpDQAQpukAvTR6Ki0mJtBCAiEowhKrXV4OJXhxDiBL0mA0uXLuWCCy7g2muv5ZxzzuHcc89l8eLFXHDBBVx11VUuLzxt2jS2b98OwN69exkzZozzWFpaGiUlJTQ1NWGxWNi1axdTp07lww8/5I9//CMAAQEBKIoi7ZLFsGYyt1HX3kiCMRaDrv97Azc2QGWFQrDRt/bp/3HdgKJw9PUrlJX0+jQhxI+cciDx1ltv5cYbb6SgoACA1NRU/P39+3ThhQsX8s0335CZmYmqqqxdu5aPP/6Y9vZ2MjIyuO+++7jppptQVZXFixcTGxvLxRdfzP3338+yZcuw2Wz89re/7fP9hBiKSprLgTNbReBwQG42gMq4ieBL+bNRF44GLU22OudjMXFwKE+lohxSR3sxOCGGEJezigaDgQkTJrBnzx7eeustJk6cyMyZM11eWKPR8Mgjj3R7LC0tzfn1ggULWLBgQbfjgYGBPPPMM32NXYghr7S5AoDk0P5PERwphvY2hcRkldCwAQpsiNAoGkL04TRZ6zA7OvHT+KPRQFIyHC5QqCxXvR2iEENCr8nAv/71L377298SFRXFL3/5S5566immTp3Kpk2byMjI4P/+3//ryTiFGHbsDjsVrdWE+Ydg9Avq1zXMnV3FcnqDStoY1+efStbh+jO7gJeE6SNpstbRbK0nxq+rCDMxGYoKVUpLwOFQ0WgGvh20EMNJr8nA008/zcsvv0xLSws33ngjn332GYmJibS0tHDNNddIMiDEGaoy1WJz2EkOje/3NQ4XgMOuMGacil4/gMENIWG6KCCPphOSAYMfxMV31VH8kFfDjPFnvn+DEMPZKWcXx40bx8yZM0lJSSExseuHLCQkhMDAQI8EJ8RwVtpcCUByP1cRtJmgshwCg1Tiz2xV4pAWoA3CoPjRZK1HVY9PCySN6Prvp98WeycwIYaQXpOBE6v4/fz8PBKMEL6ktLkCrUZLnDGmX88vzO/q1pc+xreKBn9MURRC9ZHYVAvt9lbn4yGhXcsMdx2soq6pw4sRCjH49TpNUFtb69zM58Svj30vhOg/k6WNxs5mkkMT0PVj18GmRqitUQgNU4nqXy4xrITpI6m1VNBkrSdId3zrwcTkrpUWW3YeIXPhWC9GKMTg1uvniczMzB6/7ul7IcTpKXNOEfSvXuDw0fYg6WO71tb7utAf7TdwTGw8+Bu0fLGjBIdDVhYI0ZteRwZWrFjhyTiE8CmlLUeTgX4UDzY1QmODQkSkSli46/N9gUHjR6DWSIutEbtqR6t0jbbodDBvahKf7yhhb34t08bJMIoQPfHhmUYhvMPhcFDWUoXRL5gQP6PrJ/xI8dGufKPSBziwIa5ra2IHLdaGbo9fMrurkvCf3xd7ISohhoZekwFXPQGEEP1T01aH1W4lOSQe5TTH+JuboL5OISxCRgV+LEzf1VTtx1MFo5PDGJUQwn+yq2hs7fRGaEIMer0mA8uXLwdgzZo1nopFCJ9wpLn/UwTOUYG0U5/ni0KObU1srev2uKIoLJw5ArtD5asfyr0UnRCDW681A+3t7dx99918/fXXmM3mk46vW7fOrYEJMVyVtVSiUTQkGE9vI5zWVqir7VpBEB7hpuCGsG5bE9s78NMGOI/Nm5rIXz7K4stdpVx5vmRSQvxYr8nAyy+/zI4dO9i9e3efehEIIVzrtJmpa28gPjgGvfb0tgwsLer678hUWUHQm3B9FE3WOpqs9cRqjzd/Cg32Y8b4WHZkV1FU0cyoBB/p8SxEH/WaDMTHx3PllVcybtw40tLSKCoqwm63M3r0aHQ6l/2NhBA9qGytASAx5PRGBcydUFXZtdtgZLQ7IhsejtcN1BHr370T5IU/SWZHdhXbdpVy0xWSDAhxIpfv6larlUsuuYSwsDAcDgd1dXU899xzTJ482RPxCTGsVLRWA5BgjDut55Ue6dptMGWkKqMCp+CvCcRP40+TrR5VdaAox8uiZoyPwxio518/lPHLyyeg1cpiKiGOcfnT8Pjjj/P000/z/vvv8+GHH7Jx40YeffRRT8QmxLBT0VKNTqMjOqjvk/52G5SXdnUmjOt/p2OfoCgKYfoo7KoNk62l2zG9TsO8qUk0tZrZky+7qApxIpfJQHt7e7dRgClTpvRYUCiEOLV2SweNnc3EB0ejPY0tiCvKwWZVSEoG7envXOxzjk0VNP5oVQHAghnJAGzdecSjMQkx2LlMBkJDQ9myZYvz+y1bthAWFubWoIQYjpxTBKdRL6CqUFoCGo1KUoq7IhteQnURgEKzrf6kY6OTw0iMDuY/2VW0d1o9H5wQg5TLZODRRx/lj3/8I7NmzWLWrFm8+OKLPPLII56ITYhhpbwf9QINddDRrhAbDwZpHtonOo0eoy6UVlsTNkf3N3xFUTh/aiIWm4Md2VVeilCIwcdlAeHIkSN55513aG9vx+FwEBwc7Im4hBh2Klqq8NMaiAzs+8ha6dHRbBkVOD1h+ihabU009TA6MG9aEm9+nsf2PeXMn57sheiEGHz6XE4bGBgoiYAQ/dRiNtFqaSPeGING6duPXXOznfpaCAlVCZGVcKclzNnF8OS6gcToYNKSQtmTV0NLm8XToQkxKMnaGiE8oKLl9OsF9mebAYXkEW4KahgL1oaiU/Q0WetR1ZNbF8+bkoTdofLN/govRCfE4CPJgBAeUNHaNT+d2Md6AatVJfugGb1BJeb0tiQQdNUGhOojsTg6KW89uTZg7pREALbvKfN0aEIMSi5rBnJycnjxxRdpbm7ulmG/9tprbg1MiOFCVVXKW6sJ0PsT5h/Sp+fkF1gwm1VGpoJGUvZ+CdNHUm+pYl9lDkkh3ZtCRYcHMDE1kuzD9dQ3dxAZGtDLVYTwDS6TgVWrVpGRkcHo0aNPu92qEAKaOlvosHaSHjGizz9D2Qe79vJIkPq2fju238C+qhwuH3vhScfnTU0k+3A9X+8t58rz0z0dnhCDistkwN/fn+uuu84TsQgxLFWauvoRxPexS2FDg53KKhspyToCAmQtfH/5afwJ0AaTU3sIi92K4UeNoc47O4E/fnCAr/ZIMiCEywHIOXPmsGnTJoqKiqioqHD+EUL0TbWpa+vbuOC+dRjKzu0aFZg4XjYWOFNh+kgsdiu5tQUnHQsN9mPKmGgKSpuoqDV5ITohBg+XIwN/+9vfAPjf//1f52OKorB161b3RSXEMFLVWouf1tCnegG7XeVgnhl/f4XUUQZySzwQ4DAWro+isrOEvVU5nB03/qTj509N5IfcGr7aU86Si8d6IUIhBgeXycC2bds8EYcQw1K7pYNWSxspoQl9qhcoKrbS0aEy5Ww/dFqp0TlTRl04eq2efVU5wOKTjs+eFI9Bt4/te8rIXDhG6qKEz3I5TdDQ0MDtt9/OrFmzmDFjBitWrKCu7uSNPIQQJ6s63SmCgzJFMJC0ipYJ0aMpba6gob3ppOOB/npmTIilrMZEUUVLD1cQwje4TAZWr17NWWedxdatW9m2bRuTJ0/mgQce8ERsQgx5x5KB2D4kAyaTg5JSK7ExWqIiXQ7aiT6aHDcB4OjowMnmTU0CZM8B4dtcJgOlpaXcdNNNBAcHExISws033ywFhEL0UbWpFo2iIToowuW5OblmVFVGBQbalGPJQPXBHo/PGB9LgJ+O7XvLcThO3q1QCF/gMhlQFIXKykrn9xUVFeh08qlFCFesdht17Y1EB0ag05z6Z0ZVVXJyzeh0MGa0JAMDKTEkjsiAcA5UHcThcJx03E+v5Zyz4qlt7CD/SKMXIhTC+1y+q//mN78hIyODyZMno6oq+/bt49FHH/VEbEIMaTVtdaiofZoiKKuw0dziYPxYA34GKWIbSIqiMDluPNuKvuVw4xHSI0eedM7cKYls21XK13vLGTfS9SiOEMONy5GB+fPn8+GHH7J48WIWLVrEhx9+yAUXXOCB0IQY2o7vLxDl8tzsHCkcdKfJ8aeuG5g8OprgAD3/3lchUwXCJ/WaDGzevBmAjRs38uabb5KVlUVOTg5vvfUWGzdu9FiAQgxVVaauVTeuRgbMZgcFhy2EhWlIiJcpOHc4K3YciqL0mgzodRrOOSuehpZOcorqPRydEN7XazLQU9tPIUTfOFQH1W11hPoZCdD7n/Lc3EMW7HaYOM5P1rm7SbAhiPSIkeTXF9Fu6ejxnGOdDL/eW+7J0IQYFHr9GJKZmQlAYmIiV111Vbdjb7zxhnujEmKIa+xoxmq3EhvmutNQdo4ZRYHx42SKwJ0mx03gUH0RWTV5zEyactLxs9OjCA028O3+Sm658iy0WmkXKXxHr8nAK6+8gslk4q9//Svl5cczZbvdzscff8yyZcs8EqAQQ5FzsyHjqacIampt1NbZSR2lJyhQ3nzcaUrcBN7N/jt7q3J6TAa0Wg3nnpXAp98Vk1VYz+QxfdsoSojhoNffPiNGjOjxcYPBwPr1690WkBDDQV+bE8mOg56TFjGCIH0A+6pyep0GdU4V7JOpAuFbeh0ZmD9/PvPnz+enP/0paWlpnoxJiCGvylSLv86PUD9jr+fYbCp5hywEBSqMTNH3ep4YGFqNlrNix/N92Q9UmmpI6KGl9ITUSMKNfny7v4JbF52NTqYKhI9wWbp8880391jUJF0LheiZydKOydLOiLDEUxYEFhy2YDarnDXVH41GCgc9YXJcVzKwrzKnx2RAq1E4b3ICn/y7iH2Hapk+7uRzhBiOXCYDmzZtcn5ts9n44osvsFgsbg1KiKHsdKcIJsgUgcec2Kfgp2Pm93jO3CmJfPLvIr7eWy7JgPAZLsfAEhMTnX9GjBjB//k//4ctW7Z4IjYhhqS+NCdqbrZTVm4jMUFHeJjWU6H5vKigCBKNcWTX5GO1W3s8Z9yICKJC/fn+QCVWm93DEQrhHS5HBnbu3On8WlVVDh06hNlsdmtQQgxlVaZatIqG6MDet7XNyT1aOCjLCT1uctx4/nHoS/LqCpkUO+6k4xqNwpwpiXz4VSF78mqZOTHOC1EK4Vkuk4Fnn33W+bWiKISHh8tqAiF6YbFbaWhvIjY4Cq2m50/8DodKTq4Fg0EhPc3g4QjF2XET+MehL8mqyesxGYCuqYIPvyrk673lkgwIn3BaNQNCiFOraat32ZzoSKkVU5uDSRP80OulcNDTxkeno1E0ZFXnw1k9nzM6OYzYiEB2ZFdittrx08tUjhjeXNYMVFVVsWLFCmbNmsV5553H3XffTUNDgydiE2LIqWqtAU7dnChL9hbwqgC9P+kRIyloKKbd2vPWxIqiMGdyAh1mO7sPVns4QiE8z2Uy8Nvf/pZzzz2XrVu38tlnnzFp0iTuv/9+T8QmxJBT7aI5UXu7g6JiK5ERWmJj5NOmt0yKHYNDdZBbW9DrOdKrQPgSl8lAQ0MDS5cuJTg4mODgYH75y19SVVXlidiEGFIcqoOatjrC/EPw1/X8qT8334zD0TUqIE2JvGdSTFetQFZ1Xq/npCaGkhAVxM6D1XSabZ4KTQivcJkMnH322fz97393fv/ll18yadIklxd2OBysXr2ajIwMli9fTklJSbfj27ZtY/HixWRkZPD2228DYLVaueeee1i6dClXX321bGwkhpSG9iasDluv+wuoqkr2QTNaDYwbK4WD3jQmKhW9RkdWTe/JgKIozJ2SiNliZ2eOTBWI4a3XAsJx47r6f6uqyttvv82DDz6Ioii0t7cTGhrK448/fsoLb9myBYvFwubNm9m7dy/r16/nhRdeALre9NetW8e7775LQEAAS5YsYcGCBXz11VeEhYXx5JNP0tTUxJVXXsmFF144sK9YCDdxtb9AVbWNhkYHY9INBPjLNrfeZNDqGRuVRlZNHq1mE0a/4B7Pmzslkc1b8vl6XzlzpyZ6OEohPKfXZCA3N/eMLrx7927mzp0LwJQpU8jKynIeKywsJCUlhdDQUACmT5/Ozp07ufTSS7nkkkuArk9RWq3MqYqho8rFzoPZB7t27pQdBweHSbFjyarJI7smn9nJ03o8Z0R8CMmxRnYdrKa900qgv/SQEMNTr8nA5s2bycjIYOPGjT0eX7FixSkvbDKZCA4+nm1rtVpsNhs6nQ6TyYTReLyBS1BQECaTiaCgIOdzf/3rX3P77be7fAEbNmzoNUYhPEVVVapNtQTo/Anp4VOmxaqSX2DGaNSQkuRyRa/wgEkxY4GuuoHekgHoGh1487NcdmRXMX96sqfCE8Kjeh2r7K3FZ18FBwfT1tbm/N7hcKDT6Xo81tbW5kwOKisruf766/mv//ovfv7zn7u8z8qVK8nLy+v2R2oNhKeZLO20WTuIDY7qsTDwUIEFqxUmjJPCwcEiNWIE/jq/U9YNAMydkgDIqgIxvPX6ESUzMxOA8vJy1q1bd9oXnjZtGl9++SWXXXYZe/fuZcyYMc5jaWlplJSU0NTURGBgILt27eKmm26irq6OG2+8kdWrV3POOef04+UI4R2umhM5mxKNk8LBwUKn0TIhejQ/VGbR0N5ERGBYj+clxRhJTQjlh9waWtoshATJv6EYflxWMeXn53f7FN9XCxcuxGAwkJmZybp167j//vv5+OOP2bx5M3q9nvvuu4+bbrqJzMxMFi9eTGxsLC+++CItLS08//zzLF++nOXLl9PZ2dmvFyaEJ52qXqChwU5llY2UZB0hRqmDGUwmxR6dKnAxOnD+tCTsDpVv9snogBieXE5eKorC/PnzGTVqFH5+xwufXnvttVM+T6PR8Mgjj3R7LC0tzfn1ggULWLBgQbfjDz74IA8++GCfAhdiMOlqTqQlMjD8pGPZsuPgoOWsG6jJY97IWb2ed/60RF75ezZf7i7jp+eO8lR4QniMy2Tg3nvv9UQcQgxZFpuFho4m4oNjTmpOZLOr5OSZCfBXSB0lw8uDTUpYIsGGILKq81BVtdd6jsjQAM5Ki2J/QR1V9W3ERQZ5OFIh3MvlNMFnn33GzJkzu/157733PBGbEENCdVvvWxAXHrbQ2akyfpwfOq0UDg42GkXDxJgx1LU3OP8dezN/ehIAX+0p80RoQnhUryMDDzzwAKWlpWRlZXHo0CHn4zabjdbWVo8EJ8RQcKwfQU/NibJyuqYIJk2QKYLBalLMWHaU7SGrOq/XAlCAc85K4IX39vOv3WVce+EYWRUihpVek4HbbruN8vJyHn/88W57Cmi12m5z/0L4ut52HmxsslNWbiMpUUd4mBQODlYnFhFelDan1/OCAvT8ZGIc3+yroLC8mfSknlcfCDEU9TpNkJSUxKxZs/joo48YOXIkM2fORKPRkJubi8Egc59CQNf+GTVtdYQHhOKn6/5zkS2jAkNCgjGW8IBQso/WDZzK/GldUwX/2i1TBWJ4cVkz8PDDD/PCCy9QUFDAXXfdRXZ2NqtWrfJEbEIMevUdjdgc9pOGl212lZxcM/7+CmmpkjwPZoqiMClmLM3mVkqbK0557rRxsRgD9WzfU4bdcWYbswkxmLhMBg4cOMDq1av59NNPufrqq1m7di0VFaf+gRHCV/Q2RXC4yEJHp8qEsVI4OBScuMTwVPQ6DXMmJ9LYamb/oVpPhCaER7hMBux2Ow6Hg61btzJv3jw6Ojro6OjwRGxCDHq9bTYkhYNDy8SjdQPZNfkuz73g6KqCf/0gUwVi+HCZDFx55ZXMmTOHxMREJk+ezKJFi8jIyPBEbEIMal3NieoI1AdgNBxfd97UbKe0zEZigo7wcCkcHApigiKJDookp/YQDtVxynPHj4wgJiKQ7w5U0GmxeShCIdzLZTJwww038O9//5vnnnsOgDfeeINf/OIXbg9MiMGu1dJGew/NiWRUYGiaGD2GNks7R5pOPQ2qKAoXTEuiw2znP9lVHopOCPfqdWnhQw89xKOPPsry5ct7XE/rajtiIYa7nqYI7McKB/0U0qVwcEiZEDOafxV/R3ZNHiPDk0557gXTknh7Sz5f7i5j3tRTnyvEUNBrMnBsKmDlypUeC0aIoaSnToWHi6x0dKhMneyHTieFg0PJxJiuzqrZtYe4fOyFpzw3OdZIWlIoe/JqaGo1E2aUUSAxtPU6TTBp0iQAwsLCaGhooL29nfj4eOeWxEL4uipTLTqNlsiA482JDuR0ddmcNN7fW2GJfooOiiQmKJKDNfk4HKeuGwBYMD0Zu0OVQkIxLPSaDNTX17Ns2TKuu+46/vKXv/D888+zaNEibr75ZlpaWjwZoxCDjtlmobGjmZigKDSarh+j+objhYMREVI4OBRNjBlLm7WDkmbXrYrPn5aETquwdecRl5sVCTHY9TpN8OijjzJ9+nReeeUV9Ho9ABaLhQ0bNrB27VrWr1/vsSCFGGz+U3AYAI01iKzD9QDkZgMoRMRanY+JoWVizBi+LPqW7Jo8RoUnn/Lc0GA/fjIhju8OVFJQ1sTo5JPbVwsxVPQ6MpCXl8edd97pTAQADAYDd955Jzk5OR4JTojBqsXWBIBR1/UGYLVCZQX4+6tE9d7rRgxyE2JGA33bbwBg4cwUALb854jbYhLCE3pNBvz8ei6IURTFOSwqhK9qtTUCYNSFAlBRBg67QtIIkB+PoSsqMILY4GgO1hb0qW5g2tgYIkL8+GpPORar3QMRCuEevf7aOlV7TmndKXyZzW7DZGsmUGtEp9GjqlB2BDQalYREb0cnztTE6NG0Wzsobip1ea5Wq2H+9GTaOqx8n1XpgeiEcI9eawYOHTrEhReevLxGVVVqa2VPbuG7ippKceAgRNfVwrauBjo7FBKSVPSytcCQNzFmLNuKviW75hCpESNcnn/hT1J478sCvvjPEdlzQAxZvSYDn332mSfjEGLIyK0tBMB4NBkoLel6PNn1+4YYApz7DdTk8fNxF7k8PznWyIRREezNr6Wqvo24yCCXzxFisOk1GUhMlPFOIXqSV3csGQjH1AqNDQrhkSrBRi8HJgZERGAY8cExHKwrwO6wo9W4XiZ6yeyR5BQ18PmOEq6/bIIHohRiYEmpkxCnQVVV8uoKMWj88dP4y6jAMDUhZgwd1k6Km/q2odB5kxMICtDzxX+OYLO7LjwUYrCRZECI01BtqqXZ3IpRF4bNqlBVAf4BspxwuDlxqqAv/PRaFsxIpqnVLM2LxJAkyYAQpyH36BRBiC6MinJwOBSSU0AW2Awvx/cbONTn51wyu2t46J/fFbshIiHcS5IBIU7DsWQgWBNOWQlotCrxUkA+7EQEhBFvjCG3tqtuoC9GxIUwfmQEe44WEgoxlEgyIMRpyKsrxF/nh6k+mM5OhfhEOGGTTjGMTIwZS4etk6JG1/sNHHPpOTI6IIYmSQaE6KNWs4nylipGR46itEgDqKSM9HZUwl0mHp0qyOpj3QDAnMmJhAQZ+HxHCWbZkVAMIZIMCNFHx6YIwpR4WlsUYuIgMNDLQQm3mRjdVUSY08c+BQAGvZZLZo+gtd3K13uktbEYOiQZEKKPDtYWAFBc0LU9x4hR3oxGuFtYQCiJxjhy6wqx9bFuAODSc0aiUeDjfxdJa2MxZEgyIEQf5dYWoFW05OdCeIRKSKi3IxLuNjFmDJ02M4cbSvr8nJjwQGZNiudweTN5JY1ujE6IgSPJgBB90GkzU9R4BH97BDi0MirgIyYc3W8gp7bvSwwBLj+v63+QT/5dNOAxCeEOkgwI0QeH6ouwqw6aq4NITQwlIsrbEQlPOL7fQN+LCAHOTo8iOdbIv/eVU9/c4Y7QhBhQkgwI0QfH6gXsrRFkXDRGNhnyEWH+ISSFxJNbe3p1A4qi8F/zUrE7VBkdEEOCJANC9MH+iq5PhvEBScyeFO/laIQnTYgZjdluOa26AYD505MJDTbw6XfFdJht7glOiAEiyYAQLtgcdgoai3G0B5M5fxIajQwL+JJJMWOB09tvALqWGV5+7ijaOqxs+c8Rd4QmxICRZEAIF/YcOYQDG/7WaOZOkdbevmZCdFfdQM5p9Ck45rLzRmHQafjo60LsDllmKAYvSQaEcOH9nf8B4Ly0SWi18iPja0L8jSSHJpBXV4jNfnrD/aHBfsyfkUxVfTs7sirdFKEQZ05+swlxCjUN7Rxq6Np5cNHMWV6ORnjLxOgxmO0WCk6zbgDgyvPTUBR4d9sh2YRIDFo6bwcgxGD21y/yUIIaCdaFEmuM9HY4oh/+8MWHZ3yNekvX5kF/+u5DkgPSej1vXPCUHh9PTQjlUGkTf/5bFsmxxjOO59JzRp7xNYQ4kYwMCNGLyro2tmXloOitTEkY6+1whBeF6CIAaLE29Ov508bFALA7t3rAYhJiIEkyIEQv/vpFHgR1/fKfEJPu5WiEN+k1BgK1wbTamnCojtN+fkx4ICmxRspr26isa3NDhEKcGTRMvvEAACAASURBVEkGhOjBkaoW/rW7FGNMKwDjj1aUC98VqovAgQOTrblfz58uowNiEJNkQIgevPaPgzhUFW1IA+H+oSQYY70dkvCyEH3XVEGzrX9TBQnRwcRHBlFS1UptY/tAhibEGZNkQIgfySmqZ0d2FWlpWtrtbUyMHYsi+w/7vBBdOADN/awbAPjJhK6k8j85MjogBhdJBoQ4gaqqvPr3HAAmT+v68Ti2A53wbV11A0ZMtiYcat/7FJwoKSaY+KggiitbqKqX2gExeEgyIMQJduZUk1PUwKyJcdTbygCYdLSNrRChunAcOGjtZ92AoijMmhgHdP2/JsRgIcmAEEfZ7A7+95NsNAos++lYsmvziQ6KJCZY+hWLLqH6rr0mWvpZNwCQGB1MUkwwR6pbZWWBGDQkGRDiqH9+V0xZjYlLZo9E9WuhzdIuUwSim4GoGwCcowPfZ1XKroRiUJBkQAjA1G7hzc/yCPDTsfSScWTX5ANSLyC602n0BGlDaLU1YVf735Y4LjKIEXFGKuraKKlqHcAIhegftyUDDoeD1atXk5GRwfLlyykp6b6n97Zt21i8eDEZGRm8/fbb3Y7t27eP5cuXuys0IU6yeUs+re0Wrr1oDGFGP2e72omxUi8gugvTR6KinvHowDlnxaMA3x6owCEdDYWXuS0Z2LJlCxaLhc2bN3PXXXexfv165zGr1cq6/9/enYdHVZ6NH//Omkxmsu+QBEhI2EJI2BEQCSJCFSpUQSgupW6XYl9f9cW1olZsX1usRcuvbpeKfd2tlCKtbMoeEiCEhC0EsgHZt5lJMpnl/P6IhEXWZCbr/bkuLpI55zznnjuTzD3nec7zvPoq77//PqtWreKzzz6joqICgHfeeYfnnnsOm83mqdCEOM+pcgv/2nacsCAfZk6MxeFycqg8l16+4QQZAjo6PNHJBOiax5DU2Cvb1E6wv4GBfYOorrNxuKBthYUQbeWxYmDPnj1MnDgRgOTkZLKzs1u25eXlERMTg7+/P3q9nhEjRpCeng5ATEwMK1as8FRYQpxHURT+39dZOJwKv7plCHqdhuNVBTQ6bNJFIC7KVxuAGg019oo2tzVmSARajYq0nBLsjtbdriiEO3hs1UKLxYLJZGr5XqPR4HA40Gq1WCwWfH3PrtxlNBqxWCwATJs2jeLi4qs+z4oVK3jzzTfdF7joUXYcOM2+o+WkJIRyXVIkgHQRiMtSq9T464KotpfT6GzAW2NodVtGg47k+FAyDpex70g5o38cWChEe/PYlQGTyYTVeva2GZfLhVarveg2q9V6XnFwLRYvXsyRI0fO+7dx48a2BS96hEabg3e/OYBWo+bB2UktswzmnCkGQqUYEBd3tqug7VcHUgaE4eOtZe+RMuqsTW1uT4jW8FgxMHz4cLZs2QJAZmYmCQln/7DGxcVRUFBATU0NTU1NZGRkkJKS4qlQhLiozzYcpaK2kdmT+9MrtPkqls3RxOGK4/Tx742fd9vXnRfdU8CP8w24oxjQ6zRcN7QXTpfC9qxTbW5PiNbwWDfB1KlT2b59O/PmzUNRFJYtW8aaNWuor69n7ty5PPXUUyxatAhFUZgzZw7h4bIQjGg/RaVmvvnhGGGBBm6fcnZFwpyyo9iddoZFDunA6ERnZ9AY8VIbqHVUoSguVKq2fa5KiAkg53gFx0/WUlRqJjpcClHRvjxWDKjVal566aXzHouLi2v5OjU1ldTU1IseGxUV9ZPbDYVwF0VR+Ns/mgcN3vfzoXjrz/4aZJ7OAWC4FAPiCgJ0IZTaijA7avHTBbapLZVKxcTk3ny+MZetmSeZOzUBjVqmgRHtR15tosfZtv8U+3MrGDkovGUmOGguEvadzsag8yYhJO4yLQjh3q4CgNBAHxJjg6k229h3pNwtbQpxtaQYED1KfaOdd1dno9Oquf/nQ89bmvi0pYxSawVJ4YPQqjUdGKXoCvx1wahQtXm+gXONTYzEx1tLxqFSaswy14poP1IMiB5l1bpDVNU1MmdyPJEhxvO2nekiSJEuAnEVtCotvtoALM5a7C733AXgpdcwMbk3TpfC93uLZd0C0W6kGBA9Rs7xSv617QRRYabzBg2ese9088RYyRFSDIirc7arwH1XB+J6+9M30o+T5RYO5cvMhKJ9eGwAoRDt7d878y+5zeF08en65vkDxgyJYFNG0Xnb7a4mskuOEqQLY/f+GqDGc4GKbiNQF0ZhwzGq7GWEekW6pU2VSsX1Kb05WW5he9YposN98fXRu6VtIS5FrgyIHmF3Tgm1liaGxYcQEWz8yfbTtkKcOIkyxHZAdKKr8tGY8FIbqLFX4FJcbmvX10fP+KReNNldbM4oku4C4XFSDIhur7Sqnsyj5fgZ9YwZcvFPb8UNxwGI8pZiQFw9lUpFkC4Mp+KgzuHeS/qD+wURE+FLUZmFnOPu64YQ4mKkGBDdmtPpYlNGEQqQOiIanfanL3lFUShuOI5e5UWYV+/2D1J0aUH6MAAqm8rc2q5KpSJ1RDReOg3bs05TY5G7C4TnSDEgurWMw2VU1TWSGBtM7zDTRfepdVRhcdbSy7sv6jbOJCd6Hj9tAFqVjuqmMrdfzjcadFyf0huH08X6tAKcLvd1RQhxLvnLJ7qt8poG9h4uxeSjY9zQSw/uKmo4BiDjBUSrqFRqAnWhNCk2KptK3N5+QkwgA/sEUlbdwK4D7m9fCJBiQHRTDqeLDbsLcSkweXgUet2lJxEqqD+KChUxhv7tGKHoTs50FRQ05Hqk/YkpvQnw9SIzt5z803UeOYfo2aQYEN3SzgOnW7oHYiL8LrlfvcNMWdMpIryi8db4tGOEojsJ0AWjQk3hj1eZ3E2v1XDTmD5o1Co2phdSVlXvkfOInkuKAdHtFJyuI+tYBYG+XlyX1Ovy+/74Sa6PT8Jl9xPicjQqLQG6YKrt5dQ5PDNHRWiAgQnJvWlscvLqh7tpsjs9ch7RM0kxILqV+kY7GzOKUKtVTB3T56J3D5wrv/4oADGGn85IKMS1CNI1dxUU1numqwBgSL8gBvUN4lhxLSu/ypL5B4TbSDEgug1FUdiUUUSDzcG4xEhCAwyX3b/R2UCJrZAQfSQm7aW7EoS4GoH6UAAKGo567BxnZifsH+XPhvRCvt2R77FziZ5FigHRbWTnVVJQYiY6zMSw+JAr7p9ffwQFhX4+A9ohOtHd6dVeRHjFUGorxuLw3CA/rUbN03ePxt+k5+1vDrD3iHvnNxA9kxQDolsoKKlje9YpvPUapoyKOW9p4ks5Xn8QgH4+gzwdnughYo3Nr6UT9Yc8ep6wIB+eu3cMGrWKP3yUTkGJ3GEg2kaKAdHlNTY5+OPHe3C6FCaPjMZo0F3xGKvDTImtiHCvKOkiEG7T1zAAFWqOWz1bDAAM7BvEb+amUN/o4KX30qiua/T4OUX3JcWA6NIUReGtL/eTf7qOxLhgYnv5X9VxZz65xfoM9mR4oofx1hjo7d2XSnsptXbPLz88aXgUC24eSFlVPS+8sxNrg93j5xTdkxQDokv7dvsJvt9TzIA+gUwYdvnbCM+VV38IFSr6yngB4WaxxuYCM8+a0y7nm3tjAtPH9eXEqTpefj9NbjkUrSLFgOiyDudX8e4/s/E36XnqrlFo1Ff3cq5qKqOyqYQo71gMMtGQcLM+hni0Kj3HrNntcuufSqXigdlJjE/qRc7xSv7wUQZ2h6xhIK6NFAOiS6o2N/Lqh+m4XApP/nIkIVe4jfBcudYDAMSbhnoqPNGD6dR6+vkMwOKs47StoF3OqVGreHzBcJLjQ9l9sIQ//j0Dh1MKAnH1pBgQXY7T6eK1VXuoqmtk4YzBDIsPvfpjFSd51hy81QaiZS0C4SHxxiQAci0H2u2cOq2GZ381mqFxIezIOs3y/9uLUwoCcZWkGBBdzkffHuJAXgXjhkYyZ/K1vaEXNRyj0dVAnHEIGtWlFy8Soi3CvXrjpw0kv+EoTS5bu53XW6/l+UVjGBIbzNbMk7z28R7pMhBXRYoB0aV8v6eIr78/Rq8QI7+Zm3JV8wmc64glEzj7yU0IT1CpVMQbk3AqjpZuqfZi8NLy20VjSIwLZnvWKZZ9sBubDCoUVyDFgOgyDp6o5I3PMjF6a3nuV2Ouaj6Bc9XaKznZmE+EVzRB+qvvWhCiNRJMSajRcNi8r93XEPDx1vHCr8cyfGAYGYdKefGdXXLbobgsKQZEl1BSaWXZB7txKQpL7hpFdLjvNbdxyLwPgEG+w90dnhA/YdD40M9nILWOKk6100DCc3nrtTx372iuS4rkQF4FT721jcrahnaPQ3QNUgyITq/O2sSL7+6i1tLEg7cNJWVA2DW3YXfZyLUewEdjoo+sUCjayZnC85B5b4ecX6fV8D8LRzHjur7kn67jib9slamLxUVpOzoAIS6nscnBy+/torjMws8nxTH9un6taueIZT92pYmhpjGoZeCgaCeh+khC9JEUNuRSa6/CXxfklnb/vTP/mvbvG+nH2MQIdmWX8NjrP3DjqBhie1/dbJ03j+t7zfGJrkeuDIhOy+l08ceP93C4oJpJKVHce8uQ1rWjOMg2p6NV6RhkSnFzlEJcmkqlYqjfaAAO1KV1aBwjBoZz05g+KAqs25nP7pySdh/LIDovKQZEp+RyKbzx2T7SckpIjg/lN/NSUKuv7c6BM45Zc6h3WhhoSsZLc/WTEwnhDn0MCfhpA5tfhw5zh8YSHx3AnMn98fXRk36olG935Mv0xQKQYkB0Qoqi8Nev9rN5TzEDYgJ5+p5R6LSte6m6FBcH6tJQoyHRd5SbIxXiytQqNUP9RuPCyQFzekeHQ0iAgdunxNM71ET+6Tq+3JRLjbn95kIQnZMUA6JTURSFd1Zn859dBcT29mfp/ePw8b62WwjPddSaRZ2jmgTTUHy0134HghDu0N+YiFHjy2HzXqyOjh/AZ/DSMnNiLMPiQ6g22/h841EO51dJt0EPJsWA6DRcLoWVX2WxZutxYiJ8een+cZiucS6BczlcdvbVbker0pHsP96NkQpxbTQqLSn+E3DiZF/tjo4OBwC1WsWEYb2ZOjoGFbAxo4j/pBXQ2OTo6NBEB5BiQHQKzh/HCKzbmU+/Xn4se2g8/iavNrWZY86gwWlhsO9IfDQm9wQqRCv1Nybirw0i15pFrb2yo8NpkRATyNypA4gM9iGvuJZP1x+luMzS0WGJdibFgOhwNruTP3yUzqaMIhJiAtxSCFgcdeyv24m32kCS3xg3RSpE66lVakYETEJBYVf1hk51Sd7PqOfnk/ozZkgE9Y12Vm/JY3vWKVn5sAeReQZEh6q12Hj5/TSOFFST1D+EZ+8d3aYxAmekVW/EodgZFzgVvbpthYUQ7tLHEE9v776cbMynoOEofX0GdHRILdRqFSMHhRMVZmL97kIyj5Zz/GQtvUJNDG/FRF+ia5ErA6LDFJWaeXLFVo4UVHPD8CiW3jfWLYVAYX0uBQ1HCfeKor8x0Q2RCuEeKpWKsYFTUaMhrXpju65oeLUigo3MnZpAcnwo5vomXnh7J6+tyqC6rrGjQxMeJMWA6BA7D5zi8Td+4HSFldunxPPf84ej07Z9ZsAGp5VtVf9Gg4brAm+65lUNhfA0f10Qw/zHYnWa2VW9oaPDuSi9VsP4Yb24fUoCCTEBbMk8yUN/2Mi6HSdwujpP94ZwHykGRLtyOF18uPYgyz5Ix6XAk78cwV0zBrvlTVtRFLZVrqPRVc/IgEkEysqEopMa5jeOEH0kx6zZnKg/3NHhXFJogIH/XXw9D85OQgH++lUW/7X8ezIOlXaqMQ+i7aQYEO3mVIWF/1mxlS835RIR7MNriydyfUqU29rfX7eDosY8enn3ZbDvSLe1K4S7qVUaJgX/DI1Ky7bKddTYKzo6pEvSqFX8bHw/Vi6ZwpRR0RSU1PHiu7t4duUOjhZWd3R4wk2kGBAe53QprN12nN/86Xtyi2qYPCKKN/77Bvr1urqFUq7G7uJM9tZuw6jxY1LwLdI9IDo9f10wE4KmY1ea2FD+NTZX5+6TD/Lz5r/mDecvj09m5KBwDuRV8PgbW/j9h+kUykqIXZ7cTSA8Kv90HW9+kcmRgmqMBh1PLEhm0nD3XQ0AOFB6mDd2vY9WpePG0DkYNEa3ti+Ep8QZB1NlL+NAXRrry79kWugd6NT6jg7rsvpG+vHCr8dyIK+CD/91kO1Zp9iedYpRg8OZMzmewf2CpBjvgqQYEB5RXdfI3/9zmPVpBbgUuD65N7+elUign7dbz5NTdpT/3boSRVFIDbmNYL3cAiW6lhH+12N11HG8/hAbyr/mxtDZnb4gABgaF8Jrj04kLaeErzblkn6wlPSDpQzoE8jsG/ozJjESTSsXFxPtT4oB4Va1Fhurt+SxZutxGpucRIWZWDQzkZGDwt1+rq35u1mZvgoFhcevu5+KQj+3n0MIT1Or1Fwf/DMcioPChlzWlX3CjaFzusSsmSqVirGJkYxNjOTgiUq+3nyMtJwSXv0wnbBAA6kjY5gyKpqIYLla19lJMSDc4mS5hbXbT/CfXQU02Z0E+Hrxq1uHcNOYPmg07h2a0uS0839Z3/Dt0U346Aw8Mf5+EsMH8u/CfLeeR4j2olZpSA2Zxfaqf5NrzWZNyUdMCr6VCO/ojg7tqg3uF8zgfsEUlZpZvSWPH/YW8+n6I3y6/giJccFMGRnD+GG9MHjJ205nJD8V0WqNNgdpOSV8l1ZA1rHm0dAhAQbmTO7PjaNj8Na7/+V1qDyX9/Z8RmHtSXr7RvDf4+8j2r+X288jRHtTqzRMCJqBnzaIvbVbWVf2CYm+o0j2vw5dF5pFMzrcl0duT2bRzER2ZJ1iY3oRB/IqyM6rZOVX+0mKD2X0kAhGDQonJMDQ0eGKH0kxIK5JrcXGvqPlpGWfJv1QKbYmJ9Dcf3jzuD6MG9oLndb9N6mcqC7i64PrSCveB8DUuInclfwLvLSdv29ViKulUqkY5j+OCO9ofqj4FwfMu8mrP0iS3zgSjEPRqts+Q2d7MXhpmTIqhimjYiiptLIpo4jtWafIOFRKxqFSAGJ7+TNiUBiJsSEk9Als0yqlom2kGBCXpCgK5dUN5BbXcOhEFQdPVHKsuIYzc430CjEyMaU3k1KiiA73dfv5zTYLGSez+D5/J4fKjwHQP6gv9w6/g/jgfm4/nxCdRbhXFLMjF5FVl8YBcxq7qteTWbud/sYhxBmHEKQL61Ij9iOCjcyfNpD50wZSUmkl41DzYMOsYxUcP1XLFxtzgearCgP7BDKwbxD9evnRO9TklinKxZV5rBhwuVwsXbqUI0eOoNfr+d3vfkefPn1atm/atIm33noLrVbLnDlzuOOOO654jHA/l0uhztpEVV0jlbUNnCy3UFhibv5XaqbBdnZtc41axZDYYIYPCGPkoHD6Rvq57Q9Sk6OJsvpKjlcVcqwqn7zKfI5VF7TMcjY0fAC3DriJYRGDutQfQSFaS6vWMTxgAoN8U8gxZ3DEsp9sczrZ5nSMGl/CvaIJ0UcQoo8gUB+KXuXVJX43IoKN3DIhllsmxNJgc5CdV8HhgmoO51dxtLCaolIz63cXtuwf4u9NVJgvUeEmIoONBPsbCPb3JtjfQKCfF1o3j0nqqTxWDGzYsIGmpiY+++wzMjMz+f3vf8/KlSsBsNvtvPrqq3z55ZcYDAbuvPNOUlNT2bt37yWPEVdWXddIZm45TXYXDqcLu+PH/+1O6m0OrA126hsd1DfasTbaqTbbqK5rxOH86bSiGrWK3mEmYsJ96dfLn0H9goiPDmjTOIAmp51PslZT2VBNk9OO3dmE2WalsqEGs+389dM1KjXxQf0Y2TuJ66JHEGYKafV5hejKDBojIwMmkeI/nqKGPArqj1LceJzj9Qc5Xn+wZT+tSoePxoSPxoRe7YW/NoiRATd06gLB4KVl1OAIRg2OAMDpdFFQYuZIYTWFJXUUl1koLjWTmVtOZm75T45XqZqXXzYZdJgMeow+OkwGHUaDDr1Wg16nRqdRo9Np0GnV6LVqdFo1Wo0alUqFWqViSFwwQW6+5bkr8lgxsGfPHiZOnAhAcnIy2dnZLdvy8vKIiYnB3795BroRI0aQnp5OZmbmJY+5Fk5ncz92SUlJW55Cl/PBv3LYtv/UVe2r0ajxN+mJ9PMiwNeLAFPz/2GBPvQKMRIebLzgHuFGKsrals+qhhpW7/4Wl3J2jXS9RkeAwZ9+3r0INPgT5RdJn4Aoovwi0GmaLw821TRSXFN8xfYry9vv511XX9Vu5xJdR2WDZ1+DvviTyCiGeI3E6qij2l5Opb0ci72GBpcVs7OaClfz3wCtSkeMPQ6tqm3jaoqL27c3WQ8MjdYyNDoICAKaByufrrRSVWuj2txItdlGjdlGlbmROouZqnI7RTY7zot8sLmSkYPCeXB2knufRAc483535v3vWnnsp2yxWDCZzt4nq9FocDgcaLVaLBYLvr5n+5iNRiMWi+Wyx1zKihUrePPNNy+6bcGCBW54JkII0TUdYEtHh9DpndgEX7zV0VG4T3l5eau61z1WDJhMJqxWa8v3Lper5U39wm1WqxVfX9/LHnMpixcvZvHixec91tjYyLBhw/juu+/QaNq+LK64uClTprBx48aODqNbkxy3D8mz50mOPcvpdHLTTTeRmJjYquM9VgwMHz6czZs3M2PGDDIzM0lISGjZFhcXR0FBATU1Nfj4+JCRkcGiRYtQqVSXPOZaeHs39//I4EPPi4py7zoD4qckx+1D8ux5kmPPO/P+d608VgxMnTqV7du3M2/ePBRFYdmyZaxZs4b6+nrmzp3LU089xaJFi1AUhTlz5hAeHn7RY4QQQgjhWR4rBtRqNS+99NJ5j8XFxbV8nZqaSmpq6hWPEUIIIYRnyQ2aQgghRA+nWbp06dKODsJTxowZ09EhdHuSY8+THLcPybPnSY49r7U5VilnpngTQgghRI8k3QRCCCFEDyfFgBBCCNHDSTEghBBC9HBSDAghhBA9nBQDQgghRA/XbYqBv/3tb8ydO5fZs2fzxRdfUFBQwJ133sn8+fN54YUXcLlcV25EXJLdbufxxx9n3rx5zJ8/n7y8PMmxG+3fv5+FCxcCXDKvn3/+ObNnz+aOO+5g8+bNHRlul3Rujg8dOsT8+fNZuHAhixYtoqKiApAcu8O5eT5jzZo1zJ07t+V7yXPbnJvjyspKHnroIRYsWMC8efMoLCwEWpFjpRvYtWuX8sADDyhOp1OxWCzKX/7yF+WBBx5Qdu3apSiKojz//PPKd99918FRdm3r169XHn30UUVRFGXbtm3KI488Ijl2k7ffflu55ZZblNtvv11RFOWieS0rK1NuueUWxWazKXV1dS1fi6tzYY4XLFigHDx4UFEURfnkk0+UZcuWSY7d4MI8K4qi5OTkKHfddVfLY5Lntrkwx0uWLFHWrl2rKIqi7Ny5U9m8eXOrctwtrgxs27aNhIQEHn74YR588EFuuOEGcnJyGD16NADXX389O3bs6OAou7Z+/frhdDpxuVxYLBa0Wq3k2E1iYmJYsWJFy/cXy2tWVhYpKSno9Xp8fX2JiYnh8OHDHRVyl3NhjpcvX86gQYOA5tXevLy8JMducGGeq6urWb58Oc8880zLY5Lntrkwx3v37qW0tJR77rmHNWvWMHr06FbluFsUA9XV1WRnZ/PGG2/w4osv8sQTT6AoCiqVCgCj0YjZbO7gKLs2Hx8fTp48yfTp03n++edZuHCh5NhNpk2bdt5S3RfLq8ViwdfXt2Ufo9GIxWJp91i7qgtzHBYWBjT/If3444+55557JMducG6enU4nzz77LE8//TRGo7FlH8lz21z4Wj558iR+fn588MEHREZG8s4777Qqxx5bqKg9BQQEEBsbi16vJzY2Fi8vL0pKSlq2W61W/Pz8OjDCru+DDz5gwoQJPP7445w+fZq7774bu93esl1y7D5q9dka/UxeTSYTVqv1vMfP/WUX1+7bb79l5cqVvP322wQFBUmO3SwnJ4eCggKWLl2KzWbj2LFjvPLKK4wdO1by7EYBAQEti/6lpqby+uuvk5iYeM057hZXBkaMGMHWrVtRFIXS0lIaGhoYN24caWlpAGzZsoWRI0d2cJRdm5+fX8uLyd/fH4fDweDBgyXHHnCxvCYlJbFnzx5sNhtms5m8vDwSEhI6ONKua/Xq1Xz88cesWrWK6OhoAMmxmyUlJbF27VpWrVrF8uXL6d+/P88++6zk2c1GjBjBDz/8AEB6ejr9+/dvVY67xZWByZMnk56ezi9+8QsUReG3v/0tUVFRPP/88yxfvpzY2FimTZvW0WF2affccw/PPPMM8+fPx26389hjj5GYmCg59oAlS5b8JK8ajYaFCxcyf/58FEXhsccew8vLq6ND7ZKcTievvPIKkZGRLF68GIBRo0bx6KOPSo7bQWhoqOTZjZYsWcJzzz3Hp59+islk4k9/+hP+/v7XnGNZqEgIIYTo4bpFN4EQQgghWk+KASGEEKKHk2JACCGE6OGkGBBCCCF6OCkGhBBCiB5OigEherC0tDRSUlKYNWsWM2fOZPr06axcufK8fXbu3Mkvf/lLpk2bxtSpU3n00UfPm9Tr3LYuXKCmrQYMGODW9oQQF9ct5hkQQrReYmIiq1atAppnKpsxYwZTp06lf//+ZGRk8OSTT/Lmm2+SnJwMwN///ncefvhhvvrqq44MWwjhRnJlQAjRorGxEY1G0zLb5F//+lceeuihlkIAYMGCBcyYMYOmpqZLtlNQUMC9997Lbbfdxp133snBgweprq5m/PjxLdNYHz16lFtvvRWAb775httuu41Zs2bxzDPPYLPZPPgshRAXkmJAiB4uOzubWbNmceutt5Kamsro0aNbFvLJzMxk1KhRPzlm0aJF6PX6S7a5ZMkSnnzySf7xj3/w8ssv89hjjxEYGEhSUhLbtm0DYO3atcyc5tNydwAAAgVJREFUOZPc3Fw+//xzPv30U1avXk1wcDDvvfeeZ56sEOKipJtAiB7uwm6CBx98kLfffpsHHngAoGUFxaamJm6//XYAamtrWb58OcOHD/9Je1arlezsbJ5++umWx+rr66murmbWrFmsXbuWyZMns27dOj766CM2bNhAQUEBd9xxBwB2u53Bgwd79DkLIc4nxYAQooXRaOTGG29kx44dAAwdOpS9e/cSHx+PXq9n9erVACxcuPC8VSvP5XK5ztsXoKSkpGV1tVdffZX09HQiIiKIiIjA6XQyffp0nnvuOaC5mHA6nR5+pkKIc0k3gRCihdPpZPfu3S2fzBcvXsxbb73F/v37W/Y5fPgwRUVFaDSai7bh6+tL3759W4qB7du3s2DBAgD0ej0TJ05k2bJlzJw5E4AxY8awfv16KisrURSFpUuX8uGHH3ryaQohLiBXBoTo4c6MGQBoaGhg6NCh3HfffQCMHDmS119/nT//+c9UVFSgKAoBAQEsWbLksktWv/baayxdupR3330XnU7H66+/3tLdMGvWLP75z39y8803AzBw4EAeeeQR7r77blwuF4MGDeL+++/38LMWQpxLVi0UQgghejjpJhBCCCF6OCkGhBBCiB5OigEhhBCih5NiQAghhOjhpBgQQgghejgpBoQQQogeTooBIYQQooeTYkAIIYTo4f4/EHgoPM5jLDEAAAAASUVORK5CYII=",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f055a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's visualize the data\n",
    "import matplotlib as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set(style=\"ticks\")\n",
    "\n",
    "\n",
    "bins = np.arange(70,150,6)\n",
    "\n",
    "A = df['BG 1'].values\n",
    "B = df['BG 2'].values\n",
    "\n",
    "# Show the results of a linear regression within each dataset\n",
    "ax1 = sns.distplot(A,bins=bins,label='Before treatment')\n",
    "ax2 = sns.distplot(B,bins=bins,label='After treatment')\n",
    "\n",
    "plt.pyplot.xlabel('BG level')\n",
    "plt.pyplot.ylabel('Distribution of BG levels')\n",
    "plt.pyplot.legend(bbox_to_anchor=(0.45, 0.95), loc=2, borderaxespad=0.)\n",
    "\n",
    "plt.pyplot.xlim((60,160))\n",
    "plt.pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this data, we want to know if the two BG readings taken at different times (BG1 and BG2) for the 40 persons with diabetes are significantly different.    \n",
    "\n",
    "Let's use a paired t-test as defined above.    \n",
    "First, we calculate the difference between BG readings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BG 1</th>\n",
       "      <th>BG 2</th>\n",
       "      <th>Difference</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Person ID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>89.162573</td>\n",
       "      <td>94.023517</td>\n",
       "      <td>4.860944</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>90.857629</td>\n",
       "      <td>95.273755</td>\n",
       "      <td>4.416126</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>94.912999</td>\n",
       "      <td>96.612870</td>\n",
       "      <td>1.699871</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>95.254064</td>\n",
       "      <td>96.818673</td>\n",
       "      <td>1.564609</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>97.136291</td>\n",
       "      <td>97.760342</td>\n",
       "      <td>0.624052</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                BG 1       BG 2  Difference\n",
       "Person ID                                  \n",
       "0          89.162573  94.023517    4.860944\n",
       "1          90.857629  95.273755    4.416126\n",
       "2          94.912999  96.612870    1.699871\n",
       "3          95.254064  96.818673    1.564609\n",
       "4          97.136291  97.760342    0.624052"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def difference(a,b):\n",
    "    return b-a\n",
    "df['Difference']=df.apply(lambda row: difference(row['BG 1'], row['BG 2']), axis=1)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the mean and standard deviation of the differences gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "$\\overline{d} = -2.021$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$s_d = 2.972$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "diff = df['Difference'].values\n",
    "\n",
    "d = diff.mean()\n",
    "std = diff.std(ddof=1)\n",
    "\n",
    "printmd('$\\overline{{d}} = {}$'.format(round(d,3)))\n",
    "printmd('$s_d = {}$'.format(round(std,3)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, calculate the corresponding test statistic   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "SE = 0.47"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "$\\mu_1 = 109.73, \\mu_2 = 107.71$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "SE = std/(len(df))**0.5\n",
    "DF = len(df)-1\n",
    "\n",
    "# Mean of sample 1:\n",
    "mu1 = df['BG 1'].values.mean()\n",
    "mu2 = df['BG 2'].values.mean()\n",
    "\n",
    "printmd('SE = {}'.format(round(SE,3)))\n",
    "printmd('$\\mu_1 = {}, \\mu_2 = {}$'.format(round(mu1,2),round(mu2,2)))\n",
    "\n",
    "TestStatistic = (mu2 - mu1)/SE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "t-score -4.3"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "p-value = 0.00011"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import stats\n",
    "pvalue = 2.0* stats.t.cdf(TestStatistic, DF)\n",
    "# Multiply by two, since this is two-tailed test\n",
    "\n",
    "printmd('t-score {}'.format(round(TestStatistic,2)))\n",
    "printmd(\"p-value = {}\".format(round(pvalue,5)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The probability value for t-statistic = -19.04 is virtually 0, hence we can reject the null hypothesis and tell that the two BG readings are significantly different."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='scipy'></a>   \n",
    "## t-test with scipy    \n",
    "\n",
    "Let's repeat this test with scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Ttest_relResult(statistic=-4.3004899552379108, pvalue=0.00011038369115019783)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sample 1\n",
    "a = df['BG 2']\n",
    "b = df['BG 1']\n",
    "\n",
    "stats.ttest_rel(a,b)\n"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}