regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regression\n",
    "====================\n",
    "\n",
    "Copyright 2015 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import first\n",
    "import thinkstats2\n",
    "import thinkplot\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load up the NSFG data again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9148, 244)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live, firsts, others = first.MakeFrames()\n",
    "live.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And select live, full-term births."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8001, 244)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live = live[live.prglngth>=37]\n",
    "live.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And drop rows with missing data (just for the variables we want)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7912, 244)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live = live.dropna(subset=['agepreg', 'totalwgt_lb'])\n",
    "live.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check a few rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>caseid</th>\n",
       "      <th>pregordr</th>\n",
       "      <th>howpreg_n</th>\n",
       "      <th>howpreg_p</th>\n",
       "      <th>moscurrp</th>\n",
       "      <th>nowprgdk</th>\n",
       "      <th>pregend1</th>\n",
       "      <th>pregend2</th>\n",
       "      <th>nbrnaliv</th>\n",
       "      <th>multbrth</th>\n",
       "      <th>...</th>\n",
       "      <th>laborfor_i</th>\n",
       "      <th>religion_i</th>\n",
       "      <th>metro_i</th>\n",
       "      <th>basewgt</th>\n",
       "      <th>adj_mod_basewgt</th>\n",
       "      <th>finalwgt</th>\n",
       "      <th>secu_p</th>\n",
       "      <th>sest</th>\n",
       "      <th>cmintvw</th>\n",
       "      <th>totalwgt_lb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3410.389399</td>\n",
       "      <td>3869.349602</td>\n",
       "      <td>6448.271112</td>\n",
       "      <td>2</td>\n",
       "      <td>9</td>\n",
       "      <td>NaN</td>\n",
       "      <td>8.8125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3410.389399</td>\n",
       "      <td>3869.349602</td>\n",
       "      <td>6448.271112</td>\n",
       "      <td>2</td>\n",
       "      <td>9</td>\n",
       "      <td>NaN</td>\n",
       "      <td>7.8750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>5</td>\n",
       "      <td>NaN</td>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7226.301740</td>\n",
       "      <td>8567.549110</td>\n",
       "      <td>12999.542264</td>\n",
       "      <td>2</td>\n",
       "      <td>12</td>\n",
       "      <td>NaN</td>\n",
       "      <td>9.1250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7226.301740</td>\n",
       "      <td>8567.549110</td>\n",
       "      <td>12999.542264</td>\n",
       "      <td>2</td>\n",
       "      <td>12</td>\n",
       "      <td>NaN</td>\n",
       "      <td>7.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>6</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7226.301740</td>\n",
       "      <td>8567.549110</td>\n",
       "      <td>12999.542264</td>\n",
       "      <td>2</td>\n",
       "      <td>12</td>\n",
       "      <td>NaN</td>\n",
       "      <td>6.1875</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 244 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   caseid  pregordr  howpreg_n  howpreg_p  moscurrp  nowprgdk  pregend1  \\\n",
       "0       1         1        NaN        NaN       NaN       NaN         6   \n",
       "1       1         2        NaN        NaN       NaN       NaN         6   \n",
       "2       2         1        NaN        NaN       NaN       NaN         5   \n",
       "3       2         2        NaN        NaN       NaN       NaN         6   \n",
       "4       2         3        NaN        NaN       NaN       NaN         6   \n",
       "\n",
       "   pregend2  nbrnaliv  multbrth     ...       laborfor_i  religion_i  metro_i  \\\n",
       "0       NaN         1       NaN     ...                0           0        0   \n",
       "1       NaN         1       NaN     ...                0           0        0   \n",
       "2       NaN         3         5     ...                0           0        0   \n",
       "3       NaN         1       NaN     ...                0           0        0   \n",
       "4       NaN         1       NaN     ...                0           0        0   \n",
       "\n",
       "       basewgt  adj_mod_basewgt      finalwgt  secu_p  sest  cmintvw  \\\n",
       "0  3410.389399      3869.349602   6448.271112       2     9      NaN   \n",
       "1  3410.389399      3869.349602   6448.271112       2     9      NaN   \n",
       "2  7226.301740      8567.549110  12999.542264       2    12      NaN   \n",
       "3  7226.301740      8567.549110  12999.542264       2    12      NaN   \n",
       "4  7226.301740      8567.549110  12999.542264       2    12      NaN   \n",
       "\n",
       "   totalwgt_lb  \n",
       "0       8.8125  \n",
       "1       7.8750  \n",
       "2       9.1250  \n",
       "3       7.0000  \n",
       "4       6.1875  \n",
       "\n",
       "[5 rows x 244 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And summarize a few variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>agepreg</th>\n",
       "      <th>totalwgt_lb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>7912.000000</td>\n",
       "      <td>7912.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>24.978917</td>\n",
       "      <td>7.482384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>5.523504</td>\n",
       "      <td>1.169065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>10.830000</td>\n",
       "      <td>0.125000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>20.580000</td>\n",
       "      <td>6.687500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>24.410000</td>\n",
       "      <td>7.437500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>28.910000</td>\n",
       "      <td>8.187500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>44.080000</td>\n",
       "      <td>15.437500</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           agepreg  totalwgt_lb\n",
       "count  7912.000000  7912.000000\n",
       "mean     24.978917     7.482384\n",
       "std       5.523504     1.169065\n",
       "min      10.830000     0.125000\n",
       "25%      20.580000     6.687500\n",
       "50%      24.410000     7.437500\n",
       "75%      28.910000     8.187500\n",
       "max      44.080000    15.437500"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live[['agepreg', 'totalwgt_lb']].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a scatterplot of age and birthweight, with parameters tuned to avoid saturation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXtsZHl233d+9a4iWSyyqvhoNtns5vQ8d7Wr2bUsWxLU\njhRBSCI7gALHzgN6JPonMCDDhh1ZQaTZ/BE4chTYcOIAiSzFQixFsSHIUZyHvUDakLAe7c7uzs7s\n9sz0DLub7zdZZLHej5s/Tn/m3OLwUWSz+Oi+B2g0Wayq+7u/e+95fs/3OM/zJJBAAgkkkBdPQpe9\ngEACCSSQQC5HAgMQSCCBBPKCSmAAAgkkkEBeUAkMQCCBBBLICyqBAQgkkEACeUElMACBBBJIIC+o\nRC57AYeJcy7ApgYSSCCBnEE8z3PdvvfKRgCe513bf7/6q7966WsI1n/563jR1h6s//L/nVaurAEI\nJJBAAgmktxIYgEACCSSQF1QCA9ADuXfv3mUv4ZkkWP/lyXVeu0iw/usm7ix5o16Lc867iusKJJBA\nArnK4pwT73koAgcSSCCBBNJbCQxAIIEEEsgLKoEBCCSQQAJ5QSUwAIEEEkggL6j0zAA4537TObfm\nnHv/kL/9dedc2zk33KvjBxJIIIEEcrz0MgL4LRH5yYMvOucmReTfFJG5Hh47kEACCSSQE6RnBsDz\nvD8SkZ1D/vTficjf7NVxA7k48TyRclmkVru8NTQaIqWSSKvV+2O1WnqsRqP3xwokkIuQCyWDc879\nBRFZ9DzvPee6hqoGckVlc9OU4eCgSF/fxR6/XhfZ2lJDVCyK5PMi4XBvjtVqiWxsiLTbIs6JZLMi\nsVhvjhVIIBclF2YAnHMpEfll0fTPpy8f9f633nrr05/v3bv3wnXoXXVptzs94VrtcgwA/YKsp1cG\noNHQY4joMev1wAAEcvly//59uX///pk/39NOYOfctIj8oed5n3fOfV5Evioi5ad/vikiSyLyA57n\nrR/4XNAJfA1kc1MVoYhIJiOSSl3s8RsNXYPniYRCIiMj+n8vpN3WCKDV0ggglxOJRntzrEACOauc\nthP4wgzAIX97LCJf8jxv+5C/BQbgGWVvT6TZVK88Hu/NMTxPpFpVr/uyvOFmU41QPN477x9ptTTS\nicVEIldykkYgL7pcGSoI59zvisjXRORl59yCc+7nDrwl0PA9kmJRZH9flfP2tqUuzlucE0kmLzcV\nEolo5NFr5S+ix0ilAuUfyPMjARnccyiFgqJzkNHRi1GQgQQSyOXKlYkAArk86eszhX9R3nEggQRy\n/SSIAJ5jabd7VxS9DuJ5mgKr10USCZGhocteUSCB9FaCCCCQT+VFVv4iIpWKFm09T3+uVi97RYEE\ncrXkBVcRgTzPcrDXMOg9DCSQTgnwDIE8t5JMavqHFFCv4LCBBHJdJagBBBJIIIE8JxLUAAIJJJBA\nAulKAgMQSCCBBPKCSmAAroHU60rtcJm0y4EEEsjzJ4EBuOLSbCrl8f6+YdoDCSSQQM5DAgNwxaXR\nMMpjzwuGkQQSSCDnJ4EBuOLiZ7kMhxXOGEgggQRyHhLAQK+B4PlHIkF3byCBBHK0nBYGGjSCXQNx\nLpg+FUgggZy/BP5kIIEEEsgLKoEBCCSQKy6eZ0CAQAI5TwlSQIEEcoWlUtEBPyIi6bTOeggkkPOS\nIAII5NpIo6F9ES+SFIsWARSLl72aQJ43CSKAQK6FMObSOZFMRpk+XwQJh83oBZPdAjlvCSKAQK6F\nMOPY80RKpctdy0XK0JCO9UwmRYaHL3s1gTxvEkQAgUijoRQT8bj2GlxFiUTME45GL3ctFymhkEY8\nlyWVio4WTaWCgTrPo/Q0AnDO/aZzbs05977vtb/jnPvAOfcd59zvO+cGe7mGQI6XRkNkc1Nkd1f/\nb7Uue0WHSy4n0t+vhdDB4I65ECkWRXZ29N7Y2rrs1QTSC+l1Cui3ROQnD7z2L0TkDc/zviAiD0Xk\nb/V4DYEcI/W6QQzb7atLNhcKqfLv77/slbw44mef9d8ngTw/0lMD4HneH4nIzoHX/qXnee2nv/6J\niNzs5RoCOVqaTQ3xSyV9uEOh3nQce57SWRcKx6N4ajX1OPf3z38NlyXttnrQhcLVja6OEn+hPZHo\nTQqo0dC92dsLDMxlyGVnfH9eRH73ktfwwsrWliqlREIRJvl8b5Amu7tWxK3VREZHP/uedlvprj1P\njZJzzwfmfX5eDWwyqcoun7/sFXUvfX1ab2m3e0dCuLWl3y+i1z5I712sXJoBcM79FyJS9zzvdw77\n+1tvvfXpz/fu3ZN79+5dzMJeEPE880jDYfX8ewUz9Hu+POyHvcfvAT4PeP/9fc2jU2S/jnxOvVxz\nu915P1y3COkqyP379+X+/ftn/nzP2UCdc9Mi8oee533e99rPisgviMiPeZ5XPeQzARvoBcjurnqn\noZDCDePx3hynWtXUjudZIfcw2d7W94ZCWvS9qoikboXCOumN27d1nwMxuah78EWRK88G6pz7SRH5\nGyLyo4cp/0AuTgYHVSGHQr2F+CUSmvbxvOOjjOFh9QIPW0+loorUOX3fdTAO/f1q0OJxhVEGyv+z\nclH3YCCHS08jAOfc74rIj4pITkTWRORXRVE/MRHZfvq2f+153n924HNBBPCCSLvd3YyDlRVLESUS\n16spyvNOp9zKZfWMMXbXMXUUyOXIaSOAYCBMIM8szaZ69n4l126f7PHv7KhnHw5ryue4966uWr74\nuhmA04rf2MViujeBBNKNXPkUUCDPj3ie5riZVpbLqTfvz/kPDOi/g9JoqPIX0bRPqXR0bUBE0yek\ngK4bUoSCdrdpK+fMAAQT4ALppQQGIJAzS61mQ+rpKejrU/QLCmx//3ADQM63W0UXj/cOQtloaHTR\niwIk6ZyTCuB+GR5WYxcKXT9j1wup1XQvXiQKkIuSwAAEcmYh7YMSJ4UTiVhH8VFpnXBYvfpyWd9/\nWR2+fgXdi9QSTXb83I0B6Dbt02yq8eolhPeyBWSYiHIipVKXu57nTQIDEMiZJRpVJV6pqPdMs9Dg\noHps7fbh3j+SSPSuwahbqVRMQVerpy/YniSRiEVJ5+nBNpsiGxvWwd2rJr7LFM8z5S+i1yowAOcr\ngQEI5JnkMCXuXHee7lWQWMw4byKR84ciZjL6vaSAzktqtc9yOD1vMxKcU6PZCwMaiEpgAALpWppN\nTZmEw88HTYOIRijhsFEen7c4d3wUdFoplbRoHotZ+q1XHE5XQbJZvedCocD774UEBiCQrsTzjDtI\n5OT0znWQWk3TCtHo9WAZLRZtLCTQWeY4PG/pHyQUuh7X5rpKYABeMKlW1aOKRk+nwNvtTq6W687V\n02oZ+ZyIKpqrnkIhFSKi63+eIrFALkcCA/ACSatl+Hw4d0TUC47Fjs/bh8Oa69/aUs95bKz740KJ\nzEzfVEoLxc+Sb6/X1RumL6DZ1N/DYStCHyfN5vmQz3FuRES9TMWkUpb7TySerUegVOruul9FKZfN\niQlgss8mgQF4gYTuXKRatQJova7FyuPyrH19qmQTCX1/pdKd1wzh18aGKuyRERvwclbZ2bGIxPMM\nyy9ig+OPk3hclV+9rkbjqPNut/VYKPiDBe/dXWtoazROZxhPK4mE7l27/WwF0Xpd183P0ejVj36Q\nZlPnB4jYPRtEQWeXwAC8QMKDDv1CMtk59ekoqmbE89TzqtVUeXbrfWF4MD7Npsj6uq4jnT6b8vGv\n9SCtcLcsIrnc4TQWftnbsz3a2VEFf5Dy4rjjNhoWdQ0OPjvsNRx+9nz/wet80nW/SnKd134VJTAA\nL5B4nuL2/SkSCqGx2NGeVKOh+fJy2aZ1FQrqRY+NneyBpdOqaAcG9LhEDq2Wfs9ZDEA0KrK4qMrw\nzh09t2JRv/80tY2T6BlQ6kf1B6TTujft9uERze6upZcKhe4jBFJmkYgiYU6T7mHNRxk1oLvVql73\n64SuYb2kgALv/9kkMAAvgDBqsd3WB8bvuQ8NqeI6zqssFi3dEonog8f79/b0gTwunx+NKh00k8D8\n3Z1nlUZD0yEi+l3ZbG/QIsmkyNKSKvHx8c7zbLV0Pw6bcHaYdFvzaDQM7tlu68/dGrVKxVIkmczR\nxrWbjmfP0+OfFHHghT8rbxGF7ZMkkzk5xRdIdxIYgOdctrY0916r6UNfKqkRoDlpa0tzqccN5EBx\nxeOqZGH5TCb1b6ct5g4OWtommeyeEvrgmk7ydM9DymVTltWqGoJIxJhMQyE1Pkfl5DOZzhRQN+Kc\nfqZW033u1rA1GkZrIaLG+ay5/VZLif5aLb3u2ezh79vfN5K+4wzOcdJu67GaTfXws9lgNsBFSWAA\nnmOp141Iq9VSZUaaolazDtJaTb3GYlHz4gcHlwwOmjeIsisUju4FaLXsYT7sQQbDvrWlxyyVjlei\nh8lFEabhke7sGO/O4KAVfvHQj/JII5HTk9hBTNdud8cg6nlq5Esl4xuKRJ6tVlAuW9RXqx090pKU\noOfpz2cxAJWKpcm4Hy+bIuRFkcAAXBFpNk1pnsYbBs552MNJcbOvzzpdBwbU26rVrKMUPhzy89QI\ndnf1fUNDn00ZHJZCQAns7FiqKJc73Ag0GlZcbbetIatbicV0DfV690Xfo6ReP3rweTpt+PtMphNC\nS+rDr2ihue4WWdNomJfNPoHvJzd/khGo1XReAtDWVkujhmcxjP5zcu7oe5Iu6oOfeZZjPa9NbVdR\nAgNwBaReV2/Y88xj7CYEBpMvog/7wYJYOKxKsly20Xv7++olAgNMp/X1WEwVFkZjfl7k4UN9by4n\n8uabJyu0zU0rFqN88JoPSjjcqUSPU3IYKz/2nVkEzaauN5s9GwbfvxeHpTr47mbT1koqrFTSdRMF\nVataL9jdVQMwPn48q6d/bkI0apGC/7p1U+j0N+lh8J+V1TSV0u9tNPS6H3V9hoetH+OsneGJhN4v\n9br+HHD+XJwEBuAKCCyUIp0Uv+Wyee4HPTDP64RwVquHK4p43PL65bLh90X0oY3FRG7etJC/r89y\n0FA/1Goi09MiExNHnwPrprZQq+m6j1IcoZAqj0qls/+gXtd/FFhFTEFHIlb4bTQsbYBXfhYD4C9G\n02R10PhiBEDlsM8Hj1etdvYElMvHs4v6mUiJBPB+/dftJInHNUqrVlV5nldjVze1h3D4fAqyfX0B\noucyJDAAV0D8iiQUUiXj532pVD6bRz7IlHiU8iMl4eelpwkKrz8UskEulYo++H199v5k8uQ0C/h0\nogRSEPv7RzeYsQb/Wre2rLkqmzUPVEQVPkPjobSgEH3WDlw/G2g0erSy7qbr9GBzWSJx+HB7/zkV\niwbBPSuKJhxW44wBOJjKYuJao2F/P+1+1Wr67zSGKZCrL4EBuAKSSFiagTQHA1VE9ME9zJPM5Ywp\n8bD0TLutChWURV+fvm9gQL22Vss+t71tirDREHnpJV3D/r4Vhvf2bMzjQWXlnL6vWlXPPh7XZi9/\nt+5JHh75/N1d/Z5kUhUVtQUigkJBf0cpzcycvWg4MKDfe1Y2UNbWbqvRu3XLlPrBlEi5bBDNQkEV\ndzSqn31Waoxo9OjUyeqqGdbRUV3XyEj3uXb6QDzPCvaHGRB6MaC+Pvj97bZNiwuH9X5PJgODcpkS\nGIArIgc9q0TCFHIioQ9Os6lKiveR7z1K/HnrRMJSJsnkZxWmnwuHn197TRVxOKwPNutpNs0739/X\nv0OrzHrw1hE/kVm7rU1crZYqJBRvPK7HKxRUWcZiIi+/bOkguIhENEXVbKohKpctjTYwYOms/v7j\n88lEEbHY2XsIaNgS0f0ZGTlaoTUaqoxrNb0GQ0N2HQ52tJJO8tcYziLAOYtFm37W369r6dYA+HmT\noN04zAAUCpYCq9c/G7UWCnpe7DspwHy++3nJgZyv9GzbnXO/KSL/toise573+aevDYvI74nILRF5\nIiJ/0fO8Qq/WcJ0lmVRvC2VMOqhaVaVZLJpBOKo4i1fYaGheGG/3MAWVStkxUMjO2XeXyyKzs/rA\nTkyo8iK6aLdtOlVfn34/sFFm//rXuLBgirxUEvnc5/RnkEPptP7zp6b4fKul54FxSSZVwaHAt7ft\nvbWa7tVhnnWpJPLBB7o3Q0MiU1OdtYRulO7+vkY5GGK/Ea1UrGZA6mh/X9cKpQYF53C485r4SftE\nbB+KRSuUdpsvd86uOQblqHvgKIFumhrFUdHWYU6EX2hILJfNUCQSuv/nYQCaTY1SRXTPAzTRydJL\nu/tbIvL3ReS3fa/9koj8S8/zfs05958//f2XeriGKyflsqFHMpmjw34/hYBfwZEi8XtaR82EJS1D\ncfa4HLOf6Owwr3ljQ2R5WY+7u6sPGO8rFi3FtL6uyhuIJlGB32P0p7eAYLI20lTttvURoDiSSfUW\n+/tVYe/vW7F8e9vW0t/f2bR2cI89Tz1x0hqNRmfBulbTn49DPdVqqmwiEd2PcNgUOhQXFMNp5oKK\nA3z/0JCu/2DtodXqrLm0WrrvGOhaTT/TTR4/FBK5cUPX6EeXnQZrHwppZHPSfdTfb+fdauk9E49b\nUbrRMBDC/r7tybN2ECPb253AgKOa1wIx6ZkB8Dzvj5xz0wde/vMi8qNPf/5HInJfXiAD0GpZpyYh\n+EHERqlk3ZUUaPmfzlm/sj+pXf9ggZRUAPN8/YrnMMXveap8d3ftu0IhXefwcGeKSUQfwHpdH/6F\nBX0PjUmZjBm2lRXdj7ExVRqsJZ3W1+p1NUorKyJzc/q30VEzhH191uQWCqkBWljQn+/e1d9DIU01\njY/r9+KJs1/M6221dL93d60ge1LB13/eUCp7nsJASZf195sxE1FjvLmp0Q/pqcOUeCxmXD0Y0IOz\ncVdW9PWhoZM93aEhPd78vN5b1art/VEOCNcd49gN1XU0aik5FDvpIviH+vt1/YmE3g+JhO1Po2EG\nZHDw9LWBgwSBgZwsF515G/U8b+3pz2si0iWLyvMjfs+OnynAkq6B7rdc1gdveNi8Ubp6oSVIpY7O\nc+P9tlrq8a6tGXIGJT45+dnPbG3Zgx8O67HGxlSpNxqmiNNpfYBzOY0OIDsjf12vq6KiSJpI6DHb\nbf0M51ytqmIsFPTv0FLU6xpRoPyggHbOiqnNptFa9/db1yqTsra39TO3btngenD2+/t6TiKGQKpU\nFBY7NXX8dUwk9PvwpKNRXV+hoK85p0bFD29MJHS/h4eNTiKbtUiEBjyUYCZj6bhUStdG4TuZ1HPc\n2/ts5zbSbBpRXThs/RSlkhlsai57exY1QRnCdVxZ0es1Onp8TWV31/iLoA3nnhLR+2VlRV+PRPRe\nxAiI6LUglVip6DU7jQwO6rlcp5nUly2XVnrxPM9zzh0JLnzrrbc+/fnevXty7969C1hVb4ViKYiW\ngQG90VFw1aqFseGwKqJsVhXT+ro+MEND+tCMjBjF8sE0R7NpGH7CbBTh1pY+6CMjplD8ntb+vqVo\n/IpicFDkh3/YEEogWmgkm56279ja0s8Wi/pAAxH1Q1GBu4pYymR/X41UpaKKkoJqJGKIEYq/GAC8\nZPLk/qipVLJc+sTEZ40vUZCfSRQ+m40Ne2142Brk2Gd6A9h/CsEIym55WRXj3buW6vDnu3d31fig\n+P0Q4I0NkVdesVRYLmfpMv+5tNu2NtYTCqlS3983VBW5dphLnbN0IqiqcNjoHUQM+UX0elxjG4Jx\nck7vH5yXREKNCN55LGYpsVBI75u1p+5hvX56A5BMXp+5Bucl9+/fl/v375/58xdtANacc2Oe5606\n58ZFZP2oN/oNwPMkAwOdBUa6OAmXw2FTCM6Zkt7fN0w/A7L3901JQeRGIYwiKQVUWuxh7ozHrenL\nL6FQJ0c+KRew69Af8z68XCINoKy1mhqwdNry6yg70ibw6CwuanqC6IAIKJEwxU1ksLtrmHnSZIOD\n1nMQjYrcvq3K5IMPLEpZWzPitmZT5NEjq6MMDOj34Inv7ho8sdnU1FI0amksv5IpFvW7SM8NDOjn\n/BQVpZJGOCMjuk6MUiql3w09NmkwoJSep+eI0iXdl8no66SxUJp9fZ3Ng84Z7JIGvdFRu44iZjgO\npsFu3rTPxeOW3jlO/Od2547xGUGgJ2Jr5x5cXdXXMUjh8PnWBZ53Oegcf+UrXznV5y/aAPwfIvIz\nIvLfPP3/Dy74+FdKyPVDiNbfrw8gcMhqVR8cvzKBlwWvVcTw1dWqfg8PWH+/Td4iKsjlVDkwVpC8\nLsqbprF225RAPN4JP200OlE/N27o8Z88UY+3VNLPjo9r4ZGpWiisUEh/LpcNfQRf0fq6fo60ytCQ\nRR3VqhmAaFQ9c+CtuZyug2gjl1MI6caGyCefWKHW34RFequvT1NXoG8KBd0X9mNpSZXb+Hhns5uI\nKl+M1vy8yOuvmzFnRCUpD65VMqnnvben14i9p7BO/0Qiof9TTF9ftw7rkRG9xsvLBgRYXTWDXS7r\nHkBnTVE7k+mM+NhfKL9XV82THhrSyAWqh5OYWw8S30FNTV8Axtif8uG+LhY7PxvQPV+M9BIG+rui\nBd+cc25BRH5FRP62iPzvzrn/RJ7CQHt1/IuWVssKXt16L4TawAD7+9ULxROFmK3R0L/Dgz4wYN4S\nnj4pFnK8OzuqIF95Rf+Wz38Wl81owGbTUjX1uuWUt7f1AV5f14d/ZsaULN5dNKrr3NpShba4aIao\nXNYwPpczagke+KEhYyhFOaIYUV4iZixIKWQyejwI1xYWTJnCm8PIyUZDld/8vJ1zKqUKCDQTBowp\nZzTiJZOWFonHLdpiXdBGEMGBtV9d1b0cHLSu4Olp4/Z5+ND2bWLC6hGhkHUOv/SSyOPHuoeFQmdq\niMI3xn172+odfmhoKGRkfIuLFlWWy50GgPckk3rMgQF97/Kynlerpddqd9cMUS7XHcTSn/KDpdQv\nGGR+Hh62xr9n6X0IpHvpJQroLx/xpx/v1TEvSxoNw3dHIvqAdGME8O7JcZMXz2TMIwPqiPcGyRep\nFMLrdNqKgk+eWKE4mdSGroNSLJrHPjdnDVmZjH5XqaTKpVSyFAgoIFBEImZsikVdE4aJYjaePtw4\nDDtZWlLFgvddKqkCov5AVy3fUa0a/XShYNz3UCt89JEq/qEhfa2vT9e9vKx7wj7fumU5fnoKiDoG\nBzuPA78R7KY0wO3u6ndSTxHRvUul9L2rq4ZaAu1CNzZD5Ot1XQdKmzQd9RWQP9vbnWyt1Brm543E\nT0T/Nj5uhV7ScvW63jM4KEfh7YkGh4etqN5q6b6K6OsQDh7FOyXSyZbK/S1yuMEgHeR5NpToqIJ2\no2GOUIDvPz8J+u/OQQ4jc+sGwgaTYn+/pThmZlRRkAv1pxJANzSblmslJYQnODurytXzOvH0o6Nq\nTEIhNVClknqj8/OG1Jme1igBdtKdHZF337UUDTWKdlvXSpheKllDWjRqKQ2ihe1tVYqFgirevj7z\n6IFCTk/rund27CEH+kpzl3P6eVJgzCQgz0xxGGw914Aiajis+0PqDMROJKJGaGRE94HcOZEE6Zl0\nWtdIExtKdXhYzx8DfOuWYfeTSUOmtNvWPwANBD9j7Gs1Qzb5p7BRvE0kdI+Wl41Se3RUi+axmJ1z\ns6nfh8dPzeEwz5pGPqC1+byhora29HNbWwbjPA4J5KcUGRgwEMJhndYU3U8SP1su0WxQIzgfCQzA\nOchhZG4nCUMwUil7+FE2t24Z22WlokqCVEOppJ+nSLi8bNDN3V1VtHt7NvkrHBZ5+21DsoBOGR7W\nhxU6h81No4wmBbK1ZQyi8bjVBvCmUylViCL6+0svqSKClAyqZtJL1BdgN202re+BsZHkzCcmdD3z\n87oGDCudus2mNZjt7alnDrcR1BXVqrGHbm/rMVZW9PsY9k7RdXnZPGx6Gcpl80pbLTWYNK5BmJdO\n699Ib9y8aR3CtZqeXzyux6DxDwTYkyf6eypl50Y3NMY+mzVFGovpZ4pFNQLb24ajn5mxgn67ree5\nsaHfmU7rWo+CtnJOAA1IRdVq+jkUPw1+/vsdLiTOA+UPHPfGjbMT9SGk20TM6Ab8QecjgQE4B4FH\nnofnpBC1UtEHWKQTailiKQ8/9QETu/Aw/bQNNN3s7OjfqUP4+X7ggGEYzMCArndkxNAdQ0Mir75q\nyKOlJf3c6GjnnALQPKCRQJTs7mrK5fbtztTW/r55oiJqIMjjZzI2pN45kS9/WYuqrZYaAKCvGCEK\noNAh5PNWSMxkdC9nZgxd5Jx+9927igiCGykWM1QMvQMov6Ul3R/m8NZq1jzXbut1oMErl7O1O2eU\nG9QcgHNSOIVKemBAryPXEhQP8GBQVP39NrN5eFiNO8VsoJb5vBX1kWpV37u6qucxNXV8Y1ssZoRv\nIgZPHh3V8wVl5TcggBRWVgyCy0Ajahf9/QY88EcNdANHo50EfBgTrgtePjUYIoBgXsD5SWAAzklO\nQ5N7kBgtEjEMNx4gEg7bSEa8apqj/N4qeWLSAJmMTfjiQVpbU+VCjv/zn1cltLlpBmNry7zGZlMH\nwWxsqFK4e1e/x9+YlcmocnVO0yjvvGNFzY0NSwklEvr3qSmDuO7vmwIdGFDlC+qFIijpLbD25KZv\n3NC1Fwpa5CRFtLlpHEqZjL6f5q5QSPduaMgiCKaeVSr6vaRJaJSjuYkUCAYarzwety7keFwjiaUl\nXQ9pMH8fAdTa9GKIWGfs0JBRcjN2MxIxcMB3v2s8OjdvWr7c3+UsYvBRIoVMxrz0w4Rz6uuzmgoY\n/VTKIhGQZ9QIikXde+pP8bgVybnmpDG5p9lLv8HCgAEFJkXId8ViRmkSjwfpn/OUrgyAc65PRCZF\nxBORRc/zSj1d1XMoe3vWXUuTEkqdMPvgHN1yWT0sUkUUQ0kLUSyMRlU544WOjup7SI9UqxY1kBZp\ntRSrTZ6ZQjZFx2RSlfXwsBYXUylV4HD9kKtG0fT1KR6+UFDDsrOjyoH1AQm9ccOUjIi+b2pKj/vk\niRHMTU6aIqPIzZ6xh8BTyZknEtYbQJ4eY3Tzpr7vlVcMvkohd21NlQ+GuL9fFc7amhmTUEj3NZPR\nvxEFFApZGiWdAAAgAElEQVS6b9Wq7g1eKjWdWk2NyOSkwTdrNSOi6+vTa1evq8dOox+4/EjEIiEK\n66QMczn9eXjYHAiav4hyiPRwIg4zApxfLqfrQCFTTE6l9PoUi/r9dH1DLkjNByNMyuww8riDPEd+\nskOekXr9s17+cXTXxwmINH/z3MDA6b6L2sxpSPiuixxpAJxzAyLyCyLyl0QkJ0rd4ERk1Dm3JSL/\nWET+Z8/z9o/6jkBUSIOI6M09MqI34XvvWXESegARK3C+/76lBJzTQileGvN8UynLzff363sHBgxv\nLmLKYHvbeG4gbqMfoNFQhQ00tF433vdiURUl2Pv1dX3twQPzDj/+WL83ElGFSP6Zf5lMJ0dMpaLo\nI+T991UBUrCcmDBlieIaGzMY58CAKmhSW3jjeKPUCcD9gz5qt80oQQONkS0WOxvqaOjCmx0ZMWQV\nkQPzGFCSrH91Vd/fbBr/k4gawtlZQxAlk/o9y8uW8gmH9dz89YFIRGsQIkaTQZPghx8aMiceV6M6\nO2v7AUIJNBRGgCIsyCTP0+ubzZozgQEqFrWfgk7slRW9RqTx/DxNFGrx/P0eO0OA6nXbNxGLvqh3\nnYeipf+E54loCidJRPetWDyanLFaNYZRAAPPU/3huAjgD0TkfxORn/Lx94iIiHNuTJTY7Z+JyI/1\nbnlXU8gJx2LHM3oifmIqsOMffWQe+dKSyPd/vz4AoCagiCBKIL0TiejNSyQwOGjDWkQMBz88rA8s\nKSY/lLFa1eOsP+3DrtVUYRCuT06qorpxwwrBIpZOGBkR+drX9Lv9iKOZGUvvoHyZY7C6anWPyUkr\nQHuepnvoEo7HtY6QTndSHLTbWhyHJA7vHKbK8XH9f23NsPaplJ5DsagGanlZ13PzphrToSFVZMWi\nvk5B/u23bag7fRfMNyBig+6CqC2fVyMBf1E2a4VvjFihoMqahi2GzO/v24CYmRnjBNrdVUNKLh7j\nRxRVr+vrW1tWX+C6ky7Dg4/F9H30dUxM6B6DhCLl5ycbjMcNTszcB5rCcjnd75s3rWDvl1DoaEVJ\nOgekF/ctvRip1PE1i26fv3bbDC9oMf/zWKnoXlFwD4U+e9yDpHLPG8nckQbA87wjFbvneasi8j89\n/fdCSb3eOaoR9sjjBOy0n7wNFspSSW98oJJ9faY8xsb0tXBYFZyI/j43Z3N6iRCAfNJtSz4WBUOO\nWqST06VSUSUDV7yIKsahoU62SnDsy8v6D6+JdNGdO8abMzmpipgHGnTLyooqHJRoMqnnQqGRYmSj\noQYNiCA8MSjgrS1N2ZDH53zpaRDRPXnyxIa0LC+r0tjZ0e/e2dF14t3fuGEdsTs7VmDO5/XziYQa\nIBrCQGllMtbBC4/OSy/pOVM7QaHCksrnmWuAsdze1msei+lei1hXN8NUMK4gxeidoNO7UrH6Cc14\nMzNqcGgkIyohUiTCOpgWgUwO0jzPU4XP+WSzFpHWanbPHPfs4ARkMkZ4CBQ2m7XGv6MEWDT3bjR6\n9DAfjGKp1FmsB7EF2yk/UwsiahLRe9SPhDvr5LmrKifWAJxzPywi73qet++c+49F5E0R+bue582d\n8NFAngoplXjcPIzXXtOCaSikyhOsP5h28o1TU+blzM4qLn99XQ1Bq6WG4Ytf1Ju9UDDoJ4Pd5+b0\nuysVgxBGo4r4IfXEgwBfDpBDumJp9trY0PNYWtJzYpTlyIj+DcgiKYeJCVUen3yi6yiVjKqZNEx/\nv651ZcU6ZJ3TdNTmprFfNhqqfCB/w3NdWLCaRLmsv6P8Egn1/LNZUwCkAkT0O0dG9PfFxc5Uzdyc\n7kcqpcbmzTf1O3Z3DbZKDcQ5TYeRp/c8VeQYUuoD8/P6fqKGl17S/995xxr/aHKjIO6fFUBufGvL\nlBm4fhBow8PWDCeiCuyTTyydV6vpe/b2NFJijUSXfmVKTp5Ir69PDeKXv2wIqYcP9d7Ag87lzMCC\nrAIJRbf70JAZ/Y8/NnI6ZiQAawbtRfOjSOegGX89jP6KWEwNO5G0n0LFT5a3saHXxN9oCL05XFB8\nBkfqsLkSfPd1lW6KwP+jiHyfc+4LIvLXROQ3RIe8/GgvF3ZVhbzxYTC2w6RSsUEqeBADA3rz/9iP\nmXIFTw8HDgU58qXf+54+UO++qzfr/LwdG8+7XFZP9t13tRC7taUPKPlciL5g3ZyaskhjeVnXVKup\nMiwU1DDBHUNRcnXVWEmnpizFsrOj63j40PiI+H9iQr93ddWU0NiY9RaMjlqNYntbFRYREjBYiM5Q\nNIuLxrEDMmpnR78PGmqQKyjlzU1LPQC1BAJaLhsXE+gqoKCep2tKJHSf3n9fr8HgoDWLMQ6RdBkT\n04aHDVIL1xBEaB99JPKlLxn3EXQSqZSlaqJRNRSQu5E7h1xuZET3hq7uDz/UdaE0SSHWaho1zM3Z\nezMZa/C7ccPmE1O8rVbVoH7wgRX/Wy29dxj+QgSBEWg29XzX162eRAGXonapZLMhoJcoFm2AC+kt\nkc4RkvPzltKkxtDfr9+3vGzw1Y0Nvc82NnStnCNgAmYppNNG74FBRw4q9YO/E7XTxXzWkaKXLd0Y\ngOZT6uZ/V0T+B8/zfuMpl88LK4ws7EbwmCnajo1Z6I1yA9pIQ9bkpN6MDDWhq5QwnYlcsZh54qRq\n8JArFUXlwOfTbutrTIcqFvVvpE9KJVVIFGJR7AwTWVxUZbC3p+ufmNBcfbmsKQcM2N27xlAJiR0T\n0Pjuvj5r8Nrb0z2haYoGObpv6QTGe0Rh0zhH+qNcNtghkVEspoqBxqhsVlFAdBaPjppxCYfN4+3v\n18//8R9bjWByUs/x0SOLAqBZSKX0f+Yt7O7qv1DIril7QNG4WlWl+fChNVqBvmEPRkeNduLGDb0W\n0HPQM0AxmFnD0HD09alio2aEgV1eVmMGVp95DhTBIe0bGND3zc7qfeHvLJ6d7RxXSo8E9zKpKqJZ\nItt4XJUyqbOtLatJ5fOdzxTF6nZbo4S5OZsfQN8H0c/B5kt4k/DaS6VOdB0pJojpINbjOemGh4jv\nZ4jN82wAis65XxaR/0hEfsQ5FxaRoBWjS8EzhGQNkjJC3Lk5vZHW1gxFgwdKXrlUMkKumRlVstWq\nenLhsCoHlBj/f/3rqrB3dqzmMDCgv+dyWj8oFs2YMMiDBxhe/aUl+/zenqZh8Kh2d03BVqv6gIMa\nQokziGRx0c55f1+jBwrDQDB5uEGvkH7CeMVi+n68Oj+CCsy9c9b1e/Om/o5xqNVUAdGdCp0CyJ9W\nS/c3GtV13r6ta6OZ6eFDVYYUbDHChYLONZ6YEPn2t61bGmbNSMS4hzDmKMtYTPciEjHkEimRxUVr\nMqSPoq9PvfxmU/sC6FmoVEyxbWxYZ/Zrr6mCA7K6stJJrcGAnXffVUVLNEmakdoJzYLQkEO5Qf2J\nSKhet/s4HjevHEAA12NpyfoT4nGNTkIh86wxnByPe6DZtGPxfJBapb4EWouibqVi0888T68la/KP\nAj3NCMlw2Boigdtex1kE3RiAf19E/gMR+fmnPP5TIvJ3erus50NADABn7O+39BGhM6khPw8/OeZI\nRJUYSrndVuUPGubRI32AqR+MjalyXFzszKlDvdDfbykTbnqKXvDv0KjUbuvDU6/rA4mXB8fO1JT+\nzEPNw5bNWlMaPQIi+nBiTLa27AECRVSr6bnduaPfeeuWPVwzM3o+eNI0SwFBpVt4b88YUKHMuHXL\nlC8G8uWXLS//5Imu48YNvU5TU4a+evhQ94yZw7u7eh4oRAzb/r4qT7qal5asMbDZtMhlZEQ/T/Ha\nf124LyoVy69joBikAn0H942fqbXZtFQeYIOFBd07UnArK7rGbNZYRVstkW98w9JKiYSdp78RbnTU\ncP+k6rJZvT9zOesgB4kE0R3Ke2ND95p+FAAQeOdEdcwkwOADK6bYjyGEPdXvrfuH1RARJxK6xoN0\nFN3SU+AUHWzQHBxUA0n0Vio9pwbA87wV59zviMgPOOd+SkS+4Xneb/d+addTyAHzcHPDUnwjTSFi\nnkMioTcUNBKMSiyXDarJZCzyu3SKfve7hmiheLi9bQO5o1HlqF9a6hwEs7Ki38Ws1zfe0OMODtoo\nwlTKoIgUa3kQFhf1ZwrOdIym09boxrAahrfAExOLWQMa0QPe/fKypgO+/nVDydBNzIhHPN7tbV3X\nxx9bPrtQ0OimVFJDA00BvRTJpEEbQeCIWAMVOP502uglEglT9PADZTIGj6ToOTysUUM+b5z69BHA\nyElqg45jehhCIT13vHiw8MvLqqCZkQvhHqkmajOZjKatWi1N0bBv7bZGDJmMDaf/wheM14j5B/m8\nroVIbnzcIq1k0ig8GINZr9ssA/Y9HDYqDSCe4bDePzDI8nwQkVHHmZ+3ehf05Dw/oLCgQ4ckkFQq\ntR9qPsCcT5rrfJJQN4LTamjIUED+4UIi15eeohsU0H8qyuX//z196b93zv1Xnuf9w56u7BoK3qGI\nweLgWsf7FDEvu9lUxUBhEj4VoIKrq+YxMzAGz4immZ0dU5xQ+aLoQOkkEupZ1+vWTQskFNz++LjI\nD/6gcceTf8ZwNRqqJEiXUFugX4H0RrlssFKQSKRLUDR0HU9NdXYqi9iwmYcPjQJARFFL5FuXly0X\nTrd0oaDfR0i/s6O1CfoWhoeND2hryyaSwS+0tKRKn1GX2axBAzc2bEZuuWy02Nvbhq33X9uxMfOk\nUXR4z6CfgCO+9pp66F//uqYDIxG7juwhNAmk3jCa7BfGDwgsnP7j47rutTVdY6ul+z8+roZqdtZS\nIUQkwGKZRT0wYDWPbNaMEDTOzhmktFq1ngSKxDghIyN6T+F4gHjCeED10WwaSIKmQQzQxoYNJsJQ\nrqzodcWhAF6aSulnnqWhjDSmiEWpfhhoLme1qOe5BvA3ReT7Pc/bEhFxzmVF5F+LyKUYAEJHWvYv\nAoJFoQ2StaP+Tl6U1/ykZH60ELDGSMQofcmREqJTICsU9EGF5mB2VkPpkRFLiczNWaqGPDi0C7Oz\n+nfmAmxs6HFu3FClAgvle++pAiX3/PCh3vRQReNVodAp4gLhZMj66qoem+JqPm95bf/AeNIXw8OW\nw282tfDIe1ot65D1M19SR1hYsAYf/4NYr+vegKqJxxVFtbho9NsjI0bDvLtrBfatLStIQ6lAVMD1\nYmiJiCG78GSZf0ADHQNzqDuwl0QJH39sQ2SgfCbVh4HHGcjl1FjQg5HPW7/F5qaey8iI3j+plF7P\njz4ySOburkaD7bb1PTCYvlw2egm4+b/0JV3zo0f6vmpVPwcdBv0d1Kk2Ngxp5hea8kIhNTzMzqCu\nRXFYxAbOpFKdefxUSg0Ihfp83gbx8Kzw+eFhIw/EAADRTSa7T/9wbhiXg59j2p6IFcVJVV0X6Wap\nmyLip3vYf/rapQgdpCJ6YXo9Og6PESIsZtVSfPOTW+EJ433D0c/3QEJGHpdiJd40nat37uiN1WgY\nA6SIkZ+trZkyINcMiqHZ1IcMeGCzqQYjnTaOFiZ47ezod6ysGHpoeVkV8+ysoUpIvywsGFcMDx5R\nBd2Z4bD+jaLY9LQak+1tGyZTKun75uZU6YyPGwYfcjuKfHC4kIpCYZNbhj8nFrN01Fe/ql4h6QQg\njf6ay/6+nisFb4jjqHXw/Xi8yaT1BRApULcZHrYC+MyMTSkjgspkLO8ej+t1/eQTWzP3C6myeFwV\n3J07+l2kQejDYI0otNVVAw3QKEaRe2XFILg07hH1ManN338A9HV01PZVxMALUJJQayIlB+x0ZMQY\nVylq07E8Nmb1EBGLcChGp1K6JuoIm5sW7WYyeix/iowIk9QZUYxIZ9MWDoeI7h3G/yTBkaIWdVRn\nM2lAjjU6en16A47jAvrrT3/8RET+xDn3B09//wsi8l6vF3aU+LG6/p97eTy8PUJacsfkAPl7Oq0P\n4cCA3nwrK53fQ3pnZ8foBVA2wAY3NqxoBbyQCKNYVOUMUmh4WB+qzU19nRwvofnamqY2QOxMTOhn\ngTKS5tna0vUQXczN6UNFwRjmzXzeSLuAwOF5o2RIGbRauj944RSK+QzREcqYVAepM4qU6bSeB2Rk\n7CVpJZRDJKLnSu8B8ElytyCBiBLwWkE7kaIhYqGjGtK5aNTSQqBPbt3S46ys6PUaGNB9q9Usfz4y\nosfLZlVhEsEuLdl6olE9rp+076tf1c8xbxlEEXtFBCxiRpriO/DW1VVrehLR46yuGocPcxsYjck0\nuOFhPa7nqXFmWAx5dfL7+/tq3FMprSFBnwFnVKmke7O6agR3NFslEtbXwXQ0oiVQPRDDDQzYvS+i\n37O9bdDqqSn9TuZRMKsC4+Vn3wUm3O1UMWpIJ+kIBHj3dYkCjlvmgCj756yIPHr6s4jy/3hHfajX\nAvmZcxfDzMfDz02EV+EPC0kdtFqqvHkPKAg8YwqneKYwI0ITUa2qsqLhBfbNkRHD7qdS5v3QVAOx\nFaibUEhz5hTLqAusrNgDv71tg1ReesnQFih20lagWCioQVlMcw/jIoET4mHSxwAaCO8SXDq0F8xQ\nBubHQ09qgP9h+0wmzZCJGFwUzDh5d4wSOPJUStfO/AIYQTHMID3oDiZFgydLvprRlhQ4ofkgdYNz\nQCHZj0WPx7UugVcOdBNEVDpts4ohDLx5U5Us6SHYXEGfPH5sRgSDwcCZjQ29zhgMuKaqVTPw/f16\nr0BJks+rQi0WrWHr5k39noEB3UMiiQ8+0HPBuL/+uk2fg07h8WNda7VqsFacFKg1CgW7B+BAIlJy\nzlBVoZA6J8Bvi0X9mbpbNmtken6hUQ3QxcGCLWk5YKunlVTKjC9gjesix3EBvXWB6+ha8PpO4gw5\nLyFUpjgFqyfeFwVXCr1MvIpGrb2drlIiAopbfigc834ZBUnNgeIpmOXJSf0sRFiffGKUw6R0UIoU\naOkBwBiRlqL2gGGBfgL8PDfzq6/q+REK03mbSmk6gD2hZiFioyD7+9UA9PcbTQR8OhT0VlfNKDL+\nkMlPjCQUMXplkEbkljc21FhiZOj07euzVAkeM529eJp9ffbggpphf7iG1CugOiDKY3QnXa14/1yr\nBw/0Ok1M6HeMjxv1wOSkTUijn8E/RJ2JYTRRjY1Zw1c0qh43hVO8/+lpNQB0uqLURSyVyEAfP80H\nTXzUj+hcZ85Cf79BL/1rhnIC2CiIJFJu/mMQ0dJ/QFcuETVrxCHY2DCng5QW0RjXCCMDjBZabL9w\nT1B4pqMagUJCxIYYnRbOidEDeXWd5LgU0B8e8znP87w/f9aDOuf+lmhjWVtE3heRn/M8r9bt5y96\nkykQQmIm8tmhFP41+fnO6VIk9KQ4tb9vBdxEQuSHfsiaoWjdB4e9uqoPGqmLfF7X9PixjTKs19Vz\ny+fNy8eT5CGlPkGkEI8baue11+zmf/NNS9OQh4dr6NEjWwNUFRTseLD8ueRHjwwyx8NfLBr0lTw7\nXbRQReAV0mzElCsK/xSGSW9QkJyY0ALq5KQia9bWzMt/8sRoi9Np/V6ULp/P53UdkYhGNyBoSFsA\n3ZydtWlYGBCMKwaD41DToNCLx+zPkff1qQFnVgN1JTz3uTmt7czM2D2UyRgTKeMt8bKJOkAs3bxp\n9ybnByQ2ErEIlugiEjFjJqK9FTgZOB8TEzZhDm8fh4BeAAbN12rW8Eb0vrlpTgrPFfBXaLAxVjyH\n/og8mbTaw1FDYujFAY58Erund4bcBt3qRMvXCRF0XLDy68f87cwpIOfctOicgdc8z6s5535PdObA\nPzrrd16kHHajEQLCIb+wYKmYXM5CWlgjIVZbWjKPeX5eH9J8Xh9i0hpg2NttfeByOfXGuHHh52Ek\nYzqtXalvvKHvKRbNW2ZAON4xCjeft+EiIyPqCb79tuG3/YqZWgWcK/7Zw+GwKlh4dRDqGMPDqiCI\ndEiNra6qAgIZxfdA3AUJG/MAolF9Dd4fUkXZrA3DwdPGu4ZuA0QMw1kgI8Mo0rNBgXhw0JBBg4O6\nr598Yk1a/n4BKBdIY4TDVmzH44U1lbQT12h/3ww71Nh0CqOUFhaM859akD+FSC2A2Q4Yz3hc7x0/\nORrT2Px1KQxKLmc9FMBoV1Y6G/yiURtyE4tZnebxY7vXGG4zMGDeOpEvjYqQuRFFFwr6LExO2uhS\nIgaG0/vJ406aDkYRn3rHwbQxkS0poLM0cxWLFjVyLheRnTgPOS4FdL9Hx9wTkYaIpJxzLRFJichS\nj451ZqEb9qgLSc6bghK457U1Q2tsbOiNfOeO3nwgbtpt9eaY4ytinZwUKXd29CEn9UC+GC71jQ1V\nQsyT9fPhLCxo+gGa5mJRlQt4ecJ152wAO5w2Ipbb/egjm9FLzpxiHUqAfHA4bA+55+k5UOAOh62e\n8O671rUL4RuRBMqBln7EH0XBgEl3K0qJ1BhDb+gUplsX5AzoKlBPg4O2h7WayLe+ZakSOJVQ0JzP\n4KAaSUjhQMGg+Fnj2po1vQ0MGJIondZrvbFhyCAiNthKKZgzfCaV6uxNyGatd4EeBdBiFMyHhy3t\nxRSs9XW95js7VovwD5PhPmDPQUaRhqHuRIMbe7mxYcaLvfDTLOAQrK930qBDj8LaNzbUaCwu6j2b\nz+vrUEhQV6tWrRGMe+Tg84kTxejMo4SZGv7U0GkKxX65Lugf5LgU0D8Xkf9FRP6553nlA39LichP\nicjPeJ73b53mgJ7nbTvnfl1E5kWkIiL/r+d5Xz3twnspcN9z0+OFIJubehOinPFc/ZO36DzFkxPR\nh9s56/rd3dUHIZ3W1z/5RPPEhP943zA17u7qe0CFfO97hlNnItXCgtE1M9gFPD8pCMYWep5h6d9/\n3xR3o6Eph/V1U+CJhHpmYLTxpFBgGDK4c2hcm5627l4gk6CTSAdFIgah7e+3WoK/6E0DVSajnawT\nE3qODx+qsqBbmTm71BlItRGh4e05pymh1VUjQcPgk2smRQW2mzRSLGbkZNwHFEjTaesqhtohEjGM\nPc1EoKDgBFpf132KxfR6j41Z1ytdzYmEGdKPPzbaCAr3kAdCC1Io2ACcbNaMBXORmZNA8X521mia\n33hDwQFPnuj6GC6Uz+v30YlNXYmeCcj/mD8BpTkd0DgTiYTtDTMjnLPojO7wdtt6RbgfiAKJSEUs\npUcdgnx+Lmepn8Ny9H5qj2zWwBA803z+OPHPdSBNe13kuBTQz4nIXxGRrzz11FdExInI2NPP/Z6I\n/MxpD+icmxGRvyoi0yKyKyL/xDn3H3qe949P+129EjwewnRGLYpYGzyt+ITz0Dgkk6p0aZXnhqWP\ngPBYxHjqwZDv7Nig8kbD0hWPH9sM229/294josfY37fuTxrQ5uf19ffeM4oAsOKplOWlm03N05fL\nuh6OCyyTOgCFN8Y6Hlx/sahKnYErvDY7q6mRwUHzlkWsmCdiCoyUTbttHOzFop0bfDQbG6r0FxcN\ncQO/DGkYogquAcNmSHOsr6uSoreAay5iRoniORFMPq//49mnUjbgBHpwRixC8VAqmVKikxVKbZQ8\nEE8R6072PC2+P36seX9SfX/mz9iwmHxez8PfVAjKhtQeNZxvf9u4dOCPgsaaZjLAA5ubev/gYROp\n+SG33I/sUTRqRXF4+Wka5F7LZIzjPx5XQwL0ky5kprxh+KB0rlatKxyvnj4bz9N9wMDQeV2rWWpp\nc9Put1zOIjOeRaI7xqOKmMNwEtrQ3wdx3eS4FNC6KAXErzwdAXnr6Z/mnk4EO6t8WUS+5uss/n0R\n+bOiM4Y/lbfeeuvTn+/duyf37t17hkOeTkhrgA4AISPSmVfGi261bBoSBgCaAOgb1tb0YRYxPp98\n3rhVCIVpTNrc1FQMqQtQLtAhkHZh2At4cxG9keEEghufh2d83JpaaNoi0tjdNcPHgzMyYsgaCMVI\n9UxNiXznOzagnm5IIJE8nPQ9AMEDTgkSJxq1ovTAgFEbMDKT1FalYqmLVksNAPl+YIKgrSDIgzoY\nw5vL6V7TH+Hv9ERZ0fBDamR01OY0MDwHyKYf4oqCAbtORy7DWdgL0kn0crz5pnr0kLHdvm04e2Yn\nc72Wl/VY29vWnc0+0oDFz85ZmowmKpoNs1nrEscxoSDsL67CtYNRTqeNS7/dNh4jah5EohRq4SGq\n1TQdCQrHj1YTMeNw65YN9qGRjt6V1VUjHNzasv4Z0oZwAdVqVqz2d1cDbWXeNWvHyGHYcH64J66y\n3L9/X+7fv3/mz3eFWH2q8J9F6fvlQxH5L51zSRGpisiPi8jXD77JbwAuWuDXp2DKGEjQG5mMkbpB\nxAVZFR2+jx7p5wlHCwVLDU1OWj6Wh5zPUmz78EP9DEoUjDSKYHDQED3JZGfTGVw98MaQviCNNTho\nhT68cmgBRCxPjqED+gf6hPkBUFAANQX7Dr4c6l2as27ftsY1ulBHRnRv2d9wWBFJlcpnicX8fPNE\nVtvblqrJ5fTfyop5rem07vfjx4YqGRiwmsHSko1ghOguFrNmN/b6zh09j8ePrQO1v1+ZRZn5QMoP\n0j/uHxQxUE+aquACooGsXLau5GrV+ihIMQ0OGtoEJA4NicWifh+v470vL+s1BY4LPUU+r+fC1DWg\ntUS4oII2N/V1Ij2cF/8sBK7t5z9vNQZmSKyv63fhULz2mg1Q4RrxfgTiOgYD7e+bUYFKG4eMlCZN\nX4xAhf8fw1Yum1PAPY2RpqMdgkUGI12HEZAHneOvfOUrp/r8hbcseJ73Hefcb4vIO6Iw0G/JFZwt\nzE0I9UCtZgNS4LgBi1yvG6kZTStgoVGMRAdgo/N5Y0DkYahU9IGhO5eaAR2Z5LNHRkzRtNsaXfjp\nEWBRhOOfB42u2+Fhw5NPTVlemzAYWGOzqVEI34ERKZctT7q3Zw1cYOSZaAbsk4Lh/LylLIAach7L\ny6pUm00tPmMQgTLSGATklfqJiKWQmObFIHrQL+xho2FFSJqGwK2Pj1t0Atsp1BSPHun+PHlisFD4\naG38dNoAACAASURBVAYH9bMYJu4Zej/oXcB5mJiw44+N6f7Pzup5gPh6/33DlYPAKhR0z/GsQdhA\nnQwPDTDijz/W72A05/q63r8TE0bHQA8Bho69rtVEvvY1VeiFgn0n908kYp3CODGk3ogQHj2y7mcG\nuUCA126bcyNi8NtqVa8HCC16PB48MMZPoh4gtLCq4rRVKrqnflJGKDJwCMJhqwXwXBK9+dOkL4Jc\nSs+a53m/JiK/dhnH7lbwPHko/PQHKFQKs+QOUdCkQcCUk/fFeNC9iIeC14RXzgQsCnqplCpFcOt4\nhjQq7ex0ttL7IwWQNdAowPw4OGjNQuGw8QiR6wbbDFSPQrafVrlaNX4k0icMvwHHD1Ha8rLh6RnI\nQWevf9pUq6XhP2kpIpRy2egzMhmLUMjni1jahCI70M/1dfN+QVRB1ra/bx3JcBexpo0NXfPqqn6e\nrlruhd1dI5/zs2SinCoVM4LsH93D0DcTVRUKejwiTlBVpHWaTQMUhMP6+9yczdBlfCQpHFJMdDX7\nKY393ej0ENDnQZc3RX9SQFCNM6nt4UONfh4/ti53upHfe8+uCU2O+bz2Z3CPQffMwJiFBatt0LHb\n32/DgWhKxCv3w3yZgwwUF6OJoSWtykD7/v7OvH4sZs4EFCAi1vAYCtkQmudNTjQAzrlf9Dzv7530\n2vMo5P7h1oefBuWE5wS9c6OhN1xfn75vfNwoGOAoAQUzN2fcKIT32awqG7xPYG/QKuB1LiwYPJGw\nHCXOwzE6ah5VPK6fZz4r4a2IpT2A/I2O6rofPbI6A9EBUEG8QfZIxDxzCuesD1IyPG/gnjyURA9g\n2Uk7QLbHww9iBKPLfFqMEVEC6QFyuOToKXD7kVqkCbh+KNe9PcsbQ0PBedH0RM6fwns0aqkpZjOQ\nLmLdFIcnJjTXPTpqKTRqDCJ2jaDxprhKREmXeC6nymx52WCfjx/bvGX++buXKbiCronFjB5bRM+r\n0bApcxh1aBhQqv65w6RvKCDv7RkMc3hYPXKU7uysoWvKZTUk1aoWunFk6A3BqSI16ad5wDni2uOo\ncY8Dm2UeAes5jAoCPijAHKRKGX0pontw+/b56JWrJN1EAD8rIgeV/c8d8tpzKTRz8eAzupGCEg92\nKKQPNAoBfhcUACiPclkHfHz8sXXL0g7PHNdKRfOfpBLIW29vmzdKHjOTMVpdcqngy8PhzoIm3i8D\ntEknRKOWv6d/4WAEwZg/IJEMT+Hho1iNsSiV9IGhBkI6inVTJOa90D/70zl055IzB6KXTJqXB6Eb\nChwWS+CrcChBZ5BIqMID+ghdBMXCjz+2tB00EaCDoBuAewelwzqBsULTDZsmA99p+gOeCkqGyOTG\nDYOzQiHxxS/q9VpY6Jy5ALMsXepEm0Qv4NpHRw0ZRHqRiBBjTwNWMmnR3927RkzHvQd0mTWGwzaD\nIZfT75ubM6pyxkNClbCwoPsbjaqhIoouFvU76YIm/Tc3p8ejFwAHbG3N6LGJkkAjwVmE8RCxmpSf\np6dYNMoO6jUH2T6J2nh+nkc5rg/gL4uOgrx9gBZiQES2er2wqySE3HiF3JQUjUBD0ECzsaE3DAVK\nP4cOoT6Im+1tveFBmoAJ39qy2bSkEFZXjX0TT4ViMt4yuUwQKuTz/XnwYlEfxuVlSyGBw0axgTSC\ntx3SO76/3dbIhjwu3403TMMQSBw8LJAmeGGkwTh3CsEU/cizf+5zmoOnAEt6SsQKjO22rpVUUbms\nygKFTf4a7528vIh5/bBtEkmRI+Y+mJnRY7MHnEsiYdBKf42CdAv0FUQSzEEulYxeA6I8MPoDA5r6\n8xfIx8aMSoMCNp24FFT99Ag7O1Zbwvhzv+7vq5FG4T5+bHvneRp9UGPgmuIEUb8hbYaBHh3VazY7\nq+sDMUTKjEiTaWlQhMNlxDyMDz+0hraBAU03PXig58fzQXMjXew0CqLIMVSIfzYAKVEI8A4jcGM2\nssjp5gVfJzkuAviaKPY/LyL/rWgPgIhIUUS+0+N1nVl4+JwzzPJ5COmGVMoKTLTu4+WjMOCFoTlH\nxCZzLS5atzAPFtOP6AIlKqBuMD0t8s1vmvcGYRq5fLxcf6s//PzUCjAANNu020a+BkcQBcxq1Qqy\nRAgcq9nshH2y51AR0AzEzF3w2KwDUr1q1YaCI+S6UaQoKed0v2BKLZeNxx7FDrc+g+sZTE9HKzw2\npKcoKNNTQRGfvPTeniFMWJuIvs4aq1UjeAONBfyXGQnMD2B9RFAMTkEZ4kxQY6LpiyY9lPLQkNF4\n7+2Z4QIC2WoZudvIiKUrQcDgeOC4FAr6ncCYQUGhlGMxNRLFoo4fXV62FBBoJOY1QylCnwmRHdd1\netroqZl1/cknxqHPM0R9iBGQjIqkkY5UKkYpErH108QJIsvffAlZHLBfEdsz7hf/TIGbNy3iu+po\noLPKcX0AcyIyJyI/eHHLeXZBMYgYP8x5CMiJdNpGI+LV0rzCw03Y/8kn+nCTM2Z9NJahdPGMwL6T\n5wUSyGfxkiMR6+7M5fShBLY3NGQPEUoXdAaoGrpuecjJq0IR4OdsYS08BHi9NLeBQEFR8neOT/qH\nYmc4bA1vwE9FLEdLCknEcPMYAlJgGLZ2247LfsEVg8KFSE7EvhtDV6+rkQKDThogmzUPHu756Wkb\nAeiHHDIjIhIxbieoLDC02ayNWGRv4U+ik5yCM0Po19ZUWW1u2oxfoinSGVByMOyHxkRqL7Ozdm3h\no2LPqStQ1CYSo2eDlF8yaffA4KDl41k/3bibm9bMhQGhAA8cud3W2sft22o4we5vbVnajLnDQFy5\nT4haGJQzN6ffOzam6wqFdA+2tixFS3+AiHUS08PCfcszRa8GNQQMvd8gIAA3gLfeuHF9DUQ3ReCf\nFpG/LSKjYlGA53leupcLO6v42fzOwux3lESjRlDmRwMwWJ1IQEQVhXOaH8UY0dhF6H/3rr6O90GB\nKZPphOuRQnjtNbuxZ2Y0JeKHf5JSgSPlxg3D7wN9QwGzfoqKIlYcpSBGLpkGGSiHRQzJk89rThfu\nmnbb+Pr7+vQhIdVCCg24Hik0Jkr5Sc/8o/tIg83PmzFjDSh4FDOoJv8wdiCzXB94d4pFg26SPwbu\nyEQqwn9qOsB4KXAztGTraUKU+cLUIyg49vercvrCF/T4KKX+fi0Yf/e7RosBMomi8sqKfR/RJayx\nU1P6dwrceMEPHxrHEZFaoWD3BRPNqMsw02BszCCZMzNGR87YyHpd6SGon1QqivjBwPGc0AtCtzQR\nT6GginVyUtdRKOj6k0m9/2/csPRVtWpzBXZ2LHolxQclBtGAP38P4yjGnJSev7mLyWWhkDV8ihg9\nCoYAI+HPJBSLVuOgs/nGjZN1CN9PJJpInMxT1GvpJkHyayLy73ie90GvF3Meksl0hn3nKX4l7xe8\n5VDIDAQpDEJvWtTDYfXM7t61ZpqZGWvZHx7WBxLiNuCCf/InehPOzKgXBevi3p7RRpNWCIf1e8D+\nozj9XiR0B+RyyVdTYAUBA10yOWrqBJ6nv09NaVqLIuLgoKWI/C37XJdw2HoK8O5ETBH710Rdg5TS\n7q4VjhnswhSr8XGLxogMUADJpHXPEk35eyygzabhC6ZUUj1AAJ1TpYpBwUBiTOJxVaJAEgcG1Dgx\n02FpyZrEKEaL6HqZjbC1Zcgu+kWoXVCP8Hc8wxj68svm2dKgd+eOft/HH9vAIaLYV1/VddP5Clc/\n1NGvvqqfg6WV+5PJafSRJBL6fbdu6c/j4/odUJXcuKGRMCNEcY7W1vS9DJL//Of1WvonedGVzjOH\nMeV3UkY8myL2zAGPZe/otzn43PIZfxqTdCoILZy4btLJODxHkcLBRcXPOBWXJd0YgNXrovxFdDPH\nxi72mFtbRic8NGRTl7a39YEbG1OFsrhonuGtWzZsQ0Rfo2WewlmrpQ8KOc2REXsQmWkLxHFvz9AM\nhNzptD6Ia2v2N0JxQnvSNSKdSoyUCzhqcr5453h5GDG8MYrJnqfKIhazYfAYR2CfpDIojoqYRzo+\nbs1DPMB0rWJcSXeRA2afHz3SfRoe1nXzgBMVQEPNXqDESyVLgzH6EeQQaKflZV0rw8Dxknmo+/r0\nH0qDhsFEwozOrVsWEW1sGCkd150RliCt8MyHh0W+7/vUQ4WNlbRIOKzfA6Qzn7fzxsuEpZRUGIaN\nWgd7u7mpRVjI45hzkE7rsXd2rCltdFSdktdfN4UciahBoj5C30W9bte+VlMDw3WEZK7Vsnt7fd06\nnUmh4Y1T6wO1BnEeNSOMIfDT45xBkH48AxgD/j18aPs1NqZryed1D4gWRTo72oeGDlfsBw3DZbOH\nHocC+umnP77zlLP/D0TkKeZDPM/zfr/Xi7sOAj0BDWDc7LTOE8qn03qDz8/r32j6CYX0AXrlFUuF\nQOS2tqbfy9g88uhAK8mzg5TwpwhIi0D0Rf8CP+PRgi7hZo9GLVIAuUSjVK1myKVmUz1M1pHN2uAZ\nitgUk+muvXPHvHPQIAMDtncodigkpqZ0n6JRU9qE8hS2MWr7+7pvGBIKjwwjgfYCUjAahVgLjU4Y\nRyic4fSPRNSz5hoDuaXgmUzqHty5Y+iVzU0zJEQTr7xi0SBrZSAMHbY08o2P6zGZejU1Zflpoj/S\nCBT3MdI7O/r666/r35aWjKwOxYRRZ08bDU3fMUSISM7fvyFi9AoUUEslG0wP/BLyO88z54Yh96GQ\nPgsLC9a4t7RkUNZaTY0ZUSEFXyKmZNKmkdExv7Nj8zZE9L7HuSIa66aRC+NBpAyKijRWJqPHHx21\nUawIgAYij8MMAI4XKaCD0NOLluMigJ8S+XTwS0VEfuLA318YA1AoWD7zoPjJ4UT0gQXKSQEJZQRc\nDU+JEHN/X+Sdd6xItburyI+lJZukRDeun2KY3gFmDZCegDdnb888brDvdJr6Q2uMgogN76Y5qVo1\n3D7wR/DgYLH7+42LiMIpKSA6hCHJoymNod4YGYwQ/PWtluaXITrju7e3LSLAW2M4ysFU2/Cw1l1S\nKUMBQVpGSoXmMdJMIK9AMWGc/PBRYJ1376pyg1KYh3tzU6M9hvGQc+deeu89VYbNpnq6wHVptiIn\nTxGTbmv4qYaHrWnq9u1ONJiI7iHHHRrSdM6NG9YUt7Wlyot0HiyiRA0gxEhrQq1AgRvDQ8T65IlR\njAwO6p5HIjYIZ3TU8u7OKZzzo4/0niXFWKvpvtDDwTEw2hD3AauORtWAzM5a2go0lYilkUTM2TlJ\nKPwTKQLXBcKLQ3CUcA/z81FC9HgV5DgU0M9e4DqurHz8sWGGwWjTfShi6QeKixsbxnJITpSOVxpL\ngOLB5Q5cEo/nIMFavW6D4f1MpO22MWeSe0XxMwoxmeyEREIuBxwOzx8jVa2qUmHmAROxCgVriqLB\nLJPRv7/5psgf/7GeBzxDfL8f0jg9bREGnhVdr/7CZChkaB/SFVAIUFTEq6YuIWINWfzv70ZGCeAB\nO2d8QFAdQGEN4R8erL+fAWqOoSEzpLmcURQMDup9srLSiaKhyE1NiBRONmseNbz41AioY+BkAOck\nugKiyeyHZFINSDptMyYADTQaamwgnYP07NVXLRK8dcvoM8bGrEuYVCUGFucDxV0uq9KfnDScPikY\nojZmNczNWcqQ6AveKYrypH1Am42MWNc20TTnT/8HxWsRa4Cjq5kIz28UDhOMNM8etCMiNjjpuO8A\naUTt6jpINyigvy8aCXyKABLl8X/H87x/1sO1XboQtovoTTE3Z+3m5FKBfPobs0gRTExYPWJ01Jqr\naNwi35nJqLcPHh94JukGPECUGIVA5q7SUwDmmsIda8M79BdzCWlhPuS1YlG919FRIxsDhglMkjA5\nkbDi4/d/vyptRimiSIgsMJp4lvRV+Psn4PjBq0cBMliE38kDY0ihWCAy8isFvGZGIoKS4gH1PEMt\nQeTHYJe+Pi22o2iIptJp8wbJ1aNo8MxHRmycIgRneLW3blkhGedBxNg7y2U9xo0b1jfSaGhtgy5v\nisYYWJT26KgeZ2XF5ihA8IfyxhtPJIyEjWiDjmk8aSIaUk+cP/tBzp1ZCdQciIgowGP8WTf3D/QO\nDBAiDfnGG2bgMYIg2oaHrSeG8aj5vClrhNkRpGhx3o4SEGeAHDCMRLmHZQBwmogmT4Pqge4F1N1l\nSDdF4ISIvCIi/0TUCPy0iDwWkS845/6c53l/tYfru1RBmVQqBjNEIR+EoIFmmJy0Ts2xMb2hNjcN\n0uh5+hm8XTxEcqnI2Jjh7PFyGSsIJw9wvFzOPEuOQ18BipUiLNO8GFlIREFjE15eLKb5YwpgInoe\nNF8x+erjjxXGGI+b0Xj0qBNiSJ6efSP3TwFSxCIVsOikzoAe+onw4EdCiTunr/mJv4CmvvTSZ4ef\nM1iGwjAGhhQFURsFbhSlH0dPrQL4K8gdzg+4JaM/i0Xr7va8zrw7TUowW2JUuG4Q5PX32xyA6Wm7\n/vBJ1evGAYWCLpfVcdnbM2w/8wbicfl0RCkjHV9+Wff05Zdt8hyzE5jJTL2C0aUiBiPGuICkQYly\n/C98wYy5n+QODi2eCUaf4lhxnuD22bdi0fphECJo+nQQmu4QUnp9fVbLa7UMogziDGjwQaHpjOeD\n9XUrpLZELCNw0dKNAfg+Efkhz/OaIiLOuX8gIn8sIj8sIu/3cG1XQu7e1YcAzDDKHEWJd0Q6oFLR\nG4iL+eCB5eKXlvTGX11V5UtO0+/Ni5hnhUcEKdf0tOWsSXHQlAKtAYoIdkq6UaemzHv2PGtCg2CO\nsYqgc6C/zmZV+VAz8Ec9Kyv6NxQwJG5g5SlyQcvw8KEeC2XnTy2sr9sUNfLu0FKUSrp3KF1a/MH5\nE10A/aU4vr/fWTchPIdlFE4lYIscE0OI0qHwl8+bUaLI5++gXluz7nBqHghwSVJJFL8p/JPWoq/E\nn4fHSHOv8f6JiU7qinzeeI8wro8eGf69VlOjPj9v+XmQT6QoGTAEHLlQMF4m7qO9PYu4mHPdbGoe\nnp4JqELgE6KrORzW95HOxCDm8+YAAXf1K2u6fEHbwdUEZHl31wr4jFSlWZB8vL/gCjoOUAJRHBEK\nUc5xQvSJgPbqRjD+/s9ehnRjADIi0i8iT9tXpF9Ehj3PazrnriVFEoVNYHTHSSRiOWweXLxWcPe0\n1TPohfTA7KzC6aBDYBgI3jHNVeTa+Z0bi6IabIn+oi0pHrwohoV4nkHiKE7Sag+DJB4pxceFBV3v\n3bud3ZgMWGePeJBI0aBgyfWDVsFzBdaIl8fDFgppGgRlz1QqFBeFSBQ5ShyPDmWUTut5k6qiD4Ac\nbqul3+1nkKRAjtePh+jvnHZOUzEUeyHRq9ctDUDnLfl2FCmRGemcel33MZ+3cwABBFEazWG5nFEo\nQFu8uKjHItfONcQp2NszplcQYOGwNaFR9PX3XfT3W8GTaWsorsePNWri3iEaJF9PgZrIk3tORI8x\nMmJoKsj6MLoMSALeC1yY55BU1kHFy/GWlozXCaI5EUtXUc9iPTgEpDf9+XvmGTca6pCRLoOplujl\nOCHSxjk6jffP+4EPX1YvQLeNYN92zv2rp7//qIj81865PhH5as9W1iNpNq2zkAfyuOEPNNwUCqbg\nJibsRqNQuL+v//C6V1YsdIWka2pKbxi/YiDdkctpiP3ggf4PBC8et3GPsEJSgKXQ6k+34FnxXtJH\n1arIl7+sr3FejM6rVjVlATKEKGBy0oZ0A8sj/w6aiBuYfRkasi5ZDBiGZmPDYKN0T7MXRC8TE4Ya\nYWQhOXEeNDxMhp1gzCgMkpNtt23cJl4u076gaMDY4Q1GIpbqwhiDpgEeSF1iakqv94MHnUVKcuDf\n+pZ57E+eWH0IQ/fyy1bzoQkKT3loSBUx64HyY2zMIjeYYym4fvSRDeKpVg0CSnc26Z+NDYs48nk1\nxtRHSF8CP/bzNzEjACMJrh+aiGZTrxdzA0SMZgOnhHuPyIm9xEgiNOhh/BltCrIunzd+LQw1jWwo\nZZ7twxq4eC4wONT7cCyOUv7VqjVNptNmtDEyp5HhYUtVnVSg7pWcuGTP8/6hc+7/FpEfEC0A/7Ln\nectP//w3erm4XoifP0akkxjqMNnYsOYnmou42P39xow4OmozePFsw2FTUChYFN3ysnb30nH6wQf6\ncLZa1ja/v2/omJs3dS0bGzYEfmjIPDVQJc2mKhuofwmFGw0llIN6GuNEuoMCJoiZSkXz+zwg8K8z\nbIV8vr+uAO8RSJtGQyMgogKikGZT5N13zQDxAMBeCe1BKKTHw1skxCeiQBlub1utJBKxNn9SVlwz\n9iKXM4QUdYupKd0LqIahJcbjJq1GXQjsPj0aDHx/802jiKjXjRTQOV3H1pYxUFJQTaX02O+8YyMb\n8byZn7y1pWv0Q45ZAxHd/Lz+fvu2GRcMBR4vsFoKyS+/rO9nZGk0qvcmtSZSTBDMMQ+A/cIhoEBO\nvWBvT993+7bVgLiHqCHg5JTLn23U4m8MyyFHj0HwXxvSTuwLe8u6DhOMDWulpgD77FG6A+cR6DSG\n8axy2RxCxzWCveZ53gfOuS+JKv6nvqOMOefGPM/71oWs0Cd+FsN0+mybR8gNVvw475/wlLwvaQNk\ncNCGWDN7lXxzNmvzAl591VrU4amfntaHemvLKKTpTN3c1Ie51dLvZSgNig68Md3GIvq5atX6Cngo\nWAPjKv3EYcWikZkxrJwmLjw9JntR03j0yNJYEJMRljOUBnTKo0dG4bC/bxz2dEQT6pPTxUOHex4P\nEU8Vtk2uSz5v1Nrk7ElZUWiFY2liQhX7N7+px6PQx7XiXgD/zcSqbFbXTFdqKqXXPBq13gIknbau\nZzrAoZ/wdzvDFgsqa3dX6RLoCYAB9sMP9TP0UaTTunYQY87p+eG9oxz9HDO3btn0s0pFo5VaTR0K\n0mdra0btDX1FNqvn5m+0A3YcClldACTU8LBd/3BYj8+sgFTKIoTBQTMq8bg1TxHBNJs2SIlaC8Y1\nHreB7rGYRZqxmK7Vv08nedTw9xC9Uw/C0TjMcBx0Hv35/14KsGnQbueJGDouAvhrIvILIvLrYg1h\nfvlz57eM7oQOWBG90c5K+ZDNWrfrUULHq/+BOmwsnJ/CAUrgalUVDkRazpmXyIB5biYKrv4u2VDI\nMNNgk0nxkAqh8DU0ZFEKNwcFUj4LPI+8Jh42NzCMl+GwsV2SuiINgsKi0E0nLbUL8q0YOLDioKXq\ndZtNCzIIyCR8P4ySFDGkEJ4q+f+REfN6eS+dwiB4XnlFlT1cQEtLaozITbMH0GuI6D6yDyJ6b9G3\nUaupciLVsrio6RWMIoYhHBb53vf0WgLHhKaDpq98Xnlvkkmtuayvq1JeWbHmM+hAuD9guKWGVCxa\nV+ramjWGbW7qumEbbTYtmimXO6G/7GehYEgjrtfEhNVL8JAnJmzYPPUSGiTpa/CnRkZGFMpJTYLr\n5OfzYbiPiK6TZ2t+Xq/h22/rz8mk7umrr+r7gHuijKGdFrHi/EkGgOeb5yqRsDrJUakcaiD+Gs1F\nyEG9559a9qxyXCPYLzz9/975HOp85VmZPo9T/iKW4iBCwCs5TFA8eGDj45bWyGQMBURI7b+g8bh6\nIkQHjYYqfhpiSIVAO5tKqQKiOYv0BXjzlRW9McnTlkpGRzA0pA8WEQXFZLhwQiH93uFhG+hNwZVG\nKWCmg4Pm2YE8evDAGnUGBtTLXFnRvX7lFe0VYFAIM4gxaKSTQFKABKFZDkheLmeKgtoChGiE9R98\nYPwyIIoWFvSY4+PqFYfDqqRBqxSL+jfWXqmYlwqlMZTbsZjId75j5H/cJ/AgwXj5p/6Uru3tt/X1\nV1/VfSLVRw0FvDq5dorgkMy9/77VIPzeLYoIGCqIqGxWz31+Xn/GgMVils6k3pLP27xpz7Pz53z2\n9/U6ZjJ6LIawpFIWqY2M6D22sKDfPzlpiDWk1bJidH+/FfBJqYCKQoikMRbb22aYCgWjh1hdtc9h\nTLrRDRhlET3H8fHuCr+Dg8fXCK6bdNMI1icaDUx5nvcLzrm7IvKK53n/Z89Xd0AgSGu3z5/p86D4\nL/BxFxs0RbutD/bamr4OlwrQP7DGpEHgjI/HDbkjYlA/EfN8mMrEkBY4hBh4027r+95/32bd+onW\nSJllMvow9/dbg4+/QalatUIoypmegRs3bPhHJGI87KWS1QCoIVAQT6dVyWazIj/+43qMN97QB3lp\nyYqBH32kygFEDfUSqHwhPwNSCH3A2lpnfwWzAF55Rfd5fNyKuwdJ51CGi4tGbbGyYogezp2IZmrK\nri1snByXqAuOF4am3L2r58bAe4aaUCz+4ANLh6CAqKEQDcIrBbLszh07NumyJ0/0mkxO6t7evGkR\nloghckDzjI+bA8Q1I50yNaXvL5etlgRy6LXXzOstl/VYoNkgS9vft+fBn0ohLUk6iR4VOn6ZZ+yv\n7UBm55yui8hRxNKbUIxznTBuCDW+4xw+nu9ulfpFK3/0Xqt1/sanm7r1b4nIN0Xkzz79fVlE/qmI\nnNkAOOcyIvIbIvKGaHrp5z3Pe/ukz/nbvXstjFeEFOywTaeTlyLcd7+rShXY5fi4oWIYJPPSS6ps\nZmeNkpmbFKoGvCFgcYS0wPYo3oJaoEg3N2eDPkCiUCibnDTaafKseKwiFqWAjqAJJ5fT/8kNDw9b\nGoWaCAVP+hJgpKSZaWhIUyN37+rnKbR99JHuBWE3BV8Rg0VGIuo5Ly2psiZ1QFFXxPDkpKfg6yGS\nQfmABopG9Zo9eGA1GNhEPc+GmvtpPsi3893w4oAAwbCg8ECxrK/bNWPwy82bRueNMZ+eNrZX0hnp\ntJ4zStLzdI2rq50UEzMz+pqIKs3VVb03bt7Uda+v6/6tr9v5w1oL5z4F4uVlg0ySn2eWBM2C8ACB\njIPqmp6SUkmjj81Nvd/h7qdpKpGwGgO/Dw+r4UbabU2zUXN6+WVLabLPfsXNfeUXom3SkxhAWHkD\n3gAAIABJREFUIlj6BK46bUM43Nnodp7SjQGY8TzvLzrn/pKIiOd5JffsJujvicj/5Xnev+eci4jI\nIX12lys0Z530Hop0hYIqdZQC4To5/t1dfWi/9z3r1Lx503LvTGJ6/XUtCALnJMxG6WIMbtxQJA1G\nZGHBHpBKRY/9xhv6AMAxTxpqYcGGXpDXhvAKrh9SBTw84M79nZIQZ3HO1arRBvuRUNQPCgXr5qUu\nAD8PBUXG+OG1gxyiEWtx0XK+pDoogMMxw55sb+se/8RPqMIl7/+Nb2hT2vKydaSCzuKYyaTlvvEu\nMSSk1+hWprENVBMFS1J97BsMqN/8pnrXiYSujy5mJlrRJc0sA3D9k5N6/5BHB5mTSKh3nkoZHw18\nSBT8aURrt9XwDgzY3lLUDYft3gTEQHRKRAHJHhQOpKlmZqy/AfRVs6kG79Yte2aYIUAdTMTSnf7c\nOwiwl17S3ylAM7wFTxgCvMOUuJ+dc33dom2irOd1zu9ppBsDUHPOfVp3ds7NiEjtrAd0zg2KyI94\nnvczIiJPO4x3z/p9lykUvfz0C0Am6brFI6WfYHbWcOrksPFSKPAR3lOHyGYNx72+rsqE9BKpCxHL\nhff32yxTP3R1edkKb3jvcPyQzkCJ4+kyVGZ01EjEIhFriKM4WSoZfp+cOHsB2drKihmh+XldO/0P\nNJT5seUoBaB2QPvghndODWEup+fFHN+BAVX+X/+6HudP/2mlIKDQSkoDRBeNflB6oIgw8CjvfF6/\nd3vbpsNlMlrUBUq4vq7fzWjHmzeNxRTvm30AxTY9bSk4WGeJhj73Od1HkFQUisNhG93YbquxHxy0\nYmo2a9QNNKhBZdJoWDqFiHN42Iq7ABlGRy3lQ1qy1TIqc7p1ORfmATPInuNCS47Cx9BSB4NOg74X\n7mWUPQZiauqzTtlxThpOhkgnGRww0tNKrWZQ5OdFujEAb4nI/yMiN51zvyMiPyQiP/sMx7wtIhvO\nud8SkS+Ippd+0fO88vEfu1oCVFPEumhBQ1APuH1bH3Y6Mjc3DYnCUG8GcoPBJizGGwNlQb4fBMvi\nouVZgan29+vPjKKkjb+vT/PNtZq+h6IeHimFZrxw0jt+2gI41+NxPQf6IjAydNfChEhahk7QUsm6\nU3d31QOH/TEUsgIf9QMU4cCA7iG8NPyt2VQF55wqJPh/Zmf1d4rh5N9ff93GbZKWErFCLtQQKKbB\nQSOVA9bKYJylJZs9MD6u+z0wYIXkeNyUZDarKSy85b09QxeR9qpUdN2Tk4ZsY6bDzo4ea2XF6i0Y\nEiIHumjbbV1bpaLfRW8EyCYKzMCP2+3OPD5oLBBgoM38vRjr67ov3/62NZ3BQyRixXxYWKHI4DpT\niyGtChPt1pY5S6TFaASDBBAUVLdCNESTGKm10zZsidiwFxGjyngepJtGsH/hnPuW2HD4X/Q8b+MZ\nj/mmiPwVz/O+4Zz7uyLySyLyK8/wnRcu3AwiepOPj+tNsbdn/Owo6+1tVSDc+BCa8R66F0WsgxEP\nF88/lzOsvX+YO3WBeFyVHCkcujPJBwPPg5oXJT09bZwk5EzxehMJDcFZi3PWuEaqo15XgwKrJK/5\nO3MHBw1nT6crqJO9PT0uyoS0EIiYUEjzwAyEyed1rxcWVMmRm6ZAjPdKHYMCNwoGA8dMY1AyDx4Y\n2mp01PaFcwKJND5uxgukF4VPEEvUQIAVvvqqne/srB6TXoRwWBUkCKzJSds32EDBqK+sqIHI5US+\n+EW9j0Bdcb4UVCnKU3y+eVOPS0qGecesl7y+iKHDaKTyM5YSFaLkofygXkPdCdAAxXvn9O90PItY\ndMH+YRzhOUqnLaIKh3XtRNvdCGkkBOjzWeCbdLyzP8+LdIMC+l9F5F+JyB95nvfhORxzUUQWPc/7\nxtPf/6moAeiQt95669Of7927J/fu3TuHQz+bgGAgXw7Kghb0WOyzjH7AFhkyPTpqHnNfnz5co6Oq\nDFZXjZcllTIcdTarioduxeVl/c7xcfteMMoQ0X34oSq8sTEdI/jggXGd3L1rKSpyv8D8IhH5dFA3\nQ8bhsdnYsCIhrKZ4QtPTxoSJ15bJGCIGSgnWCs0vXaXsBcU6FB5GY3PTkE2gXcj/glzBUAFRBEuP\n4dvbs65k8sHMYWDE5cyMGj1m8YK5B8U1MKBIHBEz6iKd83Yp6sNyOT5utRs/dxTMqRSwUU7JpB57\na0vfA7AAtBEkcyCCMGJPnujrYPhRlBhV7h9qKkwO83fD+7muMNYUceG5KpUMYkxaB6QUir7ZNEJC\nvxEhCvSLPx3EdSEKprGOofZbW2Y4Tyt+xU+HNs/uSUIUKHJ5tA2Hyf379+X+/ftn/rzzTgDNOuf+\nDRH5EVH2z5dE5Fvy/7d3rsFxXud9/z8gSBAEeANIkKJIW1JkV7YTRVZk14mTmnbSGU/GTZtpJ42n\nTdIkTZqZ1ta4rSdxptPyS+02bWJn2mk+xJbjS5tM6rRO0k5Ve2QzdlPVcqq7JVpSRcmiZfECgjcA\nJAji9MPZn5+D1e5iF1hg98U+vxkOFwvsu+c97/uey3P5P3ky+Niqv9TsK5L+fkrpGTM7Jmk0pfSr\nxe/TSu3qBaTxS54Cj84JDw0F1rHRwsJCjtIh5Z649UOHfFv+4ov5dwxGt9/uZiHMPDgzGSBPnPCV\nzs03S29/e37/uefyIHLrrfk7pqezA/HQIa/6RcQO2vUTE+4YRQoB/8ALL+RJhExSdhJveUs+P6QD\nkJ/AvMXD+uyzrmiJ5v3Zs3mS+8EfzMdA1pcB6eTJ3Oa5ubxKPnAg2/Kp8/v88164A/MMtuRnn83f\nTbw51wAFSVburJTx5dx5p0eI4I9A14kJDLPG/Hzuu4mJvCtg18TkxgTAhEZVLvwbO3bkvnjiCbeh\nlxLKk5PZv/DYY/lcKeJO5A6TCmG4lP+8ds0d0wyeiLcxgJO4dvKk14w+dMilsElIRCn11CmfjBG0\nQ4WW/Izbb8+fI0lv1y4vFVr6PBpx8aKHQI+MeIlT/F5ExXGflKGerUBcD9+Q5CHVkicDtjOhECZd\nFoTqN8xMKaW2W9eOCehLtQH7HknvkvQrkr5X0qonAEnvk/QfzWybpP8n6efXcKwNo9QWX1x8tR2w\njEXn57J4OKFrmENY+fGwE5XDSokkpqmpfAOfOOGy09ev54GAaKFt23KoXOlklPJDhZb8xISHOkqe\nMHbwoOv7E6NOVjNaNUSGsIravz8/9OjVUPwFmzN5Aq99ra9Yp6fzoE9Ey6FDnsRGW6T8oJ8+7buC\nnTtziCBhkKyQGfguXMgTCZLbt9ySj711q/sMiH8fGXGTB4M0uxyiShg0qTG7b5/H3h8+7Cth6gxM\nT3uxEcxFFD0x8+zUuTmve0tpRUIxkcl+4glP+WcQ3r/fbetzc7kPMMPhIKfkJxFDBw+6I5biRaUW\n0uys5xjQ55gZy5Uyph52v5KLIRKMIC3PziWZj8G3HZNLvc+lHGAnJnxXUK/qWbK46PcwZjGyfaXl\nkW7l+WE2XYle6/asB+2YgB5QDtN8ULkOwD0ppTNr+dKU0mOS3rKWY/QC7MxS8wIRpT4IK1liwkt7\nKatPUtLLUDjSzSmCPTubV6Y4ViU3haSUV4eHD+eHniQtiqCfOeP2aEIVSQgjBJUEJ5KrcGRj02XV\nvXt3HohffNEF18piNUhdUCScXcvCgrfhttty24g+2rPHk+Ukd9hhtrj55tzOQ4c8FBQfA4ljDFwv\nveQRGnfdlc+VGgxEVFE+k0IfRNgQsskq9MyZvKLduzd/nqQ0wn6JbsHZS6w2xXKI5MJBLfn1pQoU\niX1793pGMmYzxAZJ0CJLfO/efA6UWiRMFFlxBk4E6EhCnJ/37yUBjh3RgQPueK3XxqKfMOswiSAV\nTUEk4tTrnwFqL7BLbTXQNhvYmaBxGDdbfZ8/74M7zuuyLfyOKmWcWz+ZdDaadvzhjyuv/r9X0iVJ\nM2b2YEppfl1b1oeUQkyNIgkwPxBmxiSBo1bKD/nkpDs4JU/CwQ6KFAG2YgqEHDzodU7LQht33bV8\ntf7CC/mmnpzMPz/5ZD7uG96QV9wMxqdOeUjeyIjrHlGwnZ0KEwS2emoVP/+8h6mSuVsmPZGP8K1v\n5UHt1CnP9KVu8OKi50hQeo+qUtRKftOb3C+CJAOTJyGXZQ1jQlN373bF1FK1EVGy8vNS7n8K9RCZ\nIrkEhuQOUK7b7Kz/nozdUi9qfNzrKbOqZ2WPzj8ZtAgCUlqQnRX9UmrzsyJmhctuCF0mdHrIkZic\ndCmL6WmXT6CmAeqomMeAjNzS1MN9yf1aJmARCk1NarPlZRJXGz1Tqn02oxzsWZwQ/UP4MW2cmvLq\ndoNMOyagD0iSme1UDv/8pKSDktq0wm0uGPgRDiNTGFMEJhwcfpKv8pGORs+FEERWdAwUOC8JzyO+\nnHR/JgxkcoeH3aGHUBtZuw8/nB9QJCOee8416xcX3bZ+882e/DQzkwcjBNZYKRLRceRInmQwp1y5\nkicnFCGJiUd9U/KBimNgxiK+m0xnBuPDh930VQ6opaQwUVXnz7u9+tIlj9HH3k6+AzHv9N+RIy43\n8K1v+UDBzkTK/USiHUqgiKpRb1la7vAsKaUJ+N3Onb6Dwea/a5fr35CohKkJuW+SsMqB8MCBfG9Q\noAe9KAZDIpsIsWQSLbNoEf1DQba+BCIicORFpOQLEXwESGVLuV8oPnPmjGcPr7dbj7wY+ljyHRzR\nYnv35nPGUT7otGMCep+yE/gHlGsB3yfpq+vcrr7n8mWPCsBezg1OqKHkqzX+sS3HCcz709O+2mLF\ndMcd+eEhN+DcOR9ISWfH7MKAjuAWDx8KjJgjxsbc+WvmNlBW1aw4jxzx8FIG5ZTcwYmIHbuEgwfd\nrj82lmsdSJ5Yhk0Y5+W1a3nAwnGKLDJaR4h9MQGdP+/KlUQrnTmzXFJgbMx9Ek8+md97zWvczHL+\nfP4OopiYTAjJXFrK54CTnGtx9uzybF/KBU5NeX+TFdwIrjtBBNu3LxcXZJWekq+mEX9LyUuSnjuX\n23LgwPLQXOr6MrEePJjNdNPT+b3bb3e/jpSvLyYb1CVRiy0ptXxKRU/ksakV8dxzuQ1TU/n6l3+P\nD2zLlvYlF86f9/wUTG/NKOVA8FXxs+S5JlL+f26us1yCkmYy0VWmHRPQdmVJ6IdTStdX+uNBodFN\niT5KSi6URpz79u3Z/k1MPitMwvOQAWaVSqTCqVP59YsvemIUZhhWl2TRSvkGn5nxFeHFizkiRsoP\nJ6s/ahhMTuaBRXJnIiqmU1M5iujKlRy/jkbO6GienE6edKnsl15yZUtMFlQ8O3Ag+w4uXPAiOg8+\nmL/zlls8UQm/xfbtbrZhZczDOz3tYZEHDy7Xq5+dzWGc3/lOjh6anc3/3313PpcdOzyZrTShUQgG\nqQOkLoj4wfzEgAvY8LHX13PjRr5+1DAoi78wESAfgVOV8Nx6Wzwr9KEhr+lLcXsp9wE7QhLmJNfU\n373b/THXr+fzHBrKEyT28G3b3K9w7pz7bspBmPyAI0fy9X/6aY8II/KHMqbXrrlceruQYInYIe1h\n1V7fJ0iCo05b/1w2mzyIdGsnKYza39z7pUmr6rRjAvo3G9GQqrFzp2dLsj3ets315SVfmVKQfM+e\nPAnU5wqcO+fhpRSrRrSKyBoGQRKBzPxBxylH4RcKZV+65LsAEp8wQVH4pNTfwR7PwHvwoGvdnD7t\n2cUkE919t1ceY4U7M5MHe4TIcLqy4yFMsdSVx1eB6QuNe/R0yDCmLCA+hoUFrxhFX1+65OUryYRm\nJT02lq8BfUbeAvUOWBleuZIHOFasExN+jgxm7GwoeIM8BHUFkP9gh4DmPgM7AxN5BmRMkzQGfCeD\nMZE3ONvvuCNPos884zsLsqlZrV+96mG+S0suHIcd/8CB/A9KfxQif7SbSDb8Ozi1cdAySK9mkMS0\nRNsI/SSyh3sCKCTE60bFncjNwTTKLpcBnYi7VszO+gILierN4jtYRVJ0IL1aLI7JoIxSGBryknfE\nUEs+8LL6qBfBIiIEtdH9+/ONygPAVp8VMdEy3OAp5QFp926PAWfFx8NAoe16eWgeeOQRSJTBxETi\nGeGVMzMe9XPbbe4rQG5h714XFZuZ8X9bt3omaTnAXr/uxeZffjkPbrfe6pr8SP7edJPbqjGTYPue\nmsqT1P79PqgTwnryZP7+S5fy9774otuzJc/KJpmNENt6pUmigebmXBCtVCLF2Xr6tBdhwURFCKrk\nmd3S8pDRktHRPOFeuuSZ3WW9h4mJrBnEahg1WJRKS9MLK+kygavRvV1S/g1Jg/QPZS5JQlxLkRQk\nSHbt8rrUZUg1kWNlW+pVWBuxc+fyRRcRTJIvSFrlFZTn32gnUmViAugCRKYQeUBM+cSElwRE0//S\nJc8OJr6ZiBRWVRRqTykPfmShkvGJgiXKlOVKhpXf93yPm5tQAyW9HiVNzBsHDrgA2fXrvsqlmtT4\nuNvIpfz/zEyepJAARnMFQTjKF5LJTH1XdO3PnMk7hZtvzjZkzGWUrMShffq0V/wiVh2n79135+MT\n8YOE9RvfmM+nrEmwZ49fCxz4TJ5EiMzP5wmn1PZnkGEXIrlZ79o1d5Sz+p2edgVYAgSorTA15dca\nezK5BvhxWpkkUFmVPIKpDJ8dGXEfyuSkiw6WK2Puy4sXfUfCZD887IlYmDDLjGIYHc39hGzIwYPd\ni5EncY6+wdcgvbpv6tVB24U+KzWZWsGuoaz8tlmICaALkCEoLY8TRxgNM8WWLX4jUc2Ih2t8PNvo\n5+byw8mDvmdPXu1euCA9/ngeGPfsyQ/g6Gh+mJEPIENyaCiHfH796/nhRtQNmYnxcQ+XvP12zyjF\njkokT5k5uWtX/q6zZ92sceiQO3jvvNNr5ZIUhkY78gCEmd644QlVmK4IQ8RBi+Mac8PQkIeO3nab\nZ7USoXT4sO8Q0DjascMzQFlBXr68XHCPwYUV+tKS7+SoxYrzm3NhBUmB9XqNeWm5nhGTK3WZ2WWQ\nFEYRGjSSKJTTaKW5sOCJYJwbzmV2btQ2blY7g4iwpaV8PYGQTcTmuFfZTdW3p5s6+pj02OkCixp2\nryXsdDqFxQGRXBcurGwGWuu5MkaQb9EvxATQBcoVEjb78jUx3KwMWUEQCin5Vr0Mwbtxw2P/L1/O\nDyJJTS+/nFe2i4tuN8W5vLCQBwWOjY7O1avuC2CSYeXGzfm61+XPnznjkU4MJDgymeRIHBobc9E4\nyXVWkA2Q8gC4fbu35XWvy8cnbh1zy8REHvzIZmZ1fP68Z6QSIUO/kxSFzwTZXrJgiSTatStPFJjf\nkETYts3zIFBFRZKCKC1MX4jNzcy4RtD+/d4POLm//W33C+3b585XfEHUoGXlzQSBQ3x21iccYCfJ\nPTcxkfsaqWnOm1DPZnWvMVFJuT3IRHDuW7d6IXmynSkJWoKfhXDYZg5XdrfUAqgH82OZ+Uy7y4CJ\nbrJli4fwlua+9aD+um1UUat2iAmgCxC2iAmifOgo6oGwGqab4WF3IEo+aRD6hvAZK1+iOi5e9EGV\n7Tqyy6ymKfbCQFVuqcfG3CnIahF5CfT0Kf5y+XI+JiGDiIjddFN+oE+c8PKWY2PLV271D9Ply3lg\nLG3f1DHGz0A9Wgrck3CE8B6x2wy8yGxLXm2L3QShqufPuwAZyVckNWH6mppyx3G5QiOsEwc5kyom\nJJyQW7e6hs6lS169DME+nNwMsOxakDXA4VpS1iGAUoUSpVNMUeQ7lOaSRoM/RXwIHMA/c/68S3ug\nMkq9ikZOTyLTyC4+fLhx1SrCV9khT0z4JIDzeGnJd05LS76TW0+2bfP+bJVd3A3qr1s/hZPGBNAl\nMK/UUy8KJ/kgiP2V11AmkWE22bs3Z/yePZsHSqI9OBYFU6TlGcWEG5aJZ0gE4zcgIkby1Rr2Tslv\nYOQOrl3Lk8mJE777eOmlvAsgHFHK3zM15Tr4TDBEMTGAE/rKjqasoMU/8gRKjR3OdXQ0TxAUxMH2\nToST5CtLonEYYMjJYHeGCYjEpn37fFXMzgZBMCZ7ko0knyBIsCJZyswXCjt35gmYBDCK8Uh5Z4dq\nZ71tmqIx3CMM9nv3elgskg2NipYwYUjeHoTqSA4sfUqlGmz9gMWORvIM+EYwwJefw8zF4qe8Ho18\nDusBqqfsmNaT+uvWL4O/FBNAT0H8iteNoDAL8eH79nl5yNe8xp1k27blmwxxM1azU1Nuj967N6/G\nTp9enhHKDYqIGdEqOGDLtjEgS3mQomjHyIgXAQEGU1Z0rHaxNZvlAY9wTITeSEBjwhodzXb/Cxd8\nJ8MkRaEQNITKYin1/UoGMdFV/J5kNim3tSwkgqIlkg9DQ24aIvmI3QbHGx93J/uhQ+7cJWu7kVNV\ncgc++Qb1Zgl8GjjvyyQnpCO4ZxpRiqBJXr8iJRezk/w6MQk0Wo1TDYz7rX6RA6WOENeIfuZ/6kHs\n3On3+0bQzcpeOMTLwvVQf936iRXloHtBv8pBrwes1pEXQPunfiDlQcOejR0b9UXkFyQfuMh0Rc74\n0iVP3SfBh6LvZSKSlM1OrLjHxxtv71PKUhM4Mjkephu+A1sxUVIUOZ+f9+QvCtjcdNNyyQEGQOoA\nM2kxiFLg/Ny5/JnFRU/kYkJr1e9lP+OfKUXFWLXz90QFXb3q4Zj79i0ftK5f9wI2TDDtMDPjkVc7\nduQJfqUIlfrr1gpMgkyqhPbWf557qgxVbvbd+DBaDWwcrzRL4adCHoUiNmUBl6pQVgsjEq5Xq/yu\ny0EH6wsDEGGPUh5AKA0ouT4NlK8J/6u3YZbp7keO+ANNYlop6SAtX9lzXFZ1ZcGQEm72V15xzR5s\n+5hwWPUhTFaufomHJ3v66tWc0ETNgvKckG7+9rez7blU9UQBlYI37Qy4jWQb6ItGiq7l3yMzzG5g\ndnZ5lAgRQc1g8iDblfYuLbkZiryNlai/bq1gR1heh0af555q57vbiftvdDxyKwhllprfZ/1O2W40\nvPrJzNOKmAC6AElOCwt55baakLFy0FnN5gdRNF6XUIwDITQKwLTaAu/a5bb8ZudD9AQVzvAvkHVZ\nhug1UnM0c5nf69ddToJokPpBlOLq5DYQ/klWJ1E4kidoEY7aaFJodt0Qa2vUl3yuvEbltWsHVuGS\nF86hj/bv98lzPTbB7ahqrgTBAcPDue2rHez4PGGea0kig3aue7cp8z4w91SFCjW1f5mddVssmZyN\nVpitQs2IfimLyLQDxxwZ8V1Ds+/YtcuVR1eiPEarB5xBC8coTuFOhLcmJrzoOeqc6OCXMNCyc2Cw\nxDxDRBGRQui9lINsSaPrNjTkUhEI1zU65/HxxpXfOqUc5MlrIDGt/jqtZ6hiuywuurkDX8paYuRJ\nbuvGinlpqb3r3m2Gh/NOuEorf4gJoAus9FAiC0B0RqObctu25WafdqCIO6udlQpbYKscGmpcm7UE\n56CUH6pS4wiIyMHhTIhpp4yM5JU8A399VSogg/Xy5TxYTk76Sl3qfHBs9PcI3klegrHR31GMZjXs\n3u3XrTwG5rNGYCLEh9Iv2ajdmJDWOmiibruwkO+PUmhxI6na4C/FBNAViKVvVmSa0o5S98SkSAiT\n8gNw+XLrbEb+htdXrrgTsJE2evnwtHqQ6pN0xsZcmKsT2V2KdxBy2eyzt966/GdK/BGBgQOZKlr1\ngywQOUUkC9etG4NRs0IjREOVulDtQHEdyRP/6s+JKJtmEwOhvvhlVgsJi5iAumG2WSuE57ITxhSz\nlp3ZoBATQBeor3pUTzmodEtMitBOBrKVVoT1ImBDQ9nGj5YNEgCwd6/HaZemGEwUjeycRHtQdhDt\n/nYmOwZHoo3m5tqLCClF+W7c8LDYss5vPeycJNcegl27XCaBLOJ2KTWh6lUmkWtoZ/fV6BwJC+Xn\nklbXUcqfO3vWa9/WRy21A0EEmL/6YeCH8lwoG7pZ1DrXm5gANgBUDZFH6Mb2He0YYuBXMkcgs3D5\nsj/EZS1gxNYAU0MJg1izgebaNQ9Fpej52Nir1RgbgYZPWVSnU+bn3XyDmFsjGz5+Cl6Xg1m9ymsn\nIMMtucw01dpIpiNpqpNwRxQ4yfSuL9pOEhZyyvXnTL1pKfcPtanbZaXr3mtGRz2BcHQ0Bv9OiAlg\ng1iNvTgl16epf+gYDMoSdytBXD6U4nTIQTSLYKB+AO2i+AfFOiQ3L9Buznl+fuUJgIF3taYFEsMA\nAbpGkJTDa84J2ex2QiCbHRdKsxohnUwOq4kSIfGqHs6Tc2907PK6dPr99dedzPR+I8w9q6NnE4CZ\nbZH0F5JOpZT+Wq/a0c8gCy3lAbE+1pzV9mpXPJOTbj4YHs7hlKUQVwkr9CtX3OnLahazDQk9FHkv\n9Y3aocxL6AQKqsDYWG5fs4GO0EPJJ5pz57w/9+xZnTwAq2pkNJgAMPvQz6stSdiMVoqZkgvHXb3q\nwnftQl4FE2asrjcXvdwB3CvpKUkxdzehFJEqX0teNJxauM3AgVgv6SC545UwSoqZtMoPoMA2f7Ow\nsDz8DWcj6qMUFF9Pyr5plD9QT71jGM2e8nir1Ydpdi1WSg5bC6WyZTOaTa7IZFO9rVFJxX37PLN5\nvXVzugm1H9oNfR5EejIBmNlhST8u6V9K+se9aEMvQM4YJciVGBnxHUCjbXergf/KlTw4k+i0dWsu\nvlI/uCNURvJVs1BS9G+I2SZ5qj5yhhKHRNhs375+4XhU5CozMVda3SKERz1lKbcPlU5UWNthft51\nkLqpK7ORUIBdypNgo2CGteY69IJSnmF+fnnJy8Dp1Q7go5I+KKmLJSX6m1KJ8erV9m7IiYl88yLY\n1i5EuczPS88/n49j5oVjSoiWQfyr2UqpFHMbG8urQgTqgDKVaPtQWH09ILcCuzbFw1sNxGUR8fn5\nHHO/ZUt+j8Iu7a5ymVwlXx2vlIfRj5Q7n/J11SnPhRrNVYzTX282fAIws/dIOpNSesSmfp6ZAAAN\nKElEQVTMjjb7u2PHjn339dGjR3X0aNM/rQSlg7KTG3I1K8uyTirOx1aDcbt6MpSkbBbaSaWzyUlX\nH12v1T87EsCc1YryGmCqQtwNzaJ2B4l6ZU2qiFUNqmPxerOwY4fLk/RbFa5ucvz4cR0/fnzVn99w\nNVAz+7Ckn5G0KGm78i7gj1JKP1v8zaZTA712LW+3KcKyHlWOgApaJHxRX2Bqan11SlJyZ+pq4807\ngSL2Q0P5u9pRzpye9tyJyUkvrE7mc1l3txUUOiFyiopa/czFi27+Ks+RJKoqadi0w40bfq0HhU7V\nQHsqB21m75D0T+ujgDbjBCB5gtFGPGjld7HS3ajU+I38vk6/C52Y+muAH6GTCYtjbWTfrpa5OTd/\nSV79K9hcVFEOevON9E1oJLmwEd+10Su79VKybPZdndBspbuanUq3V804ZOsLo3eD+uuxCddXwSqI\ngjBBVylNI+tt6tpMzM+7U1laXqGrG1D3FxPQIF0XKuj1+y6tG1RxBxBsIq5ccQcpctCRPLQyjeLv\nu338VmKBmxUSHbds8Up0gdPnbqugatQL3w3CqqsbbN/uKrHj4/0pt1A1CEuWlteBDpyYD4NXQWbs\nauSRx8ddX6hq1ZF6AVLWIyNrK6wSvJp6FdXIBn414QMIlnH9erbhp5QH72aSysHaKYu0b9vmmkpB\n97h2zU1AgzDBhg8gWBPU5JVcJXS16phBa+bnva8XFjykNOgeqxUYHBTCBxAso4wNHxoKB+56Uvb1\n8HD/J5IFm48wAQWvgsIurSSVBxkyrKXOq4bVQ7GW0dFY/Qdrp1KZwM2ICSDoZ6iQJeVJsttJW0Gw\nWsIHEGwIi4tZP4dylFVfvV67llf1aOu3WtWXSpP1onCDyNWrWU9peHj9ah4E60NMAMGqOH/eB7+l\npWonGaWUo3Eo2WjWeiAbH/ei8t2u7lU1btzwSCaK3letdsAgExNAsCoYLOtfV5GUOjsfErXMqr/z\nWSspLdcVqvq9MGhE3EGwKnbt8kzfqq/4hoZcHnnLlvbkoIeHY/CXcj9QR6Ddvgv6h3ACB6smpc2V\nJLbZzmcjib7rDzp1AscOIFg1m+2B32zns5FE31WTmACCIAgGlJgAgiAIBpSYAIIgCAaUmACCIAgG\nlJgAgiAIBpSYAIIgCAaUyAQO1gx1gEdHQ3s9CKpET3YAZnbEzL5sZt8wsyfN7P29aEewdmZnsy7O\n3FzWB7pxo9ctCoKgXXq1A7gu6QMppUfNbFzS/zWzL6aUnu5Re4JVUg74aOqEREIQVIOe7ABSSq+k\nlB6tvb4i6WlJh3rRlmBt7NjhA/7ISFQQC4Iq0XMtIDO7RdKfSXpTbTIILaCKESv/IOgPKlUQpmb+\n+Zykexn84dixY999ffToUR09enRD2xa0T8giB0FvOH78uI4fP77qz/dsB2BmWyX9N0n/I6X0sbrf\nxQ4gCIKgQypRE9jMTNKnJE2nlD7Q4PcxAQRBEHRIVSaAH5b0FUmPS6IBH0op3V/7fUwAQRAEHVKJ\nCWAlYgIIgiDonCgIEwRBELRFTABBEAQDSkwAQRAEA0pMAEEQBANKTABBEAQDSkwAQRAEA0pMAEEQ\nBANKFIQJgk3O3Jx0+XLWa9q7N3SbAid2AEGwiVlaki5ezHUbFhbyRBAEEBNAEATBgBITQBBsYoaG\npN27peHhXLBn585etyjoJ0ILKAiCYJMQWkBBEARBW8QEEARBMKDEBBAEQdAGKUmzs9L8fK9b0j0i\nDyAIgqANzp2Trl/PrxcXN4dDPXYAQRAEK5CSD/5SzqnYDMQEEARBsAJm0vbt/nP5uspEGGgQBEGb\nXL2acyu2bet1SxrTaRho+ACCIAjaZLOs/KEnJiAze7eZnTCzZ83sV3vRhiAIgkFnwycAM9si6d9L\nerekN0p6r5m9YaPbsZ4cP368101YE9H+zNJSjvx45RXp0qWuHHJFou97S9Xb3ym92AG8VdJzKaUX\nUkrXJf2BpL/eg3asG1W/iaL9mStXcrTH0lJ+XUaBrBfR972l6u3vlF5MADdLeqn4+VTtvSDoK8xa\n/xwEVacXTuAI7wkqwfh4TvhZXJR27MiKmkGwmdjwMFAze5ukYymld9d+/pCkpZTSvy7+JiaJIAiC\nVdBJGGgvJoBhSd+U9KOSXpb0kKT3ppSe3tCGBEEQDDgbvqlNKS2a2T+S9D8lbZH0iRj8gyAINp6+\nzAQOgiAI1p+eawGZ2X1mdtrMnijemzCzL5rZM2b2BTPb08s2tqJJ+4+Z2Skze6T27929bGMzzOyI\nmX3ZzL5hZk+a2ftr71ei/1u0vyr9v93MvmZmj5rZU2b2kdr7Ven/Zu2vRP9LOS+p1sY/rf1cib6H\nBu3vqO97PgFI+qRyUljJr0n6Ykrp9ZIeqP3crzRqf5L0WymlN9f+3d+DdrXDdUkfSCm9SdLbJP3D\nWlJeVfq/Wfsr0f8ppauS3plSukvSnZLeaWY/rIr0f4v2V6L/a9wr6Sl5dGIl+r6gvv0d9X3PJ4CU\n0lclzdS9/ROSPlV7/SlJf2NDG9UBTdovSX0fNZ5SeiWl9Gjt9RVJTyvnZFSi/1u0X6pA/0tSSmmu\n9nKbsk9sRhXpf6lp+6UK9L+ZHZb045I+Lm9vZfq+SftNHfR9zyeAJhxIKZ2uvT4t6UAvG7NK3mdm\nj5nZJ/p9GylJZnaLpDdL+poq2P9F+/9P7a1K9L+ZDZnZo8r9/OWU0jdUof5v0n6pGv3/UUkflLRU\nvFeZvlfj9id10Pf9OgF8l5oudNU81b8j6VZJd0n6jqTf7G1zWmNm45L+SNK9KaXL5e+q0P+19n9O\nuf1XVKH+Tykt1UwohyX9FTN7Z93v+7r/G7T/qCrQ/2b2HklnUkqPqMmKuZ/7vkX7O+r7fp0ATpvZ\nQUkys5sknelxezoipXQm1VDenr21121qhpltVR78P5NS+nzt7cr0f9H+z9L+KvU/pJQuSvrvkn5A\nFep/KNp/T0X6/4ck/YSZnZT0+5LeZWafUXX6vlH7P91p3/frBPAnkn6u9vrnJH2+xd/2HbUbB35S\n0hPN/raXmJlJ+oSkp1JKHyt+VYn+b9b+CvX/PrboZjYq6a9KekTV6f+G7WcArdGX/Z9S+vWU0pGU\n0q2SflrSl1JKP6OK9H2T9v9sp/d+z9VNzOz3Jb1D0j4ze0nSP5f0ryT9oZn9oqQXJP1U71rYmgbt\n/xeSjprZXcrbx5OS/kEPm9iKt0v6u5IeN7NHau99SNXp/0bt/3VlifEq9P9Nkj5lZkPKi7HPpJQe\nqJ1LFfq/Wfs/XZH+L8HUU5V7v8Tk7f8NM/t+tdn3kQgWBEEwoPSrCSgIgiBYZ2ICCIIgGFBiAgiC\nIBhQYgIIgiAYUGICCIIgGFBiAgiCIBhQYgIIgiaY2feZ2X09+u4HzGxnL747GBxiAgiC5nxQWVtl\nXbBcHrUZfyDpl9bru4NAigkg2ASY2X81s7+wXBTml4r3f9HMvlkrWvK7Zvbvau/vN7PPmdlDtX8/\n1OCYI5LellL6ek3x8hkz21f73ZCZPWdmk82OZWZvNbP/bWYPm9mfm9nra+//PTP7EzN7QNIXzeyg\nmX2lVrzjiZqevpQlCX56fXsuGHR6LgURBF3gF1JKMzU9mofM7HOSRiX9M2WJ6CuSviTp0drf/7ak\nj6aU/tzMXiPpfklvrDvmmyV9U8qKl2b2WUl/p/bZH5P0SEpp2sz+U5NjPS3pR1JKN8zsxyR9WNLf\nKo79fSmlC2b2TyTdn1L6cE3baKz2nadrWjtjKaXZ7nZXEGRiAgg2A/eaGYU7Dkt6vbJOzZ+llC5I\nkpn959r7Uh7A35DHW0nSTjPbURQ3kaTXKsvpwn2S/lh5AvgF5UpwTY8laY+kT5vZ7cq6LOWz9gXa\nJekhSffVVE0/n1J6rPi705KOSDrRdk8EQQfEBBBUmpr+/I8qm2uumtmXJW3Xq3XcS8Esk/SXU0oL\nLQ6dVOisp5ROWa79/C5Jb5H03lbHMrP/IOmBlNJPmtlrJR0vfv3diSal9FUz+xFJ75H0e2b2Wyml\nzzRocxB0nfABBFVnl6SZ2uB/h3Jt4CTp65LeYWZ7as7Wv1l85guS3s8PNeXKel6UdLDuvY9L+qyk\nP0yuolh/rO8v2vVy7fXPN2t8zWx0NqX08drx7y5+fUDSqWafDYK1EhNAUHXulzRsZk9J+oikByUp\npfSyst39IUn/S1ka91LtM++XdI/lsnnfkPTLDY77mKS/VPfenyrb6D9ZvFd/LOR3f0PSR8zsYeVa\nuWXR7nJVf1TSo7W/+ylJH5Okmqb+dNj/g/Uk5KCDTQsO1NoO4L9I+kRK6Y87+PzvSfqdlNLXaj/f\nI+k3U0rvWJcGL//uX5Y0llL66Hp/VzC4xA4g2MwcqxVXeULS850M/jX+raRfkSQz+zXlusMf6m4T\nm/K3Jf3uBn1XMKDEDiAIgmBAiR1AEATBgBITQBAEwYASE0AQBMGAEhNAEATBgBITQBAEwYASE0AQ\nBMGA8v8BNarLmKTJ3tcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2151b0e150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ages = live.agepreg\n",
    "weights = live.totalwgt_lb\n",
    "thinkplot.Scatter(ages, weights, alpha=0.1, s=15)\n",
    "thinkplot.Config(xlabel='age (years)',\n",
    "                 ylabel='weight (lbs)',\n",
    "                 xlim=[10, 45],\n",
    "                 ylim=[0, 15],\n",
    "                 legend=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean of mother's age:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24.978916835186908"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live['agepreg'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean and standard deviation of birthweight:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.482384352881699, 1.1690646813035233)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live['totalwgt_lb'].mean(), live['totalwgt_lb'].std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the coefficient of correlation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.085023137198215998"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "thinkstats2.Corr(ages, weights)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Pandas `corr` function gets the same result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.085023137198215887"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live['totalwgt_lb'].corr(live['agepreg'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the relationship more clearly, we can group mother's age into 3-year bins and plot percentiles of birth weight for each bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEPCAYAAABCyrPIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8XHW5+PHPk5lM9q1Zu6Q7paUtbQGhIEtA4EqvgOwo\nqCwKP6UXfnpduCiXgqJXuXqViz8FQSkoRZGKeuECgoZNFoG2tKUt0C1daLZm3yYz8/z+mCWTZLLP\nTGaS5/165ZUzZ845881J8n3OdxdVxRhjjEkZ7wQYY4xJDBYQjDHGABYQjDHGBFhAMMYYA1hAMMYY\nE2ABwRhjDBDjgCAiN4nIZhHZIiI3RXi/QkSaRGRD4OtbsUyPMcaYgTljdWERWQJ8HvgI0A08LSL/\no6o7+xz6gqqeF6t0GGOMGZ5YlhAWAq+raqeqeoEXgAsjHCcxTIMxxphhimVA2AKcIiJTRCQT+Gdg\nRp9jFDhJRDaJyFMiclQM02OMMWYQMasyUtXtIvJ94FmgDdgA+Poc9jZQrqrtInIO8ASwIFZpMsYY\nMzCJ11xGIvJdoEpVfz7IMbuBY1X1cJ/9NuGSMcaMgqoOu1o+1r2MSgLfZwIXAI/0eb9URCSwfTz+\nAHW434UAVR3067bbbhvymET9Sta0J2u6kzntyZruZE57sqZbdeTP0TGrMgr4vYgU4u9l9CVVbRaR\n6wFU9V7gYuCLIuIB2oHLY5weY4wxA4hpQFDVUyPsuzds+6fAT2OZBmOMMcMzYUYqV1RUjHcSRi1Z\n056s6YbkTXuyphuSN+3Jmu7RiFuj8liIiCZDOo0xJpGICDqCRuVYtyEYY0xUBfqhmD6i8dBsAcEY\nk3SsxqC3aAXJCdOGYIwxZmwsIBhjjAEsIBhjjAmwgGCMMQlo9uzZPP/883H9TAsIxhgTJdu2beOM\nM84gPz+fI444gieeeCL03vPPP8/ChQvJysrijDPOoKqqKvTeVVddxa233trrWiIS9x5VFhCMMSYK\nPB4P559/Pueddx4NDQ3cd999XHnllbz//vvU1dVx4YUXcuedd9LQ0MBxxx3HZZddNt5J7scGphlj\nkkpgsNV4J6OfLVu2cOKJJ9LS0hLa90//9E+ccMIJzJgxg4ceeoiXX34ZgPb2doqKitiwYQMvvPAC\nq1evRkRwuVycccYZ/PGPf2TOnDmsXr2ahx56iL179/Lxj3+ctWvXkpaW1u+zB7onNjDNGDNpffSf\nb4nq9V558rtjOt/n87Flyxaam5s5+uijQ/szMzOZP38+7777Ltdddx2vvvoq5eXl3HHHHaFjVJXH\nHnuMZ555hrS0ND760Y/y4IMPcv31148pTYOxKiNjjImCI488kpKSEu666y66u7t59tlnefHFF2lv\nb6e1tZW8vLxex+fm5vYqTfR9whcRbrzxRsrKyigoKODcc89l48aNMf0ZLCAYY0wUpKam8sQTT/Dk\nk08ydepUfvSjH3HppZcyY8YMsrOzaW5u7nV8U1MTOTk5g16zrKwstJ2RkUFra2tM0h5kVUbGmAlj\nrFU8Y7V06VIqKytDr0866SSuvvpqANauXRva39bWxs6dO1m8eDEwvKkn4tHjyEoIxiSZpuZ2du2p\npqm5PSEbVyezzZs309nZSXt7O//5n/9JdXU1V111FZ/85CfZsmUL69evp7Ozk9tvv53ly5ezYIF/\nCfnS0lJ27do16LXj8buOaQlBRG4CPg8I8AtV/UmEY+4GzsG/YtpVqrohlmkyJpn4fD6q9texeXsV\nW7ZV8c67e6naXxt6Pz3dRWlxPqUleZQV5/u3i/MC+/IpKcwlNdUqAuLl4Ycf5v7776e7u5tTTz2V\nv/zlL6SmplJcXMzjjz/O6tWrufLKK1m5ciWPPvpo6Lxrr72WSy65hIKCAk4//XTWr1/f79rxGJcQ\ns26nIrIEWAd8BP8Smk8D/0dVd4YdswpYraqrROQE4CequjLCtazbqZkUOjrdbHtvP5u37eWdd6vY\nuqOKlpaO0V9QhMKC7ECgCASLkvxA8MijpCiP/LyspJpSOlG7nY6nZOh2uhB4XVU7AUTkBeBC4K6w\nY84D1gKo6usiki8ipapaHcN0mUmis9ONu9tDTnZGwmZ41bWNbA48+W/ZVsX7uz/E5/UNeo7D6aC0\nOJ/6wy10dbkH/wBV6g+3UH+4hXd37It4iMuVGggS/pJFSVFeKHCUFudTWpRHWlrqaH9Ek0RiGRC2\nAHeKyBSgE/hn4I0+x0wHwv9K9wMzAAsIZlCdnW6q65qorW+mpraJ2vomamqbqK5roqauidq6Zppb\n2gFIS3MFqlH6Z3glRf6qlnhkeB6Pl/d2HgxV/2zZXkVNbdOQ5+XlZrH0qJksXTiTJYtmsuiIGaSl\npaKqNLd0UF3bSE1dE4dqGqmubaK6tjHw5b8/DPE07XZ3s29/LfvCqqL6KsjP7lXCKC32B5CSIv/r\nKfnZCRt0zfDFLCCo6nYR+T7wLNAGbAAiPfr0/SuK+Ne7Zs2a0HZFRcWkWud0sunq6g7L2P3fg5l9\nMAAEM/vhXc9N1f7aXnXvfeXlZlFa0hMwyorzKQkLIoUFOTgcI+uD0dTczuZAxr/53b1se//A0E/0\nwJxZpRx91CyWBAJA+bTCiJmtiJCXm0lebiYL5k2LeC2Px+u/f3VNPcGippFDtU2BINJAe3vXkGlq\naGylobGV7e/vj/h+qssZum/hbRihIFKUR3q6a8jPMWNTWVnZq5fTSMVt6goR+S5Qpao/D9v3c6BS\nVR8NvN4OnNa3ysjaECaOYGYfyugDmX1NfXPoyb6puS0qn+VwOkh1OujsHDoTHs61iqfk9goS/pJG\nT118Q2MrmwOZ/+ZtVYMGoKD0dBdLFpYHMv9ZLFlYTk52xpjTOxKtbZ1U1zZGLGFU1zRSe7h5yGqs\n4QgPusESRk/gyGNKfjYpKUMHXWtD6C9abQgxDQgiUqKqNSIyE3gGOEFVm8PeD29UXgn82BqVo6dv\nlUJ1bZ9MuM6fEXu6PeOd1BFzOB2UFOZSXOTPjHu+cv3fi/MoCDSWtrR2BJ6Gg0/KgfsRyABrDzfj\n9Xjjku6ppQUsWTSTJQtncvRRs5g3u2zEJY9483p91B1u7lfCCA8cra1jaPgOcDgdHLtsLv91x9WD\nHmcBob9kaFQG+L2IFOLvZfQlVW0WkesBVPVeVX1KRFaJyAf4q5UG/0swIapKa1sntXXNVIdlcsGM\nvjrwtD2cKopEE/40XlKUR3FhYLswN1RnXZCXNaynSYDcnExyczKZP2dqxPe9Xh/1DS2hIFFT1xTK\n8GoC3xubRl5qcTgdLJw/nSWLZnL0opksXjiT4sLcEV9nvDkcKaGnepgV8Zi29k6qa5pCf4vhJYzq\n2kZq6ocOul6PF+lXg2ziyWY7TVDtHV09VSmBf65gRh/M9Ds6hq77TTQOp4OighxKS/L9GX1RHsVF\nuZQGn/CL8yjIy064p+bwdo1QJhcIHMGgkZHuYvHCco5eNIsli2aycP50650TEAy6vQJFKHj4Xze3\ntHP+Ocfz9dWfHPRaVkLoLymqjKJlogUEVeVwY6u/6B38qg1uN0StCA6QkZEWymiDmW5pcR4lxf5B\nSyXFeWRm9J9O15h4a+/oorvbS15u5qDHWUDozwJCAvN6fWHdAHsy+kPButeaRtzu7jF/jsuV2iuj\nD1axBPuOlxTnkZWZZt0BzYQyWQLC7NmzeeCBB/jYxz425LHJ0oYwIXV2uqmubeLDQCbfk+n7G9xq\n6ptR39h6ZThT/d34SgNVKiXhT/eBTD83J3EHXBkz2WRn9x6L0dHRwZe+9CXuvvtu9uzZw9y5c8nK\nygq9f/PNN/PNb34T8C+hWV5ezre//e3Q++OxhKYFhGF6a9NOfrXub+yuqh5VA2NfWVnplJXkU1ZS\n4P9enB94nT/iRlNjzPgLn5q6ra2NsrIyLr300l7HNDc3J/RDnAWEYXj6rxv47k/Wj6hrYuGUHH9f\n60CmX1qcF8rwy0oKyM5Kj2GKjTHj6fe//z2lpaWcfPLJvfb7fD4cDkevfffddx+PPPIIIsKPf/zj\n0BKaABs2bODLX/7ykEtoRosFhCGs+8PL3HP/U732BfvA98roS/1P+jb3izHj5/zbno3q9f54+9mj\nOm/t2rV89rOf7bd/1qxZiAhnnXUWd911F4WFhQm1hKYFhAH4fD5+9uCzPPL4i6F9c2eX8Z2bP8WM\naYUJ1y3SGJMY9u7dy4svvsivfvWr0L7i4mLefPNNli9fTl1dHTfccANXXHEFTz/9dOiYwZbQBOKy\nhKYFhAg8Hi/f+8l6nv5rz9IMyxbP5j9uvZLcnMG7xBljJreHH36YU045hVmzegbxZWVlccwxxwBQ\nUlLCPffcw9SpU2lra+vV0NxX3yU0Dx48GLuEYwGhn45ON9/63jpee3NHaN/JKxdxx9cvt2ogYxLc\naKt4oumhhx7illtuGdaxvkBvxERZQtMCQpim5na+umZtr3njz/v4R/jql863KiJjzJD+/ve/c/Dg\nQS655JJe+9944w3y8vI44ogjaGho4MYbb+T0008nJycHSJwlNC2XCzhU08j/+dq9vYLBVZefztdX\nf9KCgTFmWB566CEuuuiiftVAu3bt4pxzziE3N5elS5eSkZHBunXrQu9fe+21vPvuuxQUFHDhhRdG\nvHZSL6EZTbEeqbxrTzVfue1Bauuagh/Il6//BBefe2LMPtMYMzqTZaTySNhI5SjZtHUPX7/j4dDc\nQc5UJ7f96yWcccrScU6ZMcbE16QOCC+/vo1b/+PR0LxCmZlpfO9bV3LcsnnjnDJjjIm/SRsQ/vzs\nm3z/v58IzTk0pSCHH97+uQGXIjTGmIlu0gUEVeXhx17g3rU9IxqnTS3kv759FTOmFo5fwowxZpzF\ntPuMiPybiGwVkc0i8oiIpPV5v0JEmkRkQ+DrW7FMj8/n48f3PdkrGCyYP42f33WdBQNjzKQXsxKC\niMwGvgAsUtUuEfktcDmwts+hL6jqebFKR5Db7eE7//V7nn/xndC+Y5fN43vfuoKsTJtozhhjYlll\n1Ix/LeVMEfECmcCBCMfFfPhdW3sn3/zuI/xjwwehfWecspRbv3IJLtekqzUzxpiIYlZlpKqHgR8C\nVcBBoFFVn+t7GHCSiGwSkadE5Khop+NwYyur/+3+XsHgonNP5PavX2bBwBhjwsSyymge8H+B2UAT\n8JiIXKGqvwk77G2gXFXbReQc4AlgQaTrrVmzJrRdUVFBRUXFkGk4cOgwX771Vxw4WB/ad91nz+Kz\nl1Yk9CIVxhgzkiU0gyorK6msrBz1Z8ZspLKIXAacpaqfD7z+DLBSVW8Y5JzdwLGB0kX4/hGPVH5v\n50H+9ba1HG5o8V8jJYWvrz6f8/7pIyP8SYwZnKpS09hJS8fY18lONCkC5cXZpDoTZ/qWRB6pXFFR\nweuvv47T6X/WnjFjBtu2bQPg+eef54YbbmDfvn2ccMIJPPjgg8ycOROIvITmnDlzeOCBBzjjjDOG\n/NxkGKm8HbhVRDKATuBM4I3wA0SkFKhRVRWR4/EHqMP9LzUyb23ayc3f+TXt7V2AfzH6O75xGaes\njHqNlJmEPF4fuz5sYfu+RrZVNbJtXyMNLV3jnayYKchO418vXsrSOVPGOykJT0T46U9/yjXXXNNr\nf11dHRdddBEPPPAA5557Lt/61re47LLLePXVV8cppZHFLCCo6iYReQh4E/Dhrx76hYhcH3j/XuBi\n4Isi4gHa8fdCGpO/vbyFNf/5OzzdHgCyszP4/q1XsnzJnLFe2kxSrR3d7NjfxPaqRt6tauD9A810\ndQ9/OdVk19DaxZqH3+K6VYv4p+NmjHdyEl6kJ/X169ezZMkSLrroIsBfBV5UVMSOHTt44YUXEmYJ\nzQk1ud36J1/nhz/7EwSOLSrM5Ud3XMW82WVDnGmMn6pS3dDBtn2NbA88/VfVtDLUn19WupOygom3\neFJ1YwetYVVhq44v59qPH4lzHGcAHqzK6J8ffSuqn/Xk5ceO6PjTTz+drVu3oqoceeSR3HnnnZx2\n2mncdNNNeDwefvrTn4aOPfroo7n99tu54IILuPrqq/stoTl79mzKysp44oknQkto3nTTTRGX0EyG\nKqO4UVUe+M3z/GrdX0P7Zs4o5kd3XMXU0oJxTJlJdF6fv/pnW1Uj2/c18m7V8Kp/SvMzWDgzn0Uz\n81lYns/MkiwcKYlTzx4tdU2dfHfdRnZ+2AzAU2/so6qmlW9cuozcLNc4py7xfP/732fx4sW4XC7W\nrVsXWvayra2N4uLiXsfm5ubS0tISem1LaEbJvoP1/Obxl0KvjzqynLtu+yz5eQMvTWcmp7ZOf/XP\ntr3+p//39jcNWf3jSBHmlOWwsDwQAGbmU5Q7OQYzFuWl871rPsLdT2zl5a2HANiyp4Gv/eJ1bvnU\nCmaVZo9zChPL8ccfH9r+7Gc/y7p163jqqafIzs6mubm517FNTU2hBXIGYktojsLM6UXc8Y3L+Lc7\nH+H4FfO585ZPk5FuTy+TnapS29TJu3sbQw3Ae2tahqz+yUxzcuSMPH8JoDyfBTPyyEibEP8qo5Lm\ncvDVS5YyuyybXz/vH89zqKGDb9z/Bl++aAknLCwZ5xT2GGkVT7wsXryYtWt7Jmloa2tj586dLF68\nGLAlNKPulJVH8ZPvXM2yxbNxOh3jnRwTBx1dHuqaO6lv7qKuqZO65k7qmjqpb/G/rm/upK3TM+R1\nivPSWTgzn6NmFrBo5sSt/hkLEeGSU+cysySb/3p8Cx1uDx1uD997dCNXnDGfi0+ZM+nH9jQ1NfHa\na69x2mmn4XQ6+e1vf8tLL73Ef//3f5Ofn8/XvvY11q9fz6pVq7j99ttZvnw5Cxb4h10lyhKaEyYg\ngH9uIjMxdLo91DV1Ud/cSW0gs69v7gp892f8w8ns+0qRQPVP4Ol/0cx8ivImR/VPNJywsIT/+PxH\n+N66jRxq6EAVfv38B+ytbuVfzl9MmmvyPox1d3dz6623sn37dhwOB4sWLeKPf/wj8+fPB+Dxxx9n\n9erVXHnllaxcuZJHH300dO61117LJZdcQkFBAaeffjrr16/vd31bQjMg1ktomvjqdHsGfKr3v+6i\nrTM6g7wy05wsmJHHwvJ8jpqZzxEz8sicxNU/0dLc5ub7v9vElj0NoX3zpuZyy6eWxzzAJvLAtPES\nrV5GFhAmqbbObv7y1gH+uukgze3xG2Hr7vZFLbNPdaZQmJtOUW5a4Hs6hblpFOel+1/npZObmTrp\nqzJixeP18cDTO3jqjX2hfQXZadx82TIWzsyP2edaQOjPAoIZldrGDv78WhV/efsA7V0jr3KJl1RH\nCoV5ltkng6f/sY9f/O92PF7//2iqI4UvnruIj62YHpPPs4DQnwUEMyI7Dzbzh7/v4e9bq/H6xvde\nBjP7wpw0ivL8GX1RIJMPZvx5WS7L7JPI5t2H+cHv3qG53R3ad/6Js/jc2UdEvYHeAkJ/FhDMkHw+\n5e0P6vjDK3t61fUGzSjK4rwTZ3HsEUXEK+91OlLIyUglJcUy+4nm0OF2vvvoRvZWt4b2rZhfyFcv\nPprsjNSofY4FhP4sIJgBdXV7eeGdD/nj3/eyv66t3/tLZhfwyZNmc+wRRZYxm6hq7/Lw4/VbeH17\nTWjf9KIsvvmp5Uwvis5AUQsI/VlAMP00t7n533/s46k39tHY5u71niNFOGlxKZ88aTbzp+WOUwrN\nZODzKY9W7uS3L/T0q89Kd/KvFx/NsUcUjfn6FhD6s4BgQg7Wt/GnV6v468aD/aZhyExzctYx0zl3\n5UyK8zPGKYVmMnp5yyHufmJr6G8yRYSrzj6C806cNab2IQsI/VlAmORUlW1VjTzx9728saOm33QM\nRXnpfOKEmZx97HSy0qNXf2vMSOz6sJnvrttIbVNnaN8Zy6fxfz6xiLTU0Q1is4DQnwWEScrr8/Ha\nthqe+Pte3tvf1O/9eVNzOf+kWXx0cem4TlFsTFBDaxff/+0mtlU1hvYdOSOPmy9fzpSckc/tP1kC\nwkiW0IxWQLAcI0l0dHn482tVfPEnr/CD373TLxgce0QR3/7csfzw+hM47eipFgxMwijITuOOzx3L\nmWHjEnbsb+Kr973O+wf6P9QkK7fbzbXXXsvs2bPJzc1lxYoVPP300wDs2bOHlJQUcnJyQl933nln\n6NyrrrqKW2+9tdf14jFVRV+DjuEXkRLgEuBUYDagwF7gReAxVa0Z+GwQkX8DrsS/Ytpm4GpV7epz\nzN3AOfhXTLtKVTeM6ieZoOqbO3ny9X08/eb+fiN8Ux0pnLZsKuefOIuZJTYNsUlcLqeD1ecfxazS\nbH71zHv4VKlv7uSWX/6DGz+5hFOWJv8iVh6Ph5kzZ/Liiy8yc+ZMnnzySS699FK2bNkSOqa5uTmh\nx9cMWGUkIg8A84D/xb8W8oeAAFOB44GPAx+o6ucHOH828Fdgkap2ichvgadUdW3YMauA1aq6SkRO\nAH6iqisjXGvSVRntPtTCn17dy4ubPwyNAA3KyUzlnOPKWXVCOQXZsVtOz5hY2PBBHXc9trnXA84l\np87h06fPH1Y36GSqMlq2bBlr1qxhxYoVzJ07l+7ubhyO3m0n9913H6tXr0ZEcLlcoSU058yZw+rV\nq3nooYeGXEIz5m0IInK0qr4z6MmDHCMiU4BXgZVAC/AH/Bn+c2HH/Bz4m6r+NvB6O3Caqlb3uda4\nB4T9tW28saOGrm5fzD9rx75GNuys77d/6pRMzjtxFmcsn0q6yyZoM8nrYH0bdz6ysdc4mRMWFvN/\nL1w65OSDgwWEe7e/HtV0Xr/whFGfW11dzezZs9m0aRMul4u5c+cybdo0RISzzjqLu+66i8LCQoDE\nX0IzUkYfyORnBN8bLGCo6mER+SFQBXQAz4QHg4DpwL6w1/uBGUA1CaCjy8MrW6t5bsOBXg1i8XbU\nzHzOO2kWxx9ZbPP0mwlhWmEWP/jC8fzw95t56/06AF7fXss37n+DOz53bNKXfLu7u7niiiu46qqr\nWLBgAW1tbbz55pssX76curo6brjhBq644opQGwMkyRKaIvICcG7g2LeAWhF5RVW/PMR584D/i7/t\noQl4TESuUNXf9D20z+uIoX/NmjWh7YqKCioqKoZK+qioKu8daOK5tw/y0uZDdLjHZwK4FBFWLirh\nkyfN4sjy2M0cacx4yUpP5ZufXs6vn/+A9S/vAaAwJ43czOTuJu3z+fjMZz5Deno699xzDwBZWVkc\nc8wxAJSUlHDPPfcwdepU2trayMoaeAT3SJfQrKyspLKyctRpH069Q56qNovI54GHVPU2Edk8jPOO\nA/6uqvUAIrIeOAkIDwgHgPKw1zMC+/oJDwix0NzmpvKdD3luw4Fec7EEOR3CcQuKmV06+Bqo0ZDu\ncnDiohLKpmTG/LOMGU+OlBQ+d9YCZpZk88Qre/nXS44eUyl4LFU80aCqXHvttdTW1vLUU0/1ay/o\ny+fzV0FHawnNvg/Lt99++5DnhBtOQHCIyFTgUuBbgX3DqdDfDtwqIhlAJ3Am/sbpcH8CVgOPishK\noLFv+0Es+XzKpl31PPf2QV7fXkO3t3/7wIyiLM48ZjoVy6YmfTHWmER1+rJpnLq0LOmrRL/4xS+y\nfft2nnvuuV6Nv2+88QZ5eXkcccQRNDQ0cOONN3L66aeTk+N/wEymJTTvAJ4BXlHVNwJVQe8PdZKq\nbhKRh4A38Xc7fRv4hYhcH3j/XlV9SkRWicgHQBtw9Wh/kJGobezguQ0H+evGg9Q0dvR7Py3VwSlL\nyjjzmGksLM9P6G5ixkwUyR4M9u7dy3333Ud6enqvqp57772XlJQUbrnlFmpqasjNzeXss89m3bp1\noWNsCc0RiEYvo26Pjzd21PKXt/ezcWd9v6keABbMyOPMFdM5ZWmZLbNoTIJKpm6n8RLzXkZhF5wH\n/Bg4EX9V0d+BL6vq4OWbBFFV08pzbx/gb5s+7LV4R1BOZioVR0/lzGOmx6V9wBhjEtVwHoMfAe4B\nLgy8vgxYB4xv680gOro8vLy1mufePsD2fZG7iy6fV8iZK6ZzwqJiXM7RTbJljDETyZBVRiLyjqoe\n3WffJlVdFtOU9f68IauMVJX39jfx7NsHeGVLdcTuokV56XxsxTQ+tnw6pQU2FbQxyciqjPqLeZVR\nYBCaAP8bmJMo2AJyGf7pLBJGa0c3//bAP6iqjdxd9PgjSzjzmOksnzcl6RuujDEmVgarMnqb3t1L\nrwt8l8D+m2OVqJHKzkjFldo7oy8vzuKsY6ZTsWwaeVmucUqZMcYkj8Gmrpgdx3SM2ZnHTGd/XRsn\nLynjrBXTObI8z7qLGmPMCAw2ud1FDDIATVX7d5SNkeG0IXS5vXhVrbuoMROctSH0F49up+cy+Ijk\nuAWE4UhzWU8hY4wZi0kzMM0YMzFMhhLCSJbPhDgsoSkiV4nIYL2QXCISl6kmjDEm0d1zzz0cd9xx\npKenc/XVvbPG559/noULF5KVlcUZZ5xBVVVV6L1EWT4TBl9TORv4h4isE5GviMinReQKEflXEVkH\nvA5YZ35jjAGmT5/OrbfeyjXXXNNrf11dHRdddBF33nknDQ0NHHfccVx22WXjlMrBDRgQVPUe4Bjg\np4ALOBn4KP52h3uAY1T1/8UjkcYYk+guuOACzj///NAqaEHr169nyZIlXHTRRbhcLtasWcOmTZvY\nsWMH9913H4888gg/+MEPyMnJ4fzzzw+dt2HDBpYtW0Z+fj6XX345XV1dfT8y6gbtkhOouH858GXM\npNbe7aW23U1tu5u69m5q29zUtLupC+zr8sR+edUgR4pQnOnyf2W5KMpIpSTLFdqX7XJMym7Xmw9H\nN6taOuXkEZ/Tty5/69atLFvWM7FDZmYm8+fP59133+W6667j1Vdf7bd8pqry2GOP8cwzz4SWz3zw\nwQcjLp8ZTdZH0xjA41PqwjL32vbusMzfTW2bm9Zu73gns5faCJM1BmU4HaFgUZKZ2hM4Ml2UZLoo\nzEgl1WGj9mOhbyBua2ujuLi4177c3FxaWlpCrxNh+UywgGAmAVWlqcvT82Tf7qamrSfzr2vvpr6j\nGx3Wuk9eZuw6AAAgAElEQVTJocPjpaq5g6rm/ut9AAjClAwnRZk9pYrwEkZxlovcSVrKGKu+mXt2\ndjbNzc299jU1NYUWxxnISJfPjIbhTH89t+9U15H2DXDukcCjYbvmAreq6t1hx1QAfwSC13tcVb8z\njLQbg6rS4vZS3+F/qq8Py/TrOnq23RFWwxup1BQJZaDBJ23/U7f/CTwr1dFvgfBY6fL6ekoxbe5Q\naaYmsN01xM+rKPUd/kC4o74t4jFpjhSKM10sK83hU4unUpCR+Gsdj6aKJ9r6BtHFixezdu3a0Ou2\ntjZ27tzJ4sWLIx4/nGvGynBKCI8DK/rseww4dqgTVXVH8FwRScG/XvIfIhz6gqqeN4y0mElkqMy+\nPrA9VOY3HIJQkO7sVa0SzOiDX3npTlIS6Il5Wk56xP3B+xas7qoJCxi1gbaP4ZSIurw+9rd0sr+l\nk8q9h/n0kql84ogSnCmJcw8Sidfrpbu7G4/Hg9frpaurC6fTyQUXXMDXvvY11q9fz6pVq7j99ttZ\nvnw5CxYsABJn+UwYfLbTRcBRQJ6IXEjPpHa5QOS/xMGdCexU1X2RPm4U1zNJTFVpdXtDT/F1YZl9\nfUd3T0NtFDJ7gKxUR+jJ3p/Bp07YOnURITfNSW6ak3kFmRGPCbaZhDeSBwNHMIi0e3raTNq6vfxi\nw36e3lnH9ceUs6IsN14/TtL49re/3ath+Ne//jVr1qzh3//933n88cdZvXo1V155JStXruTRR3sq\nThJl+UwYfC6j84EL8E9h8aewt1qAR1X17yP6IJFfAm/27aoqIqfhnwZjP/4SxFdV9d0+x9hI5STT\nFnxCDau2CWb60c7s050pFGUEq3JSKQpk/EWZqRRl+OvGM1NtapORanN72VLbwi83HmB/S2ev906c\nns/nV8ygLDttgLNjZzKMVB6paI1UHs4COSeNNPOPcA0X/sz+KFWt7fNeDuBV1XYROQf4iaou6HOM\nBYQE0unxUhuWsfc84btD+8OfLseib2ZfmBmW8Wf4v2elWuNnLHV7ffz5/VrWbfmw1+/VlZLChQtL\nueSoUtLjuOqgBYT+4hkQSoAvALPpqWJSVb1mwJP6X+N84Iuq+vFhHLsbOFZVD4ft09tuuy10TEVF\nBRUVFcP9eDMCXR5fqM6+Vybf0fOE3xJhNbrRsMw+uTR0dPPgOwd4bnd9r/3FmS6uWT6dU8oL4lOt\nYQGhn+A9qayspLKyMrT/9ttvj3pAeBV4EXgLCJbxVVUfH0FiHwX+V1XXRnivFKhRVRWR44Hf9V2L\nwUoI0RHe176uvZvasIbZ4NN+U1d0MntXoIdKMHMPr84Jbltmn5x21Lfx87f28d7h3r2TlpbkcP0x\nM5iTH7ndIlosIPQXzxLCRlVdPvIkhs7PAvYCc1S1JbDvegBVvVdEbgC+CHiAduArqvpan2tYQBiC\n16cc7oxQjRNsoG1z09DpiUpfe6dIqJ6+OKzOPnzb+rBPbD5Vnttdz4ObDvR6iEgRYdX8Iq5YMo3c\nGK1NYgGhv3gGhO8Ar6rqkyNOZZRM9oDgU6Wp0xOaJiG8kba+w799uKMbbxTukUOEKRmp/Z/oM3q2\nE637pRk/rW4Pj249xJ/fq8ET9veXm+bkM0un8U9zi3BEuZuqBYT+Yh4QRKSVngVysgA30B14raoa\nt35nEz0gtHd7qQnMixPe9a+23V+lU9fu7vXPNlrho1OLMl0UZ/R/ui9IT436P7CZ+PY2dfCLt/ez\nobr3iNx5BZlcf0w5i4uzo/ZZFhD6i1sJIREkc0AITptQHRwc1Cfjr2lzR62RNj8tNTSYKtjtMny7\nMMNlg4pMzKgqrx1o4hcb9lPd1ntmzopZU7h62XSKMl1j/hyriowsXlVGx0TY3QTsVdXo5GRDSOSA\n4PUpdR09GX14hj/caQSGI8fljNgw63/SdzElI5U058QYWGWSW5fHx/rt1Ty27VCvv/10ZwqXHTWV\n8xeU2N9qnMQiILyGf5qKdwK7lgJbgTz8XUmfGWVahy0RAkJHt5dXDzRyoLmLmvYuqgMZfjTq7l0p\nKYFZKf1z44RG0oY93cezn7cx0VDT5uaXG/fz0r6GXvunZqfxhRUzOH5anj3tx1gsAsJ6/BPSbQ28\nPgr4NvB1YL2qLhvs/GgY74Dw5odN3POPqkGnGx5McNqEkqyeGSVLslyUZqVRYo20ZoJ7p6aFe9/a\nx56m3jOvHjs1j+tWzGBG7mhmwjHDEYuAsFVVF0faN9YuqcM1XgGhucvDLzbs56976gc9riA9tVdG\nXxKW+ZdkppHlsqd7M7l5fcpTH9Tymy0f9mozc4pw3pElfGrxVJteJAZiERB+B9Tjn8ZagEuBYuBK\n4GVV/cjokzvMRMY5IKgqL1U18PO39/XqY52b5uTjc4soy07rNXe81YcaMzxNXR4efucgT++s6zUm\nZkp6Klctm87ps6dYaTmKYhEQMoEv4V9PGeAV4P8BnUBWcLBZLMUzINS1u/nZW/t47UBjr/0Vs6bw\nhRUzyE9P/DnhjUl0OxvaufetfWyta+21f2FhFredOj9mg9omG+t2Oko+VZ7ZWcevNh2gLWypxKIM\nF186rpwTpufH9PONmWxUlRf2NvDLTfup7/APcVpclM33P7bAGpujJJoD0x5T1UtEZHOEt1VVjx5t\nIkcq1gHhYEsnd/+jis01vQs758wr5upl060NwJgY6uj28rt3D/HEezX855lHDriGgxm5aAaEaap6\nUERmR3pfVfeMJoGjEauA4PUpT+yo4ddbDvZaYnF6TjqrPzKTo0sGX/PUGBM9zV0eqyqKsphUGQWC\nwnxVfS7QpuCIR9tB2OdHPSDsbmznJ29U8X7YjI0OES5YWMqnF0+1hmJjTNIbaUAYMhyLyHX410OY\nAswDZgA/Az422kSOJ7fXx2+3HuL32w71mh9oXkEmN35kFvOnWHHVGDM5Dad8dgNwPPAagKq+F1g0\nJ+m8W9vK3f/Yy77mnuUAXSkpfGrJVC5cWGrz/BhjJrXhBIQuVe0KtvqLiBOiMKl+HHV0e1n7zkH+\n5/3aXn2fFxdl8y/Hz6LcRkoaY8ywAsILIvJNIFNEzsI/JuHPsU1W9LwVmHaiJmzaiQyng6uWTWfV\n/CIbBGOMMQHDGZiWAnweODuw6xng/qFaeUXkSPyjm4Pm4p8T6e4+x90NnIN/tbSrVHVDhGuNuFG5\nucvD/Rv283yfaSeOnZrH6uNmUpI19ml4jTEmkcVipPLHgL+rasegBw5+jRTgAHC8qu4L278KWK2q\nq0TkBOAnqroywvnDDgiqysv7Gvn5W/to7OoO7c9xObnumBmcPmuKDXoxxkwKUe9lBHwO+JmINAAv\nBr5eVtWGwU/r5UxgZ3gwCDgPWAugqq+LSL6IlKpq9QiuHVIfmHbi1T7TTpw2cwrXHWPTThhjzGCG\nDAiq+lnwD1QDLgZ+CkwbzrlhLgceibB/OhAeJPbj79Y6ooCgqjy7q54HNu7vNe1EYUYqNxw306ad\nMMaYYRjOOITPACcDRwO1wD3Ay8P9ABFxAecC3xjokD6vR9yDqbrNzc/f2ofb1zPa2KadMMaYkRnO\nU/6PgZ34B6NVquruEX7GOcBbqlob4b0DQHnY6xmBff2sWbMmtF1RUUFFRUXodVl2Gp9eMpUH3znA\n1Ow0bjx+lk07YYyZdCorK6msrBz1+cNpVBZgMXBK4Gs+8J6qXjmsDxB5FPhfVV0b4b3wRuWVwI9H\n26js8Sl/eq+GVfOLSbdpJ8wE41Uf9Z3tHOpopbqjBbfPy/TMXObnFpKdmjbeyTMJKha9jHLxVxmd\nGvgqAl4Lti0McW4WsBeYE5z7SESuB1DVewOv7wE+DrQBV6vq2xGuM+5rKhsTTx2ebqoDmX91Ryu1\nnW141Bfx2KkZ/sAwN3cK6Q6bHM70iEVAeAf/ojgvAS+q6v6xJXHkLCCYicynSqO7g+qOVg61+wNA\nU3fn0Cf2kYIwIyuP+bmFzMopwJVi7WeTnS2QY0yCc3u91HS2BkoAPVVAQ8lNTac0I5vSjGwcksLO\nlnoOtDX3mo4lyCkpzM4uYF5uIeXZeTjEqlEnIwsIxiQQVaWl2x2q+qnuaKW+qz1iJh7OISkUpWdS\nmpFDWSAIZDr7j65v97jZ1dLAB811VHe0RrgSpKU4mZs7hfm5hZRl5Nh0LZOIBQRjxpFXfdR1tocC\nwKGOVto97iHPy3CkUpaREyoBFKVn4UwZ2VN9s7uTD5rr2dl8mMPu9ojHZDldzMstZH5uIUVpmTZq\nf4KzgGAmJZ8qnd5uWrvdtHnctHW7afW4cXuHroqJBkVpdHdS29mGd4DG3yBBmJKWEcj8/UEgNzUt\nqplzfWd7IDjU0+LpinhMviuDeblTmJ9bRL4rfjP+qipun5fWbjetni7/927/99LMHJYUlMYtLRNd\nLBqVS/AvkDObnnELqqrXjDaRI2UBYXKLlNm3edyh162Bfb4EnZXdleKgJD2b0kx/AChJzyItTr2B\nfKrUdLT6g0PLYTq93RGPK07PYn5uEfNyppCVOraJH73qo627u1dmH/x9BYPAQG0mc3OmcNb0I8b0\n+aZHLALCq/jnL3oLCD76qKo+PupUjpAFhIkrmNkHM5A2Tzdt3V2hTL7N46bN0z3kU3ciCW/8LcvI\noSAtIyHq7b3q40BbMzub69nd2kB3hExZEKZm5jA/t5A5Of27saoqXT5v6Ine/3vyb7cEAkC7p3vI\nNpKBlKRnc8HsxaM61/QXi4CwUVWXjzllY2ABITaa3J3sbW3E7fPE7TO9Pl8gE+mmLRAAopXZpzuc\nZDldZDldZKf6v6c74jehYYbTOWDjb6Lp9nmpam3kg+Z69rU1RfwdpCDMzM4nw5naU63jcUcMJCPl\nlBSyU9PIDvyugtt5rnTKMm2WgWiJRUD4DvCqqj451sSNlgWE6Gl2d7Kr5TA7mw9T19U23skZtrQU\nZyiTD37PCr4ObKdav/tR6fR62NNymA+aD3OwPXI31pHKDPxewjP78O10h9MatOMgagFBRFrpmWgu\nC3ADwQpIVdXcsSR0JCwgjE1rdxc7mw+zq+UwNZ2RuyaOp7QUJ1mpPRl7eCZvmX18tXW72dlymJ3N\n9QP+rTglhZzUNH8G70zryegDv69Mp2vEPaRMbFgvIwP4/7F3tRxmZ0v9gP3THZJCeVYehWmZ/eec\njZEUxP/0GPaUbyNqE1Oju5N9rY2ISK8n/LQUhz3dJ4lYVBk9r6ofG2pfLFlAGJ7gIKVdzfUc6miN\nWPQPTm8wL7eQWdn5cevtYoyJv6itmCYiGUAmUCwiU8LeysW/sI1JAB2ebna3HGZny2E+bG8ZMAhM\ny8plXk4hs3MKbAI0Y0xEg+UM1wM34V8d7a2w/S34F8kx46TT62F3y2F2BRoBI/W/F4RpmTnMzSlk\nTk4BGU5bPtQYM7hBq4xExAHcoqrfjl+SIqZj0lcZdXk97GlpYFfLYfa3NQ0YBMoycpiXO4U5OQVJ\n0f3RGBM7Ng5hAnF7vexpbWBXSz37WiMHAYCyjBzm5kxhbhRGmRpjJo6otSGEeU5ELgYen5S5chz5\nVGn3uDnU0crOQQYMgX9E57xcfxCwFbOMMdEwnBJCK/7GZS8QXLVjWOMQRCQfuB//EpwKXKOqr4W9\nXwH8EdgV2PW4qn4nwnUmRCzy+Hyh0Z4t3T3zvLSEzfcy2Hw8xelZgZJAIbkuCwLGmMFFvYSgqtlj\nSM9PgKdU9WIRceIf4NbXC6p63hg+I2F0eT09Gb0nmNn3ZPztA0wsNpjCtEzm5kxhXm4heXGckdIY\nM/kM1u10kapuE5FjIr0fae3jPufnAaeo6ucCx3uApkiHjiC946rT66HJ3dk7o/f0POEPZ9WroWQ4\nUslNTaM8O595OVPIT8uIQsqNMWZog5UQvoJ/2usfQcR6jNOHuPYcoFZEfgUsw9919SZVDV+5Q4GT\nRGQTcAD4qqq+O9zEx0un18Or1Xt5v7l+TPO8BEfp5gRGfIaG/6emkeP0f7ch/8aY8TJgQFDVLwS+\nV4zh2scAq1X1HyLyY+Bm4N/DjnkbKFfVdhE5B3gCWBDpYmvWrAltV1RUUFEx2mQNn6qyq+Uwr1Tv\npWMY1T2hGRxTXeT0meMlJzWNLKcrIaZBNsZMTJWVlVRWVo76/OE0KmcAXwJOxv9E/xLwM1XtHOK8\nMvyzpM4JvD4ZuFlVPzHIObuBY1X1cJ/9cW9Ubut283L1Hva0NvTaPyUtk9ywib3Cn/IzbAZHY0wC\niUW304eAZuBu/PX9nwYeBi4Z7CRVPSQi+0Rkgaq+B5wJbO2T2FKgRlVVRI7HH6AOR7pevKgq25tq\nea2mqlebQJbTxcmls5mdUzCOqTPGmNgZTkBYrKpHhb3+q4gMt57/X4DfiIgL2AlcIyLXA6jqvcDF\nwBdFxAO0A5cPP+nR1+Tu5KVDuznQ3txr/1H5JRxfXG4TwRljJrThVBn9Gvipqr4aeL0SuEFVPxOH\n9AXTENMqI58qmw8f4s26/XjCBoLlpaZz6tQ5TMuM29IPxhgTNdGc7XRz2DGviMg+/G0IM4EdY0pl\nAqnvbOfFQ7t7LQYiCMumlHFs0Qzr9WOMmTQGqwM5d5D3kn7YsMfnY2P9QTbUH+w1OrgwLZNTy+ZQ\nkjGW8XjGGJN8JuWKaYfaW3jx0G4a3B2hfQ5J4ZjCaSwrnIpDrFRgjEl+sehlNGG4fV7erN3Plobq\nXgPMyjJyOLVsDgU2KtgYM4lNmoCwr62Jlz7cTYunK7QvNcXB8cXlHJVfYgPGjDGT3oQPCJ1eD6/V\nVLGjqbbX/vKsfE4pm02OTR1tjDHABA8Iu1oO88qhPb1mGU13ODmxZBZH5BbaqOIJSFVx+zpo97TQ\n7mnBqx5SU1ykpqT1+nJKqv3+jeljQgaENo+bVw7tZXdr70HP83IKOal0pi0tOYF41UtHIPNv9zSH\ngsBQBOkXJFwpLlIdaaRKGqmONBwyIf89jBnQhOplpKrsaKrjtZoqunw9mYJNOzEx+J/+O0MZf7un\nhU5vW8w+L0UcuPoEjZ7gkYYzxUWK9UgzCSzqayonguEEhGZ3Jy8d2sP+9t5LLizKK+GEEpt2Ihn1\nfvr3lwCG8/TvECeZzhwynTmkpqTj8blx+7roDvsaznWGwxmqjnJFDB5WNWXG06QMCM3uLh7b/U6v\naSdyU9M5tWwO07Ns2olk0PP031P1M9yn/3RHJpnO3FAQcKVkDJkJe9XbK0AEv0KBw9s1prUvgvxV\nU+FtGP3bMxxis+Sa2JiUAQHgLwfeZ1fLYQTh6CllHFs0ndQUR5xSaEbKp17aPa2hzL/D04JHh15z\nIvzpP9OZS4YzOyZ1/aqKVz10+zoDQcLdL3B4fO6ofFaKpAwQLFzWnmHGZNIGhHaPm+cP7uSE4nKb\ndmIQPvX2ydzcUXkSHi6Pr5sOTzMd3naGMwOK/+k/h4xACSBtGE//8aLqC91Ld9j9jEXVVIo4Qm0X\nqSkunCku4rX6bLYzj6zUvLh8lomuSRsQTO8MKlLm1O1zD+spfLw4xBF46g+UABw5OFKS+8nYF6qa\n6t+OEfz9+HTsa3HHljA/92gynDnjnRAzQhYQJihVpVv7Z/Dh36NVhREvaYGn/8wEfPqPF1UNCxrB\n6qi+v+fotGeMRU5qAbNzFo9rGszI2VxGSainvroLt68zwtO9m26fm+hMMhts5Oypr06R+LW1CEKG\nM3tCPP1Hg4jgECeOFCfpZEU8pnd7hjuqVVGDUfVR27kfgJbuBto9LWRaKWFCi+l/pIjkA/cDi/Hn\nZteo6mt9jrkbOAf/imlXqeqGWKZpPKgqHnXj9nb1yfSDdc+d+MJ6SI1F326Qzr4NleKadE/hyU5E\ncEoqzpRU4j39otvXSZO7DoCajiorJUxwsX5E+wnwlKpeLCJO6P0IJCKrgPmqeoSInAD8DFgZ4zRF\nnU99/TL4bl9XKABEq8jvlNTIPVGCfd5toJSJspKMmaGAYKWEiS9mAUFE8oBTVPVzAKrqAZr6HHYe\nsDbw/usiki8ipapaHat0jZaq0uZposvbEdYVMXrdD/1dD9MjDG7qeeKPZ9WOMeDv5ZXnKrJSwiQR\nyxLCHKBWRH4FLAPeAm5S1fawY6YD+8Je7wdmAAkVELp9XVS1bqfd0zLqazjE6c/sHb0zfVdKug1O\nMgnNSgmTRywDghM4Blitqv8QkR8DNwP/3ue4vrlgxLqVNWvWhLYrKiqoqKiIWkIH09rdxL7W7UN2\n1/Q/yfdk9v7MPy20z57uTbJKd2SS7yqm0e2fQt5KCYmrsrKSysrKUZ8fs26nIlIGvKqqcwKvTwZu\nVtVPhB3zc6BSVR8NvN4OnNa3ymg8up2qKnWdBzjUsZfwGJXvKsblyAhk/OmhgUJidfdmAuv0tvN+\n09uh1/NyjybTadPCJLqRdjuNWS6mqoeAfSKyILDrTGBrn8P+BHwWQERWAo2J0H7gVQ9Vbds51LGH\nYDBwSipzcpZQnn0kpRkzKUgrJTs1D5cj3YKBmfCCpYSgmo59gxxtklWsexn9C/AbEXEBO4FrROR6\nAFW9V1WfEpFVIvIB0AZcHeP0DKnT08betm24vZ2hfZnOHGZmLyQ1xVZXM5NXSUZ5qNqopbuBNk8z\nWVZKmFBspHKYxq4aDrR/0GtMQGH6VMoy5lh3TmOAfa07QkEhOzWfOTlLxjlFZjA2UnkUfOrjUPtu\n6rs+DO1LkRSmZx5BflrxIGcaM7mElxJauxutlDDBTPrHXrevi90tm3sFgzRHBvNyllkwMKaPNEcm\n+a6S0OuajqpxTI2JtkkdEFq7G9nZtLHX+IJcVyHzcpeR7ow8r4wxk11JxozQdmt3I23dzeOYGhNN\nkzIgqCq1HfvZ3bIlbHyBUJYxm5lZC20xEmMG0a+U0GmlhIli0gUEr89DVeu2QJdSP2eKv0tpccYM\nGy1szDCUZJSHtq2UMHFMqoDQ4Wnjg+aNNHcfDu3LdOYyP3c52bYilDHDlubIsFLCBDRpAkJDVw27\nWjbh9vWMLyhKn8bcnCU2vsCYUfCXEvwlaislTAwTPiD41MeBtp3sb3svNL4gRVIozzqSqZlzbZSx\nMaPkLyWEjV62UkLSm9C5odvXxa6WzRzu26U0d7l1KTUmCvqXEvrOcG+SyYQNCK3djXzQtIGOSF1K\nHZnjmDJjJo40RwYFaVZKmCgmXEBQVWo69rG7ZUvYurPC1Mw51qXUmBgoTg8vJTRZKSGJTaiAEOxS\nWt2xN7TPmeJibs4SitKnW5dSY2Kgbymh2koJSWvCBIQOT+uAXUqzrEupMTEVXkpos1JC0poQAaHL\n28GulncG6FLqGseUGTM5+EsJPeMSrJSQnCZEQHClpJPnKgIgRRzMzF5oXUqNibMSKyUkvQmRY4oI\n0zLnkecqYl7uslBwMMbEj8uR3ruUYDOhJp2YBwQR2SMi74jIBhF5I8L7FSLSFHh/g4h8azSfEywZ\nWJdSY8ZPSXo5EiwleJpotVJCUolHH0wFKlT18CDHvKCq58UhLcaYGHI50slPK6Ghy780ek1HFdmp\nS8c5VWa44lVlNFR/T+sPaswEYaWE5BWPgKDAcyLypoh8YYD3TxKRTSLylIgcFYc0GWNiJFhKCLJV\n1ZJHPKqMPqqqH4pIMfAXEdmuqi+Fvf82UK6q7SJyDvAEsKDvRdasWRParqiooKKiIrapNsaMWkl6\nOY1dNSgaKiXYFPOxV1lZSWVl5ajPF1WNXmqG+jCR24BWVf3hIMfsBo4Nb3MQEY1nOo0xY7e/7f1Q\nW0KWM5c5OUtttoA4ExFUddg3PaZVRiKSKSI5ge0s4Gxgc59jSiXwVyIix+MPUoM1QBtjkkDvtoRm\n2jzWlpDoYl1lVAr8IZDfO4HfqOqzInI9gKreC1wMfFFEPEA7cHmM02SMiYPguITDYT2Ospx5VkpI\nYHGtMhotqzIyJjm5vZ281/QWiv//d07OErJT88c5VZNHQlUZGWMmN38poTT0uqajCnu4S1wWEIwx\nMVWcPsPaEpKEBQRjTEz1LSVUWykhYVlAMMbEXHgpod1KCQnLAoIxJuaslJAcLCAYY+LCSgmJzwKC\nMSYurJSQ+CwgGGPipm8podXTOM4pMuEsIBhj4sbGJSQ2CwjGmLgqzigPKyW0WCkhgVhAMMbElSsl\njYK0stBrKyUkDgsIxpi4K86YYaWEBGQBwRgTd1ZKSEwWEIwx48JKCYnHAoIxZly4UtKYElZKsHEJ\n4y/WK6btEZF3RGSDiLwxwDF3i8j7IrJJRFbEMj3GmMQSXkro8LTQ6mkY5xRNbrEuIShQoaorVPX4\nvm+KyCpgvqoeAVwH/Gy0HzSWhaXHW7KmPVnTDcmb9mRNN0ROe2q/UsK+hCslJPM9H6l4VBkNtlrP\necBaAFV9HcgXkdJBjh9QMv/SkjXtyZpuSN60J2u6YeC0h5cSHOLAp944pmpoyXzPRyrWayor8JyI\neIF7VfUXfd6fDuwLe70fmAFUxzhdxpgEkZqSxtTMuaQ7sshKzR3v5ExqsQ4IH1XVD0WkGPiLiGxX\n1Zf6HNO3BJFY5UVjTMwVpk8d7yQYQOJVXycitwGtqvrDsH0/BypV9dHA6+3Aaapa3edcCxLGGDMK\nqjpYtX0vMSshiEgm4FDVFhHJAs4Gbu9z2J+A1cCjIrISaOwbDGBkP5AxxpjRiWWVUSnwBxEJfs5v\nVPVZEbkeQFXvVdWnRGSViHwAtAFXxzA9xhhjBhG3KiNjjDGJLSlHKovIL0WkWkQ2h+1bIyL7A4Pg\nNojIx8czjZGISLmI/E1EtorIFhG5MbB/ioj8RUTeE5FnRSR/vNPa1yBpT+j7LiLpIvK6iGwUkXdF\n5HuB/clwzwdKe0Lf8yARcQTS9+fA64S/5xAx3clyv/sNBB7pPU/KEoKInAK0Ag+p6tLAvtuAFlX9\n0bgmbhAiUgaUqepGEckG3gI+ib+qrE5VfyAi3wAKVPXm8UxrX4Ok/VIS/75nqmq7iDiBl4Gv4h8D\nk/6n+w4AAAdQSURBVND3HAZM+8dI8HsOICJfAY4FclT1PBH5Aclxz/umO+HzFgAR2Q0cq6qHw/aN\n6J4nZQkh0HU10hj3hG58VtVDqroxsN0KbMM/FiM0QC/w/ZPjk8KBDZJ2SPz73h7YdAEO/H87CX/P\nYcC0Q4LfcxGZAawC7qcnrQl/zwdIt5Dg9ztM33SO6J4nZUAYxL8E5kR6IFGLo0EiMhtYAbwOlIb1\nrqrG3yCfsMLS/lpgV0LfdxFJEZGN+O/t31R1K0lyzwdIOyT4PQf+C/ga4Avblwz3PFK6lcS/39Az\nEPhNEflCYN+I7vlECgg/A+YAy4EPgR8Ofvj4CVS5PA7cpKot4e+pvw4vYevxAmn/Pf60t5IE911V\nfaq6HP8o+FNF5PQ+7yfsPY+Q9goS/J6LyCeAGlXdwABP1ol4zwdJd0Lf7zAfVdUVwDnADYGq9ZDh\n3PMJExBUtUYD8Bf3+k2mlwhEJBV/MHhYVZ8I7K4O1NEjIlOBmvFK32DC0v7rYNqT5b4DqGoT8CT+\n+uGkuOdBYWk/Lgnu+UnAeYE67XXAGSLyMIl/zyOl+6EkuN8AqOqHge+1wB/wp3NE93zCBITADxt0\nAbB5oGPHi4gI8ADwrqr+OOytPwGfC2x/Dnii77njbaC0J/p9F5GiYBFfRDKAs4ANJMc9j5j24D94\nQMLdc1W9RVXLVXUOcDnwV1X9DAl+zwdI92cT/W8c/J0PRCQnsB0cCLyZEd7zWM9lFBMisg44DSgS\nkX3AbUCFiCzHXyTaDVw/jkkcyEeBK4F3RGRDYN+/Af8B/E5ErgX24O+5k2gipf0W4FMJft+nAmtF\nJAX/A9DDqvp84GdI9Hs+UNofSvB73lewmiIZ/s6DhJ50/0BElpHY93uggcBvMoJ7npTdTo0xxkTf\nhKkyMsYYMzYWEIwxxgAWEIwxxgRYQDDGGANYQDDGGBNgAcEYYwxgAcEkGRE5TURODHv9oIhcNMZr\n7h57ysZORJ4LDi6K8+feKCKfiffnmsRjAcEkm9PxTzEQNOqBNOKXEP8DInIGsKPv3FZR/oyBftZf\nAf8Sq881ySMh/hnM5CEis0Vku4j8SkR2iMhvRORsEXklsIjHRwLHTRGRJwIzTL4qIksDs6xeD3xZ\nRN4WkZMDlz01cP7O8NKCiHxNRN4IXGNN2OfvEJG1+If2zyAwv4uIZInIk+JfkGaziPQb1SkiXwhc\nc6OI/D4wpQQiMk9EXhP/AiXfEZGWwdIRwaeBPwaOv11Ebgo7/07pWZAo4rVE5A/in+Vyi/TMdImI\ntIrIf4p/xtQTReQ/xL/I0SYRuQsgEITqRWTxkL9AM7Gpqn3ZV9y+gNlAN7AY//QAbwIPBN47D/hD\nYPu/gVsD26cDGwLbtwFfCbveg8BvA9uLgPcD22cD9wa2U4A/A6cEPt8LHB8hbRcB94W9zo1wzJSw\n7W8DqwPb/wNcFti+Hv+CKgOmI8J1twWvDcwC3go75wOgYLBr4V/4BCADf6ALvvYBFwe2C4HtYZ+Z\nF7Z9O/DF8f77sK/x/bISghkPu1V1q6oqsBV4LrB/C/4MG/xzJz0MoKp/AwrD6tfDpyZWAhN2qeo2\neuZ7Pxs4OzBn0VvAkcD8wHt7VfWNCOl6Bzgr8BR9sqo2RzhmqYi8JCLvAFcARwX2rwQeC2yvCzt+\nsHSEm6aBla5UdS/+J/blgfPfVtWGIa51U6AU8CpQDhwR2O/FP0MtQBPQKf45/S8AgovvAByk596b\nSSopJ7czSa8rbNsHuMO2w/8mh7tKlTtsO/yc76nqfeEHBqqd2iJdRFXfF5EVwD8D3xGR51X1230O\nexA4T1U3i8jn8E+yOJR+6RiG+/EvrVoK/HKwa4l/jYSPAStVtVNE/gakB97uDAReVNUjIv+/vXt3\njSKMwjj8e1MZTIIQUoh1QAURvKS2tAoIKSwsxCoBEaNWgpDeNPkHrJIuIGJhxAuCgrZpBAlpRdhC\nJIVCimNxzsIkbDZZyOrCvk8z129uxXfmzAxnZmrdOeBujcPeYm42pJwh2KD6SN6Btzu8VuSz7h3g\nKF/ivAbuVClgJJ2RNNWtgbLM8Z+IWAOWgUsdVhsDfij/DXGrMf8L2clClk7u9Ti+S5psTD8HrgNX\nahvdtjUB/KxgcJbMVjqd30ngVES8Ah4AFxuLT5PVMG2IOUOw/2H/nWh0GF8CnknaJO/o2zXdXwLr\nkmaBewe1j4g3ks4Bn6sk8A7ZgXf7a9QF4Kmkdtay0GGdJ+RvT1s1HKv594FVSY/JjvvXIcfR2rfd\nTzQ6/4jYlfSe7OgPO6cNYF7SV+Ab+dhoz/Uo48ALSSfIjGCxsWwGeHTAdbEh4fLXZsdA0mhE/K7x\nm+QL5hs9tL9WbRZqeoR8TzAXEdt9OOTmvieAdxFxtZ/7scHnR0Zmx+NyfYq6CcwDD3tpHBEfgGlJ\n45LOA1vA234Hg3IbWPkH+7EB5wzBzMwAZwhmZlYcEMzMDHBAMDOz4oBgZmaAA4KZmRUHBDMzA+Av\nLu+EKe1w3jYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f214852d590>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.arange(10, 48, 3)\n",
    "indices = np.digitize(live.agepreg, bins)\n",
    "groups = live.groupby(indices)\n",
    "\n",
    "ages = [group.agepreg.mean() for i, group in groups][1:-1]\n",
    "cdfs = [thinkstats2.Cdf(group.totalwgt_lb) for i, group in groups][1:-1]\n",
    "\n",
    "thinkplot.PrePlot(5)\n",
    "for percent in [90, 75, 50, 25, 10]:\n",
    "    weights = [cdf.Percentile(percent) for cdf in cdfs]\n",
    "    label = '%dth' % percent\n",
    "    thinkplot.Plot(ages, weights, label=label)\n",
    "\n",
    "thinkplot.Config(xlabel=\"mother's age (years)\",\n",
    "                 ylabel='birth weight (lbs)',\n",
    "                 xlim=[14, 50],\n",
    "                 legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first and last points are not very reliable, because they represent fewer data points.\n",
    "\n",
    "It looks like there is a generally positive relationshop between birth weight and mother's age, possibly leveling or dropping for older mothers.\n",
    "\n",
    "We can get more information about the mothers by reading the respondents file, which contains one row per respondent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ReadFemResp(dct_file='2002FemResp.dct',\n",
    "                dat_file='2002FemResp.dat.gz',\n",
    "                nrows=None):\n",
    "    \"\"\"Reads the NSFG respondent data.\n",
    "\n",
    "    dct_file: string file name\n",
    "    dat_file: string file name\n",
    "\n",
    "    returns: DataFrame\n",
    "    \"\"\"\n",
    "    dct = thinkstats2.ReadStataDct(dct_file)\n",
    "    df = dct.ReadFixedWidth(dat_file, compression='gzip', nrows=nrows)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 7643 respondents and 3087 variables about each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7643, 3087)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resp = ReadFemResp()\n",
    "resp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we use the `caseid` variable as the index, we can look up respondents efficiently by id.\n",
    "\n",
    "Here's what the first few rows look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>caseid</th>\n",
       "      <th>rscrinf</th>\n",
       "      <th>rdormres</th>\n",
       "      <th>rostscrn</th>\n",
       "      <th>rscreenhisp</th>\n",
       "      <th>rscreenrace</th>\n",
       "      <th>age_a</th>\n",
       "      <th>age_r</th>\n",
       "      <th>cmbirth</th>\n",
       "      <th>agescrn</th>\n",
       "      <th>...</th>\n",
       "      <th>pubassis_i</th>\n",
       "      <th>basewgt</th>\n",
       "      <th>adj_mod_basewgt</th>\n",
       "      <th>finalwgt</th>\n",
       "      <th>secu_r</th>\n",
       "      <th>sest</th>\n",
       "      <th>cmintvw</th>\n",
       "      <th>cmlstyr</th>\n",
       "      <th>screentime</th>\n",
       "      <th>intvlngth</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>caseid</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",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2298</th>\n",
       "      <td>2298</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>27</td>\n",
       "      <td>27</td>\n",
       "      <td>902</td>\n",
       "      <td>27</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>3247.916977</td>\n",
       "      <td>5123.759559</td>\n",
       "      <td>5556.717241</td>\n",
       "      <td>2</td>\n",
       "      <td>18</td>\n",
       "      <td>1234</td>\n",
       "      <td>1222</td>\n",
       "      <td>18:26:36</td>\n",
       "      <td>110.492667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5012</th>\n",
       "      <td>5012</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>42</td>\n",
       "      <td>42</td>\n",
       "      <td>718</td>\n",
       "      <td>42</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>2335.279149</td>\n",
       "      <td>2846.799490</td>\n",
       "      <td>4744.191350</td>\n",
       "      <td>2</td>\n",
       "      <td>18</td>\n",
       "      <td>1233</td>\n",
       "      <td>1221</td>\n",
       "      <td>16:30:59</td>\n",
       "      <td>64.294000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11586</th>\n",
       "      <td>11586</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>43</td>\n",
       "      <td>43</td>\n",
       "      <td>708</td>\n",
       "      <td>43</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>2335.279149</td>\n",
       "      <td>2846.799490</td>\n",
       "      <td>4744.191350</td>\n",
       "      <td>2</td>\n",
       "      <td>18</td>\n",
       "      <td>1234</td>\n",
       "      <td>1222</td>\n",
       "      <td>18:19:09</td>\n",
       "      <td>75.149167</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6794</th>\n",
       "      <td>6794</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>15</td>\n",
       "      <td>15</td>\n",
       "      <td>1042</td>\n",
       "      <td>15</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>3783.152221</td>\n",
       "      <td>5071.464231</td>\n",
       "      <td>5923.977368</td>\n",
       "      <td>2</td>\n",
       "      <td>18</td>\n",
       "      <td>1234</td>\n",
       "      <td>1222</td>\n",
       "      <td>15:54:43</td>\n",
       "      <td>28.642833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>616</th>\n",
       "      <td>616</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>20</td>\n",
       "      <td>20</td>\n",
       "      <td>991</td>\n",
       "      <td>20</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>5341.329968</td>\n",
       "      <td>6437.335772</td>\n",
       "      <td>7229.128072</td>\n",
       "      <td>2</td>\n",
       "      <td>18</td>\n",
       "      <td>1233</td>\n",
       "      <td>1221</td>\n",
       "      <td>14:19:44</td>\n",
       "      <td>69.502667</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 3087 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        caseid  rscrinf  rdormres  rostscrn  rscreenhisp  rscreenrace  age_a  \\\n",
       "caseid                                                                         \n",
       "2298      2298        1         5         5            1            5     27   \n",
       "5012      5012        1         5         1            5            5     42   \n",
       "11586    11586        1         5         1            5            5     43   \n",
       "6794      6794        5         5         4            1            5     15   \n",
       "616        616        1         5         4            1            5     20   \n",
       "\n",
       "        age_r  cmbirth  agescrn     ...      pubassis_i      basewgt  \\\n",
       "caseid                              ...                                \n",
       "2298       27      902       27     ...               0  3247.916977   \n",
       "5012       42      718       42     ...               0  2335.279149   \n",
       "11586      43      708       43     ...               0  2335.279149   \n",
       "6794       15     1042       15     ...               0  3783.152221   \n",
       "616        20      991       20     ...               0  5341.329968   \n",
       "\n",
       "        adj_mod_basewgt     finalwgt  secu_r  sest  cmintvw  cmlstyr  \\\n",
       "caseid                                                                 \n",
       "2298        5123.759559  5556.717241       2    18     1234     1222   \n",
       "5012        2846.799490  4744.191350       2    18     1233     1221   \n",
       "11586       2846.799490  4744.191350       2    18     1234     1222   \n",
       "6794        5071.464231  5923.977368       2    18     1234     1222   \n",
       "616         6437.335772  7229.128072       2    18     1233     1221   \n",
       "\n",
       "        screentime   intvlngth  \n",
       "caseid                          \n",
       "2298      18:26:36  110.492667  \n",
       "5012      16:30:59   64.294000  \n",
       "11586     18:19:09   75.149167  \n",
       "6794      15:54:43   28.642833  \n",
       "616       14:19:44   69.502667  \n",
       "\n",
       "[5 rows x 3087 columns]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "resp.index = resp.caseid\n",
    "resp.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can join the tables, using the `caseid` from each pregnancy record to find the corresponding respondent and (abstractly) copy over the additional variables.\n",
    "\n",
    "So the joined table has one row for each pregnancy and all the columns from both tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7912, 3331)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "join = live.join(resp, on='caseid', rsuffix='_r')\n",
    "join.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The encoding for `screentime` is a colon-separated timestamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    19:56:43\n",
       "1    19:56:43\n",
       "2    14:54:03\n",
       "3    14:54:03\n",
       "4    14:54:03\n",
       "Name: screentime, dtype: object"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "join.screentime.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we convert to a datetime object, we avoid some processing problems later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0   2015-10-07 19:56:43\n",
       "1   2015-10-07 19:56:43\n",
       "2   2015-10-07 14:54:03\n",
       "3   2015-10-07 14:54:03\n",
       "4   2015-10-07 14:54:03\n",
       "Name: screentime, dtype: datetime64[ns]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "join.screentime = pd.to_datetime(join.screentime)\n",
    "join.screentime.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the effect of mother's age on birthweight, we can use a simple least squares fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.0328792906633737, 0.017995378469939118, 2.879260555190259)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ages = join.agepreg\n",
    "weights = join.totalwgt_lb\n",
    "\n",
    "inter, slope = thinkstats2.LeastSquares(ages, weights)\n",
    "inter, slope, slope*16*10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The slope is almost 3 ounces per decade.\n",
    "\n",
    "We can do the same thing using Ordinary Least Squares from statsmodels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>totalwgt_lb</td>   <th>  R-squared:         </th> <td>   0.007</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.007</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   57.60</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 07 Oct 2015</td> <th>  Prob (F-statistic):</th> <td>3.58e-14</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:40:17</td>     <th>  Log-Likelihood:    </th> <td> -12433.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  7912</td>      <th>  AIC:               </th> <td>2.487e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  7910</td>      <th>  BIC:               </th> <td>2.488e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    7.0329</td> <td>    0.061</td> <td>  115.940</td> <td> 0.000</td> <td>    6.914     7.152</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>   <td>    0.0180</td> <td>    0.002</td> <td>    7.589</td> <td> 0.000</td> <td>    0.013     0.023</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>381.292</td> <th>  Durbin-Watson:     </th> <td>   1.586</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1098.730</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.215</td>  <th>  Prob(JB):          </th> <td>2.59e-239</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.774</td>  <th>  Cond. No.          </th> <td>    119.</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            totalwgt_lb   R-squared:                       0.007\n",
       "Model:                            OLS   Adj. R-squared:                  0.007\n",
       "Method:                 Least Squares   F-statistic:                     57.60\n",
       "Date:                Wed, 07 Oct 2015   Prob (F-statistic):           3.58e-14\n",
       "Time:                        15:40:17   Log-Likelihood:                -12433.\n",
       "No. Observations:                7912   AIC:                         2.487e+04\n",
       "Df Residuals:                    7910   BIC:                         2.488e+04\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      7.0329      0.061    115.940      0.000         6.914     7.152\n",
       "agepreg        0.0180      0.002      7.589      0.000         0.013     0.023\n",
       "==============================================================================\n",
       "Omnibus:                      381.292   Durbin-Watson:                   1.586\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1098.730\n",
       "Skew:                           0.215   Prob(JB):                    2.59e-239\n",
       "Kurtosis:                       4.774   Cond. No.                         119.\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "\n",
    "formula = ('totalwgt_lb ~ agepreg')\n",
    "results = smf.ols(formula, data=join).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results object contains the parameters (and all the other info in the table):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7.032879290663371, 0.017995378469939035)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inter, slope = results.params\n",
    "inter, slope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the results are consistent with my implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.8792605551902457"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope * 16 * 10           # slope in ounces per decade"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use a boolean variable as a predictor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>totalwgt_lb</td>   <th>  R-squared:         </th> <td>   0.002</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.001</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   12.34</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 07 Oct 2015</td> <th>  Prob (F-statistic):</th> <td>0.000445</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:40:18</td>     <th>  Log-Likelihood:    </th> <td> -12456.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  7912</td>      <th>  AIC:               </th> <td>2.492e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  7910</td>      <th>  BIC:               </th> <td>2.493e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    7.5265</td> <td>    0.018</td> <td>  414.288</td> <td> 0.000</td> <td>    7.491     7.562</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>isfirst[T.True]</th> <td>   -0.0924</td> <td>    0.026</td> <td>   -3.513</td> <td> 0.000</td> <td>   -0.144    -0.041</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>364.721</td> <th>  Durbin-Watson:     </th> <td>   1.581</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1001.182</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.218</td>  <th>  Prob(JB):          </th> <td>3.95e-218</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.687</td>  <th>  Cond. No.          </th> <td>    2.57</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            totalwgt_lb   R-squared:                       0.002\n",
       "Model:                            OLS   Adj. R-squared:                  0.001\n",
       "Method:                 Least Squares   F-statistic:                     12.34\n",
       "Date:                Wed, 07 Oct 2015   Prob (F-statistic):           0.000445\n",
       "Time:                        15:40:18   Log-Likelihood:                -12456.\n",
       "No. Observations:                7912   AIC:                         2.492e+04\n",
       "Df Residuals:                    7910   BIC:                         2.493e+04\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           7.5265      0.018    414.288      0.000         7.491     7.562\n",
       "isfirst[T.True]    -0.0924      0.026     -3.513      0.000        -0.144    -0.041\n",
       "==============================================================================\n",
       "Omnibus:                      364.721   Durbin-Watson:                   1.581\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1001.182\n",
       "Skew:                           0.218   Prob(JB):                    3.95e-218\n",
       "Kurtosis:                       4.687   Cond. No.                         2.57\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "join['isfirst'] = (join.birthord == 1)\n",
    "\n",
    "formula = 'totalwgt_lb ~ isfirst'\n",
    "results = smf.ols(formula, data=join).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First babies are lighter by about 1.5 ounces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.4779760012473198"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.params['isfirst[T.True]'] * 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can make a model with multiple predictors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>totalwgt_lb</td>   <th>  R-squared:         </th> <td>   0.007</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.007</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   29.50</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 07 Oct 2015</td> <th>  Prob (F-statistic):</th> <td>1.72e-13</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:40:18</td>     <th>  Log-Likelihood:    </th> <td> -12433.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  7912</td>      <th>  AIC:               </th> <td>2.487e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  7909</td>      <th>  BIC:               </th> <td>2.489e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    7.0719</td> <td>    0.069</td> <td>  102.463</td> <td> 0.000</td> <td>    6.937     7.207</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>isfirst[T.True]</th> <td>   -0.0327</td> <td>    0.028</td> <td>   -1.185</td> <td> 0.236</td> <td>   -0.087     0.021</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>         <td>    0.0171</td> <td>    0.002</td> <td>    6.826</td> <td> 0.000</td> <td>    0.012     0.022</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>382.030</td> <th>  Durbin-Watson:     </th> <td>   1.585</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1103.586</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.214</td>  <th>  Prob(JB):          </th> <td>2.29e-240</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.779</td>  <th>  Cond. No.          </th> <td>    138.</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            totalwgt_lb   R-squared:                       0.007\n",
       "Model:                            OLS   Adj. R-squared:                  0.007\n",
       "Method:                 Least Squares   F-statistic:                     29.50\n",
       "Date:                Wed, 07 Oct 2015   Prob (F-statistic):           1.72e-13\n",
       "Time:                        15:40:18   Log-Likelihood:                -12433.\n",
       "No. Observations:                7912   AIC:                         2.487e+04\n",
       "Df Residuals:                    7909   BIC:                         2.489e+04\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           7.0719      0.069    102.463      0.000         6.937     7.207\n",
       "isfirst[T.True]    -0.0327      0.028     -1.185      0.236        -0.087     0.021\n",
       "agepreg             0.0171      0.002      6.826      0.000         0.012     0.022\n",
       "==============================================================================\n",
       "Omnibus:                      382.030   Durbin-Watson:                   1.585\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1103.586\n",
       "Skew:                           0.214   Prob(JB):                    2.29e-240\n",
       "Kurtosis:                       4.779   Cond. No.                         138.\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = 'totalwgt_lb ~ agepreg + isfirst'\n",
    "results = smf.ols(formula, data=join).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we control for mother's age, the difference in weight for first babies is cut to about 0.5 ounces (and no longer statistically significant)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.52389693286041505"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.params['isfirst[T.True]'] * 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relationship with age might be non-linear.  Adding a quadratic term helps a little, although note that the $R^2$ values for all of these models are very small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>totalwgt_lb</td>   <th>  R-squared:         </th> <td>   0.009</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.009</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   37.48</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 07 Oct 2015</td> <th>  Prob (F-statistic):</th> <td>6.30e-17</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:40:18</td>     <th>  Log-Likelihood:    </th> <td> -12425.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  7912</td>      <th>  AIC:               </th> <td>2.486e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  7909</td>      <th>  BIC:               </th> <td>2.488e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    6.0378</td> <td>    0.247</td> <td>   24.428</td> <td> 0.000</td> <td>    5.553     6.522</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>   <td>    0.0979</td> <td>    0.019</td> <td>    5.050</td> <td> 0.000</td> <td>    0.060     0.136</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age2</th>      <td>   -0.0015</td> <td>    0.000</td> <td>   -4.153</td> <td> 0.000</td> <td>   -0.002    -0.001</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>378.969</td> <th>  Durbin-Watson:     </th> <td>   1.587</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1074.412</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.219</td>  <th>  Prob(JB):          </th> <td>4.95e-234</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.751</td>  <th>  Cond. No.          </th> <td>1.36e+04</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            totalwgt_lb   R-squared:                       0.009\n",
       "Model:                            OLS   Adj. R-squared:                  0.009\n",
       "Method:                 Least Squares   F-statistic:                     37.48\n",
       "Date:                Wed, 07 Oct 2015   Prob (F-statistic):           6.30e-17\n",
       "Time:                        15:40:18   Log-Likelihood:                -12425.\n",
       "No. Observations:                7912   AIC:                         2.486e+04\n",
       "Df Residuals:                    7909   BIC:                         2.488e+04\n",
       "Df Model:                           2                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      6.0378      0.247     24.428      0.000         5.553     6.522\n",
       "agepreg        0.0979      0.019      5.050      0.000         0.060     0.136\n",
       "age2          -0.0015      0.000     -4.153      0.000        -0.002    -0.001\n",
       "==============================================================================\n",
       "Omnibus:                      378.969   Durbin-Watson:                   1.587\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1074.412\n",
       "Skew:                           0.219   Prob(JB):                    4.95e-234\n",
       "Kurtosis:                       4.751   Cond. No.                     1.36e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.36e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "join['age2'] = join.agepreg**2\n",
    "\n",
    "formula = 'totalwgt_lb ~ agepreg + age2'\n",
    "results = smf.ols(formula, data=join).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can combine the quadratic age model with `isfirst`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>totalwgt_lb</td>   <th>  R-squared:         </th> <td>   0.009</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.009</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   25.12</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 07 Oct 2015</td> <th>  Prob (F-statistic):</th> <td>3.61e-16</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:40:18</td>     <th>  Log-Likelihood:    </th> <td> -12425.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  7912</td>      <th>  AIC:               </th> <td>2.486e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  7908</td>      <th>  BIC:               </th> <td>2.488e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    6.0795</td> <td>    0.256</td> <td>   23.775</td> <td> 0.000</td> <td>    5.578     6.581</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>isfirst[T.True]</th> <td>   -0.0177</td> <td>    0.028</td> <td>   -0.637</td> <td> 0.524</td> <td>   -0.072     0.037</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>         <td>    0.0957</td> <td>    0.020</td> <td>    4.865</td> <td> 0.000</td> <td>    0.057     0.134</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age2</th>            <td>   -0.0015</td> <td>    0.000</td> <td>   -4.030</td> <td> 0.000</td> <td>   -0.002    -0.001</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>379.403</td> <th>  Durbin-Watson:     </th> <td>   1.586</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1077.436</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.219</td>  <th>  Prob(JB):          </th> <td>1.09e-234</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.754</td>  <th>  Cond. No.          </th> <td>1.41e+04</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            totalwgt_lb   R-squared:                       0.009\n",
       "Model:                            OLS   Adj. R-squared:                  0.009\n",
       "Method:                 Least Squares   F-statistic:                     25.12\n",
       "Date:                Wed, 07 Oct 2015   Prob (F-statistic):           3.61e-16\n",
       "Time:                        15:40:18   Log-Likelihood:                -12425.\n",
       "No. Observations:                7912   AIC:                         2.486e+04\n",
       "Df Residuals:                    7908   BIC:                         2.488e+04\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "===================================================================================\n",
       "                      coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "-----------------------------------------------------------------------------------\n",
       "Intercept           6.0795      0.256     23.775      0.000         5.578     6.581\n",
       "isfirst[T.True]    -0.0177      0.028     -0.637      0.524        -0.072     0.037\n",
       "agepreg             0.0957      0.020      4.865      0.000         0.057     0.134\n",
       "age2               -0.0015      0.000     -4.030      0.000        -0.002    -0.001\n",
       "==============================================================================\n",
       "Omnibus:                      379.403   Durbin-Watson:                   1.586\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1077.436\n",
       "Skew:                           0.219   Prob(JB):                    1.09e-234\n",
       "Kurtosis:                       4.754   Cond. No.                     1.41e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.41e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = 'totalwgt_lb ~ agepreg + age2 + isfirst'\n",
    "results = smf.ols(formula, data=join).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the effect is cut to less that a third of an ounce, and very plausibly due to chance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.28383607712278636"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.params['isfirst[T.True]'] * 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the best model I found, combining all variables that seemed plausibly predictive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>totalwgt_lb</td>   <th>  R-squared:         </th> <td>   0.046</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.045</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   54.63</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Wed, 07 Oct 2015</td> <th>  Prob (F-statistic):</th> <td>1.11e-76</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:40:19</td>     <th>  Log-Likelihood:    </th> <td> -12275.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>  7912</td>      <th>  AIC:               </th> <td>2.457e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>  7904</td>      <th>  BIC:               </th> <td>2.462e+04</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     7</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "            <td></td>              <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>            <td>    6.0829</td> <td>    0.243</td> <td>   25.022</td> <td> 0.000</td> <td>    5.606     6.559</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(race)[T.2]</th>         <td>    0.3625</td> <td>    0.031</td> <td>   11.649</td> <td> 0.000</td> <td>    0.302     0.424</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(race)[T.3]</th>         <td>    0.2434</td> <td>    0.050</td> <td>    4.833</td> <td> 0.000</td> <td>    0.145     0.342</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>nbrnaliv > 1[T.True]</th> <td>   -1.0779</td> <td>    0.121</td> <td>   -8.889</td> <td> 0.000</td> <td>   -1.316    -0.840</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>paydu == 1[T.True]</th>   <td>    0.0985</td> <td>    0.030</td> <td>    3.255</td> <td> 0.001</td> <td>    0.039     0.158</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>              <td>    0.0714</td> <td>    0.019</td> <td>    3.741</td> <td> 0.000</td> <td>    0.034     0.109</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>age2</th>                 <td>   -0.0012</td> <td>    0.000</td> <td>   -3.341</td> <td> 0.001</td> <td>   -0.002    -0.000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>totincr</th>              <td>    0.0132</td> <td>    0.004</td> <td>    3.452</td> <td> 0.001</td> <td>    0.006     0.021</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>453.110</td> <th>  Durbin-Watson:     </th> <td>   1.609</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td>  <th>  Jarque-Bera (JB):  </th> <td>1254.734</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 0.300</td>  <th>  Prob(JB):          </th> <td>3.45e-273</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 4.856</td>  <th>  Cond. No.          </th> <td>1.36e+04</td> \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:            totalwgt_lb   R-squared:                       0.046\n",
       "Model:                            OLS   Adj. R-squared:                  0.045\n",
       "Method:                 Least Squares   F-statistic:                     54.63\n",
       "Date:                Wed, 07 Oct 2015   Prob (F-statistic):           1.11e-76\n",
       "Time:                        15:40:19   Log-Likelihood:                -12275.\n",
       "No. Observations:                7912   AIC:                         2.457e+04\n",
       "Df Residuals:                    7904   BIC:                         2.462e+04\n",
       "Df Model:                           7                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "========================================================================================\n",
       "                           coef    std err          t      P>|t|      [95.0% Conf. Int.]\n",
       "----------------------------------------------------------------------------------------\n",
       "Intercept                6.0829      0.243     25.022      0.000         5.606     6.559\n",
       "C(race)[T.2]             0.3625      0.031     11.649      0.000         0.302     0.424\n",
       "C(race)[T.3]             0.2434      0.050      4.833      0.000         0.145     0.342\n",
       "nbrnaliv > 1[T.True]    -1.0779      0.121     -8.889      0.000        -1.316    -0.840\n",
       "paydu == 1[T.True]       0.0985      0.030      3.255      0.001         0.039     0.158\n",
       "agepreg                  0.0714      0.019      3.741      0.000         0.034     0.109\n",
       "age2                    -0.0012      0.000     -3.341      0.001        -0.002    -0.000\n",
       "totincr                  0.0132      0.004      3.452      0.001         0.006     0.021\n",
       "==============================================================================\n",
       "Omnibus:                      453.110   Durbin-Watson:                   1.609\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1254.734\n",
       "Skew:                           0.300   Prob(JB):                    3.45e-273\n",
       "Kurtosis:                       4.856   Cond. No.                     1.36e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 1.36e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = ('totalwgt_lb ~ agepreg + age2 + C(race) + '\n",
    "           'nbrnaliv>1 + paydu==1 + totincr')\n",
    "results = smf.ols(formula, data=join).fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All predictors are statistically significant, so the effects could be legit, but the $R^2$ value is still very small: this model doesn't provide much help for the office pool."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Logistic regression\n",
    "-------------------\n",
    "\n",
    "Let's say we want to predict the sex of a baby based on information about the mother.\n",
    "\n",
    "I'll start by creating a binary dependent variable, `isboy`, and checking for dependence on mother's age:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.692981\n",
      "         Iterations 3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>isboy</td>      <th>  No. Observations:  </th>  <td>  7912</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>               <td>Logit</td>      <th>  Df Residuals:      </th>  <td>  7910</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>MLE</td>       <th>  Df Model:          </th>  <td>     1</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Wed, 07 Oct 2015</td> <th>  Pseudo R-squ.:     </th> <td>1.352e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>15:40:19</td>     <th>  Log-Likelihood:    </th> <td> -5482.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>           <td>True</td>       <th>  LL-Null:           </th> <td> -5482.9</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  LLR p-value:       </th>  <td>0.7002</td>  \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -0.0038</td> <td>    0.104</td> <td>   -0.036</td> <td> 0.971</td> <td>   -0.208     0.200</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>   <td>    0.0016</td> <td>    0.004</td> <td>    0.385</td> <td> 0.700</td> <td>   -0.006     0.010</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                  isboy   No. Observations:                 7912\n",
       "Model:                          Logit   Df Residuals:                     7910\n",
       "Method:                           MLE   Df Model:                            1\n",
       "Date:                Wed, 07 Oct 2015   Pseudo R-squ.:               1.352e-05\n",
       "Time:                        15:40:19   Log-Likelihood:                -5482.9\n",
       "converged:                       True   LL-Null:                       -5482.9\n",
       "                                        LLR p-value:                    0.7002\n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -0.0038      0.104     -0.036      0.971        -0.208     0.200\n",
       "agepreg        0.0016      0.004      0.385      0.700        -0.006     0.010\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "live['isboy'] = (live.babysex==1).astype(int)\n",
    "model = smf.logit('isboy ~ agepreg', data=live)    \n",
    "results = model.fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated parameter is 0.0016, which is small and not statistically significant.  So the apparent relationship might be due to chance.\n",
    "\n",
    "But for the sake of the example, I'll take it at face value and work out the effect on the prediction.\n",
    "\n",
    "A parameter in a logistic regression is a log odds ratio, so we can compute the odds ratio for a difference of 10 years in mother's age:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0158023740467292"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_odds_ratio = results.params['agepreg'] * 10\n",
    "odds_ratio = np.exp(log_odds_ratio)\n",
    "odds_ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can use the odds ratio to update a prior probability.  A mother at the mean age has a 51% chance of having a boy.\n",
    "\n",
    "In the case a mother who is 10 years older has a 51.4% chance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.51391744178527499"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 0.51\n",
    "prior_odds = p / (1-p)\n",
    "post_odds = prior_odds * odds_ratio\n",
    "p = post_odds / (post_odds + 1)\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I searched for other factors that might be predictive.  The most likely candidates turn out not to be statistically significant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.692727\n",
      "         Iterations 3\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Logit Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>isboy</td>      <th>  No. Observations:  </th>  <td>  7826</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>               <td>Logit</td>      <th>  Df Residuals:      </th>  <td>  7820</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>MLE</td>       <th>  Df Model:          </th>  <td>     5</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Wed, 07 Oct 2015</td> <th>  Pseudo R-squ.:     </th> <td>0.0004249</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>15:40:19</td>     <th>  Log-Likelihood:    </th> <td> -5421.3</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>           <td>True</td>       <th>  LL-Null:           </th> <td> -5423.6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th> </th>                      <td> </td>        <th>  LLR p-value:       </th>  <td>0.4654</td>  \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>    <td>   -0.0639</td> <td>    0.111</td> <td>   -0.575</td> <td> 0.565</td> <td>   -0.282     0.154</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(race)[T.2]</th> <td>   -0.0135</td> <td>    0.054</td> <td>   -0.250</td> <td> 0.803</td> <td>   -0.120     0.092</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(race)[T.3]</th> <td>   -0.0158</td> <td>    0.089</td> <td>   -0.178</td> <td> 0.859</td> <td>   -0.190     0.158</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>agepreg</th>      <td>   -0.0056</td> <td>    0.006</td> <td>   -0.963</td> <td> 0.336</td> <td>   -0.017     0.006</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hpagelb</th>      <td>    0.0069</td> <td>    0.004</td> <td>    1.542</td> <td> 0.123</td> <td>   -0.002     0.016</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>birthord</th>     <td>    0.0297</td> <td>    0.024</td> <td>    1.264</td> <td> 0.206</td> <td>   -0.016     0.076</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           Logit Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                  isboy   No. Observations:                 7826\n",
       "Model:                          Logit   Df Residuals:                     7820\n",
       "Method:                           MLE   Df Model:                            5\n",
       "Date:                Wed, 07 Oct 2015   Pseudo R-squ.:               0.0004249\n",
       "Time:                        15:40:19   Log-Likelihood:                -5421.3\n",
       "converged:                       True   LL-Null:                       -5423.6\n",
       "                                        LLR p-value:                    0.4654\n",
       "================================================================================\n",
       "                   coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
       "--------------------------------------------------------------------------------\n",
       "Intercept       -0.0639      0.111     -0.575      0.565        -0.282     0.154\n",
       "C(race)[T.2]    -0.0135      0.054     -0.250      0.803        -0.120     0.092\n",
       "C(race)[T.3]    -0.0158      0.089     -0.178      0.859        -0.190     0.158\n",
       "agepreg         -0.0056      0.006     -0.963      0.336        -0.017     0.006\n",
       "hpagelb          0.0069      0.004      1.542      0.123        -0.002     0.016\n",
       "birthord         0.0297      0.024      1.264      0.206        -0.016     0.076\n",
       "================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = 'isboy ~ agepreg + hpagelb + birthord + C(race)'\n",
    "model = smf.logit(formula, data=live)\n",
    "results = model.fit()\n",
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, taking these parameters at face values, we can use the model to make predictions.\n",
    "\n",
    "The baseline strategy is to always guess boy, which yields accuracy of 50.8%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5079223102478917"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exog = pd.DataFrame(model.exog, columns=model.exog_names)\n",
    "endog = pd.DataFrame(model.endog, columns=[model.endog_names])\n",
    "\n",
    "actual = endog['isboy']\n",
    "baseline = actual.mean()\n",
    "baseline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`results.predict` uses the model to generate predictions for the data.\n",
    "\n",
    "Adding up the correct positive and negative predictions, we get accuracy 51.3%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.51316125734730389"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predict = (results.predict() >= 0.5)\n",
    "true_pos = predict * actual\n",
    "true_neg = (1 - predict) * (1 - actual)\n",
    "\n",
    "acc = (sum(true_pos) + sum(true_neg)) / len(actual)\n",
    "acc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can use the model to generate a prediction for the office pool.\n",
    "\n",
    "Suppose your hypothetical coworker is is 39 years old and white, her husband is 30, and they are expecting their first child."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.48485192])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "columns = ['agepreg', 'hpagelb', 'birthord', 'race']\n",
    "new = pd.DataFrame([[39, 30, 1, 2]], columns=columns)\n",
    "y = results.predict(new)\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is one of the few scenarios where the model predicts that the baby will be a girl.\n",
    "\n",
    "But given the low accuracy of the model (and the likelihood that it is overfit to the data) this model would not be much help in the pool."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}