sliderule_dsi_inferential_statistics_exercise_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What is the True Normal Human Body Temperature? \n",
    "\n",
    "#### Background\n",
    "\n",
    "The mean normal body temperature was held to be 37$^{\\circ}$C or 98.6$^{\\circ}$F for more than 120 years since it was first conceptualized and reported by Carl Wunderlich in a famous 1868 book. But, is this value statistically correct?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"span5 alert alert-info\">\n",
    "<h3>Exercises</h3>\n",
    "\n",
    "<p>In this exercise, you will analyze a dataset of human body temperatures and employ the concepts of hypothesis testing, confidence intervals, and statistical significance.</p>\n",
    "\n",
    "<p>Answer the following questions <b>in this notebook below and submit to your Github account</b>.</p> \n",
    "\n",
    "<ol>\n",
    "<li>  Is the distribution of body temperatures normal? \n",
    "    <ul>\n",
    "    <li> Although this is not a requirement for CLT to hold (read CLT carefully), it gives us some peace of mind that the population may also be normally distributed if we assume that this sample is representative of the population.\n",
    "    </ul>\n",
    "<li>  Is the sample size large? Are the observations independent?\n",
    "    <ul>\n",
    "    <li> Remember that this is a condition for the CLT, and hence the statistical tests we are using, to apply.\n",
    "    </ul>\n",
    "<li>  Is the true population mean really 98.6 degrees F?\n",
    "    <ul>\n",
    "    <li> Would you use a one-sample or two-sample test? Why?\n",
    "    <li> In this situation, is it appropriate to use the $t$ or $z$ statistic? \n",
    "    <li> Now try using the other test. How is the result be different? Why?\n",
    "    </ul>\n",
    "<li>  At what temperature should we consider someone's temperature to be \"abnormal\"?\n",
    "    <ul>\n",
    "    <li> Start by computing the margin of error and confidence interval.\n",
    "    </ul>\n",
    "<li>  Is there a significant difference between males and females in normal temperature?\n",
    "    <ul>\n",
    "    <li> What test did you use and why?\n",
    "    <li> Write a story with your conclusion in the context of the original problem.\n",
    "    </ul>\n",
    "</ol>\n",
    "\n",
    "You can include written notes in notebook cells using Markdown: \n",
    "   - In the control panel at the top, choose Cell > Cell Type > Markdown\n",
    "   - Markdown syntax: http://nestacms.com/docs/creating-content/markdown-cheat-sheet\n",
    "\n",
    "#### Resources\n",
    "\n",
    "+ Information and data sources: http://www.amstat.org/publications/jse/datasets/normtemp.txt, http://www.amstat.org/publications/jse/jse_data_archive.htm\n",
    "+ Markdown syntax: http://nestacms.com/docs/creating-content/markdown-cheat-sheet\n",
    "\n",
    "****\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as scs\n",
    "import seaborn as sns\n",
    "from scipy.stats import norm\n",
    "import matplotlib.mlab as mlab\n",
    "import math\n",
    "%matplotlib inline\n",
    "\n",
    "df = pd.read_csv('data/human_body_temperature.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "First, take a quick look at the basic stats on the temperature data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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>temperature</th>\n",
       "      <th>heart_rate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>130.00</td>\n",
       "      <td>130.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>98.25</td>\n",
       "      <td>73.76</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.73</td>\n",
       "      <td>7.06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>96.30</td>\n",
       "      <td>57.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>97.80</td>\n",
       "      <td>69.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>98.30</td>\n",
       "      <td>74.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>98.70</td>\n",
       "      <td>79.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>100.80</td>\n",
       "      <td>89.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       temperature  heart_rate\n",
       "count       130.00      130.00\n",
       "mean         98.25       73.76\n",
       "std           0.73        7.06\n",
       "min          96.30       57.00\n",
       "25%          97.80       69.00\n",
       "50%          98.30       74.00\n",
       "75%          98.70       79.00\n",
       "max         100.80       89.00"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Your work here.\n",
    "df.describe().round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if the sample distribution is a normal one by fitting the data and plotting the along the fit nornal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean of fit normal distribution: 98.2492\n",
      "Standard deviation of fit normal distribution: 0.7304\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1055352e8>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEFCAYAAADqujDUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4HNW5+PHv7OyqS5ZsrYptNbdjcMPYBts4YHoPTgIh\nBFIIBEjg3tzwu+HmJrkQcuGmEAghtMtNIL1QQkx1KAaDDQaMsXE9Rq6SrG7Z6tKW+f0xKyNsSV6V\n1Wx5P8/jR9qdnZl3JGvemTnnvMewLAshhBCJx+V0AEIIIZwhCUAIIRKUJAAhhEhQkgCEECJBSQIQ\nQogE5XY6gHDV17cMubtSTk4aTU3tIxlO1JNjTgxyzIlhOMfs9WYa/S1LiDsAt9t0OoRRJ8ecGOSY\nE0OkjjkhEoAQQoijSQIQQogEJQlACCESlCQAIYRIUJIAhBAiQUkCEEKIBCUJQAghEpQkACGESFAx\nMxJYCCFGU3X1fr7ylSuYNk0dfm/evAUsWXIqq1e/wdVXf51Vq15jxoyZ5OZ6P7Hurl07eeih++js\n7KSjo4NFi07ha1+7DsPod1CuIyQBiLh348pbhrTeA2f8bNT2N9R9icgqLS3j/vsfOer9qVPtpPDE\nE3+htPR7n0gALS0t/PCH3+POO++iqKiYQCDAf/3Xd1m+/CmWLbt01GIPhyQAIUTUe3xlOe9trwPA\nNA0CgeHPZLhgeh6fP2PKoNdbv34dy5c/xbnnXkh5+Q7uuONWHnzwN3g8HgBWr17FiScuoKioOBSv\nyQ9+cDsej+fwurff/mMAPv3pc3nmmX9y550/5NChQzQ3H6K4uIS5c+dx/vkX0djYwHe+8288++xy\nHn74fjZu/IBgMMjll1/JGWecNeyfQcQSgFLKBTwIzAG6gGu11uW9lp8P3AYYwPvAjVprmZ9SCBE1\n9uzZzU03XXf49W233XH4+8WLlzBlyjS+853vHT75AzQ01DN+/IRPbCctLe2Y+5o3bz6XX34lu3fv\n4he/+Bnnn38R//znC1x44cWsWrWK6uoqHnroN3R1dXH99VezYMHJZGZmDuv4InkHsAxI0VovUkot\nBO4GLgFQSmUCdwFLtdYNSqlbgFygPoLxCCFi1OfPmHL4at3rzaS+vmVU9tvXI6CKin0DrpOfX8iO\nHds/8d7+/VXU1dUe9dnec7IXF5cAUFY2iUAgQE1NNa+++jL33vsgr776PFpvP5yM/H4/NTX7ycxU\nR21zMCLZC2gJsAJAa70WmN9r2WJgE3C3UupNoFZrLSd/IURMcblcBIPBT7x3yilLeOedt6iqqgTs\nk/WvfvULdu3aSVJSMo2NjQDU1FTT3Hzo8HqG8fHp+KKLLuHBB++jtLSMzMxMJk2axNy587n//ke4\n776HOeOMs5gwYeKw44/kHUAWcKjX64BSyq219mNf7Z8OnAC0Am8qpd7WWu/ob2M5OWnDKonq9Q7v\nVikWyTFHz7YiuS/5PUdGV1c6Ho951L6ys9NITvbg9WZy0knz+clPbufRRx8lOzv7cGw///ld/Oxn\nP8ayLNra2jj99NO57rqrCQQCPP54Nt/85teYPHkyRUVFeL2ZpKR4GDMm9fC+LrtsGffddzcPPfQQ\nXm8mZ5xxBu+++y7f+tb1tLe3c9ZZZ1FSUjDsYzR634KMJKXUPcBarfXjodeVWuuJoe/Pw37mf3Ho\n9S+BNT2f7ctwJoQZzVvGaCHH/LF47gUkv+fEMJxjdmpCmDXABQChNoBNvZatB2YqpXKVUm5gIbA1\ngrEIIYQ4QiQfAT0NnK2Uegu7p8/VSqmbgXKt9TNKqf8E/hn67ONa680RjEUIIcQRIpYAtNZB4IYj\n3t7ea/lfgb9Gav9CCCEGJrWAhBAiQUkCEEKIBCWlIIQQMWGovbn6I/WXJAEIIUS/1q9fx623/iel\npWWAPajrssuu4Mwzz+7z8zU1NZSX72DJklNHM8whkwQghBADmDdv/uHibe3t7dx003UUFxcfrgja\n2/r177F37x5JAEIIEW/S0tK45JLP8uqrL/PUU49TV1dLY2MDp5xyKtdccz1//ONv6ezsZNas2aSn\nZ/DYY/9HMBiko6OD226743C9n2ghjcBCCDEIY8eOZccOzYwZs7jnnvt55JHfsXz5U5imyVVXfZWz\nzz6PJUtOY/fuXdx6639z//2PcNppp/Paa684HfpR5A5ACCEGoaamhlmzZrNt2xbWr19Heno63d2+\noz7n9Xq59967SE1No76+jlmz5jgQ7cDkDkAIIcLU1tbKs88+TXp6BhkZmdx22x184QtX0dXViWVZ\nGIaBZdnVQX/60zv53vdu4/vf/+FRU0ZGC7kDEELEhJ5um6NdDO7999dx003XYZomgUCAa665nqKi\nEm6//Qds2bIJj8fDxIlFNDTUM3nyFH7/+0eZNm065557Pt/85tdJTU0hJ2ccDQ3RV/FeEoAQQvTj\nxBPn89xzL/e57He/+8tR73m9efzlL38H4Kyzzo1obCNBHgEJIUSCkjsAkXAsC4Kt2QQPegk0j4Pg\nxxMNuVJbcY1pwBzT4GCEQowOSQAiZgy3FIAVNPDXluCvLgN/sv2mEQSX3y5YbhkEOjIJHCjEB/yi\neSOfP30yE7wZw45diGgkCUAkhMDBXHz7jsPqTAfTh+mtwMyuw5V1AMMMAPadgdWZTvBQLoED+Wza\nBZt3N3LqnPF85lOTyEpPcvgohBhZkgBEXLOCBr59xxGoKwaCmHl78Uwox/Ac3W/bMMBIbcOV2oaZ\nv5frS27h8dfKWbVhPxvLG/iXz82mrDBr9A9CiAiRRmARtyy/m+4d8wjUFWOkNZM88y2SSrf1efI/\nkmHAnCm5/Oiak7h06WQOtXbzkz+t591ttaMQuRCjQxKAiEvBrhS6ti4k2JyLK7uO5OPewZXWOujt\nmC4XFyws4V8vnY3pMnh4+Raee2vPyAcshAMkAYi4Y/k8dG9fgNWZgZm/h6Sp6w8/5x+qOVNy+f6X\n5pE7JoW/v7GLV9ZVjFC0QjhHEoCIK1bApGvHPKyudNyFu0gq2Y5hjMy2J3gz+Pcr5jImPYk/v/IR\na7fWjMyGhXCIJAARN6ygQXf5CVht2Zi5Vbgn7hjxfeRlp/Ltz88hNdnNb57bxqZdjSO+DyFGiyQA\nETd8+6YTPOTFNaYOT+nmEbvyP1JxfibfunQ2LpfBQ//YTF1Te2R2JESESQIQcSFwIJ9AXQlGagtJ\nUzZiuKyI7m9aUTZfPW86nd0BHvrHFnz+YET3J0QkSAIQMS/YlUr37png8pM0ZcOwG3zDtWhmAUtm\nF7K3toXHV5aPyj6FGEmSAERMs5/7z4GAB0/JNlypbaO6/yvPnsaE3HReXV/Juu11o7pvIYZLEoCI\naf6qKXaj77j9mLlVo77/ZI/JN5bNJMnj4rcvbudQa9eoxyDEUEWsFIRSygU8CMwBuoBrtdblvZb/\nElgC9MzscInW+lCk4hHxJ9ieib+6DCOpHU/plog1+h7L+Nx0Lj99Cn94aQd/euUjkGoRIkZEshbQ\nMiBFa71IKbUQuBu4pNfyecC5WmupuysGzbKge/cMwIWndOuoPffvz2lzJ/D21lrWba8jaaoXMyf6\nZn8S4kiRTABLgBUAWuu1Sqn5PQtCdwdTgUeUUvnAb7TWjw60sZycNNxuc6CPDMjrzRzyurEqno85\nUFvy8aOf7MhcQwz25/ftK07kW/esonvPDFKy3hxUUhrO7yqef8/9kWMeGZFMAFlA70c6AaWUW2vt\nB9KBXwH3ACbwmlJqndb6w/421jSMvtajPYdoNIjnYw52peCrnApmN57i7RHbz2B/fqmmwYWLSli+\neje+CkVS6daI7atHPP+e+yPHPPh1+xPJRuBmoPeeXaGTP0A78EutdbvWugVYid1WIMQx+fZNh6Ab\nT7HG8HQ7Hc4nXLDQHosQqCsi2CaNASK6RTIBrAEuAAi1AWzqtWwasEYpZSqlPNiPi9ZHMBYRJwLN\nOQSbCnBlNDnS6+dYPG5X6K7EwLdPYUV2PJoQwxLJR0BPA2crpd7CnnDvaqXUzUC51voZpdQfgLWA\nD/i91npLBGMRccCyQlf/gKd45Iq8jTRzTCOu7DqCB/MIHszDzJHxASI6RSwBaK2DwA1HvL291/K7\ngLsitX8RfwKN47Hax2CO248rI7p7DHuKNF0Hc/HtU7jG1Ee8NIUQQyEDwURM6PIF8FdOBSMQkSqf\nI82V2oaZV4HVlR6ajlKI6CMJQMSEl96rwOpOxV2wF1dyp9PhhMUzoRxMH76qKVh+mX5bRB9JACLq\ntXb4WPHOXnB34R6/0+lwwmZ4fHa8AQ/+mlKnwxHiKJIARNR76b19dHQFcBfudnzE72C58/aBpwt/\nTSmWz+N0OEJ8giQAEdVa2rt5eV0lY9KT7JNpjDHMIJ7CXRB0468pczocIT5BEoCIaive2UdXd4AL\nFpVgmLE56YqZVwGeTvy1xVi+JKfDEeIwSQAiah1q6+bV9ZVkZySx9ITxToczZIYriGf8Tgi68e2f\n5HQ4QhwmCUBErRfX7qXbF+SixaV4hlEIMBqY3kqMpA4CdUVY3clOhyMEIAlARKnm9m5e/6CKnMxk\nPjU7dq/+exguy+4RZJnSI0hEDUkAIiq9uq6Sbn+Q804uxuOOj/+mZm6V3RZQVyTjAkRUOOZfllLq\nO0qpgtEIRgiAji4/r75fSUaqh1PnxP7Vfw/DZeEu2GP3CKotcTocIcK6A0gFVimlnldKXRaq3ilE\nxKzasJ/2Lj9nzZ9Isie2n/0fyZ1XAWY3/toSrEB8HZuIPcdMAFrrH2mtFfBj4HRgo1LqfqXUCRGP\nTiQcnz/IS+/tIznJ5Mx5E50OZ8QZZgB3/j7wJxGoj7/jE7ElrIerSqk0oAyYBASBA8AvlVI/jmBs\nIgG9vaWGg63dLD1hPOkp8Xmz6c7fCy6/PTo4GKU1rUVCCKcN4E/ATmApcIfWeqbW+lbgHOD6yIYn\nEknQsnhx7V7cpsE5C+K3gqbh8WF6K7G6Uwk0Fjodjkhg4dwBvApM0Vpfo7VeDaCUStJadwHHRzQ6\nkVA2ljdQ29TBwhkF5GTGd195d8EewMJfWyqzhgnHhJMAvq61but5oZRyAe8DaK1rIhWYSDwvvVsB\nwDkLihyOJPJcyZ2YY2uw2rMItox1OhyRoPrtjKyUWon92AelVO8iLH7gmciGJRLN3poWdMVBZpSN\nZaI3w+lwRoW7YA+BA4X4q0udDkUkqH4TgNb6DACl1C+11t8avZBEInrpPbvSZyJc/fdwZRzCldFE\n8FAe1Y1tFI5LdzokkWAGugO4SGv9HLBeKfXlI5drrX8f0chEwmhq6eLdbXWMz01nZlliPQ5xF+yh\nuzyHl9+r4MvnTXc6HJFgBmoDWBD6uhS7/3/vf0sjGpVIKCvXVxIIWpw9fyKGkVjdIl05tRjJ7azZ\nXENLe7fT4YgEM9AjoNtCX6/ueU8plQUUaa23jEJsIgF0+QK8/kEVGakeFs1IvIojhmGPC/DtO443\nNu7nwkWlTockEkg44wCuUUo9qpTyAluBJ5VSd0Q+NJEI3tlaS1unn6Vzx5MUZ2UfwmXmVpHsMXnt\ngyoCwdic9EbEpnC6gX4T+HfgCmA5MAs4L5JBicRgWRavrKvEZRgsPWGC0+E4xnD7WTyrgAPNXXyw\no8HpcEQCCasUhNb6AHAB8LzW2o9dIE6IYfmo8hCV9a2cqLyMzUpxOhxHnXmiXRfo1fcrHY5EJJJw\nipJvUUo9h10H6BWl1OPAusiGJRLBK6GT3ZknRufV/40rbxm1fd354e24suajK+Abz/0IV1prWOs9\nfvlDEY5MxLNwEsDXgMXAZq11t1LqD8ALx1opNGL4QWAO0AVcq7Uu7+MzzwPLtdYPDzZ4EbuaWrpY\nr+uZ6M1gWlG20+FEBXf+Xrqbc/HXlpBUJv0sROSF8wgoA/u5/41KqVuBucD3w1hvGZCitV4EfBe4\nu4/P3AHkhBmriCOvfVBF0LI4KwG7fvbHlV2PkdROoLFQZgwToyKc/2VPAIeAzcBgylYtAVYAaK3X\nKqXm916olLoUu7T0inA2lpOThnsYE4N7vZlDXjdWResx+/wBVn9YTUaqhwtPnUxKkpzswO4Saubv\nw18xHX/9RDyFe8JaL1p/z5EkxzwywvnLK9Banz2EbWdhJ44eAaWUW2vtV0rNBL4IXArcGs7Gmpra\nhxCCzevNpL6+Zcjrx6JoPua1W2s42NrFuScV0XKog+iM0hnu3Cr8lVMJ1BXhLthDODdH0fp7jpRo\n/r8dKcM55oESRziPgD5QSs0ewn6bgd57doV6EAF8GZgArAS+CtyslJKupQnitfVVACydG52Nv04y\nPD7McTVYXekEm8c5HY6Ic+HcAczETgK1QCdgAJbWetIx1lsDXAw8rpRaCGzqWaC1Pty9Qin1Q6BG\nax3WoyAR2yrrWvmo8hAzysaSn5PmdDhRyZ23j0DDBPx1RZhjGp0OR8SxcBLAZ4a47aeBs5VSb2En\njauVUjcD5VprKSedoF7bYF/9ny5X//0y0g9hpB0i2JSP1Z2MkdTldEgiTh0zAWit9yqlvgjMAO4E\nLg2nEqjWOgjccMTb2/v43A/DC1XEuo4uP29vriEnM5k5U+TxRn8MA9x5Ffj2zMRfV4RnYvmxVxJi\nCMKpBfQT7FHAn8VOGFcrpfrq0inEgNZuraWzO8BpJ4zHdIU1CD1hmeOqwfThr58oE8eLiAnnr/Bc\n4EtAp9a6GTgbOD+iUYm4Y1kWr62vwnQZnDpnvNPhRD3DDGDmVoEvheDBPKfDEXEqnDaAnvKEPWMA\nknu9J0RYdu5vprK+lfnKS3ZGfE/4PlLceRUEakvtxuCxtX1+5vN/+8aQtv3AGT8bTmgiToRzB/A4\n8DdgrFLq34A3gD9HNCoRd17/QLp+DpYrtQ1X5gGCzbkEO6XHlBh5x0wAWuufAr/BHhFcDNymtf6f\nSAcm4kdbp4/3tteRl5PK9BKp/DEYZl4FAIG6iQ5HIuLRMR8BhUbtZgKvA1u11rsjHZSIL29tqsHn\nD3LaCeNxSd2fQTFzavC5j8PfMBH3xI8wXIOpxiLEwAaaFD4PeBJ7INhH2G0AKtSv/0qt9cHRCVHE\nMsuyeH1DFW7T4JRZhU6HE3MMl2WXh6gpI9BUgHtctdMhiTgy0COgXwGrgXyt9cla64VAPvAhcO9o\nBCdi346Kg1Q3tjNP5ZGVluR0ODHp48dARQ5HIuLNQI+AZmutL+/9Rmg+gO8BGyIblogXqzbsB2Dp\nCdL1c6hcKe24shrsxuCOdFypbU6HJOLEQHcAnX29qbW2kG6gIgwt7d2s03UUjkuTSV+GyR26C/DL\nXYAYQQMlgIFam6QlShzTW5tr8AcsTpszXiZ9GSZXdh24uwg0jMcKyihqMTIGegQ0Qym1q4/3DUBa\n88SALMti1Yb9uE0Xi6Xxd9gMl4XbW4W/ehKBA/m4c6UxWAzfQAlg2qhFIeLOjoqD1BxoZ+GMfDJS\nPU6HExdMb6WdAOqLJAGIEdFvAtBa7x3NQER86Wn8PU3q/owYaQwWI00mYxUjrrXDd8zG3xtX3tLn\n+2Jg7rwKuptz8dcVkVRyVHV1IQZFWpPEiHtrU7U0/kaINAaLkRTOfAAvKKUuU0rJg1xxTPbI3/24\nTUMafyOgpzGYQBKBA/lOhyNiXDiXED8BzgM+Uko9oJRaEOGYRAzrafydr/Kk8TdCTG8lAIF6GRMg\nhiecKSHfAN5QSqUClwJPKaWagV8DD2mtZcJScdiqjaHGXxn5GzF2Y3AjweZx0hgshiWsh4hKqaXA\n/cD/ACuAf8WuCySTu4vDWjt8rNteT8FYGfkbaaY3VB+oXspEi6ELpxz0XmAX8Bhwk9a6I/T+68C6\niEYnYord+GuXfZbG38gyc2rxubvwN0wIlYmW6ixi8MLpBnqh1npz7zeUUgu11muBEyMTlog1lmWx\namOo8XdmgdPhxL2Py0RPItCUL2WixZAMNB/AKYAJ/FopdQ12CQgAD/AQMlJY9PJR5SGqG9tZeHw+\nmVL2eVSYeZV2AqibKAlADMlAdwBnA6dh1/35Ua/3/cD/RjIoEXte32DP+SuNv6Pnk43BabhS250O\nScSYgUpB/BBAKfUlrfUfRi0iEXN6Gn/zpfF31JneCoLN4wjUF+Eq1k6HI2LMQI+AfhhKAmcopU4/\ncrnW+msDbVgp5QIeBOYAXcC1WuvyXstvBL6KXVr651rrx4dyAGJkDLU0wwNn/Ozjxl8Z+TvqpDFY\nDMdAj4DeD319fYjbXgakaK0XKaUWAncDlwAopXKBbwBzgRRgq1LqidBkMyKG9B75e8osafwdbVIm\nWgzHQOMANiqlioHX+vl3LEuwxwwQ6jE0v2eB1roBOEFr7QMKgE45+cem3iN/pfHXGTIyWAzVQHcA\nq7Afz/R1T28Bk46x7SzgUK/XAaWUW2vtB9Ba+5VSNwG3A/cdK9CcnDTcbvNYH+uX15s55HVj1Wgc\n89tb6wC4ZOmUhPwZR4OhlImO9d9VrMc/FJE45oEagcuGue1moHfErp6Tf6993K+UegR4USl1uta6\n3zuLpqah93DwejOpr28Z8vqxaDSO2fJ5WPPhfgrHpZGXmZRwP+NoMtgy0bH8u5K/58Gv259jNgIr\npR7ta/mxGoGBNcDFwOOhNoBNvbatgB8DnwN82I3E0noVY/wNEwgELZaeMEEafx3myq4DTxeBhglY\nRTukMViEJZxG4FVD3PbTwNlKqbewHyNdrZS6GSjXWj+jlNoIvI39OOlFrfVQ9yMcYFkQqCvC43ax\nWBp/HWePDK7EXz2ZwIEC3Ln7nQ5JxICBHgE9G/r6O6VUHnAy9tX6u1rrA8fasNY6CNxwxNvbey2/\nHfv5v4hBweaxWF3pLJiZR3qKlH2OBmZeaM7guiJJACIs4UwIcxmwAfgKcB2wQSl1XqQDE9HNX1cM\nwNITJjgciejhSu7ANaaBYGsOwfbEayQVgxdOMbgfAPO01tUASqkS7DLQKyIZmIheVncywaY8jLRm\nJk/Icjoc0Ys7bx/dh7z4a4tIKtvqdDgiyoUzH4APqOl5obXei10PSCQof/1EwIU7r0Iaf6OMK7se\nI6mDQON4rMDQu02LxDBQL6Avh77dDTyrlPod9on/CmDjKMQmopBlGfjrisDlxxwnz5mjjWGAmVeB\nv3IagYbxuPMrnA5JRLGBHgH11P9pDf27IPS6jb4Hh4kEEGzygi8FM28vhhlwOhzRB3duJf6qKfjr\nijHzKpCbNNGfgXoBXd3fstD8wCIB9TT+uvPkyjJaGUndmDm1BA4UEmzNxsw86HRIIkqFMyXk54Bb\ngQzsK38TSAXyIhuaiDbBzjSCzbm4Mg/gSmt1OhwxADNvH4EDhQRqiyUBiH6F0wj8M+DfgG3Aldhz\nA0vp5gTkr7Wv/k25+o96rswmjNQWAk0FWD4p0if6Fk4CaArV6FkLjAnNEbAoolGJqGMFTAINE8DT\nhZlTc+wVhKMMw+4SiuXCXzfR6XBElAonAXQopaZh3wEsVUolAWMiG5aINoHGQgh4cHsrMFxSuTsW\nmLn7weXHX1eMFZSWYHG0cBLAD4A7gOeAM4Fa7Do/IkFYFvhrS8AISuNvDDHMAKa3CnwpBA9Kk504\n2jEbgUNF2noKtS1QSuVorZsiG5aIJsGWHKyOTMyx1RhJXU6HIwbBnbePQG0J/toSzLG1Tocjokw4\nvYAmYk/YshToBl5RSn1ba10f4dhElAiEun6a+fscjkQMliu17ePJYtozpPeW+IRwHgE9CrwMlADT\nsMtEPxbJoET0sLqTCRzIx0htxpUhN36xyB1K3D29uIToEU4xOK/W+qFer3+hlPpKpAIS0cVfVwS4\ncOfvkxGlMcqVXfdxfaCiHRhuKeUlbOHcAbyrlPpCzwul1EXAusiFJKKFFXTZCcDsxhxX7XQ4YogM\nA8z8vRB0hwr5CWEbqBhckI8nhf+6Uuo3QAB7RHATcO2oRCgcE2gsBH8y7sJdUvcnxrm9dn2gQG0J\n7oK9TocjosRAtYDCuTsQccqywF9TAgQx86TxN9YZbj9mbhWBuhKCTdIlVNjC6QWUBtyGPQbADawE\n/ktr3Rbh2ISDgi1jsTqyMMdW40rudDocMQLc+fsI1JWEErsQ4bUB3A+kA1/DnhYyCXg4kkEJ5/Wc\nJNwFe5wNRIwYV2obrjH1BFvHsqem2elwRBQIpxfQPK31nF6vb1JKyVxzcSzYmUrwYB5G+kFcGYec\nDkeMIHf+XroPeXn5vUq+fvHxTocjHBbOHYBLKZXd8yL0vfQji2P+2hLAkMbCOOQa04CR0sq722pp\napFR3YkunARwD3ZX0LuVUncD7wH3RjYs4RTL7yZQPxEjqUOqfsYhw7Af6wWCFq+8L3WdEl04CeBZ\n4LPALmAP8Fmt9aORDEo4x19XDEE3Zv5eqfoZp8zc/WSlJ/H6B/vp6JKb+UQWThvAm1rr44DNkQ5G\nOMsKGvbjH9MnVT/jmOEKcua8iTz9xi7e3Lifc06SEhGJKpwEsFEp9SXgXaCj502ttXQOjzOBxvHg\nS8ZdIAO/4t3pcyfw/Nt7eHldBWfMm4jblGE/iSicBHBy6F9vFjBpoJWUUi7gQWAO0AVcq7Uu77X8\n20BPiYkXtNa3hxu0GHmWBf7qMrvmvzT+xr2MVA+fmj2eV9+vZN32OhbOKHA6JOGAcOYDKBvitpcB\nKVrrRUqphcDdwCUASqlJ2PMLnwwEgdVKqae11h8OcV9imIIHvVidGZi5VYOq+X/jylsiGJWIpLMX\nFLFyfSUr3tnHycfnY0i1v4QzUC2g8diDwKYCq4H/1FofHMS2lwArALTWa5VS83stqwDO01oHQvvy\nAAMON83JScPtNgex+0/yejOHvG6sGswx+6vtPO8u2B2pcEQU8Xoz8XozOWX2eFZv3E9lUycnqtgp\nESF/zyNjoDuAx7Br/z8CXA78Arh6ENvOAnqPIgoopdxaa7/W2gc0KKUM4C7gA631joE21tTUPohd\nf5LXm0l9fcuQ149FgznmHRUHCbaOxTWmTiYMSRA9/zfOnDuB1Rv38+cXt1E0NtXhqMIjf8+DX7c/\nAyWACVrdaL2NAAAXbUlEQVTrcwGUUq8CGwa532ag955dWuvDfc6UUinYk820AN8c5LbFCHrurT0A\nuMfvcjYQMepKCjKZNWkcm3Y18lHlQaZOzD72SiJuDNT0393zTeiKvXuAz/ZlDXABQKgNYFPPgtCV\n/3Jgo9b6+p5HQWL07a5uZvPuA7gyGzEzB/OET8SLixbbdZ+ef1sa/xNNOL2Aegx2VNDTwNlKqbew\n5xS4Wil1M1AOmMBpQLJS6vzQ5/9Ta/32IPchhumF0B+9XP0nrqkTs5lWlM2HOxvZV9tCcX7iPV9P\nVAMlgBlKqd5nhQmh1wZgaa0H7AaqtQ4CNxzx9vZe36cMKlIx4qoa2nh/Rz1lhZlUZzU6HY5w0EWL\nSrin4iDPv72Xbyyb6XQ4YpQMlACmjVoUwhEvvL0HgIsWlfLrSkdDEQ6bUTaWkvxM1m2vo7qxjcJx\n6U6HJEbBQDOCyQPBOFbd2MbarbVM8KYzZ2ouSAJIaIZhcPEppdz/9008s2YP1396htMhiVEg478T\n1LNv7cGyYNmSMlwyAEgAc6fmUpyfwbtba6lqkAn/EoEkgAS0v6GNd7bUUpyXwdxpXqfDEVHCMAwu\nWVKGBTyzWgYEJgJJAAnomTW7sYBL5OpfHOGEKbmUFNhtAZX1Migw3kkCSDCVda28t62OkvxMTpia\n63Q4IsrIXUBikQSQYJaHrv6XfapMin+JPs2ZPI6ywkzW6Xr21SZWyYVEIwkggeyubuZ9XU9ZYRaz\nJ49zOhwRpQzD4DOn2sN8nly10+FoRCRJAkgQlmXxxGv2dAyXLZ0sV/9iQDNKx3JcSQ6bdx1g294m\np8MRESIJIEFs3n2A7fsOMnvyOKaX5DgdjohyhmFw6dLJADz5ejmWJfNDxyNJAAkgaFk88dpODOBz\np012OhwRI8oKs1gwPY/d1S28r+udDkdEgCSABLB2Sw2V9a0smllAUV6G0+GIGPLZUydhugyeWrUT\nfyDodDhihEkCiHPdvgBPv7Ebt2mw7FNDnd1TJKr8sWmcOmc8tU0drNqw3+lwxAiTBBDn/vnuPhqb\nOzlrfhG5Y2JjxicRXS5ZUkZqssk/3txFa4fP6XDECJIEEMcONHfy/Nq9ZKV5uHhxqdPhiBiVlZ7E\nxYvLaOv0s/xNGRwWTyQBxLEnV+2k2xfkc6dNJjV5MHP/CPFJZ82fSH5OKq99UEWVlIiIG0asdO+q\nr28ZcqCJNon0jStvIdCSTfe2hRhph0ie8TbS7V8M19eLb+G+Jz9kRmkON19+gmNjSRLt7xmGPSl8\nv78ouQOIQ5YFvn3HAeAp2SYnfzEi5kwex8yysWzZ08QHHzU4HY4YAZIA4lCgrgirbQzm2P0y0bsY\nMYZhcMVZUzFdBn9+ZQed3X6nQxLDJAkgzjS1dOGrnAamD0/x9mOvIMQgFI5L5/yFJRxo7uIf0iAc\n8yQBxJm/vPoRBDx4inZgJHU7HY6IQxctKiEvO5VX1lWytyaxnsXHG0kAcWRjeQPrttfhymjC9FY4\nHY6IU0keky+dqwhaFr//53aCwdjoSCKOJgkgTnR2+/njSzswXQae0i3S8CsiakbZWE4+Pp/d1S28\n+n6l0+GIIZIEECeeeG0njc2dnHtSMa406actIu8LZ04lI9XDU6t2UnOg3elwxBBIAogDm3c38toH\nVUzITeeSJVLvR4yOMelJXHXONLr9QX7z/FZ5FBSDIjY8VCnlAh4E5gBdwLVa6/IjPuMF1gCztdad\nkYolnrV3+njshe2YLoNrLzoej1tyuhg9Jx2Xz/u6nve21/HSexWcd3Kx0yGJQYjk2WIZkKK1XgR8\nF7i790Kl1LnAS0BBBGOIe3955SOaWrq4eHEpJQWZTocjEtBV50wjK83D39/YRVVDm9PhiEGIZAJY\nAqwA0FqvBeYfsTwInAUciGAMce297XWs2VxDSUEmFywqcTockaAy05L48nnT8QeC/O/yLXT7Ak6H\nJMIUyQphWcChXq8DSim31toPoLV+GUApFdbGcnLScLvNIQfj9cbX1XF1Qxu/W7GdlCST//jyAgrz\n4+v4RPQZ6G/oXG8m5dUtrHh7D8vf3ss3PzfHsVjiVSSOOZIJoBnoHbGr5+Q/FE1NQ+9lEG/Fo3z+\nIP/zx/dp7/Rz7UXHkeIiro5PRKdj/R9btriEzeX1vPjWHsryMpg/PS8iccTb33M4hlkMrt9lkXwE\ntAa4AEAptRDYFMF9JZQnX9/J3poWTplVwOKZhU6HIwRgDxC74ZKZJHlcPPbiduoPdjgdkjiGSCaA\np4FOpdRbwC+AbyulblZKfTqC+4x767bX8fK6CgrHpXHV2eE9PhNitIzPTefKs6fR0eXngac30SXt\nAVEtYo+AtNZB4IYj3j6qOpnWujRSMcSbfbUt/Pr5rSR7TL6xbCbJSUNvExEiUpbMKmRnVTNvbNzP\nYy9s4/pPz3Bs7gAxMOk0HiOa27v51VOb6PYFufai45nozXA6JCH6ZBgGV50zjSkTxvDutjpWvLPP\n6ZBEPyQBxAB/IMiDT2+msbmTZUvKmKe8TockxIDcposbPzOTnMxknnx9JxvLZQKZaCQJIMoFLYvH\nXtjGjoqDzFdeLjql1OmQhAjLmIxkbvrsLNxuFw8v38Lu6manQxJHkAQQ5Z58bSdvb6ll0vgsrrnw\neFzyLFXEkLLCLK7/9Ay6/QHufWIjtcPozi1GniSAKPbPd/ex4t19FIxN41uXzpZGXxGTTpzm5apz\nFC3tPu752wYOtclERdFCEkCUevPD/fxtZTnZGUncfPkcMtOSnA5JiCE7fe4ELl5cSv3BTu7+6wZa\n2iUJRINIjgSOSzeuvCXi+/DXT8C3eyaYPjpKV3Pb+89EfJ9CRNqyT5XR2unjtfVV/PyvG/jOFXPJ\nSPX0+/mh/K09cMbPhhNiwpE7gCjjr5uIb/cscPtIPu49mdxFxA3DMLjy7GmcPncCFXWt3PWXD2jt\n8DkdVkKTBBBF/LXF+PbMBHc3ydPfxZWWWPVORPxzGQZXnvNxEvjpn9fT1NLldFgJSxJAFLAs8FVM\nw7f3ePB0hU7+cuUv4lNPEjhr3kSq6tu48w/rZB4Bh0gCcJgVNPDtmoW/ehJGShvJx6+Vk7+Iey7D\n4IqzpnLp0skcaO7ix394nx0VB50OK+FIAnCQ1Z1M9/YFBBon4MpoIvm4tbiSpYKiSAyGYXDBwhKu\nveg4unwB7vrLB7y+ocrpsBKK9AJySKAlm+7yE8CXgjm2Gk/ZJgwz6HRYQoy6xTMLyc5I5uHlW/j9\nCs2e6maulEq3o0ISwCizLAjUluCrUGCBu2g77oI9yABfkciOLx3LrV+Zz/1Pb+KNjdXsq20l6E3D\nlSojhyNJHgGNomBXCt16Ab59x4HpI2n6OjyFcvIXAiA3O5XvXTWPU2YVsKemha4ti/HXTcSynI4s\nfkkCGAWWBf768XRtPoVg8zhc2XWkzFqDmXXA6dCEiCpJHpNrLjyeGy6ZAYaFb89MunecSLArxenQ\n4pIkgAgLtmfQvf0kfLtng2XgKd1M0tT1GB4ZCi9Ef046Lp+UWatxZTUQPJRH16Yl+KpLsYJyuzyS\npA0gQiyfB9/+yQTqisFy4cquxVOyDVdyp9OhCTEkQy2DMtTyDEZSF0lqHYGG8fgqpuOvmE6gYQKe\nIo1rTENUPDod7Z/JSJMEMMKsgIm/tgR/dRkEPBhJ7XhKtmHm1DsdmhAxxzDA7d2PmV2Pr3IagfqJ\ndO+YjyvzgJ0IMg45HWJMkwQwQiyfB39dMf7aEvAngbsbT/E2zLwKDJd07xRiOAyPj6SyLQTz9+Kr\nnEbwYB5dWxfhymrAXbgbV1ZjVNwRxBpJAMMUbM/AX1dEoGECBN1gduMeX467cDeGGXA6PCHiiiut\nleRp6wm05OCvmkywOZfu5lyMtEO48/fR1R2QeTMGQRLAEFh+N4GmfPx1RVht2QAYSR24Cz7C9FbK\niV+ICDMzmzCnryPYloW/uozAgQJ8u2dx8wOrWTSjgEUzCpg0PgtDbgsGJAkgTK0dPjbtbKRrx4kE\nD+WC5QIsXGPqcOdV4squxzCkw7IQo8mV3kzSlI0EuzSB+okkNx/PyvVVrFxfRe6YFE46Lp+503Ip\nK8yS6VT7IAmgH/5AkN3VzWzf28SHuxrZtb85NCAlDyOtGXNsNea4aunVI0QUcCV34ppYzl1Lr2Xz\nrgO8u62W9R818MLavbywdi9ZaR5mTRrH8aVjUcXZjM2ScQUgCQAAy7I42NrN7upmdlc3s2t/Mzur\nDtHttxtvDQOmTBjD7MnjeO7A73ClSulaIaKR6XIxZ0ouc6bk0uULsHnXAT7c2cCHOxtZs7mGNZtr\nAMjLTmXyhDFMGp/FpPFZTPSm43EnXttBQiWAQDDIgeYuapvaqT3QQU1jO5X1rVTWt9LW6f/EZyd4\n05lelIMqzmZ6Sc7hqeteWCknfyFiQbLHZJ7yMk95CVoWFbWtbN/XxPa9TeyoPMjbW2p4e4udEFyG\nQf7YVIryMigcl07+2FQKxqaRl51KWkr/01bGuoglAKWUC3gQmAN0Addqrct7Lf86cD3gB+7QWj8X\niTj21rRw75MfUlXXSlNLF8EjCosYgDcnFVWcQ2lBJmXjsygryIzrX7oQicZlGJQUZFJSkMm5JxUT\ntCxqGtvtO/7qZirr7AvB6saji8+lJrvJHZPC2MxksjOTyc5IZkxGEpmpSQRasjHcPgx3N5h+DFds\ntQNG8g5gGZCitV6klFoI3A1cAqCUKgD+FZgPpACrlVIva61HfG646gNtbN19gMw0D5MmZJGblUL+\n2DTyc1LJH5vG+HHp0m1MiATjMgzG56YzPjedU2YVAvaj4MbmTmoOfPyEoP5QB42HOqlr6qCirq+J\nmhYesWG/nQjMnq8BcAUwXPbXnu///sYuPG4XHtOF2zRwmy5M08B0GZguV+irgctl4HG7WJyTFpGf\nQyQTwBJgBYDWeq1San6vZScBa0In/C6lVDkwG3hvpINYeHwBF506hYYGmWVLCNE/wzDIHZNK7phU\nZpZ9cpllWXR0BWhq7eJgSxfNbd20tHfz5LaXsPweLL8H/B6sgAcCpv26Kw3L6rvc2nM1ewYV25eb\nu1g6u3CIR9a/SCaALKD3OO2AUsqttfb3sawFGDPQxrzezGH14fJ6M4ez+mGPX/7QiGxHCDGwaPxb\nKzni9ZUXznAkjpESyWqgzUDvs64rdPLva1kmIBOCCiHEKIpkAlgDXAAQagPY1GvZu8CnlFIpSqkx\nwHHA5gjGIoQQ4giGFaHpdnr1ApqN3dnmauyEUK61fibUC+g67CT0P1rrpyISiBBCiD5FLAEIIYSI\nbjIjmBBCJChJAEIIkaAkAQghRIKKu1pASqlk4DFgEnZ30xsBC3gYSMIuS/EFrXWjY0GOsH6O+f96\nfWQ68Fut9XcdCC8i+jnmEuAn2OVFXtFa/8C5CEdeP8dcBvwUaANWaK3vcC7CkaWUOhn4qdZ6qVJq\nCvBb7L/lzcCNWuugUuo24ELs3/m/aa3fdSzgYQrneEOfmwI8rbWeNdx9xuMdwNeBVq31QuBfgPuB\nR4AfaK1PxU4E0xyMLxKOOmat9VKt9VLga0AlEDcnhpC+fs93AV8GFgFLlVLD/gOJMkce84PAr4HP\naa2XANOVUkucDHCkKKVuwT62nrrN92D/DX8Ku1fhJUqpE4HTgJOBLwAPOBHrSAjneEOf+xLwV8A7\nEvuNxwRwPPAigNZaAycCecDFSqnXsU8OMXuV0I8jj/m4XsvuBf5Dax1vtTD6OuYPgLGAB/sPKd6m\nZjvymE8BmrTWu0LL12CXYIkHO4HP9no9D1gV+v5F4CzsY31Ja21prfcBbqXUiJwYHRDO8QI0YSe9\nERGPCWADcJFSyggNQMsFZgCvAKcDOcBXHIwvEo485glKKVMpNRvI0lq/6nB8kXDUMQNbgOeAbUAF\nsN3B+CLhyGNOBtKUUtOVUib2OJt0RyMcIaFxQb5ebxla654+6z2lYwZdUiZahXm8aK2f01qPWE36\neEwAj2I/H30T+AywDmjRWr8W+oE+h12FNJ4ceczva60DwFV8si0gnhx5zDuB/wBmaK0nAx8B/8+5\n8CLiqN8z8CXgIeB5QAMNjkUXWcFe3/eUjonnkjJ9He+Ii8cEsAB4NfRM9AnsE8MOpdSnQstPxb5S\njCdHHnPPI4EzCVVkjUNHHvMWoDX0D6Aa+24vnvT1ez439O98YDL2nW48+kAptTT0/fnYSXANcK5S\nyqWUKsauNxYvCbCv4x1xcdcLCPvK77+VUt/HzprXYDeYPKCUcgO7sa8U40lfxwxQEE+9nY7Q1zGf\nDLyklOoMvfdV58KLiL6O+ULsNq0O4E9a63i7uOnx/4D/U0olYT/ie1JrHVBKvQm8jX0xe6OTAY6w\no443EjuRUhBCCJGg4vERkBBCiDBIAhBCiAQlCUAIIRKUJAAhhEhQkgCEECJBxWM3UBGlQv2anwPK\nseubJAF/1FrfOYhtlAKva61Lw/js94HLQi/nABtD3z8xmH2OJqXUncDzWuu3IrDtzwIF2LWxfHz8\n8+hxLXZxuXFa64dHev8i+kgCEKNtXahIHUqpDGCbUupprfXWkd5R6CR/Z2hfltb6hJHeRwScRqje\nz0hSSqUC/449EBIg0M/PY51S6nWl1BNxPIZEhEgCEE5Kwy7YdgggVN/ml9iF3BqA67XW5UqpucBv\nQutsDH02E3tQ3yStdXPozuB5rfWMcHYcGhT4c+wTohv4tdb6PqXUWcB3sAvKlQF/wx5k9Wns0rzn\nYxfkqgT+CcwNxX+l1npfqKTvPUAqUA9cp7Xeq5RaDdRh16W6FDgD+CKQgV3K+PPAYuAE4DGl1CXY\nV+rf1VqvDpUAXqG1nqKU+iN2bZgp2AOGGvva5xGH/CVgpdbaHzr2gfwD+Cbw3+H8LEXskjYAMdrm\nK6U2KKU+BPYArwP7QyMe/wrcpLWeg122+y+hdX4P3KK1PpFQmQutdQt2/ZtLQ5/5cuhz4boB6A5t\ncwFwmVJqUWjZydgnzFnAt4BKrfU87BGZnw99Jh+7EuVs4O/AvaF6/f8HXB7a7n3A//ba53qttcIu\nVHcRcFooYT0P3KC1fgy74NvVYdwR1WqtjwNePcY+e3waeKPXazP0e+j5d1evZW+EPi/inNwBiNF2\n5COg54DvAs9ilzZ+D0Br/YRS6hGlVAkwXmvdU+Pmt3xc6uJR4Iehr1/EvqoO11nATKXUOaHXmdgn\n/F3Ah1rrqlCMjdgnWYC9fFxfqE1r/afQ978DbsMuST0JeE4pBXY7R2qvfb4TOraDobruX1RKTQPO\nA94bROyHtxXGPntMxb5r6dHfIyCwj3PqIOMRMUgSgHCM1rpVKfU0cDb2VfCRDOzibkav9/y9vn8D\nu/T1Z4HdWuv9g9i9CdystX4GQCmVG9rXEqD7iM/6OVrvuQZcoc+YwI7QlTihEs15vT7XEXq/FDup\n/Ap4AfvR0PF97MPi42P3HLGso9dxDLTPHsF+jqMvPj5ZjVLEKXkEJBwTOlktBdZjlzIep5RaEFr2\neWBvqCFyr1LqwtBqX+xZP1Te+3fYjz1+O8jdrwS+rpRyh9oT3mZwZcKzlFLnh76/GrvhditQoJRa\nHHr/OuAPfax7ErBda30vdiG3C7BP5GCfpHsuzBqw2wwAlvUTR7j73Ik9ZWY4yrB7aok4JwlAjLae\nNoAN2BO2tGPPg9oFXA7cr5TaDNwUeg32vAa3KaU+wC553NtfsRuT/zHIOB4A9mE3Kq8D/ldrvXoQ\n6weBK5VSm7AnGrpZa92B3Ubwy1AbxxXY0zge6UUgRSm1BTvx7MI+6fYs+3WoMfknwLeUUus5+g4A\ngEHs89lQnOE4HVge5mdFDJNqoCJmKaVc2I2507XW/zqK+3UDnVrrmHmEqpRKw55icJHWesBHQUqp\nt4BPx1FtfdEPuQMQsezv2IOXpLviMWit24GfAdcP9Dml1BeAP8vJPzHIHYAQQiQouQMQQogEJQlA\nCCESlCQAIYRIUJIAhBAiQUkCEEKIBPX/AQ47O2HwpP3qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108a602e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Fit the temperature data to a normal distribution, obtain its mean and standard deviation\n",
    "n_mu, n_sigma = scs.norm.fit(df.temperature)\n",
    "print('Mean of fit normal distribution: %.4f' % n_mu)\n",
    "print('Standard deviation of fit normal distribution: %.4f' % n_sigma)\n",
    "\n",
    "# Plot the distribution with fit values\n",
    "x = np.linspace(96, 101, 100)\n",
    "plt.plot(x, norm.pdf(x, loc=n_mu, scale=n_sigma), label='Fit Curve')\n",
    "plt.hist(df.temperature, bins=20, normed=True, label='Data')\n",
    "plt.xlabel('Body Temperature (F)')\n",
    "plt.ylabel('Probability Density')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is difficult to assess just by a visual inspection, so conduct a normality test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.70380143332\n",
      "p-value = 0.2587\n"
     ]
    }
   ],
   "source": [
    "# Normality test\n",
    "stats, pvalue = scs.normaltest(df.temperature)\n",
    "print(stats)\n",
    "print('p-value = %.4f' % pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sample distribution is could be distributed normally considering high p-value of ~0.26. In this situation, I would use one-sample test with t-statistic because: \n",
    "1. We only have one sample to test (temperature) and we want to estimate the population mean. Thus, use one-sample test\n",
    "2. Sample size is relatively large (n > 30) \n",
    "3. Exact population variance is unknown\n",
    "4. t-ststistic is a more conservative than z-statistic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now actually conduct the one-sample t-test to see if 98.6 is a good estimate of the population mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t-statistic = -5.4548\n",
      "p-value = 0.0000002411\n"
     ]
    }
   ],
   "source": [
    "# Test if the population mean for body temperature is 98.6F with t-statistic\n",
    "t_statistic, pvalue = scs.ttest_1samp(df.temperature, 98.6)\n",
    "print('t-statistic = %.4f' % t_statistic)\n",
    "print('p-value = %.10f' % pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "p-value is very small (~2.0e-7), thus it is quite likely that the real  population mean is <i>not</i> 98.6F.<br>\n",
    "Now try using the Z-statistic assuming normal distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z-statistic = 0.4803\n",
      "p-val = 0.3270\n"
     ]
    }
   ],
   "source": [
    "z = (98.6 - n_mu) /n_sigma \n",
    "pval_z = (1-  scs.norm.cdf(0.5 + z)) * 2\n",
    "print('z-statistic = %.4f' % z)\n",
    "print('p-val = %.4f' % pval_z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming a normally distributed sample distribution, the p-vale is very high and the value 96.8F is <i>not unlikely</i> to be the population mean. The sample distribution is not so well described by the normal distribution, so t-test should be employed as shown above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now calculate margin of error and interval for 95% confidene level  using t-statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Basic stats for the temperature data ---\n",
      "Sample Mean 98.2492\n",
      "Sample Standard Deviation 0.7332\n",
      "\n",
      "--- Parameters to obtain 95% Confidence Interval ---\n",
      "t-critical value: 1.9785\n",
      "Margin of Error 0.1272\n",
      "95% Confidence Interval 98.1220F ~ 98.3765F\n"
     ]
    }
   ],
   "source": [
    "# Set confidence level and calculate critical value for t statistic\n",
    "confidence = 0.95\n",
    "t_critical = scs.t.ppf(q = ((1 + confidence)/2), df=(len(df.temperature) - 1))\n",
    "\n",
    "# calculate mean and standard deviation of the temperature data\n",
    "mean = df.temperature.mean()\n",
    "stdev = df.temperature.std()\n",
    "\n",
    "# calculate margin of error\n",
    "sigma = stdev/math.sqrt(len(df.temperature))  # Standard deviation estimate\n",
    "m_of_e = t_critical * sigma\n",
    "\n",
    "# calculate confidence interval\n",
    "CI = (mean - m_of_e, mean + m_of_e) \n",
    "\n",
    "print('--- Basic stats for the temperature data ---')\n",
    "print('Sample Mean %.4f' % mean)\n",
    "print('Sample Standard Deviation %.4f' % stdev)\n",
    "print('\\n--- Parameters to obtain 95% Confidence Interval ---')\n",
    "print(\"t-critical value: %.4f\" % t_critical)\n",
    "print('Margin of Error %.4f' % m_of_e)\n",
    "print('%d%% Confidence Interval %.4fF ~ %.4fF' % (100*confidence, CI[0], CI[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the 95% confidence interval, if the body temperature is below 98.1F or above 98.4F, it could be considered abnormal (only 5% probability of having that temperature in a healthy/normal state). This is consistent with the above analysis that 98.6F is unlikely to be the real mean of the population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compare the temperatures between males and females."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>\n",
       "      <th></th>\n",
       "      <th colspan=\"8\" halign=\"left\">heart_rate</th>\n",
       "      <th colspan=\"8\" halign=\"left\">temperature</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gender</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>F</th>\n",
       "      <td>65.0</td>\n",
       "      <td>74.153846</td>\n",
       "      <td>8.105227</td>\n",
       "      <td>57.0</td>\n",
       "      <td>68.0</td>\n",
       "      <td>76.0</td>\n",
       "      <td>80.0</td>\n",
       "      <td>89.0</td>\n",
       "      <td>65.0</td>\n",
       "      <td>98.393846</td>\n",
       "      <td>0.743488</td>\n",
       "      <td>96.4</td>\n",
       "      <td>98.0</td>\n",
       "      <td>98.4</td>\n",
       "      <td>98.8</td>\n",
       "      <td>100.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M</th>\n",
       "      <td>65.0</td>\n",
       "      <td>73.369231</td>\n",
       "      <td>5.875184</td>\n",
       "      <td>58.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>73.0</td>\n",
       "      <td>78.0</td>\n",
       "      <td>86.0</td>\n",
       "      <td>65.0</td>\n",
       "      <td>98.104615</td>\n",
       "      <td>0.698756</td>\n",
       "      <td>96.3</td>\n",
       "      <td>97.6</td>\n",
       "      <td>98.1</td>\n",
       "      <td>98.6</td>\n",
       "      <td>99.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       heart_rate                                                     \\\n",
       "            count       mean       std   min   25%   50%   75%   max   \n",
       "gender                                                                 \n",
       "F            65.0  74.153846  8.105227  57.0  68.0  76.0  80.0  89.0   \n",
       "M            65.0  73.369231  5.875184  58.0  70.0  73.0  78.0  86.0   \n",
       "\n",
       "       temperature                                                      \n",
       "             count       mean       std   min   25%   50%   75%    max  \n",
       "gender                                                                  \n",
       "F             65.0  98.393846  0.743488  96.4  98.0  98.4  98.8  100.8  \n",
       "M             65.0  98.104615  0.698756  96.3  97.6  98.1  98.6   99.5  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check for gender difference\n",
    "df.groupby('gender').describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Slice the dataset into two\n",
    "temp_male = df[df.gender == 'M'].temperature\n",
    "temp_female = df[df.gender == 'F'].temperature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conduct the t-test assuming that male and femal temperature samples are independent from each other.\n",
    "Now check if the two samples are normally distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p-value for male: 0.6404\n",
      "p-value for female: 0.0937\n"
     ]
    }
   ],
   "source": [
    "stat_m, pval_m = scs.normaltest(temp_male)\n",
    "stat_f, pval_f = scs.normaltest(temp_female)\n",
    "print('p-value for male: %.4f' % pval_m)\n",
    "print('p-value for female: %.4f' % pval_f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the p-values, there is a 10% chance that the female samples are not drawn from normally distributed population. Thus, as the above analysis, use t-test to compare the means of male and female body temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Male Stats : Mean = 98.105F, Standard Deviation = 0.699\n",
      "Female Stats : Mean = 98.394F, Standard Deviation = 0.743\n",
      "\n",
      "\n",
      "t-test result: p-value = 0.0239\n"
     ]
    }
   ],
   "source": [
    "# calcuate means and standard deviations for each dataset\n",
    "male_mean = temp_male.mean()\n",
    "male_stdev = temp_male.std()\n",
    "female_mean = temp_female.mean()\n",
    "female_stdev = temp_female.std()\n",
    "\n",
    "# conduct t-test to see if the means are significantly different than expected by chance\n",
    "t_statistic2, pvalue2 = scs.ttest_ind(temp_male, temp_female)\n",
    "\n",
    "# print outputs\n",
    "print('Male Stats : Mean = %.3fF, Standard Deviation = %.3f' % (male_mean, male_stdev))\n",
    "print('Female Stats : Mean = %.3fF, Standard Deviation = %.3f' % (female_mean, female_stdev))\n",
    "print('\\n')\n",
    "print('t-test result: p-value = %.4f' % pvalue2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Low p-value of 0.02 indicates that the difference in mean between the males and the females are statistically significant (there's ~2% chance that they're difference by chance)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now calculate the 95% confidence interval for female mean to see if 98.6F falls within in the interval. I picked female population for this test because their mean is higher than that of males."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Parameters to obtain 95% Confidence Interval for mean of female body temperature ---\n",
      "t-critical value: 1.9785\n",
      "Margin of Error 0.1272\n",
      "95% Confidence Interval 98.2114F ~ 98.5763F\n"
     ]
    }
   ],
   "source": [
    "# calculate margin of error for female temperature\n",
    "sigma_f = female_stdev/math.sqrt(len(temp_female))  # Standard deviation estimate\n",
    "m_of_e_f = t_critical * sigma_f\n",
    "\n",
    "# calculate confidence interval\n",
    "CI = (female_mean - m_of_e_f, female_mean + m_of_e_f) \n",
    "\n",
    "print('--- Parameters to obtain 95% Confidence Interval for mean of female body temperature ---')\n",
    "print(\"t-critical value: %.4f\" % t_critical)\n",
    "print('Margin of Error %.4f' % m_of_e)\n",
    "print('%d%% Confidence Interval %.4fF ~ %.4fF' % (100*confidence, CI[0], CI[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on my analysis of the sample population, mean body temperature of 98.6F suggested by Wunderlich is statistically incorrect. The actual mean lies between 98.1F ~ 98.4F at 95% confidence level. Even for the females, who have the higher average body temeprature than the males, 98.6F is slightly above the normal range at 95% confidence level (98.2F1 ~ 98.58F). But there is a slight chance (~2.5%) that Wunderlich obtained the 98.6F value by chance if he had sampled only the female population! In any case, the real average body temeprature of the human population is quite likely below 98.6F based on this sample of 130 people."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}