Kawasaki Non-responder Model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import and clean data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kawasaki_raw = pd.read_csv('kd_treatment_response.csv', na_values=[' '], \n",
    "                           parse_dates=['dob', 'Date_admission'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Race2\n",
       "1    0.839465\n",
       "2    0.862385\n",
       "3    0.695652\n",
       "Name: responder, dtype: float64"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kawasaki_raw.assign(responder=(kawasaki_raw.Responder==1)).groupby('Race2').responder.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>Responder</th>\n",
       "      <th>StudyID</th>\n",
       "      <th>sex</th>\n",
       "      <th>dob</th>\n",
       "      <th>year_dob</th>\n",
       "      <th>month_dob</th>\n",
       "      <th>race</th>\n",
       "      <th>KD_complete</th>\n",
       "      <th>Race2</th>\n",
       "      <th>ethnicity</th>\n",
       "      <th>...</th>\n",
       "      <th>Ectasiaoraneurysm7</th>\n",
       "      <th>Ectasiaoraneurysm8</th>\n",
       "      <th>Ectasiaoraneurysm9</th>\n",
       "      <th>Ectasiaoraneurysm10</th>\n",
       "      <th>Ectasiaoraneurysm11</th>\n",
       "      <th>Totalectasia/aneurysm</th>\n",
       "      <th>age_months</th>\n",
       "      <th>age_years</th>\n",
       "      <th>age_norm</th>\n",
       "      <th>age_norm2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2</td>\n",
       "      <td>KD001</td>\n",
       "      <td>0</td>\n",
       "      <td>2004-12-06</td>\n",
       "      <td>2004</td>\n",
       "      <td>12</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>74.033333</td>\n",
       "      <td>6</td>\n",
       "      <td>0.958584</td>\n",
       "      <td>0.918883</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>KD002</td>\n",
       "      <td>1</td>\n",
       "      <td>2003-03-18</td>\n",
       "      <td>2003</td>\n",
       "      <td>3</td>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>19.200000</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.730244</td>\n",
       "      <td>0.533256</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>KD003</td>\n",
       "      <td>1</td>\n",
       "      <td>2009-01-09</td>\n",
       "      <td>2009</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "      <td>11.800000</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.958159</td>\n",
       "      <td>0.918068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>KD004a</td>\n",
       "      <td>1</td>\n",
       "      <td>2003-05-08</td>\n",
       "      <td>2003</td>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>17.733333</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.775416</td>\n",
       "      <td>0.601270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>KD005</td>\n",
       "      <td>1</td>\n",
       "      <td>1997-09-10</td>\n",
       "      <td>1997</td>\n",
       "      <td>9</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>...</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>79.033333</td>\n",
       "      <td>6</td>\n",
       "      <td>1.112580</td>\n",
       "      <td>1.237834</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 42 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   Responder StudyID  sex        dob  year_dob  month_dob  race  KD_complete  \\\n",
       "0          2   KD001    0 2004-12-06      2004         12     0            1   \n",
       "1          2   KD002    1 2003-03-18      2003          3     8            2   \n",
       "2          2   KD003    1 2009-01-09      2009          1     4            0   \n",
       "3          2  KD004a    1 2003-05-08      2003          5     3            0   \n",
       "4          2   KD005    1 1997-09-10      1997          9     4            0   \n",
       "\n",
       "   Race2  ethnicity    ...     Ectasiaoraneurysm7  Ectasiaoraneurysm8  \\\n",
       "0      3          1    ...                      0                   0   \n",
       "1      3          0    ...                      0                   0   \n",
       "2      1          1    ...                      1                   1   \n",
       "3      2          1    ...                      0                   0   \n",
       "4      1          2    ...                      0                   0   \n",
       "\n",
       "   Ectasiaoraneurysm9  Ectasiaoraneurysm10  Ectasiaoraneurysm11  \\\n",
       "0                   0                    0                    0   \n",
       "1                   0                    0                    0   \n",
       "2                   1                    0                    0   \n",
       "3                   0                    0                    0   \n",
       "4                   0                    0                    0   \n",
       "\n",
       "   Totalectasia/aneurysm  age_months  age_years  age_norm  age_norm2  \n",
       "0                      1   74.033333          6  0.958584   0.918883  \n",
       "1                      1   19.200000          1 -0.730244   0.533256  \n",
       "2                      8   11.800000          0 -0.958159   0.918068  \n",
       "3                      0   17.733333          1 -0.775416   0.601270  \n",
       "4                      2   79.033333          6  1.112580   1.237834  \n",
       "\n",
       "[5 rows x 42 columns]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kawasaki_raw.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate age from admission year and birth year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    74.033333\n",
       "1    19.200000\n",
       "2    11.800000\n",
       "3    17.733333\n",
       "4    79.033333\n",
       "Name: age_months, dtype: float64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kawasaki_raw['age_months'] = (kawasaki_raw.Date_admission - kawasaki_raw.dob).dt.days/30.\n",
    "kawasaki_raw['age_months'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kawasaki_raw['age_years'] = (kawasaki_raw.age_months/12.).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kawasaki_raw['age_norm'] = ((kawasaki_raw.age_months - kawasaki_raw.age_months.mean())\n",
    "                            /kawasaki_raw.age_months.std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kawasaki_raw['age_norm2'] = kawasaki_raw['age_norm']**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "relevant_cols = ['StudyID', 'Responder', 'sex', 'age_norm', 'age_norm2', 'age_years',\n",
    "                 'Race2', 'ethnicity', 'KD_complete', 'Tier']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>StudyID</th>\n",
       "      <th>male</th>\n",
       "      <th>age_norm</th>\n",
       "      <th>age_norm2</th>\n",
       "      <th>age_years</th>\n",
       "      <th>ethnicity</th>\n",
       "      <th>KD_complete</th>\n",
       "      <th>Tier</th>\n",
       "      <th>black</th>\n",
       "      <th>other</th>\n",
       "      <th>responder</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>KD001</td>\n",
       "      <td>0</td>\n",
       "      <td>0.958584</td>\n",
       "      <td>0.918883</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>KD002</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.730244</td>\n",
       "      <td>0.533256</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>KD003</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.958159</td>\n",
       "      <td>0.918068</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>KD004a</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.775416</td>\n",
       "      <td>0.601270</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>KD005</td>\n",
       "      <td>1</td>\n",
       "      <td>1.112580</td>\n",
       "      <td>1.237834</td>\n",
       "      <td>6</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  StudyID  male  age_norm  age_norm2  age_years  ethnicity  KD_complete  Tier  \\\n",
       "0   KD001     0  0.958584   0.918883          6          1            1     1   \n",
       "1   KD002     1 -0.730244   0.533256          1          0            2     1   \n",
       "2   KD003     1 -0.958159   0.918068          0          1            0     4   \n",
       "3  KD004a     1 -0.775416   0.601270          1          1            0     3   \n",
       "4   KD005     1  1.112580   1.237834          6          2            0     1   \n",
       "\n",
       "   black  other  responder  \n",
       "0      0      1          0  \n",
       "1      0      1          0  \n",
       "2      0      0          0  \n",
       "3      1      0          0  \n",
       "4      0      0          0  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kawasaki_subset = (kawasaki_raw[relevant_cols]\n",
    "                           .rename(columns={'sex':'male'})\n",
    "                           .assign(responder=(kawasaki_raw.Responder==1).astype(int),\n",
    "                                  black=(kawasaki_raw.Race2==2).astype(int),\n",
    "                                  other=(kawasaki_raw.Race2==3).astype(int))\n",
    "                           .drop(['Responder', 'Race2'], axis=1))\n",
    "kawasaki_subset.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "black\n",
       "0    0.820290\n",
       "1    0.862385\n",
       "Name: responder, dtype: float64"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kawasaki_subset.groupby('black').responder.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Probability of responder does not seem to vary with age. Wide variation at the end is mostly a sample size artefact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x115491b00>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAFkCAYAAADv13iSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4lfWd///nfZ8l2znZQxIgJGEJW0hIoC5grNbSatUW\nkVCkdWmp7bTTzV7OjLa/b8XpQMa2U62dobN0caRW6lKX0rG0KCqLiAaSECALkLBlIftycnK2+/79\ncSCCAieBhHPfJ+/HdXnpWe7k/fEked335/4siq7rOkIIIYQwJDXcBQghhBDiwiSohRBCCAOToBZC\nCCEMTIJaCCGEMDAJaiGEEMLAJKiFEEIIA7OGeoOu66xZs4ba2lrsdjtr164lKytr6PWXX36Z3/zm\nN8THx7N06VKWL18OwLJly3A4HABMnjyZdevWjVEThBBCiMgVMqi3bNmC1+tl48aNVFZWUlZWxvr1\n6wHo6uriySef5JVXXsHhcHDfffexaNEiUlNTAXj66afHtnohhBAiwoXs+i4vL6ekpASAwsJCqqur\nh147fvw4s2fPxul0oigK8+bNo6KigpqaGgYGBli9ejX33XcflZWVY9cCIYQQIoKFvKLu7+/H6XR+\ncIDViqZpqKpKTk4Ohw4dorOzk5iYGN555x1yc3OJiYlh9erVlJaW0tjYyP3338/mzZtRVbklLoQQ\nQoxEyKB2OBy4XK6hx2dCGiA+Pp6HHnqIb33rWyQmJjJ37lySkpLIzs5mypQpAOTk5JCYmEhbWxvp\n6ekX/D66rqMoyuW2RwghhIgoIYO6uLiYrVu3cvPNN1NRUUFeXt7Qa4FAgP379/PMM8/g9XpZvXo1\n3/ve93jxxRepq6vjkUceobW1FZfLRVpa2kW/j6IotLX1XX6LwigtzWn6NoC0w0gioQ0QGe2IhDaA\ntMNI0tKcod/EMIJ6yZIl7Nixg5UrVwJQVlbGpk2bcLvdlJaWAnDHHXcQFRXFl7/8ZRITE1m+fDkP\nP/wwq1atQlVV1q1bJ93eQgghxCVQjLR7ViScHZm9DSDtMJJIaANERjsioQ0g7TCS4V5Ry2WuEEII\nYWAS1EIIIYSBSVALIYQQBiZBLYQQQhiYBLUQQghhYBLUQgghhIFJUAshhBAGJkEthBBCGJgEtRBC\nCGFgEtRCCCGEgUlQCyGEEAYmQS2EEEIYmAS1EEIIYWAS1EIIIYSBSVALIYQQBiZBLYQQQhiYBLUQ\nQghhYBLUQgghhIFJUAshhBAGJkEthBBCGJgEtRBCCGFg1nAXEAkqTu1ld/MuPJZ+ogIOrsq8hvkT\nisJdlhBCiAhgmKDeVLeJwICF3ISpZDomhrucYas4tZfna5+l29ONatPRfApHexsAJKyFEEJcNsME\n9e6Tuwl4VI72NnJD1idME9ZbGjdT1VbJKVcrfnxYsTEhLp2kqGQJaiGEEJfNMEFd116H5ldoGziF\nw+bg9ulLw13SsLzTvIPGniP49QAoOugKAz0u7BZ7uEsbsS1H/8rrR/9Gv96NQ0nkpuwlfDL7U+Eu\na1ySz0IIcYZhgvpk30n0gEKPp4doa7RpgvpE/0l6vb14A150dBQU7BY7J/pPhru0Edly9K/8vPzf\naHW1EsCHBRvV7VUApgsIs48Z2HL0rzxV/Wtcvn4Ui44eaOJ43zHAfJ9Fc38TDT1HsHQFTHlrC8z/\n8yTMzzBBfbT7KAoWnHYHDpsj3OUMW+9gF+6A+5zn3AE3vYNdYaro0jy171c09h7Br33QM+DudfHU\nvl+ZKhwqTu3ltYY/AxAXa6djsGPosVn+uL566CW6Pd0A2C0WfJqfbk83rx56yVSfRXN/E28ef4PW\ngVYsUZopb20Fx6BspNvThWrV0fwKR3sbAfP8PJ0RCSdNEDntGAnDBHW7ux1Q6fP2kB6THu5yhq3H\n2zOi543qYOcBOgY68OMfes6KlYOdB8JY1cjtbt7FgY5qDnfV4dE9RClRTEvKIyU6xTR/WE/0HcOr\nefD4B8Gvg6YQZY3mxOmrarN4v2U3f6x/gfquWryaB7saxYykmaa6tbWlcTO1XQfpdHeiKX5U3Upy\nTDJbGjeb5ucJIuOkCSKnHSNlmKD2BXwA6FqA1oHWMFczfH7dP6LnjarV1XpOSAP48dPqMs9nAbD9\n5NtUtJYPjRkY0N1UtJajoPDVwq+Hu7xhURQLp1ytuHwuNAKoWIizxeFMnBHu0kZkY83v2NP6XvCz\nQMfFAHta38OuWk0T1HtOvc/x3mN4NR+qApoOLl8/dtVcY1Deb9lNdfu+c3oGzDYeCE63o2Mf3YNn\ntcNtvnaMlGGC+kywabpGh6cjzNWMP17dM6LnjepQV/15xwwc6qoPd2nDZlMtdHk6CWiBoee8mgeb\nagljVSNX0VaFZ+hzAB0IEKCirSrcpQ3bqYFW+ry9eDQfoAEqUaqNUya6mACoaNt73p4BM40HAqhs\n2zt08RBtteIJ+Gh1tVLZttdU7Ti7+/62vNtCvt8wQX2Ghka3uzvcZQiT6hhsP++YgY7B9jBVNHKd\ng11oAQ2f7hs62VB0hU6TjXsY8Lnwa350tKHnFFQGfK4wVjUy3oAPl89FQP/gpMmvePGe7gE0i0Nd\n9RzvOxbsuTw9BsXl7yfOap7xQACdg50jet6ImvubqGqvBCBBjxnWMSGDWtd11qxZQ21tLXa7nbVr\n15KVlTX0+ssvv8xvfvMb4uPjWbp0KcuXLw95TCg+zTvs9wpxtl7P+ccGXOh5IzrRdwyvHrwSBdDR\n8epe092jBh2NwIeeCXaDm4Vf86HpGn7tg9tCiqrg18wV1J3uDga8rnOmkfoCXjrd5uq9TIpO5tTA\nqeCMCK+OHlCIsznIjs8Od2nD1tBzhCM9R6jrrMGvDHL7zNtDHhMyqLds2YLX62Xjxo1UVlZSVlbG\n+vXrAejq6uLJJ5/klVdeweFwcN9997Fo0SL2799/wWOGI/ChX24hhutCPztm+pnq8fQOhfQZOjo9\nnt4wVXRpbOr5/7xc6Hkj8gS8BLQACspQ70ZAC+AJmOtiwq/78GpeBgNeznThR1vs+HVznXBkO7N5\nt/ldTg00E8CPBSsTYjPJdponqKvaKylveQ+A6Kjh/S6EfFd5eTklJSUAFBYWUl1dPfTa8ePHmT17\nNk6nE4B58+ZRUVFBVVXVBY8RQlycTz9/CFzoeaOyW6KGAu6M4JiBqDBWNTJneve00933OjoWLKbr\n9QtoGrquY1VUUFTQg72lAU0LfbCB9Hi66fZ04PYN4MePFSvdng56POa5XXq0p5FuT+fQeIHhCBnU\n/f39Q0EMYLVa0TQNVVXJycnh0KFDdHZ2EhMTwzvvvENubu5FjxmutDRn6DcZXCS0AaQdRmKmNuho\n5+0Z0NFM1A4d7UPtCIa2bqI2QFxMDFaXFZ/fF/xbrKhEWaKIi4kxVTuquveioxFli8Km21AVFR2N\nqu69pmmHS+/hWH8jbp/7nAGjFxMyqB0OBy7XB4M/zg7c+Ph4HnroIb71rW+RmJjI3LlzSUpKwul0\nXvCY4Wpr6xvR+40oEtoA0g4jMVMb/FoAFfV0OJ8eFIeCXwuYph2DgcHznmwMBgZN0wYAVQteLOmn\nm6Lrp/8ua1ZTteNo53HQVaLUGKxWFb9fAz34vFnacbzzJP2D/Xg170d+ti4kZHoWFxfz1ltvAVBR\nUUFeXt7Qa4FAgP379/PMM8/w+OOP09DQQHFxMUVFRRc8RghxcdYLnD9f6HmjUlAAzhkUd/bzZqDr\n+kfqVVDQdfMMiINgKGu6hqYFCGgBNC1w+rG5ur5jLOcfJX2h542ox9tDQNewKNZh7wkR8jd/yZIl\n7Nixg5UrVwJQVlbGpk2bcLvdlJaWAnDHHXcQFRXFl7/8ZRITE897jBBXggXLeQeOWTDPHGS71Y7f\n/9F7V3aruRbZsKrWYNfkWaGmKipWEw0msyjBNpx94aMqKhbFPG0A6Pf149N8aDqgKKAo+DQf/b7+\ncJc2IvMnFPNO8w48gUFAx6paiLJEM39CcbhLGzarasVhczAYcA+NfQh5TKg3KIrCo48+es5zubm5\nQ//9zW9+k29+85shjxkJ1UR/VD88WObs58WVp6qW8973UU20WIhNjUJlEA2dYEIoqCjYVPMMwgKI\nscag6RoMLXeioOkaMVbzXP0k2ONx+wfOmUdtUSwk2OPDWNXIuf0DoAfvUYOG3Q5WJSr4vInckbec\nHm83Tf0nCSg+LLqNiY5J3JG3PNylDVtm3ERcPhfR1mD3/XAY8rTQrtrCXcKwWRTLeZcLtSjmCYZI\nEm+Pp2MwODf0TDyced4skqIS8WtePIEPVoWLskSRFJUYxqpGzqrYsKk2fJoPTp+42lQbVsU8v98z\nU+bQ5+vDo3k5c7IRpdqZmTIn3KWNiMfnZ9DrRcECWFA0Fb/iN9GM9qD5E4q4L/8rpt7N7KbsJbh8\n/XQOdqIxSqO+r5TgFaiCTbWSEpsW7nKGzW65QDelyfajtik2fOeZU2kz0R9VgILUQna37MLtd6Oj\no6IQY42hILUw3KUN25yUuexo6kJRVBQlOPDHoliYkzI33KWNiEVVibbGYNPt57TDMsKBpeF03aTr\ncfvc1HeftbFI4kyum3R9uEsbtqZ2F709KhbFjtWm4fH50XRwWKJNtVPhGfMnFDF/QhFpaU7TDCA7\n25LsT+PyuajvqsWnDg7rGMMEtcPu4MzZ6jUZ14a7nGFLiU7F62oioAenbICCRVFJiU4Nd2kjEmuN\no8f30bmIsda4MFRz6fLTChjwD9A60IofL1bspMemk59WEO7Shu0T2Us42ttIs6sJn+7DrtrIjJvI\nJ7KXhLu0EUmJSWXQ72YgMMiZRTZiLdGkxJjnd2Nq4jQSohNZmHE10VFWBj3+oefNoKffw+PPVRLt\nTyc6wYvFptEz4MHn04i2OslJyA39RcSoynRMZOn0ZcG1vmNHaXrWlZIYnTi0JeHX5v99uMsZtvzU\nArwBD/2+fjQ0VFQcNgf5qeYJBoAp8TnUd9ecs1SiVbUyJT4nfEVdgvlpRQwGBs/ZXScxOon5aebp\nGrNb7BSmF5PaN4GA4sWi25nknGy6Xpr5aUX0eHrAp6ApAVTdgsPmNNVnEWuNZUH6wuDVjzJIvN3J\njKSZxFpjw11aSG6Pnyeer6Kjd5CPz1rEYMwhGnsbsCn9dHbppJDDgvSPhbvMcSnTMZFMx8Rhz/02\nTFDfN/8+c95vmPJJXH4Xne72oT9GyTGp3DTlk+EubUQ+nXsz/XW9dA92ESCABQuJ0Ul8OvfmcJc2\nIgszrqLf76LV1YLFrhHwqqTHZbAw46pwlzZsvZ5e0mMzSI/NIC7WjmvAO/S8mZRk3cD+jmoURUVT\nfKi6jZToZEqybgh3acPW7+tjauK04JV1fAw9ve6h543MH9D45SvVHG3t4/rCTAqLP8NfGv+P3MTp\nxMRYeeXtI3h7dIrSJKjNwDBB/c83/rM57zfk3IzLP0BdVy1+dRCrFk1e0kyW5Jgr4K6bdD06UHVq\nL25cxBBHwYQiU92Lg+CZ6g2TbxzqVgoMWMhNmGqqTeXjoz4YEPfh580k1hrLTdmfGroXZ9OiTXM1\neobD5qTvPKHssBl3FSxd13l6cy3VRzopmJbC3Z+eiUVVURQlOAhL7WfOxMk0HcrA3zUZzPOrMW4Z\nJqjN6sP3G8wYDAC5CVPp8/Uxf0LROVcOuQlTw1zZyJ3drWTGk78pzmx0oNXVgkJwOlN6XAZTTLTx\nAJj3avRsuQlTh7Yk/PDzRvWnHY1sr2omO8PJ331u7tDgvbMHYVVPbOXhul1sq2zi2rkZYa5YhCJB\nPQrMHgzA0IlFQ88RFCWA0+Y05QlHJDhz0pQcnWzqkyYzXo1+mNl+L7ZVNfHy9gZSE6L5bmkh0fbz\n/4lPT44lLyuRmmPdnOoaYEKSeXo5xiMJajEkEk44IoHZwuFCzHg1ej5m+b2oPtLB/75WS1y0lQdW\nFJIQd/HBhyUFmdQd72b7vhaWXW+uz2S8kaAWwoDMEg4XEyknHGZwtKWP/3i5GlVV+PbyAjJTQk+r\nXDhrAr/fUseOfc0svS4XVZXVFI1KgloIMWYyHRNx9zrIy0jB7/nogjri8rX3uHni+Uq83gBfX5rP\njMnDW8Euymbh6tnpvFnRRHVDcOCZMCbzLBEkhDCVXpeX/3ylmn95+n2+/bM3aWwx1/QyM+h3+3j8\nuUp6XF5WfnIGC2dNGNHxJYXB3o1tVU1jUZ4YJRLUQohRpes6O/Y184P/2cXug6eYmBpHT7+Hx57Z\ny74jH512Ji6Nzx/g31+sorljgE99LIslC7NG/DVyMpxMSoujor6d3tPz9YXxSFALIUZNW7ebnz1X\nya//fBB/QOeuT87gn798FQ/d8zE0Xefnz1exvao53GWanqbr/GrTQepO9PCxWRNY8Ynpl/R1FEWh\npGAiAU1nV3XLKFcpRosEtRDismmazl93H+P//fpd9jd0kj81mR995SqWLMxCVRUWFUzkwZXziYmy\n8Jv/O8ifdjScs0+1GJnntx7ivZpT5E1O4Cu3zUZVLn0g2LVz07GoCtuqmuUzMSgJaiHEZTlxqp+1\nG8rZ+MYh7FYL9982hwdKC0lNOHff6RmTE/n+3QtIiY/mpW0NPL25loCmhalq8/rbe8fZvPs4mSmx\nfPPOAmzWy9tS1xlrpygvjZPtLo40yzgCI5KgFkJcEp9f449vH+HRp96jobmXa+ak8y/3X821+Rko\nF7jCy0yJ4wf3LGBKuoO3Kpr49xf34fEObwchAe/XnGLj6/UkxNl5oLQQR8zobEN7fUEmANsq5baE\nEUlQCyFGrP5EN2t+u5tNOxtJcNj5zvICvvrZucTHht7hK9ERxT+tKmZubjKVhzv48bN7ZSDTMNSf\n6Oa//3QAu93Cd0sLSU2MCX3QMM3JSSY5PordB1vlxMmAJKiFEMPm9vjZ8Ndayn63h5aOAW4qnsyP\nVl9N4fSR7TEdE2XlO8sLWJSfQUNzL+s2lHOqa2CMqja/5g4XT75Qhabp/P3SfLIzRncZVlVVWJyf\nyaA3wHs1p0b1a4vLJ0EthBiWykPt/H+/epete06SmRLLw19cwBc+lUdM1KWtm2S1qKy+dTa3Lcrm\nVJebtRvKaZB7pB/R0+/h8ecqcQ36ufeWmeRPHZuFSa473f29XeZUG46sTCaEuKhel5ffb6lj98FT\nWFSFzy7O4dZrc7BZL/88X1EUll0/jSRnNL/7ay2P/X4PX/9c/oiv0CPVoNfPEy9U0d4zyOeuy6Wk\nYOyWX01LjGF2dhIHj3bR0jlARrJs1GEUckUthDgvXdfZWf3BwiVTJ8bzyJc+xtKSqaMS0me7sWgS\n31w2D3T4xYv7eLtSruoCmsYvX97P0ZY+rivI5LOLc8b8e5YUnh5UJlfVhiJBLYT4iPZuN48/V8mv\nNh3EF9C466YZfP+LC5ic5hiz71k0I41/uKuI2GgrT71Ww8vbjozbeb26rrNhcy37jnSQPzWZez49\n84Ij6UfTgrw04qKt7NzXIlPnDESCWggxRNN0/vbecf7fr3dT3dDJ3Nxk/mX11Sz5WNYV2V1p2qQE\nvn/3AlITonl1RyNPvVaDPzD+AmPTzkbermxmSrqDr38uH6vlyvyptlktXDMngx6Xl6rDstyrUUhQ\nCyEAONHWz7rflfPs6/VYLQpfuW0231sxutOAhiMjOZYf3LOQ7Awn26qa+cWL+xj0+q9oDeG0vaqZ\nl7Y1kBIfzXdLCy95sN6luk7mVBuOBLUQ45zPr/HytiM8+tv3ONLUy9Vz0ll7/zUsys+8It2t55MQ\nZ+efVhUxb2oK+4508Njv99Ljivy51tUNHfzvX2qIi7bywIpCEh1RV7yG7AwnU9IdVB3uoKffc8W/\nv/goCWohxrFDJ3pY89vdvLqjkfg4O99eXsDXPjuX+LjQC5eMtWi7lW/dOY/r5mVytKWPdRvep6Uz\ncudaH2vt4z9eqkZRFL51ZwETU+PCVktJwUQ0XWenbNRhCBLUQoxDbo+fZ/5aR9nvymnuGOATxZP4\nl69czXyDTYuyWlS+9JlZfHZxDm3dg6zbUM7hkz3hLmvUdfQM8vjzlXi9Ae6/fQ55WYlhreeauelY\nLSpvy0YdhiBBLcQ4U3W4nf/363d5fc8JMlJiefiLxXzxUzOv+L3Q4VIUhaUlU7n35pkMDPr5ybN7\n2VvfFu6yRo1r0Mfjz1fS0+/l85+YzsdmTQh3ScRF21g4M43WzgHqT0TeiZHZSFCLc+i6jqbJGXQk\n6h3w8t+v7ueJ56vo6fdy+6Ic1nzpY8yYHN6rt+H6+PxJfOvOeaDAv/9xH2/uPRnuki6bz6/xixf3\n0dTuYsnCLD511ZRwlzSkpEDmVBtFyFNoXddZs2YNtbW12O121q5dS1ZW1tDrr776Kk899RQWi4Vl\ny5Zx1113AbBs2TIcjuCcy8mTJ7Nu3boxaoIYDcda+9hZ3cKu/S0EdCiYmkzRjDTmTU0hyn552+iJ\n8NJ1nV37W3n29Xr63T5yM+P50i2zmDxh7OZEj5XC6an8413F/PyFSp7eXEtn3yB3lEwN26C3y6Hp\nOr/+8wHqjnezcGYan79perhLOsfM7CRSE6J5r+YUqz556UvFissX8v/8li1b8Hq9bNy4kcrKSsrK\nyli/fv3Q6z/+8Y957bXXiI6O5tZbb+W2224jKio4UvHpp58eu8rFZevp9/DO/lZ2Vrdwoq0fAEeM\njRi7hXf2t/LO/lZsVpW5OckU5aUyf3oqzmHsjiSMo73HzdOba6k+0ondprLyphl8csHkKzIneqxM\nnRjP9+9ewON/qGTTzqN09Xq495ZZV2yu8Wh54c3D7D54iumTE7j/9jmoBjvZUBWF6woyeXlbA+/V\nnOL6wrFbvlRcXMigLi8vp6SkBIDCwkKqq6vPeX3WrFn09PQMndEqikJNTQ0DAwOsXr2aQCDAAw88\nQGFh4RiUL0bK6wtQcaidHftaqG7oQNfBoioUzUhl8bxMCqalkJEez/vVTeypa2dvXRsVh9qpONSO\nokDe5ESK8tIonpF6xefXiuHTNJ3X95zgj28dweMLMDcniXtunkVahHxm6UmxfP/uBfz8hUp2VLfQ\n7fLyjaX5prnq2/L+cf7y7jEykmP59p0F2KzG7LW6bl4mr2xrYFtlkwR1GIX8qe7v78fp/GBLNavV\niqZpqGrw7HXGjBnceeedxMbGsmTJEhwOB9HR0axevZrS0lIaGxu5//772bx589Ax4srSdZ36Ez3s\nrG7hvZpTuD3BxSNyM50sys/kqtkTzrlSVhSFnIx4cjLiWXb9VFo7B9hT38beunbqjndTe7ybja/X\nM2WCg+K8NIry0picFmfK7sdIdLKtn6deq+FwUy9x0Va++KnZLMrPiLjPJz7Ozj/eVcwvX6mm6nAH\nj/1+Dw+UFpIQhrnHI1Fe28azW+qJj7PzwIpCHDG2cJd0Qcnx0cydmkz1kU5OtruYFMYpY+NZyKB2\nOBy4XK6hx2eHdG1tLW+++SZvvPEGsbGxPPjgg2zevJkbb7yR7OxsAHJyckhMTKStrY309PSLfq+0\ntNHdYzUcjNSGlg4XW8tPsPX94zR3BD/DlIRobl2cyycWZpGVfuFaz25HWpqT/Jnp3AN09Q7y7v4W\n3qlupqq+jWPb+3l5ewMZKbFck5/JNfmZzMpJxmKQrlUjfR6Xarht8PkDPP96Pc+/Xoc/oFMyfxL3\nL80nyRk9xhUOz1h9Fv/8tUX88o9VbN51lLLf7+XR+69h8oSx+V6X24aDDZ38z5/2E2W38Oj91zI9\nTNOwRtKO266bRvWRTsrr25k/O2MMqxq5SPj9Ho6QQV1cXMzWrVu5+eabqaioIC8vb+g1p9NJTEwM\ndrsdRVFITk6mt7eXF198kbq6Oh555BFaW1txuVykpaWFLKatre/yWhNmaWnOsLfB7fHzXs0pdla3\nUHe8GwC7TeXaueksmpfJ7ClJQ/cnL1RrqHYsmJ7CgukpuD1+qg53sLe+jarDHbz81mFefusw8bE2\n5s9IpWhGGnNyksLWrWeEz+NyDbcNh0728NRrNTS1u0hyRnH3p2Yyf0Yq/kEfbYO+K1DpxY31Z7Hi\n41OJsam8vK2BB3/+Nt9ZXsj0yQmj+j0utw0tnQOs21COP6Dz7TvmkRBtCcvP50jbkTshDkeMjS27\nj/GZq7IMMxYgUn6/hyNkUC9ZsoQdO3awcuVKAMrKyti0aRNut5vS0lJWrFjBqlWrsNvtTJkyhTvu\nuANd13n44YdZtWoVqqqybt066fYeQ5qmc6Cxkx3VLeypa8PnD25iMGtKIovyM1kwM21M7t3FRFm5\nek46V89Jx+fXOHi0iz11bVTUt/F2ZTNvVzYTZbcwb2oKxXmpFExNJTbaHPcQzcLt8fPHt4/wRvkJ\ndILbRS6/YZpp7tWOFkVR+OziXJKcUfzva7X8ZONevnr7XBbMDH2BcCX0uLz87A8V9Lt93HfLLAqm\npYS7pGGzWVWunZvB394/TkV9OwsNMM97vFF0Ay07EwlnR1eyDSfa+tlZ3cI7+1vo6Q+ug5yeFMOi\n/Ayuzc8gNeHSBg5dbjs0TedwUw9769rZU9fGqW43EBy0Njs7iaK8NIpmpI75OsaRcsZ9oTZUHe5g\nw+YaOno9ZCTHct8ts8K+otWFXMnPYt+RDta/VI3XF2DVkjxuWjB5VL7upbbB4w3w2O/30NjSx2cX\n57C0ZOqo1HOpLqUdJ9r6+eGvd1MwLYXvlhpjYHCk/H4Px/g67Y4AvQNe3j09pepoa/CHNDbKyg1F\nk1iUn8G0ifFhHzSkqgozJicyY3IipTdO42S7i711beypa6e6oZPqhk42bK5l2sT4ocFoGcmxYa3Z\nTHoHvGx8vZ5d+1uxqAq3Lcrm9kU5hh05fKXNm5rCP32hiCeeq+SZv9XR2TfInR+fFpbpTwFN45ev\nVNPY0sfieRl87rrcK17DaJic5iA3M559Rzro6vOQ5DT2gL1II0FtAj6/RuWhdnZWt7DvSAcBTUdV\nFAqnpbB4XiaF01MM+0daURQmpzmYnObg9sW5dPQMnh5B3kbd8R4ON/Xy/JuHmZgaR9GMVIrz0sjJ\ncIb9ZMOJ8Ip/AAAgAElEQVSIdF1n14FWnt1yZuESJ/fdMpssEy5cMtZyMuL5/j0LefwPFby26xjd\nfR6+9JnZV/T+qq7r/O6vdVQd7mBubjL33jzL1D/XJYWZNDT3sn1fM7cvygl3OeOKBLVB6brOkeZe\ndu5rYffBVlyDwSlVU9IdLMrP5Jo56YbY4WikUhKiWbIwiyULs+h3+6g8FOwer27o5M/vHOXP7xwl\nyRlF8Yw0ivJSyctKNMzglXBq73GzYXMd+450YLeqfP4T01myMMvUC5eMtQmJMXz/7gU8+UIV7+xv\npcfl5e/vmHfF7t//+Z2jvFXRxJQJDr6xNN/0P8dXz05n45Z6tlc1ceu12YZboCWSSVAbTEfPIO/s\nb2FndcvQln4JcXZuvmoKi/IzTLns44U4YmwsnpfJ4nmZeLwBqhs62VvfRuWhdl7fc4LX95wgLtpK\nwbTglXZ+bvK4W840oOlsef84L55euGTO6YVLJkTIwiVjzRlr58G7ivjvV/ezt76df31mD98tLRzz\nrtsd+5r549tHSImP4julhRExuC8mysrCWRPYWd1C7bFuZmcnhbukccP8Pz0RYNDrp7y2jZ3VLdQc\n7UInONLyqtkTWJSfydzcJCwRPmo+ym5hwcw0FsxMwx/QqD/ezZ66dvbUt/HO/uCAObtVZW5ucA3y\n+TNSDb1QxHDpuo7Xr+H2+BkY9OP2BP8ZOP3P7oOnqDnaFdELl4y1KJuFv79jHs/8rY6te0+ybsP7\nfHfF/DFbvGN/YydPvVZDbJSV766YH1H3c0sKMtlZ3cK2qiYJ6itIgjpMNE2n5lgXO/a1UF53Cq8v\nOKUqb3ICi+ZlsnDmhHE7lclqUZmdk8zsnGRWLZlBY0sfe+uDg9H21gf/URSYmZVI0eku8ksd4X65\n/AHtnHB1D/oZ8AQY8PhwewIMDAb/PfT6Oe8LPg6E2K3sqtkTuOuTeSSY8FaHUaiqwhc/lUdyfBQv\nvnWEsg3lfHt5waiPkj/W2sd//HEfigLfunNexK3klZeVSHpSDOW1bQws8REbbf6TZTMYn0kQRs0d\nLnZWB7u2u/o8AKQmRLMoP4NF+RlMSJLRz2dTFIXczHhyM+NZdv00WjoHgiPI69uoOdZNzbFunn29\nnux0J0V5wS7ySanDW85U03UGQ4TomX+ffcU7cNbzZ06wRsJuU4mJsuKMtZGeFENMlJXYaCsxUcF/\nYs/6d97UFFJi5Y/haFAUhVuvzSHREcVTr9Xw040VfPX2OaM2L7izd5Annq9k0Bvg7z43l5lTIu+K\nUzm9UceLbx3h3QOt3Fg8OlPfxMVJUF8B/W4fuw+2smNfCw3NvQDERFm4vjCTRfmZzJicIN2Zw5SR\nHMst12RzyzXZdPd7qKgPdo8fbOziaGsfL29rYEJiDIXTU0lMiKaja+BDIRvA7fEx4Akw6PEz0kUE\nLKoyFKyJjihizwrWcwPXQmyUjdgoCzHR575nJIOKImGuqNEsnpdJgsPOf7xUzS9frmblJ2ewZGFW\n6AMvYmDQx+PPVdLd72XFjdO5avbFl0s2s8XzMnnp7QbermqWoL5CJKjHiD+gse9wBzurW6g41E5A\n01EUyJ+azOL8TIpmpGK3ja+BUaMt0RHFDUWTuKFoEgODfvYd6WBPXRtVRzr42/vHP/J+BYbCMiU+\nmtihAA2Gaky05SNXtGeHbGyUFZtVlZOqCJCfm8JDq4p54vlKnt1ST1evh+U3Xtpca59f49//uI+T\n7S4+uWAyn77q8kLf6BIdURRMS6HiUDvHWvuYcpE9A8TokKAeRbqu09jSy459Lbx7oJV+d3CN5Ulp\ncSzOz+TqOekRNbDESGKjz13O9EhTD0lJcXjd3qGr3Ci7RaaUiCHZGU5+cPcCHn++kr/sPkZn3yCr\nb52DzTr8Hg9N1/nN/x2k5lg3xXlprLxpxrg4kbuuIJOKQ+1sq2rmC0skqMeaBPUoqT3WxZqn3uNY\nS7Cb0hlr45MLJ7M4P5Mp6Y5x8ctrFDaryswpSdJtLEJKTYzh4S8u4MkXq9h98BS9Li/fXDZv2IOk\nXnzrMO8eaGXapHi+evuccTOvvWBaCvFxdnbtb2HFjdMMu+BSpJCgHgXvHmjl138+gK7DwplpLMrP\nJH9qsukXOBBiPHDE2Hjw8/P57z8dYE9dG2XPBPe1To6/+Pagr5ef4LVdx0hPjuXbdxaMq1tZVovK\novwM/vLuMfbWt0f0PXkjkCS5DLqu85d3j/Ffr+7HZlV59P5r+cYd85g/I1VCWggTsdssfGNpPjcV\nT+Zkm4u1G8o50dZ/wffvqWvj93+rIz7WxgMrCnHGjr+pcyUFmQBsq2wKcyWRT9LkEmmazu+31PPc\n1kMkOaN46AsLKMwzxpZ6QoiRU1WFVUtmUHrDNLr6PJT9bg+1x7o+8r7DJ3uCJ+c2le+UFo7bVeIy\nU+KYPjmBA41dtPe4w11ORJOgvgReX4D1L1fzevkJJqXF8YO7F8jGCEJEAEVRuOWabL56+xy8vgD/\n9ocKdh9sHXq9tXOAn79QhT+g8fXP5ZObGR/GasOvpCATHdhe1RzuUiKaBPUI9Q14+cnGveypa2PW\nlEQe/kJxyHtZQghzuWZuBg+sKMRmVfnPV/azeXdwB67Hn6uk3+3j7k/PpHB6arjLDLuPzZpAlN3C\njn3NaCFW2BOXToJ6BE51u1m3oZzDJ3u5Zk46D6yYL0voCRGh5uQk89AXFpDosPOHNw7xrX/byqlu\nN7ctyuGG+ZPCXZ4hRNutXDVrAh29Hg4e/ehtAjE6JKiHqaG5l3VPv09rl5vPXJPNV24f2XxLIYT5\nZE1w8IO7FzIxNY7uPg+L8jO4oyQ33GUZSknhRAC2VcmgsrEi07OGofJQO798pRqfX+OLn8rjE7Js\nnhDjRkpCNN//4gKaugfJSYuVNRE+ZNrEeDJTYtlT10a/2xcRu9oZjVwShvBmxUmefLEKdPjmHfMk\npIUYh2KjrVw7L1OmXZ6HoiiUFEzEH9B5Z39LuMuJSPJTdwG6rvPHt4/w9F9qiYu28Q93FVEk06+E\nEOIjFuVnYFEVtlU2o+syqGy0SVCfhz+g8es/H2TTzkYmJMbwg7sXMG1SQrjLEkIIQ4qPs1M4PZUT\nbf00tsiyvaNNgvpD3B4/P3++kp3VLeRmxvP9uxeQnix7RAshxMWcWalM5lSPPgnqs3T1efjXZ/aw\nv7GL+dNT+ce7ioiPG39LAwohxEjlT00m0WFn14FWvL5AuMuJKBLUp51s62fthvc5fqqfG4om8ffL\n8omyj59F9oUQ4nJYVJXF8zJxe/yU17aFu5yIIkFNcIvKst/tobPXw50fn8rdn8rDosr/GiGEGInr\nzmzUIXOqR9W4n0e9+2Arv9oU3KLy/tvmcG1+RrhLEkIIU0pPimVmViI1x7o51TXAhCQZ3zMaxu1l\n45ktKv/zleAWlQ+sKJSQFkKIy1RSeOaqWgaVjZZxGdSapvPsh7aonJOTHO6yhBDC9BbMnEBMlGzU\nMZrGXVB7fQF++XI1W8pPMClVtqgUQojRFGWzcPWcDLr7vVQ3dIS7nIgQMqh1XeeRRx5h5cqV3HPP\nPRw/fvyc11999VWWLVtGaWkpzz777LCOCZd+t4+fbqyg/MwWlV+ULSqFEGK0nZlTva1Sur9HQ8jB\nZFu2bMHr9bJx40YqKyspKytj/fr1Q6//+Mc/5rXXXiM6Oppbb72V2267jV27dl30mHBo63bzs+cq\nae0c4Oo56Xz5M7Nl9yshhBgDORlOJqc5qDjUTq/LK+tRXKaQSVVeXk5JSQkAhYWFVFdXn/P6rFmz\n6OnpwePxAMEF2kMdc6U1NPey9un3ae0c4JZrpnC/bFEphBBjJrhRRyYBTWdntWzUcblCplV/fz9O\np3PosdVqRdO0occzZszgzjvv5Pbbb+eGG27A4XCEPOZKqjrczo9/v5e+AR9fWJJH6Q3TUWWbOiGE\nGFPX5mdgtShsq2qSjTrO47VdR4f93pBd3w6HA5fLNfRY0zTU04uB1NbW8uabb/LGG28QGxvLgw8+\nyF/+8hecTucFj7mYtDRnyPeMxOZdjax/cR9WVeHh+67i2nmZo/r1z2e02xAu0g7jiIQ2QGS0IxLa\nAFemHWnANfmZbK9sotPtZ1b26M+sMevn8drOBp5/8zD33J4/rPeHDOri4mK2bt3KzTffTEVFBXl5\neUOvOZ1OYmJisNvtKIpCcnIyfX19FBcX88Ybb5z3mItpaxudXVd0XeflbQ38aWcjjhgb315ewPQM\nx6h9/QtJS3OO+fe4EqQdxhEJbYDIaEcktAGubDuumpXG9som/vTWIVJumT2qX9usn8feujZ++dI+\nnLG2YR8TMqiXLFnCjh07WLlyJQBlZWVs2rQJt9tNaWkpK1asYNWqVdjtdqZMmcIdd9yBxWJh+/bt\n5xxzpfgDGv/7Wg07qltIS4zmgRXzyZDdr4QQ4oqbk5NMSnwU7x48xcqbZhBtH9+LYR462cN/vhpc\nZOu7pYXDPk7RDXTz4HLPjtweP+tf2sf+xi5yM518Z3nhFR1taNYzvA+TdhhHJLQBIqMdkdAGuPLt\neHnbEV7d0ciXPjOLkoKJo/Z1zfZ5NHe4KPvdHgYG/Xx7+TwKpqUOu+s+YoY+d/V5eOz0FpWF01L4\nx7uKZUqAEEKE2XXzMlEY30uK9vR7ePy5SvrdPu65eSYF01JHdHxE9EOcbHfxxHMVdPR6+Pj8iXxR\ndr8SQghDSE2MYXZOEgcau2jucJGZEhfukq4ot8fPE89X0d4zyOeuy+X6wpH3Kpg+zWqPdVG2oZyO\nXg/Lrp/KPZ+eKSEthBAGcqbLe/s4u6r2BzR++XI1R1v7uL4wk88uzrmkr2PqRNt9sJV/+0MFHl+A\n1bfO5rZFOSgyR1oIIQylOC+VuGgrO6pb8AfCs6bGlabrOv/7Wg3VDZ0UTEvh7k/PvOR8MmVQ67rO\n5t3BLSqtluDoucVXYI60EEKIkbNZLVwzN4Nel5d9h8fHRh0vbTvCjuoWcjOdfP1z+ZfV02u6oNY0\nnWdfr+cPbxwi0WHnoS8UMzdXtqgUQggjG9qoYxx0f2/de5JNO48yISmG7ywvJMpuuayvZ6rBZF5f\ngP/ZdIDy2jYmpcbxwIpC2f1KCCFMYEq6k+x0J1WHO+ju95DoiAp3SWNib10bv/trLc5YG99bMTpT\nhE1zRd3v9vHTP1RQXitbVAohhBmVFGai6To79kXmVfWHFzSZkDQ6i22ZIqjbut2s21DOoRM9XDV7\nAg+smE9s9PCXXxNCCBF+18xJx2ZV2V7VHHEbdTR3uHjyhSoCAZ1vLM0nNzN+1L624YO6saWXtRvK\naekc4Jarp/DVz86VLSqFEMKEYqNtLJiZRmuXm/oTPeEuZ9Rc7oImoRg68aoOd/DYM3vpc3mDW1Te\nKFtUCiGEmZ2ZU72tsinMlYyO0VjQJBTDBvXblU08+UIVmq7zjTvmcdOCyeEuSQghxGWaOSWRtMRo\n3qs9hdvjD3c5l2W0FjQJxXBBHdyi8ghPvVZDbLSVf7iriAUz08JdlhBCiFGgKgrXzcvE69N492Br\nuMu5ZKO5oEkohgpqf0Djt/9Xw6s7GklNiOb7dy9g+qSEcJclhBBiFC2el4miwLZK847+/mBBk/jL\nXtAkFMME9cCgjydfqGL7vmZyMpz84J6Fso+0EEJEoOT4aPJzU2ho7uVEW3+4yxmxcxc0KbjsBU1C\nMUxQP7x+x1AXwj+tKiZBtqgUQoiIdWalMrNt1DEWC5qEYpigPnKyh+sLJ/KtO+eN+dmJEEKI8Jo/\nIxVHjI2dJtqoY6wWNAnFMEG9+rP53HuzbFEphBDjgdWisig/g363j4r69nCXE9JYLmgSimFScenH\np8kWlUIIMY6c6f5+u8rYc6rHekGTUAwT1EIIIcaXSWkOpk6MZ/+RTjp7B8NdznmdvaDJ0jFa0CQU\nCWohhBBhU1KQiQ6G3Kjj3AVNJnL7GC1oEooEtRBCiLC5anY6dpvKtqpmNANt1PHRBU3ywnZ7VoJa\nCCFE2MREWfnYrAm09wxSe7Qr3OUMuZILmoQiQS2EECKshjbqMMic6iu9oEkoEtRCCCHCasbkBNKT\nY3m/tg3XoC+stYRjQZNQJKiFEEKElaIolBRk4g9ovHsgfBt1hGtBk1AkqIUQQoTd4vwMVEUJ20Yd\n4VzQJBQJaiGEEGGX4IiiYFoKR1v7ONbad0W/99kLmtwbhgVNQpGgFkIIYQglhcGVyq7kVfWHFzQp\nCcOCJqFIUAshhDCEeVNTiI+zs+tACz5/YMy/n1EWNAlFgloIIYQhWC0qi/MzcA36Ka9rG9PvZaQF\nTUKxhnqDruusWbOG2tpa7HY7a9euJSsrC4D29nYeeOABFEVB13Vqamp48MEH+fznP8+yZctwOBwA\nTJ48mXXr1o1tS4QQQpjedQWZvPbuMbZXNXPNnIwx+z5GWtAklJBBvWXLFrxeLxs3bqSyspKysjLW\nr18PQGpqKhs2bACgoqKCJ554ghUrVuD1egF4+umnx7B0IYQQkSYzJY4ZkxM40NhFe7eb1MSYUf8e\nRlvQJJSQpxDl5eWUlJQAUFhYSHV19Xnf96Mf/YhHH30URVGoqalhYGCA1atXc99991FZWTm6VQsh\nhIhYZ1Yq2z4GG3UYcUGTUEJeUff39+N0Oj84wGpF0zTUs7oJ3njjDfLy8sjOzgYgOjqa1atXU1pa\nSmNjI/fffz+bN28+55jzSUtzXvR1M4iENoC0w0gioQ0QGe2IhDaA8dtxy3UxPPt6HTv3t/LlpQVY\n1PPfOx5pO2oaO/mvPx3AbrOw5v5ryZuSNBrljrmQQe1wOHC5XEOPPxzSAK+++ir33nvv0OOcnJyh\n0M7JySExMZG2tjbS09Mv+r3a2q7s3LnRlpbmNH0bQNphJJHQBoiMdkRCG8A87fjYrAm8XdnM2+8d\nJX9qykdeH2k7mjtclP1uD36/xreXzyMpxhr2/w/DPdEI2fVdXFzMW2+9BQTvQ+fl5X3kPdXV1RQV\nFQ09fvHFF/nXf/1XAFpbW3G5XKSlpQ2rICGEEOJM9/fbo7BRh9EXNAkl5BX1kiVL2LFjBytXrgSg\nrKyMTZs24Xa7KS0tpbOz85yucYDly5fz8MMPs2rVKlRVZd26dSG7vYUQQogzpk6MZ2JqHHvr2ugb\n8OKMvbR7yWZY0CSUkEGtKAqPPvroOc/l5uYO/XdycjIvvfTSOa/bbDZ++tOfjlKJQgghxpszG3X8\n4Y1D7NrfypKPZY34a5hlQZNQ5DJXCCGEIV2bn4FFVdhW1YSu6yM61kwLmoQiQS2EEMKQ4mPtzJ+R\nyok2F40tIxv4ZaYFTUIxb+VCCCEiXknBmY06moZ9zDkLmpQaf0GTUCSohRBCGFZ+bgpJzijePdiK\nxxd6o46PLGhyiYPQjESCWgghhGGpqsLieRm4PQHKa09d9L2HTvbwX6/ux2ZV+W5pIROSYq9QlWNL\ngloIIYShXXd6TvXF9qlu7nDx5AtV+AM631iaT25m/JUqb8xJUAshhDC0CYkxzJqSSO3xblq7Bj7y\nutkXNAlFgloIIYThDW3U8aGVyiJhQZNQJKiFEEIY3oKZacREWdm+r5mApgGRs6BJKBLUQgghDM9u\ns3DNnHR6+r3sO9J5zoImhSZf0CSUkEuICiGEEEZQUpjJ1r0n2V7VTEv34NCCJn9n8gVNQpGgFkII\nYQrZ6U6yJjjYW9fGnrq2iFnQJJTIPQURQggRUc5s1KEDCQ57xCxoEopcUQshhDCNkoKJnOpy89kb\npuOwjY9rzfHRSiGEEBEhym5h1ZI8cicmhLuUK0aCWgghhDAwCWohhBDCwCSohRBCCAOToBZCCCEM\nTIJaCCGEMDAJaiGEEMLAJKiFEEIIA5OgFkIIIQxMgloIIYQwMAlqIYQQwsAkqIUQQggDk6AWQggh\nDEyCWgghhDAwCWohhBDCwCSohRBCCAOzhnqDruusWbOG2tpa7HY7a9euJSsrC4D29nYeeOABFEVB\n13Vqamp48MEHWbFixQWPEUIIIcTwhQzqLVu24PV62bhxI5WVlZSVlbF+/XoAUlNT2bBhAwAVFRU8\n8cQTrFix4qLHCCGEEGL4QgZ1eXk5JSUlABQWFlJdXX3e9/3oRz/iZz/7GYqiDPsYIYQQQlxcyHvU\n/f39OJ3OocdWqxVN0855zxtvvEFeXh7Z2dnDPkYIIYQQoYW8onY4HLhcrqHHmqahqufm+6uvvsq9\n9947omPOJy3NGfI9RhcJbQBph5FEQhsgMtoRCW0AaYfZhAzq4uJitm7dys0330xFRQV5eXkfeU91\ndTVFRUUjOuZ82tr6RlC68aSlOU3fBpB2GEkktAEiox2R0AaQdhjJcE80Qgb1kiVL2LFjBytXrgSg\nrKyMTZs24Xa7KS0tpbOz85xu7gsdI4QQQoiRU3Rd18NdxBmRcHZk9jaAtMNIIqENEBntiIQ2gLTD\nSIZ7RS0LngghhBAGJkEthBBCGJgEtRBCCGFgEtRCCCGEgUlQCyGEEAYmQS2EEEIYmAS1EEIIYWAS\n1EIIIYSBSVALIYQQBiZBLYQQQhiYBLUQQghhYBLUQgghhIFJUAshhBAGJkEthBBCGJgEtRBCCGFg\nEtRCCCGEgUlQCyGEEAYmQS2EEEIYmAS1EEIIYWAS1EIIIYSBSVALIYQQBiZBLYQQQhiYBLUQQghh\nYBLUQgghhIFJUAshhBAGJkEthBBCGJgEtRBCCGFgEtRCCCGEgUlQCyGEEAYmQS2EEEIYmDXUG3Rd\nZ82aNdTW1mK321m7di1ZWVlDr1dVVfHYY48BkJqayk9+8hPsdjvLli3D4XAAMHnyZNatWzdGTRBC\nCCEiV8ig3rJlC16vl40bN1JZWUlZWRnr168fev2HP/whv/jFL8jKyuKFF16gqamJiRMnAvD000+P\nXeVCCCHEOBCy67u8vJySkhIACgsLqa6uHnqtoaGBxMREfvvb33L33XfT09NDTk4ONTU1DAwMsHr1\nau677z4qKyvHrgVCCCFEBAt5Rd3f34/T6fzgAKsVTdNQVZWuri4qKip45JFHyMrK4mtf+xr5+fkk\nJSWxevVqSktLaWxs5P7772fz5s2oqtwSF0IIIUYiZFA7HA5cLtfQ4zMhDZCYmMiUKVPIzc0FoKSk\nhOrqau655x6ys7MByMnJITExkba2NtLT0y/6vdLSnBd93QwioQ0g7TCSSGgDREY7IqENIO0wm5BB\nXVxczNatW7n55pupqKggLy9v6LWsrCwGBgY4fvw4WVlZlJeXs3z5cl544QXq6up45JFHaG1txeVy\nkZaWFrKYtra+y2tNmKWlOU3fBpB2GEkktAEiox2R0AaQdhjJcE80Qgb1kiVL2LFjBytXrgSgrKyM\nTZs24Xa7KS0tZe3atXzve98DoKioiI9//OP4fD4efvhhVq1ahaqqrFu3Trq9hRBCiEug6Lquh7uI\nMyLh7MjsbQBph5FEQhsgMtoRCW0AaYeRDPeKWi5zhRBCCAOToBZCCCEMTIJaCCGEMDAJaiGEEMLA\nJKiFEEIIA5OgFkIIIQxMgloIIYQwMAlqIYQQwsAkqIUQQggDk6AWQgghDEyCWgghhDAwCWohhBDC\nwCSohRBCCAOToBZCCCEMTIJaCCGEMDAJaiGEEMLAJKiFEEIIA5OgFkIIIQxMgloIIYQwMAlqIYQQ\nwsAkqIUQQggDk6AWQgghDMwa7gIigdrchKXhCFgC2AIWArlT0TInhrssIYQQEUCC+jKpzU3Y3tyK\n2toCFg1rQEU9ehTfDTdKWAshhLhsxun63rQJ287tqM1N4a5kRKzv78bS2IDidoOuo7jdWBobsL6/\nO9ylCSGEiADGCWpdR+nrw1pVaaqwthw5PKLnhRBCiJEwZNe3peGIqbqNlf4+lO4uUHVUTUFPTEKP\njQl3WcLEZNyDGG3yM2Vehgxqpb8v3CUMm5aUhHXfvuCDaCuKx4fS2op/5szwFnYJIuUX2VqxF+vu\nXeDpJzrKgf+qa/DPLwp3WcOmNjdhraoMPkiIGepp8oMpPw8RfpH0MxUpf6dGwpBBrTuc4S5h2PSk\nZLT0CSjd3aDo6FF29MRE9KTkcJc2IpHyi2yt2Iv9tT8HH8TZUTs6hh6bJawtDUdQOjuHBihaAipa\neobpepqEcVgajlzweTP9TEXK36mRMmRQB3KnhruE4YuKQps8BdXtBt8g+unHREWFu7IRiZRwsO7e\nhdLfGzxxGroVkYh19y7TBLV67Ci2rVuwNh4BzyBRUdH4c6bi+8QnYdF14S5vRCLh6icS2nChXkoz\n9V5C5JxwjJRxgvqhh4h2JuC5/XNon74l3NUMn8eDZW958AfI58FqiyLQ309g8uRwVzYi6rGj2N7d\nidp0EvxebFY72sRJ+BRMFQ6WY0ex1NagdnSC5sOq2tBSkkExzrjJUKzv78ZWuRd8PkBH6evH1tOD\nHh+PZ+UXwl3esKnNTUS9/EfU+lrwDWK3RaPNmIln6TLT/FGNlOmXusOJ0vfRUDZT7yVEzgnHSIUM\nal3XWbNmDbW1tdjtdtauXUtWVtbQ61VVVTz22GMApKam8pOf/ASbzXbRY86rtxdLv4uY3/8OPTEJ\n74q7Lq9lV4ilqgJr+XvBK2pdw6KoKJ0dBLKz4fal4S5v2KwHD2Ct2Bv8gdcCWFQL6qlT6AmJeMJd\n3Agorc2ozc3BB1YVvB7U5ma0pKTwFjYClpqDMHj6/7pFgYAGAU/weROx/20z1vL3gw+irai9fajl\n76PHxTF4z5fCW9wwnZl+CUCcfWj6pf5+HF4T/X4Hcqd+0GX8oefNJFJOOEYqZFBv2bIFr9fLxo0b\nqayspKysjPXr1w+9/sMf/pBf/OIXZGVl8cILL9DU1ER9ff1FjzkvtxtFV4Jdfc9uME1Q297ZieIe\nAI8PFB10BQUd2zs7cYe7uBFQaw6gnmoBrzc4VU5RUAZcqDUHwl3ayHi9KB43ysAg6AFUxYIeGx1s\nl6PBAzgAABBDSURBVEkobhd6dBSK1wcKoCrodhuK2xXu0kbEUrH3ws+bJKgtRw6jtjYHe5oCPqwW\nG9rESVhMNqtDy5yIn9Ndx0oA3ek0ZRd+IHcqtjffQG1tPesWXTr+gsJwlzamQgZ1eXk5JSUlABQW\nFlJdXT30WkNDA4mJifz2t7+lvr6eG264gZycHDZu3HjBYy5oYAAAxWfDeuYM1gQsLU0ofj+ggQ6g\no/j9WFrMMxccQD1xAjwe8Ps/eFLXg8+biOLzoQ96wO0CTUNXVXRVQfH5wl3asGmJKVj7XcGQ1nVQ\nFBQgkJgS7tJGRPEOjuh5I1KbmrAcOhR8YLegDLixHDqEHhUd3sLGNSXEY+M7e9wDt90W8v0hg7q/\nvx+n84NuBavViqZpqKpKV1cXFRUVPPLII2RlZfG1r32NuXPnXvSYCwoEsABoGviiSUszSVeGogQD\n7nQQWABsNlAU87QBwOMO/r8/7cxnYfG4zdUOVx8E/BAd/ENqBQj4sbn6iDZLO667Fp4/AQEFNB2r\nqoBFxXrdteZpA8CcWVBeDn190Okn2moFpxPmFBBjlnboPrBbhh7aT/+3XfcRa5Y2AJw8CUfrgktc\n6ZCo+oOPUxwwaVK4qxu+Ay2QOyn4D+A483xXCxSYZErsyZOw5x1obga3e3SC2uFw4HJ90OV2duAm\nJiYyZcoUcnNzASgp+f/bu9vYKKq+j+Pf3W53W9rSAhaFG2iLppeACU+KICJoQBHUgG0NggpJbxNB\nEqJosEgCSHg2+EAgFBEJWCkBSyC+EAIIaNUbJBQBg0psFUpFBQQKbfdh5n6xdS+8FLat4ZoHfp9X\ny+5M+J90d34z55w5M4gjR46QlpZ21X2uKhwmAuD1YvgD/P6rMyYHZHh9eEOhaKhBtA2hEIbX55g2\nAKQnJJJgAmbjsGjj60hCIued1I5gCK/PhycUIcFjEjE9mIkJGMGQY9rhv60b/jt6klD1A/5gA0F/\ngEh2V4K3dSPokDYAJObeQVL5/+EJRgh4oSEYwawPU597ByGHtCPF34oErw/v2bMkRkKEEhIx2rYl\n4m/FJYe0ASDxq6/xXIwOxqWnJ3P+fPS1+dXXhPytrSytWfwnTzf2XP65HVyoc8xvw79tF74j0fkm\nKSn+Ju0TN6j79OnDJ598wvDhw6moqCA3Nzf2WefOnbl8+TInTpygc+fOHDhwgPz8fLp06XLVfeJK\nSMBMTY2/nU2YXi94vdEuysZuSrze6PsOYma2h/Pnor0DhgEJXggEou87iJGWjuemerhYC2YE05OA\nmZaKkZZudWlNFwgQun8okdM/408wCDbeKue0W/5o1YrQvfeR8N0xCNcT8SURyb0dWrWyurImMzIy\nSEhKxkhPB4+JYXogKRkjI8Pq0prFLbOl3TCZrCXLS8cN6mHDhlFeXs6YMWMAmD9/Ph999BF1dXUU\nFBQwd+5cXnzxRQB69+7N4MGDMU3zL/vEr8QXDbnkZEh0zgHJY5rRA2g4QvRUzwO+hOj7DmLk5GDW\nVENCQmxsl6RkjMbeEqcwunXHe/kyZkoa+BMwgpHY+05hpqZFezPatoX0ZCJ/XP046GAE0RAwunbF\n6NoV0pMJNbbDSeFgtG+PmZQMScmQ5MOsD8fedxI3BBy4Z/Z6c8UNao/Hw+zZs//0Xs4VB++7776b\njRs3xt0nrkAg1qVhJjtnooaZkop56SKeUCQ269tMTMBMcU6vAEAkKxtvTQ2es7+BEcHwJmC2vYlI\nVrbVpTVLw+h8POd/j90PbrZKxuj4PzSMzre6tCZzy8HIDeFg3tKBUN+7Yr0CRloakdzbMW/pYHVp\nzeKW75TRoSPG6dOxJYK9jUsEO2n2eqTrrfiaMsH6CvZZ8CQQwMSDmZzsqK7vSPfueGsvYDYuTmHi\ngcREIt2dcwUHYHTqTORf52IPF4k0PlzE6BTn/nebCffqTf2E/8W370sCDbUEHbjWt5tupXF6OJip\naX/bK+Ckkw1wz3fKW3MK7+mfMbKyIT0Z43wd3tM/46055Zi2hO/sh6f2UnQRHY8RfwfsFNSZmdGQ\nCyRhdnDO2WrwweF4GhrwVv0AwQYMfwAjuyvBB4dbXVqzGF2yCJvEVmCKNI6LGl2yrC6t2cK9ehPu\n1Zu0zDTqHTLB5D8ZHTpGDzyZaY6ZePWf3BAObjjZ+IMbvlNuWELU6NCR0JD7/317VhPYJ6hvugmj\ncTGBcLc7rK6myaJnR7V4T99xxcSfmwnf2c/q0polktMVz8WLfxkXdeIBSezD6eHghpMNN3HLpLgr\nfxdNYZ+gHjyY8KXo6lHhfv0tLqbpomdHD8TOjsIOXbTfDWM/IteD00823MQN8x5awj5B7fFgtGvn\nuPFEcMcP2Q1jPyLibm4aimgO+wT1a685djzRDdww9iMi7najDkXYJ6jFUm4Z+xERd3NDD2ZzOWv5\nLLlurjbG4/axHxERu1NQC3D1MR63j/2IiNidur4FuHHHfkRE7E5BLTE34tiPiIjdqetbRETExhTU\nIiIiNqagFhERsTEFtYiIiI0pqEVERGxMQS0iImJjCmoREREbU1CLiIjYmIJaRETExhTUIiIiNqag\nFhERsTEFtYiIiI0pqEVERGxMQS0iImJjCmoREREbU1CLiIjYmIJaRETExhTUIiIiNqagFhERsTFf\nvA1M02TWrFl8++23+P1+5s6dS+fOnWOfr1mzhk2bNtG2bVsAXnvtNbKzs3n88cdJTU0FoFOnTsyb\nN+86NUFERMS94gb1jh07CAaDlJaWcujQIebPn8/y5ctjnx89epRFixbRvXv32HvBYBCAtWvXXoeS\nRUREbhxxu74PHDjAoEGDAOjZsydHjhz50+dHjx6luLiYsWPHsnLlSgCOHTvG5cuXKSwsZMKECRw6\ndOg6lC4iIuJ+ca+oa2trSUtL+/cOPh+GYeD1RjN+5MiRjBs3jtTUVJ5//nn27NlDx44dKSwspKCg\ngKqqKp599lm2bdsW20dERESaJm5Qp6amcunSpdi/rwxpgPHjx8fGogcPHsw333zDPffcQ5cuXQDI\nzs4mIyODX3/9lZtvvvma/1dmZto1P3cCN7QB1A47cUMbwB3tcEMbQO1wmriXuH369GHPnj0AVFRU\nkJubG/ustraWRx55hLq6OkzT5Msvv6RHjx5s2rSJBQsWAHD69GkuXbpEZmbmdWqCiIiIe3lM0zSv\ntcGVs74B5s+fz9GjR6mrq6OgoICtW7eydu1aAoEAAwYMYPLkyYRCIYqKijh16hRer5eXXnqJXr16\n/VcaJCIi4iZxg1pERESso9ldIiIiNqagFhERsTEFtYiIiI0pqEVERGws7n3U11O8dcSd5tChQ7z+\n+uusW7fO6lJaJBwOM336dKqrqwmFQjz33HM88MADVpfVLIZhMGPGDCorK/F6vcyePZvbbrvN6rJa\n7MyZM+Tl5fHee++Rk5NjdTnN5pY1/1euXMmuXbsIhUKMHTuWvLw8q0tqts2bN1NWVobH46GhoYFj\nx45RXl4e+/s4QTgcZtq0aVRXV+Pz+ZgzZ44jfxfBYJCioiJOnjxJamoqM2fOjK098ncsDep464g7\nyapVq9iyZQspKSlWl9JiW7dupU2bNixatIjz588zatQoxwX1rl278Hg8rF+/nn379rFkyRLHfqfC\n4TAzZ84kKSnJ6lJaxC1r/u/bt4+DBw9SWlrK5cuXWb16tdUltcjo0aMZPXo0EH14Un5+vqNCGmDP\nnj0YhkFpaSmff/45b7zxBm+//bbVZTXbxo0bSUlJYcOGDVRWVjJ79mzefffdq25vadd3vHXEnSQr\nK4tly5ZZXcY/8vDDDzNlyhQgemXq81l6HtciQ4cOZc6cOQBUV1eTnp5ucUUtt3DhQp588knat29v\ndSkt4pY1/z/77DNyc3OZNGkSEydO5P7777e6pH/k8OHDHD9+nIKCAqtLabbs7GwikQimaXLx4kUS\nExOtLqlFjh8/zn333QdATk4OP/zwwzW3t/RIHG8dcScZNmwY1dXVVpfxjyQnJwPRv8uUKVN44YUX\nLK6oZbxeL6+88go7duxw5Nk2QFlZGe3atWPgwIGsWLHC6nJaJCkpyRVr/p87d45Tp05RXFzMiRMn\nmDhxIh9//LHVZbXYypUrmTx5stVltEhKSgonT55k+PDh/P777xQXF1tdUot069aN3bt3M3ToUCoq\nKvjll18wTROPx/O321v6i4m3jrj899XU1DB+/HhGjx7NiBEjrC6nxRYsWMC2bduYMWMG9fX1VpfT\nbGVlZZSXl/P0009z7Ngxpk2bxpkzZ6wuq1mys7N57LHHYq//WPPfaTIyMhg0aBA+n4+cnBwCgQBn\nz561uqwWuXjxIlVVVfTr18/qUlpkzZo1DBo0iG3btrF161amTZsWG2Jxkry8PFJSUhg3bhw7d+6k\nR48eVw1psDior7WOuFM5eaG33377jcLCQl5++eXYWJbTbNmyJfa41UAggNfrdeTJ3/vvv8+6detY\nt24dt99+OwsXLqRdu3ZWl9UsH374oSvW/O/bty+ffvopEG1HfX09bdq0sbiqltm/fz/9+/e3uowW\nS09Pj42rp6WlEQ6HMQzD4qqa7/DhwwwYMICSkhIeeuihuJOoLe36HjZsGOXl5YwZMwaIriPudNc6\nK7K74uJiLly4wPLly1m2bBkej4dVq1bh9/utLq3JHnzwQYqKinjqqacIh8O8+uqrjqr/7zj1O5Wf\nn09RURFjx47F6/Uyb948R540DRkyhK+++or8/HxM02TmzJmO/ZtUVlY6+s6a8ePHM336dMaNG0c4\nHGbq1KmOnGyZlZXFW2+9xYoVK2jdujVz58695vZa61tERMTGnHd6KyIicgNRUIuIiNiYglpERMTG\nFNQiIiI2pqAWERGxMQW1iIiIjSmoRUREbExBLSIiYmPOezySiDRZJBJh1qxZfP/995w5c4acnByW\nLl3Khg0bKCkpoXXr1uTk5NClSxcmT57M3r17Wbp0KZFIhE6dOjFnzhxHP4FMxA10RS3iYgcPHsTv\n91NaWsr27dupq6vjnXfeYf369WzevJmSkhJ+/PFHAM6ePcuSJUtYvXo1ZWVlDBw4kMWLF1vcAhHR\nFbWIi915551kZGRQUlJCZWUlP/30E/3792fIkCG0atUKgJEjR3LhwgW+/vprampqeOaZZzBNE8Mw\nyMjIsLgFIqKgFnGxnTt3snTpUiZMmEBeXh7nzp2jdevWXLhw4S/bRiIR+vbty/LlywEIBoN/egyt\niFhDXd8iLvbFF18wYsQIRo0aRdu2bdm/fz+mabJ3715qa2sJBoNs374dj8dDz549qaiooKqqCoBl\ny5axaNEiaxsgInp6loibfffdd0ydOpXExET8fj/t27fn1ltvJTMzkw8++ICUlBTatGnDXXfdRWFh\nIbt37+bNN9/EMAxuueUWFi9erMlkIhZTUIvcYKqqqti9ezcTJkwAYNKkSTzxxBMMGTLE0rpE5O9p\njFrkBtOxY0cOHz7Mo48+isfj4d5771VIi9iYrqhFRERsTJPJREREbExBLSIiYmMKahERERtTUIuI\niNiYglpERMTG/h8CYXNPtVPzfgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11538e128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = (kawasaki_subset[kawasaki_subset.age_years<10]\n",
    "              .assign(age_int=kawasaki_subset.age_years.astype(int))\n",
    "              .groupby('age_int')).responder.mean().plot()\n",
    "ax.set_xlabel('age')\n",
    "responders = kawasaki_subset[kawasaki_subset.responder==1]\n",
    "ax.plot(responders.age_years, np.random.normal(0.9, 0.005, size=len(responders)), \n",
    "        'go', alpha=0.3)\n",
    "\n",
    "nonresponders = kawasaki_subset[kawasaki_subset.responder==0]\n",
    "ax.plot(nonresponders.age_years, np.random.normal(0.6, 0.005, size=len(nonresponders)), \n",
    "        'ro', alpha=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "StudyID        0\n",
       "male           0\n",
       "age_norm       0\n",
       "age_norm2      0\n",
       "age_years      0\n",
       "ethnicity      0\n",
       "KD_complete    0\n",
       "Tier           0\n",
       "black          0\n",
       "other          0\n",
       "responder      0\n",
       "dtype: int64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kawasaki_subset.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def interpolate(x0,y0, x):\n",
    "    x = np.array(x)\n",
    "\n",
    "    idx = np.searchsorted(x0, x)\n",
    "    dl = np.array(x - x0[idx - 1])\n",
    "    dr = np.array(x0[idx] - x)\n",
    "    d=dl+dr\n",
    "    wl = dr/d\n",
    "\n",
    "    return wl*y0[idx-1] + (1-wl)*y0[idx]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO (theano.gof.compilelock): Waiting for existing lock by process '48770' (I am process '50575')\n",
      "INFO (theano.gof.compilelock): To manually release the lock, delete /Users/fonnescj/.theano/compiledir_Darwin-15.5.0-x86_64-i386-64bit-i386-3.5.1-64/lock_dir\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Applied log-transform to σ and added transformed σ_log_ to model.\n"
     ]
    }
   ],
   "source": [
    "import theano.tensor as tt\n",
    "GaussianRandomWalk = pm.distributions.timeseries.GaussianRandomWalk\n",
    "\n",
    "covs = ['male', 'black', 'other']\n",
    "k = len(covs)\n",
    "nknots = 7\n",
    "\n",
    "knots = np.linspace(kawasaki_subset.age_norm.min(), kawasaki_subset.age_norm.max(), nknots)\n",
    "\n",
    "with pm.Model() as model:\n",
    "    \n",
    "    # Baseline probability of response\n",
    "    θ = pm.Normal('θ', 0, sd=10)\n",
    "    # Coefficients for covariates\n",
    "    β = pm.Normal('β', 0, sd=10, shape=k)\n",
    "    \n",
    "    odds = pm.Deterministic('odds', tt.exp(β))\n",
    "    \n",
    "    # Age effect (non-linear)\n",
    "    σ = pm.HalfCauchy('σ', 2.5)\n",
    "    y = GaussianRandomWalk('y', sd=σ, shape=nknots)\n",
    "\n",
    "    α = pm.Deterministic('α', interpolate(knots, y, kawasaki_subset.age_norm))\n",
    "    \n",
    "    \n",
    "    # Probabilities\n",
    "    p_race_female = pm.Deterministic('p_race_female', pm.invlogit(θ))\n",
    "    p_race_male = pm.Deterministic('p_race_male', pm.invlogit(θ + β[1]))\n",
    "    p_race_black_female = pm.Deterministic('p_race_black_female', \n",
    "                                           pm.invlogit(θ + β[1]))\n",
    "    p_race_black_male = pm.Deterministic('p_race_black_male', \n",
    "                                         pm.invlogit(θ + β[0] + β[1]))\n",
    "        \n",
    "    π = pm.invlogit(θ + tt.dot(kawasaki_subset[covs].values, β)+ α*kawasaki_subset.age_norm.values)\n",
    "    \n",
    "    responder = pm.Bernoulli('responder', π, observed=kawasaki_subset.responder)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0 [0%]: ELBO = -277.97\n",
      "Iteration 10000 [10%]: Average ELBO = -267.98\n",
      "Iteration 20000 [20%]: Average ELBO = -224.99\n",
      "Iteration 30000 [30%]: Average ELBO = -222.47\n",
      "Iteration 40000 [40%]: Average ELBO = -222.0\n",
      "Iteration 50000 [50%]: Average ELBO = -221.95\n",
      "Iteration 60000 [60%]: Average ELBO = -221.89\n",
      "Iteration 70000 [70%]: Average ELBO = -221.98\n",
      "Iteration 80000 [80%]: Average ELBO = -221.93\n",
      "Iteration 90000 [90%]: Average ELBO = -221.92\n",
      "Finished [100%]: Average ELBO = -221.91\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    \n",
    "    advi_fit = pm.variational.advi(n=100000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "iterations = 2000\n",
    "burn = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Assigned NUTS to θ\n",
      "Assigned NUTS to β\n",
      "Assigned NUTS to σ_log_\n",
      "Assigned NUTS to y\n",
      " [-----------------100%-----------------] 2000 of 2000 complete in 57.7 sec"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    \n",
    "    trace = pm.sample(iterations, start=advi_fit[0], njobs=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Odds ratios for responder. Tabulated values below plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.gridspec.GridSpec at 0x11b8de5c0>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAF7CAYAAADiw5DPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtYVHXix/HPDCM+oKaoaO0+pmkqKksZ7Fbm5RfpbuTd\nNNEEU0t7yvLSZrmSlsF6q9w03dWsdK3ULDWitlwv3XS94NOKWoarj5cnrTBxFVAR5/v7w3WSBBlB\n+M6B9+ufcIaZ8+XMOcybM+f0dRljjAAAACxy2x4AAAAAQQIAAKwjSAAAgHUECQAAsI4gAQAA1hEk\nAADAOoIEZbZkyRL17NlT3bp1U/fu3fXUU0/pyJEjxX5/27Ztdfjw4Utu/+STT5SQkFAuY9yxY4cm\nTZokSdq5c6dGjRpVLsu52MqVK3XnnXfqwQcfLPdl+aN3797Kycm57Pd89tlnmjVrVgWNCKi8IiIi\n1KNHD/Xq1Uu9e/fW3XffrX79+mnnzp1Ffv8rr7yi5OTkK17OsGHDdPz48bIONyB4bA8AzjZt2jRl\nZmZq/vz5atiwoSRp1apV6t+/v5YvX+677WIul6vY57vcfWWxZ88e/fDDD5KkyMhIvfzyy+WynIut\nWrVKY8eOVffu3ct9Wf5YuXJlid+zY8cOnThxogJGA1RuLpdLixcvVu3atX23vf7660pOTtbSpUuv\n2nI2bNhw1Z7LNoIEpfbDDz9o6dKl+uKLL1SzZk3f7b169dLXX3+t+fPn65lnnlF6erqSk5PldrsV\nGRmpi/9ffC+//LLS0tIUFham66+/3nd7enq6pk2bJq/XK5fLpREjRqhLly6Flr9lyxalpKQoJCRE\np0+f1jvvvKPp06drx44dys3NlTFGycnJuu666zR79mzl5OToT3/6k3r16qXnn39eH3zwgXJycvTc\nc89p9+7dcrlc6tChg5544gm53W7NmjVLa9euVbVq1VSnTh1NnTpV9evXLzSGXz6+Y8eOGjNmjKZP\nn66MjAx99913OnbsmAYPHlzoce+++64WLlyooKAghYWFadq0aWrYsKGWLVumN998U0FBQapXr54m\nTpyoevXqqVOnTlq9erXq1asnSerfv79GjhypRo0aafLkycrLy9OPP/6oVq1aaebMmQoODtZvfvMb\n3XXXXfr22281Y8YM9e3bV5s2bVL16tX17LPP6sCBAzp+/Lhq1KihF198USdOnNDSpUvl9XpVs2ZN\njR49WsuXL9eSJUskSXXq1FFSUpKaNm3q1+sDVGXGmEK/686dO6fDhw+rTp06xT5m7969SkxMVFZW\nlurXr6+ZM2eqfv36Wr9+vebNm6eCggIdO3ZMvXr10uOPP67x48dLkhITE/Xqq68W+QegoxiglD75\n5BPTt2/fIu9bt26d6dmzp8nPzzd33HGH2bRpkzHGmLS0NBMREWG+++47889//tN069bN5OXlmXPn\nzpkRI0aYhIQEY4wxgwcPNh9++KExxpjdu3ebyZMnX7KMzZs3m9atW5sjR44YY4z56quvzKhRo3z3\nz5s3zzz88MPGGGNWrFhhRowY4Xtct27djDHGjBs3zqSkpBhjjMnPzzdDhw418+fPN0eOHDHR0dEm\nPz/fGGPMG2+8YdasWXPJGJ566qkiH2+MMYMGDTKffPLJJY/55ptvzG233Wa+//57Y4wxixYtMpMm\nTTL/+te/zO9//3uTnZ3tG/M999xjjDHm6aefNq+//roxxpj//Oc/5s477zTGGDNt2jSTmppqjDHm\n7Nmzpnv37mb16tXGGGNatmzpu88YYyIiIkx2drb5+OOPTXJysu/2iRMnmueff94YY8zs2bN9X2/Z\nssXcf//95vTp08YYY7788kvfePx5fYCqrGXLlqZ79+6mR48epn379uauu+4yycnJ5qeffiry+2fP\nnm06d+7s2/8feeQRM3fuXGOMMYmJiebAgQPGGGN++OEH07p1a9/3tWzZ0hw/frwCfqLyxxESlElB\nQUGRt+fn58vlcikzM1PVqlXTrbfeKknq2rWr71yOTZs2qUuXLgoJCZEk3XvvvVq8eLEkKS4uTpMn\nT9a6devUrl07jRkzpsjlXHvttbr22mslSTfffLNGjRqlJUuW6ODBg9qyZUuhIzdF+eKLL3yHT6tV\nq6YBAwZo0aJFeuihh9SqVSv17t1bHTp0UMeOHXX77bdf8vjPP/+82McXZ9OmTerQoYPvr5nExERJ\n0owZMxQXF+f7C6p3795KSUnRd999p759++q5557TkCFDtGLFCvXp00eS9OSTT2rDhg1asGCB9u/f\nr6ysLOXm5vqWFR0d7fva/O+vtT/84Q9q1KiR3nzzTR04cEBbtmxR27ZtLxnnp59+qoMHDyo+Pt73\n2BMnTujEiRN+vz5AVXbhI5tvvvlGDz30kNq2bau6desW+/3t2rXz7f8RERH66aefJEl//etf9emn\nnyo1NVX79u2TJJ06dcr3vaaSzADDSa0otZtuukn79+/37TQX27x5s9q2bSuXyyWv11vovqCgIN/X\nF+9IF9/ev39/ffDBB2rfvr2+/PJL9ejRo8gTMkNDQ31ff/rppxoxYoRcLpc6d+5c6I20OL8cm9fr\n9UXW4sWLNXXqVIWFhWnKlClKSUm55PG/fP6LH1+coKCgQufKnDlzRvv27btkLBeev6CgQNHR0Tp3\n7pwyMjKUlpamfv36SZLGjBmjd955R7/+9a81ZMgQtW7dutDjL14/F5b59ttva8KECQoJCVH37t3V\ntWvXIteT1+tVz549tXLlSq1atUqrVq3Su+++q2uuucbv1weoyi7sV61atdL48eM1YcIE3wn9w4cP\n953wun79eknn/6i54ML+eurUKd/H4JGRkRo3bpyCgoIqTYRcjCBBqTVs2FCJiYkaO3as74RRSXrv\nvfe0evVqDR8+XC1atJB0/kiCJK1du9Z30mSHDh308ccf6+TJk/J6vXr//fd9zxEfH6+vv/5avXr1\n0uTJk3Xy5MkST7bcuHGjYmNjFR8fr8jISK1du9b3Jh8UFFRkKLRv315vvfWWpPNHdZYtW6Y77rhD\nu3fvVrdu3dSsWTMNHz5cDzzwgL799lu/H385t956qzZu3KijR49KOn+V0gsvvKCOHTvqH//4h44d\nO+Zbj2FhYWrcuLEkqW/fvkpOTlZERITv6MrGjRv16KOPKi4uTsYYbd++XefOnStyuRd+gW3YsEF9\n+vTRvffeqyZNmmj9+vWF1tPZs2clSXfccYc+/PBDZWVlSZLeeustPfDAA5JK9/oAVVnXrl11yy23\n+P6wmT9/vlatWuW7Gq84Bw4cUF5enkaPHq3/+7//0+bNm3X27Fnffu7xeEr8I8gp+MgGZTJmzBi9\n9957euSRR5Sfn6/8/HxFRUVp2bJlvo9S5syZo4kTJ2rmzJmKiIjwnZjZqVMn7dmzR/fee69q166t\niIgIZWdnSzr/UURKSopefvlluVwujRw5Ur/61a8uO5b4+Hj98Y9/VM+ePRUUFKSYmBitXr1a0vlL\njf/yl7/oscceK3RpcVJSkp5//nl1795dZ8+eVceOHfXwww/L4/EoLi5Offr0UWhoqEJCQpSUlHTJ\nMidMmFDk46Xirxhq0aKFxo0bp2HDhsnlcik8PFx//vOfFR4ersGDB/tOgA0LC9O8efN8j+vVq5dm\nzpypl156qdD6f/TRR1WnTh2FhITod7/7nQ4ePFjk8i/8e+jQoZo4caJWrFght9utNm3aKDMzU5J0\n++2367HHHlO1atWUlJSkBx98UEOHDpXb7VbNmjX1yiuvSJLGjRun5OTkK3p9gKqkqP0/KSlJPXv2\n1IYNG0r8w+WCiIgIderUSXfffbeuueYaNW7cWDfeeKMOHjyoRo0aqXPnzho4cKDmzp2rG2+88Wr/\nGBXKZSrjcR8AAOAofGQDAACsI0gAAIB1BAkAALAuYE9qzco6aXsIASEsLFTZ2Xm2hxHQoqMj5Xa7\ntHXrDttDCWhsS/6p7OspPLyW7SGUir/vCU5+/Zw8dsn/8Re3DXKEJMB5PEElf1MVt23bTu3fv9/2\nMAIe25J/WE/O5uTXz8ljl8o+foIEAABYR5AAAADrCBIAAGAdQQIAAKwjSAAAgHUECRwvOjpSTZo0\nsT0MAEAZECQAAMA6ggQAAFhHkAAAAOsIEgAAYB1BAgAArCNI4HjMZQMAzkeQAAAA6wgSAABgHUEC\nAACsI0gAAIB1BAkAALCOIIHjMZcNADgfQQIAAKwjSAAAgHUECQAAsI4gAQAA1hEkAADAOoIEjsdc\nNgDgfAQJAACwjiABAADWESQAAMA6ggQAAFhHkAAAAOsIEjgec9kAgPMRJAAAwDqCBAAAWEeQAAAA\n6wgSAABgHUECAACsI0jgeMxlAwDOR5AAAADrCBIAAGAdQQIAAKwjSAAAgHUECQAAsI4ggeMxlw0A\nOB9BAgAArCNIAACAdQQJAACwjiABAADWESQAAMA6gsRh0tPdmjUrWOnpvHQXMJcNADifx/YAqpqB\nA0O0Zs2VrvZaRdxWvVTL79y5QG+/fapUjwUAoLxU6SDp2DFUu3cH2R5GhVqzxqMGDYoKnPIREXFO\nn3+eV2HLAwA4U5UOEie8UYaH11JW1klJ5z+u6dEjVAUFLnk8RqmpeYqJ8VoeIQAAZVduQbJy5Urt\n27dPTzzxRHktosqJifEqNTVPGzd61K5dATECAKg0yvUIicvlKs+nr3LS093ECAAUIz3drYwMKSrK\nze9IB/IrSFauXKn169fr9OnTOnr0qBISErR27Vrt2bNH48aN0/fff6/Vq1fr9OnTCgsL0yuvvFLo\n8W+++abS0tLkcrnUtWtXDRo0qFx+GCe5spNbf3nOR+lOaP2lynKCa3R0pNxul7Zu3WF7KAAs+fkj\nbcnjCeUjbQfy+whJbm6uXnvtNX300UdatGiRli1bps2bN2vhwoWKjIzUokWLJEnDhg3Tjh0/vzHs\n3btXH330kZYsWSJjjIYMGaL27duXOBlaWFioPJ6KO+E0MlLatavCFhcQKvoE19Jq00baubP4+93u\n80fiwsMD/2exjXXkH9ZT4CnpPSEjQyooOP91QYFLGRk1FBdXQYO7ipy+7ZVl/H4HSevWrSVJtWrV\nUtOmTSVJtWvX1tmzZ1WtWjWNHTtWISEh+vHHH1VwYauQlJmZqcOHD2vw4MEyxujkyZM6cOBAiUGS\nnV2xJ5yuX1+hi/PbhZNaq/oJrVlZxd/n9Rq53S7fyb8o2sUnSKN4lX09OfUNr6T3hKgotzyen39H\nRkXlKSvLWb8jnb7t+Tv+4rZBv4OkuPNBzp49q7Vr12rZsmU6ffq0+vTpI2OM7/4bbrhBzZs316uv\nvipJWrhwoVq2bOnvYvE/nNAKAMW78DsyI6OGoqKq1h9slUWZT2r1eDwKCQnRgAEDJEkNGjTQjz/+\n6Ls/IiJCt912mwYMGKD8/HzddNNNatiwYVkXWyXFxHgVE5NvexgAEJBiYryKi5PjjozgPJe5+HBG\nAHHyYauryemH8CoCJ7X6h23JP5V9PTn1Ixt/XxMnv35OHrtU9o9smBAFjsdcNgDgfAQJAACwjiAB\nAADWESQAAMA6ggQAAFhHkAAAAOsIEjhedHRkif/nXwBAYCNIAACAdQQJAACwjiABAADWESQAAMA6\nggQAAFhHkMDxmMsGAJyPIAEAANYRJAAAwDqCBAAAWEeQAAAA6wgSAABgHUECx2MuGwBwPoIEAABY\nR5AAAADrCBIAAGAdQQIAAKwjSAAAgHUECRyPuWwAwPkIEgAAYB1BAgAArCNIAACAdQQJAACwjiAB\nAADWESRwPOayAQDnI0gAAIB1BAkAALCOIAEAANYRJAAAwDqCBAAAWEeQwPGYywYAnI8gAQAA1hEk\nAADAOoIEAABYR5AAAADrCBIAAGAdQQLHYy4bAHA+ggQAAFhHkAAAAOsIEgAAYB1BAgAArCNIAACA\ndQQJHI+5bADA+QgSAABgHUECAACsI0gAAIB1BAkAALCOIAEAANYRJHA85rIBAOcjSAAAgHUECQAA\nsI4gAQAA1hEkAADAOoIEAABYR5DA8ZjLBgCcjyABAADWESQAAMA6ggQAAFhHkAAAAOsIEgAAYB1B\nAsdjLhsAcD6CBAAAWEeQAAAA6wgSAABgHUECAACsI0gAAIB1BAkcj7lsAMD5CBIAAGAdQQIAAKwj\nSAAAgHUECQAAsI4gAQAA1hEkcDzmsgEA5yNIAACAdQQJAACwjiABAADWESQAAMA6ggQAAFhHkMDx\nmMsGAJyPIAEAANYRJAAAwDqCBAAAWEeQAAAA6wgSAABgHUECx2MuGwBwPoIEAABYR5AAAADrCBIA\nAGAdQQIAAKwjSAAAgHUECRyPuWwAwPkIEgAAYB1BAgAArCNIAACAdQQJcIXS092aNStY6ensPgBw\ntXhsDwDwx8CBIVqzpqTNtVaFjOVn1St4eVfDpeuoc+cCvf32KQtjAYCfESRVQMeOodq9O8j2MMpR\nk//9d7/FMTjXmjUeNWhQ0TEX6Cp+fUREnNPnn+dV+HKBQEGQVAGV/ZdcdLSR2+3S1q0ny31Z6elu\n9egRqoIClzweo9TUPMXEeMt9uVdDeHgtZWWV/zpyOtYTYEeJH4KvXLlSL730UqHbYmNjlZ+ff0UL\nat++/ZWNDAhAMTFepabmKSnpjKNiBAACXamOkLhcrqs9DsAxYmK8iom5siAHAFyeX0Hy1Vdf6YEH\nHlBubq5Gjhzpu33Pnj2aOnWqvF6vsrOz9eyzz+rmm2/W8uXLtXTpUhljFBsbW+gxM2fOVE5Ojp55\n5pmr/9MAV1l6ulsbN3rUrl0BR0OAAJee7lZGhhQV5WZ/dSC/giQ0NFTz5s3TsWPH1K9fPxljJJ0P\nkqefflrNmzdXWlqaVqxYoeuvv14LFizQBx98oODgYL300kvKyzt/DsO0adMUFBREjMAv/l1ZI0nn\nj9iV74mZzruihqtnUJX8fH6X5PGE8pGqA/kVJNHR0ZKkunXrqlatWjpw4IAkqWHDhpozZ45CQkKU\nk5OjmjVr6tChQ2rRooWCg4MlSWPHjpUkHT16VJmZmWrcuLFfAwsLC5XHU5mvDPFfePiVvdFGRkq7\ndpXTYALSftsDCEhFXz3D1TT+Cez11KaNtHOn7VFUrJLeEzIypIKC818XFLiUkVFDcXEVNLir6Ep/\n3weasozfryDJyMiQJGVlZSkvL09hYWGSpJSUFL3wwgtq2rSpZs+ercOHD6tRo0bat2+fzp49q2rV\nqunxxx/XhAkTVL9+fb322mtKSEjQF198oQ4dOlx2mdnZlfvKEH+V5oz/9evLaTABrDyujHDyFTVF\n4eoR/zhlPWVlle5xTn3DK+k9ISrKLY/n5/01KipPWVnO2l+dsu0Vx9/xF7cN+hUkZ86c0eDBg3Xq\n1ClNnjxZEyZMkCT16NFDo0aNUu3atdWwYUMdP35cdevW1YMPPqhBgwbJ5XIpNjZWDRs29D1XSkqK\nHnroIb3zzjuqXbu2P4sHrLhwRQ3nkACB78L+mpFRQ1FRzv7joapymQsnhAQYJ1fi1eT0Yq4orKeS\nsY78U9nXk1OPkPj7mjj59XPy2KWyHyFhMg4AAGAdQQIAAKwjSOB40dGRatKkie1hAADKgCABAADW\nESQAAMA6ggQAAFhHkAAAAOsIEgAAYB1BAsfbtm2n9u/fb3sYAIAyIEgAAIB1BAkAALCOIAEAANYR\nJAAAwDqCBAAAWEeQwPGYywYAnI8gAQAA1hEkAADAOoIEAABYR5AAAADrCBIAAGAdQQLHYy4bAHA+\nggQAAFhHkAAAAOsIEgAAYB1BAgAArCNIAACAdQQJHI+5bADA+QgSAABgHUECAACsI0gAAIB1BAkA\nALCOIAEAANYRJHA85rIBAOcjSAAAgHUECQAAsI4gAQAA1hEkAADAOoIEAABYR5DA8ZjLBgCcjyAB\nAADWESQAAMA6ggQAAFhHkAAAAOsIEgAAYB1BAsdjLhsAcD6CBAAAWEeQAAAA6wgSAABgHUECAACs\nI0gAAIB1BAkcj7lsAMD5CBIAAGAdQQIAAKwjSAAAgHUECQAAsI4gAQAA1hEkcDzmsgEA5yNIAACA\ndQQJAACwjiABAADWESQAAMA6ggQAAFhHkMDxmMsGAJyPIAEAANYRJAAAwDqCBAAAWEeQAAAA6wgS\nAABgHUECx2MuGwBwPoIEAABYR5AAAADrCBIAAGAdQQIAAKwjSAAAgHUECRyPuWwAwPkIEgAAYB1B\nAgAArCNIAACAdQQJAACwjiABAADWESRwPOayAQDnI0gAAIB1BAkAALCOIAEAANYRJAAAwDqCBAAA\nWEeQwPGYywYAnI8gAQAA1hEkAADAOoIEAABYR5AAAADrCBIAAGAdQQLHYy4bAHA+ggQAAFhHkAAA\nAOsIEgAAYB1BAgAArCNIAACAdQQJHI+5bADA+QgSAABgHUFSjtLT3Zo1K1jp6axmAAAux2N7AE4x\ncGCI1qwp7eqqXurl3nOPtHBhqR8OAIAjVMog6dgxVLt3B9kexlXx0UdSgwa1bA+jkIiIc/r88zzb\nwwAAVCKVMkgC4c0yPd2tHj1CVVDgksdjlJqap5gY7xU/T3h4LWVlnSyHEQIAEDhKdXJDZmam0tPT\nJUmxsbHKz8+/qoOqDGJivEpNzVNS0plSxwj8w1w2AOB8pTpCsnr1aoWHhysmJkYul+tqj6nSiInx\nKiaGWAMAoCQlBklBQYHGjx+vQ4cOyRijAQMGaMWKFQoODlarVq1kjNGzzz6rQ4cOyeVyac6cOQoJ\nCdGkSZN08OBBeb1ejR49Wr/97W/VvXt3NWnSRMHBwXrxxRcr4uezKj3drY0bPWrXroAjJABQztLT\n3crIkKKi3PzOdaASg2TZsmWqV6+eZsyYodzcXPXp00exsbFq3ry5oqKiJEn9+vVT27ZtNX78eG3Y\nsEHZ2dmqW7euUlJSdPz4cQ0aNEhpaWnKzc3Vo48+qoiIiHL/wcpD6a+04SobAChPP5+3J3k8oXxU\n7kAlvrvu3btX7dq1kyTVqFFDzZo108GDB9W8eXPf97Rp00aSVL9+fZ0+fVqZmZnatm2btm/fLmOM\nzp07p+zsbEnSDTfc4NfAwsJC5fGU/kqZyEhp165SPzxgBOJVNhdr00baudP2KM4LDw/c9RQoWEf+\nYT0FnpLeEzIypIKC818XFLiUkVFDcXEVNLiryOnbXlnGX2KQNGvWTOnp6ercubNycnKUmZmpPn36\nyOstvjybNWum6667TsOHD9eZM2f0t7/9TXXq1JEkv885yc4u25Uy69eX6eFlVpWussnKsj0CZ6wn\n21hH/qns68mpb3glvSdERbnl8fz8OzcqKk9ZWc46QuL0bc/f8Re3DZYYJPfdd5+eeeYZDRw4UGfO\nnNHIkSMVFhamGTNmqGnTpoUC48LX/fv3V1JSkhISEpSbm6sBAwbI5XJVqRNgL1xlwzkk5S86OlJu\nt0tbt+6wPRQAllz4nZuRUUNRUXxc40QuY4yxPYiiOLkSryanF3NFIEj8w7bkn8q+npx6hMTf18TJ\nr5+Txy6V/QgJk6wAAADrCBIAAGAdQQIAAKwjSAAAgHUECRyPuWwAwPkIEgAAYB1BAgAArCNIAACA\ndQQJAACwjiABAADWESRwvOjoSDVp0sT2MAAAZUCQAAAA6wgSAABgHUECAACsI0gAAIB1BAkAALCO\nIIHjMZcNADgfQQIAAKwjSAAAgHUECQAAsI4gAQAA1hEkAADAOoIEjsdcNgDgfAQJAACwjiABAADW\nESQAAMA6ggQAAFhHkAAAAOsIEjgec9kAgPMRJAAAwDqCBAAAWEeQAAAA6wgSAABgHUECAACsI0jg\neMxlAwDOR5AAAADrCBIAAGAdQQIAAKwjSAAAgHUECQAAsI4ggeMxlw0AOB9BAgAArCNIAACAdQQJ\nAACwjiABAADWESQAAMA6ggSOx1w2AOB8BAkAALCOIAEAANYRJAAAwDqCBAAAWEeQAAAA61zGGGN7\nEAAAoGrjCAkAALCOIAEAANYRJAAAwDqCBAAAWEeQAAAA6wgSAABgHUECAACsI0gCgDFGkyZNUnx8\nvBITE3Xo0KFC9y9cuFDdunVTYmKiEhMTtX//fjsDDRDbt29XQkLCJbevW7dOffv2VXx8vJYvX25h\nZIGluPXE9iQVFBRo3Lhxuv/++3Xfffdp3bp1he5nW3Km4rb5QOT0/bOkfag0PFdhXCijNWvWKD8/\nX0uXLtX27ds1ZcoUzZ0713f/rl27NH36dLVu3driKAPDggUL9P7776tGjRqFbi8oKNDUqVO1YsUK\nVa9eXQMGDNBdd92lunXrWhqpXcWtJ4ntSZJSU1MVFham6dOn67///a969eql2NhYSWxLTnW5bT7Q\nVIb983L7UGlxhCQAbNu2TR06dJAk3XTTTdq5c2eh+3ft2qV58+Zp4MCBmj9/vo0hBozGjRtrzpw5\nl9y+d+9eNW7cWDVr1lS1atUUHR2trVu3WhhhYChuPUlsT5IUFxenUaNGSZK8Xq88np//NmNbcqbL\nbfOBpjLsn5fbh2JjY5Wfn3/Fz0mQBICcnBzVqlXL92+PxyOv1+v7d9euXfXcc8/p73//u7Zt26bP\nPvvMxjADQpcuXRQUFHTJ7b9chzVq1NDJkycrcmgBpbj1JLE9SVJISIhCQ0OVk5OjUaNGacyYMb77\n2Jac6XLbfKCpDPvnL/eh0aNHa/z48UpISNBPP/2kYcOGaciQIVf0nARJAKhZs6Zyc3N9//Z6vXK7\nf35pBg8erDp16sjj8ahTp076+uuvbQwzoNWsWVM5OTm+f+fm5uqaa66xOKLAxfZ03pEjRzR48GD1\n7t1b99xzj+92tiXY5KT98+J9qGvXrpoyZYoWL16sevXq6fXXX9cbb7xxRc9HkASAW265xVfB//73\nv9WiRQvffTk5OerWrZtOnTolY4w2bdqkNm3a2BpqwPjlnJDNmjXTgQMHdOLECeXn52vr1q26+eab\nLY0ucPxyPbE9nXf06FENGzZMTz75pHr37l3oPrYlZ3PSfLFO3j8vtw+5XK5SvQ6c1BoAunTpog0b\nNig+Pl4WUkF0AAAAsElEQVSSNGXKFKWlpenUqVPq16+fxo4dq4SEBFWvXl233367OnbsaHnE9rlc\nLkkqtJ7Gjx+voUOHyhijfv36qUGDBpZHaV9R64ntSZo3b55OnDihuXPnas6cOXK5XLrvvvvYliqB\nC9u8Ezh5/yxqH1qwYIGCg4O1du3aUj2nyzgpJwEAQKXERzYAAMA6ggQAAFhHkAAAAOsIEgAAYB1B\nAgAArCNIAACAdQQJAACw7v8BNEEDn+Xm/BgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b8de4e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(trace[burn:], varnames=['odds'], ylabels=covs, vline=1, \n",
    "              main='Odds ratios of covariates')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Age effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.gridspec.GridSpec at 0x11c8fa320>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiAAAAF7CAYAAADrKDC1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VOW9xvFnhgEkXIRKYk+7XFSlEoWOCuNBEUEuIUYS\nYOQioICW4+UUCxoUpHCwhqCirdIEEBqraWlpvaJRlgcbBVGQy+RU08BBPUigXhtp0MCgyTDv+SMS\nCIRkIDvvzCTfz1quFWbv2e8vb/a797PfPc52GWOMAAAALHJHuwAAANDyEEAAAIB1BBAAAGAdAQQA\nAFhHAAEAANYRQAAAgHUEEDhi5cqVuuaaa+T3+zVz5kx99dVXNcsuv/xy+f3+mv9eeeUVSdJf/vIX\nDR8+XOPHj9cnn3xSs/6tt96qjz76qN72/ud//kf/8R//Ib/fr4yMDN1+++368MMPG/17LFiwQEuW\nLJEk3Xbbbdq1a5e2bt2qjIyMOtefM2eOnnrqqVNqY8mSJcrOzm5wvaVLl+qNN944pW2fjvp+PwAN\nS05O1ogRIzRq1Cj5/X5dc801Gjt2rEpKSupcP9JjwPGmTp2q/fv3N7bcmOGJdgGIf5s3b9bvfvc7\nPfPMM0pKStJLL72k//qv/1JOTo52796tzp07a/Xq1Se8Ly8vT2vXrtVrr72mP/3pT5o1a5b++7//\nW927d9d555130va2bdumWbNmadmyZbrwwgslSS+//LImTZqkV199VV26dHHk91qxYoUkad++fY5s\n71Rt3rxZP/7xj6PSNoDIuVwurVy5UmeeeWbNa08++aSys7P1l7/8xbF2Nm7c6Ni2YgEBBI22Y8cO\nXXHFFUpKSpIkDRs2TPPmzVMoFNLf/vY3ud1uTZ48Wfv371dqaqp+9rOfyeVyqXXr1vrmm28UDAZr\nfn7qqacanFHIzc3VtGnTasKHJGVkZOiMM85QOBzW1q1btXDhQrVr107ffPONnn32Wb311ltavny5\nQqGQzjjjDM2aNUuXXHKJDhw4oHnz5un9999XYmKiWrVqpT59+kiSBg8erNzcXEnSwYMHNX36dO3d\nu1edOnXSggUL1K1bt1p17dq1Sw888ID279+vcDisSZMm6brrrqv3d5kzZ47at2+vDz74QJ9//rnO\nO+88PfbYY3rhhRdUUlKihx9+WG63WwMHDtSvfvUrbdu2TeFwWBdeeKHmzZun9u3ba/Dgwbr44ov1\nwQcfaNq0aXr88cf18ssvS5IqKio0ZMgQvf766woEAlqxYoVCoZD+9a9/aeTIkZoxY0ategKBgBYt\nWqRwOCyXy6XbbrtNKSkpEewFQMtljNGx3+l5+PBhffrpp+rcufNJ37Nr1y5NnjxZZWVl6tq1qx57\n7DF17dpV69atqzVOR40apenTp2vOnDmSpMmTJysvL09nn312k/9eTc4AjbRt2zYzaNAg8+mnnxpj\njFm5cqVJTk42ZWVl5plnnjHZ2dmmqqrKVFRUmPHjx5vf//73xhhj1q5da/x+v7nlllvMvn37zGOP\nPWZeeumlBtu79NJLzf/93/+ddPmWLVvMRRddZD777DNjjDGlpaUmPT3d7N+/3xhjzIcffmiuvPJK\nc+jQIbNw4UJz7733GmOM2bdvnxk4cKDJzc01xhgzaNAgU1JSUrO9d9991xhjzNNPP23Gjh1rjDHm\n3nvvNU8++aQJhUJm+PDhZseOHcYYYyoqKsy1115r3nvvvRPqy83NNQsWLKh5/4QJE0xVVZWpqqoy\nfr/fvPDCC8YYY2688Ubz2muvGWOMWbJkiXn44YdrtvHoo4+a+++/v6bOZcuW1SwbMmSIKSkpMcYY\ns2rVKnPPPfcYY4yZPHmy2bNnjzHGmC+++MJcdNFFpry83GzZssWkp6cbY4yZMmWKWbNmjTHGmJ07\nd5qsrKz6/xgATI8ePUxGRoYZMWKE6d+/vxkyZIjJzs42+/btq3P93NxcM3ToUFNeXm6MMeZnP/tZ\nzRg+2Tg90s6R41hzwAwIGs3n82natGmaNm2a3G63Ro8erTPPPFOtW7fW2LFja9br0KGDbr75Zq1c\nuVKTJ0/WsGHDNGzYMEnS3r179d5772nGjBl64IEH9Mknn+iyyy7TTTfddEJ7bre71tVGXb7//e/r\n+9//vqTqacsvv/xSN910U837PB6PSktL9c4772ju3LmSpO9973saOnRondvr0aOHLr74YkmS3+/X\n/fffrwMHDtQsLy0t1d69e/WLX/yipo1vv/1WO3bskNfrrbfWq666Sh5P9VC84IILan1+5si21q9f\nr4qKipop2FAopLPOOqtmPZ/PV/Pz6NGjtXr1avXs2VMvvPCCZs2aJUl6/PHHtX79ehUUFNR8xubQ\noUO1aklLS1NWVpbeeOMN9evXT3fddVe9tQOoduQWzP/+7//qlltu0aWXXqrvfe97J12/X79+NTMk\nycnJNbd6TzZOj6zb0LEvnhBA0GgHDx7UZZddptGjR0uq/szEb37zG5155pl66aWXlJycrB49ekiq\nHjytW7c+YRsPPfSQZs+erU2bNikYDGrp0qWaOnWqhgwZonPOOafWupdccon+9re/qXv37rVez8rK\nUkpKilq1aqWEhISa18PhsK644go9+uijNa99/vnnSkpKksvlqjWgjwSB47ndRz+vbYypuYV0xOHD\nh9WpU6dan3XZt2+fOnbsePKO+84ZZ5xR8/Px9Ry7/blz5+qqq66SVH1A+vbbb2uWH/v7Xnfddbru\nuus0ZswYVVRU6LLLLtOhQ4fk9/uVkpIin8+nMWPGqLCw8IS2rr/+eg0ePFgbN27Uhg0btGTJEhUU\nFKhDhw4N/h5AS3ZkLF144YWaM2eO5s6dq0suuUQ/+MEPdOutt+qf//ynXC6Xpk+fLkm1jh8ul0tS\n9bgeNWqUhg0bVu84bS74v2DQaP/85z81adKkmhmBZcuWKT09XZL04YcfKjc3V+FwWN98843++Mc/\n6tprr631/nXr1unss89WcnKyKisra4WAb7755oT2br/9di1btkw7duyoee2FF17Qa6+9VhN0jnX5\n5Zdr48aNNVcTb775pkaOHKnKykpdddVVeu6552SM0VdffaXXX3+9zt9x586d2rlzpyTp6aefVu/e\nvdW2bdua5eeee67atm2rgoICSdJnn32m9PR0bd++veEOPAmPx6NQKCSpepbkT3/6k6qqqhQOhzV3\n7txagepYZ599tn7yk59o/vz5NTNQe/bs0cGDB3XnnXfq6quv1pYtW1RVVaXDhw/Xeu/48eO1Y8cO\njRo1SllZWaqoqNDXX3992r8D0BINHz5cvXv31sKFCyVJv/3tb/Xiiy9q9erVGjRo0Enft2fPHgWD\nwZOO02OPCc0BMyBotHPPPVe33nqrxo0bJ2OM+vTpo/nz50uS7rjjDi1YsEAZGRkKhUJKS0vTmDFj\nat5bWVmp5cuXKy8vT5J05ZVXatWqVRo1apQuvvjiOv8vEJ/Pp+zsbGVnZ+vQoUOqqqrSOeecoz/8\n4Q91Tnl2795dWVlZyszMlCS1atVKjz/+uM444wz9/Oc/13333ae0tDSdddZZtQLMkasSSTr//PO1\ndOlS7d27V127dtWiRYtqtdG6dWstW7ZM2dnZeuKJJ3T48GHddddduvTSS0+7XwcNGqRFixapsrJS\n06ZN00MPPSS/31/zIdTZs2efUOcR48aN04wZM7R8+XJJ1beQrr76al1zzTXq1KmTunXrpu7du2vv\n3r21rsTuueceLVy4UL/5zW/kcrl0xx136Ac/+MFp/w5AS1DXGJw3b55GjhypjRs36sorr4xoO8nJ\nyRo4cGCd4/Scc87R0KFDNXHiRC1btuyEGeB45DLNdW4HAADELG7BAAAA6wggAADAOgIIAACwLiY+\nhFpWVuH4Nrt0SVB5edDx7TZH9FVk+vTpJbfbpW3b/h7tUuIC+9Wpaar+Skxs+H8Fj0WRnBfifR9r\nCfXXt/812xkQj6dVtEuIG/RVZIqKSlRaWhrtMuIG+9Wpob9OXbz3WUuvv9kGEAAAELsIIAAAwDoC\nCAAAsI4AAgAArCOAAAAA6wggQIT69OmlH/3oR9EuAwCaBQIIAACwjgACAACsI4AAAADrCCAAAMA6\nAggAALCOAAJEiGfBAIBzCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCBAhngUDAM4hgADN\nTCDgVk5OGwUCDG8AscvTmDfPmTNHO3bs0C9+8Qv17dtXhw8f1l133aVx48apf//+Wrx4sZ577jk9\n9NBD6t+/v1M1AzFr4sR2Kixs1LByUNtoF6ChQ9tp1apD0S4DQAxq9JHynnvuUd++ffWPf/xDs2bN\n0hdffKFx48ZJku6880598cUXjS4SzcuAAQnaubNVtMs4DS5JUlJSxyjXET8KCz30l6Tk5MPasCEY\n7TKalUDAreJiyet1y+cLR7scnIZ6A8jMmTM1YsQIDRw4ULt27dLDDz+sFStW1LluMBjUwoULlZeX\n1ySFovmI1wNxnz5GbrdL27ZVnNL7YmtWJPqGDg0xK4JGCQTcGjEiQaGQ5PEkqKAgSAiJQ/UeFceN\nG6c///nPGjhwoJ5//nmNHTv2pOv26NHjtIvo0iVBHo/zV8SJiVx5RSoafdWrl7R9u/VmG4EZECcw\nK1KfE/ulZ0+ppCQKpURZfeeF4mIpFKr+ORRyqbi4vdLSLBbnoHg/TzWm/noDSN++fZWdna1//etf\n2rhxo2bOnHnaDdWnvNz5K+LExI4qKzu1K9WWKlp9tW6d9SYb6e9xuV8dvVp0yeMx1q4W47Gvoqm+\n/iora9x241F95wWv1y2P5+g+7fUGVVYWfzMg8T5GIqm/vv2vwXnhkSNHauHCherfv79atYrH+/ZA\ny+bzhVVQENSmTR716xdiqhpx78g+XVzcXl4vt1/iVYMBxO/3a/HixXrllVds1AOgCfh8Yfl8ldEu\nA3CMzxdWWpricuYD1RoMIKFQSJdddtlJv4DJGFPr3w8++KAjhQEAgOar3m8q+utf/6pbbrlF06dP\nP+k6v/rVr7Rly5Y6ly1evFhvvfVW4yoEAADNjsscP4URBU3xIZx4/3CPTfRV5OiryNFXp6ap+ite\nP4QaSV/E+z7WEuqvb//ju5qBCPEsGABwDgEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEi\nVFRUotLS0miXAQDNAgEEAABYRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEixLNgAMA5BBAAAGAd\nAQQAAFhHAAEAANYRQAAAgHUEEAAAYB0BBIgQz4IBAOcQQAAAgHUEEAAAYB0BBAAAWEcAAQAA1hFA\nAACAdQQQIEI8CwYAnEMAAQAA1hFAAACAdQQQAABgHQEEAABYRwABAADWEUCACPEsGABwDgEEAABY\nRwABAADWEUAAAIB1BBAAAGAdAQQAAFhHAAEixLNgAMA5BBDgNAQCbuXktFEgwBACgNPhaewGVq9e\nrZycHI0ZM0aBQEChUEiSlJWVpQ8++ECLFy9WSkqKMjMzG10sIEkTJ7ZTYWGjd93T4JIkJSV1POa1\nttZaHzo0pFWrDllrDwCakiNH8YyMDH388ce68cYbNWTIEL399tt69NFHlZubq2AwqN27dzvRTJMZ\nMCBBO3e2inYZUdax4VUQVYWFnuPCTzxofL3JyYe1YUPQgVrQnAQCbhUXS16vWz5fONrl4DQ4dhk5\ne/ZsdexYfbAJhUJq29belWFjtfSDW2JiR5WVVUS7jBNEb6YjfsTyrEis7leIf4GAWyNGJCgUkjye\nBBUUBAkhccixo3vnzp0lSR999JEeeeQRLVu2LOL3dumSII/H+RmIxMSTX3316iVt3+54k3Es3q6s\nIcXDrIid2nr2lEpKrDTVpOo7ZrU09Z0Xioul7+72KxRyqbi4vdLSLBbnoHj/mzemfkcvLzdv3qwF\nCxbokUceUbdu3SJ+X3m58zMQDV19rVvneJNxiyvVSP1diYkd9eqrB7+7+nLJ4zFcfZ2E7f2qrMxa\nU02iqforXk9w9Z0XvF63PJ6jY9DrDaqsLP7GYLwfeyOpv779z7EAsmXLFj3wwAN64okn9G//9m9O\nbRaIOT5fWAUFQW3a5FG/fiHCB2DZkTFYXNxeXi8XAPHKkQBijNGDDz6oUCik2bNnyxij8847T/ff\nf78Tmwdijs8Xls9XGe0ygBbL5wsrLU1xOfOBao7NgLz44otObQoAADRzjnyL0po1a5Sfn3/C62vX\nrlVeXp4TTQAAgGak0TMgfr9ffr+/zmWpqalKTU1tbBMAAKCZ4XukgQjxLBgAcA4BBAAAWEcAAQAA\n1hFAAACAdQQQAABgHQEEAABYRwABIlRUVKLS0tJolwEAzQIBBAAAWEcAAQAA1hFAAACAdQQQAABg\nHQEEAABYRwABIsSzYADAOQQQAABgHQEEAABYRwABAADWEUAAAIB1BBAAAGAdAQSIEM+CAQDnEEAA\nAIB1BBAAAGAdAQQAAFhHAAEAANYRQAAAgHUEECBCPAsGAJxDAAEAANYRQAAAgHUEEAAAYB0BBAAA\nWEcAAQAA1hFAgAjxLBgAcA4BBAAAWEcAAQAA1hFAAACAdQQQAABgHQEEAABYRwABIsSzYADAOQQQ\nwGGBgFs5OW0UCDC8AOBkPI3dwOrVq5WTk6Prr79emzdvVlVVlTp37qxHHnlEb731lhYvXqyUlBRl\nZmY6US9wSiZObKfCwkbv5t9xSZKSkjpGuH5bh9o9uaFDQ1q16lCTtwMATnPkyJyRkaHy8nL5/X6N\nHDlSS5Ys0bPPPqspU6YoGAxq9+7dTjSDOgwYkKCdO1s5sKVIT6qIJYWFnlMIRNEQy7U1neTkw9qw\nIRjtMpq1QMCt4mLJ63XL5wtHuxycBqcuDTVnzhxJUjgc1meffaYf/vCHTm0a9XDiIJeY2FFlZRUO\nVBMbnJ31iD+xMivS3PYrxI5AwK0RIxIUCkkeT4IKCoKEkDjk6FE6FApp5MiRqqys1B133BHx+7p0\nSZDH48RVfG2JifFx9dWrl7R9e7SriI++QsNia1YkVuqwq2dPqaTk1N8XL8csG+o7LxQXS6FQ9c+h\nkEvFxe2VlmaxOAfF+9+8MfU7GkA8Ho/WrFmjd955R7NmzdLKlSsjel95ufNTlfF09bVuXXTbj6e+\niq6/19tXR6/KXPJ4TIu/Kmvp+1VZ2amt31T9Fa8nuPrOC16vWx7P0bHm9QZVVhZ/Yy3ex0gk9de3\n/zkSQIwxysrKUmpqqvr27auEhAS53fwfAGhZfL6wCgqC2rTJo379Qi06fABN6chYKy5uL6+3ZQf9\neOZIAHG5XJo0aZLmz5+vZcuWye1267777nNi00Bc8fnC8vkqo10G0Oz5fGGlpSkuZz5QzbEZkHPP\nPTfiWy4AAKBlc+Q+yZo1a5Sfn3/C62vXrlVeXp4TTQAAgGak0TMgfr9ffr+/zmWpqalKTU1tbBMA\nAKCZ4ZOiQIR4FgwAOIcAAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAESoqKlFpaWm0ywCA\nZoEAAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAEeJZMADgHAIIAACwjgACAACsI4AAAADr\nCCAAAMA6AggAALCOAAJEiGfBAIBzCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCBAhngUD\nAM4hgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggQIR4FgwAOIcAAgAArCOAAAAA6wggAADA\nOgIIAACwjgACAACsI4AAEeJZMADgHAIIcJoCAbdyctooEGAYAcCp8jR2A6tXr1ZOTo6mTJmim266\nSVu3btWsWbO0fv16rV27VosXL1ZKSooyMzOdqBc4JRMntlNhYaN38++4JElJSR2Pe72tQ9uv39Ch\nIa1adchKWwDQ1Bw5MmdkZOimm27S559/rvz8fIVCIUlSamqqgsGgdu/e7UQzcNCAAQnaubPVMa8c\nf1JFrCks9NQRfmJdw/UmJx/Whg1BC7UAiCVOXRqqsrJSv/zlL5WVlaXrrrvOqc2iiRx7wE9M7Kiy\nsoooVhMf+vQxcrtd2ratQoGAWyNGJCgUcsnjMXr44W/0r3+51a9fSD5fONqlxgT2KzSlQMCt4mLJ\n63Uz5uKUIwHEGKOsrCz99Kc/VVJSkhObBBzj3G2Yum/BhEIuZWa2++5fTXs7htswgI65AJA8ngQV\nFAQJIXHIkQDy1VdfqaioSHv37pUxRvv379fMmTP161//OqL3d+mSII+nVcMrnqLExHibrm46vXpJ\n27fXtwZ91bDSaBcQh7dhmrbWnj2lkpImbcIqjllH1XdeKC6WvrvTr1DIpeLi9kpLs1icg+L9b96Y\n+h0JIJ07d9arr75a8+/+/ftHHD4kqbzc+fu/TP/Wtm7dyZfRV5Grq6+Ovx3D1Vg1W/tVWVmTN2FF\nU/VXvJ7g6jsveL1ueTxHx5zXG1RZWfyNuXg/9kZSf337n2OfAQFaKp8vrIKCoDZt8vAZEMCCI2Ou\nuLi9vF4Cf7xy7DMgx3r77bed2CwQN3y+sHy+ymiXAbQYPl9YaWmKy5kPVHPkG5TWrFmj/Pz8E15f\nu3at8vLynGgCAAA0I42eAfH7/fL7/XUuS01NVWpqamObAAAAzQzfIQ1EiGfBAIBzCCAAAMA6AggA\nALCOAAIAAKwjgAAAAOsIIAAAwDoCCBChoqISlZaWRrsMAGgWCCAAAMA6AggAALCOAAIAAKwjgAAA\nAOsIIAAAwDoCCBAhngUDAM4hgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggQIR4FgwAOIcA\nAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAEeJZMADgHAIIAACwjgACAACsI4AAAADrCCAA\nAMA6AggAALCOAAJEiGfBAIBzCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCBAhngUDAM4h\ngACnIRBwKyenjQIBhhAAnA5PYzewevVq5eTkaMqUKVq+fLkuuOACSVJKSoqSkpK0ePFipaSkKDMz\ns9HFApGaOLGdCgsbvXsfxyVJSkrqeMxrbR1u46ihQ0NatepQk20fAKLJkSN0RkaGevToofT0dM2b\nN6/WsmAwqN27dzvRDBwyYECCdu5sddyrHetcF9FTWOg5LuzEruTkw9qwIRjtMgDEEccuEUtKSlRS\nUqJJkybprLPO0rx589S1a1enNg8HHX+iSEzsqLKyiihVEz/69DFyu116/PGDGjEiQaGQS61aGU2a\nVKVx46rk84WjXSLQYgQCbhUXS16vm7EXpxwLIOeff7569eqlK664Qi+//LKysrKUk5Pj1OaB0+bc\n7ZjqWzDXXtu+5pXDh13Kz2+j/Pw2Dmz/5LgdAxwVCLi/uwiQPJ4EFRQECSFxyLEA0rdvX7Vr105S\n9ec/cnNzI35vly4J8niOvyXQeImJ8TF9HQvq66tevaTt2y0WE7NKo9ZyPN2OqS3ymnv2lEpKmrCU\nOMAx66j6zgvFxVIoVP1zKORScXF7paVZLM5B8f43b0z9jgQQY4zmzZunYcOGKS0tTZs2bVLPnj0j\nfn95ufP3jrmtELmG+mrdOovFxLhj++roVZhLHo/hKuw4pzMGy8qaqJg40FTHrHg9wdV3XvB63fJ4\njo49rzeosrL4G3vxfp6KpP769j9HAojL5dLdd9+tOXPm6M9//rMSEhKUnZ3txKaBmOXzhVVQENSm\nTR716xcifACWHBl7xcXt5fUS/OOVYzMgP/zhD/WHP/zBic0BccPnC8vnq4x2GUCL4/OFlZamuJz5\nQDVHvkVpzZo1ys/PP+H1tWvXKi8vz4kmAABAM9LoGRC/3y+/31/nstTUVKWmpja2CQAA0MzwPdJA\nhHgWDAA4hwACAACsI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAARKioqUWlpabTLAIBmgQACAACs\nI4AAAADrCCAAAMA6AggAALCOAAIAAKwjgAAR4lkwAOAcAggAALCOAAIAAKwjgAAAAOsIIAAAwDoC\nCAAAsI4AAkSIZ8EAgHMIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIECGeBQMAziGAAAAA\n6wggAADAOgIIAACwjgACAACsI4AAAADrCCBAhHgWDAA4hwACAACsI4AAAADrCCAAAMA6AggAALCO\nAAIAAKwjgAAR4lkwAOAcAgjggEDArZycNgoEGFIAEAlPYzewevVq5eTk6Prrr9fu3bv18ccfKxQK\nad68efr000+1ePFipaSkKDMz04l6gVM2cWI7FRY2eleX5JIkJSV1rGedtg60E7mhQ0NateqQ1TYB\nwAlOHJWVkZGhUCikCy64QIsWLdL777+v999/XyNGjFAwGNTu3budaAZRMGBAgnbubBXtMnAShYWe\nBgJRtNmvLTn5sDZsCFpvF8CpcSSAGGP09ttvKy0tTVOnTlXHjh01f/58JzaNKONAflSfPkZut0vb\ntlXUej0QcGvEiASFQi55PEYFBUH5fGEFAm5t2uRRv34h+XzhKFUdPYmJHVVWVtHwisBpCATcKi6W\nvF53ixxfzYEjAUSSysvLVVFRod/97nd68cUXtWjRIi1atMipzQOOOf1bMg3fggmFXLr22vbHvWr3\ntozErRk0b0dDv+TxJNSEfsQXxwJIly5dNHjwYEnS4MGD9cQTT5zCexPk8Tg/zZ+YGMtT07Glob7q\n1Uvavt1SMTGrNNoFRCx2bs3Yq6FnT6mkxFpzTYJj1lH1nReKi6VQqPrnUMil4uL2SkuzWJyD4v1v\n3pj6HQsgvXv31vr163XRRRdp69at6t69e8TvLS93fpqf6d/IRdJX69ZZKibGRbpfney2TEsSjTFY\nVma1OUc1VX/F6wmuvvOC1+uWx3N0fHm9QZWVxd/4ivfzVCT117f/ORJAXC6Xbr/9ds2dO1fjx49X\n69atuf2CFs3nC6ugINiiPwMCNJUj46u4uL283pYX7psLxz6E2qlTJ+Xm5jqxOaBZ8PnC8vkqo10G\n0Cz5fGH8SLkcAAAJn0lEQVSlpSkuZz5QzZFvTVqzZo3y8/NPeH3t2rXKy8tzogkAANCMNHoGxO/3\ny+/317ksNTVVqampjW0CAAA0M3xvNBAhngUDAM4hgAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA\n6wggQISKikpUWloa7TIAoFkggAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggQIR4FgwAOIcA\nAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsI4AAEeJZMADgHAIIAACwjgACAACsI4AAAADrCCAA\nAMA6AggAALCOAAJEiGfBAIBzCCAAAMA6AggAALCOAAIAAKwjgAAAAOsIIAAAwDoCCBAhngUDAM4h\ngAAAAOsIIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggQIR4FgwAOIcAgmYjEHArJ6eNAgF2awCIdZ7G\nbmD16tXKyclRRUWFLrroIhlj9OWXX+rMM8/UzTffrMWLFyslJUWZmZlO1IsYNXFiOxUWNnp3ckjb\nJtquS5KUlNSx3rWGDg1p1apDTVQDADQPjpwxMjIyagJGKBTSDTfcoOzsbHXv3l3BYFC7d+92opnT\nNmBAgnbubBXVGmJf/SdVRK6w0NNgSIlHycmHtWFDMNplAGgmHL9kXblypa688kp1797d6U2fNg6a\n9UtM7Kiysopol9Eof/yjR3fffYbCYZc8HqOCgqB8vrCjbfTpY+R2u7RtW3z3FdAcBAJuFRdLXq/b\n8bEOOxwNIFVVVXr66af13HPPOblZxIDYusVSv1DIpWuvbd8EW274Fgy3X4CmFwi4NWJEgkIhyeNJ\naJILDjQ9R88omzZt0r//+7+rQ4cOp/S+Ll0S5PGc/i2SXr2k7dvrWtL8psGbDn3VsNIG12iut19O\nX+2+6NlTKimJUilxIDGRfeeI+s4LxcVSKFT9cyjkUnFxe6WlWSzOQfH+N29M/Y4GkHfeeUcDBgw4\n5feVlzfuFsm6dSe+1hxuK9gS73119Gqo6W6/HBHvfWXTyfqqrCwKxcSBptq34vUEV995wet1y+M5\nOua93qDKyuJvBiTejyeR1F/f/udoACktLdWoUaOc3CTQIJ8vrIKCoDZt8qhfvxBTsUAzd2TMFxe3\nl9fL7Zd45UgAMcZIkpYvX+7E5oBT5vOF5fNVRrsMAJb4fGGlpSkuZz5QzZFvbFqzZo3y8/NPeH3t\n2rXKy8tzogkAANCMNHoGxO/3y+/317ksNTVVqampjW0CAAA0M3xnNRAhngUDAM4hgAAAAOsIIAAA\nwDoCCAAAsI4AAgAArCOAAAAA6wggQISKikpUWloa7TIAoFkggAAAAOsIIAAAwDoCCAAAsI4AAgAA\nrCOAAAAA6wggQIR4FgwAOIcAAgAArCOAAAAA6wggAADAOgIIAACwjgACAACscxljTLSLAAAALQsz\nIAAAwDoCCAAAsI4AAgAArCOAAAAA6wggAADAOgIIAACwjgACAACsa1YB5K9//atmzpxZ57JnnnlG\no0eP1vjx47V+/Xq7hcWYb7/9VtOnT9cNN9yg2267TeXl5Sess3DhQo0ePVqTJ0/W5MmTdeDAgShU\nGj3GGN13330aP368Jk+erH/84x+1lr/xxhsaM2aMxo8fr2effTZKVcaGhvoqPz9f6enpNftSaWlp\ndAqNIe+9954mTZp0wuvsV5E7WR/GqpPVGw/jIxQKadasWbrhhhs0btw4vfHGG85s2DQT2dnZJi0t\nzWRmZp6wrKyszKSnp5uqqipTUVFh0tPTTWVlZRSqjA1PPfWUyc3NNcYYs2bNGpOdnX3COhMmTDDl\n5eW2S4sZr732mrn33nuNMca8++675j//8z9rllVVVZmUlBRTUVFhKisrzejRo82+ffuiVWrU1ddX\nxhhz9913m+3bt0ejtJiUl5dn0tPTzfXXX1/rdfaryJ2sD2NVffXGw/h4/vnnzQMPPGCMMWb//v3m\n6quvdmS7zWYGpHfv3vrlL39Z57Li4mL16dNHHo9HHTp00I9+9CO9//77dguMIUVFRRowYIAkacCA\nAXrnnXdqLTfGaM+ePZo/f74mTJig559/PhplRlVRUZGuuuoqSdLFF1+skpKSmmW7du1St27d1KFD\nB7Vu3Vp9+vTRtm3bolVq1NXXV5K0fft2rVixQhMnTtRvf/vbaJQYU7p166alS5ee8Dr7VeRO1oex\nqr5642F8pKWlacaMGZKkcDgsj8dTs2zw4MGqrKw8re16Gl4ltjz33HP6/e9/X+u1Bx98UGlpadq6\ndWud7zlw4IA6duxY8++EhARVVFQ0aZ2xoq7+6tq1qzp06CBJat++/Qm3V4LBoCZNmqSbb75ZoVBI\nkydP1k9+8hNdcMEF1uqOtuP3GY/Ho3A4LLfbfcKy9u3bt5j9qS719ZUkDR8+XDfccIM6dOigadOm\n6c0339TAgQOjVW7UpaSk6JNPPjnhdfaryJ2sD2NVffXGw/ho166dpOp9dMaMGbrzzjs1Z84cffzx\nx9q3b5+mTp0qj8ejp5566pS2G3cBZMyYMRozZswpvadDhw61TrIHDx5Up06dnC4tJtXVXz//+c91\n8OBBSdV9cexBT6re2SZNmqS2bduqbdu2uvzyy7Vz584WFUA6dOhQ00eSap1QW/L+VJf6+kqSpkyZ\nUhN4Bw4cqB07dsTcATYWsF+1TPEyPj777DPdcccduvHGGzV8+HANHz5cUvUMyJNPPqnWrVuf8jab\nzS2Y+ni9XhUVFamyslIVFRX66KOP9OMf/zjaZUVN79699eabb0qS3nzzTfl8vlrLd+/erQkTJsgY\no6qqKhUVFalnz57RKDVqju2jd999t1b4Ov/887Vnzx59/fXXqqys1LZt23TJJZdEq9Soq6+vDhw4\noPT0dB06dEjGGG3evLnF7UsnY457Dij71ak7vg9j3fH1xsv4+PLLLzV16lTdc8898vv9tZa5XK7T\n/jvE3QzIqcjPz1e3bt00aNAgTZo0SRMnTpQxRpmZmWrTpk20y4uaCRMmaPbs2Zo4caLatGmjX//6\n15Jq99eoUaM0duxYtW7dWn6/X+eff36Uq7YrJSVFGzdu1Pjx4yVV3+Z75ZVXdOjQIY0dO1Zz5szR\nT3/6UxljNHbsWCUlJUW54uhpqK8yMzNrZtSuuOKKms8ftXQul0uS2K8a4Ugfxou6/ubxMD5WrFih\nr7/+WsuWLdPSpUvlcrn0xBNPqE2bNnr99ddPe7suE28REgAAxL0WcQsGAADEFgIIAACwjgACAACs\nI4AAAADrCCAAAMA6AggAALCOAAIAAKz7fycJhyYOMuVWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b8de6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(trace[burn:], varnames=['y'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "odds:\n",
      "\n",
      "  Mean             SD               MC Error         95% HPD interval\n",
      "  -------------------------------------------------------------------\n",
      "  \n",
      "  0.552            0.175            0.009            [0.246, 0.898]\n",
      "  1.260            0.441            0.015            [0.569, 2.203]\n",
      "  0.422            0.168            0.006            [0.168, 0.780]\n",
      "\n",
      "  Posterior quantiles:\n",
      "  2.5            25             50             75             97.5\n",
      "  |--------------|==============|==============|--------------|\n",
      "  \n",
      "  0.277          0.432          0.528          0.648          0.967\n",
      "  0.634          0.945          1.179          1.498          2.370\n",
      "  0.192          0.307          0.393          0.498          0.857\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pm.summary(trace[burn:], varnames=['odds'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}