_8_zero_inflated_nb.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An example with zero inflated negative binomial regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here: https://en.wikipedia.org/wiki/Negative_binomial_distribution\n",
    "\n",
    "A negative binomial RV models the number of failures in a sequence of i.i.d Bernouilli trials before the a specified (non-random) number of successes (denoted $r$) occurrs. E.g. we can call \"success\" rolling a 6 on a 6-sided die and ask: \"How many failures (non-6s) will occur before we roll the third ($r=3$) 6?\" In this case, the probability of the number of non-6s will be negative binomial.\n",
    "\n",
    "\n",
    "A generalization is to allow for $r$ to assume any positve real value. \n",
    "\n",
    "NB allows for letting mean $\\neq$ variance, unlike Poisson distribution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1) Theory\n",
    "again, https://en.wikipedia.org/wiki/Negative_binomial_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that a **Binomial** RV describes the distribution of the number of successes in $n$ independent Bernoulli trials with success probability $p$. \n",
    "\n",
    "**Example**: probability of getting exactly 2 heads in 6 tosses of a fair coin:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}[X=2] &= \\text{(number of ways of distributing 2 successes into a sequence of 6 trials)}\\times(\\mathbb{P}[\\text{sucess}])^2\\times(\\mathbb{P}[\\text{failure}])^{6-2}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the first part of this, the success can happen on trial 1 and 5, 1 and 2, 3 and 4 etc. So it is a problem of ${n \\choose x}$. But, in each of these 2-out-of-6 sucesses, I need to get exactly a string of 2 sucesses and 4 losses. This happens with probability $(\\mathbb{P}[\\text{sucess}])^2\\times(\\mathbb{P}[\\text{failure}])^{6-2}$. Hence the required probability is\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}[X=2] &= {6 \\choose 2}\\left(\\frac{1}{2}\\right)^2\\left(\\frac{1}{2}\\right)^4 \\\\\n",
    "&= \\frac{6!}{2!4!}\\frac{1}{64} \\\\\n",
    "&= \\frac{6\\times 5}{2}\\frac{1}{64}\\\\\n",
    "&= \\frac{15}{64}.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a **Negative Binomial** RV we are instead asking: how many times ($k$) will we fail until the $r$th success? So with $k$ failures and $r$ successes, there will be a total of $k+r$ Bernouilli trials. But we are interested in the number of failures **BEFORE** the $r$th success. So the total number of trials is $k+r-1$. We are now asking about the way of distributing $k$ failures among these. We therefore have (success probability $p$, $r$=number of successes, $k$=number of failures)\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(k;r,p) &= \\mathbb{P}[X=k] = {k+r-1 \\choose k}p^r(1-p)^k.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now one can extend the negative binomial distribution to the case where $r$ is any positive real, not necessarily integer. In that case:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "{k+r-1 \\choose k}p^r(1-p)^k &= \\frac{(k+r-1)(k+r-2)...(r)(r-1)...(1))}{k!(r-1)!} \\\\\n",
    "&= \\frac{(k+r-1)(k+r-2)...(r)}{k!} \\\\\n",
    "&= \\frac{\\Gamma(k+r)}{k!\\Gamma(r)},\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where the binomial coefficient is extended to its real-values counterpart, the Gamma function (see wikipedia). Hence, the probabiluty mass function becomes:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(k;r,p) &= \\mathbb{P}[X=k] = \\frac{\\Gamma(k+r)}{k!\\Gamma(r)}p^r(1-p)^k.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In negative binomial regression, the distribution is specified in terms of the mean $m = \\frac{pr}{1-p}$, hence $p=\\frac{m}{m+r}$ and $1-p=\\frac{r}{m+r}$. Substituting this:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "f(k;r,p) &= \\mathbb{P}[X=k] = \\frac{\\Gamma(k+r)}{k!\\Gamma(r)}\\left(\\frac{m}{m+r}\\right)^r\\left(\\frac{r}{m+r}\\right)^k.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For maximum likelihood estimation, see: https://en.wikipedia.org/wiki/Negative_binomial_distribution#Maximum_likelihood_estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Zero Inflated Negative Binomial** regression model is then analogous to the ZIP model we saw in `_3_poisson_statsmodesl_zero_inflation`: we use a binary outcome model such as a logistic regression to model whether something does or does not happen and then we use a count model (in this case a negative binomial) to model the expected count of whatever we are interested in conditional on the thing happening. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2) Example: fish caught in state park again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider the same example as we did in 2.2) in `_3_poisson_statsmodesl_zero_inflation`. The example is taken from here: https://stats.idre.ucla.edu/r/dae/zinb/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
      "  import pandas.util.testing as tm\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "plt.style.use('seaborn')\n",
    "mpl.rcParams['font.family'] = 'serif'\n",
    "%matplotlib inline\n",
    "\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.discrete.count_model import ZeroInflatedNegativeBinomialP as ZINB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('https://stats.idre.ucla.edu/stat/data/fish.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nofish</th>\n",
       "      <th>livebait</th>\n",
       "      <th>camper</th>\n",
       "      <th>persons</th>\n",
       "      <th>child</th>\n",
       "      <th>xb</th>\n",
       "      <th>zg</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.896315</td>\n",
       "      <td>3.050405</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.558345</td>\n",
       "      <td>1.746149</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>-0.401731</td>\n",
       "      <td>0.279939</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.956298</td>\n",
       "      <td>-0.601526</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0.436891</td>\n",
       "      <td>0.527709</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nofish  livebait  camper  persons  child        xb        zg  count\n",
       "0       1         0       0        1      0 -0.896315  3.050405      0\n",
       "1       0         1       1        1      0 -0.558345  1.746149      0\n",
       "2       0         1       0        1      0 -0.401731  0.279939      0\n",
       "3       0         1       1        2      1 -0.956298 -0.601526      0\n",
       "4       0         1       0        1      0  0.436891  0.527709      1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nofish</th>\n",
       "      <th>livebait</th>\n",
       "      <th>camper</th>\n",
       "      <th>persons</th>\n",
       "      <th>child</th>\n",
       "      <th>xb</th>\n",
       "      <th>zg</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>250.000000</td>\n",
       "      <td>250.000000</td>\n",
       "      <td>250.000000</td>\n",
       "      <td>250.00000</td>\n",
       "      <td>250.000000</td>\n",
       "      <td>250.000000</td>\n",
       "      <td>250.000000</td>\n",
       "      <td>250.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.296000</td>\n",
       "      <td>0.864000</td>\n",
       "      <td>0.588000</td>\n",
       "      <td>2.52800</td>\n",
       "      <td>0.684000</td>\n",
       "      <td>0.973796</td>\n",
       "      <td>0.252323</td>\n",
       "      <td>3.296000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.457407</td>\n",
       "      <td>0.343476</td>\n",
       "      <td>0.493182</td>\n",
       "      <td>1.11273</td>\n",
       "      <td>0.850315</td>\n",
       "      <td>1.440277</td>\n",
       "      <td>2.102391</td>\n",
       "      <td>11.635028</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>-3.275050</td>\n",
       "      <td>-5.625944</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.008267</td>\n",
       "      <td>-1.252724</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>2.00000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.954550</td>\n",
       "      <td>0.605079</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>4.00000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.963855</td>\n",
       "      <td>1.993237</td>\n",
       "      <td>2.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>4.00000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>5.352674</td>\n",
       "      <td>4.263185</td>\n",
       "      <td>149.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           nofish    livebait      camper    persons       child          xb  \\\n",
       "count  250.000000  250.000000  250.000000  250.00000  250.000000  250.000000   \n",
       "mean     0.296000    0.864000    0.588000    2.52800    0.684000    0.973796   \n",
       "std      0.457407    0.343476    0.493182    1.11273    0.850315    1.440277   \n",
       "min      0.000000    0.000000    0.000000    1.00000    0.000000   -3.275050   \n",
       "25%      0.000000    1.000000    0.000000    2.00000    0.000000    0.008267   \n",
       "50%      0.000000    1.000000    1.000000    2.00000    0.000000    0.954550   \n",
       "75%      1.000000    1.000000    1.000000    4.00000    1.000000    1.963855   \n",
       "max      1.000000    1.000000    1.000000    4.00000    3.000000    5.352674   \n",
       "\n",
       "               zg       count  \n",
       "count  250.000000  250.000000  \n",
       "mean     0.252323    3.296000  \n",
       "std      2.102391   11.635028  \n",
       "min     -5.625944    0.000000  \n",
       "25%     -1.252724    0.000000  \n",
       "50%      0.605079    0.000000  \n",
       "75%      1.993237    2.000000  \n",
       "max      4.263185  149.000000  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5ClZXnn//cMM/bQM9NtD9UwVjSjorn8CSSAm1lAUMbE1GysSPwWJgRlN5aWZinltxokiGhAxGTNGuBbyy7ZdWo1AgF3qajENRSSKcHBEuVbXJaCYBCmx6HpHjjYzIz9/eM87d603dLD+fGcc/r9qqKY5z6nz7nONafu/sx97uc5K2ZnZ5EkSZLUtLLuAiRJkqReYkCWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKq+ouQJLU2yJiA3AdsIHm7433ZebdtRYlSR3kCrIk6dn8LfDtzDwe+BBwc0QM1VyTJHWMAVmStKhq9fj/Aa4FyMyvA08D/67OuiSpk9xiIUltEBGHA/8ZGAaeB3wVuAK4CtgIrAG+B5yZmfsi4iLgfcD1wPOB1wH/ANwAfBg4Ajg7M2+KiLcAnwR+Wj3Ga4Eh4N1zWx2q57+qGl8JXJCZ/1L87E7gTuAE4LDMfMkSX9pvATOZ+VAx9v8Bx1S1StLAMSBLUosi4iDgfwGXZ+bfRcQocC/wn4CvZObnqvtdB7wTuDYzL4mIlwK/DfxbYBR4GNiTmX8QEacAlwE3ZeaXqpXc/wK8NzO/GxF/DNwYES8HZoH/DVyRmf81Io4A/k9EvKT42b8F3pOZ50XEFVU9n6cZ3hfy5cy8DDgMmJp32+PAoa32TZJ6lQFZklr328DLgG0AmTkVEacCk8CmiPgGsA94MdCg2q5QuS0znwZ2RcQE8J1q/B7gpfOe597M/G715y/QPHFuM7AfOBz4H9Xz3xMRD9PcBvE/q/tnZt5X/eG86v9vX+Lrm11gbMUSf1aS+o57kCWpdS8EJjNz39xAZt5Bc7X4PcBbMvMkmoF2eN7P7in+vK843gesnnffyeLx99NcyX1B9fyzwK0R8c8R8c80t1qMFj87fxV4qSbmPQ40t4RMPMfHk6Se5wqyJLXux8BYRKyaC8kR8UrgRODOzHysut/8wHugNsz9ISJW0Qyqj9AM03urED53+1rg57/qwZa4xeJuYE1EvCgzf1zd9iqaYV+SBpIBWZJa903gB8AfA/+92vP798DngHdVl0TbD5wM3NfC8/xGRLy22mbxduAnwHaaQfihiDglM2+swvNNwLn83y0bv2QpWywyc3dEfBH4D8BHI+JEmich3tLC65CknrZidnahrWWSpAMx7yoWK4ELgW/RDMmvoHn1iVU09wx/svqx84GfAWcBb6IZQu8DTgX+Dvg3wK2Z+TsRcQbw74Fv07yCxME0r2Kxo3j+z1bjK4H/Vp2w90aaJ+htpLma/TvP4bWVXxRyEPBnflGIpEFmQJakPlAF5DPKbRSSpM7wJD1JkiSpYECWpB5XfdnHB4GjIuJv6q5HkgadWywkSZKkgivIkiRJUqEnLvO2b9/+2cnJRt1lDLyxsWHsc2fZ4+6wz91hnzvPHneHfe6Ofuzz+Pj6Bb8VtCdWkFetOqjuEpYF+9x59rg77HN32OfOs8fdYZ+7Y5D63BMBWZIkSeoVBmRJkiSp0BN7kKempnj00Z0L3jY8PMzIyGiXK5IkSdJy1RMB+cUveSmPTz624G2jo2Ps2HGPIVmSJEld0RMB+fHJxzjhtCsZWjv2jPGZJye5fds5NBoNA7IkSZK6oicCMsDQ2jHWrNtQdxmSJEla5jxJT5IkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSqsqruApZiY2Lng+PDwMCMjo12uRpIkSYOspwPyvqcbrFh5EFu2vH7B20dHx9ix4x5DsiRJktqmpwPy/r0zzP58PyecdiVDa8eecdvMk5Pcvu0cGo2GAVmSJElt09MBec7Q2jHWrNtQdxmSJElaBjxJT5IkSSp0ZAU5IlYBFwCbMvPdnXgOSZIkqRM6tYK8FvhyBx9fkiRJ6oglryBHxEbgUuDIzDy2GN8CnAJMALOZ+dHMnIqI3W2vVpIkSeqwA9licTxwM3DU3EBEDANXA6/OzJmIuCEiTs7Mr7W5zkUdcsg6xsfXd+vp+p696jx73B32uTvsc+fZ4+6wz90xKH1eckDOzOsj4qR5w5uBBzNzpjq+A9gKdC0g7979BKtX7+nW0/W18fH17NplrzrJHneHfe4O+9x59rg77HN39GOfFwv0re4RPhQoOzENHBoRK4BTgYiI32rxOSRJkqSuafUqFhNAGb1HgInMnAUur/6TJEmS+karK8jbgU0RMVQdHwfc0uJjSpIkSbVZckCOiBOB04EXRMSFEXFwZjaA9wKfiYhLgXu6eYKeJEmS1G4HcpLebcBtC4zfCtzazqIkSZKkuvhFHpIkSVLBgCxJkiQVDMiSJElSwYAsSZIkFQzIkiRJUsGALEmSJBUMyJIkSVLBgCxJkiQVDMiSJElSwYAsSZIkFZb8VdP9Znp6ikajseBtw8PDjIyMdrkiSZIk9YOBDMjT01McffQRTE1NLnj76OgYO3bcY0iWJEnSLxnIgNxoNJiamuSE065kaO3YM26beXKS27edQ6PRMCBLkiTplwxkQJ4ztHaMNes21F2GJEmS+ogn6UmSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSYVVdRfQa6anp2g0GgveNjw8zMjIaJcrkiRJUjcZkAvT01McffQRTE1NLnj76OgYO3bcY0iWJEkaYAbkQqPRYGpqkhNOu5KhtWPPuG3myUlu33YOjUbDgCxJkjTADMgLGFo7xpp1G+ouQ5IkSTXwJD1JkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCm3/opCIGAYuBh4CdmbmF9v9HJIkSVKnLCkgR8RG4FLgyMw8thjfApwCTACzmfnR6viuzPxiRNwEGJAlSZLUN5a6gnw8cDNw1NxAtVJ8NfDqzJyJiBsi4mTgRcD26m4Ht7PYQTQ9PUWj0VjwtuHhYUZGRrtcUXf8qte9f/8+Djpo4bfmIPdEkiT1hiUF5My8PiJOmje8GXgwM2eq4zuArcDdwHg19lQ7ivxVDjlkHePj658xtnfvnq7+3HM1NTXFMcceyeOTjy14+/PHNvCjB+5ndLR9gbBdtbfi2V73yoNW8fP9+xa8rRM9abde6PFyYJ+7wz53nj3uDvvcHYPS51b2IB8KlIlyuhq7Ebg4Ig4DtrXw+Euye/cTrF6955fGuvlzz9Wjj+7k8cnHOOG0KxlaO/aM22aenOT2befw0EM72bixPedSjo+vZ9eu9tTeil/1uqd3/Yg7/+GSrvWk3Xqlx4POPneHfe48e9wd9rk7+rHPiwX6VgLyBFA+6ggwkZkN4PwWHnfZGVo7xpp1G+ouo+sWet0zT04uepskSVI3tLIMtx3YFBFD1fFxwC2tlyRJkiTVZ8Xs7Oyz3ikiTgTeAbwZuAq4MjOfiog3AW8DdgF7q6tYSJIkSX1rSQFZkiRJWi5680wnSZIkqSYGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqdDKN+m1LCK2AKfQ/Fa+Wa+j3B4RcThwKXA38EJgd2ZeEhEbgMuA+4GXAx/OzJ31Vdr/IuJg4JvAVzPz3IhYA3wKeJhmjy/LzO/XWWO/i4gA/gh4CjgRuJjmnPER4AfAi4FzMvPZvytei4qI82j28qc037t/ChyMc0ZLImIjzfn4yMw8thpbdJ6IiD8BfhPYD/wwM6+ppfA+s0ifLwA2Ao8AxwAXZeZ91W32+QAt1OPittOAzwHr5+bifs94ta0gR8QwcDVwVmZeDBwRESfXVc+A2QB8PjOvyMz3A2+PiKOBTwD/lJmXATfRnKDVmkuBbxfHHwAeysy/BP4KuLaWqgZERBwEfBq4JDMvpxnaHqA5d1xT9fl7wAX1Vdn/ql98HwLOzMy/ANbS/MXmnNG644GbgRXF2ILzRES8EDgXODczzwfeFREv73K9/WqhPq8Dzs7MTwI3AFeAfW7BQj0mIl4JvGreWN9nvDq3WGwGHszMmer4DmBrjfUMjMy8KzNvLoZWAk/S7O/2asx+tygiTqfZxweK4V/0ODO/CxwZESM1lDcojqU5GZ8ZER8Cfh94HHgDcFd1H9/LrWsATwNz79V1wL04Z7QsM68H9swbXmye+F1gR2bOfYPXduD3ulVrP1uoz5n5kaKXK4G5T5ns83OwUI+rIHw+MH91uO8zXp0B+VCe2ejpakxtFBFvBb5SfaxU9nwaGIuIWrfZ9KuIeBXwysy8cd5Nvq/baxPNifa6arXt9TRXfp4qfrnZ4xZl5jRwHvCFiLgO+Fea21ecMzpjsXnC+aMDIuJ5wDuBC6sh+9w+H6f5Cd/T88b7vsd1BuQJYH1xPFKNqU0i4g00V9rOqobKno8Ak5m5r47aBsBbgZ9FxAdpfuz0uoj4AL6v220auC8zp6rjbwCvAQ6OiLmP+exxiyLiKJoBeWtmnkFzH/JFOGd0ymLzhPNHm1Xh+CrgzzPzh9WwfW6DiHgRMAacWv0uBDg7Io5hAHpcZ0DeDmyKiKHq+DjglhrrGSgRsZXmx0jvBzZGxGaa/d1c3cV+tyAzP56Zl1R7M78B3JmZf03R44h4LfCdanVOz803gUOqvcjQXFG+F/g6ze0X4Hu5HX4NeKwIv48Aa3DO6JTF5omvAEcX//jbDPxjPSX2v+rj/2uAT2fmjoj4w+om+9wGmfnjzDwjMy+rfhdCs9ffYgAy3orZ2dlnv1eHRMSbgLcBu4C9/XaGY6+qTsi7DfhWNbQW+CzwJeBy4EHgcOCDnpHemmrC/TPgeTR7PHci0yPAy4BPeBWL1lTbhN5Ic574deBM4DCaK5z3V2NnexWL5676B8hngJ/R3OP9Gponks3gnNGSiDgReAfwZpormVdWNy04T1RXVziG5tUVvu/VFZZmkT5vo/le/kl1t7XFFS7s8wFaqMeZ+VREjAPvAT5W/XdNZj7c7xmv1oAsSZIk9Rq/KESSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmS9Kwi4s0RsTMiTqq7FknqtFV1FyBJ6m0R8R+BVwD76q5FkrrBFWRJ0rO5KTP/I7C37kIkqRtcQZakNoiIw4H/DAwDzwO+ClwBXAVsBNYA3wPOzMx9EXER8D7geuD5wOuAfwBuAD4MHAGcnZk3RcRbgE8CP60e47XAEPDuzLy7eP6rqvGVwAWZ+S/Fz+4E7gROAA7LzJcs9bVl5r8+175IUj8yIEtSiyLiIOB/AZdn5t9FxChwL/CfgK9k5ueq+10HvBO4NjMviYiXAr8N/FtgFHgY2JOZfxARpwCX0Vy9/VJEbAD+C/DezPxuRPwxcGNEvByYBf43cEVm/teIOAL4PxHxkuJn/xZ4T2aeFxFXVPV8nmZ4X8iXM/OydvdKkvqBAVmSWvfbwMuAbQCZORURpwKTwKaI+AbN/bsvBhrAtcXP3paZTwO7ImIC+E41fg/w0nnPc29mfrf68xeA64DNwH7gcOB/VM9/T0Q8DPw74H9W98/MvK/6w3nV/9/e8iuXpAHkHmRJat0LgcnM/MVJbJl5B83V4vcAb8nMk2gG2uF5P7un+PO+4ngfsHrefSeLx98PPA68oHr+WeDWiPjniPhnmlstRoufnXoOr0uSliVXkCWpdT8GxiJi1VxIjohXAicCd2bmY9X95gfeA7Vh7g8RsYrm3uVHaIbpvVUIn7t9LfDzX/VgbrGQpIUZkCWpdd8EfgD8MfDfqz2/fw98DnhXRAzR3AZxMnBfC8/zGxHx2mqbxduBnwDbaQbhhyLilMy8sQrPNwHn8n+3bPwSt1hI0sJWzM7O1l2DJPW9eVexWAlcCHyLZkh+Bc2rT6yiuWf4k9WPnQ/8DDgLeBPwH2gG6FOBvwP+DXBrZv5ORJwB/Hvg28AxwME0r2Kxo3j+z1bjK4H/Vp2w90aaJ+htpLma/TvP4bVtBc6judf6PuBfMvN9B/o4ktQvDMiS1AeqgHxGuY1CktQZnqQnSZIkFQzIktTjqi/7+CBwVET8Td31SNKgc4uFJEmSVHAFWZIkSSr0xGXe9u3bPzs52ai7jIE3NjaMfe4se9wd9rk77HPn2ePusM/d0Y99Hh9fv2Kh8Z5YQV616qC6S1gW7HPn2ePusM/dYZ87zx53h33ujkHqc08EZEmSJKlX9MQWi6mpKR59dGetNQwPDzMyMlprDZIkSapfTwTkF7/kpTw++VitNYyOjrFjxz2GZEmSpGWuJwLy45OPccJpVzK0dqyW5595cpLbt51Do9EwIEuSJC1zbQ/I1QXtXwGsBr6fmV9cys8NrR1jzboN7S5HkiRJOiBLCsgRsRG4FDgyM48txrcApwATwGxmfhTYkZlfiohR4FpgSQFZkiRJ6gVLvYrF8cDNwC+uFRcRw8DVwFmZeTFwREScnJkPV3d5K/CpNtYqSZIkddySVpAz8/qIOGne8GbgwcycqY7vALYCX4uIrcD9wMP0kUMOWcf4+Pq6y+ioQX99vcAed4d97g773Hn2uDvsc3cMSp9b2YN8KLCnOJ4GDo2IPwAuAL4DrAdOa+E5umr37idYvXrPs9+xT42Pr2fXrsF9fb3AHneHfe4O+9x59rg77HN39GOfFwv0rQTkCZoBeM4IMJGZNwE3tfC4kiRJUm1a+Sa97cCmiBiqjo8Dbmm9JEmSJKk+SwrIEXEicDrwgoi4MCIOzswG8F7gMxFxKXBPZn6tg7VKkiRJHbfUk/RuA25bYPxW4NZ2FyVJkiTVpZUtFpIkSdLAMSBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVFjViQeNiFXABcCmzHx3J55DkiRJ6oROrSCvBb7cwceXJEmSOmLJK8gRsRG4FDgyM48txrcApwATwGxmfjQzpyJid9urlSRJkjrsQLZYHA/cDBw1NxARw8DVwKszcyYiboiIkzPza22usysOOWQd4+Pr6y6jowb99fUCe9wd9rk77HPn2ePusM/dMSh9XnJAzszrI+KkecObgQczc6Y6vgPYCvRlQN69+wlWr95TdxkdMz6+nl27Bvf19QJ73B32uTvsc+fZ4+6wz93Rj31eLNC3ukf4UKDsxDRwaESsAE4FIiJ+q8XnkCRJkrqm1atYTABl9B4BJjJzFri8+k+SJEnqG62uIG8HNkXEUHV8HHBLi48pSZIk1WbJATkiTgROB14QERdGxMGZ2QDeC3wmIi4F7unXE/QkSZIkOLCT9G4Dbltg/Fbg1nYWJUmSJNXFL/KQJEmSCh35qul+NTGxs+4SGB4eZmRktO4yJEmSli0DMrDv6QYrVh7Eli2vr7sURkfH2LHjHkOyJElSTQzIwP69M8z+fD8nnHYlQ2vHaqtj5slJbt92Do1Gw4AsSZJUEwNyYWjtGGvWbai7DEmSJNXIgNyDOrUXeu/ePeze/cSz3s990JIkaTkzIPeQXtkL7T5oSZK0nBmQe0gv7IV2H7QkSVruDMg9yL3QkiRJ9fGLQiRJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIJXsVDPmp6eotFo1FqDX5oiSdLyY0BWT5qenuLoo49gamqy1jr80hRJkpYfA7J6UqPRYGpq0i9NkSRJXWdAVk/zS1MkSVK3eZKeJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEkFA7IkSZJUMCBLkiRJBQOyJEmSVDAgS5IkSYVV7X7AiBgGLgYeAnZm5hfb/RySJElSpywpIEfERuBS4MjMPLYY3wKcAkwAs5n50er4rsz8YkTcBBiQJUmS1DeWuoJ8PHAzcNTcQLVSfDXw6syciYgbIuJk4EXA9upuB7ezWKkOExM7l3zfvXv3sHv3E219/uHhYUZGRtv6mP1qenqKRqPRkT4fCP9OJOnAzM3fdTqQuXtJATkzr4+Ik+YNbwYezMyZ6vgOYCtwNzBejT21pCrUcw45ZB3j4+tre/69e/fU9txz9j3dYMXKg9iy5fW11vH8sQ386IH7GR1d3oFsamqKY449kscnH6u7lGX1d1LnPLBc2OPusM/dsVCfe2X+PpC5u5U9yIcCZYqZrsZuBC6OiMOAbS08vmq0e/cTrF5dX0itc3Vwzv69M8z+fD8nnHYlQ2vHaqlh5slJbt92Dg89tJONG5f3ObWPPrqTxycfq/XvA5bX38n4+Hp27ar/H6uDzB53h33ujsX63Avz92Jz92L/cGolIE8A5aOOABOZ2QDOb+FxpZ4ytHaMNes21F2GKv59SFJ/6qf5u5Xlj+3ApogYqo6PA25pvSRJkiSpPitmZ2ef9U4RcSLwDuDNwFXAlZn5VES8CXgbsAvYW13FQpIkSepbSwrIkiRJ0nIx2GeYSJIkSQfIgCxJkiQVDMiSJElSwYAsSZIkFQzIkiRJUsGALEmSJBVa+Sa9lkXEFuAUmt/KN+t1lNsjIg4HLgXuBl4I7M7MSyJiA3AZcD/wcuDDmbmzvkr7X0QcDHwT+GpmnhsRa4BPAQ/T7PFlmfn9OmvsdxERwB8BTwEnAhfTnDM+AvwAeDFwTmbW//3kfSwizqPZy5/SfO/+KXAwzhktiYiNNOfjIzPz2Gps0XkiIv4E+E1gP/DDzLymlsL7zCJ9vgDYCDwCHANclJn3VbfZ5wO0UI+L204DPgesn5uL+z3j1baCHBHDwNXAWZl5MXBERJxcVz0DZgPw+cy8IjPfD7w9Io4GPgH8U2ZeBtxEc4JWay4Fvl0cfwB4KDP/Evgr4NpaqhoQEXEQ8Gngksy8nGZoe4Dm3HFN1efvARfUV2X/q37xfQg4MzP/AlhL8xebc0brjgduBlYUYwvOExHxQuBc4NzMPB94V0S8vMv19quF+rwOODszPwncAFwB9rkFC/WYiHgl8Kp5Y32f8ercYrEZeDAzZ6rjO4CtNdYzMDLzrsy8uRhaCTxJs7/bqzH73aKIOJ1mHx8ohn/R48z8LnBkRIzUUN6gOJbmZHxmRHwI+H3gceANwF3VfXwvt64BPA3MvVfXAffinNGyzLwe2DNveLF54neBHZk59w1e24Hf61at/WyhPmfmR4pergTmPmWyz8/BQj2ugvD5wPzV4b7PeHUG5EN5ZqOnqzG1UUS8FfhK9bFS2fNpYCwiat1m068i4lXAKzPzxnk3+b5ur000J9rrqtW219Nc+Xmq+OVmj1uUmdPAecAXIuI64F9pbl9xzuiMxeYJ548OiIjnAe8ELqyG7HP7fJzmJ3xPzxvv+x7XGZAngPXF8Ug1pjaJiDfQXGk7qxoqez4CTGbmvjpqGwBvBX4WER+k+bHT6yLiA/i+brdp4L7MnKqOvwG8Bjg4IuY+5rPHLYqIo2gG5K2ZeQbNfcgX4ZzRKYvNE84fbVaF46uAP8/MH1bD9rkNIuJFwBhwavW7EODsiDiGAehxnQF5O7ApIoaq4+OAW2qsZ6BExFaaHyO9H9gYEZtp9ndzdRf73YLM/HhmXlLtzfwGcGdm/jVFjyPitcB3qtU5PTffBA6p9iJDc0X5XuDrNLdfgO/ldvg14LEi/D4CrME5o1MWmye+Ahxd/ONvM/CP9ZTY/6qP/68BPp2ZOyLiD6ub7HMbZOaPM/OMzLys+l0IzV5/iwHIeCtmZ2ef/V4dEhFvAt4G7AL29tsZjr2qOiHvNuBb1dBa4LPAl4DLgQeBw4EPekZ6a6oJ98+A59Hs8dyJTI8ALwM+4VUsWlNtE3ojzXni14EzgcNornDeX42d7VUsnrvqHyCfAX5Gc4/3a2ieSDaDc0ZLIuJE4B3Am2muZF5Z3bTgPFFdXeEYmldX+L5XV1iaRfq8jeZ7+SfV3dYWV7iwzwdooR5n5lMRMQ68B/hY9d81mflwv2e8WgOyJEmS1Gv8ohBJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJ6lMR8e6I+FFEXFcdXxQRj0bExfVWJkn9zYAsSX0qM/9f4Lri+BLgy7UVJEkDwoAsSZIkFVbVXYAk6dlFxEXASdVhA3h3Zv6kvookaXC5gixJ/WESODkz3whcD1xecz2SNLBcQZak/vBj4OsRsRIYAZ5Xcz2SNLBcQZakHhcRLwf+HjgvM18PfAAYrrcqSRpcBmRJ6n2/CUxn5l3V8eo6i5GkQWdAlqTe9wNgLCJ+ozp+c53FSNKgWzE7O1t3DZKkZxERHwNOB74DPAq8E3gMeBpYA/wNcBDwPuBnwMcy89p6qpWk/mZAliRJkgpusZAkSZIKBmRJkiSpYECWJEmSCgZkSZIkqdAT36S3b9/+2YKNlkcAAAwzSURBVMnJRt1lDLyxsWHsc2fZ4+6wz91hnzvPHneHfe6Ofuzz+Pj6FQuN98QK8qpVB9VdwrJgnzvPHneHfe4O+9x59rg77HN3DFKfeyIgS5IkSb3CgCxJkiQVemIP8tTUFI8+urPWGoaHhxkZGa21BkmSJNWvJwLyi1/yUh6ffKzWGkZHx9ix4x5DsiRJ0jLXEwH58cnHOOG0KxlaO1bL8888Ocnt286h0WgYkCVJkpa5ngjIAENrx1izbkPdZUiSJGmZ8yQ9SZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqdD26yBHxFuAVwCrge9n5hfb/RySJElSpywpIEfERuBS4MjMPLYY3wKcAkwAs5n5UWBHZn4pIkaBawEDsiRJkvrGUrdYHA/cDKyYG4iIYeBq4KzMvBg4IiJOzsyHq7u8FfhUG2uVJEmSOm5JK8iZeX1EnDRveDPwYGbOVMd3AFuBr0XEVuB+4GH6yCGHrGN8fH3dZXTUoL++XmCPu8M+d4d97jx73B32uTsGpc+t7EE+FNhTHE8Dh0bEHwAXAN8B1gOntfAcXbV79xOsXr3n2e/Yp8bH17Nr1+C+vl5gj7vDPneHfe48e9wd9rk7+rHPiwX6VgLyBM0APGcEmMjMm4CbWnhcSZIkqTatXOZtO7ApIoaq4+OAW1ovSZIkSarPkgJyRJwInA68ICIujIiDM7MBvBf4TERcCtyTmV/rYK2SJElSxy31JL3bgNsWGL8VuLXdRUmSJEl18Zv0JEmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqrOrEg0bEKuACYFNmvrsTzyFJkiR1QqdWkNcCX+7g40uSJEkdseQV5IjYCFwKHJmZxxbjW4BTgAlgNjM/mplTEbG77dVKkiRJHXYgK7zHAzcDK+YGImIYuBo4KzMvBo6IiJPbWqEkSZLURUteQc7M6yPipHnDm4EHM3OmOr4D2Ap8rT3lddchh6xjfHx93WV01KC/vl5gj7vDPneHfe48e9wd9rk7BqXPrZ6kdyiwpzieBg6NiBXAqUBExG9l5t0tPk9X7N79BKtX73n2O/ap8fH17No1uK+vF9jj7rDP3WGfO88ed4d97o5+7PNigb7VgDwBlI88Akxk5ixwefVf35iY2Fl3CQwPDzMyMlp3GZIkSctWqwF5O7ApIoaqbRbHAX/belndte/pBitWHsSWLa+vuxRGR8fYseMeQ7IkSVJNDuQqFicCpwMviIgLgSszsxER7wU+ExG7gHsys+/2H+/fO8Psz/dzwmlXMrR2rLY6Zp6c5PZt59BoNAzIkiRJNTmQk/RuA25bYPxW4NZ2FlWXobVjrFm3oe4yJEmSVCO/yEOSJEkqdOSrptWaTp0suHfvHnbvfuJZ7+eJgpIkaTkzIPeQXjlZ0BMFJUnScmZA7iG9cLKgJwpKkqTlzoDcgzxZUJIkqT6epCdJkiQVDMiSJElSwYAsSZIkFdyDrJ41PT1Fo9GotQYveSdJ0vJjQFZPmp6e4uijj2BqarLWOrzknSRJy48BWT2p0WgwNTXpJe8kSVLXGZDV07zknSRJ6jZP0pMkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgoGZEmSJKlgQJYkSZIKBmRJkiSpYECWJEmSCgZkSZIkqWBAliRJkgqr2v2AETEMXAw8BOzMzC+2+zkkSZKkTllSQI6IjcClwJGZeWwxvgU4BZgAZjPzo9XxXZn5xYi4CTAg96GJiZ3L+vlLB1LL3r172L37ibY+//DwMCMjo219zH41PT1Fo9HoSJ8PhH8nknRg5ubvOh3I3L3UFeTjgZuBo+YGqpXiq4FXZ+ZMRNwQEScDLwK2V3c7eMlVqyfse7rBipUHsWXL6+supXa90ovnj23gRw/cz+jo8g5kU1NTHHPskTw++VjdpSyrv5Px8fV1lzDw7HF32OfuWKjPvTJ/H8jcvaSAnJnXR8RJ84Y3Aw9m5kx1fAewFbgbGK/GnlpSxeoZ+/fOMPvz/Zxw2pUMrR2rrY7pXT/izn+4pLbnh97oxcyTk9y+7RweemgnGzcu71MGHn10J49PPlb7e3M5/Z2Mj69n1649dZcx0Oxxd9jn7lisz70wfy82dy/2D6dW9iAfCpRdmK7GbgQujojDgG0tPL5qNLR2jDXrNtT2/DNPTtb23PPV3Qs9k38fktSf+mn+biUgTwBl7B4BJjKzAZzfUlWSJElSTVr5fHA7sCkihqrj44BbWi9JkiRJqs+K2dnZZ71TRJwIvAN4M3AVcGVmPhURbwLeBuwC9lZXsZAkSZL61pICsiRJkrRcDPYp2JIkSdIBMiBLkiRJBQOyJEmSVDAgS5IkSQUDsiRJklQwIEuSJEmFVr5Jr2URsQU4hea38s16HeX2iIjDgUuBu4EXArsz85KI2ABcBtwPvBz4cGburK/S/hcRBwPfBL6amedGxBrgU8DDNHt8WWZ+v84a+11EBPBHwFPAicDFNOeMjwA/AF4MnJOZT9RU4kCIiPNo9vKnNN+7fwocjHNGSyJiI835+MjMPLYaW3SeiIg/AX4T2A/8MDOvqaXwPrNIny8ANgKPAMcAF2XmfdVt9vkALdTj4rbTgM8B6+fm4n7PeLWtIEfEMHA1cFZmXgwcEREn11XPgNkAfD4zr8jM9wNvj4ijgU8A/5SZlwE30Zyg1ZpLgW8Xxx8AHsrMvwT+Cri2lqoGREQcBHwauCQzL6cZ2h6gOXdcU/X5e8AF9VXZ/6pffB8CzszMvwDW0vzF5pzRuuOBm4EVxdiC80REvBA4Fzg3M88H3hURL+9yvf1qoT6vA87OzE8CNwBXgH1uwUI9JiJeCbxq3ljfZ7w6t1hsBh7MzJnq+A5ga431DIzMvCszby6GVgJP0uzv9mrMfrcoIk6n2ccHiuFf9DgzvwscGREjNZQ3KI6lORmfGREfAn4feBx4A3BXdR/fy61rAE8Dc+/VdcC9OGe0LDOvB/bMG15snvhdYEdmzn2D13bg97pVaz9bqM+Z+ZGilyuBuU+Z7PNzsFCPqyB8PjB/dbjvM16dAflQntno6WpMbRQRbwW+Un2sVPZ8GhiLiFq32fSriHgV8MrMvHHeTb6v22sTzYn2umq17fU0V36eKn652eMWZeY0cB7whYi4DvhXmttXnDM6Y7F5wvmjAyLiecA7gQurIfvcPh+n+Qnf0/PG+77HdQbkCWB9cTxSjalNIuINNFfazqqGyp6PAJOZua+O2gbAW4GfRcQHaX7s9LqI+AC+r9ttGrgvM6eq428ArwEOjoi5j/nscYsi4iiaAXlrZp5Bcx/yRThndMpi84TzR5tV4fgq4M8z84fVsH1ug4h4ETAGnFr9LgQ4OyKOYQB6XGdA3g5sioih6vg44JYa6xkoEbGV5sdI7wc2RsRmmv3dXN3FfrcgMz+emZdUezO/AdyZmX9N0eOIeC3wnWp1Ts/NN4FDqr3I0FxRvhf4Os3tF+B7uR1+DXisCL+PAGtwzuiUxeaJrwBHF//42wz8Yz0l9r/q4/9rgE9n5o6I+MPqJvvcBpn548w8IzMvq34XQrPX32IAMt6K2dnZZ79Xh0TEm4C3AbuAvf12hmOvqk7Iuw34VjW0Fvgs8CXgcuBB4HDgg56R3ppqwv0z4Hk0ezx3ItMjwMuAT3gVi9ZU24TeSHOe+HXgTOAwmiuc91djZ3sVi+eu+gfIZ4Cf0dzj/RqaJ5LN4JzRkog4EXgH8GaaK5lXVjctOE9UV1c4hubVFb7v1RWWZpE+b6P5Xv5Jdbe1xRUu7PMBWqjHmflURIwD7wE+Vv13TWY+3O8Zr9aALEmSJPUavyhEkiRJKhiQJUmSpIIBWZIkSSoYkCVJkqSCAVmSJEkqGJAlSZKkggFZkiRJKvz/jwxfEZwk1I8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(3,1,figsize=(10, 6))\n",
    "\n",
    "ax[0].hist(df[df.camper==0]['count'], bins=20, edgecolor=\"k\", linewidth=1.2)\n",
    "ax[0].set_yscale('log')\n",
    "ax[0].set_xlim([0,150])\n",
    "ax[0].set_title('camper=0')\n",
    "ax[1].hist(df[df.camper==1]['count'], bins=20, edgecolor=\"k\", linewidth=1.2)\n",
    "ax[1].set_yscale('log')\n",
    "ax[1].set_xlim([0,150])\n",
    "ax[1].set_title('camper=1')\n",
    "ax[2].hist(df['count'], bins=20, edgecolor=\"k\", linewidth=1.2)\n",
    "ax[2].set_yscale('log')\n",
    "ax[2].set_xlim([0,150])\n",
    "ax[2].set_title('all')\n",
    "\n",
    "f.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the $p=2$ is for the choice between **NB1** model (variance = mean + $\\alpha\\times$mean) and the **NB2** model (variance = mean + $\\alpha\\times$mean$^2$). (see here: https://towardsdatascience.com/negative-binomial-regression-f99031bb25b4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Maximum number of iterations has been exceeded.\n",
      "         Current function value: 1.731564\n",
      "         Iterations: 35\n",
      "         Function evaluations: 36\n",
      "         Gradient evaluations: 36\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/base/model.py:512: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "endog_var = df['count'].values\n",
    "exog_vars = np.concatenate([np.ones(df.shape[0]).reshape(-1,1),\n",
    "                            df.child.values.reshape(-1,1),\n",
    "                            df.camper.values.reshape(-1,1)],\n",
    "                           axis=1)\n",
    "exog_vars_infl = np.concatenate([np.ones(df.shape[0]).reshape(-1,1),\n",
    "                                 df.persons.values.reshape(-1,1)],\n",
    "                                axis=1)\n",
    "m = ZINB(endog = endog_var,\n",
    "         exog = exog_vars,\n",
    "         exog_infl=exog_vars_infl,\n",
    "         inflation='logit',\n",
    "         p=2).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/base/model.py:492: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available\n",
      "  'available', HessianInversionWarning)\n",
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/base/model.py:492: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available\n",
      "  'available', HessianInversionWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ZeroInflatedNegativeBinomialP Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>                 <td>y</td>               <th>  No. Observations:  </th>  <td>   250</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>           <td>ZeroInflatedNegativeBinomialP</td> <th>  Df Residuals:      </th>  <td>   247</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                       <td>MLE</td>              <th>  Df Model:          </th>  <td>     2</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                  <td>Fri, 17 Apr 2020</td>        <th>  Pseudo R-squ.:     </th>  <td>0.06793</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                      <td>00:38:47</td>            <th>  Log-Likelihood:    </th> <td> -432.89</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>                   <td>False</td>             <th>  LL-Null:           </th> <td> -464.44</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>           <td>nonrobust</td>           <th>  LLR p-value:       </th> <td>1.989e-14</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>           <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_const</th> <td>    1.6027</td> <td>    0.836</td> <td>    1.917</td> <td> 0.055</td> <td>   -0.036</td> <td>    3.242</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_x1</th>    <td>   -1.6661</td> <td>    0.679</td> <td>   -2.454</td> <td> 0.014</td> <td>   -2.997</td> <td>   -0.336</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>         <td>    1.3709</td> <td>    0.256</td> <td>    5.353</td> <td> 0.000</td> <td>    0.869</td> <td>    1.873</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>            <td>   -1.5152</td> <td>    0.196</td> <td>   -7.747</td> <td> 0.000</td> <td>   -1.898</td> <td>   -1.132</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>            <td>    0.8791</td> <td>    0.269</td> <td>    3.265</td> <td> 0.001</td> <td>    0.351</td> <td>    1.407</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>alpha</th>         <td>    2.6784</td> <td>    0.471</td> <td>    5.684</td> <td> 0.000</td> <td>    1.755</td> <td>    3.602</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                     ZeroInflatedNegativeBinomialP Regression Results                    \n",
       "=========================================================================================\n",
       "Dep. Variable:                                 y   No. Observations:                  250\n",
       "Model:             ZeroInflatedNegativeBinomialP   Df Residuals:                      247\n",
       "Method:                                      MLE   Df Model:                            2\n",
       "Date:                           Fri, 17 Apr 2020   Pseudo R-squ.:                 0.06793\n",
       "Time:                                   00:38:47   Log-Likelihood:                -432.89\n",
       "converged:                                 False   LL-Null:                       -464.44\n",
       "Covariance Type:                       nonrobust   LLR p-value:                 1.989e-14\n",
       "=================================================================================\n",
       "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "inflate_const     1.6027      0.836      1.917      0.055      -0.036       3.242\n",
       "inflate_x1       -1.6661      0.679     -2.454      0.014      -2.997      -0.336\n",
       "const             1.3709      0.256      5.353      0.000       0.869       1.873\n",
       "x1               -1.5152      0.196     -7.747      0.000      -1.898      -1.132\n",
       "x2                0.8791      0.269      3.265      0.001       0.351       1.407\n",
       "alpha             2.6784      0.471      5.684      0.000       1.755       3.602\n",
       "=================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative: run the model with robust standard errors -- these are larger for the parameters of the count component of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: Maximum number of iterations has been exceeded.\n",
      "         Current function value: 1.731564\n",
      "         Iterations: 35\n",
      "         Function evaluations: 36\n",
      "         Gradient evaluations: 36\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/base/model.py:512: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "m = ZINB(endog = endog_var,\n",
    "         exog = exog_vars,\n",
    "         exog_infl=exog_vars_infl,\n",
    "         inflation='logit',\n",
    "         p=2).fit(cov_type='HC0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/base/model.py:492: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available\n",
      "  'available', HessianInversionWarning)\n",
      "/home/seb/py36/lib/python3.6/site-packages/statsmodels/base/model.py:492: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available\n",
      "  'available', HessianInversionWarning)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ZeroInflatedNegativeBinomialP Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>                 <td>y</td>               <th>  No. Observations:  </th>  <td>   250</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>           <td>ZeroInflatedNegativeBinomialP</td> <th>  Df Residuals:      </th>  <td>   247</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                       <td>MLE</td>              <th>  Df Model:          </th>  <td>     2</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                  <td>Fri, 17 Apr 2020</td>        <th>  Pseudo R-squ.:     </th>  <td>0.06793</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                      <td>00:38:47</td>            <th>  Log-Likelihood:    </th> <td> -432.89</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>                   <td>False</td>             <th>  LL-Null:           </th> <td> -464.44</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>              <td>HC0</td>              <th>  LLR p-value:       </th> <td>1.989e-14</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>           <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_const</th> <td>    1.6027</td> <td>    0.665</td> <td>    2.410</td> <td> 0.016</td> <td>    0.299</td> <td>    2.906</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_x1</th>    <td>   -1.6661</td> <td>    0.430</td> <td>   -3.873</td> <td> 0.000</td> <td>   -2.509</td> <td>   -0.823</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>         <td>    1.3709</td> <td>    0.389</td> <td>    3.520</td> <td> 0.000</td> <td>    0.608</td> <td>    2.134</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>            <td>   -1.5152</td> <td>    0.241</td> <td>   -6.280</td> <td> 0.000</td> <td>   -1.988</td> <td>   -1.042</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>            <td>    0.8791</td> <td>    0.470</td> <td>    1.869</td> <td> 0.062</td> <td>   -0.043</td> <td>    1.801</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>alpha</th>         <td>    2.6784</td> <td>    0.577</td> <td>    4.645</td> <td> 0.000</td> <td>    1.548</td> <td>    3.809</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                     ZeroInflatedNegativeBinomialP Regression Results                    \n",
       "=========================================================================================\n",
       "Dep. Variable:                                 y   No. Observations:                  250\n",
       "Model:             ZeroInflatedNegativeBinomialP   Df Residuals:                      247\n",
       "Method:                                      MLE   Df Model:                            2\n",
       "Date:                           Fri, 17 Apr 2020   Pseudo R-squ.:                 0.06793\n",
       "Time:                                   00:38:47   Log-Likelihood:                -432.89\n",
       "converged:                                 False   LL-Null:                       -464.44\n",
       "Covariance Type:                             HC0   LLR p-value:                 1.989e-14\n",
       "=================================================================================\n",
       "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "inflate_const     1.6027      0.665      2.410      0.016       0.299       2.906\n",
       "inflate_x1       -1.6661      0.430     -3.873      0.000      -2.509      -0.823\n",
       "const             1.3709      0.389      3.520      0.000       0.608       2.134\n",
       "x1               -1.5152      0.241     -6.280      0.000      -1.988      -1.042\n",
       "x2                0.8791      0.470      1.869      0.062      -0.043       1.801\n",
       "alpha             2.6784      0.577      4.645      0.000       1.548       3.809\n",
       "=================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = m.predict(exog=exog_vars, exog_infl=exog_vars_infl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAFyCAYAAABbdsanAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3jcV33n8bcutmRFkpHtMQl2LhBHByghcRIoIYDBcdunTdMupM8TcLnkgdItBZYAAQoJEGjShOVaum1hS1mgKZA+sKWwlMsmhgC5cEnMdclxHAjETkCyrVhSZI0u1v4xI6PYM6OLNXPmN/N+PU+eeM4ZzXylMxp95vzO7/xaZmZmkCRJUjqtqQuQJElqdgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSqw9dQHHa3BwpGr7dvT1dTE0NFath9cyc7yyxfHKFscrWxyv+pXL9bSUaneGrIL29rbUJWgRHK9scbyyxfHKFscrewxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkqanlJ6cZGBojPzmdrIbMXzpJkiRpKaYPH+bGHbvZuWuQA8N51vR2sLk/x6VbN9HWWts5KwOZJElqSjfu2M1N39tz5Pb+4fyR29u39de0Fg9ZSpKkppOfnGbnrsGSfTt37av54UsDmSRJajoHR/McGM6X7BsaGefgaOm+ajGQSZKkprO6u4M1vR0l+/p6OlndXbqvWgxkkiSp6XSsaGNzf65k3+b+dXSsaKtpPS7qlyRJTenSrZuAwpqxoZFx+no62dy/7kh7LRnIJElSU2prbWX7tn4u2XI6B0fzrO7uqPnM2CwDmSRJamodK9pY39eVtAbXkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlVrWd+kMIJwLXAGfFGJ9yVN+fAjcAPTHG0WLbNuB5wAAwE2N8R7VqkyRJqifVnCF7BvAfQMvcxhDCE4AnHtXWBXwIeG2M8WrgySGEC6tYmyRJUt2oWiCLMX4GGJnbVgxebwSOnv06H/hFjDFfvH0rcFG1apMkSaontb64+LXAO2OMEyGEue3reWR4Gy62zauvr4v29updmT2X66naY2v5OV7Z4nhli+OVLY5XttQskIUQTgb6gEvnhLHXhRD+k8K6sbmvnN5i27yGhsaWs8xHyOV6GBwcmf+OqguOV7Y4XtnieGWL41W/ygXlmgWyGOP9wGWzt0MI1wHvizGOFg9lnhpC6CgetrwA+Ida1SZJkpRS1daQhRC2AC8CTgohXBVCWFVsz4UQrire7Y0hhA0xxjHgFcAHQwjXAD+MMd5crdokSZLqScvMzEzqGo7L4OBI1b4Bp3yzxfHKFscrWxyvbHG86lcu19NSqt2NYSVJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJtVfzwUMIJwLXAGfFGJ9SbHsTcCLwIHAe8LYY493FvhcCm4Fp4N4Y44erWZ8kSVI9qPYM2TOA/wBa5rR1A6+LMf534LPAuwFCCBuBK4ArYoxvBP4shHBGleuTJElKrqqBLMb4GWDkqLa3xhhn5jz/aPHfvwfcOafvduD3q1mfJElSPajqIctKQggrgZcAryw2reeR4W242FZRX18X7e1ty19gUS7XU7XH1vJzvLLF8coWxytbHK9sSRLIimHsH4ErY4z3FpsHgE1z7tYL7J7vsYaGxpa/wKJcrofBwZH576i64Hhli+OVLY5Xtjhe9atcUK75WZYhhC7gw8D7Yox3hhAuKXZ9BTg3hDC73ux84Eu1rk+SJKnWqn2W5RbgRcBJIYSrgPcCNwBPAh4bQgA4AfhsjHFPCOE9wPtDCNPAR2KM91SzPkmSpHpQ1UAWY7wFuOWo5udVuP8NFAKbJElS03BjWEmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJi7dV88BDCicA1wFkxxqcU2zqB9wB7gTOA62OMu4p9LwQ2A9PAvTHGD1ezPkmSpHpQ7RmyZwD/AbTMabsc+GWM8Trg/cA/A4QQNgJXAFfEGN8I/FkI4Ywq1ydJkpRcVQNZjPEzwMhRzRcBtxf7fwScFULoBX4PuDPGOFO83+3A71ezPkmSpHpQ1UOWZaznkSFtuNhWrr2ivr4u2tvblrXAuXK5nqo9tpaf45Utjle2OF7Z4nhlS4pANgDMfZX0FtsGgE1Hte+e78GGhsaWtbi5crkeBgePnuBTvXK8ssXxyhbHK1scr/pVLiinOMvyi8D5ACGEM4EfxBiHga8A54YQZtebnQ98KUF9kiRJNVXVQBZC2AK8CDgphHBVCGEV8LfAqSGEq4DXAy8DiDHuoXD25ftDCO8FPhJjvKea9UmSJNWDqh6yjDHeAtxSouuVZe5/A3BDNWuSJEmqN24MK0mSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUqsPcWThhDeAJwG7APOAF4GrAKuB35WbHtLjPHXKeqTJEmqpZrPkIUQTgTeDLw6xvh24ATgecDfADfFGK8HPge8p9a1SZIkpbCkQBZCePtxPOcYMAH0Fm93Az8BLgJuL7bdWrwtSZLU8FpmZmYq3iGE8OfA1cB6oKX430yMsW2pTxpCeBHwQuDB4uO9EtgPPDrG+FAIoR2YBFbEGKcqPdbU1PRMe/uSS5EkSaqlllKNC1lDdjmwBbg3xngYIITw5qVWEUI4G3gDcE6McSqE8F7gbcAA0AM8RGH2bGi+MAYwNDS21FLmlcv1MDg4UrXH1/JyvLLF8coWxytbHK/6lcv1lGxfSCD7QYzxnqPavnQctWwADswJWw8CpwBfBM4H7gcuKN6WJElqeGUDWQjhbcV/PhxCuBm4A8gX2/4AeNoSn/PLwB8UZ8YeAp5EYRYuD7wrhNAPnA5cscTHlyRJypRKM2QXA/+HwozV/cW2lqP+v2gxxmkKa8ZKeflSH1eSJCmrKgWyN8YYv1aqI4TwlSrVI0mS1HTKbnsxG8ZCCH85tz2EcBnwguqWJUmS1DwWsg/Z4+feiDF+DOiqSjWSJElNqNKi/p8DM8CaEMIfzulqA35c7cIkSZKaRaU1ZM+msHj/amDuzvzjXmNSkiRp+ZQNZDHGXxT/ednRfSGEJ8UYnSWTJElaBvNuDBtCeHGJ5hcCv7v85UiSJDWfhezU/1fAt4v/XgGcBdxXrYIkSZKazUIC2atijDtmb4QQWoC/rl5JkiRJzWXebS/mhrE5X3NmdcqRJElqPgtZQza7/QUUtrxYDfx9NYuSJElqJgs5ZHkLv9n2YgYYjDEeql5JkiRJzWUhgawP+J0Y40eqXYwkSVIzWsilk04FPlrtQiRJkprVQgLZ14F1cxtCCNdUpRpJkqQmtJBDlucAu0IIPwHyxbYzgKuqVpUkSVITWUgg6wL+eM7tjcDzqlOOJElS81lIIHsgxnjL7I0QwoUUduyXJEnSMigbyEIIpxT/uT+EcDLQUrz9c2Cy2oVJkiQ1i0ozZLOzYmuAZ89pPwR8qloFSZIkNZuygSzG+FiAEMJ/izF+sHYlSZIkNZeFXMvSMCZJkhpWfnKagaEx8pPTyWpYyKJ+SZKkhjN9+DA37tjNzl2DHBjOs6a3g839OS7duom21oVs1bp8DGSSJKkp3bhjNzd9b8+R2/uH80dub9/WX9Naahv/JEmS6kB+cpqduwZL9u3cta/mhy8NZJIkqekcHM1zYDhfsm9oZJyDo6X7qsVAJkmSms7q7g7W9HaU7Ovr6WR1d+m+ajGQSZKkptOxoo3N/bmSfZv719Gxoq2m9bioX5IkNaVLt24CCmvGhkbG6evpZHP/uiPttWQgkyRJTamttZXt2/q5ZMvpHBzNs7q7o+YzY7MMZJIkqal1rGhjfV9X0hpcQyZJkpRYshmyEEIAXkDhYuVbgKuBAeCtwG7gNOD1McbRRCVKkiTVRJIZshBCG/A+4J0xxncBLwN+DnwI+HCM8Trgx8CbUtQnSZJUS6kOWT4FaAFeHUJ4M3Ax8BDwHOC7xfvcClyUpjxJkqTaSXXI8lTgfOAFMcaDIYQbgLXAoRjjTPE+w8D6+R6or6+L9vbqnRGRy/VU7bG1/ByvbHG8ssXxyhbHK1tSBbJh4O4Y48Hi7W8BzwRWhRBaiqGsl8KasoqGhsaqVmQu18Pg4EjVHl/Ly/HKFscrWxyvbHG86le5oJzqkOW3gbXFtWRQmDH7CfA1CoczAS4AvpigNkmSpJpKMkMWYzwQQngT8IEQwiCQA94JfBJ4Wwjhd4FTgNelqE+SJKmWkm17EWP8d+Dfj2q+D3hp7auRJElKx41hJUmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMTaUz1xCGEV8G3gqzHGK0IIncB7gL3AGcD1McZdqeqTJEmqlZQzZNcAO+fcvhz4ZYzxOuD9wD8nqUqSJKnGkgSyEMKLgFuBn89pvgi4HSDG+CPgrBBCb4LyJEmSaqrmhyxDCE8EnhBjfEsI4clzutYDI3NuDxfbhis9Xl9fF+3tbctfaFEu11O1x9byc7yyxfHKFscrWxyvhRufmGJoOE9fbwedK9Os5krxrM8FxkMIfwU8A1gZQrgcGADmvnp6i20VDQ2NVaVIKLyYBwdH5r+j6oLjlS2OV7Y4XtnieC3M9OHD3LhjNzt3DXJgOM+a3g429+e4dOsm2lqrcxCxXFCueSCLMV47++/iQv7uGOMHiv8+H/hmCOFM4AcxxoqzY5IkSUt1447d3PS9PUdu7x/OH7m9fVt/TWtJtqg/hHAJ8CzgaSGEFwB/C5waQrgKeD3wslS1SZKkxpafnGbnrsGSfTt37SM/OV3TepJtexFj/Czw2aOaX5miFkmS1FwOjuY5MJwv2Tc0Ms7B0Tzr+7pqVo8bw0qSpKazuruDNb0dJfv6ejpZ3V26r1oMZJIkqel0rGhjc3+uZN/m/nV0rKjeDg6lJDtkKUmSlNKlWzcBhTVjQyPj9PV0srl/3ZH2WjKQSZKkptTW2sr2bf1csuV0Do7mWd3dUfOZsVkGMkmS1NQ6VrTVdAF/Ka4hkyRJSsxAJkmSlJiBTJIkNbWRsQl+et8BRsYmktXgGjJJktSUJqamuPYTd7F3cJTDM9DaAhty3Vz54nNY2V7biOQMWQPJT04zMDRW88s9SJKURdd+4i7uHyiEMYDDM3D/wCjXfuKumtfiDFkDSHG1ekmSsmxkbIK9g6Ml+/YOjjIyNkFP18qa1eNf6wYwe7X6/cN5ZvjN1epv3LE7dWmSJNXMYo4U7ZkzM3a0wzOF/lpyhizj5rta/SVbTk+2yZ0kSbWwlCNFG9d309pCyVDW2lLoryVnyDJuIVerl1JybaOkalvKkaKerpVsyJUOXRty3TU9XAnOkGXe7NXq95cIZSmuVi/Ncm2jpFo4niNFV774nLJnWdaagSzjZq9Wf9P39hzTl+Jq9dKs2U+ss2Y/sQJs39afqixJDWYhR4rKXRZpZXs773jpUxkZm2DPwCgb19d+ZmyWH1MbwKVbN7HtvI2s7e2ktQXW9nay7byNSa5WL8H8n1g9fClpucweKSploUeKerpW8oTT1iQLY+AMWUOop6vVS7CwT6wba1yTpMbUKEeKDGQNpB6uVi+Baxsl1dbsEaGdu/YxNDJOX08nm/vXZepIkYFM0rJrlE+skrKhEY4UGcgkVUUjfGKVlC1ZPlJkIJNUFY3wiVWSasVAJqmqsvyJVZJqxW0vJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkNbX85DQDQ2NJL+vmWZaSJKkpTR8+zI07drNz1yAHhvOs6e1gc3+OS7duoq21tnNWBjJJkhrUyNgEewZG2bi+O+mFs+vVjTt2P+KKIvuH80dub9/WX9NaDGSSJDWYiYkp3v7R77B3cJTDM9DaAhty3Vz54nNY2e6ffigcpty5a7Bk385d+7hky+k13czaNWSSJDWYN/zdN7l/oBDGAA7PwP0Do1z7ibvSFlZHDo7mOTCcL9k3NDLOwdHSfdWSJCaHEE4HrgHuAjYC+2OM7wwhrAGuB34GnAG8Jcb46xQ1SpKURSNjE9z3q+GSfXsHRxkZm/DwJbC6u4M1vR3sLxHK+no6Wd3dUdN6Us2QrQE+HWN8d4zxNcDzQwjnAn8D3BRjvB74HPCeRPVJkpRJewZGOXy4dN/hmUK/Cpd129yfK9m3uX9dza+9m2SGLMb43aOaWoGHgYuAa4tttwIfr2VdkiRl3cb13bS2UjKUtbYU+lVw6dZNQGHN2NDIOH09nWzuX3ekvZaSr+wLITwX+EqM8e4QwnpgpNg1DPSFENpjjFPlvr6vr4v29uql2Fyup2qPreXneGWL45Utjlc25IDTTuzlZw8ce9jytJN6edypa2tfVB17zQvOZXxiiqHhPH29HXSuTBONkgayEMJzgOcAlxebBoAe4CGgFxiqFMYAhobGqlZfLtfD4ODI/HdUXXC8ssXxyhbHK1ve/epncvn7bznmLMs3bj/bcSyjHRg5eIhq/3TKfbBJFshCCBcBzwReA5wUQjgV+CJwPnA/cEHxtiRJWoSVK9t5x0uf6j5kGZLqLMtzgRuB7wFfA04A/h54C/CuEEI/cDpwRYr6JElqBD1dK3nCaWtSl6EFSLWo/06g3KrCl9eyFjWm/OQ0B0fzrO7uqPmZMpIkLVbyRf3Scqqn65JJkrRQ/oVSQ5m9Ltn+4Twz/Oa6ZDfu2J26tKaQn5xmYGiM/OR06lIkNaEsvwc5Q6aGMT4xVVfXJWsmzkxKSqkR3oMMZGoYQ8PzX5dsfV9XjatqDrMzk7NmZyYBtm/rT1WWpCbRCO9B2YiN0gL09RauS1ayr6ej5tclaxb5yemKM5NZPHQgKTsa5T3IQKaG0bmyna7OFSX7ujpXeLiySg6Ozj8zKUnVshzvQfWw9sxDlmoY4xNTPHxoomTfw4cmyU9OG8qqYHV3YWZyf4k3xL6eTmcmJVXV8bwH1dPaM2fI1DCGhvMMjZQOZA+N5p2pqZKOFW1s7s+V7Nvcv84QLKmqjuc9qJ7OzDeQqWFUXkPmTE01Xbp1E9vO28ja3k5aW2BtbyfbztvIpVs3pS5NUhNYyntQva0985ClGkbnynY29+cecabNLGdqqquttZXt2/q5ZMvpXiFBUs0t5T1oIWvPanlmvoFMDWX209DOXfsYGhmnr6eTzf3rnKmpkY4VbW4tIimZxbwHre7uoK9nJQdKLHV5VHftz8w3kKmhOFMjSc1rMdcx7ljRxgmrSgeyE1bV/sx8A5kakjM1ktQ8lnK2ZH5ymrHxyZJ9Y+O1PzPfRf1Sk6jFPjv1sJePpOazlLMlK68hq/2Z+c6QSQ2uFvvs1NNePpKay3xnS5a7jnG97aHoO6UakjM1v1GLfXbqaS8fSc1lqTv119seis6QqaE4U/NIS/3kWG/PIUnlHM9M1588+3HEXz7E3sFRDs9AawtsyHXzJ89+XDVLLqn5/kKpoTlT80i1uM5kM17L0hlYqX4cz0zXZ77+M+4fKIQxgMMzcP/AKJ/5+s+qUWpFzpCpYYxPTFWYqRlsypmaWqyRqLd1GNXkDKxUn5Yy01Vvs/u+g6hhDA3nS4YCKMyUNeJMzXxqsUZivucAGmY2yRlYqT4tZaar3mb3nSFTw+jqbKe1hSO/kHO1tsCqjuZ8udfi6gWlnuPsM9ZyeGaGq/7pjoaYTaq3T9OSCvKT09x5969L9t1590DFsyzdqV/A4nYUbjZL+dmMjU+VDGNQCGmH8lP0dK1cxiqzoRZXL5h9jouffhp7BkbZuL6bL9x2HzfPua7o7GwSwPZt/cv6/LVQb9e9k1RwcDTP0GjpDV6HRifK/m66U79ch1LB8fxs+no76OteUfIXs697RUOtZVqKal694Ohx6+tZyeih0m+Qs59Ys6aZ1spJWdLW2rKkfnfql+tQKjjen01+8nDJ9snpGWchq+jocTswMsHEVOnpytlPrFlTb3sWSQvRDGcEDwwdWlL/wdH6WnfsDFmNuQ6lvPzkNHfFgZJ9d8XyZ0nOzs58f/c+xvKl33TGxqcYGZtoykOW1VbpNV3OfJ9o69Vzn/VY7vjJrxg9NHWkrXtVO8991mMTViUda3r6MJ+8aVdTHIl5VHfl9/Vy/as62mkBSn10bKH2644ba1QyoNpndWT509DB0XzJY/kABypcV2x2dmbfQ+NlH/vwDOwZGF2WOrOqWq+NSq/pcub7RFuvrr9h5yPCGMDooSmuv2Fnooqk0j76hZ80zZGYh0ZL/92Yr/9QfqpkGINCSDuUnyrTWx3OkNVYtdahNMK6tFUdiz9LcjGzMxvXdx9viZlU7ddGpdd0Ka0t2RyLkbGJsqF+z8CoM7CqG/nJae748YMl+xpxT8b53k/K9S/lb041ZeMvdQNZrnUoR892NMK6tEP5+c+SPFqlNQBztUDT/rGs9muj0mu6lA257kyOxZ6B0Yqfppt9Blb14+BovuwsdCPuydjTtZLuVaXDU/eq9rLvN0v5m1NNzpAlcDz7QpWa7Xjy6Wv54b37S94/S5+GVncv/izJSp9w5poB9h88xNrVq5ap2myovD/Pr5fttVFqbVVXZxv5iWmm55xn0d7Wwhu2n33cz5fCUtepSLVW6X0xxdqoastPTtNeZnqpvZWyZ0uu7u5gbZnZ/bW97kPWFNpaW7lky+k87YnrGT00xWNP6l3wjMHsbMes/cN5vrbzgbL3n/00lIX9kTpWtDE5XTpZlTtLstInnKPFXz7E089srkBWeX+eyWV7bZRaWzU2fuxatanpGd79ye/zjpc+9bifs9YWsk7lpHU1KkaqoNL74uzaqCzOUpdzcDTPQw+Xns166OGpivuQbe7PPeJv6qzN/Tn3IWt004cP86mb7+G2Hz3I+ERh6qBzZRsXnHkiz7/wjIpreirNdpSTpR3qR8YmGBsv/UtV7izJ1d0drO5q4+DY/AvVwymPWpY6s+RQvnQYW2j/QoyMTbB3cOGH6/YOFtZbLfwgZ31whkxZUYvf+3qy1H3IYGnXwKwW15DN48F9o3z5jl/w4L5Rdu8Z4uNf+im79wwd6f/cN+7hv33gFj73jXsW9Hg37tjNjjv3HgljAOMT09x859551/RUmu0oJ8Vx8KXaM+c6ZEcrd5Zkx4o2pmfm30Khva2l6Q5XAvz0voeOq38hKo1bKVk943WpZ3JJtVaL3/t6stR9yGBp18CslrqbOgkhbAOeBwwAMzHGd6SoY3R8gr940xeYmCoEp3/7+r1H+m75wbFnr3z+tvv5/G338/rnn81vnbam5GPmJ6e57YfHTo3OurPCXluw9L2bsjJD1tVZeXq4VP/I2MQxh8pK6VnVXvNdl+vBpo29x9W/EOv7Fhd0s3qWZbkdvRfaL9XKY9ZVXoYwX3/WbFzfXfFsyXLvN/W2L2hdzZCFELqADwGvjTFeDTw5hHBhilpe93e3Hglji/HeT3+/bN/B0TxjFT5ED1XYawuWvndTVmbI9g6OLbp/oTMts+ulms1kmd3yF9q/ENOLmR4ju2dZfvfu0psWL7RfqpUH9lV+L52vP2t6ulbSWmbCorW1pez7TbX3BV2sugpkwPnAL2KMsz+FW4GLal3Eg/tGmSqzuHwhyh2+/H8/L30m5FzT0+VD4P6DSwtkv9r/8JK+rtbue2Bo0f0LnZXI0lq65TTfa2apr6m5Fjpz29oCJ6/v5soXn3Pcz5nC2KHKhyTn65dqpdnWO46MTTBd5m/29PQMI2VmQmb3UCwlxfVp6+0v1HpgZM7t4WJbWX19XbS3L++U4jd/vLiF80fb8f0Hefklx/7RueVHpTfqm2vfw1M8+Qk9Jfu+s8jL08y6c/d+Ljy/9gsUF+uOuG/e/te88JE/mx9++e4FPfbhGVh1Qie5dScsub4smu81851dgzx32+OP6zkeqHCFBIA3v+QpnLBqBaed1HvMG1wuV/q1Xo/u3nNw3v4sfT9L0ejfX6P4xb575+kf4+IGGssH7hmsuEfgyMRhHndq6e/3grM28PlvHrte7IKzHsPGx9T2RLB6C2QDwNyfWm+xrayhoeWfet104vH90d569kkMDo4c037h2Rv46IOx4tc+5lEdJb8WYMuZJ/Hjew8sup6nPT5X9jHrydazT+Lzt91fsf/o7+Npj8/xze+X3/Zj1pqeDqYnJjPxc1hO871mtpxZ+szKO5YAAAgHSURBVLW6GD0rWyuu3zhxdQc9XSuZODTB4JxZpFyuJ1Pj8fQnrucbPyr/dvT0J67P1PezWFkbr2Z27qY1fOWOX1Tsb6SxnO89qGdla9nv9+LzT2Hs0MQx+4JefP4pVfsZlftgU2+HLG8HTg0hzH6MvgD4Yq2LOGldN+1tS7/48X951hkl259x1oaKX9e9qr3imYDnPv7RS6rnCaetXdLX1Vq5n1ul/oV+b+eE2u8pUw/me80s9TU1V0/XSjbkSi+azep6sVIuu+hJx9Uv1cqmjX3H1Z81x/Me1NbayvZt/Vzz8t/mb/78aVzz8t9m+7b+JJccrKtAFmMcA14BfDCEcA3wwxjjzSlqed+rL2Blua1/K3j98yvvQv7Wlz6lZHtHO7zrFefP+/h//ee/vah63vWX8z9mPSn386v0c630PXa0wbbzNi7oKgiNqtxrZrGvpUqufPE5nFw80wmyv16snFf80W8tql1K5eNXlT4f7j2venqNK6mN430P6ljRxvq+rqQf3FtmZo7/LKuUBgdHqvYN5HI9/PCnD/KD3fs5a9NaHh6f5NYf/YoLzjzxyCeMz33jHnbc9QBbz3nMvDM8c33rB3v50h33senkR/FHFzxu0Xtk3Xn3r7npzj1sO3fjkVmO3XuG+Pyt9zE1dZiLLzgtMzNjpSzl5/rA0CH+x7/dxYGRPNvO3cCzzj6Z1d0dTTkzVkqp18xym70A98b188+MZfkQ2Me++GNu+8kAT/+t9U0zM5bl8WpGs+O1e8/QMX+3Gtli3oNSyeV6Sh6CM5BV4BtQtjhe2eJ4ZYvjlS2OV/0qF8jq6pClJElSMzKQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhLL/E79kiRJWecMmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUqsPXUB9SiEsA14HjAAzMQY35G4JJURQjgduAa4C9gI7I8xvjNtVZpPCGEV8G3gqzHGK1LXo/JCCAF4AXAI2AJcHWP8TtqqVE4I4Q3AacA+4AzgZTHGQ0mL0oI4Q3aUEEIX8CHgtTHGq4EnhxAuTFuVKlgDfDrG+O4Y42uA54cQzk1dlOZ1DbAzdRGqLITQBrwPeGeM8V3Ay4Cfp61K5YQQTgTeDLw6xvh24AQKkwvKAGfIjnU+8IsYY754+1bgIuDmdCWpnBjjd49qagUeTlGLFiaE8CIKv1dPBroTl6PKngK0AK8ufljdD/xT2pJUwRgwAfQCD1H4/fpJ0oq0YM6QHWs9MDLn9nCxTXUuhPBc4CsxxrtT16LSQghPBJ4QY/zfqWvRgpxK4UPqx2KM1wHPAl6StiSVE2McBt4A3BhC+BiwB9idtCgtmIHsWANAz5zbvcU21bEQwnOA5wCvTV2LKnouMB5C+CvgGcBTQwiXJ65J5Q0Dd8cYDxZvfwt4drpyVEkI4WwKgeyiGONlFNaRvS1pUVowD1ke63bg1BBCR/Gw5QXAPySuSRWEEC4Cngm8BjgphHBqjPH2xGWphBjjtbP/DiF0At0xxg8kLEmVfRtYG0JoizFOU5gx25W4JpW3ATgQY5wq3n4QOCVhPVoELy5eQgjhd4A/AQaBSc+yrF/FBfy3AN8rNp0A/H2M8WPJitK8QgiXAK8EVlIYr08lLkllFJcCbKXwfngKhQXjnrVXh4onYXwQGKewhuxJwOUxxgeTFqYFMZBJkiQl5hoySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsx9yCQ1pRDCWcAHKFwq7bIQQitwD3BejHGowtddBjy7uPEmC/06SarEGTJJTSnG+APg43NuHwa2LTZULfXrJGku9yGTVPdCCH8BvBX4TyBPYcPLzwJTxfZ/BR5H4bI+LwN2Au+gsFP5ycAnY4z/WXys64DzgJ8Bh4FVxRmylwNXA38aY/x6COE04HrglxR2qI/AvwD/SGFH9JuBL1DYLPXor7sOuL/Y988xxv8bQnhFsdZ/K7Y/Cbgixvj55f55ScoeZ8gk1b0Y44coHBZ8MMb4KuAPgCspXFvxq8BpMcZLgIuB+4AbgI/EGN8MvAL4XyGER4UQ/hD4feD3Yoz/lTnXrY0x/lPxOWb9K/CpGOMbgRcBZ8UY7yk+9rdjjK+KMX6lzNd9uvh1rwD+JYTw6BjjPxZr7YsxPg94FfDmZfwxScowA5mkLLkVIMY4RuFyWc8ptt9UbL8NuJfChcsvCyF8CHh3se1kCpcA+lbxMCPAN0s9SQihB3j6nOebiDH+8XzFlfi6IWAPhWutzrqt+P97gUfP+x1Lagou6pfUCPIl2q6MMf4aIISwCpiobUllzdY6jR+KJRX5ZiApS84HCCF0UVgH9rWj7xBjHKEwQ/W7xfu2Al8COoAdwDOLbVCYSTvGnMe4oPgYnSGEG4vd40BbCKElhPCSeb6uD9hImZk4SZrlon5JmRBC+DrwXWA1EIDPAD+ncEhyL/DB2QXyIYTHAu+lsLD+BOAzMcYvF/uuA84FdgN9wNnA24BHUVic/x0K67tWUljUvxdYC/zPGOOtIYQNFE4o2E0h4LVV+LqNFNayfbW4fm221j8D3gI8H3hrjPFvl/nHJSljDGSSMqEYyK6OMX49cSmStOw8ZCmp7oUQ/hLoB14XQtiYuh5JWm7OkEmSJCXmDJkkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElK7P8Dzp2JLk7R0twAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(1,1,figsize=(10, 6))\n",
    "ax.scatter(yhat, df['count'])\n",
    "ax.set_xlabel('prediction')\n",
    "ax.set_ylabel('truth')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TODO** \n",
    "\n",
    "- add Vuong tests comparing ZINB to ordinary NB regression\n",
    "\n",
    "- add boostrapped standard errors"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}