ELSO Survival.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:da3baff1a9862c95f030ab6d00170f86f723aa1bdf8f9a2377c5dbb5807ff380"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu = pd.read_csv('data/flu.csv', index_col=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "(1406, 58)"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Turn into table of unique patients, with appropriate organism columns"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_only = flu.groupby('PatientID')['OrganismName'].apply(lambda s: \n",
      "                        len([x for x in s if str(x).startswith('Influenza')])==len(s)).astype(int)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coinfections = flu.groupby('PatientID')['OrganismName'].apply(lambda s: \n",
      "                        [x for x in s.dropna() if not str(x).startswith('Influenza')])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "MRSA coinfection"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mrsa = coinfections.apply(lambda x: x.count('Staphylococcus aureus, meth resist') > 0).astype(int)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Any coinfection"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "any_coinf = coinfections.apply(lambda x: int(len(x)>0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_unique = flu.drop_duplicates(cols=['PatientID']).set_index('PatientID')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_unique['flu_only'] = flu_only\n",
      "flu_unique['mrsa'] = mrsa\n",
      "flu_unique['any_coinf'] = any_coinf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Add time from admission to time on ECMO"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "admit_through_emco = flu_unique.AdmitToTimeOnHours.add(flu_unique.HoursECMO, fill_value=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See how many have no time information from admission, and drop these"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "missing = admit_through_emco.isnull()\n",
      "missing.sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "1"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete = flu_unique[missing-True]\n",
      "admit_through_emco.dropna(inplace=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Confirm no null here\n",
      "admit_through_emco.isnull().sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Time off ECMO through to event (death or discharge)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "off_ecmo_to_event = flu_complete.TimeOffToDCDateHours.add(flu_complete.TimeOffToDeathDateHours, \n",
      "                                                          fill_value=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "off_ecmo_to_event.isnull().sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "51"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Not sure why add() does not work here\n",
      "flu_complete['time_to_event'] = (admit_through_emco.values + off_ecmo_to_event.fillna(0).values)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete.time_to_event.isnull().sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Clean covariates"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete.Race.value_counts()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "W    554\n",
        "A    102\n",
        "B     83\n",
        "H     82\n",
        "O     56\n",
        "dtype: int64"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete['non_white'] = (flu_complete.Race!='W').astype(int)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete['male'] = (flu_complete.Sex=='M').astype(int)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "covariates = ['male', 'AgeYears', 'non_white', 'HoursECMO', 'pH', 'PCO2', 'PO2', 'HCO3', \n",
      "              'SaO2', 'mrsa', 'any_coinf']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Get counts of missing values"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete[covariates].isnull().sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "male           0\n",
        "AgeYears       0\n",
        "non_white      0\n",
        "HoursECMO      1\n",
        "pH            95\n",
        "PCO2          85\n",
        "PO2           81\n",
        "HCO3         167\n",
        "SaO2         147\n",
        "mrsa           0\n",
        "any_coinf      0\n",
        "dtype: int64"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to drop the missing values for some of the covariates, and not use some covariates with many missing values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "drop_subset = ['HoursECMO', 'pH', 'PCO2', 'PO2']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete = flu_complete.dropna(subset=drop_subset)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "flu_complete[['pH', 'PCO2', 'PO2', 'HCO3', \n",
      "              'SaO2']].hist(bins=25);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAENCAYAAADjW7WQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXtYVVX6x7+Hm9w8IBgyRoTKQdCIUiwJTbxko9KMDTPO\niP5EGpu8lhdCJ820Mi9IYoiYOb9f3mbKHsPRxmxGEANyJgQ0E1BBLQwBFYSD3M5l/f5gzu4A57bP\n/fJ+nofnYe/zrnXW2u8+66zz7nd9l4AxxkAQBEE4BE6WbgBBEARhPmjQJwiCcCBo0CcIgnAgaNAn\nCIJwIGjQJwiCcCBo0CcIgnAgaNAnCIJwIBxm0J8/fz6ee+45la85OTnhr3/9K3f8008/YcGCBYiK\nioK3tzdCQkLwm9/8Bvn5+T3K3bx5E3PnzkV4eDj69++PUaNGYfXq1Whra+thd/DgQYwePRp+fn7w\n9PTEiBEj8P777xu9j0S3n52cnODk5ARXV1eEhIRg0aJFaGxs5GwqKysxe/ZsjBgxAl5eXggNDcX/\n/M//oKysrEdd3333HV588UUMGzYMPj4+GDt2LDZv3gypVMrZyGQyvPHGGxg1ahSEQiEeeugh/PKX\nv8S3335rtj4Tuvmd/NmNwwz6AoEAAoFAq11hYSEiIiLwww8/YNu2bfjuu+9w+PBhPPHEE0hKSuLs\nzp49i8ceewwtLS3Izs7GxYsXsX79epw4cQKRkZG4d+8eZzto0CC89dZbOHfuHL799lssWLAAKSkp\nyMjIMElfHZ1nn30WdXV1qKqqQmpqKj755BPMmzcPAPDZZ58hKioKALBr1y58//332Lt3LwYNGoSF\nCxdydXzyyScYM2YM/Pz8cODAAZSUlGDZsmXIysrCM888g66uLgBAR0cH/v3vfyMlJQX/+c9/8Omn\nn8Ld3R1TpkzB9evXzd95B0aT38mfSjAHISkpiT333HMqXxMIBOzw4cOsvb2diUQiNn36dJV2d+/e\nZYwx1tHRwUQiEZsxY0Yfm5aWFjZo0CA2b948je0ZP348S0xM5NkLQhtJSUlsypQpPc4tWrSIOTs7\ns5qaGubn58cWLVqksqzCvw0NDczPz48tXry4j83NmzeZh4cHW79+vdo23Lt3j3l6erJdu3YZ0BOC\nD7r4nfzZjcPM9AGAaVGc+Oabb1BVVYW1a9eqfN3f3x8AUFRUhKqqKrzxxht9bPr374+lS5fib3/7\nG2QyWZ/X29racOrUKRQXFyMhIUGPXhDa6P2LztvbG3K5HJ9//jmampqwbt06leUU/j1x4gSamppU\n3gePPvooEhMTsX//frXv7+TkBIlEgoEDBxrQC4Iv2vxO/uzGoQb9/Px89O/fv88f0P2FUFFRAQAY\nMWKExnoUdiNHjlT5+ogRIyCVSnH16lXuXHNzM7y9veHj44P4+Hjs378fv/nNb4zRLaIXii93qVSK\nf/7znzh48CDGjh2LW7duQSgUYvDgwRrLV1RUaLQbMWIEampq+jy7UbB48WKEhoZi5syZhnWE4IW+\nfnc0f7pYugHmZOzYsX2+0RljEIlEOsX7DUEoFOK7775DfX09/vWvf2HBggV48OABkpOTTfq+joji\ny10qlaKrqwtTpkzBnj17sHv3bpO/9+rVq3H69GmcOXMG/fr1M/n7ET9jCr/boz8datB3d3fH0KFD\n1b6umOFfvnwZsbGxau0iIiIAAJcuXcK4ceP6vH758mW4uroiLCyMOycQCDB06FAMHToUMTExuH37\nNt5//30a9E2A4svdxcUFgwcPhotL920eERGBlpYW/PTTT3j44YfVllfY3bp1C0FBQX1ev3z5MoKD\ng+Hp6cmdY4zh1VdfxZEjR5Cbm6v2VyBhOrT5nfzZjUOFd7TN5p955hmEhoZi06ZNKl+/e/cuACA2\nNhbDhg3D5s2b+9i0tLRg165dmD17NpydndW+V1BQEFpaWni0ntAVxZd7cHAw98EHgBdeeAF+fn54\n9913VZZT+PeFF17AgAED8N577/WxuXnzJv7617/2yOSSyWR46aWX8PnnnyM/Px+RkZFG7hGhC5r8\nTv5UwoIPkc2Kqqf7ChTZO4wxdvbsWda/f382ZcoU9uWXX7KqqipWWFjI3n33XRYcHMyVOXPmDPPy\n8mLx8fEsLy+PXbt2jeXk5LDw8HA2bNgwLhOEMcbWr1/PTp8+zaqqqlhBQQF77733mFAo1JgxQOiH\nJj8zxtgnn3zC3Nzc2B/+8AeWm5vLqqur2ZkzZ1hqaip7+umnObvDhw8zV1dXlpyczAoLC9mVK1fY\noUOH2ODBg9lTTz3Furq6GGOMSaVS9tvf/pb5+/uzgoICdvv2be6vtbXV5P0lutHmd/LnzzjMoD9/\n/nytKZsKbt26xV566SX22GOPMU9PTxYcHMwSEhJYQUFBj3LXr19niYmJLCwsjHl7e7MnnniCpaam\nsgcPHvSwW7FiBQsNDWUeHh5swIABLDo6mmVnZzO5XG78jjo4mvysoLy8nP3+979nw4cPZx4eHmzo\n0KEsKSmJXbx4sYddWVkZmzlzJhsyZAgTCoXs6aefZu+99x6TSqWczY0bN5hAIGBOTk5MIBD0+Nu4\ncaNJ+kj0RRe/kz+7ETCmPo9x9+7dKCsrg1AoRHp6OgCgvb0dmZmZaGhowKBBg7Bs2TK4u7sDAE6e\nPIm8vDw4OzsjOTkZ4eHh5vm5QuiNKh8rOHHiBA4dOoS//OUv8Pb2BkA+tlXkcjnWrFkDPz8/rFmz\nBkeOHEFeXh6EQiEAYPbs2XjyyScBkI/tHY0x/YkTJ/bJRT969CiGDx+O7du3QyQS4ejRowCAW7du\n4cyZM9iyZQtWrVqFrKwsyOVy07WcMAqqfAx0x7e/++67HrnJ5GPb5eTJkwgKCuKeawkEAsTHx2Pb\ntm3Ytm0bN+CTj+0fjYN+REQEvLy8epw7f/48JkyYAACIi4tDcXExAKC4uBixsbFwcXFBQEAAAgMD\nUVVVZaJmE8ZClY8B4MCBA5g7d26Pc+Rj2+TevXsoKyvDpEmTuFx21h3a7WNLPrZ/eGfvNDc3w9fX\nFwDg4+OD5uZmAEBTUxO3ohHoXt2oLHZE2A7FxcXw8/PDo48+2uM8+dg22b9/P+bOnQsnp58/7gKB\nAKdOncKKFSuQnZ2NBw8eACAfOwIGpWxqS4E09YInwvh0dnYiJycHs2bN4s5peOxDPrZySkpKIBQK\nMWTIkB5+nDp1Knbt2oVNmzbByckJBw4cUFsH+di+4L04y8fHB/fv34evry+amprg4+MDAPDz8+uh\nLHnv3j34+fmprSc3N1eP5hKGMnnyZI2v19fX486dO3j99dcBAI2NjVizZg02bdpEPrYRlH185coV\nlJSUoKysDBKJBO3t7di1axeWLl0KAPD09MTzzz+PzMxMAPQ5thW0fY41wXvQj46ORn5+PmbOnImz\nZ89izJgx3PmdO3ciPj4ejY2NqKurQ2hoqMa6Ro0apV+r0S2BrGo1rKnL2vJ7l5aWarUJDg7GRx99\nxB0vWbIEW7duhbe3t9l9rIyh180UdVljm3r7ODExEYmJiQCA8vJyHD9+HEuXLkVTUxMGDBgAmUyG\nwsJCBAcHAzD955hvP/nY22rdfO11+RxrQmPKZkZGBioqKiAWi+Hj44NZs2Zh7NixGlM2c3NzuVQv\nhVyBKnJzc402IOhLnbgT9eJuHe1B/d0Q2N8+tDXUUVpa2meGoMrHEydO5F5funQptmzZ0iNl05Z8\n7Gio8rGCy5cv44svvsDq1auRmZmJH374AS4uLoiIiMCvf/1r7lkd+di60eRjXdA46JsSa7hZLtaK\n8frJ7syEtOmhiBrc36LtMTWG3ix8sQYfOxrkY/vHUB/brPZOYWGhRcra8nvbMsbst7HqssY2WTt8\n+8nH3lbr1sfeEBxKZVMRzpH7h+BirRhdMlp0QhCEY+FQ4R3lcA4AvDVlCDaevgGAwjumgH76mx91\nPu4tw2AsORXysflx2PAOQRC601uGwVRyKnXiTlysFeNirRh14k6T9IUwDJsd9C0ZA6WYvvmxxvi5\nNbZJFapkGEwlp1Iv7sLrJ6vw+skqLjNOGVuNu1NMn7AbVKlsHjx4EKWlpXBzc0NERARmzZrF7Spk\nCQXGOnEn9xzGEVJrjY1ChqG9vZ07p0lORSQScXYkw2B/2OxM31iLYsz93oa229j9VqWyGRUVhfT0\ndGzevBkdHR3IyckBYDkFxnpxF94tuqN29sgXY11DY/rCVPezOhkGZQyVU1GepSq+PJRf6z2LVT7u\n/boh9uPGjdNan772Cv9Yg72h0EzfwYmIiEBDQ0OPc48//jj3/xNPPIFvv/0WgPqf/sp7ARPWhSoZ\nhszMTKPJqQA9v7C667mj8jU6Ns6xoStybXamb+wYmJuzgHsApe0hlCPF9HNzczmpDXtRYHSkmH5i\nYiKys7ORlZWF5cuXY+TIkVi2bBknpwKgj5xKUVERpFIpGhoadJJh4IOtxt0ppm+HNLZLufRNoDuF\n09Fjx59//jnc3d0RExOj1kaXn/7KP10B8D7uPzSKq6+5uRn4b2qtvvUpt02f8orjS5cuGVTeFMeK\nZy/qUPgrISEBmZmZSElJ4VI2ASAoKAgTJ07E6tWr4ezsjMWLF5PKpp1Befr/HeiV/wfsM29fXX5v\nQ0MDtm7d2mO7xPz8fOTm5uLNN9+Em5sbAODYsWMAgJkzZwIANm3ahFmzZvV48KeMsXzsaHIZhmDp\ntRjkK9NjqI/1numfPn0a+fn5kEgkiIiIwPz58zUu+CBshwsXLuD48ePYuHEjN+AD+ikwEgRhXegV\n029tbUVOTg7WrVuHzZs34/bt27hw4YLaBR+mgPL0jUNGRgbefPNN1NbWYtGiRcjLy8P//u//oqOj\nA++88w5SU1Oxb98+AD1/+m/fvt1mf/o7Ukzf2rDVuLvDx/QVs7+2tjYA3bsteXl54fz589iwYQOA\n7gUfGzZswJw5c4zTUsIkLF++vM+5SZMmqbWfPn06pk+fbsomEUamq6sLGzZsgEQigZubG2JiYhAf\nH48jR44gLy8PQqEQADB79mxug3RLrMcgzIPeg/6CBQuwZMkSuLq6Ytq0aRCJRGoXfJgCytN3TBRZ\nVoD+eyA4Up4+0P15feutt9CvXz9IJBKsWbMGo0ePhkAgQHx8POLj43vYK6/HaGxsxDvvvIOdO3f2\n2GNXX/j2k4+9rdatj70h6OXFlpYW7Nu3Dzt27EBWVhauXr2KkpKSHja6/Ozns0jDGMe9v4SkUqnK\n/4HuLBFzt8/Ux/ZAY7tU4zJ/QjX9+nV/OXZ0dEAmk8HV1RWA6v2PDZViIKwbvWb6VVVVEIlECAwM\nBADExMSgoqJC7YIPdSh/u/FdpKBP+e4Z4s8LR1xcXFT+D3T/UomKUF2/choi3/YrylqivKGLOuwF\nVf6zZD3GrksVcrkcq1evRk1NDebPn4+BAwcCAE6dOoW8vDyEhYVh3rx58PLyMqkUA99+8rG31br1\nsTcEvWb64eHhqK6uRmtrKyQSCcrKyhAVFaV2wQdBEJbFyckJaWlp+OCDD/DVV1/hxo0bmDp1Knbt\n2oVNmzbByckJBw4cUFte0y935V+R2mQYLl26xOsXKl97Rzg2FL1m+p6enkhISEBaWhq6uroQFRWF\nkSNHIjQ0VOWCD1NAMX1CXxwtpq9MQEAAnnzySZSXl2PGjBkAuj/Pzz//PDIzMwHwl2LgI8OwaNEi\ntWUNtdcWDTDU3lqODf3FrneeflxcHOLi4nqc8/DwQGpqqkENIsyLKpVNY22wQVgHLS0tcHZ2hpeX\nF8RiMS5cuIDk5GQuFCuTyVBYWIjg4GAAtB7D3iHtHTO/t7Xl6atS2TTVBhvWgqPl6d+/fx9vv/02\nXn/9dXzwwQeIj49HZGQkDh06hJSUFKxduxYymQxJSUkATLsew1Zz6R0+T5+wH1SpbKpbb2FqlU3F\nHsaA/umYRF+Cg4OxdevWPueXLl2qtgytx7BfbHbQN3UMVFM+uL3H9C21wYZi1yUA2PmCiPsCMPYG\n9o4c07c0tppL7/B5+o4A5YN3Y8wNNrRlJihnfihf/y7Zz7nkyusp7HEthaUzOwj7x2Zn+ubMazXm\nexvabnP021QbbGjLTOid+aEK5fUUmtZSaDo2dK2EunOG1GfI2g/lY1WZHepkGCzxwN5Wc+kdPk+f\nsG8stcEGYRoUMgxpaWnYsGEDzpw5g9u3b9v9A3tCNTY76FOevnFQqGzevn0bixYtwpkzZ5CQkICr\nV68iJSUF165dQ0JCAgD7Udl0xJi+sgyDXC6Hq6srzp8/jwkTJgDofmBfXFwMwLQyDLYad7enmL7N\nhncI46BKZROA2vUWlNVhm6iSYbDUA3vCstjsTJ/y9Al9cbQ8fUC1DIMyhjyw5yPDkJ2dzevhNB97\nxf+msO9dxpL2hkIzfYJwIJRlGIz1wJ6PDENkZCSvh9987U11rBhsrcHeUBkGvWf6HR0d2LVrF1JT\nU7FixQpcu3YN7e3t2LZtG1JSUpCWloaOjg6DGqcJiukT+uJoMf2WlhY8ePAAADgZhuDgYIs8sLfV\nuDvF9AHs27cPI0aMwNKlSyGTydDZ2cllA6SmpuLYsWM4evSoXeycpbxQC6DVooRtcf/+fS4Dx9fX\nl5NhUCeQqPzA3tnZ2WYf2BOq0Wum39bWhsrKSm5bPWdnZ3h6eqrNBjAF5oxtKy8Uev1kFa7e0pxL\nrgmK6VseR4vpK2QY0tLSsHbtWu5zqxBI3L59O15//XUuRx/ofmCfnp6Obdu2ISIiwmhtsVV9HIfX\n3mloaIBQKERWVhauX78OkUiE5ORks26XqCvKei7GXs5v75w+fRr5+fmQSCSIiIjA/PnzNS7oMTfG\n2DqRIBwNvQZ9mUyG6upq/OY3v8HLL7+MvXv34ty5cz1sdN0uUd2DDG3HupbvPzSK03NZG/dIj/fX\ntF1i72NlfHx8eLfXGo49PT3V9qk3ra2tyMnJQXp6Otzc3LB161ZcuHAB33//vdWE8Brbpdh4ujsL\nJW16qM6DvqPF9O/evYusrCw0NzdDKBRysuiW2BjdVuPuDh/T9/f3h7e3N6KjowEAsbGxOHv2LHx9\nfc22XaKux8qx+N5bImraLrH3saa228oxn6f+bm5uALpDeQDQ2dkJLy8vtQqchPXi4uKCpKQkhISE\noKWlBatWrUJoaKhFNkYnLI9eXvT19UVgYCCuXbsGuVyO0tJSREZGYvTo0WbbLtGSsW1Dwla2EtN3\nc3PDggULsGTJEvzpT3/C8OHDIRKJrDKExxdHi+n7+voiJCQEACAUCjFs2DBusZW5N0a31bi7w8f0\nAWDJkiXIyspCS0sLgoODMWfOHDDGzLZdImFaWlpasG/fPuzYsQNeXl54//33UVJS0sPG2CE8Xb5A\nNIXkdA158bVXd3zp0iWDylsihFdXV4dbt24hLCwMV65cMfvG6ITl0XvQHzx4MDZt2tTnvLm2S7Rk\nvrq2sJUmbCVPv6qqCiKRCIGBgQCAmJgYVFRUqF3Qow4+ITy+Kpu9Q3DmDpnx3e9V07Eqvxo7hNfR\n0YGMjAwkJSXB3d0dU6dOxW9/+1u0t7fj4MGDOHDgQJ8+KTBWyqatxt3tKaZPQTpCJeHh4aiurkZr\nayskEgnKysoQFRWldkEPYd1IpVKkp6dj/PjxnM98fHwgEAi4jdEVIRy+K3L5yDDQseX3TLBZGQbl\nsIG5aW5uBgb316usoe02V789PT2RkJCAtLQ0dHV1ISoqCiNHjlS7oMeWMNY1NKYvTOlXxhj27NmD\noKAgzJgxgzvf1NSEAQMGGLwxunK7tckw9D6n7dcLH3tV19BY9gpba7A3VIbBZgd9wvQoUvuUUSzo\nIWyHK1euoKCgAMHBwZzvZs+ejaKiIty8eRMuLi6IiIhQuTE6rci1P2x20KeYPqEvjpanHx4ejk8/\n/bTPeUVOvipMJaFtq3F3iukTBEEQNonNDvqUp0/oi6Pl6VsTtppLT3n6BEHYDOpkGCyxMTpheWx2\n0KeYPqEvjhbTVyfDkJ+fr1JHyZQyDLYad6eYvoOjUHdU/NWJOy3dJJNg6Y1yCOOgTobBEhujE5bH\nZgd9S8ZAG8QdPfT1FdLNumBLMX3FRjnbtm3D9u3b8fDDD3Mb5Wzfvh0ikQhHjx41W3uMhSPH9JVl\nGDRtjO7v78+VMaYMg63G3e0ppm/QoC+Xy5GamootW7YAAM0C7Qhr2CiHD8q/vuz1l5eh9JZhUMaQ\njdEJ28KgQf/kyZMICgribghzzgItGdvWJrusCVuJ6StvlLNq1Srs2bMHnZ2dVquyqby7mbZfXo4W\n0wfUyzDcv38fAAzaGF2dDIObswBfV9zC1xW3enwZ85Ed4GM/btw4XrIGfOyVV85a2t5Q9B697t27\nh7KyMrz44ov44osvAIC01u0IS2yUY6jKpqr31vR+9nisSmVTnQyDQkdp5syZfTZGN4YMQ2O7FO8W\nKSQZ7iBteqjZRfHs8dhiMgz79+/H3Llz0d7ezp0z5yzQkto72gYbTdiK9o4lNsoxVGVT03srH/PV\nRVF33PucIfXx1Y1Rd6xqQFAlw5CYmIiEhASzb4zO9/7lY2+rdetjbwh6DfolJSUQCoUYMmQILl++\nrNKGYoC2jfJGOcOGDeM2yhk4cKDK2SFhvaiTYQDUS6GbSoaBsDx6DfpXrlxBSUkJysrKIJFIuEUe\nfLXW+fz0732sa/n+Q6M4e0374PLZI1fVDNMaftprO+azRy5gvxvlOGJM31qw1Vx6e8rT12vQT0xM\nRGJiIgCgvLwcx48fx7Jly3Do0CFes0Bj/TTWdGyqPXJ7Yw2xPm3HfGOBlt4ohyAI42OUPH1FKCch\nIQFXr15FSkoKrl27hoSEBGNUrxJL5ukbGtM3BEfRaDEljpanv3v3brz88stYtWoVd+7IkSNYuHAh\nUlNTkZqairKyMu61kydPIiUlBatXr0ZlZaVR22KrufT2lKdvsAzDiBEjMGLECACktU4Q1sjEiRMx\nbdo07Nq1izsnEAgQHx+P+Pj4HramlGAgrAOb9STl6RP64mgx/YiICHh5efU5zxjrc87UEgy2Gnd3\n+Ji+NVMn7uyxOKdLJrdgawjCejl16hTy8vIQFhaGefPmwcvLC01NTRCJRJyNMSUYCOvAZmf66mJg\n9eKuHro4XbK+sxlDoZi+beNoMX1VTJ06Fbt27cKmTZvg5OSEAwcOqLXVln6tbkWuKrKzs3mtyOVj\nr/jfFPa9y1jS3lDsbqZPGBe5XI41a9bAz88Pa9as0ajBTtgOinRqT09PPP/888jMzATAX4IBUL8i\nVxWRkZG8svb42pvqWFXKuKXsDV2Ra7MzfWuN6deJOzXKLttaTN+S+kqmwtFi+qpoamoC0C23UVhY\niODgYADdEgxFRUWQSqVoaGjQKsHAF1uNu1NMn1CLIrykIG16KAL797Ngi/SH9JXsg4yMDFRUVKCl\npQWLFi3C7373O5SXl+PmzZtwcXFBREQEkpKSAJhWgoGwDmx20Lcm7R2FrC+g/cGxrWjvAJbXVzIV\nxrqGxvSFKf26fPnyPucUktmqMKUEg63q49iT9o7NhnesCWVZX1M8OLYEyvpKqlL7AN1VNpX/13Rs\nCpVNPu/P9/jSpUsmrV+fY4LQhs3O9K01pq8NW4npG0tfic9DOHOpbBrreNGiRUarT5VfzSG1YW5s\nNe5OMX3C7jGWvpIlUA63DervZrPPVIzF7t27UVZWBqFQiPT0dADQmIV18uRJ5OXlwdnZGcnJyQgP\nD7dk8wkjY7PhHdLeMS+W0FfSF227aDlanv7EiRPxxhtv9DinLgtLWYZh1apVyMrKglxuvAWOtqqP\n4/DaO3fv3kVWVhaam5shFAoRFxeHuLg4yuG2U0hfybaJiIhAQ0NDj3PqsrDUyTCEhYVZoOWEKdBr\n0HdxcUFSUhJCQkLQ0tKCVatWITQ0FPn5+Rg+fDhSU1Nx7NgxHD161GTpfBTTtw+UZTPMJZlBefrq\ns7BMLcNgq3F3e4rp6xXe8fX1RUhICABAKBRi2LBhaGxsxPnz5zFhwgQA3bOH4uJiozWUsE+UZTPs\nJfPJ1tCWhWVMGQZryHCy9WNDMfhBbl1dHW7duoWwsDDaI1cHbClP315xtDx9VajLwtJHhuGh4U/i\nyp0HcHES6JSBxSeji4+9qmtoLHuFrTXYW1SGoaOjAxkZGUhKSuoTuzd2DjffHGlVGGu7RH2whxmC\nLaLI5FEliaEsmdH7NXsnOjoa+fn5ANAjC0sfGYZlf7+KZX+/in3f1pq62YQR0HumL5VKkZ6ejvHj\nx3M3jClzuHXNkVbeHrE3xtouUZ+YvqWEopSP+cwQ7OVhfWO7FBtP3wDwsySG4tooS2boI5dhKzF9\nhQyDWCzGokWLMGvWLCQkJKjc69jUMgy2Gne3p5i+XoM+Ywx79uxBUFAQZsyYwZ1XzB6sOYeb0A1r\neFhPGAdVMgyA+r2OTSnDQFgevcI7V65cQUFBAb7//ntuj80LFy7QHrkq6B1esJU8fXt+WO9oefrW\nhK3m0jt8nn54eDg+/fRTla+ZO4fb2nfK6h1esEUs9bCeMD1LliyBh4cHnJyc4OzsjM2bN9tcCI/g\nh83KMKiKywLAW1OGmPy9HSlP39CH9dYI5en3ZMOGDfD29uaOFat1TRHCs9W4uz3F9G1WhoEwPZoe\n1gPQ6WG9towi5V8KuoTN+Khs9rGRS/F1xS0u3Nba3jNjxxoyrCyRodVbRdUeQniEemx2pm+Lefpu\nzgJ8XXGLGyiF7i5o6eiuS1dhMHP121gP67VlaHVnW3XndevyC4qPymZvm1apAO8W3eHer/evQr4Z\nUb3PGZJhxTfHXN0x3xxugUCAt99+GwKBAFOnTsWUKVNMGsLje//ysbfVuvWxNwSbHfRtkcZ2aZ9B\np3c6obWgeFgfHBzMPadJTExUm+pH2CbvvPMOBgwYgFu3bmHz5s14+OGHe7xuqyE8Qj02O+jbqvaO\noZir39b0sN5asYeY/oABAwB05+c/9dRTqKqq4r3eRpnWBw+02ijPahUhKXXHfOzHjRuntT5D7K3l\n2NPTs++GI3yIAAAgAElEQVRF5YFNDvqWEOkiCHujs7MTcrkcHh4eaGlpQVlZGZKTkw1ab+Pt5aXx\ndWtYpGjrxxaVYbAUlhbpMrZEA9Azn1+TLICj5HObG01yDeqw9Tz95uZmrF+/Hq+//joyMjIwY8YM\nREVFmWy9jeKZlqmuMeXp64ZNzvTtEeV8fsD6Yvz2jiq5BnsnICAAaWlpfc6bas8E5WdajnKNrREa\n9PXAEWL6xqRe3IkbTR0AgCCffgjysY+FPvYQ07d2rCWX3p7y9GnQJ0zO/Q4p1v/zOgBg3aQQuxn0\nCcIWMXpMv7y8HKtXr0ZKSgq+/PJLY1dvFZgipq8r1hDTN9THypLG1vggXtf4vq3H9DVhLZ9ja4m7\nU0xfDXK5HNnZ2XjzzTfh5+eHP//5z4iMjERQUJBB9Vq7vo656X09dF3YZQyM4WNl6QxzyGbwxRHj\n+8qY6nNMWAdGHfSrqqoQGBiIgIAAAEBsbCzOnz+v083S0NrFDebuLk6QylmPtMy1X13nbC09UFg6\npn+xVtxDb0jTwNT7C8JQDPGxLaKY9QN9v1ztNaZvDh9ruq7KWEvcnWL6amhsbIS/vz937Ofnh6qq\nKg0lfuZ0VSM+Pn8bAJDybDAGebtZ9WzQVlCeVW8ZZXh9hvjYFlGe9e98QcR9gSpLaOgjp2HNmMPH\n6n5NWfJXrKNgNQ9yY4J98Mh/H/AN8XPHvQcSC7dIPeaI6SvPhICfB5bm5mZ4ePdXa6s8AAHWEQrz\ncXfBm5O7v7it3beaUB6olCU01MlpKA9g2r4kFKtOlcvYyoCn8K2wn7Ne5ZXv396/6tN+GaL1Gip/\nNnx8fNR+Cfe+tlUXi7XOsBVl+NQtdHdBTf09biWzOj8qlxG03cfjoY/wum76ImC9JfYM4OrVq/js\ns8+wdu1aAEBOTg4EAgFmzpzZxzY3N9dYb0vwYPLkyQaVJx9bP+Rj+8cQHxt1pj9s2DDU1dWhoaEB\nfn5++Oabb/Daa6+ptDX0xiQsA/nY/iEf2zdGnekD3aleH3/8MWQyGSZPnkx7bdoh5GP7h3xsvxh9\n0CcIgiCsF5sUXCMIgiD0gwZ9giAIB8IiKZvl5eXYv38/Fy+cNm2aRvvdu3ejrKwMQqEQ6enpAID2\n9nZkZmaioaGB28Gp98bdAHD37l1kZWWhubkZQqEQcXFxiIuL06l8V1cXNmzYAIlEAjc3N8TExCA+\nPl7n91Ygl8uxZs0a+Pn5Yc2aNTqXX7JkCTw8PODk5ARnZ2ds3ryZ13t3dHRg3759+PHHHyGRSLB4\n8WIEBQXxaru+8PWxMob0m++9cvLkSeTl5cHZ2RnJyckIDw/XWNeRI0eQl5cHoVAIAJg9ezaefPJJ\njXXpcw/yrUufdhlKbx/fuHGD17X/7LPPcPz4ccjlcjz88MOYPn262uvi5OSEDRs24N69e3jw4AG8\nvb2xfPlyPProoxqvYW5uLurr6xESEoJ3331XY3uSk5PR1dUFgUCAhx56CDt27FBr39HRgY0bN+LH\nH38EAMyfPx/jxo1TaVtbW4u3334bYrEYAoEAAoEAs2fPxsSJE1Xanz59Gjk5Obh//z48PDyQkpKi\ntZ+8/MvMjEwmY0uXLmX19fVMIpGwlJQUVlNTo7FMeXk5u379Olu5ciV37uDBg+zYsWOMMcZycnLY\noUOHVJZtampiN27cYIwx1tzczBYsWMBqamp0Lt/R0cEYY6yrq4utXLmS1dbW6lxWwYkTJ9jOnTvZ\nli1beLV98eLFTCwW9zjH570zMzNZbm4uY4wxqVTKHjx4wLvt+qCPj5UxpN987pWamhqWkpLCJBIJ\nq6+vZ0uXLmUymUxjXUeOHGEnTpzo876a6uJ7D+pTlz7tMgRVPj579iyva79ixQp27do1Vl9fzxYv\nXqz1ulRXV7OUlBTW1tbGXn31VbZw4UJ24MABjdfw2LFjbOvWrWzOnDlMJpNpbM8f/vAH1tTU1OM6\nqbPfunUrW7hwIZNIJKy2tpYtXrxYa1skEgm7ffs2+/3vf8/q6+tV1i0Wi9krr7zCVq5cyTo7O9lb\nb73FXn75ZZ3q1tW/Zg/vKC/xdnFx4ZZ4ayIiIgJevXbkOX/+PCZMmAAAiIuLQ3Fxscqyvr6+CAkJ\nAQAIhUIMGzYMjY2NOpfv1697UUVHRwfkcjlcXV11LgsA9+7dQ1lZGSZNmgT232fmfMqzXs/ZdS3b\n1taGyspKTJo0CQDg7OwMT09PXu+tL/r4uDf69pvPvVJcXIzY2Fi4uLggICAAgYGBPVaeqqpLVdu0\n1cX3HtSnLn3aZQiqfNzY2Mjr2j/77LMIDQ1FQEAABg8ejMDAQI3X5eLFi4iNjYVUKoWTkxMCAgJw\n7tw5tfU/+eST+O677zB9+nS4u7ujqqpKY3vc3d3h7Ozc4zqpsm9ra0NFRQWef/55uLi44Be/+AUG\nDx6ssS0KH9TX18PDwwP3799XWbebmxs6OzsxevRoyOXdCysHDhyoU926+tfs4R1jLfFubm6Gr68v\nAMDHxwfNzc1ay9TV1eHWrVsICwvTubxcLsfq1atRU1OD+fPnY+DAgbzee//+/Zg7dy7a29t5t10g\nEODtt9+GQCDA1KlTMWXKFJ3LNjQ0QCgUIisrC9evX4dIJEJycrJe140vhvrYkH6rQl3ZpqYmiEQi\nzs7f358bQDVx6tQp5OXlISwsDPPmzYOXl5fOdelyD+pT15UrVwxqF1909bGuffTw8MDNmzc1XhfF\nF8Inn3yC+fPn4+bNm7h+/bra+quqqvDyyy+jvb0drq6uaGxs1NgeZ2dn7r7r16+fWvuGhga4uLig\ntLQUBQUFEIlE8PX1RWVlpda+FhUV4eGHH1Zbt5ubG8LDw3H8+HGcOnUK06ZNQ3NzM2pqaox2D9vF\ng1yBQKDVpqOjAxkZGUhKSuoTC9ZU3snJCWlpafjggw/w1Vdf4caNGz1e11S2pKQEQqEQQ4YMUTkL\n01b+nXfeQVpaGl599VXk5OSgoqJC57IymQzV1dV4+umnsXnzZkilUpw7d07n8pbEkH5rQ1tZba9P\nnToVu3btwqZNm+Dk5IQDBw7oXJe+96AudRnSLnOh7n07Ojpw+fJlxMXFabwuAoEAiYmJ3GdRLBar\ntW1oaICHh0ePz17v9+99/Nvf/pa773744Qf89NNPKu1lMhmam5sRGhrKfbYaGhq09lUqlaKkpAQB\nAQFq29LS0oLLly9jzpw5yMrKwtWrV3H37l2tdfN53eyDvp+fH+7du8cd37t3D35+frzr8fHxwf37\n9wF0f9spdC5UIZVKkZ6ejvHjx3ObPPMpD3RvLffkk0+ivLxc57JXrlxBSUkJlixZgp07d+Ly5cvI\nzMzUufyAAQMAAEFBQXjqqadQVVWlc1l/f394e3sjOjoabm5uiI2NxYULF+Dr68ur3/pgqI8N6bcq\n1JXVp50+Pj4QCATw9PTE888/z81utdXF5x7Uty592qUvutarrY+KvvTv3x9PPfWUTmUUn8WffvoJ\n/fv3V2nb0dGB6upq7rPX1NSE3NxcjXV3dPx3d7egIHh7e6OlpUWlvb+/P/r16wehUMh9tm7fvg2h\nUKix3WVlZRg6dCjEYjH8/PxU1l1VVYWAgADIZDL0798fMTExuH37ttp+6uNfsw/6yku8pVIpvvnm\nG0RHR/OuJzo6Gvn5+QCAs2fPcjd/bxhj2LNnD4KCgjBjxgxe5VtaWvDgwQMAgFgsxoULFxAcHKzz\neycmJiI7OxtZWVlYvnw5Ro4ciWXLlulUvrOzkwsJtbS0oKysjNd7+/r6IjAwENeuXYNcLkdpaSki\nIyMxevRoncobgiE+NrTfqlBXNjo6GkVFRdxMra6uDqGhoRrrampqAtA92yssLERwcLDWuvjeg/rU\npU+7DEFXH2vqY2FhIXbv3g0/Pz/IZDKuXarKtLS0YOTIkSgqKkJTUxPOnz+PtrY2xMTEqKx/wYIF\neOihh7Bz504kJSXBxcUFf/7zn9W25/HHH0dBQQGkUimqq6vR1NSEJ598UqW94rOVl5eHrq4uFBUV\nQSaTYezYsRr9WVBQgMjISM4HquoODw/HgwcP8PXXX6O9vR3//ve/IZFI1PZTH/9aZEUu3yXeGRkZ\nqKiogFgsho+PD2bNmoWxY8fqlMJXWVmJt956C8HBwdzPnsTERAwfPlxr+R9//BFZWVmQy+Xw9fVF\nTEwMJk2axDtlU9HnEydOYPXq1TqVb2ho4DatVnzjP/fcc7zeu7a2FllZWWhpaUFwcDCWLVsGxpjZ\nUjb1WcZvaL/53iuKtD5FultERESfulpaWuDr64vf/e53KC8vx82bN+Hi4oKIiAj8+te/5mKt6urS\n5x7kU9fs2bNRVFTEu12G0tvHV69e5XXt/+///g9ffvklXF1d4e/vD3d3d7XXpaGhAVlZWbh//z7a\n29vh5eWFV199FSEhIVqvoUQigY+PD9555x2191FDQwPefPNNtLa2wtnZGVOmTMG8efPU2tfW1uLd\nd9/F/fv34ebmhhUrViAsLExtW/7+97/jb3/7G4KCgvDHP/4RERERauvOz8/H0aNH0djYCE9PTyxf\nvhxDhw7V6x5WiW4JWvZBW1sbW7duHQsNDWUeHh7Mz8+PjRkzhn3wwQe867pz5w575ZVXWGRkJPP2\n9maPPfYYe+WVV9idO3d62P3lL39hcXFxbODAgax///5s9OjR7PDhw8bqEmECkpKS2JQpU1S+JhAI\nyH+ETWM1evrmYNGiRfjiiy+wZcsWLvWrtLQUNTU1vOqprKzk0s3S0tIgEolw7do1bNy4ESNHjsTX\nX3+N4cOHAwDOnDmDF198Edu3b4eLiwv+9re/Ye7cuXBxccGsWbNM0U3CQBQLaAjCLrH0t4458fHx\nYRs3btRoU1JSwn75y1+ygIAA5u3tzcaMGcNOnTrFvS6Xy9mzzz7LoqKi+iyCkEqlLDIykk2YMEHj\nezz11FMsISFB734QhjFhwgS2YMEC9t5777HQ0FD2i1/8gq1YsYJJJBLGWPdM/7nnnlNZlmb6hK1j\nFymbujJ48GDk5uZyD71UIRaLMXv2bOTn5yM3NxejRo3Cr371K1y7dg0AcOPGDRQUFCA1NRVOTj0v\nn7OzM1JTU/H111/jhx9+UPsecrkcDz30kHE6RejF559/jps3b+KLL77A3r17cfjwYaSmpnKvMxKf\nJewVS3/rmJOioiL26KOPMmdnZ/b444+zP/3pT9zSZnVIJBIWFBTENm3axBhj7IsvvmACgYBduHBB\npX1JSQkTCATs5MmTKl/fs2cP8/DwYBUVFYZ1htCbCRMmMA8PD05igzHG1q5dy9zd3dmDBw9YUlIS\nc3FxYd7e3n3+aKZP2DoONdN/5plnUF1djYKCAiQlJaGiogIvvvgifvWrX3E2d+7cweLFixEREYEB\nAwZgwIABqKur44SVDOHvf/87VqxYgY8++shooleEfsTGxnISGwDw3HPPobOzE9XV1QCAsWPH4uLF\niz3+Lly4YKnmEoTRcKgHuUB3CCYmJgYxMTFYuXIlMjIysHLlShQUFGD8+PGYP38+Ll26hN27d2PI\nkCFwd3fHr371K3R1dW9grEiHunTpEqKiovrUf/nyZQDAiBEjepz/5JNPkJycjH379mHOnDkm7iWh\nDaYlfOPu7o6hQ4eaqTUEYT4caqavCoVuRWtrKwCgoKAACxcuRHx8PEaOHAlXV1dUVlZy9kOHDsW4\nceOQlpYGmUzWoy6pVIq0tDRMmDABjz76KHf+o48+QnJyMg4cOEADvpXwzTffoLOzkzv+17/+BXd3\ndwwbNsyCrSII0+NQg/6ECRPw4Ycf4vz587h69SoOHTqE1NRUDB48GM888wwAYPjw4fjss89QWFiI\nCxcuYOnSpX2W/H/44Yeora3F+PHjcerUKVRXV+PUqVN49tlnUV9fjw8//JCz3bFjBxYvXoydO3di\n/PjxqKurQ11dnVFErwj98fDwwGuvvYbKykr84x//wL59+7Bw4UJ4enpaumkEYVIcKrwzffp0HD58\nGOvXr0dLSwsCAgIwYcIErFu3jhvYDxw4gK1bt+LFF1/EI488gtTU1B4KmUB3iKe8vBzr1q1DSkoK\nbt68iZCQEIwfPx4nTpzooT74wQcfQC6XY+HChVi4cCF3Pi4uDnl5eebpONEDgUCAhIQEPPLII5g+\nfTra29sxe/ZsbNu2jXud8vQJe0WjDIOqnYMOHjyI0tJSuLm5ISIiArNmzeJmR6baoYcwHY7o44kT\nJ0IkEmHv3r2WbgpBmB2N4Z2JEyfijTfe6HEuKioK6enp2Lx5Mzo6OpCTkwMAuHXrFs6cOYMtW7Zg\n1apVnGYNYd04oo8ZY5SHTzgsGgd9VTsHPf7443BycoKTkxOeeOIJTtbTVDv0EKbFEX1M4RvCkTEo\npp+bm8ttx2eqHXoIy2KPPj5z5oylm0AQFkPv7J3PP/8c7u7uiImJUWtDsynbhnxMEPaHXjP9/Px8\nlJWV4c033+TO8d3BJTc3V5+3Jgxk8uTJOtmRj20XXX1MOCa8B/0LFy7g+PHj2LhxI9zc3Ljz0dHR\n2LlzJ+Lj49HY2KjTDi6jRo3i32IVFBYWYty4cVZVlzW2qbS0VCc7S/mYbz/52Ntq3XztdfUx4bho\nTNlUtXPQsWPHIJVK4e3tDQAICwvDggULAPDbwUWhYEmYj9LS0j6zQPKxfaHKxwShjEW2SwRoQLAE\n5h4QyMfmhwZ9Qht2IcNQWFhodXVZY5usHb795GNv6brrxJ24WCvGxVoxvqvit1Obo/ifMA8OJcNA\nEKakTtyJenG3Guug/m4I7P+zdHO9uAuvn+xe07AuljbQISyHXcz0jfXA1Jh1WWObrB2+/eRjb466\nFQP76yeruMFfFb0F/IzZFoLQhl0M+gRBEIRu2MWgb43xc2tsk7Vj6bi7uepubm42WVsIQhsU07dB\nlGPHQN/4MUEQhDrsYtC3xvi5Kduk/FAQANKmh9rFoG/rMX1doZg+YUnsYtB3dNycBbhYKwZAs36C\nIDSjMaa/e/duvPzyy1i1ahV3rr29Hdu2bUNKSgrS0tLQ0dHBvXby5EmkpKRg9erVPfaVNTXWGD83\nZ5sa26U6ZY2owpp8bE1xd13t68Sd+LriFi7WitEl021vAYrpE5aE9yYqR48exfDhw7F9+3aIRCIc\nPXoUgP1ssOFokI8No17chXeL7uD1k1XoktHGLIT1w3sTlfPnz2PChAkAuvd5LS4uBmDZDTYcLaZv\nTKzJx9YUdzflNaeYPmFJeKdsNjc3w9fXF0D3zav4qdrU1NRjQ3Bb3WCDIB8ThD1j0INcbRtoaHtd\nWTJWEbfU51g55mlofb3r1Le+7OxsREZGmqx/2lBVn2Jzcz6Y08d8rxkfe773iK72cv8Qzk4qlfbp\nu7K9gubmZmBwf5Pck/r4mHAseA/6Pj4+uH//Pnx9fdHU1MT9VOW7wQbQ88PQ+4NhqWN1H1S+x8qD\nkSnaqw1V5XXVWreUj/leM1NfY12Ou7Om7gAAXFxcNNr3xhT3JOnpE9rgHd6Jjo5Gfn4+AODs2bMY\nM2YMd76oqAhSqRQNDQ06bbBhLKwxfm6NbdIVS/nYnmL6ijRaxZ9yZg/F9AlLonGmr9hgQywWY9Gi\nRZg1axYSEhKQmZmJlJQUDBo0CMuWLQMABAUFYeLEiVi9ejWcnZ2xePFi2j/VBiAfm4bGdik2nr7B\nHb81ZYgFW0MQP6Nx0F++fLnK86mpqSrPT58+HdOnTze8VTyxxq0JrbFNqrAmH9vyloZ8UMT0raEt\nhONhF4JrBEEQhG7YxaBvjfFza2yTtWNPMX1NUEyfsCR2MegTBEEQumEXg76ja+/YC7aovaMPpL1D\nWBK7GPQJgiAI3bALaWVrjJ8bu03KG6foquZoa1BMXzUU0yeMiV0M+o6A8sYplPNNEIS+6B3eOX36\nNNatW4fVq1fj448/BqBZh92UWGP83BrbxBdz+5hi+qqhmD5hTPQa9FtbW5GTk4N169Zh8+bNuH37\nNi5cuKBWh52wPcjH6qkTd6qUVyAIW0CvQd/NzQ0A0NbWhq6uLnR2dsLLy0utDrupcYSYvq701nyp\nE3fq9Z6W8LGtxPQVoTZ9N06hmD5hSfSK6bu5uWHBggVYsmQJXF1dMW3aNIhEIrU67IT56K35ou+m\n6eRjgrBP9Br0W1pasG/fPuzYsQNeXl54//33UVJS0sNGFyEu0tPXvX/Kuu36wFdr3RI+thU9/f5D\no7jXlDX0e+vpq4P09AlLotegX1VVBZFIhMDAQABATEwMKioq1Oqwq8PSWuiqjq1VT7/7Wt6BvvDV\nWreEj21FT79bQ78bZQ393nr62iA9fcIS6BXTDw8PR3V1NVpbWyGRSFBWVoaoqCi1OuymhmL6xscS\nPraVmL6hUEyfsCR6zfQ9PT2RkJCAtLQ0dHV1ISoqCiNHjkRoaKhKHXbC9iAfE4R9ovfirLi4OMTF\nxfU45+HhoVaH3ZRYo3a9sdukHEc2F+b2Menpm78thONB2jsEQRAOhF0M+tYYP7fGNlk7FNNXjaP4\nnzAPpL1j5ygWaxEEQQB2MtO3Rp0ba2lTY7uUWz1q7TiK9g7kUl4rpkl7hzAmNNMnCDNzv0OGTfk3\nAei/Ypog9MUuZvrWGD+3xjZZO44S0zd0ERdBGALN9K0U5U1TAPvdOIUgCPOi96Df0dGBffv24ccf\nf4REIsHixYsRFBSEzMxMNDQ0cAt33N3djdleldhjnr7ypimAZTZOMbePHSVPX1eNHnO0hXA89B70\n9+3bhxEjRmDp0qWQyWTo7OzktNZTU1Nx7NgxHD16FHPmzDFmewkzQj4mCPtDr5h+W1sbKisrMWnS\nJACAs7MzPD09SU/fBPVYCkv4mGL6qrH1e4mwLvSa6Tc0NEAoFCIrKwvXr1+HSCRCcnIyaa3bEeRj\ngrBP9Jrpy2QyVFdX4+mnn8bmzZshlUpx7ty5Hja6aq0r/6/vseJ/Y9TXu05968vOzjaofO/BlG8c\n2FAs4WO+14yPPd97RBd7QD89fWU7ZT8b454kCG3oNdP39/eHt7c3oqOjAQCxsbE4e/YsfH19SU//\nvxiq9d5bP59vSMBQLOFjR9PTB3pKMpCePmEO9Jrp+/r6IjAwENeuXYNcLkdpaSkiIyMxevRo0tM3\ncj2WwhI+ppi+amz9XiKsC72nj0uWLEFWVhZaWloQHByMOXPmgDFGWut2BPmYIOwPvQf9wYMHY9Om\nTX3Ok56+ceuxJOb2MeXpm78thONhFzIMBEEQhG7YxaBPMX37gGL6qqF7iTAmdjHoEwRBELphF4O+\ntWjXm6IeR8JR9PT1iekThLGwi0GfIAiC0A27GPQppm8fUExfNXQvEcbELgZ9giAIQjcMGvTlcjlS\nU1OxZcsWAEB7ezu2bduGlJQUpKWloaOjwyiN1AbF9E2HOX1MMX3V2Mu9RFgHBg36J0+eRFBQECe8\npdBa3759O0QiEY4ePWqURhKWg3xMEPaF3oP+vXv3UFZWhkmTJoExBgCkp8+znjpxJy7WinGxVow6\ncadR3teYmNvH1hzTV/aVoVtXUkyfsCR6yzDs378fc+fORXt7O3eOtNb5obwl4s4XRFa3Jy75+GeU\nfWWJrSsJwljoNeiXlJRAKBRiyJAhuHz5skobXbXW1cnH8jlWjnkaWl/vOvWtLzs7u4f0ryp7uX8I\n954N4g5syq/hjtfGPdLjWplbT98SPtblmulrz/ceUWWvCkP19CGXclLNgrb7aKn7waB70tPTU6c2\nEI6LXoP+lStXUFJSgrKyMkgkErS3tyMzMxM+Pj6kp/9fdNF67/6wd2vm9/7Jr+3Y1FjCx9aup68K\nQ/X0W6UCbDzd/QsibXooxo37+cue9PQJU6DXSJKYmIjExEQAQHl5OY4fP45ly5bh0KFDyM/Px8yZ\nM0lP38bjsJbwsTXF9EOjxvTYLMWY4TaK6ROWxCh5+oqf+QkJCbh69SpSUlJw7do1JCQkGKN6wgpw\nNB8rYviKvy4Zs3STCMIoGBwzGDFiBEaMGAEA8PDwID19I9djDZjLx9akp2/KB9Skp09YElqRSxAE\n4UDYxaBPMX37wJpi+toeUBsCxfQJS2IXgz5BEAShG3Yx6JP2jn1gTdo71hbTJwhjYReDPkEQBKEb\ndjHoU0zfPqCYvmroXiKMiV0M+gRBEIRu6JWnf/fuXWRlZaG5uRlCoRBxcXGIi4vjluo3NDRg0KBB\nWLZsGdzd3Y3d5j5Qnr7xsYSPLZ2nXyfu5ETvWttNp3pKefqEJdFr0HdxcUFSUhJCQkLQ0tKCVatW\nITQ0FPn5+Rg+fDhSU1Nx7NgxHD16FHPmzDF2mwkz4Ig+VlbS7C14RxD2gl7hHV9fX4SEhAAAhEIh\nhg0bhsbGRtLTN0E9lsISPrammL4pBe4opk9YEoPv7Lq6Oty6dQthYWEOq7WuK8rhA8A6NPN1gXxM\nEPaDQQ9yOzo6kJGRgaSkpD5xXV201o2FreTp26KIlzl9bE15+qbcv4Dy9AlLovegL5VKkZ6ejvHj\nx3PyugqtdQA6aa0r38yFhYV2dXzp0qU+r/eeFWvagEPbsTkwt49VXTNj2ms71udXiz6bqKgr39zc\nbHB/CEIbeoV3GGPYs2cPgoKCMGPGDO58dHQ0L611Y22AoSrmaelNWRYtWgTg55BO/6FR/w3n3OFs\nNG3AYelNVCzhY8U1M4W9LvdI9xeY6k1t1KHPJirqyvj4+CAqwrDPBG2iQmhD752zCgoKEBwczMns\nJiYmIiEhAZmZmUhJSeHS+RwdW91blXxMEPaJXoN+eHg4Pv30U5WvkZ6+ceuxFJbwsaXz9JWxtpi+\nLd9LhHVh3pgBQVgZyhlVtpJNRRCGYBeDvjXl6SvH8C/Wimkg4YEl8vTVhd8oT5+wV+xi0LcmlAcR\nwLbi+I4Cze4JR8YuBNdMnadfJ+7ExVox93ejsb3HcZ3YdDotjoS58vSV10uoWythbTF9gjAWNNPX\nAQwuCZcAAAtxSURBVFWz942nb3DHadNDEdi/nyWaRtgRbs4CXKwVAwAG9Xeje4owCXYx6Fs6pq/8\nYaVwgf44uvZOY7uUm0woTyQopk8YE7sY9JXprW+jPGPS9JohKH9YKYZvfdiq5hFBmAKjD/rl5eXY\nv38/ZDIZJk+ejGnTphn7LfqgnMfcOxSjPGPq/drOF0TcYCB0d0FLhxTNzc3w8fHhjgEaJHpjKh+b\nKk9fn4fr1hbTp9k+YSyMOujL5XJkZ2fjzTffhJ+fH/785z8jMjISQUFBxnwbo9F7hv5znP5Oj2Oa\nvf+MrfnYVlEOGcr9Q1An7qQYP2EUjDroV1VVITAwEAEBAQCA2NhYnD9/3ugDgvLP9UH93WgWZEZM\n6WNjxvQNTcu0dJ6+8oQEANKm+9CgTxgFo97ZjY2N8Pf35479/PxQVVWlocTPVDQ8wKXb3TObMY/4\nwMPVqU/oRUGXTI61X10HQJkz5sYQHyuQs+40SQFMJ8Ftq5pH6qDMHsJYWM2DXCcAzk7dywacBH0/\ntMqznt4fYk0xT8qssS7qxJ2ob+3+Mn9Y2A8Dvdy419T5Ud0D+MLCQoRGjenxq89Yg6E1xfQB9Zk9\nBMEXAWPMaDt5XL16FZ999hnWrl0LAMjJyYFAIMDMmTP72Obm5hrrbQkeTJ482aDy5GPrx1AfE/aN\nUWf6w4YNQ11dHRoaGuDn54dvvvkGr732mkpbujFtE/IxQdg2Rp3pA93pfB9//DGXzjd9+nRjVk9Y\nAeRjgrBdjD7oEwRBENaLXQiuEQRBELpBgz5BEIQDYZGUTX2X8d+9exdZWVlobm6GUChEXFwc4uLi\n0N7ejszMTDQ0NHD7trq7u+tUp1wux5o1a+Dn54c1a9boXVdHRwf27duHH3/8ERKJBIsXL0ZQUBDv\nuk6fPo38/HxIJBJERERg/vz5Ordp9+7dKCsrg1AoRHp6OgBoLHvy5Enk5eXB2dkZycnJCA8P1+ma\naev7okWLEBYWxr1eUFCA48ePAwCCgoIwffp0fPXVV2rtFVRVVWHt2rWIiIhAW1ubWtuqqirs378f\nHR0d8PDwQEBAgNq6u7q6sHfvXvz4449wcnJCa2srPD09AQD19fX4/e9/3+cZxV//+lf85z//QWNj\nIwYOHAhXV1eVtsr99PPzQ319Pdzc3DTWDQDnzp3Djh078NBDD8HLy0utraKfYrEYjY2NCAwMVFu3\ncj89PDwQHx+vdRN7wkFgZkYmk7GlS5ey+vp6JpFIWEpKCqupqdGpbFNTE7tx4wZjjLHm5ma2YMEC\nVlNTww4ePMiOHTvGGGMsJyeHHTp0SOf2nDhxgu3cuZNt2bKFMcb0riszM5Pl5uYyxhiTSqXswYMH\nvOsSi8Vs8eLFrL29nclkMvbee++xsrIynespLy9n169fZytXruTOqStbU1PDUlJSmEQiYfX19Wzp\n0qVMJpPp1Fdd+q7MlStXuHNnzpxhL730kkZ7xrrvkw0bNrBXXnmFffTRR2ptW1tb2YoVK9jdu3cZ\nY4y9//77Guv+6quvuPoaGhrY0qVLmVwuZzKZjL388svszp07PexLSkrYe++9xxhj7OrVq+yNN95Q\na9u7n2+88QbXF1X2yv3cvHkzO3funFrb3v1sbm7WWLe6fhKE2cM7ysv4XVxcuGX8uuDr64uQkBAA\ngFAoxLBhw9DY2Ijz589jwoQJAIC4uDgUFxfrVN+9e/dQVlaGSZMmgf33ebY+dbW1taGyshKTJk0C\nADg7O8PT05N3XYpZYVtbG7q6utDZ2QkvLy+d64mIiICXl1ePc+rKFhcXIzY2Fi4uLggICEBgYCDv\nlbWa+q5MWFgYdy48PBwPHjzQaA8AX375JUaNGoW2tjY89thjam0LCwvx9NNPw9/fH21tbaiurtZY\nt6enJ9rb2yGVStHa2go3NzcIBAJcunQJgwYNwsCBA9VeP5FIhAcPHuDcuXMqbZX7OWrUKNy7dw8A\n1Nat6OfYsWMhFAo12ir3E4BWe3X9JAizh3eMsYwfAOrq6nDr1i2EhYWhubkZvr6+AAAfHx80Nzfr\nVMf+/fsxd+5ctLe3c+f0qauhoQFCoRBZWVm4fv06RCIRkpOTedfl5uaGBQsWYMmSJXB1dcW0adMg\nEon07p+m/jQ1NUEkEnF2/v7+aGxs1LleBar6/tJLL3FfYL35xz/+AW9vb432ii/yefPm4dixY/jq\nq6/w2WefqbS9ffs2ZDIZ1q9fj+bmZggEAo11jxs3DiUlJfjjH/8IuVyOd999FwBQVFSkcjVw7/vV\n399fJ9XL06dPIzo6Wmvd58+fx/r165GdnQ2BQKDWVrmfHR0deOGFFzB+/Hi19ur6SRA2+SC3o6MD\nGRkZSEpK6hPb1nU2U1JSAqFQiCFDhnCz/N7oWpdMJkN1dTWefvppbN68GVKpFOfOneNdV0tLC/bt\n24cdO3YgKysLV69eRUlJiV5tUoW2svrUrUvfFXz//fcoLS2FWCzWaP/xxx8jMTERcrkcYrEYoaGh\nam1lMhnKy8uxcuVKzJ07F/X19Rg9erRa+1OnTsHZ2Rl79+7F+vXrsWXLFnR1daGkpAQxMTEq2618\nf8jlclRWVqq1VfSzoKAAs2fPhlQqVVu3op8CgQCMMY22yv1ctWoVjhw5gra2NrX2qvopl5MMCWGB\nmb6fnx/3sxfoDrH4+fnpXF4qlSI9PR3jx4/nHkz5+Pjg/v378PX1RVNTE3x8fLTWc+XKFZSUlKCs\nrAwSiYR74KlPXf7+/vD29uZmdrGxsTh79ix8fX151VVVVQWRSMQ9oIuJiUFFRYVebVKgrqyhflCg\nru+KkIiCH374AXv37sVrr72GrVu3arS/fv06MjIyuEEqPz8fw4cPV2nr7++PJ554Ar6+vhCJRHBx\ncYG7uzvc3NxU2peXl2PSpEno168fRCIRBgwYgLy8PAwdOpQLmSjT+zrV1tYiJCREpa1yP9944w14\neXmhuLhYbd2KfgKAWCxGcXExBg0apNJWuZ8AMHToUHzxxRdq61bVz9u3b+Phhx9W2W7CcTD7TF95\nGb9UKsU333zDDQDaYIxhz549CAoKwowZM7jz0dHRyM/PBwCcPXtWpyyFxMREZGdnIysrC8uXL8fI\nkSOxbNkyvery9fVFYGAgrl27BrlcjtLSUkRGRmL06NG86goPD0d1dTVaW1shkUhQVlaGqKgovdqk\nQF3Z6OhoFBUVQSqVoqGhAXV1dQgNDdW5Xk19f/zxx3vY3L17F+np6Vi2bBnCwsK02u/atQtZWVnI\nzs6GUCjEjBkzMGrUKJW2Y8aMQXl5OTo7O+Hi4gKBQABXV1e1dUdGRqKkpARyuRz19fVobW1FZWUl\nYmNj1V6/r7/+GkC37pBEIunzhaaqn4ov7qKiIrV1K/qZlZWFsWPH4pFHHlGbyabcz9bWVty8eRM1\nNTVq61bVTxrwCcBCK3L1XcZfWVmJt956C8HBwVwoIjExEcOHD9c7ZVPRnhMnTmD16tV6p2zW1tYi\nKysLLS0tCA4OxrJly8AY411Xfn4+zpw5g66uLkRFRWHWrFno7OzUqZ6MjAxUVFRALBbDx8cHs2bN\nwtixYzWmbObm5nIpmxERETpfM219LygoAAA899xz2LNnD7799lvuYaNMJoO7u7tae2W2b9+Ompoa\nyOVytbb//Oc/8eWXX0IikeDZZ5/FxYsX1dbd1taGTz/9FJWVlRAKhZg8eTI++ugj7Nq1Cx4eHgCA\nf/3rXz3qP3z4MEpLS+Hm5oa6ujrs3r1bpW3vfgoEAty9e1dj3QoyMzNRXFyMDz/8UK2tcj9/+ctf\nIicnR23dvfs5bdo0jBo1Sg/vEvYGyTAQBEE4EDb5IJcgCILQDxr0CYIgHAga9AmCIBwIGvQJgiAc\nCBr0CYIgHAga9AmCIBwIGvQJgiAcCBr0CYIgHIj/B7b49EmZXT4zAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x113acc8d0>"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Survival Model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Surivival time, in days"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "obs_t = flu_complete.time_to_event.values/24."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(obs_t, bins=25);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAECCAYAAAASDQdFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFXRJREFUeJzt3W9MVHe+x/HPMKOiLDPj0ADJEmLjiJDNZuNm6F1DjIN2\nNUXubZOmTdRNKUmTm6LsupULZvugbomRhCWBHVAfbt19UvcaTdxwbbKOmEqfDC3dB8V1w1a3kmac\nODMMXAPqwLkPTCe1yt+ZYej9vV/POHOO8+UXeDOcORxtlmVZAgAYIS/XAwAAVg7RBwCDEH0AMAjR\nBwCDEH0AMAjRBwCDOOZ78OHDhzp+/LgePXqktWvXavv27aqvr9fU1JQCgYAikYhKSkrU3Nys/Px8\nSVJ/f7+CwaDsdrsaGxtVWVm5Ip8IAGBhtoWu03/w4IHWrVunR48e6dixY2ppadGVK1dUWFiol19+\nWRcvXtT9+/d18OBBjY2NqaenRydPnlQsFlN7e7t6enqUl8cvFACwGixY43Xr1kmSpqenNTs7qzVr\n1mhoaEg7d+6UJPn9foVCIUlSKBRSTU2NHA6HiouLVVpaqtHR0SyODwBYinlP70jS7Oys2tradOfO\nHb355pt67rnnlEgk5Ha7JUkul0uJREKSFI/HtWXLltSxRUVFisViWRodALBUC0Y/Ly9PnZ2dikQi\nOnnypLZu3frE4zabbd7jF3ocALByFoz+N4qLi7Vt2zaNjIzI5XJpfHxcbrdb8XhcLpdLkuTxeBSN\nRlPHRKNReTyeOf/NK1eupDE6AJhr9+7dyzpu3uhPTEzIbreroKBAk5OT+vzzz9XY2Cifz6eBgQG9\n8sorunbtmqqrqyVJPp9PPT09qq+vVywWUzgcltfrnXeAn/70p8saHABM9dlnny372HmjPz4+rr6+\nPs3Ozsrtdqu+vl4//vGP5fV6FQgE1NLSkrpkU5LKyspUW1urtrY22e12NTU1cXoHAFaRBS/ZzKYr\nV67wSh8Aluizzz5b9ukdLqAHAIMQfQAwCNEHAIMQfQAwCNEHAIMQfQAwCNEHAIMQfQAwCNEHAIMQ\nfQAwCNEHAIMQfQAwCNEHAIMQfQAwCNEHAIMs+r9LzJavxqcXtZ9NUmnhWq2x83MKAJYr59F/679v\nLGq/cne+uv9jC9EHgDRQUAAwCNEHAIMQfQAwCNEHAIMQfQAwCNEHAIMQfQAwCNEHAIMQfQAwCNEH\nAIMQfQAwyLz33rl37576+vqUSCTkdDrl9/vl9/t17tw5BYNBOZ1OSdL+/fu1bds2SVJ/f7+CwaDs\ndrsaGxtVWVmZ/c8CALAo80bf4XCooaFBmzZt0sTEhI4ePSqv1yubzab6+nrV19c/sf/Y2JiuXr2q\njo4OxWIxtbe3q6enR3l5/EIBAKvBvDV2u93atGmTJMnpdGrz5s2KxWKSJMuynto/FAqppqZGDodD\nxcXFKi0t1ejoaOanBgAsy6JfgofDYY2NjamiokKSdPnyZf3617/W6dOndf/+fUlSPB5XUVFR6pii\noqLUDwkAQO4tKvrT09Pq7u5WQ0OD8vPztWfPHvX29urEiRPKy8vT2bNn5zzWZrNlbFgAQHoWjH4y\nmVRXV5d27Nih6upqSZLL5ZLNZtOGDRu0d+/e1Ckcj8ejaDSaOjYajcrj8WRpdADAUs0bfcuydObM\nGZWVlWnfvn2p7fF4XJI0MzOj69evq7y8XJLk8/k0ODioZDKpSCSicDgsr9ebxfEBAEsx79U7N2/e\n1Mcff6zy8nK1trZKenx55uDgoG7fvi2Hw6Gqqio1NDRIksrKylRbW6u2tjbZ7XY1NTVxegcAVpF5\no19ZWakPP/zwqe3fXJP/LHV1daqrq0t/MgBAxnEBPQAYhOgDgEGIPgAYhOgDgEHmfSN3tRmfSuqf\n96YWvX9J4VqVFq7L4kQA8P3yvYp+5H8f6tj//HPR+3fWeYk+AHwLp3cAwCBEHwAMQvQBwCBEHwAM\nQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQB\nwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCCO+R68d++e+vr6lEgk5HQ65ff75ff7NTU1pUAg\noEgkopKSEjU3Nys/P1+S1N/fr2AwKLvdrsbGRlVWVq7IJwIAWNi80Xc4HGpoaNCmTZs0MTGho0eP\nyuv1amBgQFu3blVra6suXryo8+fP6+DBgxobG9PVq1fV0dGhWCym9vZ29fT0KC+PXygAYDWYt8Zu\nt1ubNm2SJDmdTm3evFmxWExDQ0PauXOnJMnv9ysUCkmSQqGQampq5HA4VFxcrNLSUo2Ojmb3MwAA\nLNqiX4KHw2GNjY2poqJCiURCbrdbkuRyuZRIJCRJ8XhcRUVFqWOKiooUi8UyPDIAYLkWFf3p6Wl1\nd3eroaEhde7+Gzabbd5jF3ocALByFox+MplUV1eXduzYoerqakmPX92Pj49Levzq3uVySZI8Ho+i\n0Wjq2Gg0Ko/Hk425AQDLMG/0LcvSmTNnVFZWpn379qW2+3w+DQwMSJKuXbuW+mHg8/k0ODioZDKp\nSCSicDgsr9ebvekBAEsy79U7N2/e1Mcff6zy8nK1trZKkg4cOKBXX31VgUBALS0tqUs2JamsrEy1\ntbVqa2uT3W5XU1MTp3cAYBWZN/qVlZX68MMPn/nYNz8Evquurk51dXXpTwYAyDguoAcAgxB9ADAI\n0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcA\ngxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9\nADCIY6EdTp06peHhYTmdTnV1dUmSzp07p2AwKKfTKUnav3+/tm3bJknq7+9XMBiU3W5XY2OjKisr\nszg+AGApFox+bW2tXnrpJfX29qa22Ww21dfXq76+/ol9x8bGdPXqVXV0dCgWi6m9vV09PT3Ky+MX\nCgBYDRascVVVlQoKCp7ablnWU9tCoZBqamrkcDhUXFys0tJSjY6OZmZSAEDaFnylP5fLly8rGAyq\noqJCb7zxhgoKChSPx7Vly5bUPkVFRYrFYhkZFACQvmWdd9mzZ496e3t14sQJ5eXl6ezZs3Pua7PZ\nlj0cACCzlhV9l8slm82mDRs2aO/evalTOB6PR9FoNLVfNBqVx+PJzKQAgLQtK/rxeFySNDMzo+vX\nr6u8vFyS5PP5NDg4qGQyqUgkonA4LK/Xm7lpAQBpWfCcfnd3t27cuKGJiQm9/fbbeu211zQyMqLb\nt2/L4XCoqqpKDQ0NkqSysjLV1taqra1NdrtdTU1NnN4BgFVkwegfOXLkqW27du2ac/+6ujrV1dWl\nNxUAICu4gB4ADEL0AcAgRB8ADEL0AcAgRB8ADLLs2zB8H6y12/S3ryeXdExJ4VqVFq7L0kQAkFv/\nr6Mfm0rqt3+9taRjOuu8RB/A/1uc3gEAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI\n0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcAgxB9ADAI0QcA\ngzgW2uHUqVMaHh6W0+lUV1eXJGlqakqBQECRSEQlJSVqbm5Wfn6+JKm/v1/BYFB2u12NjY2qrKzM\n7mcAAFi0BV/p19bW6je/+c0T286fP6+tW7fqd7/7nbZs2aLz589LksbGxnT16lV1dHTo6NGj6uvr\n0+zsbHYmBwAs2YLRr6qqUkFBwRPbhoaGtHPnTkmS3+9XKBSSJIVCIdXU1MjhcKi4uFilpaUaHR3N\nwtgAgOVY1jn9RCIht9stSXK5XEokEpKkeDyuoqKi1H5FRUWKxWIZGBMAkAlpv5Frs9nSehwAsHKW\nFX2Xy6Xx8XFJj1/du1wuSZLH41E0Gk3tF41G5fF4MjAmACATlhV9n8+ngYEBSdK1a9dUXV2d2j44\nOKhkMqlIJKJwOCyv15uxYQEA6Vnwks3u7m7duHFDk5OTevvtt/X666/r1VdfVSAQUEtLS+qSTUkq\nKytTbW2t2traZLfb1dTUxOkdAFhFFoz+kSNHnrm9tbX1mdvr6upUV1eX3lQAgKzgL3IBwCBEHwAM\nQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQB\nwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAM4sj1AKvN\nWrtNf/t6ctH7lxSuVWnhuixOBACZQ/S/IzaV1G//emvR+3fWeYk+gO+NtKJ/6NAhrV+/Xnl5ebLb\n7Tp58qSmpqYUCAQUiURUUlKi5uZm5efnZ2peAEAa0n6lf/z4cf3gBz9IfXz+/Hlt3bpVra2tunjx\nos6fP6+DBw+m+zQAgAxI+41cy7Ke+HhoaEg7d+6UJPn9foVCoXSfAgCQIWm90rfZbHr//fdls9m0\nZ88evfjii0okEnK73ZIkl8ulRCKRkUEBAOlLK/rt7e3auHGjxsbGdPLkSf3whz984nGbzZbWcACA\nzErr9M7GjRslSWVlZXrhhRc0Ojoql8ul8fFxSVI8HpfL5Up/SgBARiw7+g8ePNDU1JQkaWJiQsPD\nwyovL5fP59PAwIAk6dq1a6qurs7IoACA9C379E4ikVBnZ6ckqbCwUPv27dNPfvITVVRUKBAIqKWl\nJXXJJgBgdVh29IuLi1PR/7b169ertbU1raEAANnBvXcAwCBEHwAMQvQBwCBEHwAMQvQBwCBEHwAM\nQvQBwCBEHwAMQvQBwCBEHwAMQvQBwCD8x+hpWmu36W9fTy56/5LCtfxH6gByhuinKTaV1G//emvR\n+3fWeYk+gJzh9A4AGIToA4BBiD4AGIToA4BBiD4AGIToA4BBiD4AGIToA4BBiD4AGIToA4BBiD4A\nGIR776xy4ckHujv5cEnHOPMdmphOLnp/bgIHmIPor3J3Jx/qv/pHl3TMey8+z03gADwT0Qe3hwYM\nQvRX2FID+3BmNovTPMbtoQFzZCX6IyMj+uCDDzQzM6Pdu3frpZdeysbTfC8tNbDvvfh8FqcBYJqM\nX70zOzur06dP6+jRo+ro6FAwGNTY2FimnwYAsAwZj/7o6KhKS0tVXFwsh8OhmpoaDQ0NZfppAADL\nkPHTO7FYTEVFRamPPR6PRkeXdvUJVreVeON3qZeqZvsy1eVcOssb3liNcv5G7n/+2w8XtZ8z3y5b\nlmfB4iz1fYmef9+y5GA+nJnVux99uej9l3qZ6lJnWuo80up7w3upP7hW4ocWM608m2VZVib/wX/8\n4x/685//rHfffVeSdOHCBdlsNr3yyitP7XvlypVMPjUAGGP37t3LOi7jr/Q3b96scDisSCQij8ej\nTz75RL/61a+eue9yhwYALE/GX+lLjy/Z/MMf/pC6ZLOuri7TTwEAWIasRB8AsDpxl00AMAjRBwCD\n5OSSzdVym4ZDhw5p/fr1ysvLk91u18mTJzU1NaVAIKBIJKKSkhI1NzcrPz8/q3OcOnVKw8PDcjqd\n6urqkqR55+jv71cwGJTdbldjY6MqKytXZKZz584pGAzK6XRKkvbv369t27at2Ez37t1TX1+fEomE\nnE6n/H6//H5/TtdqrplyuVYPHz7U8ePH9ejRI61du1bbt29XfX19Ttdprply/TUlPb6LwLFjx+Tx\neHTs2LGcf+/NNVfG1spaYTMzM9bhw4etu3fvWo8ePbJaWlqsO3furPQYlmVZVlNTkzU5OfnEtj/+\n8Y/WxYsXLcuyrAsXLlh/+tOfsj7HyMiI9eWXX1rvvPPOgnPcuXPHamlpsR49emTdvXvXOnz4sDUz\nM7MiM507d866dOnSU/uu1EzxeNy6deuWZVmWlUgkrLfeesu6c+dOTtdqrplyvVbT09OWZVnWw4cP\nrXfeecf6+uuvc/419ayZcr1OlmVZly5dsnp6eqyOjg7LsnL/vTfXXJlaqxU/vbPabtNgfed97KGh\nIe3cuVOS5Pf7FQqFsj5DVVWVCgoKFjVHKBRSTU2NHA6HiouLVVpampW/eH7WTNLT67WSM7ndbm3a\ntEmS5HQ6tXnzZsVisZyu1VwzSbldq3XrHv+x0PT0tGZnZ7VmzZqcf019e6aZmRmtWbNGUm7XKRqN\nanh4WLt27UrNket1mmsuy7IyslYrfnpnNd2mwWaz6f3335fNZtOePXv04osvKpFIyO12S5JcLpcS\niUROZptrjng8ri1btqT2KyoqSkVmJVy+fFnBYFAVFRV64403VFBQkJOZwuGwxsbGVFFRsWrW6tsz\n3bx5M6drNTs7q7a2Nt25c0dvvvmmnnvuuZyv07NmknL7NfXBBx/oF7/4haamplLbcr1Oc81ls9ky\nslZGv5Hb3t6uzs5O/fKXv9SFCxd048aNJx632VbHjR8WmmOl5tyzZ496e3t14sQJ5eXl6ezZszmZ\naXp6Wt3d3WpoaHjq/ZZcrdV3Z8r1WuXl5amzs1O///3v9dFHH+nWrSdvUZGLdXrWTLlcp08//VRO\np1PPP//8M19BL+Y5s7FOc82VqbVa8eh7PB5Fo9HUx9FoVB6PZ6XHkCRt3LhRklRWVqYXXnhBo6Oj\ncrlcGh8fl/T4J7vL5crJbHPNkcv1c7lcstls2rBhg/bu3Zv6DW0lZ0omk+rq6tKOHTtUXV2dmiuX\nazXXTLleK0kqLi7Wtm3bNDIykvN1mmumXK3TzZs39emnn+rQoUPq6enRF198oUAgkPN1etZcvb29\nGVurFY/+t2/TkEwm9cknn8jn8630GHrw4EHqV6eJiQkNDw+rvLxcPp9PAwMDkqRr166lvolX2lxz\n+Hw+DQ4OKplMKhKJKBwOy+v1rshM8XhckjQzM6Pr16+rvLx8RWeyLEtnzpxRWVmZ9u3bl9qey7Wa\na6ZcrtXExITu378vSZqcnNTnn38+79d2Lmf6Jq65WKcDBw7o9OnT6uvr05EjR/SjH/1Izc3NOf/e\ne9Zchw8fztjXVE7+Inc13KYhEomos7NTklRYWKjt27fr5z//eU4u2ezu7taNGzc0OTkpl8ul119/\nXT/72c/mvWzsypUrqcuzqqqqsjbTxMSE3G63XnvtNY2MjOj27dtyOByqqqrSyy+/nDr3uRIz/f3v\nf9d7772n8vLy1K+vBw4c0NatW3O2Vs+aaf/+/RocHMzZWn311Vfq6+vT7Oys3G63tm/frl27di14\nKWIuZurt7c3p19Q3RkZGdOnSJbW1teV0nb7riy++0F/+8he1tbUpEAjoX//6V9prxW0YAMAgRr+R\nCwCmIfoAYBCiDwAGIfoAYBCiDwAGIfoAYBCiDwAGIfoAYJD/A1bpFSAoDAOiAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x113acb990>"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Fate of each patient"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "died = (flu_complete.fate=='Died').astype(int).values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Unique failure times"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "times = np.unique((flu_complete.time_to_event[flu_complete.fate=='Died'].values/24).astype(int))\n",
      "times"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,\n",
        "        13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,\n",
        "        26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,\n",
        "        39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  50,  53,  54,\n",
        "        55,  59,  61,  62,  63,  67,  68,  73,  77, 108, 121, 124, 125,\n",
        "       126, 133, 142, 150, 201, 227, 235, 301])"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N_obs = len(obs_t)\n",
      "T = len(times) - 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate risk set"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = np.array([[int(obs >= t) for t in times] for obs in obs_t])\n",
      "Y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 31,
       "text": [
        "array([[1, 1, 1, ..., 0, 0, 0],\n",
        "       [1, 1, 1, ..., 0, 0, 0],\n",
        "       [1, 1, 1, ..., 0, 0, 0],\n",
        "       ..., \n",
        "       [1, 0, 0, ..., 0, 0, 0],\n",
        "       [1, 1, 1, ..., 0, 0, 0],\n",
        "       [1, 1, 1, ..., 0, 0, 0]])"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Counting process. Jump = 1 if $\\text{obs}_t \\in [ t_j, t_{j+1} )$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dN = np.array([[Y[i,j]*(times[j+1] >= obs_t[i])*died[i] for i in range(N_obs)] for j in range(T)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Process for one patient\n",
      "dN[:,1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 33,
       "text": [
        "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,\n",
        "       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
        "       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
        "       0, 0, 0])"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Covariates"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "standardize = lambda x: (x - x[np.isnan(x)-True].mean()) / x[np.isnan(x)-True].std()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "male, AgeYears, non_white, HoursECMO, pH, PCO2, PO2, HCO3, SaO2, mrsa, any_coinf = \\\n",
      "    flu_complete[covariates].values.T\n",
      "SaO2 /= 100.\n",
      "SaO2 -= 1e-6\n",
      "\n",
      "HoursECMO_std = standardize(HoursECMO)\n",
      "AgeYears_std = standardize(AgeYears)\n",
      "pH_std = standardize(pH)\n",
      "PCO2_std = standardize(PCO2)\n",
      "PO2_std = standardize(PO2)\n",
      "HCO3_std = standardize(HCO3)\n",
      "SaO2_std = standardize(SaO2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Model"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import theano.tensor as tt\n",
      "from pymc import Model, Normal, Poisson, DensityDist, T as StudentT, Uniform, Deterministic, HalfCauchy, Laplace, Exponential\n",
      "from pymc.distributions.timeseries import GaussianRandomWalk\n",
      "from pymc import sample, psample, NUTS, Slice, Metropolis, forestplot, find_MAP"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with Model() as model:\n",
      "    \n",
      "    X = np.c_[[male, pH_std, PCO2_std, PO2_std, AgeYears_std, HoursECMO_std, non_white, mrsa, any_coinf]].T\n",
      "    \n",
      "    # Intercept for survival rate\n",
      "    beta0 = Normal('beta0', 0.0, 0.001)\n",
      "    # Treatment effect\n",
      "    beta = Normal('beta', 0, 0.001, shape=X.shape[1])\n",
      "\n",
      "    # Survival rates\n",
      "    lam = tt.exp(beta0 + tt.dot(X, beta))\n",
      "\n",
      "    def exp_logp(failure, value):\n",
      "        # Exponential survival log-likelihood\n",
      "        return tt.sum(failure * tt.log(lam) - lam * value)\n",
      "\n",
      "    surv_like = DensityDist('surv_like', exp_logp, observed=(died, obs_t))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 79
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with model:\n",
      "    trace = psample(1000, step=NUTS())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\r",
        " [--                6%                  ] 66 of 1000 complete in 0.5 sec\r",
        " [-----            14%                  ] 140 of 1000 complete in 1.0 sec\r",
        " [-------          20%                  ] 207 of 1000 complete in 1.5 sec\r",
        " [----------       28%                  ] 286 of 1000 complete in 2.0 sec\r",
        " [--------------   37%                  ] 377 of 1000 complete in 2.5 sec\r",
        " [---------------- 44%                  ] 445 of 1000 complete in 3.0 sec\r",
        " [-----------------52%                  ] 525 of 1000 complete in 3.5 sec\r",
        " [-----------------61%---               ] 612 of 1000 complete in 4.0 sec\r",
        " [-----------------69%------            ] 697 of 1000 complete in 4.5 sec\r",
        " [-----------------78%---------         ] 783 of 1000 complete in 5.0 sec\r",
        " [-----------------87%-------------     ] 870 of 1000 complete in 5.5 sec\r",
        " [-----------------96%----------------  ] 969 of 1000 complete in 6.0 sec\r",
        " [-----------------100%-----------------] 1000 of 1000 complete in 6.2 sec"
       ]
      }
     ],
     "prompt_number": 80
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "forestplot(trace, vars=['beta'], ylabels=['male', 'pH', 'PCO2', 'PO2', 'Age', 'HoursECMO', 'non_white', 'mrsa', 'any_coinf'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 82,
       "text": [
        "<matplotlib.gridspec.GridSpec at 0x125108850>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEjCAYAAACB7F6fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVGX/P/D3GSZEQMABERU3EIQBlVBDU1QgSbCszDJL\nRW17Mu2XS4RWRqaBuORalo+4pLZpj+bToxmguWspmMzkSppoRIpiuLHdvz/I+YowOMjAzC3v13Vx\nXZwzZ859n5nhvM99zmcOihBCgIiISAIqS3eAiIjIVAwtIiKSBkOLiIikwdAiIiJpMLSIiEgaDC0i\nIpIGQ4ukNmLECPTt29fk6crEx8fDx8en1vooI5VKhTVr1li6G2QG27Ztg0qlwrlz5yzdFbNgaEnq\nypUriIuLg7e3N5ycnNCjRw+kpKSUW2b58uVQqVQVftLS0gzL6HQ6dOnSBU2aNMHLL7+M7OzscuuY\nOHEiRo8ebXK/9u3bh8ceeww+Pj5wcHCAn58fXnnlFRw/frxmG1wFRVEMvy9YsABr164t99itj5uy\njrvVpk0bTJ8+vVrPeeihhzBy5Mgat031w4gRIwx/x2q1Gs2bN8dTTz1Vq39fN61atQoqleUjw/I9\noLvy0ksvYc2aNXjnnXewY8cOhIWFISoqCjt37iy3nI2NDf7880/k5OQYfkJDQw2PT58+HYGBgdiy\nZQt+++03jB8/3vDY/v378Z///AczZ840qU8ffvghevbsiVatWiE5ORmZmZmYO3cuioqKEBcXZ/R5\nxcXF1dz68m79fnyjRo3g7Oxc6WOmruNumRqQtaWmryPJoVevXsjJycGpU6ewYMEC/PTTT3jkkUcs\n3a26I0g6165dE/fdd5+YPXt2ufmdO3cWUVFRhully5YJtVpd5bo8PDzEgQMHhBBCfPXVVyI4OFgI\nIcSNGzdEYGCg2LJli0l90uv1wtbWVsyYMaPSx8+fPy+EEGLr1q1CURSxc+dO0a9fP9G4cWOxePFi\nIYQQK1asEPfff79wdHQUPj4+Yvr06aK4uNiwjr///luMGDFCNGnSRLRv317MmjVLxMTEiIceesiw\njLHppUuXCq1WK9zc3MTw4cPF33//bVjm3XffFe3atSvX371794rw8HCh0WhEixYtxMiRI8WFCxeq\nfA3atGkjpk+fbphu3bq1eP/998XkyZNFmzZtRNOmTcW4ceMM2xQTEyMURSn38+OPPwohhLh06ZJ4\n5ZVXhJeXl3BychI9evQQ27dvN6y7stdx3rx5wt7eXqxZs6Zcv86ePStsbGxEamqqEEKI1atXiwce\neEA4OzsLNzc30b9/f3Hs2LFyz1EURaxevdowvWrVKjFgwADRuHFjodFoRK9evUR2dnaVrweZX0xM\njOjbt2+5efHx8UJRFHH58uVKn3Pzs7Jnzx4RFRUlXFxchFarFZs2bSq33OTJk4W/v7+wt7cXLVu2\nFP/6179Efn5+uXXc+jNy5Mja2cg74EhLQkVFRSgpKYGjo2O5+Y0aNcKOHTvKzSspKYG3tzeaN2+O\nsLAwfPfdd+Uef/jhh7Fu3Trk5+djw4YNeOyxxwAAU6dORUhIyB2vB930+eefw9bWttxI7Vaurq7l\npl9++WWMGTMG6enpeOSRRzBnzhy8+eabeP7557F3715MnDgRCxYswHvvvWd4zrBhw7Bz50589dVX\n+Oabb5CRkYGNGzeWG91UNto5cOAANm7ciK+++gpr167F7t27MXz4cKPb8tNPPyE0NBShoaHYsmUL\nFi5cCJ1Oh4EDB5r0Wtxq7ty5yM/Px9dff424uDjMmzcPK1asAADMnz8foaGhGDx4sGEU3L17d5SU\nlCAsLAy//fYbPvzwQ2zduhWdO3dG3759ceTIEaOv45NPPonHH38cn332WbllVq1ahRYtWiA8PBwA\nUFhYiClTpuDgwYNYuXIlrl69iv79+6OoqKjSbThw4ACGDRuGPn36YM+ePdi+fTtiYmKq/VqQeYh/\nzgqUlpbi0KFD+Prrr9G1a1c0atSoyue9/PLL6N+/P7Zs2YJ27dph8ODBuHTpkuFxe3t7LFmyBIcP\nH8a0adOwefNmvPbaawCAHj16YOHChQBg+KzOmzevlrbwDiwSlVRjoaGhIiAgQOzYsUPk5+eL5ORk\noVKphEqlEpcuXRJCCLFnzx6xfPlykZ6eLtatWycGDx4sFEURX331lWE9WVlZ4l//+pfQarXijTfe\nEOfPnxfp6emiTZs24sKFC2LcuHHCy8tLhIWFiePHjxvtz6BBg0SnTp3u2O+bR2zvv/++YV5hYaFw\ndHQUS5YsKbfstGnThIuLixBCiOPHjwtFUcTSpUsNj1+6dEmo1epyR56VjbQURRGnTp0yzFuyZIlQ\nFEWcPHlSCFFxpBUZGSmee+65cn3ZuXOnUBRFZGRkGN22ykZagYGB5Zbp1KmTGDJkiGH6oYceqnDE\numbNGuHs7GwYnQohRGlpqWjfvr14/fXXhRCVv45CCLF582ahVqtFTk6OYV5gYKCYPHmy0X7rdDqh\nKIrYtWuXYd6tI620tDRhZ2d3x5Em1b6YmBihVquFo6OjsLOzE4qiiH79+lX53tz8rMyaNcswLz09\nXSiKUuWZlI8++kg0aNDAMP3ZZ58JRVHMsyE1wJGWpFatWoXmzZujd+/e0Gg0+PjjjzF27FgIIaBW\nqwEA3bp1Q0xMDIKCgjBw4EB88cUXiImJKVcs0LZtW3z88cfQ6XRISkqCs7MzRo0ahYULF2Lfvn3Y\nuXMndu3ahdDQULzzzjtV9klU47pQv379DL8fP34cV65cweuvv45GjRoZfqZPn47Lly/jwoUL0Ov1\nAIDIyEjD85ydndGtW7c7tuXt7Y3WrVsbpm+OHm+u83bp6elYt25dub7069cPiqLgxIkTJm+joigV\nrjUEBwfjzz//rPJ56enpuHLlCtq0aWNo38nJCVlZWRXav/V1BMoKO9zd3Q2VfwcPHoROpys3sszI\nyMATTzwBLy8vODk5ISQkBADw+++/V9qfPn36IDAwEL6+vhgyZAiWLFmCCxcumPYikNl169YNhw4d\nwvbt2zF+/Hhs377d8Fl2dHQ0fGb69+9f7nm3fhY7duwIlUpV7rP4zTffoFevXmjRogUaNWqEiRMn\noqioCDk5OXWzYSZSW7oDdHdatWqFLVu24Pr167h+/TpcXFwwbtw4NGvWDA4ODkafFxISUq667naJ\niYkICAhA//79MXr0aAwcOBAeHh546aWXEBgYiJKSEtjY2FR4nlarxaZNm1BUVIT77rvvjv1v3ry5\n4feSkhIAwNq1a+Hr61th2caNGxtdjylBWZ0wBcpOu8TFxWHYsGEVHmvatGm11nX7KRuVSoXS0tIq\nn1NSUgJ/f3+sX7++wmP29vblpm99HYGywpvnnnsOK1euxLhx47By5Uo88MADaN++PQDg6tWriIyM\nRMeOHbF8+XI0bdoUV69eRZcuXVBYWFhpfxRFwb59+7Br1y6kpKRg8eLFiI2NRUZGRrmDAaobdnZ2\n8PLygpeXF7p27YqcnByMHj0aGRkZ+OWXXwzLNWzYsNzzbv0sqlQqKIpi+Czu27cPTz31FF599VXM\nnj0bjRs3RmpqKl555RWjnwtL4UhLcnZ2dnBxccGNGzfwxRdfYNCgQVUuf/DgQXh7e1f6mF6vx5Il\nSzB//nwAZdfOrl+/DgC4du0a8vPzcfXq1Uqf+8wzz6C4uBhz5syp9PGqjszbt28PBwcHnDx50vDH\neOuPSqWCVqsFAGzZssXwvPz8fOzbt6/K7QWArKwsnD592jD9ww8/AIBhnbcLDg5GZmZmpX2p6oDg\nbtja2lao+uvcuTOysrLQqFGjCu17eHjccZ3Dhw/HoUOHkJGRgc8//7zcKOvXX3/F+fPnMXXqVPTq\n1Qvt27dHZmbmHYNdpVIhNDQU7733Hn766Sc0atTIcG2O6tbt12zj4+Oh1+vx5ZdflvusNGvWzOR1\n7ty5Ew0bNsScOXPQtWtXtGvXDvv37y+3jK2tLQDzVNrWBEdakkpJScGNGzfg7++PM2fOYMqUKWjQ\noAHi4+MNy8THxyMkJAQ+Pj44deoUvv76ayQnJ1e6sykpKcGoUaPw4YcfGkY2jz76KGJjYxEWFoYV\nK1YgKCjI6MVef39/TJ8+HXFxcThz5gwGDx6MFi1a4MSJE1i7di3y8vKMjvBsbW3xwQcfYPLkyVAU\nBRERESguLsbhw4eRkZGBxMREtGvXDgMGDEBiYiK8vLzg7u6OxMREODk53fGP6OYo9P3338eFCxeQ\nlJSExx9/HF5eXpUun5iYiO7du2PChAkYNmwYGjVqhOPHj2Pt2rVYuHAh7OzsKn3e7f0w5Y+7bdu2\n2Lp1KzIyMuDp6QkXFxc888wzmDdvHvr374/p06fDx8cHf/75J9LS0qDVag3FMsYEBgbi/vvvx8iR\nI3H58mUMGTLE8Fjr1q3RoEEDLFy4EI6Ojvjrr7/w0UcfVVmqv2HDBvz2228IDQ2Fk5MTduzYgXPn\nzqFly5Z33D4yv9s/Vzf/NpKSksq919Xh5+eHq1evYtasWXjiiSewb98+pKamllumbdu2AIDVq1ej\nb9++cHR0NPtBnEksdTGNambdunXCx8dH2NnZiXbt2okXXnhBXLx4sdwy48ePF23bthUNGzYUGo1G\n9OjRQ3zzzTeVrm/mzJniqaeeKjfvxo0bIiYmRjg7O4vg4GBx+PDhO/Zrz549YsCAAcLLy0vY29uL\n9u3bi1dffdVQ9LB161ahUqnE2bNnKzz366+/Ft26dRPOzs6icePGolu3boZyeCHKSt5jYmKEm5ub\naNeunUhKShIjRowoV4hhbPrf//638PPzE66urhVK3uPj44WPj0+5vmRkZIj+/fsLDw8P4eDgIPz9\n/cuVq1fm9kKM26eFEOKFF14QYWFhhumsrCzRq1cv4ejoKFQqlaHk/e+//xYTJ04U7du3F3Z2dqJF\nixZi4MCBhkKQql5HIYSYN2+eUBRFDBw4sMJj//vf/0RUVJRo3Lix6N27tzh06JBQq9VixYoVhmVu\nLcTYvn27CA8PF02aNBF2dnbC19fX6FcbqHbd/vm+affu3UKlUonvv/++wmPGPiu3v+dz584VnTt3\nFk2bNhWDBw8WGzduFCqVSpw+fdqwzOuvvy7c3d0tWvKuCMH/XExERHLgNS0iIpIGQ4uIiKTB0CIi\nImkwtIiISBoseTfi9nJPIpJPRESExdrmPqTmKnv/WD1oxMWLFy3dBTJCo9EgLy/P0t0gK1fVnVTq\nAvchNWPs/ePpQSIikgZDi4iIpMHQIiIiaTC0SDqxsbGW7gIRWQgLMYzgRVQiubEQQ24sxCAiIukx\ntIiISBr3fGjpdDokJiZauhtkxYQQ0Ol00Ol0Fv8Hd0RUtXs+tIiqIoTAyJEfICrqFKKiTmHUqAQG\nF5EVk+I2Trm5uUhISICfnx90Oh2CgoLQo0cPrFq1CgDw4osvorCwEMuXL0dhYSFatWqFQYMGVfjX\n5NevX0dycjKOHz8OIQSGDx+O4OBgS2wS1UBiYiLi4uKqXEajqc5F+FmG3zZsGApX1+r1Jy+PF9zJ\nNEII6PV6AIBWq63yP0ZT5aQILQA4d+4cXnzxRTz//POYMGECLl26hHfffRc7d+5ESkoKnn32WUyd\nOhUqlQq7d+/G6tWrMWHChHLr+Oabb9CmTRuMHj0aly5dQlJSEkOrllQvNKorCUlJM2px/dVTu9tq\nHMNSLjdH9WlpQQCAiIgNSE6exOCqJmlCS6PRQKvVAgC8vb0RGBgItVoNX19fbN68GYWFhVizZg1+\n/fVXCCFQUlJSYR2//PILioqKsG3bNgDAlStXkJubC3d397rclHqhNneoGs2d1296kAgAbwEI/Gda\nB2AaANN2JAwOMpVer0daWhAKCoYCAFJTy+YFBARYuGdykSa07O3tDb+r1WrDtFqtRlFREb7//ns0\natQICQkJyM7OxsyZMytdz/PPP28IP7IOdzNSMd/oRgEwHcDhf6aHwNTAMm8/qsZwJCojTWjdSV5e\nHvz8/AAAKSkplS7TqVMnpKSkoG3btmjYsCF+++03tG3bti67SZWo7g7ZlJGWqSqesvmMp2yoVmi1\nWoSHr0daWtl0RMQhaLXRlu2UhKQJrdt3IrdOK4qCqKgoLFmyBP/973/Ro0ePCo8DwJNPPonly5dj\n4sSJsLOzg7u7O95888262QCySoqiYNmyybdcHI9mYFGt4GfNPHgbJyN4CxbrZUr1IBFv4yQ3Y+8f\nQ8sIfuCI5MbQkhvvPUhERNJjaBERkTQYWkREJA2GFhERSYOhRdLhXfuJ6i9WDxrByh/rpdFokJeX\nZ+lukJVj9aDcWD1IRETSY2gREZE0GFpERCQNhhYREUmDoUXSiY2NtXQXiMhCWD1oBCt/iOTG6kG5\nsXqQiIikx9AiIiJpMLSIiEgaDC0iIpIGQ4ukw3sPEtVf9Sa04uPjkZWVZZjOzc3FhAkTLNgj+Qgh\noNPpoNPpYMmi06SkJIu1TUSWVW9CS1EUS3dBakIIjBz5AaKiTiEq6hRGjUqwaHARUf2ktnQHzC03\nNxcJCQnw9fVFZmYmvL29MWbMGEt3SxoaTVXfbZll+G3DhqFwdTW+ZF4ev6NCdDshBPR6PQBAq9Xy\nYPou3HOhBQDnzp3DkCFD8MILL2DevHk4ePAgAGD+/PmwtbUFABQXF0OlurcGmlUHTt2qSV8YeHQv\nunm2Ii0tCAAQEbEBycmTGFzVdE+Glr29PR544AEAQM+ePZGRkQEAeO211+Dl5QUA+Ouvv+65C/rm\n2tlXHjgCwFsAAv+Z1gGYBqDiHxxDh6givV6PtLQgFBQMBQCkppbNCwgIsHDP5HJPhpYp6tP1GPOM\nwBQA0wEc/md6CCoLrJq2Z0rg8d6DRPXXPRlaV69exf79+3H//fdj9+7dePDBB/H9999bulsWY46R\nT8VTG59Z7NRGXFxcnbdJVFNarRbh4euRllY2HRFxCFpttGU7JaF7MrSaN2+On3/+GStWrIC3tzeC\ng4MrDS2eSzadoihYtmzyLReRo/n6EVUD/4bM4567y3tubi5mzJiB2bNn12g9vEMzkdx4l3e51au7\nvPPohYjo3nTPjbTMhUdJRHLjSEtu9WqkRfe2e+2rCkRkOo60jOBRkvXSaDTIy8uzdDfIynGkJTeO\ntIiISHoMLSIikgZDi4iIpMHQIiIiaTC0SDq89yBR/cXqQSNY+UMkN1YPyo3Vg0REJD2GFhERSYOh\nRURE0mBoERGRNBhaJB3ee5Co/mL1oBGs/LFevPcgmYLVg3Jj9SAREUmPoUVERNJQW7oDVRk8eDBa\nt24NRVEQGBiIIUOGQK1Wo6CgACtXrsRvv/2GkpISNGnSBCNGjECzZs2Qm5uLpUuXIi8vDyqVCr17\n90Z0dDQAYMeOHfj2228BAJ6ennjiiSfQqlUrS24iERFVg1WHVoMGDZCUlITi4mLMnDkThw4dQnBw\nMBISEhAWFobRo0cDAI4dO4aLFy/C1dUV7733Hl5++WV07NgRhYWFmDlzJmxsbPDwww+jadOmeO+9\n92Bvb49t27bhk08+wfTp0y28lWQKIQT0er2lu0FEFibF6UG1Wo3AwEAcOXIEOp0OarUaDz30kOFx\nX19faLVa7Ny5Ey4uLujYsSMAwNbWFtHR0Vi/fr1hOXt7ewBAcHAwLly4UPcbQ9UmhMDIkR8gMjIL\nkZFZaN++O1g/RDISQkCn00Gn0/EzfJeseqR105UrV3Dw4EE8+eSTOHPmDLy8vCpdLjs7G506dSo3\nT6vV4vLly7h+/Trs7OwM81NSUtClS5da7TdVn0ZjrOJrluG3o0eHw9XV+Dry8li1Rdbn5sFXWloQ\nACAiYgOSkydBURQL90wuVh1ahYWFiI2NRcOGDdG1a1dotVqcOXOmxuvNzMzEjh07MG3aNDP0kipj\nPHysu20GHtUWvV6PtLQgFBQMBQCkppbNCwgIsHDP5GLVoWVra4ukpKRy81q2bIm9e/dWurynpydS\nU1Px9NNPG+bp9Xo4OzsbRlmnT5/Gp59+ismTJ8PBwaH2Ol/P3brzr3mACQBvAQj8Z1oHYBqAmh+h\nMqSI5GLVoVWZwMBAfP7550hJSTFc1zpx4gQKCwvRs2dPrFu3DocPH0aHDh1QWFiIzZs3Y8CAAQCA\n8+fPY/bs2Rg7diw8PDwsuRn1Sk2DQQiBESNuICXlFACgb98bWLbsIk+rkFS0Wi3Cw9cjLa1sOiLi\nELTaaMt2SkJWfUeMmJgYrFixosL8goICrFixAqdOnUJpaamh5N3Dw6NcybuiKOjVqxceeeQRAMDi\nxYuxf/9+uLm5AQBsbGyQkJBQadv8Nrt1ubV6UKvVMrDojqzxjhj8HJvO2Ptn1aFlSQwt65WYmIi4\nuDhLd4OsnDWGFpmOoVVN/MBZL957kEzB0JIb7z1IRETSY2gREZE0GFpERCQNhhYREUmDoUXSiY2N\ntXQXiMhCWD1oBCt/iOTG6kG5sXqQiIikx9AiIiJpMLSIiEgaDC0iIpIGQ4ukk5iYaOkuEJGFsHrQ\nCFb+WC/ee5BMwepBubF6kIiIpMfQIiIiaTC0iIhIGgwtIiKSBkOLpMN7DxLVX9JWDw4ePBitW7eG\noigIDAzEkCFDoFarcfnyZSQnJ+Ps2bMAgM6dO+Ppp5+GSqXCL7/8gjVr1qC4uBju7u4YMGAA/Pz8\nKl0/K3+I5MbqQbndc9WDDRo0QFJSEqZPn47s7GwcOnQIQggkJCSgU6dOmDlzJmbMmIG8vDysWbMG\nAODk5IS4uDjMmjULjzzyCObMmWPhrbh3CSGg0+mg0+kg6XEREVkhaUPrJrVajcDAQBw5cgSZmZnI\ny8tDWFgYAEClUmHAgAHYtGkTCgsL0aZNG7i4uAAA/P39UVRUhOLiYkt2/54khMDIkR8gKuoUIiOz\nMGjQJJSWllq6W0QWx4O5mlNbugM1deXKFRw8eBBPPvkkzpw5gw4dOpR73NPTE05OTsjJyUGrVq0M\n83ft2gVfX1+o1dK/BBal0Rg7BTPL8NvWrcPh5lb5Unl5PIVC9cPNg7m0tCAAQETEBiQnT4KiKBbu\nmVyk3WMXFhYiNjYWDRs2RNeuXaHVanHmzJlKl739Q3HmzBl8+eWXeOedd+qiq1bLeODI1QcGH8lA\nr9cjLS0IBQVDAQCpqWXzAgICLNwzuUgbWra2tkhKSio3z9PTE+vXry83Lzs7G/n5+fDw8AAAXLhw\nAbNmzcLYsWPh7u5eZ/21Rqbs7O8uVASAtwAE/jOtAzANgLmOKOMBxDOsiOohaUOrMh06dIBGo8G2\nbdvQp08flJaW4ttvv0W/fv1ga2uLK1euIDExEc899xx8fX0t3V0p3G0wlJZOwJNPjsf+/e2hUvnh\noYdizXYqRKN5D3l5/6/G6yGqS1qtFuHh65GWVjYdEXEIWm20ZTslIWlL3mNiYrBixYoK828veb//\n/vvxzDPPQKVSYd26dVi/fj2aNWtmWP7tt9+Gk5NThfWwXLXmhBDQ6/UAyv5gzXXunjfMJVNYY8l7\nbf1N3IuMvX/ShlZtY2hZL4YWmcIaQ4tMd899T4uIiOofhhYREUmDoUXS4b0HieovXtMyguejieTG\na1py4zUtIiKSHkOLiIikwdAiIiJpMLSIiEgaDC2STmJioqW7QEQWwupBI1j5Y714RwwyBasH5cbq\nQSIikh5Di4iIpMHQIiIiaTC0iIhIGgwtkg7vPUhUf7F60AhW/hDJjdWDcmP1IBERSY+hRURE0mBo\nERGRNKQPrf3792Pw4ME4d+6cpbtCVk4IAZ1OB51OB17KJZKT9KG1a9cuBAcHY+fOnZbuCtWRu7n3\noBACI0d+gMjILERGZmHUqAQGF9U5HjjVnNTVg9evX8eECRPw3nvvYfr06fjwww8hhMDSpUuRnp6O\n5s2bw9bWFqGhoejWrRuysrKwdOlSXLlyBW5ubnjttdfg5ORU6bpZ+WO9TLn3oEZz95VjeXl87+8F\n1lY9ePPAKS0tCAAQEXEIycmToCiKJbpn9Yy9f+o67odZ/fTTT+jUqRPc3Nzg5OSErKwslJaW4syZ\nM5g1axby8vIQGxuLXr16obi4GB999BEmTZoEV1dXfP/990hLS8Pjjz9u6c0gI6oKnpqEUk3avRMG\nHhmj1+uRlhaEgoKhAIDU1LJ5AQEBFu6ZXKQOrV27dqF///4AgG7dumHXrl1wcHBA165d0bBhQ7Ro\n0QK+vr4AgHPnzuGvv/7CjBkzAAClpaVo0qSJxfpOd5aXd9GM4SQAvAUg8J9pHYBpAO7uKJfhRGQZ\n0oZWQUEBdDodzpw5A6AshBRFQUREhNHnODo6Iikpqa66SGZQWThoNNUPDSEERoy4gZSUUwCAvn1v\nYNmyizw1Q3VGq9UiPHw90tLKpiMiDkGrjbZspyQkbWjt3bsXvXr1wosvvmiYFx8fD39/f3z11VeI\niIjAxYsXcfz4cfTr1w/Nmzc3PC8kJAQlJSXIycmBp6enpTaB6pCiKFi+/C3o9XoAZTsQBhbVJUVR\nsGzZ5Fs+g9H8DN4FaUNr165dFa5HhYSEYPfu3WjRogUmTpyI5s2bw8fHB/b29lCr1XjjjTewdOlS\nfP7551Cr1ejfvz9DS0J3e+9BRVF4/YAsip/BmpO6etCY69evw87ODrm5uZg0aRLmzJkDZ2fnaq2D\n1YNEcrO26kGqnnuyetCYGTNm4PLly7h+/TqGDx9e7cAiIiLrdE+OtMyBR0lEcuNIS268yzsREUmP\noUVERNJgaJF07ubeg0R0b+A1LSN4Ptp6mXLvQSJe05Ibr2kREZH0GFpERCQNhhYREUmDoUVERNJg\naJF07vbeg0QkP1YPGsHKHyK5sXpQbqweJCIi6TG0iIhIGgwtIiKSBkOLiIikwdAi6fDeg0T1F6sH\njWDlj/XivQfJFKwelNtd/efiYcOG4bPPPjNMb9u2DVlZWRg1apR5e3eL+Ph4XLp0Cba2tgAADw8P\njB8/HgCwf/9+bNiwAYWFhSgpKUFYWBgeffRRLFq0CHv37sWSJUtgZ2cHAFi+fDk2bdqEpUuXwtHR\nEZcvX0bTnnRgAAAZTklEQVRycjLOnj0LAOjcuTOefvppqFT1c7AphIBerwcAaLVaKIpi4R4REd1Z\nlaFVWzuykpIS2NjYGG3ztddeg5eXV7n56enpWLduHSZNmgQXFxcUFhYiJSXF8LiHhwd++uknhIaG\norS0FJmZmdBoNADKdtAJCQmIjIzE66+/jtLSUixevBhr1qzB0KFDa2UbrZkQAiNHfoC0tCAAQETE\nBiQnT2JwEZHVqzK0qnL+/HksXLgQubm5aNq0KV599VW4ublh0aJF6Ny5M7p16wbg/0ZrOp0Oa9eu\nhYODA86ePYsPP/wQH330EU6fPo3i4mIMGjQI3bt3N9re+vXrMXz4cLi4uAAAbG1tER0dDaAs6B58\n8EHs3r0boaGh0Ov18PPzQ3p6OgAgMzMTeXl5CAsLAwCoVCoMGDAAb775Jp5++mnDqO5epNEYO0Uy\ny/Dbhg1D4epa9Xry8niqg6imeIaj5qoMrcLCwnK3zCkoKECXLl0AAN999x2Cg4MxYMAA/Oc//8Gy\nZcvwxhtvVHgTbp3W6/VITExE27ZtkZmZiZKSEsyYMQMAcPXqVQBlb+r8+fMNQdKxY0cMHToUZ86c\nqTD6ulWzZs3w888/48qVK9i1axdCQ0ORnp4OIQTOnDmDDh06lFve09MTTk5OyMnJQatWre74QlmK\n8dCpW7XRDwYh1Sc8w2EeVYaWra0tkpKSDNM3r2kBQEZGBqZOnQoACAsLw3//+987NtamTRu0bdsW\nQFlonDhxAitXrkSfPn0MwWHs9KApQkJCsGvXLpw4cQIvvfSSYb6xD4UMHxZz7NgrBo4A8BaAwH+m\ndQCmAaj89bC2cOG9B0lGer0eaWlBKCgouySRmlo2LyAgwMI9k8tdnx4Eyo4cgPI7//vuuw/FxcUA\ngBs3bqCoqMjw2K3VIC4uLpg5cyb27NmDTz75BL169cLDDz9stK2WLVvi5MmTCAwMrPTxm6cI33zz\nTfTp08fQJ0VR0KJFC6xfv77c8tnZ2cjPz4eHh0c1t9pyzDfaUQBMB3D4n+khMBZY1Wm3rsItLi6u\nTtohIutz16F1//3348cff0T//v2xdetW+Pv7AwB8fX2h1+vRs2dP/PjjjygtLa30+RcvXoSDgwN6\n9+4NGxsb6HS6Ktt7/PHHsWrVKsTFxcHFxQVFRUVISUlBVFQUgLIAdXNzw5AhQ9CxY8dyz+3YsSM0\nGg22bduGPn36oLS0FN9++y369esn1fUsc4ZC2amKlf+cqvgFERGHeKqCqBZptVqEh69HWlrZdETE\nIWi10ZbtlITuunowOjoaixYtwqZNm+Dh4YFXX30VQFkpeUZGBsaNG4fu3bsbStBvX9/vv/+OVatW\nQaVSoXHjxhgxYoThsVuvaTk5OeHtt9/G/fffj8LCQiQlJaGoqAilpaUIDw+vsO6HHnqo0vYmTZqE\n5ORkfPfddwDKQveZZ56pavPvaYqiYNmyybdcFI5mYBHVIv7NmQe/XGwEvxhIJDd+uVhu/NckREQk\nPYYWSYf3HiSqv3h60AgO7a0X7z1IpuDpQbnx9CAREUmPoUVERNJgaBERkTQYWkREJA2GFkmH9x4k\nqr9YPWgEK3+I5MbqQbmxepCIiKTH0CIiImkwtIiISBoMLSIikgZDi6TDew8S1V+sHjSClT/Wi/ce\nJFOwelBurB4kIiLpMbSIiEgaDC0iIpKG2tIdIMsQQkCv1wMAtFotFEWxcI+IiO5M2pGWTqczWkW2\nePFinD17FgDwzTff1GW3pCCEwMiRHyAq6hSiok5h1KgEyFSPw3sPEtVf0lYP6nQ6bNy4EXFxcVUu\nN3z4cKxcubLa65ep8kejqfsqqbw8eV4fqp+ssXqQZzhMZ+z9q9HpwdzcXCQkJCAgIAA6nQ7+/v4Y\nNWoUsrOz8fHHH6OgoABt2rTB6NGj4eDggPj4eAQEBODAgQNo0KABRowYgbZt21a67sTERDz77LNo\n1aoVYmNj8cADD2DQoEH48ssv4ebmhmbNmqGoqAhz587F6dOn8cADD2DIkCEAgPj4eAwfPhx79uxB\nYWEhYmNj0bJlS4wdOxbbt2/Hhg0bUFRUhA4dOuDFF1+syUsAwDKhYWnWss0MT5LFzTMcaWlBAICI\niA1ITp7E4KqmGl/TOnfuHEaNGoVRo0YhISEBx44dw//+9z889thjCAkJwbJly/D1119jxIgRUBQF\nf/31Fz744APs2rULmzZtwujRoytdr5+fH3799Ve4ubnBxsYGx44dAwAcOXIEL730EvLy8qDX6zF7\n9mw0adIEcXFxiIyMhKurq+FD8Nxzz+H7779HUlISACA7OxupqalITEyEjY0NFixYgOPHj8PHx6dG\nr4E17ThNCxMB4C0Agf9M6wBMA1D9Px5r2nYia6bX65GWFoSCgqEAgNTUsnkBAQEW7plcahxaGo0G\nHTp0AFA23NXr9Th58iQmTJgARVHQp08fLFmyxLB8aGgoVCoVAgICqrze5O/vj02bNsHd3R3BwcE4\nfPgwCgsLkZubi2bNmiEvLw/t2rVD8+bNAQDt27fH0aNH8eCDDxpdZ2ZmJnJycvDWW28BAIqKiqDT\n6WocWqawlpFJGQXAdACH/5kegrsJLKButovBSEQ31Ti07O3t/29lajWuX78OoGwoXNmw18HBwbBs\nYWGh0fV6e3vj5MmTcHd3R8eOHfH3338jJSUF3t7eFdZ1c31FRUV37G+nTp2Mju5qkzXueIXwvO10\nxSGeriCqJVqtFuHh65GWVjYdEXEIWm20ZTslIbNXD6rVanh7e2Pfvn0oKSnBtm3b7mr4q1ar4erq\nir1796J9+/bw8/PDxo0b4e/vX631qFQq3LhxAwAQGBiIQ4cOITs7GwBQUFCA8+fPV7tv9wpFUbBs\n2WRs2tQGmza1kSaweO9BkpGsf2/WpsYjrdtfdEVR8NRTT+Hjjz/GqlWr0LZtW6Mjmzu9Yf7+/sjM\nzMR9990HPz8/5OXlGUJLURST3vBHH30UU6ZMgaenJ8aOHYuYmBjMmjULKpUKtra2eOGFF+Dm5mbi\n1t57FEWR7px6UlLSHatGiayRjH9v1kbakvfaJlPJe33DG+aSKayx5J1MxxvmEhGR9Cx+G6eMjAys\nWbOm3Dx3d3dMnDjRQj0iIiJrZfHQCgoKQlBQkKW7QUREEuDpQZIO7z1IVH+xEMMIXkQlkhsLMeTG\nQgwiIpIeQ4uIiKTB0CIiImkwtIiISBoMLZIO7z1IVH+xetAIVv5YL97GiUzB6kG5sXqQiIikx9Ai\nIiJpMLSIiEgaDC0iIpIGQ4ukw3sPEtVfrB40gpU/RHJj9aDcWD1IRETSuydDq6SkxNJdILojIQR0\nOh10Oh14woPINBb/J5CmyM3NRUJCAvz8/KDT6RAUFIQePXpg1apVAIAXX3wRe/fuxcWLF5GdnY0m\nTZrgueeewyeffGIYok+YMAEeHh6YOXMmzp8/D3t7e0RGRqJ79+6W3DSqh26G1ZQpn+Hnn3sBACIi\nNiA5eRIURbFw76i2CCGg1+sBAFqtlu/1XZLimlZubi7Gjh2Ld999F76+vpgwYQJat26N1157DTt3\n7kRWVhYaNWqEXbt2YerUqXBycsJXX30FNzc3hIeHo6SkBCUlJbC1tUVBQQEcHR1x9epVTJw4EUlJ\nSXB0dKzQJs9H0600GsteH7kbeXn1+zNsTde0hBAYOfIDpKWV/Zf2iIhDPEi5A2PvnxQjLaDs1j1a\nrRYA4O3tjcDAQKjVavj6+mLz5s3o3LkzgoOD4eTkBABo164dVq9ejcuXLyMsLAzOzs4AgN27d2Pv\n3r3Iz8/HtWvX8Mcff8DHx8di23Wvqt2dfPw/P1QVGYPWFDKGsV6vR1paEAoKhgIAUlPL5gUEBFi4\nZ/KRJrTs7e0Nv6vVasO0Wq1GUVERgPLJHBwcDC8vL2zfvh3vvPMOxo8fj4YNG2LLli2Ij4+Ho6Mj\nYmNjDc8lmbwH+UNLAHgLQOA/0zoA0wDc+0feMoYOWQ9pQqu6cnNz4e7ujgEDBiAnJwfZ2dlwc3OD\nk5MTHB0dceTIEZw+fdrS3bxn1eaOSaORd8d362kiIfzRtesWTJ36PAIC+kFRLlm6e1RLtFotwsPX\nIy2tbDoi4hC02mjLdkpS0oTW7ed+K5u+dd7u3buxY8cO2NraokWLFujevTtsbGzg5uaGcePGoWXL\nlujQoUOd9J3oJkVRsGzZ5FsuyM/hdY16oOL7Hs33/S5JUYhhCSzEsF781yRkCmsqxKDq45eLiYhI\negwtkg7vPUhUf/H0oBEc2hPJjacH5cbTg0REJD2GFhERSYOhRURE0mBoERGRNBhaJJ3ExERLd4GI\nLITVg0aw8sd68cvFZApWD8qN1YNERCQ9hhYREUmDoUVERNJgaBERkTQYWiQd3nuQqP5i9aARrPwh\nkhurB+XG6kEiIpIeQ4uIiKTB0CIiImkwtIiISBr3VGjl5eVhzpw5d1zul19+wdtvv42pU6fWQa/I\nFEII6HQ66HQ63Kk2iPceJKq/6mX14KxZsxAVFYWAgACjy7Dyp+4IITBy5Af44YeOAIDIyMNITp4E\nRVEqXZ73HiRTWGP1oBACer0eAKDVao1+xsn4+2e20Jo5cybOnz8Pe3t7REZGonv37hg2bBgee+wx\n7N69G61bt8bzzz8PGxsbvPHGG5g3bx5sbGxw9epVxMbGYv78+VCpKg78cnJy8Mknn+DChQto2LAh\nJkyYAHd3d3z22WfYt28fFEXBkCFD8OCDDyI3NxczZszA7NmzsW3bNmRkZODatWs4f/48IiIiEB0d\njbVr1+Lbb7+FRqNBly5dMHTo0Eq3h6FVMxpNbe4wFAAVP7Z5eXzP6P9YW2jdPDhLSwsCAEREHKry\n4Ky+M/b+qc3VwCuvvAJHR0dcvXoVEydORIcOHVBYWAiNRoM5c+Zg8eLFOHDgAHr37g2tVouDBw+i\na9eu2L17N0JCQioNLACYP38+evfujYcffhjFxcUoLS3F6dOncfToUSQlJeHSpUuYMmUKtFpthefq\ndDokJSXBzs4O48ePR2RkJAYNGgSdTodhw4bBy8vLXJtvdWo3NKyTpbaZYUmm0Ov1SEsLQkFB2YFy\namrZvKrO+FBFZgut3bt3Y+/evcjPz8e1a9fwxx9/QKVSoWfPngCAwMBA/Prrr+jduzciIiKwYcMG\ndO3aFdu2bcO//vWvStd57do1nD17FhEREWWdVZd198CBA+jWrRvs7e1hb28Pb29vnDhxAq1atSr3\n/I4dOxrS2tPTE6dOnUK7du3MtclWzdQdqXWEmwDwFoDAf6Z1AKahbERlfRhSRJZjltD6888/sWXL\nFsTHx8PR0RGxsbEoKirCfffdB1tb27KG1GoUFRUBANq3b4+//voLOp0OpaWl8PT0rFZ7iqJUuFhf\n2RDbwcHB8LuNjQ0KCwuru2n3PGvYAQshMGLEDaSknAIA9O17A8uWXazimpZ19JuoOrRaLcLD1yMt\nrWw6IuIQtNpoy3ZKQmYJrYsXL8LJyQmOjo44cuQITp8+fcfn9O7dG/Pnz8egQYOMLtOwYUN4enoi\nJSUFkZGRKCkpgRACnTt3xpIlSxAREYFLly7h5MmTaNeuHW7cuGGOzaE6pigKli9/y+QL1Lz3IMlI\nURQsWzb5ls95NK9n3QWzhJafnx/c3Nwwbtw4tGzZEh06dABQ+ejnpp49e+KLL75Ajx49qlz32LFj\n8cknn+C7776Dg4MDxo8fj1atWqF9+/aIjY2Foih4/vnn4ezsjNzcXH4IJKUoisnn9uPi4mq5N0S1\nozqfc6qcxUred+7ciYyMDIwZM8YSzd8RqweJ5GZt1YNUPbVePVgdycnJOHr0KMaNG2eJ5omISFJW\n8+XipUuX4ujRo+XmRUdHo0+fPhbpD4+SiOTGkZbcav3LxfcafuCI5MbQkhv/nxbdM3jvQaL6iyMt\nI3iUZL1470EyBUdacuNIi4iIpMfQIiIiaTC0iIhIGgwtIiKSBgsxjEhNTbV0F4iohm7+hwhL4D6k\n5ip7/xhaREQkDZ4eJCIiaTC0iIhIGgwtIiKSBkOLiIikYZF/TSKba9euYcGCBcjNzUXTpk0xduxY\n2NnZVVju+vXr+Pe//43ff/8dRUVFeOWVV+Dr61tn7QNAaWkp4uLioNFozPbPEk1p//z581i0aBHy\n8/Ph5OSEPn361OgO/Xq9HitWrEBJSQkiIiIQFRVVYZk1a9bg4MGDaNCgAUaPHo0WLVrcdXvVbX/H\njh349ttvAQCenp544okn0KpVqzpr/6YTJ07g7bffxrhx4xASElKn7Z84cQIrVqzA9evX4eDggPj4\n+Dprv7CwEJ9++il+//13NGzYEI888gi6du1qtvZrw0cffYT09HQ4OTlh9uzZVtGmTqdDUlISmjZt\nCgAICQnBk08+adY+mHvfAEF39Nlnn4n169cLIYT4z3/+I1atWlXpcgsWLBCpqalCCCGKi4vFlStX\n6rR9IYTYuHGjmDdvnkhMTDRL26a2f/HiRfHbb78JIYTIz88XL7zwgjhz5sxdtVdSUiLGjBkj/vzz\nT1FUVCQmTpxYYV0HDhwQH3zwgRBCiGPHjonJkyffVVt32/7Ro0cN7+/WrVvrvP2by8XHx4uEhASx\nZ8+eOm2/oKBAjBs3Tpw/f14IUfae12X733//vViyZIkQQojc3FwxZswYUVpaarY+1Aa9Xi+ysrLE\n+PHjrabNzMxMs+4rKmPKvmH06NEmr4+nB03w888/o3fv3gCAPn364KeffqqwzNWrV3HkyBGEh4cD\nAGxsbGBvb19n7QPAhQsXkJ6ejvDwcAgzfpPBlPZdXFzQpk0bAICTkxO8vb3v+oahJ06cgIeHB9zd\n3aFWq9GjRw/8/PPPRvvk4+ODK1eu4NKlS3fV3t207+vra3h/g4ODceHCBbO0bWr7ALBp0yZ069YN\nTk5OZmvb1PZ37tyJkJAQuLq6AoBZ+2BK+/b29rh27RqKi4tRUFAAW1tbKIpitj7UBn9/fzg4OFhd\nm+bcV1TGnPsGgNe0TJKfnw8XFxcAgLOzM/Lz8yssk5ubCycnJyxatAgTJkzA4sWLUVhYWGftA8CK\nFSswdOhQqFTmfVtNbf+mnJwcZGdnw8fH567ay8vLM+wMgcrv6n77Mq6urma787sp7d8qJSUFXbp0\nMUvbprafl5eHn3/+GZGRkQBg1h22Ke3/8ccfKCgowJQpUxAbG4sdO3bUafs9e/ZEaWkpnn/+eUyZ\nMgWvvfaa2dqvTxRFwbFjxzB+/HgkJCQgOzu7Vtur6b4B4DUtg/fff7/SI/UhQ4aUmza2cygpKcHJ\nkycxcOBAvPjii/j000+xZ88ew2igtts/cOAAnJyc0LZtW+h0OpPaNGf7N12/fh1z585FTEyM0etu\n5lLbR4imyMzMxI4dOzBt2rQ6bXf58uV49tlnoSgKhBB1/lqUlJRAr9fjnXfewY0bNzBt2jSEhITA\n1ta2TtrfvHkzbGxsDNe1EhMTsWjRIrMfsN3r2rZti48//hg2Njb48ccfMWPGDCxYsKBW2rp93/DN\nN99g7969AMr+jUtsbCwAwM/PD6NGjTK6HobWP9555x2jjzk7O+PSpUtwcXHBxYsX4ezsXGEZV1dX\nODo6Go64e/TogR9//NHk0Kpp+0ePHsWBAweQnp6OoqIiXLt2DQsXLsSYMWPqpH0AKC4uxuzZsxEa\nGlqji+Iajabc6bYLFy5Ao9FUe5nabB8ATp8+jU8//RSTJ08262kfU9rPysrC3LlzAQB///03MjIy\noFarzTLiM6V9V1dXBAUFGUbgXl5e0Ov1CAoKqpP29Xo9wsPD0aBBA/j4+KBx48b4448/zFqMUx80\nbNjQ8Ht4eDhWr16NgoICODo6mrWdyvYNAwcOxMCBAwEAr776KpKSkkxaFw9LTNClSxds27YNAPDj\njz9WukN2cXGBh4cHjh8/jtLSUhw8eBAdO3ass/afffZZfPzxx1i0aBFef/11BAQEmBxY5mhfCIHF\nixfD09MT/fv3r1F73t7eyMnJQW5uLoqLi7F79+4KO+MuXbpg+/btAIBjx47BwcHBsAOtKVPaP3/+\nPGbPno2xY8fCw8PDLO1Wp/2FCxdi0aJFWLRoEbp164YXXnjBbKcoTWm/a9eu0Ov1uHHjBgoKCnDq\n1Cn4+fnVWfsdOnTAgQMHUFpaij///BMFBQUMrLtw6dIlwyj9wIEDsLW1NXtgmXPfAPDegyYxVvKd\nl5eHTz75BJMmTQIAnDt3DosWLcLly5fRqlWrKkvTa6P9m/R6PTZu3Ig333yzxm2b2v6RI0fw7rvv\nolWrVoZTiM8+++xdH3nr9XosX77cUPIcHR2NH374AQDQt29fAMDq1atx8OBB2NnZ4ZVXXoGnp6dZ\ntteU9hcvXoz9+/fDzc0NQFnhTUJCQp21f6uPPvoInTt3NnvJ+53a37JlCzZt2oSioiI8+uijePjh\nh+us/atXr+LLL7/EkSNH4OTkhKioKAQHB5ut/dowd+5c/Prrr/j777/h7OyMp59+GmFhYXXS5uXL\nl+Hi4oKnnnoKJSUlAMpex82bN+OHH36ASqVC69atER0dDS8vL7P2wZR9w5gxY7Bw4UKT1sfQIiIi\nafD0IBERSYOhRURE0mBoERGRNBhaREQkDYYWERFJg6FFRETSYGgREZE0GFpERCSN/w85j9572TW8\ngAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x12528d250>"
       ]
      }
     ],
     "prompt_number": 82
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Same model above, but with a Weibull hazard function"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with Model() as model_weibull:\n",
      "    \n",
      "    X = np.c_[[male, pH_std, PCO2_std, PO2_std, AgeYears_std, HoursECMO_std, non_white, mrsa, any_coinf]].T\n",
      "    \n",
      "    # Intercept for survival rate\n",
      "    beta0 = Normal('beta0', 0.0, 0.001)\n",
      "    # Treatment effect\n",
      "    beta = Normal('beta', 0, 0.001, shape=X.shape[1])\n",
      "\n",
      "    # Survival rates\n",
      "    lam = tt.exp(beta0 + tt.dot(X, beta))\n",
      "\n",
      "    alpha = Exponential('alpha', 1)\n",
      "    def weibull_logp(v, t):\n",
      "        # Weibull survival log-likelihood\n",
      "        return tt.sum(v*(tt.log(alpha) + (alpha-1)*tt.log(t*lam) + tt.log(lam)) - (lam*t)**alpha)\n",
      "\n",
      "    surv_like = DensityDist('surv_like', weibull_logp, observed=(died, obs_t))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with model_weibull:\n",
      "    trace_weibull = sample(1000, step=NUTS())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Library/Python/2.7/site-packages/theano/scan_module/scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility\n",
        "  from scan_perform.scan_perform import *\n"
       ]
      },
      {
       "ename": "PositiveDefiniteError",
       "evalue": "Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [ 0  1  2  3  4  5  6  7  8  9 10]",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mPositiveDefiniteError\u001b[0m                     Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-39-e8c4d31d3e9f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mmodel_weibull\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mtrace_weibull\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNUTS\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;32m/Users/fonnescj/Code/pymc/pymc/step_methods/nuts.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model)\u001b[0m\n\u001b[1;32m     66\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 68\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpotential\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mquad_potential\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscaling\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_cov\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mas_cov\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     69\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     70\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mstate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fonnescj/Code/pymc/pymc/step_methods/quadpotential.pyc\u001b[0m in \u001b[0;36mquad_potential\u001b[0;34m(C, is_cov, as_cov)\u001b[0m\n\u001b[1;32m     33\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mQuadPotential_SparseInv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m     \u001b[0mpartial_check_positive_definite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     36\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mis_cov\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mas_cov\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fonnescj/Code/pymc/pymc/step_methods/quadpotential.pyc\u001b[0m in \u001b[0;36mpartial_check_positive_definite\u001b[0;34m(C)\u001b[0m\n\u001b[1;32m     56\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     57\u001b[0m         raise PositiveDefiniteError(\n\u001b[0;32m---> 58\u001b[0;31m             \"Simple check failed. Diagonal contains negatives\", i)\n\u001b[0m\u001b[1;32m     59\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     60\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mPositiveDefiniteError\u001b[0m: Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [ 0  1  2  3  4  5  6  7  8  9 10]"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following incorporates a spline for the age effect."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def interpolate(x0, y0, x):\n",
      "    \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]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nknots = 10\n",
      "\n",
      "with Model() as model_spline:\n",
      "    \n",
      "    X = np.c_[[male, non_white, mrsa, any_coinf]].T\n",
      "    \n",
      "    coeff_sd = HalfCauchy('coeff_sd', 100)\n",
      "    spline = GaussianRandomWalk('spline', sd=coeff_sd, shape=nknots)\n",
      "    knots = np.linspace(AgeYears.min(), AgeYears.max(), nknots)\n",
      "    beta1 = Deterministic('beta1', interpolate(knots, spline, AgeYears))\n",
      "    \n",
      "    # Intercept for survival rate\n",
      "    beta0 = Normal('beta0', 0.0, 0.001)\n",
      "    # Treatment effect\n",
      "    beta = Normal('beta', 0, 0.001, shape=X.shape[1])\n",
      "\n",
      "    # Survival rates\n",
      "    lam = tt.exp(beta0 + tt.dot(X, beta) + beta1)\n",
      "    \n",
      "    def exp_logp(failure, value):\n",
      "        # Exponential survival log-likelihood\n",
      "        return tt.sum(failure * tt.log(lam) - lam * value)\n",
      "\n",
      "    surv_like = DensityDist('surv_like', exp_logp, observed=(died, obs_t))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with model_spline:\n",
      "    trace_spline = sample(1000, step=NUTS())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "PositiveDefiniteError",
       "evalue": "Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [0]",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mPositiveDefiniteError\u001b[0m                     Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-43-2a112baa9e33>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mmodel_spline\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mtrace_spline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNUTS\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;32m/Users/fonnescj/Code/pymc/pymc/step_methods/nuts.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model)\u001b[0m\n\u001b[1;32m     66\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 68\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpotential\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mquad_potential\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscaling\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_cov\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mas_cov\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     69\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     70\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mstate\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fonnescj/Code/pymc/pymc/step_methods/quadpotential.pyc\u001b[0m in \u001b[0;36mquad_potential\u001b[0;34m(C, is_cov, as_cov)\u001b[0m\n\u001b[1;32m     33\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mQuadPotential_SparseInv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m     \u001b[0mpartial_check_positive_definite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     36\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mis_cov\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mas_cov\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/fonnescj/Code/pymc/pymc/step_methods/quadpotential.pyc\u001b[0m in \u001b[0;36mpartial_check_positive_definite\u001b[0;34m(C)\u001b[0m\n\u001b[1;32m     56\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     57\u001b[0m         raise PositiveDefiniteError(\n\u001b[0;32m---> 58\u001b[0;31m             \"Simple check failed. Diagonal contains negatives\", i)\n\u001b[0m\u001b[1;32m     59\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     60\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mPositiveDefiniteError\u001b[0m: Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [0]"
       ]
      }
     ],
     "prompt_number": 43
    }
   ],
   "metadata": {}
  }
 ]
}