31_BasicModelEvaluation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=float:right><img src=\"assets/img/appliedai-logo.png\" width=100 style=\"margin: 0px 20px\"></img></div>\n",
    "\n",
    "\n",
    "##### Jonathan Sedar Personal Project\n",
    "## PyMC3 vs PyStan Comparison\n",
    "_Spring 2016_\n",
    "\n",
    "This set of Notebooks and scripts comprise the **pymc3_vs_pystan** personal project by Jonathan Sedar of Applied AI Ltd, written primarily for presentation at the PyData London 2016 Conference.\n",
    "\n",
    "The project demonstrates hierarchical linear regression using two Bayesian inference frameworks: PyMC3 and PyStan. The project borrows heavily from code written for Applied AI Ltd and is supplied here for educational purposes only. No copyright or license is extended to users.\n",
    "\n",
    "\n",
    "    \n",
    "# 31_BasicModelEvaluation\n",
    "\n",
    "#### Demonstrate posterior predictive checks and info criteria\n",
    "\n",
    "      \n",
    "+ [Setup](#Setup)\n",
    "    + [Local Functions](#Local-Functions)\n",
    "    + [Load Data](#Load-Data)\n",
    "    + [Describe Dataset](#Describe-Dataset)\n",
    "\n",
    "\n",
    "+ [Reload Models and Traces](#Reload-Models-and-Traces)\n",
    "    + [Reload PyMC3 model and traces](#Reload-PyMC3-model-and-traces)\n",
    "    + [Reload PyStan fitted model](#Reload-PyStan-fitted-model)\n",
    "\n",
    "\n",
    "+ [Posterior Predictive Checks on Training Data](#Posterior-Predictive-Checks-on-Training-Data)\n",
    "    + [Calculate Posterior Samples](#Calculate-Posterior-Samples)\n",
    "    + [Eyeballing: Plot the yhat Distributions](#Eyeballing:-Plot-the-yhat-Distributions)\n",
    "    + [Eyeballing: Calculate Coverage](#Eyeballing:-Calculate-Coverage)\n",
    "    + [MSE and R-squared](#MSE-and-R-squared)\n",
    "    + [KS test](#KS-test)  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Interactive magics\n",
    "%matplotlib inline\n",
    "%qtconsole --colors=linux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# filter warnings for presentation's sake\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# general packages\n",
    "import sqlite3\n",
    "from convenience_functions import *\n",
    "from ipywidgets import interactive, fixed\n",
    "from itertools import combinations\n",
    "\n",
    "# scientific packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import patsy as pt\n",
    "from scipy import optimize\n",
    "from scipy.stats import ks_2samp\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# pymc3 libraries\n",
    "import pymc3 as pm\n",
    "import theano as thno\n",
    "import theano.tensor as T \n",
    "import pystan\n",
    "\n",
    "sns.set(style=\"darkgrid\", palette=\"muted\")\n",
    "pd.set_option('display.mpl_style', 'default')\n",
    "plt.rcParams['figure.figsize'] = 12, 4\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# none"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cnxsql = sqlite3.connect('data/car_emissions.db')\n",
    "dfs = pd.read_sql('select * from cars_post_exclusions_2sd', cnxsql, index_col=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## convert sqlite bool storage (as ints) back to bools\n",
    "for ft in ['parent_is_vw', 'mfr_is_vw', 'is_tdi']:\n",
    "    dfs[ft] = dfs[ft].astype(bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2593, 13)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>1653</th>\n",
       "      <th>835</th>\n",
       "      <th>763</th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "      <th>dtype</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>emissions_nox_mgkm</th>\n",
       "      <td>12</td>\n",
       "      <td>36</td>\n",
       "      <td>32</td>\n",
       "      <td>2593</td>\n",
       "      <td>37.32</td>\n",
       "      <td>17.9</td>\n",
       "      <td>1</td>\n",
       "      <td>23.000000</td>\n",
       "      <td>35.000000</td>\n",
       "      <td>51.000000</td>\n",
       "      <td>76</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parent_is_vw</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>bool</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mfr_is_vw</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>bool</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>parent</th>\n",
       "      <td>daimler-ag</td>\n",
       "      <td>bmw</td>\n",
       "      <td>bmw</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>aston</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>volksw</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mfr</th>\n",
       "      <td>mercedes-benz</td>\n",
       "      <td>bmw</td>\n",
       "      <td>bmw</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>abarth</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>volvo</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trans</th>\n",
       "      <td>auto</td>\n",
       "      <td>semiauto</td>\n",
       "      <td>auto</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>auto</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>semiau</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fuel_type</th>\n",
       "      <td>petrol</td>\n",
       "      <td>petrol</td>\n",
       "      <td>petrol</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>diesel</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>petrol</td>\n",
       "      <td>object</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>is_tdi</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>2593</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>bool</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_combined</th>\n",
       "      <td>-0.0728208</td>\n",
       "      <td>0.80692</td>\n",
       "      <td>0.220426</td>\n",
       "      <td>2593</td>\n",
       "      <td>-0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.685973</td>\n",
       "      <td>-0.339409</td>\n",
       "      <td>-0.152797</td>\n",
       "      <td>0.167108</td>\n",
       "      <td>2.75301</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_extra_urban</th>\n",
       "      <td>-0.0751821</td>\n",
       "      <td>0.47462</td>\n",
       "      <td>0.0728415</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.47696</td>\n",
       "      <td>-0.180913</td>\n",
       "      <td>-0.075182</td>\n",
       "      <td>0.093988</td>\n",
       "      <td>21.5997</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_urban_cold</th>\n",
       "      <td>-0.0249258</td>\n",
       "      <td>0.724487</td>\n",
       "      <td>0.291493</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.657763</td>\n",
       "      <td>-0.308037</td>\n",
       "      <td>-0.158155</td>\n",
       "      <td>0.158264</td>\n",
       "      <td>2.92276</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>engine_capacity</th>\n",
       "      <td>-0.259591</td>\n",
       "      <td>0.502214</td>\n",
       "      <td>-0.0383156</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.920116</td>\n",
       "      <td>-0.278857</td>\n",
       "      <td>-0.045471</td>\n",
       "      <td>-0.037765</td>\n",
       "      <td>2.57901</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>emissions_co_mgkm</th>\n",
       "      <td>-0.325803</td>\n",
       "      <td>-0.456874</td>\n",
       "      <td>0.0567836</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.5</td>\n",
       "      <td>-0.846545</td>\n",
       "      <td>-0.382482</td>\n",
       "      <td>-0.088458</td>\n",
       "      <td>0.272874</td>\n",
       "      <td>2.28145</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                             1653       835        763  count   mean   std  \\\n",
       "emissions_nox_mgkm             12        36         32   2593  37.32  17.9   \n",
       "parent_is_vw                False     False      False   2593    NaN   NaN   \n",
       "mfr_is_vw                   False     False      False   2593    NaN   NaN   \n",
       "parent                 daimler-ag       bmw        bmw   2593    NaN   NaN   \n",
       "mfr                 mercedes-benz       bmw        bmw   2593    NaN   NaN   \n",
       "trans                        auto  semiauto       auto   2593    NaN   NaN   \n",
       "fuel_type                  petrol    petrol     petrol   2593    NaN   NaN   \n",
       "is_tdi                      False     False      False   2593    NaN   NaN   \n",
       "metric_combined        -0.0728208   0.80692   0.220426   2593  -0.00   0.5   \n",
       "metric_extra_urban     -0.0751821   0.47462  0.0728415   2593   0.00   0.5   \n",
       "metric_urban_cold      -0.0249258  0.724487   0.291493   2593   0.00   0.5   \n",
       "engine_capacity         -0.259591  0.502214 -0.0383156   2593   0.00   0.5   \n",
       "emissions_co_mgkm       -0.325803 -0.456874  0.0567836   2593   0.00   0.5   \n",
       "\n",
       "                         min        25%        50%        75%      max  \\\n",
       "emissions_nox_mgkm         1  23.000000  35.000000  51.000000       76   \n",
       "parent_is_vw           False        NaN        NaN        NaN     True   \n",
       "mfr_is_vw              False        NaN        NaN        NaN     True   \n",
       "parent                aston         NaN        NaN        NaN   volksw   \n",
       "mfr                   abarth        NaN        NaN        NaN    volvo   \n",
       "trans                   auto        NaN        NaN        NaN   semiau   \n",
       "fuel_type             diesel        NaN        NaN        NaN   petrol   \n",
       "is_tdi                 False        NaN        NaN        NaN     True   \n",
       "metric_combined    -0.685973  -0.339409  -0.152797   0.167108  2.75301   \n",
       "metric_extra_urban  -0.47696  -0.180913  -0.075182   0.093988  21.5997   \n",
       "metric_urban_cold  -0.657763  -0.308037  -0.158155   0.158264  2.92276   \n",
       "engine_capacity    -0.920116  -0.278857  -0.045471  -0.037765  2.57901   \n",
       "emissions_co_mgkm  -0.846545  -0.382482  -0.088458   0.272874  2.28145   \n",
       "\n",
       "                      dtype  \n",
       "emissions_nox_mgkm  float64  \n",
       "parent_is_vw           bool  \n",
       "mfr_is_vw              bool  \n",
       "parent               object  \n",
       "mfr                  object  \n",
       "trans                object  \n",
       "fuel_type            object  \n",
       "is_tdi                 bool  \n",
       "metric_combined     float64  \n",
       "metric_extra_urban  float64  \n",
       "metric_urban_cold   float64  \n",
       "engine_capacity     float64  \n",
       "emissions_co_mgkm   float64  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "custom_describe(dfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Declare feats for use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fts_cat = ['parent_is_vw', 'mfr_is_vw', 'parent', 'mfr', 'trans', 'fuel_type', 'is_tdi']\n",
    "fts_cat_smp = ['mfr_is_vw','trans','fuel_type','is_tdi']\n",
    "fts_num = ['metric_combined', 'metric_extra_urban', 'metric_urban_cold'\n",
    "           ,'engine_capacity', 'emissions_co_mgkm']\n",
    "fts_num_smp = ['metric_combined', 'engine_capacity', 'emissions_co_mgkm']\n",
    "ft_endog = 'emissions_nox_mgkm'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Describe dataset\n",
    "\n",
    "+ The dataset is 2593 rows, with 12 exog features, 1 endog feature.\n",
    "+ These are observations of car emissions tests, one row per car.\n",
    "+ You can read off the basic distributional statistics of the features in the table above. Numeric features have been standardized according to [Gelman's 2sd principle](http://www.stat.columbia.edu/~gelman/research/published/standardizing7.pdf).\n",
    "+ I have selected these particular 12 features to work with. Some are derivatives of original features.\n",
    "\n",
    "We have the following features to choose from:\n",
    "\n",
    "```\n",
    "+ Categoricals:\n",
    "    + `trans`     - the car transmission, simplified to 'auto', 'semiauto', 'manual'\n",
    "    + `fuel_type` - the car power supply, simplified to 'petrol', 'diesel'\n",
    "    + `parent`    - the parent company of the car manufacturer, 20 values\n",
    "    + `mfr`       - the car manufacturer, 38 values\n",
    "\n",
    "+ Booleans:\n",
    "    + `parent_is_vw` - if the parent company of the car manufacturer is Volkswagen\n",
    "    + `mfr_is_vw`    - if the car manufacturer is Volkswagen\n",
    "    + `is_tdi`       - (processed feature) if the car engine type is a turbo diesel\n",
    "    \n",
    "+ Numerics:\n",
    "    + `metric_combined`    - a score for fuel efficiency in combined driving\n",
    "    + `metric_extra_urban` - a score for fuel efficiency in an extra-urban driving\n",
    "    + `metric_urban_cold`  - a score for fuel efficiency in an urban setting, cold start\n",
    "    + `emissions_co_mgkm`  - a count of CO particulates emitted mg/km\n",
    "    \n",
    "+ Numeric endogenous feature:\n",
    "    + `emissions_nox_mgkm` - a count of NOx particulates emitted mg/km    \n",
    "```\n",
    "\n",
    "For the purposes of this Notebook, the final feature mentioned `emissions_nox_mgkm` will be used as the _endogenous_ / _dependent_ / _output_ feature of the linear models. All other features may be used as _exogenous_ / _independent_ / _input_ features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create design matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'emissions_nox_mgkm ~ metric_combined + engine_capacity + emissions_co_mgkm + mfr_is_vw + trans + fuel_type + is_tdi'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fml_all = '{} ~ '.format(ft_endog) + ' + '.join(fts_num_smp + fts_cat_smp)\n",
    "fml_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2593, 9)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>1731</th>\n",
       "      <th>1033</th>\n",
       "      <th>count</th>\n",
       "      <th>mean</th>\n",
       "      <th>std</th>\n",
       "      <th>min</th>\n",
       "      <th>25%</th>\n",
       "      <th>50%</th>\n",
       "      <th>75%</th>\n",
       "      <th>max</th>\n",
       "      <th>dtype</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Intercept</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>1.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mfr_is_vw[T.True]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.20</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trans[T.manual]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.49</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>trans[T.semiauto]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.12</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fuel_type[T.petrol]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>is_tdi[T.True]</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.12</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>metric_combined</th>\n",
       "      <td>-0.312750</td>\n",
       "      <td>-0.472703</td>\n",
       "      <td>2593</td>\n",
       "      <td>-0.00</td>\n",
       "      <td>0.50</td>\n",
       "      <td>-0.685973</td>\n",
       "      <td>-0.339409</td>\n",
       "      <td>-0.152797</td>\n",
       "      <td>0.167108</td>\n",
       "      <td>2.753013</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>engine_capacity</th>\n",
       "      <td>0.042048</td>\n",
       "      <td>-0.278857</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.50</td>\n",
       "      <td>-0.920116</td>\n",
       "      <td>-0.278857</td>\n",
       "      <td>-0.045471</td>\n",
       "      <td>-0.037765</td>\n",
       "      <td>2.579014</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>emissions_co_mgkm</th>\n",
       "      <td>0.099293</td>\n",
       "      <td>-0.421449</td>\n",
       "      <td>2593</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.50</td>\n",
       "      <td>-0.846545</td>\n",
       "      <td>-0.382482</td>\n",
       "      <td>-0.088458</td>\n",
       "      <td>0.272874</td>\n",
       "      <td>2.281452</td>\n",
       "      <td>float64</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         1731      1033  count  mean   std       min  \\\n",
       "Intercept            1.000000  1.000000   2593  1.00  0.00  1.000000   \n",
       "mfr_is_vw[T.True]    0.000000  0.000000   2593  0.04  0.20  0.000000   \n",
       "trans[T.manual]      0.000000  1.000000   2593  0.49  0.50  0.000000   \n",
       "trans[T.semiauto]    0.000000  0.000000   2593  0.12  0.32  0.000000   \n",
       "fuel_type[T.petrol]  0.000000  0.000000   2593  0.50  0.50  0.000000   \n",
       "is_tdi[T.True]       0.000000  0.000000   2593  0.12  0.32  0.000000   \n",
       "metric_combined     -0.312750 -0.472703   2593 -0.00  0.50 -0.685973   \n",
       "engine_capacity      0.042048 -0.278857   2593  0.00  0.50 -0.920116   \n",
       "emissions_co_mgkm    0.099293 -0.421449   2593  0.00  0.50 -0.846545   \n",
       "\n",
       "                          25%       50%       75%       max    dtype  \n",
       "Intercept            1.000000  1.000000  1.000000  1.000000  float64  \n",
       "mfr_is_vw[T.True]    0.000000  0.000000  0.000000  1.000000  float64  \n",
       "trans[T.manual]      0.000000  0.000000  1.000000  1.000000  float64  \n",
       "trans[T.semiauto]    0.000000  0.000000  0.000000  1.000000  float64  \n",
       "fuel_type[T.petrol]  0.000000  1.000000  1.000000  1.000000  float64  \n",
       "is_tdi[T.True]       0.000000  0.000000  0.000000  1.000000  float64  \n",
       "metric_combined     -0.339409 -0.152797  0.167108  2.753013  float64  \n",
       "engine_capacity     -0.278857 -0.045471 -0.037765  2.579014  float64  \n",
       "emissions_co_mgkm   -0.382482 -0.088458  0.272874  2.281452  float64  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(mx_en, mx_ex) = pt.dmatrices(fml_all, dfs, return_type='dataframe', NA_action='raise')\n",
    "custom_describe(mx_ex, 2, )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reload Models and Traces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reload PyMC3 model and traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded object from disk at models/mdl_robust_pymc.pkl\n"
     ]
    }
   ],
   "source": [
    "with pm.Model():\n",
    "    mdl_robust_pymc = read_pickle(relnm='models/mdl_robust_pymc.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## load traces within model context\n",
    "with mdl_robust_pymc:\n",
    "    trc_robust_pymc = pm.backends.text.load('traces/trc_robust_pymc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# alternatively, load with a ref to reloaded model\n",
    "trc_robust_pymc = pm.backends.text.load('traces/trc_robust_pymc', model=mdl_robust_pymc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reload PyStan fitted model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Fit object contains traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded object from disk at models/mdl_robust_stan_model.pkl\n",
      "Loaded object from disk at models/mdl_robust_stan_fit.pkl\n"
     ]
    }
   ],
   "source": [
    "## load the model first, and then the fit:\n",
    "\n",
    "mdl_robust_pystan_model = read_pickle(relnm='models/mdl_robust_stan_model.pkl')\n",
    "\n",
    "mdl_robust_pystan_fit = read_pickle(relnm='models/mdl_robust_stan_fit.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Posterior Predictive Checks on Training Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate Posterior Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyMC3: native method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Feed native method with model and traces, output ppc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 19.5 s, sys: 419 ms, total: 19.9 s\n",
      "Wall time: 20.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "ppc_pymc = pm.sample_ppc(trc_robust_pymc[-1000:], samples=1000, \n",
    "                           model=mdl_robust_pymc, size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Return is a dict of predicted values (traces * calcs * posterior * datapoints)\n",
    "\n",
    "    for each datapoint:              2593x\n",
    "        use each trace value:        1000x\n",
    "            calc the posterior         10x\n",
    "            \n",
    "    creates 10,000 estimates for each data point (25,930,000 values)\n",
    "    \n",
    "    (probably a little overkill)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 10, 1, 2593)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ppc_pymc['y'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2593, 10000)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yhat_pymc = ppc_pymc['y'].reshape(1000*10*1, 2593).T\n",
    "yhat_pymc.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PyStan: remodel using Stan 'generated quantities'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stan can generate a posterior prob given datapoints within the modelsepc itself.\n",
    "\n",
    "Very handy, but appears to be only an option at modelling time, so I'll rerun the model from earlier, this time speccing the `generated_quantities{}`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mdl_robustppc_stan_spec = \"\"\"\n",
    "\n",
    "    /* Spec for robust regression linear model */\n",
    "\n",
    "    data {                                 \n",
    "        int<lower=0> N;                     // count of observations\n",
    "        int<lower=0> M_ppc;                 // count of ppc samples per observation\n",
    "        int<lower=0> N_ppc;                 // count of ppc observations\n",
    "        int<lower=0> K;                     // count of exog features\n",
    "        matrix[N, K] X;                     // exog features\n",
    "        matrix[N_ppc, K] X_ppc;             // exog features for posterior pred\n",
    "        vector[N] y;                        // endog feature\n",
    "    }\n",
    "    parameters {\n",
    "        vector[K] beta;                     // exog coeffs\n",
    "        real<lower=0> sigma;                // linear model error\n",
    "    }\n",
    "    transformed parameters {}\n",
    "    model {  \n",
    "        sigma ~ cauchy(0, 10);              // explicit prior for error (half-cauchy)\n",
    "        y ~ student_t(1, X * beta, sigma);  // student-T likelihood\n",
    "    }\n",
    "    generated quantities {\n",
    "        matrix[M_ppc, N_ppc] y_ppc;\n",
    "        \n",
    "        // needs loops (the *_rng functions only support 1D reals inputs)\n",
    "        \n",
    "        for (m in 1:M_ppc)\n",
    "            for (n in 1:N_ppc)\n",
    "                y_ppc[m, n] <- student_t_rng(1, X_ppc[n] * beta, sigma);  // posterior pred\n",
    "    }\n",
    "    \"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stan_datadictppc = {}\n",
    "stan_datadictppc['N'] = mx_ex.shape[0]\n",
    "stan_datadictppc['K'] = mx_ex.shape[1]\n",
    "stan_datadictppc['X'] = mx_ex.values\n",
    "stan_datadictppc['y'] = mx_en[ft_endog].values\n",
    "\n",
    "stan_datadictppc['M_ppc'] = 10\n",
    "stan_datadictppc['N_ppc'] = mx_ex.shape[0]\n",
    "stan_datadictppc['X_ppc'] = mx_ex.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mdl_robustppc_stan_fit = pystan.stan(\n",
    "        model_code  = mdl_robustppc_stan_spec,\n",
    "        data        = stan_datadictppc,\n",
    "        sample_file = ensure_dir('traces/trc_robustppc_stan/chain'),\n",
    "        model_name  = 'mdl_robustppc_stan',\n",
    "        iter        = 1000,\n",
    "        warmup      = 750, \n",
    "        chains      = 4, \n",
    "        n_jobs      = 2, \n",
    "        verbose     = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## don't print summary this time, since it now also contains 2593 estimates of y_ppc\n",
    "\n",
    "# print(mdl_robustppc_stan_fit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## similarly it's a terrible idea to simply call plot\n",
    "\n",
    "# mdl_robustppc_stan_fit.plot()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Extract ppc from model fit object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ppc_pystan = mdl_robustppc_stan_fit.extract(pars='y_ppc', permuted=True, inc_warmup=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Return array is datapoints * traces\n",
    "\n",
    "    for each datapoint:              2593x\n",
    "        use each trace value:        1000x\n",
    "            calc the posterior         10x\n",
    "            \n",
    "    creates 10,000 estimates for each data point (25,930,000 values)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 10, 2593)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ppc_pystan['y_ppc'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2593, 10000)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yhat_pystan = ppc_pystan['y_ppc'].reshape(1000*10, 2593).T\n",
    "yhat_pystan.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eyeballing: Plot the yhat Distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df_yhat_mean = pd.DataFrame({'observed': mx_en[ft_endog]\n",
    "                       ,'pymc3': yhat_pymc.mean(1).ravel()\n",
    "                       ,'pystan': yhat_pystan.mean(1).ravel()}\n",
    "                       ,index=mx_en.index)\n",
    "df_yhat_mean_melt = pd.melt(df_yhat_mean, var_name='method', value_name='y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Boxplot the distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3wAAAFfCAYAAAAVs86TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGyBJREFUeJzt3X+Q3XV97/HXhrDZ/HQPSZYQEGr8EYmI1cgomsmQADFm\n7BQvDgq3UrSFjkAcU/xBICQkKJYyLVdBpyXtVdMBScpYdbTVTGiiQxmqxNKkUqAMEOiEjSTZSDab\ndSE59w9lL4FsupI9+9397OPxD+zZs+e89/uZc873mfM9+22q1+v1AAAAUJxRVQ8AAABAYwg+AACA\nQgk+AACAQgk+AACAQgk+AACAQo2ueoDh5p577sk73vGOqsfgNyZNmpTnnnuu6jH4DesxtFiPocV6\nDC3WY2ixHkOL9RhaarXaUf28d/gY1o455piqR+AlrMfQYj2GFusxtFiPocV6DC3WoyyCDwAAoFCC\nDwAAoFCCDwAAoFCCDwAAoFCCDwAAoFCCDwAAoFCCDwAAoFCCDwAAoFCCDwAAoFCjqx4AgKHnwxf+\n73TserbqMRgmapOnZu0376h6DAAOQ/AB8Aodu55N7X9dXfUYg6bjW382on7fgdbxrT+regQA+uCQ\nTgAAgEIJPgAAgEIJPgAAgEIJPgAAgEIJPqChLr300qpHAGCAeW6H4UPwAQ21bdu2qkcAYIB5bofh\nQ/ABAAAUSvABAAAUyonXAaAf1l/zoT6/t+DGuwdxEgDoP+/wAQAAFErwAQAAFMohnUDDLViwoOoR\ngAbzOAcYmgQf0HDr16+veoQRqVarpaOj41X9rJ13flvD7XF+NI8PPEfAcOKQTgAAgEIJPgAAgEI5\npBMA+sGpFwAYjrzDBwAAUCjBBwAAUCjBBzTUKaecUvUIAAwwz+0wfAg+oKFWr15d9QgADDDP7TB8\nCD4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCja56AACGntrkqen41p9VPcagGmm/70CqTZ5a\n9QgA9EHwAfAKa795R9UjFKdWq6Wjo6PqMQAYYRzSCQAAUCjBBwAAUCjBBwAAUCjBBwAAUCjBBwAA\nUCjBBwAAUCjBBwAAUCjBBwAAUKhKgm/p0qVV3G2/DOXZAAAAfhve4QMAACjU6EbfwaZNm7Jhw4Yk\nybx583L22Wenu7s7f/VXf5Wnnnoqv/u7v5sLLrggSfLDH/4wGzduTHNzc+bOnZtzzjknSfLggw/m\n7rvvTlNTUz7wgQ/kXe96V+9t/9d//Ve2b9+eer2e5cuXp729PevWrcunPvWpJMmKFSuydOnStLS0\n9Hk7//AP/5Cf/OQnOemkk9Ld3d3oTQIAADAoGhp8e/fuzYYNG7Jq1aokyQ033JC3v/3t6enpySWX\nXJIxY8Zk+fLlOffcc1Or1XL//fdn2bJlmTBhQu9t1Ov13HHHHbnhhhvS3NyclStXZvbs2Rk9+tej\n79mzJ9ddd11Gjfr1m5XTp09PZ2dn9u/fn127duWEE05IS0tLn7fT2dmZzZs358Ybb0x3d3c+/elP\nN3KTAAAADJqGBt+OHTsyY8aM3hh7wxvekGeeeSaTJk1KS0tLkuR1r3tddu3alVqtlsWLF2fDhg3p\n7OzMO9/5zrz5zW/O3r17s3v37tx0001Jkq6uruzevTttbW1JktNPP7339l80Z86c3HfffdmxY0fm\nz5+fJH3eznPPPZcZM2akqakpY8eOzaRJkxq5SQAAAAZNQ4Pv+OOPz+OPP54DBw4kSR599NEsWrQo\nu3fvTmdnZ8aOHZvHH3+895DO4447Luedd16ef/75LFu2LDfddFMmTZqUk046KZ/97GczduzYft3v\nu9/97vzlX/5lkuSiiy5Kkj5vp6WlJU8++WTq9Xr27t2bjo6OgdwEAAAAlWlo8E2cODHnnntuVqxY\nkSSZP39+arVaJkyYkDVr1uTpp5/O3Llzew/h/PrXv55t27Zl//79WbhwYe/tXHjhhbnpppvS1NSU\nWq2WT37yk0e835aWlhx33HE58cQTD7n8cLczadKkzJ49O9dcc01e+9rXZvz48QO8FQAAAKrRVK/X\n61UPMZzcc889ecc73lH1GPxGrVbzruwQYj2GFusxtFiPocV6DC3WY2ixHkNLrVY7qp93WgYAAIBC\nCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4A\nAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBC\nCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4A\nAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBC\nCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4A\nAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBC\nCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4A\nAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBCCT4AAIBC\nCT4AAIBCCT4AAIBCCT4AAIBCjT7SNx966KEj/vCsWbMGdBgAAAAGzhGD77vf/W6SZP/+/dm5c2de\n+9rXJkmeeOKJTJs2LStXrmz8hAAAALwqRwy+q6++Okny5S9/OUuWLElra2uSpL29PevWrWv8dAAA\nALxq/foM37PPPtsbe0kybdq07Nixo2FDAQAAcPSO+A7fi8aNG5d169blrLPOSpLce++9mTBhQiPn\nAgAA4Cj1K/gWL16cdevW5eabb86oUaNy2mmn5ZOf/GSjZwOgAS688KLs2rWz6jEYoSZPnpJvfvPO\nqscAGDH6FXwTJkzIxz/+8UbPAsAg2LVrZ94/5/KqxyjCP937Vdvyt/RP93616hEARpR+BV+SdHZ2\n5tFHH01TU1Pe9KY3Zfz48Y2cCwAAgKPUr+D7+c9/nq9+9auZOXNm6vV6/uZv/iZXXHGF8/ABAAAM\nYf0KvrvuuivXX399pk6dmiTZsWNHbr311nz+859v6HAAAAC8ev06LUOS3thLkuOPPz5NTU0NGQgA\nAICB0a/gGzduXB588MHerx944IGMGzeuYUMxvFx66aVVjwAAwACyf1eOfh3S+Ud/9Ee59dZb89d/\n/ddJkra2tlx55ZUNHYzhY9u2bVWPAADAALJ/V45+BV9bW1tuuOGGdHd3J0laWloaOhQAAABHr9+n\nZUjKDb277rorjzzySEaNGpU/+ZM/SVtbW9UjAQAAHLV+Bd+jjz6a9evXZ9++fYdc/rnPfa4hQw22\nj3zkI0mShx9+ON/+9rdz2WWXVTwRACVZ9X8+3uf3ln/q/w7iJACMNP0Kvq985Sv54Ac/WNk7X0uW\nLMnMmTPz1FNPZfbs2Tn//POzffv2rF27NkuWLEmSrFixIkuXLs3999+frVu35sknn8yCBQvygx/8\nINdee22mTJmSxx57LHfeeWfq9XqmTJmSK6644pD7eeyxx3LiiSdW8SsCAAAMuH5/hu+ss85q8Ch9\n6+npySWXXJLm5uYsX748Z599dqZPn559+/alq6sru3btyrRp03oPOW1ra8vrXve6dHd3573vfW+e\nfPLJTJkyJatXr84111yT17zmNa+4jxUrVmTPnj1ZtWrVYP96AAAADdGv4HvXu96Vn/70pznjjDMa\nPc9hTZo0qTfmZsyYkZ07d6a1tTVz5szJfffdl/b29syfP7/3+q2trUmSiRMnpqOjIwcOHMhzzz2X\n1tbWw8ZekqxcuTKPPfZYbrvttlx77bWN/6UKs2DBgqpHAGCY8JoBMHiOGHwXX3xxmpqaUq/X8/zz\nz6e5uTlJUq/X09TUlG984xuDMuTu3bvT2dmZcePG5YknnsgFF1yQJDnzzDNz8803p16v5w/+4A96\nr1+v1w/5b/LraNyzZ0927dqVyZMnH/Z+WltbnVD+VVq/fn0l91ur1dLR0VHJffNK1mNo6Ws97GxT\ntapeM17K89XQYj2GllqtVtkbPQy8IwbfmjVrBmuOIxo/fnzWrFmTp59+OnPnzs2ECROSJGPGjMnU\nqVNzwgknHHL9F6OtqanpkIC77LLLctttt6Ver2fSpEn50z/90yTJLbfckr1796a5uTkf/3jfH6wH\nAAAYTn6r0zJUZcyYMbn88ssP+73Ozs7Mmzev9+sjfdbw9a9/fVasWPGKy1/8wy8AAAAl6Vfwvfww\nyHq9no0bNx7yubnB9tBDD+Wuu+7Ke97znkycOLGyOQDgf+LUCwBUZVR/rvTlL3/5kK+bmpryL//y\nLw0Z6HC++MUvvuKyWbNmZdWqVVm4cOGgzQEAADCcHPEdvp6envT09OTAgQPp7OzsvXzPnj159tln\nGz4cAAAAr94Rg2/Dhg35/ve/n1/+8pe5+uqre//qZXNzs3fW6HXKKadUPQIAAAPI/l05jhh8ixYt\nyqJFi7J8+XInJKdPq1evrnoEAAAGkP27cvTrM3zLli1r9BwAAAAMsH4F34snXAcAAGD46FfwJckv\nfvGL/Nu//Vvv193d3Q0ZCAAAgIHRr+D78Y9/nC996Uu58847ey873KkSAAAAGDr6deL19evXZ+XK\nlfnCF77Qe9mLf7ETgOFl8uQp+ad7v1r1GMWwLX87kydPqXoEgBGlX8F3zDHHZPTo/3/V7u7u9PT0\nNGwoABrnm9+883++EgOuVqulo6Oj6jEAGGH6FXxvfOMbc8cdd6SrqysPPPBAvvvd72bOnDmNng0A\nAICj0K/P8F100UVpa2vL1KlTc++99+Z973tfPvCBDzR6NgAAAI5Cv4Kvvb099Xo9M2fOzOtf//p0\ndHTke9/7XqNnAwAA4Cj0K/i+8IUv5Jlnnkl3d3e6u7uzf//+7N+/v9GzAQAAcBT69Rm+3/u938tJ\nJ52U3/md32nwOAAAAAyUfgXfgw8+mLVr12b8+PG9p2NoamrKbbfd1tDhAAAAePX6fVqG22+/Pcce\ne2yj5wEAAGCA9Cv4enp68ulPfzrTp08/5PLPfe5zDRkKAACAo9ev4PvgBz/Y6DkAAAAYYP0Kvlmz\nZjV6DgAAAAZYv07LAAAAwPAj+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAA\nAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol\n+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAA\nAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol\n+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAA\nAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol\n+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAA\nAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAApV\nbPBt2LCh6hEAAAAqVWzw3XPPPVWPAAAAUKnRVQ/wUkuWLMnMmTPz1FNPZfbs2Tn//POzffv2rF27\nNkuWLEmSrFixIkuXLk1LS0t++MMfZuPGjWlubs7cuXNzzjnnJEluu+22bN++PStXrsxpp52W888/\nP0mycePG/PSnP82zzz6b+fPn5/3vf3+S5DOf+Uze+ta35pFHHslpp52WCy+8sJoNAAAAMICGVPD1\n9PTkkksuSXNzc5YvX56zzz4706dPz759+9LV1ZVdu3Zl2rRpaWlpSZLcf//9WbZsWSZMmHDI7Vx5\n5ZVZunRpVqxYccjlc+fOzbx58/L888/nmmuu6Q2+zs7O/P7v/34mTpyYz3zmM4IPAAAowpAKvkmT\nJvXG3IwZM7Jz5860trZmzpw5ue+++9Le3p758+f3Xn/x4sXZsGFDOjs78853vjNvfvObj3j7//mf\n/5nNmzdnzJgx6enp6b28tbU1r3nNa5Ikzc3NDfjNAAAABt+Q+gzf7t2709nZmYMHD+aJJ57ItGnT\nkiRnnnlm7r///jzxxBOZOXNm7/WPO+64nHfeefnwhz+cr33ta4fc1gsvvJB6vX7IZV/72tfyh3/4\nh72HfgIAAJRsSL3DN378+KxZsyZPP/105s6d23uo5pgxYzJ16tSccMIJh1z/61//erZt25b9+/dn\n4cKFh3zv9NNPz4033pi2trZceumlSZKZM2fmuuuuyymnnPKKw0ABAABKM6SCb8yYMbn88ssP+73O\nzs7MmzfvkMsuueSSPm/rox/96Csuu+yyyw573S9+8YuH/X8AAIDhbEgF3+E89NBDueuuu/Ke97wn\nEydOrHocAACAYWNIBd/h3l2bNWtWVq1aVcE0AAAAw9uQ+qMtAAAADBzBBwAAUCjBBwAAUCjBBwAA\nUCjBBwAAUCjBBwAAUCjBBwAAUKghdR4+ABrvogs+nJ17OqoegxFoSmstd65bW/UYACOK4AMYYXbu\n6cjF2+tVjzHsrJneZLsdpTXxDw0Ag80hnQAAAIUSfAAAAIUSfAAAAIUSfAAAAIUSfAAAAIUSfAyI\nSy+9tOoRAAAYIPbtyiH4GBDbtm2regQAAAaIfbtyCD4AAIBCCT4AAIBCja56AABg6Fq8fXOf37t1\n+uxBnASAV8M7fAAAAIUSfAAAAIVySCcDZsGCBVWPAMAQ57UCYHAJPgbM+vXrB/0+a7VaOjo6Bv1+\nOTzrMbT0tR52uKlSFa8Vh+P5amixHkNLrVbLGWecUfUYDBCHdAIAABRK8AEAABTKIZ0AQJ+cegFg\nePMOHwAAQKEEHwAAQKEEHwPilFNOqXoEAAAGiH27cgg+BsTq1aurHgEAgAFi364cgg8AAKBQgg8A\nAKBQgg8AAKBQgg8AAKBQTrwOMMJMaa1lTTqqHmNYWjO9qeoRhrUprbWqRwAYcQQfwAhz57q1VY8w\nItVqtXR0CG0ABpdDOgEAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAAAAol+AAA\nAAol+AAAAAo1uuoBhqOf/exnVY8AAACMEGefffar/tmmer1eH8BZAAAAGCIc0gkAAFAowQcAAFAo\nwQcAAFAowQcAAFAowQcAAFCoY66//vrrqx5iuHjhhReyePHijB49Om94wxuSJFu3bs1XvvKVbNq0\nKVOnTk1bW1vFU44Mt99+e773ve9l06ZNmTVrVsaPH5/EelTJtq/W4R4T1qR6L3/d2LJlizWpyO7d\nu/Pnf/7n2bhxY55++um87W1vsx4V2rhxY1avXp1Nmzalra0tbW1t1mMQPfzww/mLv/iLtLe35/TT\nT0/S9+u415LGO9x6DOi+bp1++/73v1+/+eab6z/4wQ/q9Xq9fvDgwfqyZcvqv/rVr+q/+tWv6suX\nL694wpFn69at9dtvv71er1uPKtn2Q8fWrVvrq1evtiZDxEtfN6xJtW655Zb6I4880vu19ajWVVdd\nVT9w4EB937599WuvvdZ6DLItW7bU//Vf/7X+d3/3d/V6ve/Hg3UZHC9fj5caiH1dh3T2U09PT7Zs\n2ZIzzjij97JnnnkmJ5xwQpqbm9Pc3Jzjjz8+7e3tFU458owdOzbHHntsEutRJdt+6Bg7dmxGjx5t\nTYaAl79uWJPqHDx4MDt27Mib3vSm3susR7VOPvnkbN26NQ888EDe9ra3WY9B9ta3vjUTJkzo/bqv\n7W9dBsfL1+OlBmJfd/SATluALVu25Dvf+U6amppSr9fT1NSUiy++OD/72c+ycOHC7Nmzp/e6nZ2d\nGTduXL7xjW+kXq9n3Lhx2bt3b6ZNm1bhb1CWvtbj5JNPTpL88z//cxYtWpTEelTJth86XnxMWJPq\n/eM//uMhrxvWpDrPPfdcenp6cvPNN2f//v1ZuHBhWltbrUeFTj311PzoRz/KwYMHM2fOHI+PivW1\n/V/8f+tSnYHY1xV8L3P66af3Hjv7oq6urjz88MM577zzsmnTpt7LJ0yYkK6urvzxH/9xkmT16tWZ\nOHHiYI5bvMOtx4s2b96c6dOn58QTT0xiPapk2w8NL31MbN++3ZpU6OWvG/V63eOkQhMnTsz48eNz\n1VVX5eDBg7nuuuvyiU98wnpUpL29Pf/xH/+RJUuWJElWrlyZj33sY9ajQn09Px08eNC6VGig9nUF\nXz88/PDDef755/OlL30pv/jFL3Lw4MG85S1vyfTp0/PMM88kSer1etrb2/2LxyB5/PHH8/Of/zwX\nX3xx72XTpk2zHhWx7av38seENanW4V43Tj31VGtSkWOOOSaTJ0/Onj17ctxxx+XYY4/1GKlQvV5P\nV1dXkl//YaN9+/ZZj4rU6/Ukfb9mHDx40LoMohfXIxnYfd2m+ktvmf/Rj370o3R3d+d973tfkl8f\ncvj3f//3aWpqyoc+9KE+341iYF155ZWZPHlyRo0alZNPPjkf+9jHkliPKtn21TrcY+Lf//3fc/fd\nd1uTir30dcOaVGfnzp1ZvXp1urq6cuaZZ2bRokXWo0Lf+ta3snnz5iTJueeem7POOst6DKJvf/vb\nefDBB/PLX/4yp556ai677LI+t7/X98Y73HoM5L6u4AMAACiUv9IJAABQKMEHAABQKMEHAABQKMEH\nAABQKMEHAABQKMEHAABQKMEHAABQKMEHAABQKMEHAA3S3d2dT3ziE3nhhReSJAcOHMgVV1yRrq6u\niicDYKQQfADQIC0tLXn729+en/zkJ0mSzZs35y1veUvGjRtX8WQAjBSCDwAa6JxzzsnGjRuTJBs3\nbsy5555b8UQAjCSCDwAaaMaMGdm3b18ee+yxdHR05I1vfGPVIwEwggg+AGiw+fPn55Zbbsn8+fOr\nHgWAEUbwAUCDvfe9701XV1fmzp1b9SgAjDCCDwAabMuWLZk3b15aWlqqHgWAEWZ01QMAQKn++7//\nO3/7t3+b0aNH56qrrqp6HABGoKZ6vV6veggAAAAGnkM6AQAACiX4AAAACiX4AAAACiX4AAAACiX4\nAAAACvX/AH3WUWUxw6xvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11e3e66a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = sns.factorplot(x='y', y='method', hue='method',\n",
    "                   data=df_yhat_mean_melt, kind='box', size=5, aspect=2.5,\n",
    "                   legend=False, showmeans=True, sym='')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eyeballing: Calculate Coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ppc_coverage(yhat, mx_en, cr=.5):\n",
    "    \"\"\" calc the proportion of coverage from full yhat ppc (nobservations x nsamples)\"\"\"\n",
    "\n",
    "    lower_bound = np.percentile(yhat, 50 - (50 * cr), axis=1)\n",
    "    upper_bound = np.percentile(yhat, 50 + (50 * cr), axis=1)\n",
    "    ninliers = np.sum((mx_en.values.T >= lower_bound).T & (mx_en.values.T <= upper_bound).T)\n",
    "    \n",
    "    return ninliers / mx_en.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dfcov = pd.DataFrame({'cr': np.arange(0, 1.01, .05)})\n",
    "dfcov['coverage_pymc'] = dfcov['cr'].apply(lambda x: ppc_coverage(yhat_pymc, mx_en, cr=x))\n",
    "dfcov['coverage_pystan'] = dfcov['cr'].apply(lambda x: ppc_coverage(yhat_pystan, mx_en, cr=x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dfcov_melt = pd.melt(dfcov, id_vars='cr', var_name='model', value_name='coverage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGnCAYAAACqxPZNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VPW9//H3mZnMJJMJZDMJawOKggYiyNrYeu2vWqX3\nh1CLBLd7vVLbapfr1uVaevVHse29Fatw26tUK9aFxfZ2t9ZbKkqkLKKBiKACIwghgSyQySSZ7fz+\nQFIg22QyycnMvJ6Ph48HzDln5vOdIB/eZ/l+DdM0TQEAAAAA+sxmdQEAAAAAkCwIWAAAAAAQJwQs\nAAAAAIgTAhYAAAAAxAkBCwAAAADihIAFAAAAAHFCwEJSe/LJJ7VixYqo973++us1efLkfq4KAJBI\n6CWpbffu3dqwYYPVZSCBELCAj/zLv/yLnnvuORmGYXUpAIAERS9JPu+8845effVVq8tAAnFYXQBw\nti1btuixxx5Tenq6wuGwLr30Uj355JN64oknNGbMGLW2tmrp0qV67733FA6HNWfOHN10003tx3/v\ne9/T1q1bVVRUpNzcXI0YMaJ929tvv60f/vCHMk1TQ4cO1ZIlS5STk2PFMCVJra2tWrZsmXbs2CG7\n3a5x48bp/vvvb9/W1TgfeeQRZWZmatGiRZKkX/ziF6qpqdE999wjqftxHjp0SF/60pf06U9/Whs3\nblRmZqaeeuopSZLP59PSpUtVU1Ojw4cP6zOf+YzuvPPO9nr/4z/+Q1u3bpUkpaena9q0afra177W\n42d2Z8WKFTp06JCOHj2qo0ePaurUqVq8eLEkafny5XK5XLrtttskSU899ZTq6+t11113acWKFTpy\n5IjeeecdXXHFFTp8+LDefPNN/eY3v5HNZuv2uwWQ/Ogl97dvo5d03UskaePGjVqxYoUcDodCoZAe\neuih9p/3s88+q6efflotLS3avXu3Zs2apa985Ss9jnPLli366U9/qvHjx6uqqkp+v19PPPGEsrOz\ne/vjRSIygUFm8+bN5pVXXmm2tbWZ06ZNMzdt2mQ++OCD5nPPPWeapmkuW7bM/OEPf2iapmm2traa\n1113nblp0ybTNE3zT3/6k3nzzTebkUjEDIVC5m233WYuX77cNE3TDAQC5pw5c8yampr2fb/97W93\n+PzJkyf3qt5HH33UvOmmm87476677orq2AceeMBctmxZp9u6G+eBAwfMOXPmtO977bXXmvv27TNN\ns+dxfvjhh2ZJSYn54osvdvq5DQ0NpmmaZktLi3nppZe2v8+uXbvMG2+80TRN09y6dau5cOHC9mOi\n/W47s3z5cvO2224zw+GwGQqFzOuvv9585ZVXTNM0zerqanP27Nnt+15zzTXmhx9+2H7cN77xDfOd\nd94xp0+fbjY2NprXXHON+f7775um2f13CyD50UvMHseZyr3k4MGD7b+fO3euuWvXri7f+1e/+pW5\nZMmSXo1z8+bN5mWXXWbu3bvXNE3T/OY3v2muW7cuqrEg8XEFC4PSmDFj5HQ6lZWVpfHjx2vbtm1q\nbW2VJL322mv68Y9/LElyuVy69tpr9eqrr2rmzJnatm2b5syZI8MwZLfbNWPGDPn9fknSvn37VF1d\nrXvuuUemaUqSnE5nn2v96le/GvOxf/7zn7V+/fpOt3U3zlGjRik7O1t79uxRWlqaXC6XxowZIym6\ncRYXF+uqq67q9HMdDodeeeUVHTp0SC6XS8eOHVNBQYFcLpdaW1sVCoXU2Nio/Pz89mP6+t3OnDlT\nNtvJO5Y/85nPaPv27brssstUVFSkUaNG6a233lJaWpry8/PPOIs8fvx4ZWVlqaCgQEOHDlVWVlb7\nn5PuvlsAqYFeQi/pqpeMHDmy/bjrrrtO3/nOd3T55Zfr6quv1rnnnhv1Z3Y1Tulkjxo7dqwkaeTI\nkWpqaor6fZHYCFgY1E79Bdvd66Zptt/rbrfbu3wvh8OhkSNH6umnn45rjY888og2b958xmsFBQXt\nDa0noVCo0wZiGEaX45Ska6+9Vr/+9a/ldDp17bXXtr/el3Hu2bNH9957rxYuXKjx48crJyenvYax\nY8eqtLRUc+bM0dixY/Wd73wnLp95NtM0lZaW1v77BQsW6IUXXpDT6VR5eXmv3qur7xZAaqGX0Eu6\n6yULFy7U3Llz9dprr+nuu+/W7bffriuvvLJP40RqY5ILJJxPfvKTev755yVJLS0teuGFF/TJT35S\nkjRjxgy9+OKLMk1Tzc3NZzyUOmbMGAUCAb388svtr3X2F2Fv/3L8+te/rueee+6M/6JtiFdccYUe\neeSRTj+zu3FKJ8/ObdiwQevXr9fVV1/d/nq04+zMpk2bdPnll2vhwoXyeDz68MMP249977335PV6\n9bvf/U4rVqxQUVFRXD5Tkv7yl78oGAwqEAjo97//vWbOnNm+7bLLLtObb76prVu36lOf+lTU79nd\ndwsA9JKT6CVSJBJRRkaGrrzySs2ePVuVlZVnbD91ZerserobJ1IbV7AwqHU2C9OXvvQlLV26VAsW\nLFA4HNa8efM0Y8YMSdLll1+u119/XfPmzVNeXp6GDRvWfpzNZtNPfvITLVmyRE888YRsNptmz56t\nG2+88Yz3DwQCuuGGGzRnzhwtWLCgX8f3jW98Qw899JCuu+46uVwujR49Wg8++KAk6Ytf/GKX45RO\n/oU/bdq09sbQ23F2Zvbs2brjjju0efNmjRkzRtOmTWtvKkOHDlV1dbVuvPFG2Ww2OZ1OXXPNNZo7\nd26fPlM6eUbz9ttvV01Nja644gpNnTr1jPH8wz/8g1wuV/utH105/c9Ld98tgNRCL6GXdNdLfvCD\nH2jnzp0yTVN5eXlasmTJGdtnzZqllStXtgep5cuXKz09vdtxIrUZJlEbQBQ2bNigiooKfetb35LN\nZtP69ev1+OOPa/Xq1X163xUrVigzM1O33HJLp9tDoZBuvvlmLV++XHl5eX36LACAteglSAWWXMHa\nvXu3Vq1apYsuuqjHsxKPP/64qqurZZqmbr/99vYHBwEMrBEjRqiqqko333xz+1m+H/zgBz0ed/31\n13c4e3zqGYClS5d2e+xDDz2kiooK3XLLLTREAEgC9BKkAkuuYO3cuVMtLS169913o77sW1VVpU2b\nNukLX/hCP1cHAAAAALGxZJKLiRMnyuPx9OqYjIyMM2aDAQAAAIDBJmFmEVy/fr2uuOKKHvf7y1/+\nMgDVAABSFX0GANCdhJhF8I033tDw4cPPWGC0Ow0NDf1c0cDLyclhXAmEcSUWxpVYcnJyrC4hab9X\nxpU4GFdiYVyJpa99xtIrWNE8/rVv3z69/fbb+uxnPzsAFQEAAABA7Cy5gvXrX/9ab731lo4fP66W\nlhbddtttkk4u2OZyuTRlypT2fZctW6a8vDw98MADGj16dJfTbwIAAACA1SwJWHPnztXcuXM7vD5r\n1qwOr61YsWIgSgIAAACAPkuYSS4AAAAAYLAjYAEAAABAnBCwAAAAACBOCFgAAAAAECcELAAAAACI\nEwIWAAAAAMQJAQsAAAAA4oSABQAAAABxQsACAAAAgDghYAEAAABAnBCwAAAAACBOCFgAAAAAECcE\nLAAAAACIEwIWAAAAAMQJAQsAAAAA4oSABQAAAABxQsACAAAAgDghYAEAAABAnBCwAAAAACBOCFgA\nAAAAECcELAAAAACIEwIWAAAAAMQJAQsAAAAA4oSABQAAAABxQsACAAAAgDghYAEAAABAnBCwAAAA\nACBOCFgAAAAAECcELAAAAACIEwIWAAAAAMQJAQsAAAAA4oSABQAAAABxQsACAAAAgDghYAEAAABA\nnBCwAAAAACBOCFgAAAAAECcELAAAAACIEwIWAAAAAMQJAQsAAAAA4oSABQAAAABxQsACAAAAgDgh\nYAEAAABAnBCwAAAAACBOCFgAAAAAECcOKz509+7dWrVqlS666CLdeOON3e67c+dOrVu3ToZhaP78\n+SopKRmgKgEgfmp9AVV4G1XnDyrPnaay4mwVeJz9ctypY5ojdcq0haP+LABA4upLv4i1N9FnOmdJ\nwAoGg5o3b57efffdbvczTVNr167V4sWLJUlLly4lYAGwXG8b0kZvo9ZU1ihimu2vrd/boPLSQpUV\nZ8f1uNOPcdgdCoVDUX0WAGBwiCX09LVfRHvM2cfRZzpnScCaOHGidu3a1eN+1dXVGjZsmJzOk3+o\nCgsLdeTIERUVFfV3iQDQqd42pFpfoMP+khQxTa2urNG4fHenjTOW42L9LADA4BBL6BnIfkGfic6g\nfgbL5/PJ7XZr1apVeuqpp+R2u9XU1GR1WQBSVE+NpdYX6HBMhbexw/6nH1fhbex0WyzHnX7MVNv7\nGqXaqD4LAGC9WHqM1Pd+Ee0xZx5n6ir7dmWoLarjUo0lV7Ci5fF45Pf7tWjRIknSypUrlZWV1eNx\nOTk5/V2aJRhXYmFciSWacb2076BsNnuXZ6berA2ofFThGa81R+rksHf9V21zxN7pZ8dy3KljZhm7\ntND2in4UvvaM9+jqs9B7yfo9Mq7EwrgSS0/jiqXHSH3rF7055tRxaXa7bra9rI8ZtXo9PIk+0wlL\nA5bZRXI+paioSNXV1e37Rnt7YENDQ1zqG0xycnIYVwJhXIkl2nEdPHZcoXCoy+0fHjuuhgbPGa9l\n2sLdHpNpC3f62bEcl2kL6+Nmpa6zV+jB4Od0xFZ4xnt09VmJZjA072T4Hs+W6v9/JxrGlViiGVcs\nPUaKvV/E0ps8RlC3GC+qSA36fmCegvZ0+kwnLLlF8Ne//rXWrl2r7du36/HHH29/fdOmTdq+ffvf\ni7PZNH/+fC1ZskRLly7V/PnzrSgXACRJee60brfndrK9rDhbNsPodH+bYXR5T30sx5V9bKjOsx3R\n0uDn9aGZH/VnAQCsF0uPkWLsFzH2pk+MSJMh6YfBz6lFrqiPSzWWXMGaO3eu5s6d2+H1WbNmdXht\n0qRJmjRp0kCUBQDdKivO1vq9DZ3et95VYynwOFVeWqjVZ91XbzMMlZcWdvkwcCzHFWS59F7JDaqp\nrJEU/WcBAKwXS4+RYuwXMfamvNw87S65WcHKGqkXx6WaQf0MFgAMJrE2pLLibI3Ld6vC26h6f1C5\nUU67G8txpx/THLGzPgkAJIhYe4zU934Ra2+iz3SOgAUAvRBrQyrwODWvpKDXnxfLcaeOSdZnGQAg\nWcXaY6S+9Yveos90j4AFAL0Ua0OKK9OU8+BrCgyfLjnSra0FABA3g6LHSFKoTc7DmxUY9Qmpi+e1\n0LlBvQ4WAKATpqmM3S/IdeBVGZGuZ4ECACAmwRZ5tj0qx3GvTn+mF9HhChYAJBIzooxdq+U4/oF8\nM+6UmZZpdUUAgCRiBJvl2fqoQkM/ppYLyyWD6zG9RcACgERhRuSuelb25mo1TftXKS3D6ooAAEnE\nCPjk2fpjhXLPV8v4+dwaGCMCFgAkCOeB12Tz16pp6td47goAEHfuqmcVzL9IrefPJVz1AQELABJE\nYNSlCoycJdmZChcAEkGtL/DRdOZ1CTGdefPEm0+ewCNc9QkBC0DKSrTGJ5tdkt3qKgAAUdjobdSa\nj9a0ctgdCoVDWr+3QeWlhV0uGmw5bj2PCwIWgJSUkI0PAJAQan2B9h5zuohpanVljcbluwf3CT30\nCdOCAEg5PTW+Wl/AospOE2qTgi1WVwEAiEGFt7FDjzklYpqq8DYOcEUdGa3HpS5qRN8QsACknEHf\n+EKt8mx7VK4DG6ytAwAQkzp/sNvt9T1s72+2psMa8vqDsjfut7SOZMUtggBSzmBufEbQL8+2RxXO\nGqm2sVdaVgcAIHZ57rRut+f2sL0/2U8clGfbo2q54PMK54y1rI5kxhUsAClnsDY+I+CTZ8vDCg0d\nI/9FN7C4IwAkqLLibNm6mInPZhiWPetrP+6VZ+uj8k8oV2DEDEtqSAV0bwApZzA2PiPYfDJc5Y1X\ny4TrmCIXABJYgcep8tLCDr3GZhgqLy20ZIILe6NXnm0r5C+5UcFhlwz456cSbhEEkHJONb7VZ010\nYWXjMx0Zaj13toJFUwhXAJAEyoqzNS7f/dFyIHbLlwOJZOSq+eIvKJR3gSWfn0oIWABS0mBrfDJs\nnFEEgCRT4HFqXkmBcnJy1NDQYGktpmuIQq4hltaQKghYAFLWYGp8AAAgOfAMFgAAAADECQELAAaY\nzVetzDf+S4qErS4FAJCE0qq3KmPPr6wuI2URsABgANlPfKisLQ8rWHSJZLNbXQ4AIMk4D22S+511\nCgybbnUpKYtnsABggNiPf3ByitwLFyg4bKrV5QAAkozzwKvKeP8Papp+pyKeYVaXk7IIWAAwAOwN\n++TZ/pOT648UXmx1OQCAJOPyrpdr/8tqmnG3IpkFVpeT0ghYADAAnDVvqnniPylUMNHqUgAAySYc\nlKN+j3wz7lbEnW91NSmPgAUAA6Bl/LVWlwAASFb2NDVP+bLVVeAjTHIBAAAAAHFCwAIAAACAOCFg\nAUCcOWp3ymhrsroMAEAyMk05D/1NMiNWV4IuELAAII6ch/6mzKqnZQROWF0KACDZmBG5335WrgMb\npHDA6mrQBSa5AIA4cR7cqIz3fqemaXcqkjXc6nIAAMnEjMi9c5Vs/jo1Tfu65Ei3uiJ0gYAFAHHg\n+uAVpe/7k5pm3KVIZqHV5QAAkkkkrMwdP5cR8Mk37WuS3Wl1RegGAQsA+shRt0eu/X9W04x7WH8E\nABB36e//XkaoRb5L7pDsaVaXgx4QsAAkhVpfQBXeRtX5g8pzp6msOFsFnoE5wxfKPV9NH/83mU7P\ngHweAGBgWdljJKl1zJWSzUG4ShAELAAJb6O3UWsqaxQxzfbX1u9tUHlpocqKs/u/AMMgXAFAkrK8\nx0hSWsbAfA7iglkEASS0Wl+gQ+OTpIhpanVljWp9zLIEAIgNPQaxIGABSGgV3sYOje+UiGmqwtsY\n3w80TRlBf3zfEwAwKA14j5GkUKsUCcf/fTFgCFgAElqdP9jt9voetveKGZH77eeU8c6a+L0nAGDQ\nGtAeI8kINCtryzI5q7fE9X0xsAhYABJanrv7B35ze9geNTMi985fyOY7LP+FC+PzngCAQW3Aeowk\no61Jni0PK5R7vgLDZ8btfTHwCFgAElpZcbZshtHpNpthxOcB5EhY7sqfy9ZaJ9/Ur7G4IwCkiAHp\nMZKMtuPK2rJMwYIStVxwrdTFZyIxELAAJLQCj1PlpYUdGqDNMFReWtj3aXQjYWW+9TPZgs3yXfIV\nyeHq2/sBABJGv/cYSUZro7I2P6TAsKlqPX8u4SoJME07gIRXVpytcfluVXgbVe8PKjeea5QYNgXP\nKVFg+HTWHwGAFNSvPUaS6UhX67mzFRjBbYHJgoAFICkUeJyaV1IQ/zc2DAVGlcX/fQEACaPfeowk\nOdIJV0mGWwQBAAAAIE4IWAAAAAAQJ5bdIrhz506tW7dOhmFo/vz5Kikp6XLfv/71r3r55Zdlt9u1\nYMGCbvcFgFgZwWa5335e/gsXynRmWl0OACDJ2I8fkOvga/JfdD2TWSQxS65gmaaptWvX6jvf+Y7u\nu+8+rVu3rtv9//CHP+h73/uevv3tb+v5558foCoBpBIj4JNny8OKOLNkprmtLgcAkGTsjfvl2fao\ngvkTCFdJzpIrWNXV1Ro2bJiczpOzrxQWFurIkSMqKirqdP/Ro0dr586dOn78uC6++OKBLBVAKmg5\nOUVusGCiWs6fR+MDAMTX0d3yvPFf8k+8WcGCSVZXg35mScDy+Xxyu91atWqVTNOU2+1WU1NTlwFr\nwoQJ2rBhgyKRiC699NIBrhZAMjNaG6WKRxUomqLW8/6RcAUAiCtH3R6pcqWaJ/2LQudcaHU5GACW\nBCyPxyO/369FixZJklauXKmsrKxO9z1y5Iiqqqp05513SpIeeOABTZw4US5X14t95uTkxL/oQYBx\nJRbGlSB2V0hjPqmMi+Ypw+pa+kHS/bwGiWT9XhlXYmFcCeLtTdKldyqrMDnnEEi6n1ccWBKwioqK\nVF1dLenk81jd3R5omqb8fr8kKRQKqbm5WUYPZ5gbGhriW/AgkJOTw7gSCONKIIVlyTkuJenPS4Oj\nmSfr98q4EgfjSiAX/XNyjktJ+vNS3/uMJQHLZrNp/vz5WrJkSfssgqds2rRJLpdLU6ZMkSQNGzZM\nEyZM0H333SdJmj17dvuzWwAAAAAwmFg2TfukSZM0aVLHh/xmzZrV4bXPfe5z+tznPjcQZQEAAABA\nzCwLWADQmVpfQBXeRtX5g8pzp6msOFsFnvhctbafOCjTkaGIOz8u7wcAwOkcR6sUyr1AsqdZXQos\nRMACMGhs9DZqTWWNIqbZ/tr6vQ0qLy1UWXF2n97bftwrz7YV8pfcRMACgBTWXyfynAdeVcbeP6pp\nxj30mRRHwAIwKNT6Ah3ClSRFTFOrK2s0Lt8dcwO0N+yVZ/tP5S+5ScHC0niUCwBIQP11Is/l/Ytc\n3r+oafrdhCvIZnUBACBJFd7GDuHqlIhpqsLbGNP7OurelWf7T9U86Z8JVwCQwno6kVfrC8T0vq59\nL8nlXS/fjLsVyTwnDpUi0RGwAAwKdf5gt9vre9jeGZv/mDLfelzNFy9S6JzkXH8EABCd/jiRl3Z4\ni1wfVqhp5j2KZOT1tUQkCW4RBDAo5Lm7fyA4t4ftnYm489X08W/T9AAA/XIiL1h4sUJ542W6hsRa\nFpIQV7AADAplxdmydbGIuM0wYr43nnAFAJD650Se7E7CFTogYAEYFAo8TpWXFnYIWTbDUHlpYdym\nagcApKb+OpEHnI1bBAEMGmXF2RqX71aFt1H1/qByezt9bqhNcrj6t0gAQEI6dSJv9VkTXUR9Is+M\nSJGQZOeEH7pHwAIwqBR4nJpXUtDr45wHX/voQeNvSl2coQQApLaYT+RFwnLvfFqR9Gy1XjBvYIpF\nwiJgAUh4rg/+Kte+P8s3/U7CFQCgW70+kRcJK7PyCRmhFvlLbui/wpA0CFgAEppr35/lOrDh5Poj\nLO4IAIincFCZb62UzIh8U26X7DFMhIGUQ8ACkLDS3/+jnIc2qWnGPTIzcqwuBwCQTMIBed58TKYt\nTc2Tb5Ns/LMZ0eFPCoCEFXENORmu0odaXQoAINkYNgXzxqvtY5+SbHarq0ECIWABSFiBUZdaXQIA\nIFnZHGobc4XVVSABsQ4WAAAAAMQJAQsAAAAA4oSABWDwMyPKeGetbM1Hra4EAJCEjLYTclc9I4WD\nVpeCJEDAAjC4RcLKrHxS9qbDijCZBQAgzozW48ravEwR1xBmCkRc8KcIwOAVCSnzrZUyIiH5LrmD\n9UcAAHFltNQra8vDCoycpdZzZ1tdDpIEAQvA4BQOKvPNxySbXb7JXyJcAQDiyuY/Js+Wh9X2sX9g\ntkDEFQELwKCUdrRKcrjUPOlfWH8EABB3Lu//qm3Mp9X2scutLgVJhoAFYFAKFk1WsPBiyTCsLgUA\nkIRaJiygx6BfMMkFgMGLxgcA6C/0GPQTAhYAAAAAxAkBC4DljLYm2ZoOW10GACBJ2Y8fkIItVpeB\nFMEzWAD6Ta0voApvo+r8QeW501RWnK0Cj/OMfYy248ra8mO1DZ+htqzhFlUKAEg00fQYSXLUv6fM\nNx9T8+TbFMo934JKkWoIWAD6xUZvo9ZU1ihimu2vrd/boPLSQpUVZ0uSjNaGk+uPDJ+htnOvsqpU\nAECCiabHSJLj2DvKfOtnar54EeEKA4ZbBAHEXa0v0KHxSVLENLW6ska1voBsLXXK+tuP1DayTK3n\nfdaiSgEAiSaaHiNJjqNVJ8PV5C8qlD/BilKRoghYAOKuwtvYofGdEjFNbd5XK8/mZWor/j9qG/uZ\nAa4OAJDIeuoxFd5G2Ru9ytzxlHyX3K5QHleuMLC4RRBA3NX5g91ur2015LvkDkV45goA0Es99Zh6\nf1DhIaPUNOMeRTxFA1QV8HcELABxl+dO63Z7rjtNkayCAaoGAJBMoukxstkJV7AMtwgCiLuy4mzZ\nuljA0WYYZzyADABAb9BjMNgRsADEXYHHqfLSwvYGaFdY0snGV15a2Ok0ugAAROPsHiOd7DP0GAwW\n3CIIoF+UFWdrXL5b7+2p1Kxj/6PfD7tDM8eeQ+MDAPTZqR5T4W3UmLrXVGweUmDql+gxGBQIWAD6\nzfC2/RrXsFrNU27VnHNGWF0OACCJFHicKs/YLpfeVNOMu2RmEK4wOBCwAPQLx9G3lVn5pJon36ZQ\n3gVWlwMASCamqfT3fy9n9TY1zbhHZjrPXWHwIGABiLu02h1y73xavku+rHDOeVaXAwBIJqap9Hd/\nLWftDjXNuFuma4jVFQFnIGABiDsj6JfvkjsUzh5jdSkAgKRjSnbnyXDl9FhdDNABAQtA3AVGzLS6\nBABAsjJsaj3vs1ZXAXSJadoBAAAAIE4IWAAAAAAQJwQsAH3i+uAV2Ru9VpcBAEhGkZAy9vxKRqDZ\n6kqAqBGwAMTMte8lufa/LNOVZXUpAIBkEw4q883HZPNVy7SzxhUSB5NcAOg901T6+3+Qs3qLmmbe\nIzM9x+qKAADJJByQZ/tPZTrS1Vx6q2Tjn6xIHJb9ad25c6fWrVsnwzA0f/58lZSUdLlvfX29li9f\nrkgkonPPPVc333zzAFYKoNYXUIW3Uc2ROmUaIV1jbJSzoeqj9UeGWl0eACDBnd5nhqpF1/pWK+LO\nkX/SP0s2u9XlAb1iScAyTVNr167V4sWLJUlLly7tNmA9/fTTWrhwoc4///yBKhHARzZ6G7WmskYR\n05TD7tCIyGEdd1Rq2wVf1HTCFQCgj87uM5/WFm01nPINm6cywhUSkCUBq7q6WsOGDZPTefJ+2sLC\nQh05ckRFRUUd9o1EIqqpqSFcARao9QXam94pH5gFWhwsl/G2T8VFARV4uC8eABCbzvrMS+HJkiRj\nx1GNO8dDn0HCsSRg+Xw+ud1urVq1SqZpyu12q6mpqdOAdeLECQUCAf3nf/6nWlpadNVVV2n69Ond\nvn9OTnI+D8K4EksyjOulfQdls9nPmA3HYf/7Xxtv1gZUPqpw4AvrB8nw8+pMso7Lasn6vTKuxJIM\n4+qsz9iDldJGAAAgAElEQVTtae2/ps8Mfsk6rr6IKmA1NTXpmWeeUX19ve677z6ZpqmXXnpJV111\nVUwf6vF45Pf7tWjRIknSypUrlZXV+SxkWVlZyszM1N13361IJKLFixfr4osvbr/61ZmGhoaY6hrM\ncnJyGFcCSZZxHTx2XKFwqP33DrvjjN9/eOy4Gho8VpQWV8ny8zpbMo/Lasn6vTKuxJEs46LPJLZk\nHldfRDVN+2OPPabJkycrGAxKkgzD0Ouvvx7zhxYVFam6ulrSyeexuro9UJLsdrvy8vLU2Ngoh8Oh\ntLS0TvcDEH/5GTaNMo51uT3Xzf+PAIDYjXS1aKi6XuOKPoNEFNUVLJ/Pp5kzZ+qll15qf8087V7Z\n3rLZbJo/f76WLFnSPovgKZs2bZLL5dKUKVPaX7vhhhv02GOPye/3a9asWd1evQIQJ+GgPt/yS1XZ\nA1oemt1hs80wVFacbUFhAIBkYLTU69r6n6vZPlnrwx0nO6PPIFFFFbAMwzjj8t+WLVvkdrv79MGT\nJk3SpEmTOrw+a9asDq/l5+fr29/+dp8+D0AvhAPyvPmYTIdD9RfdLNuOujMeQLYZhspLC3nwGAAQ\nE5v/qDxbHlZb8adUZEyV7ayJLugzSGRRBaybbrpJDz74oI4dO6ZvfvObCoVCuvfee/u7NgBWCLXJ\ns/0nMp1Zap50iz5us+u8c4Z8tD6JXZm2sMqKs2l6AICY2JprlLXlYbWMvVqBj12mMknj8t30GSSN\nqALW2LFj9f3vf1+HDx+WzWbT8OHDZbNF9fgWgERimvJs/4ki6TnyT7xZMk7+f17gcWpeSUHSPswK\nABgYRutxZW1eppZxcxQYVdb+On0GySTqadodDodGjx7dn7UAsJphqOWCeQoPGd0ergAAiBfTNUS+\nKV9WOLvY6lKAfhPVv6DWrl3b4bUNGzZo8eLF7bMBAkgO4aHFhCsAQP8wDMIVkl5U/4ratWuXvvvd\n72rZsmWqq6uTJL366qu65ppr9Nxzz/VrgQAAAACQKKIKWIFAQAsXLtQnPvEJ/fznP5ckhUIhTZ06\nVSdOnOjXAgH0IzNidQUAgGRGn0EKivo+oAkTJmjKlCk6dOiQIpFI+38AEpPR0qCsiqUyWuqtLgUA\nkIQcx95R1uYfEbKQcqKa5OLCCy/UAw88oLa2Nk2aNEk/+MEP1NzcrKeffrpPCw4DGHi1voAq39+v\nK2ue1BtDZio77FGB1UUBAJKKo3anMneuUvPkL/JcL1JOVAHrxhtvlNfrldPp1PDhw3Xw4EF5PB7t\n379fn/70p/u7RgBxstHbqFcqd+tbaS/of0JT9XLNBNlq96u8tFBlxdlWlwcASHC1voAO7npds+p/\nqz+cc4POd36Mk3hIOVFP015cXNz+61GjRkmScnJy4l4QgP5R6wvotcpd+re0X+qXoZl6JTJRkhQx\nTa2urNG4fDeLOgIAYrbR2yjvzld1k/2vWhqcK+/BTNk+5CQeUk/U12wbGhq0f/9+7du3T/v27dP2\n7dv7sy4AcVbhbVSecVxrQmXt4eqUiGmqwttoUWUAgERX6wtoTWWNzlGjvh+8Vl6zUNLfT+LV+gIW\nVwgMnKiuYD3//PN65ZVX5HQ6lZWVpdra2vZJLwAkhjp/UJWRMV1ur/cHB7AaAEAyqfA2KmKa+m14\neodtp07izSvhZkGkhqgC1ubNm7V8+XKtX79excXFcrvd+v3vf9/ftQGIozx3Wrfbc3vYDgBAV+p6\nOEnHSTykkqhuEczPz5fT6VRBQYEOHDig0aNH6/Dhw/1dG4A4KivOls0wOt1mMwzujwcAxIyTeMDf\nRRWwcnNz5fP5NGHCBL388st6/vnnWQMLGOQcdbvlOFrV/vsCj1PlpYUdQpbNMFReWsgEFwCAXnF9\n8Nf2tRQ5iQf8XVS3CN5yyy3KyMiQJH31q1/Vjh07dO+99/ZrYQBi5zhapcwdT6l58m1nvF5WnK1x\n+W5VeBtV7w8q152msuJswhUAIHqmqfT3fifnke0KFJ58Hv/USbzVlTWKnLZGKifxkIqiClinwpUk\njR49WqNHj+63ggD0TVrNW3JXPSPflC8rnHNuh+0FHicPGgMAYmOaytjzKzmO7VLTjLtkuoa0b+Ik\nHnBS1OtgARj80qrfkHvX8/JN/YrCQ4utLgcAkExMUxnvrJGjYa980++U6fR02IWTeECUz2AtWbKk\nv+sA0EdGW5My9vxKvmlfJ1wBAOLOcbRKjuPeLsMVgJOiuoLV0tLS33UA6CPTlaUTn7hfsjNTEwAg\n/kLnlKgpbzx9BuhBVFewSktL9be//a2/awHQVzQ9AEB/MQz6DBCFqK5gbdmyRb/97W/1wgsvSJJM\n05RhGPrRj37Ur8UBAAAAQCKJKmB961vf6u86APSSrblGkcxCq8sAACSjcFBGoElmRq7VlQAJJ6qA\ndc455/R3HQCiZZpKf/93SqupVFPZfZIR1Z2+AABEJxyQ542fKOwZppYLF1hdDZBwmKYdSCSmqYx3\n/0eOo1XyXniHXn37mOr8QeWx1ggAIB5CrfK88V+KpOfqg1FzVFFVS58BeimqgNXU1KRnnnlG9fX1\nuu+++2Sapl566SVdddVV/V0fgFNMUxnvrJWj4X39b9Gteua1Y4qYZvvm9XsbVF5aqLLibAuLBAAk\nrGCLsrYtV9gzTH/2fFZr1h+gzwAxiOreoscee0yTJ09WMBiUJBmGoddff71fCwNwpozd6+Ro3Kf9\nF92hZ972ndH0JClimlpdWaNaX8CiCgEACSscVNbWHys0ZJS8xfO1ZsdR+gwQo6gCls/n08yZM2UY\nRvtr5ln/0wHoX4Fh09Q07V/12oeBDk3vlIhpqsLbOMCVAQASnj1NLefPVcuF5ar44AR9BuiDqG4R\nNAxDDQ0N7b/fsmWL3G53vxUFoKNw9hhJUp2/vtv96v3BgSgHAJBkQvkTJEl1PfQR+gzQvagC1k03\n3aQHH3xQx44d0ze/+U2FQiHde++9/V0bgE7kubtf5DG3h+0AAHSHPgP0TVQBa+zYsfr+97+vw4cP\ny2azafjw4bLZmBoa6DemKZ12S+7pyoqztX5vQ6e3b9gMg4ePAQA9o88A/SbqlORwODR69GiNHDmS\ncAX0p3BAmW/8l+wN73e6ucDjVHlpoWxnNUabYai8tJApdAEA3bI11yrrbz+Ugi2dbqfPAH0T1RWs\nZ599VpdffrmGDx/e3/UAqS3UdnL9EddQhYeO6XK3suJsjct3q8LbqHp/ULmsTwIAiILNd0RZW3+s\nlnNnS2kZXe5HnwFiF1XAys3N1fLly+V0OnX55Zdr1qxZcrlc/V0bkFqCLfK8sUIRd4H8E2+SjO6v\nFBd4nJpXUjBAxQEAEp2t6ZCytj6ilvPnKjDy4z3uT58BYhNVwLr66qt19dVX6/Dhw3r11Vd13333\nady4cfriF7/Y3/UBKcEINsuz9VGFh4yW/6KFPYYrAAB6w378gDzblss/Yb6Cw6dbXQ6Q1KIKWKeE\nQiEFg0GFw2GewwLiyOY/qmDeeLWeP7fLh44BAIiVvXG//BctVLBoitWlAEkvqoD129/+Vq+99pqy\nsrJ0+eWXa8GCBXI6uQcXiJfw0GKFhxZbXQYAIEkFPnaZ1SUAKSOqgNXS0qJ7771XBQXchwsAAAAA\nXYkqYC1YsKC/6wAAAACAhBf1M1hVVVV66623JEmTJ0/WRRdd1G9FAcnM5j8q+/EDCg67xOpSAABJ\nyFG7U5GMXEWyRlhdCpCSopqp4sUXX9Tzzz+voqIiFRUV6dlnn9Wf/vSn/q4NSDq25hplbX5IRtBn\ndSkAgCSUdmS7MneukhEOWF0KkLKiuoL12muv6f7772+f2OITn/iE7r//fl111VX9WhyQTGxNh0+u\nPzJujgKjys7YVusLqMLbqDp/UHks5ggAiEHa4S1yv7NOvqlfU3jo6PbX6THAwIoqYNlstjNmDXS5\nXEzTDvSC/cRBebY9qpYLPq/AiBlnbNvobdSayhpFTLP9tfV7G1ReWqiy4uyBLhUAkICcH76ujHd/\no6bp/3rGrYH0GGDgRZWSRo0apWeeeUZHjhzRkSNH9PTTT2v06NE9HwhAioSV+dZK+SeUdwhXtb5A\nh8YnSRHT1OrKGtX6uMUDANA9m++IMt77rZqm33lGuKLHANaIKmDdcsstcjgc+vGPf6xHHnlELpdL\nt9xyS58+eOfOnfrud7+rf//3f1dVVVWP+4dCId1xxx166aWX+vS5wICz2XXi4//W6aQWFd7GDo3v\nlIhpqsLb2N/VAQASXMRTpOOfuF8RT9EZr9NjAGtEdYug0+lUeXm5ysvL4/Khpmlq7dq1Wrx4sSRp\n6dKlKikp6faYP//5zxozZkxcPh8YcI70Tl+u8we7Pay+h+0AAEjqtM/QYwBrWPIgVXV1tYYNGyan\n0ymn06nCwkIdOXKky/0DgYB27NihadOmDWCVQP/Lc6d1uz23h+0AAHSFHgNYI6qA9bOf/azDaz/9\n6U9j/lCfzye3261Vq1bpqaeektvtVlNTU5f7//GPf9RVV10ls4vL3MBgYrQ2RL1vWXG2bIbR6Tab\nYfAAMgDgTKYZdZ+hxwDWiOoWwYMHD3Z4rbsrTj3xeDzy+/1atGiRJGnlypXKysrqdF+/36/du3dr\n7ty5euWVV6IKWTk5OTHXNpgxrgRwcIu09WfSPy6Lalw5OdKiMpt+/jevTv+jbRjSLTOLdcGoc/qx\n2Ngk1c/rNIwLvZGs3yvjGuRMU3rzF1LDB9KwxT2OKxF7jJREP6+zMK7UEVXACofDZ/zeNE0Fg7Hf\nt1tUVKTq6ur29zpy5IiKioo63Xf37t0KBoN65JFHVFtbq0gkopKSEo0cObLL929oiP4KQqLIyclh\nXINcWvVWuXetlW/qVzXE6Yl6XKX5Dv3bP4xWhbdR9f6gctvXKHEMuu8mmX5ep2NciWUwNPNk/V4Z\n1yBmRpSxa40cx73yTf2ashXdn8NE6jFSEv28zsK4Ektf+0xUAeu8887T6tWrde211yocDmvNmjUa\nN25czB9qs9k0f/58LVmyRIZhaP78+e3bNm3aJJfLpSlTpkiSpkyZ0v7rDRs2qLW1tdtwBVjBeWiT\nMvb8T4f1R6JV4HFqXklBP1QGAEh4ZkTuqmdl91Wradq/SmkZvTqcHgMMrKgC1vXXX69f/OIXuuOO\nO2QYhqZNm6abbrqpTx88adIkTZo0qcPrs2bN6vKYyy67rE+fCfQH54eblPHur0+uP+IZZnU5AIBk\nYppy73xatpY6NU37Wpez0gIYPKKepv3WW2/Vrbfe2t/1AAknnDVCTTPuUiSz0OpSAADJxjAUKLpE\nobwLJLvT6moARCGqgAWga+Gho60uAQCQxEIFE60uAUAvWLIOFgAAAAAkI65gAXFU6wuowtuo5kid\nMm3hj2Zq4pYOAACAVEHAAqJlmsrY8yuFhoxWcPi0Dps3ehu1prJGEdOUw+5QKBzS+r0NKi8tZDFH\nAEDPQq3KrHxCLePnK5LZ+ax/nMgDBj8CFhAN01TGO2vkaNir1rGf6bC51hdoD1eni5imVlfWaFy+\nmwYIAOiSEfTLs+1RhbNGKuLO73QfTuQBiYFnsICemBG533725OKO0++U6fR02KXC29ghXJ0SMU1V\neBv7u0oAQIIyAj55tjys0NAx8l90g2R0/OdZTyfyan2BgSoXQA8IWEB3zMjJ9Uc+WtzRTHN3ulud\nP9jt29T3sB0AkJqMthMnw1XeeLVMuE4yjE7340QekDi4RRDohq21QUYkJN/Ur0kOV5f75bnTun2f\n3B62AwBSk6Nxn4KFpWo97/92Ga4kTuQBiYQrWEA3Ihl5ar54UbfhSpLKirNl66Ix2gyDe+MBAJ0K\nFl6s1nFzug1XEifygERCwALioMDjVHlpYYeQZTMMlZcWMsEFAKBPOJEHJA5uEQTipKw4W+Py3R9N\nn2tn+lwAQNycOpG3+qyJLjiRBww+BCzglFCr0o5WKThsasxvUeBxal5JgXJyctTQ0BDH4gAAic7m\nq5YR9Cucc25Mx3MiD0gMBCxAkoItytq2XGHPMAWLLunxXngAAHrDfuJDebY9qpYLro05YEmcyAMS\nAQELKc8INMuz7RGFho5Ry4ULCFcAgLiyHz8gz7bl8l94nYLDplldDoB+RsBCSjPamuTZ+ohCeReo\nZfznCVcAgLiyN+6X543/kv+iGxQsmmx1OQAGAAELKS1zx5MKFpSoddw1hCsAQHwFW+TZ/lM1T/wn\nhQomWl0NgAFCwEJKay69VabTY3UZAIBklJahE5d+lz4DpBgCFlJad02v1hdQhbdRdf6g8txpzNQE\nAOi1rvoMPQZIXgQsoBMbvY1ac9ZaI+v3Nqi8tJDFHAEAfUKPAZKbzeoCgIFiBHxR7VfrC3RofJIU\nMU2trqxRrS/QH+UBABJcNH2GHgMkPwIWUoKt6ZCGbPx/sjUd7nHfCm9jh8Z3SsQ0VeFtjHd5AIAE\n5zz0N2Vt+oEUCXW7Hz0GSH4ELCQ9+/EDytryY/nHf16RrOE97l/nD3a7vb6H7QCA1OI8WKGMPf8j\n35TbJVv3T1/QY4DkxzNYSGp/X3/kegWLpkR1TJ47rdvtuT1sBwCkDucHG5Sx70U1zbhLkczCHven\nxwDJjytYSFr2hr0nw9XEm6MOV5JUVpwtWxdrYtkMgweQAQCSJJd3vdL3v6SmGXdHFa4kegyQCghY\nSFqmc4iaS29VsGBSr44r8DhVXlrYoQHaDEPlpYVMowsAkCSFMwvUNOMeRdznRH0MPQZIftwiiKQV\nyTxHkczom97pyoqzNS7frQpvo+r9QeWyRgkA4Cyhc0piOo4eAyQ3AhbQhQKPU/NKCqwuAwCQhOgx\nQPLiFkEAAAAAiBMCFpJC2uEtcu1/2eoyAADJyIwoY9ca2Rvet7oSAAmAgIWE5/zwdbl3/1LB/Aut\nLgUAkGzMiNxVz8h+4oDCnhFWVwMgAfAMFhKa88Crytj7RzVNv1MRT5HV5QAAkkkkLPfOp2RrOy7f\n1K9JDpfVFQFIAAQsJCyXd71c3v9V0/S7Y54tEACATkXCyqz8mYxQm3yXfEWydz3DX60voApvo+r8\nQeUxIyCQ8ghYSEyhVqXV7pBvxt2KZOR1uyuNDwDQW/YTByUZ8k35smRP63K/jd5GramsUcQ0219b\nv7dB5aWFLBoMpCgCFhKTI12+6f/a4240PgBALMLZxWqefFu3+9T6Ah16jCRFTFOrK2s0Lt/NCT0g\nBTHJBZJWT42v1hewqDIAQDKo8DZ26DGnRExTFd7GAa4IwGBAwELSovEBAPpTnT/Y7fb6HrYDSE4E\nLAx+ZkRpNW9JXYSlrtD4AADRMIJ+Oer29Pq4PHfXz2ZJUm4P2wEkJwIWBjczInfVs0rf95IU6V0g\novEBAHpiBHzybFmmtKNVvT62rDhbNsPodJvNMHjWF0hRBCwMXpGw3DtWyeavVdO0r3c7RW5naHwA\ngO4YbSfk2bJMofwL1XLB53p9fIHHqfLSwg69xmYYKi8tZIILIEUxiyAGp0hYmZVPyAj65Zv61V6H\nK+nvjW/1WRNd0PgAAEZro7K2PKzAsKlqPe8fpS5OyPWkrDhb4/LdqvA2qt4fVC7LgQApj4CFQSlj\n9wtSOCDfJXd0u/5IT2h8AIAOzIg825YrMGKmWs+9us9vV+Bxal5JQRwKA5AMCFgYlFrPvUpmWqZk\n6/sfURofAOAMhk2+S+6QmZFrdSUAkhABC4OS6RpqdQkAgCRGuALQXywLWDt37tS6detkGIbmz5+v\nkpKSLvd9/PHHVV1dLdM0dfvtt6uggKsRAAAAAAYfSwKWaZpau3atFi9eLElaunRptwHrtttukyRV\nVVXpN7/5jb7whS8MSJ0YIKG2k5NYxPiAMQAA3Qq1SQ6X1VUASBGWBKzq6moNGzZMTufJiQYKCwt1\n5MgRFRUVdXtcRkaG0tJYuyiptPmUtflHaj13toJFk3vcvdYXUIW3UXX+oPKYsAIA0JO6vRr66g91\nouzferz9nB4DIB4sCVg+n09ut1urVq2SaZpyu91qamrqMWCtX79es2fPHqAq0d+MtiZp03KF8sYr\nWHhxj/tv9DZqzVlTrq/f26Dy0kLWtAIAdGBv2Cu99Zj8F13fY7iixwCIF0sClsfjkd/v16JFiyRJ\nK1euVFZWVrfHvPHGGxo+fLhGjBjR4/vn5OTEpc7BJqnG1dIgVTwsjZ6p9InXKb2H2wOPnGjVC1Xv\ny2azd1gde13VMU09d5iKhqT3X70xSKqf12kYV2JJ1nFZLVm/16QaV+0u6c3/lmZ9RZ7h3Z/ES8Qe\nIyXZz+s0jCuxJOu4+sKSgFVUVKTq6mpJJ5/H6un2wH379untt9/WzTffHNX7NzQ0xKXOwSQnJydp\nxmW0NChryzIFRsxSxqQFUY3rj1W1CoZCXW+v/GBQTcWeTD+v0zGuxJLM47Jasn6vyTIux7F3lPnW\nz9R88SJlDb+4x3ElWo+RkuvndTrGlViSeVx9YUnAstlsmj9/vpYsWdI+i+ApmzZtksvl0pQpU9pf\nW7ZsmfLy8vTAAw9o9OjRuuWWW6woG/Fis6t17GcUGHWpMqI8pM4f7HZ7fQ/bAQApxOZQ85QvKZQ7\nLqrd6TEA4smyadonTZqkSZMmdXh91qxZHV5bsWLFQJSEAWK6higw6tJeHZPn7n5yk9wetgMAUke0\nweoUegyAeDr7VmNgUCorzpati+e0bIbBA8gAgJjRYwDEEwELCaHA41R5aWGHBmgzDJWXFjKNLgAg\nZvQYAPFk2S2CSA324weUVrNdrefP7fN7lRVna1y+WxXeRtX7g8pljRIASHnOgxtlOjIUHHZJn96H\nHgMgXghY6Df2hn3ybP+J/CU3xu09CzzOQTeTEwDAGq4P/irXvj/LN/3OuLwfPQZAPBCw0C8c9e8p\n883H1DzxnxQqmGh1OQCAJOPa/7JcH7wi34y7FXHnW10OALQjYCHuTl9/JJQ/wepyAABJJv39P8p5\naJOaZtwtMyPX6nIA4AwELMSXaSp930tqnvJFhXLP73K3Wl9AFd5GNUfqlGkLc587ACAqRmuD0o7u\nUNOMe2SmD+1yP/oMAKsQsBBfhiHftK9LXUx3K0kbvY1aU1mjiGnKYXcoFA5p/d4GlZcWMhUuAKBb\nZnqOmmZ+kz4DYNBimnbEXzdNr9YXaG96p4uYplZX1qjWF+jv6gAAiY4+A2AQI2BhQFV4Gzs0vVMi\npqkKb+MAVwQASCb0GQBWI2ChTxx1eyQzEvX+df5gt9vre9gOAEghZkSOund7dQh9BoDVCFiImfOD\nDcrc8XMZbcejPibPndbt9twetgMAUkQkrMzKJ5W+78VencijzwCwGgELMXHt/1+l73/p5BS56TlR\nH1dWnC1bF/fO2wyDh48BAFIkpMy3VsoItcg35XbJiP6fK/QZAFYjYKHXXHv/JNeBV9Q04x5F3Of0\n6tgCj1PlpYUdmp/NMFReWsgUugCQ6sJBZW7/b0mmfJO/JNl7d8WJPgPAakzTjl5xHnxNrkObPlp/\n5ORZwFNrjdT5g8pzp/W41khZcbbG5bs/Wp/EzvokAIB2mTtXSXanmktvlWz2mN6DPgPASgQs9Eqw\naIqCBaUyXUMknbnWyCnRrDVS4HFqXkmBcnJy1NDQ0O91AwASQ8u4OYpk5LWHq96exDuFPgPAKgQs\n9IqZltn+657WGhmX7+ZsIQCgVyKZBe2/jvUkHgBYiWewEDPWGgEA9BcWDAaQqAhY6JoZkSLhLjez\n1ggAoE8ioS43cRIPQKIiYKFzkbDcO34ul/d/u9yFtUYAALEy2o5rSMVS2Ru9nW7nJB6AREXAQkeR\nsDIrfyZboFltH7u8y91YawQAEAujtUFZmx9SYNhUhbOLO92Hk3gAEhUBC2cKB5X55mNSJCTflC9L\n9q4nqWCtEQBAb9la6pS1+SG1jSxT63mf7XI/TuIBSFTMIoi/Cwfk2f7fMh2uj9Yf6fmPx+lrjdT7\ng8rtxRS6AIDUYms+Ks/Wh9VW/Gm1FX+q231PncRbfdZEF5zEAzDYEbDwd2ZEoZxz1Tr2ql4t7nhq\nrREAALoVCaj13M8qMKosqt05iQcgERGw0K621aaK1mmq23akV4s5AgAQjUjWCAWyRvTqGE7iAUg0\nBCxIYjFHAED/q/UFVOFtVJ0/yIk8AEmLgIUeF3Mcl++mAQIA+oQTeQBSBbMIpiij7YQydv9SMiMs\n5ggAiDt7w/tyHtggqecTebW+gBUlAkC/IGClIKO1UVmbH5Jpd0oyWMwRABBXjro98rzxU0Uy8iWJ\nE3kAUgq3CKYYo6VeWVseVmDkx9V67tWSWMwRABA/jqO7lFn5hJon36ZQ3gWSxIk8ACmFK1gpxOY/\npqzNP1Lb6Mvaw5XEYo4AgPhIq92hzB1PynfJl9vDlcSJPACphYCVQtLf+41ax3xGbWM+fcbrpxZz\nPDtksZgjACBq4aDS3/2NfJfcoXDOeWds4kQegFTCLYIpxD/pFsnoPFOzmCMAoE/saWoqu6/TPnPq\nRN7qsya64EQegGREwEolXYSrU1jMEQDQJ930GU7kAUgVBCwAADAgOJEHIBUQsJJUY41XG2qdOtpi\nKo+zhACAOKr1BbT73bf1bqiAHgMAZyFgJaF33t6ukg+e1vvBudpnFkmS1u9tUHlpIQ8SAwD6ZKO3\nUb6dv9en7Dv0m8CNapWTHgMAp2EWwSTjO7hDJR88reXB2e3hSjq5kOPqyhrV+gIWVgcASGS1voDa\ndv5Gl9mq9P8C16lVJ69a0WMA4O8IWEnEUbtT+bue0sPBf9Tb5ugO2yOmqQpvowWVAQASnmmq9a11\nmmHboyXB+WqQ54zN9BgAOImAlSQc9e8qc+cq/dKzQHvMkV3uV+8PDmBVAIBkkb7vRRX539X3gvN1\nXJmd7kOPAQCewUoaoaHFapp+lwIfOKSj9V3ul+tOG8CqAADJIlA0TS/7LlTT/tYu96HHAABXsJKH\n3RPQgvYAAAvYSURBVKlI1nCVFWfLZhid7mIzDB5ABgDEJJJ5jqadO5weAwA9IGAlmQKPU+WlhR0a\noM0wVF5ayDS6AICY0WMAoGfcIpiozIhkdJ6Py4qzNS7frQpvo+r9QeWyRgkAoDfoMQAQMwLWIFfr\nC6jC26jmSJ0ybWGVFWdr1NENsjcfkb/kpi6PK/A4Na+kYAArBQAkorP7zKWjM1X83ioFhk1TcPi0\nTo+hxwBA1ywLWDt37tS6detkGIbmz5+vkpKSuOybTDZ6G7WmskYR05TD7lAoHJJn/580wr1HgY/f\nbXV5AIAEd3afMcKtmvLB/2/v/kKjvNc8gD+TmMlxTHYRi0bdelMoSDV4AhaELnTploqFUtTQvWgL\nQirUunSLB3poa6G0XixelFL2xlwFQU4tp1SWloJgm4tdEWopiYJlu7IXrTM99I+tOq3RzLsXHtN4\nzqnJO/ll/vXzuTPzy+R5fMff4zfzzvv+Z3zz98viNwO/bXZ5AG2pKQEry7I4evRo7N+/PyIiDhw4\n8IuhKc/aTvKny1MzQ++GLIa7/zs2d/1P/P7yzvjX6b7wu0MA6vWXc6YYU/FvPe/GD1kpfv/NP8cL\n1Vqs7JvjSQD4K025yEW5XI7Vq1dHsViMYrEYq1atikqlsuC1neS//u/iLeHqsa7xGOr633jt2nB8\nmy1zM0cAFmT2nPlNTMXvuv8Y32T98R/Xt8a1rMucAahTU97Bunz5cpRKpRgbG4ssy6JUKsWlS5di\nYGBgQWtvWr58+WKW3xBXat/Eku4bh6c7pqMQWfx77V+i2r00lkTElVp3R/QZ0RnH62/RV3vRF3l0\nwt/r7Dnzd1GNz7J/iHey+6K7u/Dnx82ZVqev9qKvX4+mBKy+vr6oVqsxMjISERGjo6PR39+/4LU3\nfffdd2kLboJlXdNxffp6RERcj4g/xD/9+c/XZx7vhD6XL1/eEX38JX21F321l1YY5p3w9zp7zvwp\nSvHH+MeZP998vBP67OR/B/pqH/pqLwudM005RXBgYCDK5XJE3PiMVaVS+cV3pPKs7SRuGAzAYjJn\nABZHU97B6urqiuHh4Xj11Vdnrgx408mTJ6O3tzeGhobmXNvJbt7M8Q+3XOjCzRwBSMOcAVgcTbtM\n++DgYAwODv7V17ds2TLvtZ1u9s0cr9S6Z+6DZegBkII5A5CeGw23uJs3c+zUc1wBaC5zBiCtpnwG\nCwAAoBMJWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAA\nAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkI\nWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAA\nAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkI\nWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIkIWAAAAIksacYPnZycjLfffjsKhUIMDw/Hhg0bbrv+0KFD\nUS6XI8uy2LNnT6xcubJBlQIAAMxfwwNWlmVx9OjR2L9/f0REHDhwYM6AtXv37oiIOHPmTBw7diye\neuqpRa8TAAAgr4afIlgul2P16tVRLBajWCzGqlWrolKpzOt7ly5dGj09PYtcIQAAQH0W9R2siYmJ\nOHbsWBQKhciyLAqFQuzYsSNKpVKMjY1FlmVRKpXi0qVLMTAwMOfznThxIrZt27aYJQMAANRtUQPW\n4OBgDA4O3vK1CxcuRLVajZGRkYiIGB0djf7+/jmf6/Tp07FmzZpYu3btnGs/+eST+goGoC088MAD\nTf355gxAZ1vInGn4Z7AGBgaiXC5HxI3PY1UqlTnfvTp//nycPXs2nnzyyTmfv9lDF4DOZs4AcDuF\nLMuyRv/QiYmJmasI7ty585Z3uU6ePBm9vb0xNDQ087W9e/fGihUroqurK9atWxe7du1qdMkAAABz\nakrAAgAA6ERuNAwAAJCIgAUAAJCIgAUAAJBIw68imNLk5OTMxTKGh4djw4YNt11/6NChKJfLkWVZ\n7NmzJ1auXNmgSucnTz95e2+mPLW2+jGaLe8xuH79ejz77LPxyCOPxEMPPdSgKvPL09e3334bb775\nZtRqtbjrrrvmdaXPZsnT14cffhjHjx+P7u7ueOyxx1r639e5c+dibGws7rnnnnj88cdvu7ad9o08\nfS3mvmHOtMfrxZy5wZxprk6cM2ZMnXtG1qZqtVr20ksvZVevXs2uXr2avfzyy/P+3snJyezQoUOL\nWF1+efpZSO+NVm+trXiMZqunr/feey87ePBg9sEHHzSgwvrk7ev111/PPvvsswZVV7+8fe3bty+b\nnp7Orly5kr3wwgsNqrI+ExMT2alTp7LDhw/fdl077RtZNv++Zku9b5gz7fF6MWd+Zs40T6fOGTPm\nZ3n2jLY9RbBcLsfq1aujWCxGsViMVatWRaVSmdf3Ll26NHp6eha5wnzy9LOQ3hut3lpb8RjNlrev\nqampmJiYiM2bNzewyvzy9FWr1eKrr76Ku+++u8FV5pf3eK1bty4mJyfj448/jk2bNjWw0vw2btwY\nfX19c65rp30jYv59zZZ63zBn2uP1Ys7cYM40V6fOGTPmZ3n2jLY4RXBiYiKOHTsWhUIhsiyLQqEQ\nO3bsiFKpFGNjY5FlWZRKpbh06dKcNy2OiDhx4kRs27atAZXP3+XLl+fdT561zVZvra14jGbL29f7\n778fW7dujYsXLza40nzy9PXDDz/E1NRUHDx4MH788cfYunVr3HvvvU2oem55j9f69etjfHw8arVa\n3HfffQ2udnG0075Rr4XsG+ZM/WubzZy5wZxprl/7nGmnPaNeefaMtghYg4ODt9yMOCLiwoULUa1W\nY2RkJCIiRkdHo7+/f87nOn36dKxZsybWrl27KLXWq6+vb9795FnbbPXU2qrHaLY8fVWr1Th37lw8\n+uij8dFHH0XWwreey9NXf39/LFu2LPbt2xe1Wi32798fmzZtimKx2MiS5yVPX5VKJc6cORPPPfdc\nRES88sorsXHjxujt7W1YvYuhnfaNeix03zBn6l/bbOaMOdMKfu1zpp32jHrk3TPa9hTBgYGBKJfL\nERGRZVlUKpU5U/L58+fj7Nmz8fDDDzeixFzy9FNP782St9ZWPkaz5enr3Llzce3atXjjjTfi+PHj\nMT4+Hl988UUjy523PH11d3fHihUr4uLFi7FkyZKWPtUmT19ZlkW1Wo2IGx8Yv3LlShQKhYbVWq+5\n/kPVTvvGbPP5j+Ji7RvmTHu8XswZc6YVdPqcMWPy7RmFrJV/zTGHiYmJmauV7Ny585bfPp48eTJ6\ne3tjaGho5mt79+6NFStWRFdXV6xbty527drVjLJ/0S/187d6uV3vrSZPX61+jGbL09dN4+Pj8dNP\nP7X01Z3y9PX111/H6OhoVKvV2LJlS0ufbpOnr3feeSdOnz4dEREPPvhg3H///c0oeV7efffd+PTT\nT+P777+P9evXx+7duyOi/feNPH0t5r5hzrTH68Wc+Zk50zydOGfMmPr2jLYOWAAAAK2kbU8RBAAA\naDUCFgAAQCICFgAAQCICFgAAQCICFgAAQCICFgAAQCICFgAAQCICFgAAQCJLml0AMD9TU1Nx5MiR\n+Pzzz6O7uzvuvPPOGBkZmXn8mWeeie3bt8eJEydiamoqnn/++bjjjjuaWDEA7cScgTQELGgThw8f\njmXLlsVrr732i2u+/PLLOHDgQAOrAqBTmDOQhlMEoU2cOnUqduzYcds127dvb1A1AHQacwbSELCg\nTRQKhZienm52GQB0KHMG0hCwoE1s3rw53nrrrciyrNmlANCBzBlIw2ewoE088cQTceTIkXjxxRej\np6cnBgYG4umnn555vFAoNLE6ANqdOQNpFDK/pgAAAEjCKYIAAACJCFgAAACJCFgAAACJCFgAAACJ\nCFgAAACJCFgAAACJCFgAAACJ/D/wPul0OwEv1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1257c9ac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = sns.lmplot(x='cr', y='coverage', col='model', data=dfcov_melt, fit_reg=False,\n",
    "               size=6, scatter_kws={'s':70})\n",
    "_ = [g.axes[0][i].plot((0,1),(0,1),ls='--',color='#FFA555') for i in range(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observe:**\n",
    "\n",
    "The coverage seems fairly poor:\n",
    "+ Slight under-coverage at narrow CR \n",
    "+ Strong over-coverage at wide CR, particularly above 0.8\n",
    "+ The over-coverage at wide CR indicates the model might be overly influenced by outliers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Compute Mean-Squared Error (MSE)\n",
    "\n",
    "Mean-Squared Error of prediction vs observed $$\\frac{1}{n}\\sum_{i=1}^{i=n}(\\hat{y}_{i}-y_{i})^{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_mse(observed, yhat_dist):\n",
    "    return np.sum((observed - yhat_dist.mean(1))**2) / len(observed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### PyMC3 model MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "26140.064509482203"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_mse(mx_en[ft_endog], yhat_pymc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### PyStan model MSE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "143789.00082573877"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_mse(mx_en[ft_endog], yhat_pystan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute R-squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have the MSE anyway, lets calculate the (unadjusted) r-squared ($R^{2}$) value for each model. It's a ratio of the sum of squared errors for the model and the sum of squared errors for a hypothetical model that just predicts the mean of the dataset. \n",
    "\n",
    "$$R^{2} = 1 - \\frac{\\sum e_{model}^{2}}{\\sum e_{mean}^{2}}$$\n",
    "\n",
    "Ideally the $R^{2}$ lives in the range $[0,1]$, but in practice, it's quite easy for a model to perform arbitrily worse than the mean guess, particularly for oddly shaped data, resulting in a possible range (-inf, 1]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_r2(observed, yhat_dist):\n",
    "    \n",
    "    sse_model = np.sum((observed - yhat_dist.mean(1))**2)\n",
    "    sse_mean = np.sum((observed - observed.mean())**2)\n",
    "    \n",
    "    return 1 - (sse_model / sse_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### PyMC3 model R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-80.591101952395036"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_r2(mx_en[ft_endog], yhat_pymc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### PyStan model R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-447.80926065619184"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_r2(mx_en[ft_endog], yhat_pystan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Kolmogorov–Smirnov (KS) Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean values / Max Likelihoods are a clumsy way to summarise data, wasting lots of rich information present in our carefully gathered Bayesian samples.\n",
    "\n",
    "Instead of summarising to calcuate MSE and $R^{2}$, lets compare the _distributions_ of the predicted values $\\hat{y}$ vs the actual values $y$\n",
    "\n",
    "I'll use the [Kolmogorov-Smirnov (KS) Test](https://en.wikipedia.org/wiki/Kolmogorov–Smirnov_test), specifically the 2-sample version, which lets us quantify the distance between the empirical distribution functions of two 'samples': the model predictions in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Using `scipy.stats.ks_2samp`\n",
    "\n",
    "As per the [docs](http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.stats.ks_2samp.html), this is a two-sided test for the null hypothesis that 2 independent samples are drawn from the same continuous distribution.\n",
    "\n",
    "If the K-S statistic is small or the p-value is high, then we cannot reject the hypothesis that the distributions of the two samples are the same. For my purposes here, I'm simply looking for mutual differences between the posterior predictive distibutions for the `pymc3` and `pystan` yhats, and their individual differences with the `observed` distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "observed vs pymc3:\n",
      "Ks_2sampResult(statistic=0.083686849209409941, pvalue=2.2920448851288229e-08)\n",
      "\n",
      "observed vs pystan:\n",
      "Ks_2sampResult(statistic=0.09409949865021211, pvalue=1.8248854176675423e-10)\n",
      "\n",
      "pymc3 vs pystan:\n",
      "Ks_2sampResult(statistic=0.018897030466640929, pvalue=0.73972718941801741)\n"
     ]
    }
   ],
   "source": [
    "for (a,b) in combinations(['observed','pymc3','pystan'], 2):\n",
    "    print('\\n{} vs {}:'.format(a,b))\n",
    "    print(ks_2samp(df_yhat[a].values, df_yhat[b].values))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observe:**\n",
    "\n",
    "The first 3 results printed above show the difference between the `observed` distribution and the posterior predictive distribution for `pymc3`, and `pystan` models\n",
    "\n",
    "+ The KS statistics for the two pairs `observed` vs `pymc3`, and `observed` vs `pystan` are fairly similar, which agrees with what we saw in the mean distrobution plots earlier\n",
    "+ The KS statistic for the final pair `pymc3` vs `pystan` is a little smaller and the p-value far higher, also suggesting that the `pymc3` and `pystan` models have fitted the data quite similarly - which is what we would hope to see!\n",
    "+ This is amore detailled look than was possible when just using the MSE and R2 checks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**&copy; Applied AI Ltd 2016**  \n",
    "<a href='http://www.applied.ai'>applied.ai</a>"
   ]
  }
 ],
 "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": 0
}