_3_poisson_statsmodesl_zero_inflation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Zero-inflated Poisson regression and some examples with `statsmodels`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The zero-inflated poisson model is a natural choice to model phenomena that are counts but that are characterized by an over-abundance of zeros compared what one would expect from a poisson random-variable. In general one can model this situation is with a mixture between two distributions. In particular (see Cameron and Trivedi (2013), section 4.6), given a base count density $f_2(y_i)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\mathbb{P}[y_i=j] = \\begin{cases} \\pi_i + (1-\\pi_i)f_2(0) \\;\\; \\text{if} \\; j = 0, \\\\\n",
    "(1-\\pi_i)f_2(j) \\;\\;\\;\\;\\;\\;\\;\\;\\; \\text{if} \\; j > 0\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence we add a proportion of zeros $\\pi_i$ to the baseline distribution and the probabilities from the baseline model are decreased by the proportion $(1-\\pi_i)$ to make the probabilities add to one. One can either set $\\pi_i$ to a constant (just estimate an intercept) or let it depend on the $\\mathbf{x}_i$ e.g. via a logistic regression model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model is estimated by maximum likelihood. Define a censoring indicator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "d_i = \\begin{cases}1 \\;\\;\\;\\; \\text{if } y_i > 0, \\\\\n",
    "0 \\;\\;\\;\\; \\text{if } y_i = 0, \n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $d_i= 0$ with probability $\\pi_i$ and $d_i=1$ with probability $1-\\pi_i$. Hence the density becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(y_i) = [\\pi_i + (1-\\pi_i)f_2(0)]^{1-d_i}\\times[(1-\\pi)f_2(y_i)]^{d_i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we let the zero-inflation component $\\pi_i$ and the base count density $f_2(y_i)$ be functions of the parameter vector $\\mathbf{x}_i$, the likelihood can be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "L(\\mathbf{\\theta}_1,\\mathbf{\\theta}_2|\\mathbf{X},\\mathbf{Y}) &=  \\prod_{i=1}^m[\\pi_i(\\mathbf{x}_i, \\mathbf{\\theta}_1) + (1-\\pi_i(\\mathbf{x}_i, \\mathbf{\\theta}_1))f_2(0|\\mathbf{x}_i,\\mathbf{\\theta}_2)]^{1-d_i}\\times[(1-\\pi(\\mathbf{x}_i, \\mathbf{\\theta}_1))f_2(y_i|\\mathbf{x}_i,\\mathbf{\\theta}_2)]^{d_i}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One particular version of such a model is to make $f_2(y_i)$ poisson and model $\\pi_i$ with a logit. This is the model considered by Lambert (1992). In particular"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\mathbf{\\pi}_i(\\mathbf{x}_i,\\mathbf{\\theta}_1) &= \\frac{\\mathrm{e}^{\\mathbf{\\theta_1}'\\mathbf{x}_i}}{1+\\mathrm{e}^{\\mathbf{\\theta_1}'\\mathbf{x}_i}} \\\\\n",
    "f_2(y_i|\\mathbf{x}_i, \\mathbf{\\theta}_2) &= \\frac{\\mathrm{e}^{y_i\\mathbf{\\theta}_2'\\mathbf{x}_i}\\mathrm{e}^{-\\mathrm{e}^{\\mathbf{\\theta}_2'\\mathbf{x}_i}}}{y_i!} \n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For details on gradient and Hessian of the likelihood, see https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Zero-Inflated_Poisson_Regression.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cameron and Trivedi (2013) suggest that economic data still suffer from over-dispersion and therefore recommend to use the Negative Binomial model for the base count density. The problem is that this more heavily parameterized model is more difficult to estimate. Hence the convergence of the ML estimation may be slow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1) Example one: weekend drinking\n",
    "\n",
    "The is a nice example since it gives an explicit motivation for using a mixture. Ask a group of students: \"How many drinks did you have this weekend?\" The group may contain drinkers and non-drinkers. The non-drinkers never drink (so they are part of the $\\pi_i$), but even some of the drinkers may not have had a drink this weekend, so they may also report a zero (their poisson RV had a draw of zero). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://bookdown.org/roback/bookdown-bysh/ch-poissonreg.html#cs:drinking"
   ]
  },
  {
   "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 ZeroInflatedPoisson as ZIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('./weekendDrinks.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZUAAAGqCAYAAABzrqrBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfbidZ10n+m+bHV4CyUw67grIOZ0K9edAC+OgMqW8CK2AVM9gQdBRxIHROQjX4WWw1LYyRVNIaxFBPRbPwAhyRqnUN6jAlDIglEKhVdpcTm+gFuhgx2whnMTGljTk/LGeTXfSJPtOsvdae+/1+VxXr6zn/bfW82Tlfr691/0ct2/fvgAAAAAAQI/jJ10AAAAAAACrh1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmWAJFX1f1TVLVX1kUnX0quq/rCq7qqqH5h0LQAAMAlV9edruT1cVY+rqr+qqi92rPv9vesCHCuhMkCS1tqfJdk66TqORGvtx5L8r0nXAQAAE/TjST466SKWS2vthiSv6Fz3+t51AY7VzKQLAAAAADgarbWdk64BYBoJlYGpU1WnJPmNJBsy+h78ldba+4fFx1XVm5M8LskDkvxca+3Gqjo+yW8lOS3JN5N8LsnLW2t3VtWJSS5P8m3D/i5trf1JVX1/kt9J8k+Hbc9O8uQk25KckuQPWmv/rqpeleT8JO9qrb2iqp6R5KIk30iyM8l/aK397VD7DyW5JMnXkvy3ZfuQAABgDKrqtUl+Psn7MmpPf0eS7Ule2Fr7+2GdZyZ5bZJ7kvxjkpe11j5fVb+Q5D8muby1dlFVfXuS382oHb8+yXtba5dU1YZh/rcnWZfk0621Vw77Pui9wUHa8j+cZHOS57bWPneY9/GeYZvvT/LHSa7MqK3/mCSvaq39ybD+xiS/nqSGmv4ko/uIfVW1Lskbk5yZ5PYkf3HAsR6c5C1JviujX6C/s7V2+ZF98gDHxvAXwFSpqpkk780o0H1ykp9O8u6qesSwyuOT/OfW2hOT/FqSP6qq9UmemeSft9aeOGz3z5LMDtv8v0m2DfOfk+RtVfXPF/z87GFJ/rK19gMZNQ6fOmz3S8Ofb0nysSFQPjmjhujPtNaekuQDSd451P5tSf4wo5D5BzIKp79jaT8hAAAYn9baL2fU5v2BJC9orX1fkq9k1EZOVX1n7m0fPznJ7yV5b1XNtNZ+ddh23n9M8pHW2lOTPCPJjwzzfybJV4f29VOSPHHY9yHvDQ5oy183bPvx4RiHex//OsmLkpyR5JVJnt5ae3aSV2f/4fZ+Pcm64b7jaUl+MslPDcv+Q0b3DN+XUceUf3nA4d60YNtnJDm3qp540A8YYJkIlYFp8/gk35nkXUnSWvtCkk9l1IgbZrWbh9fvzqgReXqSHUlOq6ofHHot/0SSL1fVdyQ5K8nbho3vyKix+RMLjnlna+1Dw/JfaK19NckHk7xgWP7M3NsY/rdJPtNaa8P0f01yZlU9NKMG5d+11q4b9vVnGfXUAACA1e6q1tqu4fXvJXnu0GP3J5Jcv6B38O8nOSnJEw6yj68l+aGqenRr7c4kT18w/0lV9fjW2t6MguVk8XuDJPmH1trHh9c3JTl5kffx0dbaN1prcxn1uP7sgm2/M0mG+4mfTPL24bj/mNG9x78b1v2xJO9prd3VWtuX5Ir5nQ/bvmDBtrsyCsbn7y0AxkKoDEybhyfZ0Vq7Z8G8uWF+MgqPkyRDg/PrSR46BLk/l+Q1Sb6UUU+D4xZs986q+khVfSTJI5I8aMH+/7+D1PHOjHpCJMnzMmpEztf3qAX7+uPheN+e5KFJ/v6A/Xxt8bcMAAAr3o4Fr7+a0fAV35ZR+3hufsHQRt+Re9vhC/1qRsNNvLuq/iqjThlprf1BRr2E31xVn0/ywmH9xe4NktFwdPPuSnK/Rd7HrgWv71kwfc/wnpLRLx7vv/B9HXDcA9v9C9v889teuuCe4SnDPICxMaYyMG1uT7J5+LncfONxNsktw+sT5lccfg73T5PcUVX/JKOf0v35MFTGBzL6Wd4Hh9WfO/RGSFU9IIt/v74vyf9TVU9Pktba1xfU95nW2tkL6ticUWP2sbl3yI15JwQAAFa/he3ab0uyJ6Ng9faMxh1Okgy9lzcn+Z8H2ceJrbXfSPIbVXVWkvdV1Y0ZdfJ4d2vtXVX1PUk+VFW3ZPF7g+Uyl+Tu4Vj/Y8Fx59/THdm/3f/PDrLty1prn06SYbi+DctZMMCB9FQGps2nknwho2Em5sdoe3xG4yIno17Cpw2vfzzJ3ya5LsmPZtRTOa21WzNq8K0bHqB3dfb/udnluXfc5INqrX0jo97J/yWjcZLn/X6Sx1fVSUN9Jyb5SEbf11clObGqzhiW/ZskDz6idw8AACvT04eH1yWjX/S9Z+iV/PtJvreqHjkse35Gv+T7xEH28Yaqmh9/+FMZPfj6uCQvy9BrOcnNGfX8XZfF7w2WRWvtmxn9cvGFw3EfmNGvF//LsMoVGQ3/8YCqOi4LhtZbsO3C+48Lc++vIAHG4rh9+/ZNugaAsRp6Gv9mRkNUzCT5lYx+inZpRr0h/keSR2X01Oifa63dUFWV0YP7HphRkHtTkp9vrX1jCH7/74yGqEiS97fWXl9Vj8poTOTvTvLJJOe01r7107WqenxG4589bOFP7obey/8po94Z30xyfmvtk8OyZyW5JKOf/H08owbw14c6r1/SDwoAAMagqn43yZ0ZtadPyqg37k+31v5+WD7fPt6b0TNFXtZa+3xV/UJGD867K6M2/f/KaLi6e5L8kyTvaK29par+dZKLMwqYNyW5urX2i8O+73Nv0Fp7/wFt+Xdk9AyVdyZ5SJLfaa2de8B7eFWSc4daXpnkBzN6YN8tGQXh78gosL66tfb0qnpwRg/r++7huH+c5NLW2r7hF5NvTHJmRp1cPjns+89ba+cs2PZRGd0z/OXwOTwuye8M+3xva+3HjuqEAHQQKgMAAAATM4TKX2ytXTThUgDoZPgLAAAAAAC6CZUBAACAiaiq1yZ5ZpKfqaoXT7oeAPoY/gIAAAAAgG56KgMAAAAA0G1mnAebm9s19m7RmzdvyI4du8d9WFYg1wLzXAvMcy0wz7Vw7GZnNx436RqYjEm08Xv5uz2dnPfp5dxPJ+d9ejn3y+9wbfw131N5ZmbdpEtghXAtMM+1wDzXAvNcC7A2+bs9nZz36eXcTyfnfXo595O15kNlAAAAAACWjlAZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALrNTLqA5XThha9Ja3+dPXv2TrqU+zj11NOyZcslky4DAAAAAOCIrOlQedu2m3P9DTdm0+zJky5lPzvnbpt0CQAAwIS9aOuHl3yfbz/vaUu+TwCAA63pUDlJNs2enCc87+JJl7GfT1xxwaRLAAAAAAA4KsZUBgAAAACg25rvqQwAAHCklnpoCsNSAABriZ7KAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0m5l0AQAAwNKpqrOSnJNke5J9rbXXHbD8NUkekuSOJN+b5LWttVuGZV9M8sVh1a+01n5yPFUDALCaCJUBAGCNqKoNSS5P8ujW2t1VdWVVndlau2bBag9O8qrW2r6qen6SX03yI8Oy322tXTTeqgEAWG2EygAAsHacnuRLrbW7h+lrk5yd5Fuhcmvtlxasf3ySf1gw/aSqOjfJxiTvb619YpnrBQBgFRIqAwDA2nFikl0LpncO8+6jqu6X5IVJXrpg9i+21q4fejzfWFU/3Fr7wuEOuHnzhszMrDvGspfP7OzGSZeQZHx1rJT3O2k+h+nl3E8n5316OfeTI1QGAIC1Y3tGvYznbRrm7WcIlH87yQWttVvn57fWrh/+3F1Vf5XkjCSHDZV37Ni9BGUvj9nZjZmb27X4imMwrjpWyvudpJV03hkv5346Oe/Ty7lffocL7Y8fYx0AAMDyui7JSVV1/2H6jCRXVdUJVbUp+da4y29N8muttRuq6jnD/DOr6pkL9vXIJLcGAAAOoKcyAACsEUMP45ckeUtVzSW5qbV2TVVdmuRrSbYmeVeSU5OcXFVJ8qAkV2bUo/miqvpXSR6W5I9aax+fxPsAAGBlEyoDAMAa0lq7OsnVB8w7d8Hrcw6x3c1JnrO81QEAsBYY/gIAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKDbTO+KVfXAJJ9K8t9aa6+uqgckuSzJV5KckmRra+1zy1MmAAAAAAArwZH0VN6S5C8XTL8iyZdba29I8qYkb1vKwgAAAAAAWHm6QuWqekGSa5PctmD22UmuS5LW2s1JHltVm5a8QgAAAAAAVoxFh7+oqkcl+RettfOr6jELFp2YZNeC6Z3DvJ2H2tfmzRsyM7PuaGs9YuvXj+9YR2r9+nWZnd046TKmjs+cea4F5rkWmOdaAAAA6NMzpvKPJrmrqs5L8sQk96uqVyTZnmTh3demYd4h7dix+2jrPCp79uwd6/GOxJ49ezM3t2vxFVkys7MbfeYkcS1wL9cC81wLx04oDwAA02PRULm1dvH86+HhfA9urf368Pr0JB+rqtOSfLa1dsheygAAAAAArH49PZWTJFX1nCRPzqin8k8keXOSy6rqwiSPTPLi5SkRAAAAAICVojtUbq1dmeTKA2a/dGnLAQAAAABgJTt+0gUAAAAAALB6CJUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoNvMpAsAAACWTlWdleScJNuT7Gutve6A5a9J8pAkdyT53iSvba3dMiz7qSTfk2Rvkltba28dZ+0AAKwOeioDAMAaUVUbklye5JWttYuSPKaqzjxgtQcneVVr7dIkVyb51WHbhyd5dZJXt9bOTfLvq+qUsRUPAMCqoacyAACsHacn+VJr7e5h+tokZye5Zn6F1tovLVj/+CT/MLx+RpIbWmv7hunrkvxQks8va8UAAKw6QmUAAFg7Tkyya8H0zmHefVTV/ZK8MMlLj3TbhTZv3pCZmXVHVew4zM5unHQJScZXx0p5v5Pmc5hezv10ct6nl3M/OUJlAABYO7YnWXh3tWmYt58hUP7tJBe01m5dsO0jD9j2C4sdcMeO3Udd7HKbnd2Yubldi684BuOqY6W830laSeed8XLup5PzPr2c++V3uNDemMoAALB2XJfkpKq6/zB9RpKrquqEqtqUfGvc5bcm+bXW2g1V9Zxh3Q8meVxVHTdMn57k/WOsHQCAVUJPZQAAWCNaa7ur6iVJ3lJVc0luaq1dU1WXJvlakq1J3pXk1CQnV1WSPCjJla21/1lVlyV5U1XtTfKfW2vGUwYA4D6EygAAsIa01q5OcvUB885d8Pqcw2z7roxCZwAAOCTDXwAAAAAA0E2oDAAAAABAN6EyAAAAAADdhMoAAAAAAHQTKgMAAAAA0E2oDAAAAABAN6EyAAAAAADdhMoAAAAAAHQTKgMAAAAA0E2oDAAAAABAN6EyAAAAAADdhMoAAAAAAHQTKgMAAAAA0E2oDAAAAABAN6EyAAAAAADdhMoAAAAAAHQTKgMAAAAA0E2oDAAAAABAN6EyAAAAAADdhMoAAAAAAHSbmXQBAAAALJ8Xbf3wku7v7ec9bUn3BwCsPnoqAwAAAADQTagMAAAAAEA3oTIAAAAAAN0WHVO5qo5P8t4kn0pyvySPSPKiJA9MsjXJ3yQ5Jcn5rbW/W75SAQAAAACYtN6eyte11n65tXZhkg1Jzkny+iQfaq1tTfInSS5bphoBAAAAAFghFu2p3Fr7ZpItSVJVM0kenqRl1Ev54mG1a5O8Y7F9bd68ITMz64662CO1fv34jnWk1q9fl9nZjZMuY+r4zJnnWmCea4F5rgUAAIA+i4bK86rqGUlemeR9rbXPVNWJSXYNi3cm2VxVM621ew61jx07dh9TsUdqz569Yz3ekdizZ2/m5nYtviJLZnZ2o8+cJK4F7uVaYJ5r4dgJ5QEAYHp0P6ivtfbB1tozk5xcVT+fZHuS+buHTUl2HC5QBgAAAABg9Vs0VK6qR1XV2Qtm3ZbkO5NcleT0Yd4ZwzQAAAAAAGtYz/AXdyd5cVV9T5L1Sf5Fkv8ryTeSXFJV35XkEUlevWxVAgAAAACwIvQ8qO/WJOccYvHPLm05AAAAAACsZN1jKgMAAAAAgFAZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALrNTLoAAABg6VTVWUnOSbI9yb7W2usOss7zk7w+yctba+9bMP+TSe4aJve21s4cQ8kAAKwyQmUAAFgjqmpDksuTPLq1dndVXVlVZ7bWrlmwzskZBc63H2QXH2itXTSeagEAWK2EygAAsHacnuRLrbW7h+lrk5yd5FuhcmvttiS3VdV/Osj2p1XVa5I8MMmnW2tXLXfBAACsPkJlAABYO05MsmvB9M5hXq9LWmvXV9W6JH9RVbtaa39xuA02b96QmZl1R1HqeMzObpx0CUnGV8c4jrNSPtPDWQ01sjyc++nkvE8v535yhMoAALB2bE+y8O5q0zCvS2vt+uHPvVX1sSRPTXLYUHnHjt1HUeZ4zM5uzNzcrsVXHINx1TGO46yUz/RQVtJ5Z7yc++nkvE8v5375HS60P36MdQAAAMvruiQnVdX9h+kzklxVVSdU1abDbVhV311VL14w65Qkty5TnQAArGJ6KgMAwBrRWttdVS9J8paqmktyU2vtmqq6NMnXkmytquOSXJDkpCTPr6o9rbUPZjRUxtlV9bCMejjfnuS/TuadAACwkgmVAQBgDWmtXZ3k6gPmnbvg9b4kW4b/Fq7zt0nOGUeNAACsboa/AAAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACg28ykCwAAADgSL9r64SXd39vPe9qS7g8AYK3TUxkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG5CZQAAAAAAus0stkJVPSLJliQ3Jnl4kq+21n65qk5IsjXJ3yQ5Jcn5rbW/W85iAQAAAACYrEVD5SQnJPmD1tqfJklV/XVVXZXkZ5N8qLV2RVX9SJLLkrxg+UoFAAAAAGDSFg2VW2ufPmDW8UnuTHJ2kouHedcmecdi+9q8eUNmZtYdaY1Hbf368R3rSK1fvy6zsxsnXcbU8Zkzz7XAPNcC81wLAAAAfXp6Kn9LVf1okg+21m6pqhOT7BoW7UyyuapmWmv3HGr7HTt2H32lR2HPnr1jPd6R2LNnb+bmdi2+Iktmdnajz5wkrgXu5Vpgnmvh2AnlAQBgenQ/qK+qnprkqUleOczanmT+7mFTkh2HC5QBAAAAAFj9unoqV9XZSZ6U5OVJHlpVJyW5KsnpSW5PcsYwDQAAAADAGrZoqFxVj0vy7iSfSfLfkzwoyW8lOT/JJVX1XUkekeTVy1gnAAAAAAArQM+D+m5I8uBDLP7ZpS0HAAAAAICVrHtMZQAAAAAAECoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQLeZSRcAAAAsnao6K8k5SbYn2ddae91B1nl+ktcneXlr7X1Hsi0AAOipDAAAa0RVbUhyeZJXttYuSvKYqjrzgHVOzig0vv1ItwUAgESoDAAAa8npSb7UWrt7mL42ydkLV2it3dZa++9Hsy0AACSGvwAAgLXkxCS7FkzvHOYt27abN2/IzMy67gJXotnZjWviGOM6zrjey7FYDTWyPJz76eS8Ty/nfnKEygAAsHZsT7Lw7mrTMG/Ztt2xY3d3cSvV3NyuxVdaBccY13HG9V6O1uzsxhVfI8vDuZ9Ozvv0cu6X3+FCe8NfAADA2nFdkpOq6v7D9BlJrqqqE6pq09Fsu0x1AgCwigmVAQBgjWit7U7ykiRvqaotSW5qrV2T5LwkP58kVXVcVV2Y5KQkz6+qZyyyLQAA7MfwFwAAsIa01q5OcvUB885d8Hpfki3Df4tuCwAAB9JTGQAAAACAbkJlAAAAAAC6CZUBAAAAAOgmVAYAAAAAoJtQGQAAAACAbkJlAAAAAAC6zUy6AAAAAFa3F2398JLu7+3nPW1J9wcALC09lQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6CZUBgAAAACg28ykC2BlufDC12TbtpsnXcZBnXrqadmy5ZJJlwEAAAAAU02ozH62bbs5199wYzbNnjzpUvazc+62SZcAAAAAAESozEFsmj05T3jexZMuYz+fuOKCSZcAAAAAAMSYygAAAAAAHAGhMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0m1lshap6SJItSR7bWvu+Yd4DklyW5CtJTkmytbX2ueUsFAAAAACAyevpqfzEJH+a5LgF816R5MuttTckeVOSty1DbQAAAAAArDCLhsqttfck2XXA7LOTXDcsvznJY6tq09KXBwAAAADASrLo8BeHcGL2D5p3DvN2Hm6jzZs3ZGZm3VEe8sitXz++Yx2p9evXZXZ246TLuI+1/pmtxM+cyXAtMM+1wDzXAgAAQJ+jDZW3J1l457VpmHdYO3bsPsrDHZ09e/aO9XhHYs+evZmbO7AD+OSt5c9sdnbjivzMGT/XAvNcC8xzLRw7oTwAAEyPnjGVD+aqJKcnSVWdluSzrbXD9lIGAAAAAGD1W7SnclU9JckLkjy0qi5M8sYkb05y2TD9yCQvXtYqAQAAAABYERYNlVtrH03y0YMseunSlwMAAAAAwEp2tMNfAAAAAAAwhYTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQLdFH9QHAACsHlV1VpJzkmxPsq+19roDlj8gyWVJvpLklCRbW2ufG5Z9MckXh1W/0lr7yfFUDQDAaiJUBgCANaKqNiS5PMmjW2t3V9WVVXVma+2aBau9IsmXW2uXVtVpSd6W5EnDst9trV003qoBAFhtDH8BAABrx+lJvtRau3uYvjbJ2Qesc3aS65KktXZzksdW1aZh2ZOq6tyq+pWqesJYKgYAYNXRUxkAANaOE5PsWjC9c5jXs87OJL/YWrt+6PF8Y1X9cGvtC4c74ObNGzIzs+7YK5+g2dmNa+IY4zrOajjGuD5vVh7nfjo579PLuZ8coTIAAKwd25MsvLvaNMzrWqe1dv3w5+6q+qskZyQ5bKi8Y8fuYyx58ubmdi2+0io4xriOs9KPMTu7cWyfNyuLcz+dnPfp5dwvv8OF9oa/AACAteO6JCdV1f2H6TOSXFVVJywY4uKqjIbJyDCm8mdbazur6syqeuaCfT0yya3jKhwAgNVDT2UAAFgjhh7GL0nylqqaS3JTa+2aqro0ydeSbE3y5iSXVdWFGQXHLx42357koqr6V0keluSPWmsfH/+7AABgpRMqAwDAGtJauzrJ1b9Xr5QAAAleSURBVAfMO3fB639M8tKDbHdzkucse4EAAKx6hr8AAAAAAKCbUBkAAAAAgG5CZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG4zky5gGt359TuybeftefaznzXpUu5j27abMrPpf5t0GQAAAADACiVUnoC937gru49LbvnyjkmXch87d/1DTtg06SoAAAAAgJVKqDwhm2ZPzhOed/Gky7iPD/zmv510CQAAAADACmZMZQAAAAAAugmVAQAAAADoJlQGAAAAAKCbUBkAAAAAgG4e1AcAAMCK96KtH17S/b39vKct6f4AYJoIlVkV7vz6Hdm28/Y8+9nPOup9rF+/Lnv27F3Cqu516qmnZcuWS5Zl3wAAAACwkgiVWRX2fuOu7D4uueXLOyZdyn3snLtt0iUAAAAAwNgIlVk1Ns2enCc87+JJl3Efn7jigkmXAAAAAABj40F9AAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3WYmXQCsdnd+/Y5s23l7nv3sZ026lPs49dTTsmXLJZMuAwAAAIA1RKgMx2jvN+7K7uOSW768Y9Kl7Gfn3G2TLgEAAACANUioDEtg0+zJecLzLp50Gfv5xBUXTLoEAAAAANYgYyoDAAAAANBNqAwAAAAAQDehMgAAAAAA3YTKAAAAAAB0EyoDAAAAANBtZtIFANPlwgtfk23bbp5oDevXr8uePXsPuuzUU0/Lli2XjLkiAABWihdt/fCS7u/t5z1tSfcHACuBUBkYq23bbs71N9yYTbMnT7qU+9g5d9ukSwAAAABY8YTKwNhtmj05T3jexZMu4z4+ccUFky4BAAAAYMUzpjIAAAAAAN2EygAAAAAAdBMqAwAAAADQTagMAAAAAEA3oTIAAAAAAN1mJl0AAMBqdOGFr8m2bTdPuoyDOvXU07JlyyWTLgMAAFijhMoAAEdh27abc/0NN2bT7MmTLmU/O+dum3QJAADAGidUBgA4SptmT84TnnfxpMvYzyeuuGDSJQAAAGucMZUBAAAAAOimpzIAAACM0Yu2fnhJ9/f28562pPsDgMXoqQwAAAAAQLdj6qlcVWclOSfJ9iT7WmuvW5KqAACAo7JYG72qHpDksiRfSXJKkq2ttc8Ny34qyfck2Zvk1tbaW8dZOwAAq8NR91Suqg1JLk/yytbaRUkeU1VnLlVhAADAkelso78iyZdba29I8qYkbxu2fXiSVyd5dWvt3CT/vqpOGVvxAACsGsfSU/n0JF9qrd09TF+b5Owk1xxzVUto59xtK+4p6PfsuWtF1pWs3NpWal3Jyq1t59xt2bbz9jz72c+adCn72bbtpuy++54V93klK/czY/msX78ue/bsnXQZrABHcy2s1O+znXO3Jf/75kmXweT0tNHPTnJ+krTWbq6qx1bVpiTPSHJDa23fsN51SX4oyefHUjmwpIzbDMByOpZQ+cQkuxZM7xzmHdLs7MbjjuF4R+zaaz82zsMBAMCk9bTRD7XOEbfvk/G38ZPkvW/8N46xwo6zVo4xruOslWOsRLOzGyddAhPgvE8v535yjuVBfduTLDxzm4Z5AADAZPS00Q+1jvY9AABdjiVUvi7JSVV1/2H6jCRXHXtJAADAUTpoG72qThiGuEhGbfbTk6SqTkvy2dbaziQfTPK4qprveXx6kvePr3QAAFaL4/bt27f4WodQVT+Y5LlJ5pLsOfDJ0gAAwHgdrI1eVZcm+VprbWtVPTDJZUnuSPLIJK9vrX1u2Panknxvkr1JPtdae+tE3gQAACvaMYXKAAAAAABMl2MZ/gIAAAAAgCkjVAYAAAAAoJtQGQAAAACAbjOTLmC5VNVZSc5Jsj3JPg8RnG5V9ckkdw2Te1trZ06yHsanqh6SZEuSx7bWvm+Y94CMHlD0lSSnJNk6/4Ai1q5DXAs/k+T/zL3fD29rrf3eZCpkHKrqERldBzcmeXiSr7bWfrmqTkiyNcnfZPS9cH5r7e8mVylwLNwLTC/t/umhnT+dtOmnkzb8yrQmQ+Wq2pDk8iSPbq3dXVVXVtWZrbVrJl0bE/OB1tpFky6CiXhikj9N8i8XzHtFki+31i6tqtOSvC3JkyZRHGN1sGshSX68tfbF8ZfDhJyQ5A9aa3+aJFX111V1VZKfTfKh1toVVfUjGd2QvmCCdQJHyb3A1NPunx7a+dNJm346acOvQGt1+IvTk3yptXb3MH1tkrMnWA+Td1pVvaaqLqoq18IUaa29J8muA2afneS6YfnNSR5bVZvGXRvjdYhrIUleVlWvrqrXDv+nmzWstfbp+cbo4Pgkd2bB90K0G2C1cy8w3bT7p4R2/nTSpp9O2vAr01oNlU/M/l8yO4d5TK9LWmuXJPmVJOdX1ZMnXRAT5TuCeR/N6PvhsiSfSfKHE66HMaqqH03ywdbaLdn/e2Fnks1VtSZ/0QVTwL/z0027f7r5+z+dtOmniDb8yrFWQ+XtSTYumN40zGNKtdauH/7cm+RjSZ462YqYMN8RJElaa7e11uaGyQ8neUpVrZtkTYxHVT01o38LXjnMWvi9sCnJjtbaPZOoDThm/p2fYtr9U8/f/ymkTT89tOFXlrUaKl+X5KSquv8wfUaSqyZYDxNUVd9dVS9eMOuUJLdOqh5WhKsy+mlshrHWPtta2znZkpiEqnrDgv+TfUqSLw43oaxhw8+hn5Hk5UkeUlWnZ8H3QrQbYLVzLzCltPuJdv5U0qafDtrwK89x+/btm3QNy6KqfjDJc5PMJdnjic/Tq6oeluQ3k/xlRv/nan2SV7XWvjnRwhiLqnpKkp9O8swkv53kjcOiy5LckeSRSV7vqdBr3yGuhZ9LcmqS25KcluTNrbVPTqxIll1VPS6jn0h+Zpj1oCS/leTPklyS5EtJHpHkPE+OhtXLvcB00u6fLtr500mbfjppw69MazZUBgAAAABg6a3V4S8AAAAAAFgGQmUAAAAAALoJlQEAAAAA6CZUBgAAAACgm1AZAAAAAIBuQmUAAAAAALoJlQEAAAAA6Pb/A94EVgjyLdMUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1800x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots(1,2, figsize=(25,7))\n",
    "ax[0].hist(df.drinks, bins=15, edgecolor=\"k\", linewidth=1.5)\n",
    "ax[0].set_title('observed')\n",
    "\n",
    "lamb = df.drinks.mean()\n",
    "theoretical = np.array([(lamb**k)*np.exp(-lamb) / np.math.factorial(k) for k in range(df.drinks.max()+1)])\n",
    "\n",
    "ax[1].bar(np.arange(0,df.drinks.max()+1), theoretical)\n",
    "ax[1].set_title('poisson model')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we can fit a simple Poisson model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['off_campus'] = (df.dorm == 'off campus')\n",
    "df['male'] = (df.sex == 'm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "endog_var = df.drinks.values\n",
    "exog_vars = np.concatenate([np.ones(df.shape[0]).reshape(-1,1),\n",
    "                            df.off_campus.values.reshape(-1,1),\n",
    "                            df.male.values.reshape(-1,1)],\n",
    "                           axis=1)\n",
    "m = sm.GLM(endog = endog_var,\n",
    "           exog = exog_vars,\n",
    "           family = sm.families.Poisson()).fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>           <td>y</td>        <th>  No. Observations:  </th>  <td>    77</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    74</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>         <td>Poisson</td>     <th>  Df Model:          </th>  <td>     2</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>          <td>log</td>       <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -175.55</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Fri, 17 Apr 2020</td> <th>  Deviance:          </th> <td>  230.54</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>00:26:38</td>     <th>  Pearson chi2:      </th>  <td>  254.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>5</td>        <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>    0.1293</td> <td>    0.124</td> <td>    1.041</td> <td> 0.298</td> <td>   -0.114</td> <td>    0.373</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    0.8976</td> <td>    0.201</td> <td>    4.470</td> <td> 0.000</td> <td>    0.504</td> <td>    1.291</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>    <td>    1.1154</td> <td>    0.161</td> <td>    6.925</td> <td> 0.000</td> <td>    0.800</td> <td>    1.431</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   No. Observations:                   77\n",
       "Model:                            GLM   Df Residuals:                       74\n",
       "Model Family:                 Poisson   Df Model:                            2\n",
       "Link Function:                    log   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -175.55\n",
       "Date:                Fri, 17 Apr 2020   Deviance:                       230.54\n",
       "Time:                        00:26:38   Pearson chi2:                     254.\n",
       "No. Iterations:                     5                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.1293      0.124      1.041      0.298      -0.114       0.373\n",
       "x1             0.8976      0.201      4.470      0.000       0.504       1.291\n",
       "x2             1.1154      0.161      6.925      0.000       0.800       1.431\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "230.5439891062946\n"
     ]
    }
   ],
   "source": [
    "print(m.deviance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the code here https://github.com/proback/BYSH/blob/master/04-Poisson-Regression.Rmd, the `first_year` variable is generated as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['first_year'] = df.dorm.isin(['kildahl', 'mohn', 'kittlesby'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 1.828103\n",
      "         Iterations: 19\n",
      "         Function evaluations: 20\n",
      "         Gradient evaluations: 20\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.drinks.values\n",
    "exog_vars = np.concatenate([np.ones(df.shape[0]).reshape(-1,1),\n",
    "                            df.off_campus.values.reshape(-1,1),\n",
    "                            df.male.values.reshape(-1,1)],\n",
    "                           axis=1)\n",
    "exog_vars_infl = np.concatenate([np.ones(df.shape[0]).reshape(-1,1),\n",
    "                                 df.first_year.values.reshape(-1,1)],\n",
    "                                axis=1)\n",
    "m = ZIP(endog = endog_var,\n",
    "        exog = exog_vars,\n",
    "        exog_infl=exog_vars_infl,\n",
    "        inflation='logit').fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ZeroInflatedPoisson Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>          <th>  No. Observations:  </th>  <td>    77</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>           <td>ZeroInflatedPoisson</td> <th>  Df Residuals:      </th>  <td>    74</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.1315</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>00:26:38</td>       <th>  Log-Likelihood:    </th> <td> -140.76</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>              <td>True</td>         <th>  LL-Null:           </th> <td> -162.08</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>      <th>  LLR p-value:       </th> <td>5.540e-10</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>   -0.6036</td> <td>    0.303</td> <td>   -1.992</td> <td> 0.046</td> <td>   -1.198</td> <td>   -0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_x1</th>    <td>    1.1364</td> <td>    0.609</td> <td>    1.865</td> <td> 0.062</td> <td>   -0.058</td> <td>    2.330</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>         <td>    0.7543</td> <td>    0.139</td> <td>    5.408</td> <td> 0.000</td> <td>    0.481</td> <td>    1.028</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>            <td>    0.4159</td> <td>    0.205</td> <td>    2.026</td> <td> 0.043</td> <td>    0.014</td> <td>    0.818</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>            <td>    1.0209</td> <td>    0.172</td> <td>    5.926</td> <td> 0.000</td> <td>    0.683</td> <td>    1.359</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                     ZeroInflatedPoisson Regression Results                    \n",
       "===============================================================================\n",
       "Dep. Variable:                       y   No. Observations:                   77\n",
       "Model:             ZeroInflatedPoisson   Df Residuals:                       74\n",
       "Method:                            MLE   Df Model:                            2\n",
       "Date:                 Fri, 17 Apr 2020   Pseudo R-squ.:                  0.1315\n",
       "Time:                         00:26:38   Log-Likelihood:                -140.76\n",
       "converged:                        True   LL-Null:                       -162.08\n",
       "Covariance Type:             nonrobust   LLR p-value:                 5.540e-10\n",
       "=================================================================================\n",
       "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "inflate_const    -0.6036      0.303     -1.992      0.046      -1.198      -0.010\n",
       "inflate_x1        1.1364      0.609      1.865      0.062      -0.058       2.330\n",
       "const             0.7543      0.139      5.408      0.000       0.481       1.028\n",
       "x1                0.4159      0.205      2.026      0.043       0.014       0.818\n",
       "x2                1.0209      0.172      5.926      0.000       0.683       1.359\n",
       "=================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = m.predict(exog=exog_vars, exog_infl=exog_vars_infl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAFyCAYAAAAkp9RWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAf20lEQVR4nO3de3Cdd33n8bcuthxFspHtY5LaSYAY/cgOwTgkTBNDkxi33SVlOsAOLmxJMwvsNky4DlNaAl1TkqFpQxYoDMwWdktYLmazu5O24dbgBNqQK/GkYZv+HEEutnGwbAtdIvtIlrR/nCNHkY9k2db5Pc85z/s1k4mf53cu36Pf8xx99Hsuv5apqSkkSZJUf61ZFyBJklQUBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKpD3rAhaiv3+4oe950dPTycDAaNZlaAb7JH/sk3yyX/LHPsmf2X1SKnW3zPVYR7wSaG9vy7oEzWKf5I99kk/2S/7YJ/lzMn1i8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSU2kPD7B/oFRyuMTWZciqYaGmKtRkjS/iclJtu/oY+eufg4NlVm5vIONvSW2bl5PW6t/Y0t5YfCSpCawfUcfdz6059jywaHyseW3benNqixJs/hnkCQ1uPL4BDt39dds27nrgIcdpRwxeElSgxscKXNoqFyzbWD4CIMjtdskpWfwkqQGt6Krg5XLO2q29XQvY0VX7TZJ6Rm8JKnBdSxpY2NvqWbbxt7VdCxpS1yRpLl4cr0kNYGtm9cDlXO6BoaP0NO9jI29q4+tl5QPBi9JagJtra28bUsvb778fAZHyqzo6nCkS8ohg5ckNZGOJW2s6enMugxJc/AcL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSZKUiMFLkiQpEYOXJElSIgYvSZKkRAxekiRJibTX64VDCOcDNwAPA+uAgzHGPwshrAT+HPg58FLgIzHGX9arDkmSpLyoW/ACVgLfjDHeDhBC+JcQwh3Au4A7Y4zfCiG8AbgZeHsd65AkScqFugWvGOODs1a1As8CVwE3VtfdA3ylXjVIkiTlScvU1FTd3ySE8Ebgihjj+0IIZeCFMcZfhRDagXFgSYzx6FzPP3p0Yqq9va3udUqSJC2Clrka6nmoEYAQwpXAlcD7q6v2A93Ar4DlwMB8oQtgYGC0rjXWW6nUTX//cNZlaAb7JH/sk3yyX/LHPsmf2X1SKnXP+di6Bq8QwlXAa4H3AWeHEM4D7gAuBXYDm6rLkiRJTa+eVzW+CtgOPATcBZwJfB74CHBTCKEXOB/4UL1qkCRJypN6nlz/E6BrjuZ31et9JUmS8sobqEqSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSZKUiMFLkiQpEYOXJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwkiRJSsTgJUmSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSZKUiMFLkiQpEYOXJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwkiRJSsTgJUmSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSZKUSHu9XjiEcBZwA7AhxnhJdd01wB8CR6oP+3KM8av1qkGSJClP6ha8gNcAtwOvnLX+92KMT9bxfSVJknKpbsErxnhbCOGKGk3XhRCeATqBz8UYD9WrBkmSpDyp54hXLT8E7ogx9ocQXg/8L+B1J3pST08n7e1tdS+unkql7qxL0Cz2Sf7YJ/lkv+SPfZI/C+2TpMErxvjEjMUdwN+GENpijBPzPW9gYLS+hdVZqdRNf/9w1mVoBvskf+yTfLJf8sc+yZ/ZfTJfCEt6VWMI4ZMhhOmw91LgyROFLkmSpGZRz6saLwfeDpwdQvgo8CngGeALIYQngAuB36/X+0uSJOVNPU+u/yGVc7pm+ky93k+SJCnvvIGqJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwkiRJSsTgJUmSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSZKUiMFLkiQpEYOXJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwkiRJSsTgJUmSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVySsErhPBfFrsQSZKkZtd+ogeEEP4TsA1YA7RU/5sCPl7XyiRJkprMQka83g9cDiyNMbbFGFuBj9a3LEmSpOZzwhEv4JEY4+Oz1n2nHsVIkiQ1szmDVwjhT6v/fDaE8APgPqBcXfd64NfrXJskSVJTme9Q4xuonM+1G/gRMMZz53i11L80SZKk5jLfocY/ijHeVashhPC9OtUjSZLUtOYc8ZoOXSGEd89cH0K4BnhrfcuSJElqPgu5qvFlMxdijH8DdNalGkmSpCY238n1T1C5X9fKEMLvzGhqA35a78IkSZKazXzneF1B5ST6bcDMO9UfiTH+so41SZIkNaU5g1eM8anqP6+Z3RZCeHmM0VEvSZKkk7CQKYOurrH694HfWvxyJEmSmtdC7lz/x8D91X8vATYAT9arIEmSpGa1kOB1XYxxx/RCCKEF+ET9SpIkSWpOJ7ydxMzQNeM5F9anHEmSpOa1kHO8pm8rAZVbSawAPl/PoiRJkprRQg41/pDnbicxBfTHGA/XryRJkqTmtJDg1QP8ZozxS/UuRpIkqZktZMqg84D/Xu9CJEmSmt1CgtfdwOqZK0IIN9SlGkmSpCa2kEONFwG7Qgj/DyhX170U+GjdqpIknZLy+ASDI2VWdHXQsaQt63Kk3MjLvrGQ4NUJ/O6M5XXAm070pBDCWcANwIYY4yXVdcuAm4G9VMLbn8cYd51s0ZKk55uYnGT7jj527urn0FCZlcs72NhbYuvm9bS1LuTghtSc8rZvLCR4/SLG+MPphRDC66jcwf5EXgPcDrxyxrr3A0/HGP8ihHAh8GXgtSdRrySphu07+rjzoT3Hlg8OlY8tv21Lb1ZlSZnL274xZ9QLIZwbQjgXOBhCOGfG8hPA+IleOMZ4GzA8a/VVwL3V9keBDSGE5adcvSSJ8vgEO3f112zbuesA5fGJxBVJ+ZDHfWO+Ea/pUa6VwBUz1h8GvnGK77eG54exoeq6ofme1NPTSXt7Y5+rUCp1Z12CZrFP8sc+OTX7DjzLoeFyzbaB4SO0LV1CafWZp/z69kv+2CcLU+99Y6aF9smcwSvG+GKAEMJ7Y4yfXZSqYD8ws7Ll1XXzGhgYXaS3z0ap1E1//+zBP2XJPskf++TUTYxPsLK7g4NDx/+C6elexsTY+Cn/bO2X/LFPFq6e+8ZMs/tkvhC2kLkaFyt0AdwBXApQPcfrkRjjvKNdkqT5dSxpY2NvqWbbxt7VXt2owsrjvrGQk+tPSQjhcuDtwNkhhI8CnwI+A9xcXV4PvKNe7y9JRbJ183qgct7KwPARerqXsbF39bH1UlHlbd9omZqaOvGjMtbfP5z/IufhsHD+2Cf5Y58sjsW+V5H9kj/2yamp5328ahxqbJnrsXUb8ZIkpdexpI01PZ1ZlyHlTl72De+qJ0mSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqREDF6SJEmJGLwkSZISMXhJkiQlYvCSJElKxOAlSZKUiMFLkiQpEYOXJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwkiRJSsTgJUmSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhIxeEmSJCVi8JKkJlIen2D/wCjl8YmsS5FUQ3vWBUiSTt/E5CTbd/Sxc1c/h4bKrFzewcbeEls3r6et1b+xpbwweElSE9i+o487H9pzbPngUPnY8tu29GZVlqRZ/DNIkhpceXyCnbv6a7bt3HXAw45Sjhi8JKnBDY6UOTRUrtk2MHyEwZHabZLSM3hJUoNb0dXByuUdNdt6upexoqt2m6T0DF6S1OA6lrSxsbdUs21j72o6lrQlrkjSXDy5XpKawNbN64HKOV0Dw0fo6V7Gxt7Vx9ZLygeDlyQ1gbbWVt62pZc3X34+gyNlVnR1ONIl5ZDBS5KaSMeSNtb0dGZdhqQ5eI6XJElSIgYvSZKkRAxekiRJiRi8JKmJOEm2VFte9g1PrpekJuAk2VJteds3DF6S1AScJFuqLW/7hn8GSVKDc5JsqbY87hsGL0lqcE6SLdWWx33D4CVJDc5JsqXa8rhvGLwkqcE5SbZUWx73DU+ul6Qm4CTZUm152zcMXpLUBJwkW6otb/uGwUuSmoiTZEu15WXf8BwvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwqsrLrOWpFO3znozh0TEee/IQw6NjWZcinTS3XynfMruqMYRwH3CkujgRY3xdFnXkbdbyeiva5z0ZY0ePcuOtD7O3f4TJKWhtgbWlLq6/+iKWtnsBsPLN7VdqDFn+pv1ujPGK6n+ZhC54btbyg0Nlpnhu1vLtO/qyKqmuivZ5T8aNtz7M7v2VX1oAk1Owe/8IN976cLaFSQvg9is1hiyD14UhhA+HELaFEK7KooA8zlpeT0X7vCdjeHSMvf0jNdv29o942Ea55vYrNY4sx59vijE+EEJoA34UQhiOMf6o1gN7ejppb1/8u8zuO/Ash4bnnrW8bekSSqvPXJT3KpW6F+V1TkfKz9sIZvbJLx7vPzZSMNvkFAyPTfKS87Lvw2aXh/2kEdV7+7Vf8sc+yZ+F9klmwSvG+ED1/xMhhH8ErgRqBq+BgdG61DAxPsHK7g4ODh0fRnq6lzExNk5///Bpv0+p1L0or3O6Un3eRjC7T7qXttLaQs1fXq0tlfai/Gyykpf9pBHVc/u1X/LHPsmf2X0yXwjL5FBjCOFlIYR3zFj1UuBnqevI46zl9VS0z3syujuXsrbUVbNtbamL7s6liSuSFs7tV2ocWZ3jNQRcFUL4WAjhL4HdwNezKGTr5vVsuXgdq5Yvo7UFVi1fxpaL12U2a3m9Fe3znozrr76Ic9Z00dpSWW5tgXPWVK4Kk/LO7VdqDC1TU3OcGJAj/f3DdS+yPD5Rt1nL8zgsXM/P2wjm65Ph0TH27B9h3RpHClLK437SiBZ7+7Vf8sc+yZ8ahxpb5nqsN3epysus5akU7fOejO7OpVzwopVZlyGdErdfKd+KfcdMSZKkhAxekiRJiRi8JEmSEjF4FZQT6c7NCcTVyNx+pXzz5PqCcSLduTmBuBqZ26/UGNwbC8aJdOfmBOJqZG6/UmMweBWIE+nOzQnE1cjcfqXGYfAqkD0zRrpmm5yqtBfV4EiZQzXmsITKBOKDI7XbpDxw+5Uah8GrQNbNmE5kttaWSntRrejqYOXyjpptPd3LWNFVu03KA7dfqXEYvArEiXTn5gTiamRuv1LjMHgVjBPpzs0JxNXI3H6lxuAk2QnkcULTok8EPV+fFH0C8azkcT9pRIu9/dov+WOf5I+TZOuEnEh3bk4grkbm9ivlm4caJUmSEjF4SZIkJWLwkiRJSsTgVVW0iWX79gzwle88Rt+egaxLkbSIivZdJjWawp9cX7SJZQ+NHOFDn/vxseUfPrIPgJuvu4yVXcuyKkvSaSrad5nUqAq/NxZtYtmZoWsh6yU1hqJ9l0mNqtDBq2gTy57osKKHHaXGVLTvMqmRFTp4FW1i2Xsefea02iXlU9G+y6RGVujgVbSJZTddeNZptUvKp6J9l0mNrNDBq2gTy65f13Na7ZLyqWjfZVIjK3TwguJNLHvzdZed1HpJjaFo32VSo3KS7Kp6ToycxwlN+/YMcM+jz7DpwrMKOdKVxz4pOvtkcThJdvOzT/LHSbJPQdEmll2/rqeQgUtqdkX7LpMaTeEPNUqSJKVi8JIkSUrE4CVJkpSIwUuSJCkRg1dVeXyC/QOjhZlaY9+BEb5731PsOzCSdSmSFtHw6BiPPXmI4dGxrEuRVEPhr2qcmJxk+44+du7q59BQmZXLO9jYW2Lr5vW0tTZfLh05MsYH/+oejk5U7tDxrbt/RntbC7e8ZxNdy5ZmXJ2kUzV29Cg33vowe/tHmJyC1hZYW+ri+qsvYml74b/qpdxovmRxkrbv6OPOh/ZwcKjMFHBwqMydD+1h+46+rEuri5mha9rRiSk++Ff3ZFSRpMVw460Ps3t/JXQBTE7B7v0j3Hjrw9kWJul5Ch28yuMT7NzVX7Nt564DTXfYcd+BkeNC17SjE1MedpQa1PDoGHv7a++/e/tHPOwo5Uihg9fgSJlDQ+WabQPDRxgcqd3WqB7pO3ha7ZLyac+Mka7ZJqcq7ZLyodDBa0VXByuXd9Rs6+lexoqu2m2NasP6VafVLimf1q3ponWOCUpaWyrtkvKh0MGrY0kbG3tLNds29q5e9Dkbs3b26i7a22p/O7e3tXD2ar+cpUbU3bmUtaXa++/aUhfdnV44I+VFoYMXwNbN69ly8TpWLV9GawusWr6MLRevY+vm9VmXVhe3vGfTceFr+qpGSY3r+qsv4pwZI1+tLXDOmspVjZLyo2Vqao4TA3Kkv3+47kWWxycYHCmzoqtj0Ue68jiT/L4DIzzSd5AN61cVcqQrj31SdPbJ4hgeHWPP/hHWrVmckS77JX/sk/yZ3SelUvccB/+9j9cxHUvaWNPTmXUZyZy9uquQgUtqdt2dS7ngRSuzLkPSHAp/qFGSJCkVg5ckSVIiBi9JkqREDF5VBwcP8+NH93Fw8HDWpSTx1DOD3HZXH089M5h1KblTtAnT1VzcfqV8K/zJ9YfHxvnwF+5l5PDRY+u6zmjnpmsv5YylSzKsrD4GR8t88LP3MH2Z6Lfvf5oW4Jb3bmJFZ3PdMPZkFW3CdDUXt1+pMRR+b5wdugBGDh/lw1+4N6OK6mtm6Jo2VV1fdEWbMF3Nxe1XagyFDl4HBw8fF7qmjRw+2nSHHZ96ZvC40DVtqtpeVEWbMF3Nxe1XahyFDl7x6V+dVnujefCx2l/MC21vZkWbMF3Nxe1XahyFDl7h3BecVnujueSC2vNSLrS9mRVtwnQ1F7dfqXEUOnitWnEGXWfUvr6g64x2Vq04I3FF9XXeWSuYaw6Dlmp7URVtwnQ1F7dfqXEUOngB3HTtpceFr+mrGpvRLe/ddFz4mr6qseiKNmG6movbr9QYnCS76uDgYeLTvyKc+4JFH+nK44SmTz0zyIOP9XPJBaVCjnTN1yf1nDBdc8vjftKIFnv7tV/yxz7JHyfJPgWrVpzBZRc216HF+Zx31opCBq6FKNqE6Woubr9SvhX+UKMkSVIqBi9JkqREDF6SJEmJGLwkSZISMXhV3ffTfXziKw9w30/3ZV1KErfdFbnulru57a6YdSm5Ux6fYP/A6AmnWRkeHeOxJw8xPDqWqDLpxPr2DPCV7zxG356BrEuRVEMmt5MIIWwB3gTsB6ZijB+f7/H1vJ3E0wdG2PalB45bv+2dr+bc1V2L8h55uvT30Z8f5L9+65Hj1n/gLRu48CWrMqgoG7X6ZGJyku07+ti5q59DQ2VWLu9gY2+JrZvX09b63N8oY0ePcuOtD7O3f4TJKWhtgbWlLq6/+iKWtnuh8KnK037SiA6NHOFDn/vxcetvvu4yVnYtO+XXtV/yxz7Jn5O5nUTyEa8QQifwReADMcZtwCtCCK9LXce0WqFrvvWNrlbomm99kWzf0cedD+3h4FCZKeDgUJk7H9rD9h19z3vcjbc+zO79ldAFMDkFu/ePcOOtD6cvWqqqFbrmWy8pG1kcarwUeCrGOD1r6z3AVRnUccLDis122PFEhxWLfNixPD7Bzl21JwnfuevAscOOw6Nj7O0fqfm4vf0jHnZUJk50WNHDjlJ+ZHFcZA0wc4x0qLpuTj09nbS3L/4dxO/655+coP0XvOHK3kV5r1Kpe1Fe53Tc/cgzJ2y/9i0XJ6omezP7ZN+BZzk0XK75uIHhI7QtXUJp9Zn84vH+YyNds01OwfDYJC85L/u+blR52E8a0fa7fzZv+0/6DnHpxnNP+fXtl/yxT/JnoX2SRfDaD8ysbnl13ZwGBkbrUsiVr/g1Hn96cN72xTiOnpfj8VdsOItv37933vY81JnC7D6ZGJ9gZXcHB4eOD1893cuYGBunv3+Y7qWttLZQM3y1tkD30tbC/AwXW172k0b0qvUr+d59T83bfqo/W/slf+yT/Klxjtecj83iUOO9wHkhhI7q8ibgjgzq4NdffvZptTeaf39lOK32ZtaxpI2NvaWabRt7Vx+b8667cylrS7Uvulhb6qK7c2ndapTmsn5dz2m1S0onefCKMY4C1wKfDSHcAPxzjPEHqeuYtu2drz6p9Y3uA2/ZcFLri2Tr5vVsuXgdq5Yvo7UFVi1fxpaL17F18/rnPe76qy/inDVdtFavWWltgXPWVK5qlLJy83WXndR6SdnI5HYSJ6uet5OYdt9P9/EPP9nNb77qnEUf6crjsPBtd0Xu3rmPKzaeXciRrvn6pDw+weBImRVdHcdGumoZHh1jz/4R1q1xpGsx5HE/aUR9ewa459Fn2HThWYsy0mW/5I99kj8nczsJg1cC7iT5Y5/kj32ST/ZL/tgn+ZPr+3hJkiQVlcFLkiQpEYOXJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJWLwkiRJSsTgJUmSlEhD3LlekiSpGTjiJUmSlIjBS5IkKRGDlyRJUiIGL0mSpEQMXpIkSYkYvCRJkhJpz7qAZhZCOAu4AdgQY7wk63qKLoRwPpX+eBhYBxyMMf5ZtlUphNAK/B1wP7AUOB/4jzHGw5kWJkIIZ1Dpl+/HGD+UdT1FF0K4DzhSXZyIMb4uy3pUEUIIwFuBw8DlwLYY4wNzPd7gVV+vAW4HXpl1IQJgJfDNGOPtACGEfwkh3BFj/EnGdQnujTHeABBCuB14E/C1bEsSlT9UdmZdhI75boxxW9ZF6DkhhDbgFuANMcbJEMKtwNH5nmPwqqMY420hhCuyrkMVMcYHZ61qBZ7NohY9J8Y4SeUXPCGEdiqjkTHTokQI4e3APcArgK6My1HFhSGEDwNnAA/GGO/IuiBxCdACvCeE0AkcBP56vid4jpcKKYTwRuB7McZ/zboWVYQQfhv4e+DvY4wPZV1PkYUQ/g1wQYzx/2Rdi57nphjjTcAngI+EEH4j64LEecClwN/EGD8J/AbwB/M9weClwgkhXAlcCXwg61r0nBjj92KM/xZ4cQjh3VnXU3BvBI6EEP6YyikTrw4hvD/jmgpv+ryhGOME8I9UvseUrSHgX2OMg9XlfwKumO8JHmpUoYQQrgJeC7wPODuEcF6M8d6Myyq06ujKi2ccNnkCeEmGJRVejPHG6X+HEJYBXTHGT2dYUuGFEF4GbIoxfrm66qXA/82wJFXcD6wKIbRVA/F5wK75nmDwqqMQwuXA26n8gv8o8Cmv1MpOCOFVwHbgIeAu4Ezg84DBK1tl4B0hhI3AEuAC4L3ZliSAEMKbqRw6WRpCeGuM8RtZ11RgQ8BVIYRfA5YDu4GvZ1uSYoyHqufdfTqE0A+UgHmvlm+ZmppKUpwkSVLReY6XJElSIgYvSZKkRAxekiRJiRi8JEmSEjF4SZIkJeLtJCQ1tRDCBuDTwFMxxmuqk3I/DlwcYxyY53nXAFfEGK+pLi/oeZI0H0e8JDW1GOMjwFdmLE8CW042PJ3q8yRpJu/jJSk3Qgh/CHwM+DaVG6u+HPjfwNHq+q9Ruav9FcA7gJ3Ax4F9wDnA12OM366+1ieBi4GfA5PAGdURr3cB24D/EGO8O4TwIuDPgaep3HU6Al8FvgCsBX4A/B1wbo3nfZLKjSzPBb4cY/yHEMK11Vq/VV3/cuBDMca/Xeyfl6TG44iXpNyIMX6RyuG8fTHG64DXA9dTmf/s+8CLYoxvBt4APAn8T+BLMcY/Aa4F/kcI4QUhhN8B/h3w2zHG/wx0z3iPv66+x7SvAd+IMf4RlZkmNsQYH6++9v0xxuuq80jWet43q8+7FvhqCOGFMcYvVGvtiTG+CbgO+JNF/DFJamAGL0l5dA9AjHGUyhRP05MB31ld/2PgZ1QmcL4mhPBF4C+r684BNgP/VD08CJUJhY8TQugGLpvxfmMxxt89UXE1njcA7KEyD+i0H1f//zPghSf8xJIKwZPrJTWSco1118cYfwkQQjgDGEtb0pyma53AP3IlVfllICmPLgUIIXRSOU/rrtkPiDEOUxlx+q3qY1uB7wAdwA7gtdV1UBkZO86M19hUfY1lIYTt1eYjQFsIoSWE8AcneF4PsI45RtYkaZon10vKlRDC3cCDwAogALcBT1A5lLgX+Oz0ieohhBcDn6JygvuZwG0xxu9W2z4JvAroA3qAVwJ/CryAyknyD1A5/2oplZPr9wKrgP8WY7wnhLCWyon9fVSCXNs8z1tH5Vyz71fPL5uu9Z3AR4DfAz4WY/zMIv+4JDUYg5ekXKkGr20xxrszLkWSFp2HGiXlRgjh3UAv8MEQwrqs65GkxeaIlyRJUiKOeEmSJCVi8JIkSUrE4CVJkpSIwUuSJCkRg5ckSVIiBi9JkqRE/j/kd7dW9Q1/+AAAAABJRU5ErkJggg==\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.drinks)\n",
    "ax.set_xlabel('prediction')\n",
    "ax.set_ylabel('truth')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2) Second example: fish caught in state park"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example is taken from here: https://stats.idre.ucla.edu/r/dae/zip/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('https://stats.idre.ucla.edu/stat/data/fish.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAFkCAYAAADi/3UkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASEElEQVR4nO3de4yl9V3H8c/C4sLszYEspbEVk9r87MVSBUwIt1bWqCEaQZNqEG2i0TSm6YW2tEorIK3UlmpIjOWPJv1DEk1phSZEKxJDKCG9bE3RJv01lhYUy866DLvAoctC1z/OWTNnO9vOMt8zt/N6Jfxxnpl9nt98meec9z7znNlNR44cCQAAy3fSai8AAGCjEFYAAEWEFQBAEWEFAFBEWAEAFBFWAABFNq/2ApJk376nJv47H2ZnZzI/P5j0YdYN8/h+ZjLOPMaZxzjzGGce4zb6PHbt2r7peB+bmitWmzefvNpLWFPM4/uZyTjzGGce48xjnHmMm+Z5TE1YAQBMmrACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKrIl/hHklHDhwII8/vnfix5mZmcmOHTsnfhwAYO2ZirA6ePBAzjv/nDw5/8TEj7Vz52z27HlIXAHAFJqKsBoMBnly/olcfNUt2bJ1dmLHOfTMfO6//ZoMBgNhBQBTaCrC6qgtW2dz6rbTV3sZAMAG5eZ1AIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAosrl6h621X03yU0lOSfKN3vunqo8BALAWLSmsWmtnJbkpyTm99/MXbN+d5Mokc0mO9N5vSLKn9/7Z1trOJJ9IIqwAgKmw1B8FXpTkriSbjm5orc0k+XiSd/Ter0/yutbaZb33x0afckWSjxauFQBgTVvSFave+x2ttTccs/mCJI/03g+NHj+Q5PIk97bWLk/ycJLHsgSzszPZvPnkpa34RTh8+KmJ7XsxZ5yxLbt2bV/RY74Y62GNK81MxpnHOPMYZx7jzGPctM5jOfdYnZlkYbEcTHJma+3Xklyb5KtJtie56oftaH5+sIxl/HD79z890f0vdrxTTlnZmDtRu3Ztz759a3uNK81MxpnHOPMYZx7jzGPcRp/HD4rG5YTVXIbhdNSOJHO99zuT3LmM/QIArEvL+XULDyY5u7W2ZfT4wiR3L39JAADr05LCqrV2aZKrk7y0tXZda+203vsgyVuS3NpauynJQ733eye4VgCANW2pN6/fl+S+Rbbfk+Se6kUBAKxHfvM6AEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQBFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAEWEFAFBEWAEAFBFWAABFhBUAQJHNk9hpa21zkmuTnN17/4NJHAMAYK2Z1BWrrUn+aYL7BwBYc5Z8xaq1dlaSm5Kc03s/f8H23UmuTDKX5Ejv/Ybe+4HW2v7y1QIArGEn8qPAi5LcleT1Rze01maSfDzJa3rvh1prn26tXdZ7v/dEFjE7O5PNm08+kT9yQg4ffmpi+17MGWdsy65d21f0mC/GeljjSjOTceYxzjzGmcc48xg3rfNYclj13u9orb3hmM0XJHmk935o9PiBJJcnOaGwmp8fnMinn7D9+5+e6P4XO94pp6xszJ2oXbu2Z9++tb3GlWYm48xjnHmMM49x5jFuo8/jB0Xjcu+BOjPJwskdTHJma21Tkjclaa21n13mMQAA1oXlvitwLsnCbNuRZK73fiTJh0f/AQBMheVesXowydmttS2jxxcmuXuZ+wQAWJeWHFattUuTXJ3kpa2161prp/XeB0nekuTW1tpNSR460RvXAQA2ihO5ef2+JPctsv2eJPdULgoAYD3yCzwBAIoIKwCAIsIKAKCIsAIAKLLc32PFIubm9k50/zMzM9mxY+dEjwEAnDhhVej55wbZdNLJ2b37kokeZ+fO2ezZ85C4AoA1RlgVeuHwoRz53gu5+KpbsmXr7ESOceiZ+dx/+zUZDAbCCgDWGGE1AVu2zubUbaev9jIAgBXm5nUAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACK+EeY16m5ub3L+vOHDz+V/fufPu7HZ2ZmsmPHzmUdAwCmjbBaZ55/bpBNJ52c3bsvmehxdu6czZ49D4krADgBwmqdeeHwoRz53gu5+KpbsmXr7ESOceiZ+dx/+zUZDAbCCgBOgLBap7Zsnc2p205f7WUAAAu4eR0AoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKLJ5tRfAdDt48EAGg8FEjzEzM5MdO3ZO9BgAkAgrVtHBgwdy7rmvy4ED8xM9zs6ds9mz5yFxBcDECStWzWAwyIED87n4qluyZevsRI5x6Jn53H/7NRkMBsIKgIkTVqy6LVtnc+q201d7GQCwbG5eBwAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAoIqwAAIoIKwCAIsIKAKCIsAIAKCKsAACKCCsAgCLCCgCgiLACACgirAAAimyu3mFrbSbJ9UkeTbK39/6p6mMAAKxFSwqr1tpZSW5Kck7v/fwF23cnuTLJXJIjvfcbRo+/1Hv/VGvtziTCCgCYCku9YnVRkruSvP7ohtGVqY8neU3v/VBr7dOttcuSvDzJg6NPO61ysaysubm963r/J3qsw4efyv79T7/oY8zMzGTHjp0v+s+vFQcPHshgMFj2PH6YjTIvYGmOPrdM2mo/tywprHrvd7TW3nDM5guSPNJ7PzR6/ECSy5N8Jcmu0bZnl7L/2dmZbN588lI+9UU5fPipie17I3r+uUE2nXRydu++ZLWXsmwr+bX86Ozp+fa3Hs7Ones3Fg4cOJDzzj8nT84/MfFjrcd57dq1fbWXsKaYxzjzGLdwHtP03LKce6zOTLKwWA6Otn0myfWttZckuX0pO5qfn2zBTvJv3RvRC4cP5cj3XsjFV92SLVtnJ3acg/u+nS/+w40T23+ycl/LoWfmc//t1+TRR/fmrLPW73tCHn98b56cf8K8FrFr1/bs2+cvaUeZxzjzGHfsPDbac8sPiujlhNVckoV73pFkrvc+SPKeZeyXNWLL1tmcuu30ie3/0DPzE9v3sSb9tWw05gVMwjQ8tywn5x5McnZrbcvo8YVJ7l7+kgAA1qclhVVr7dIkVyd5aWvtutbaaaMrU29Jcmtr7aYkD/Xe753gWgEA1rSl3rx+X5L7Ftl+T5J7qhcFALAerY+7RgEA1gFhBQBQRFgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEWEFQBAkU1HjhxZ7TUAAGwIrlgBABQRVgAARYQVAEARYQUAUERYAQAUEVYAAEU2r/YCVkJrbXeSK5PMJTnSe79hlZe0olprr0hyU5KvJHlZkv299xtba6cnuTnJw0lemeSPe+97V2+lK6e1dlqSLyT55977u1prpyb5aJLHMpzFzb33b6zmGldSa60l+a0kzya5NMn1GZ4v70/yn0l+Isk1vfenV2mJK6q19u4Mv+b/zfD74feSnJYpOV9aa2dl+JxxTu/9/NG2454jrbXfTvIzSV5I8s3e+22rsvAJOs5Mrk1yVpLvJDkvyQd6718ffWxDz2SxeSz42FVJ/jbJ9qPPGdP0Orzhr1i11maSfDzJO3rv1yd5XWvtstVd1Yo7Pcnf9d4/0nt/W5LfbK2dm+RDSf6l935zkjszfNKcFjcl+bcFj9+e5NHe+58n+cskn1iVVa2C1trJST6W5Mbe+4czjIhvZXje3DaayX8kuXb1VrlyRi8Y70vy1t77nybZmuELwjSdLxcluSvJpgXbFj1HWmsvS/KuJO/qvb8nye+31l65wutdCYvNZFuSd/be/yLJp5N8JJmamSw2j7TWXpXk1cdsm6rX4Q0fVkkuSPJI7/3Q6PEDSS5fxfWsuN77l3rvdy3YdFKSZzKcw4OjbVMzl9ba1Rl+vd9asPn/Z9F7//ck57TWdqzC8lbD+Rk+Ob61tfa+JL+S5Mkkb0zypdHnTM33R5JBkueSHP3/vy3J1zJF50vv/Y4kTx2z+XjnyC8m2dN7P/rbph9M8ssrtdaVsthMeu/vX/B1n5Tk6BXdDT+TxeYxCqj3JDn2atRUvQ5PQ1idmfH/+QdH26ZSa+2KJJ8bXa5eOJuDSWZbaxv6x8OttVcneVXv/TPHfGiav0/OzvCJ75OjqxGXZPi37WcXvDBMzTx67weTvDvJ37fWPpnkvzP8cejUnS/HON45Ms3nTpKktfYjSX43yXWjTdM6kw9meOX7uWO2T9U8piGs5pJsX/B4x2jb1GmtvTHDqxDvGG1aOJsdSeZ778+vxtpW0BVJvttae2+Gl7J/rrX29kz398nBJF/vvR8YPf58ktcmOa21dvQy/9TMo7X2+gzD6vLe+5szvM/qA5nO82Wh450j03zuHI2qv0nyJ733b442T91MWmsvTzKb5E2j59ckeWdr7bxM2TymIaweTHJ2a23L6PGFSe5exfWsitba5Rlenn5bkrNaaxdkOIcLRp8yFXPpvX+w937j6D6Zzyf5Yu/9r7JgFq21n07y1dGVi2nwhSRnjO61SoZXsL6W5F8z/DFhMiXfHyM/luSJBdH0nSSnZgrPl2Mc7xz5XJJzF0T4BUn+cXWWuLJGP/q6LcnHeu97Wmu/PvrQ1M2k9/5fvfc3995vHj2/JsO5fDlT9jo8Ff8Ic2vtF5L8RpJ9SQ5v5HcjLGZ0o/p9Sb482rQ1yV8n+WySDyd5JMkrkrx3o77L6VijJ8A/SvIjGc7i6M3I30nyk0k+NGXvCrwiyc9neI78eJK3JnlJhldqHh5te+c0vCtwFJi3JvluhveavTbDG7cPZUrOl9bapUl+J8kvZXg15pbRhxY9R0bvgDsvw3fAfWOjvQMuOe5Mbs/w++N/Rp+2dcE7Bjf0TBabR+/92dbariR/mOTPRv/d1nt/bJpeh6cirAAAVsI0/CgQAGBFCCsAgCLCCgCgiLACACgirAAAiggrAIAiwgoAoIiwAgAo8n/9LE85NyiddgAAAABJRU5ErkJggg==\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",
    "\n",
    "ax.hist(df['count'], bins=20, edgecolor=\"k\", linewidth=1.2)\n",
    "ax.set_yscale('log')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 4.126434\n",
      "         Iterations: 16\n",
      "         Function evaluations: 18\n",
      "         Gradient evaluations: 18\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 = ZIP(endog = endog_var,\n",
    "        exog = exog_vars,\n",
    "        exog_infl=exog_vars_infl,\n",
    "        inflation='logit').fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ZeroInflatedPoisson 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>ZeroInflatedPoisson</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.08466</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>00:26:40</td>       <th>  Log-Likelihood:    </th> <td> -1031.6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>              <td>True</td>         <th>  LL-Null:           </th> <td> -1127.0</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>      <th>  LLR p-value:       </th> <td>3.647e-42</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.2974</td> <td>    0.365</td> <td>    3.554</td> <td> 0.000</td> <td>    0.582</td> <td>    2.013</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_x1</th>    <td>   -0.5643</td> <td>    0.151</td> <td>   -3.735</td> <td> 0.000</td> <td>   -0.860</td> <td>   -0.268</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>         <td>    1.5979</td> <td>    0.085</td> <td>   18.797</td> <td> 0.000</td> <td>    1.431</td> <td>    1.765</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>            <td>   -1.0428</td> <td>    0.089</td> <td>  -11.660</td> <td> 0.000</td> <td>   -1.218</td> <td>   -0.868</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>            <td>    0.8340</td> <td>    0.093</td> <td>    9.000</td> <td> 0.000</td> <td>    0.652</td> <td>    1.016</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                     ZeroInflatedPoisson Regression Results                    \n",
       "===============================================================================\n",
       "Dep. Variable:                       y   No. Observations:                  250\n",
       "Model:             ZeroInflatedPoisson   Df Residuals:                      247\n",
       "Method:                            MLE   Df Model:                            2\n",
       "Date:                 Fri, 17 Apr 2020   Pseudo R-squ.:                 0.08466\n",
       "Time:                         00:26:40   Log-Likelihood:                -1031.6\n",
       "converged:                        True   LL-Null:                       -1127.0\n",
       "Covariance Type:             nonrobust   LLR p-value:                 3.647e-42\n",
       "=================================================================================\n",
       "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "inflate_const     1.2974      0.365      3.554      0.000       0.582       2.013\n",
       "inflate_x1       -0.5643      0.151     -3.735      0.000      -0.860      -0.268\n",
       "const             1.5979      0.085     18.797      0.000       1.431       1.765\n",
       "x1               -1.0428      0.089    -11.660      0.000      -1.218      -0.868\n",
       "x2                0.8340      0.093      9.000      0.000       0.652       1.016\n",
       "=================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative: run the model with robust standard errors -- these are MUCH larger for the parameters of the count component of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 4.126434\n",
      "         Iterations: 16\n",
      "         Function evaluations: 18\n",
      "         Gradient evaluations: 18\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 = ZIP(endog = endog_var,\n",
    "        exog = exog_vars,\n",
    "        exog_infl=exog_vars_infl,\n",
    "        inflation='logit').fit(cov_type='HC0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ZeroInflatedPoisson 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>ZeroInflatedPoisson</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.08466</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>00:26:40</td>       <th>  Log-Likelihood:    </th> <td> -1031.6</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>converged:</th>              <td>True</td>         <th>  LL-Null:           </th> <td> -1127.0</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC0</td>         <th>  LLR p-value:       </th> <td>3.647e-42</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.2974</td> <td>    0.363</td> <td>    3.570</td> <td> 0.000</td> <td>    0.585</td> <td>    2.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>inflate_x1</th>    <td>   -0.5643</td> <td>    0.147</td> <td>   -3.843</td> <td> 0.000</td> <td>   -0.852</td> <td>   -0.277</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>         <td>    1.5979</td> <td>    0.283</td> <td>    5.647</td> <td> 0.000</td> <td>    1.043</td> <td>    2.152</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>            <td>   -1.0428</td> <td>    0.309</td> <td>   -3.374</td> <td> 0.001</td> <td>   -1.649</td> <td>   -0.437</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th>            <td>    0.8340</td> <td>    0.388</td> <td>    2.148</td> <td> 0.032</td> <td>    0.073</td> <td>    1.595</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                     ZeroInflatedPoisson Regression Results                    \n",
       "===============================================================================\n",
       "Dep. Variable:                       y   No. Observations:                  250\n",
       "Model:             ZeroInflatedPoisson   Df Residuals:                      247\n",
       "Method:                            MLE   Df Model:                            2\n",
       "Date:                 Fri, 17 Apr 2020   Pseudo R-squ.:                 0.08466\n",
       "Time:                         00:26:40   Log-Likelihood:                -1031.6\n",
       "converged:                        True   LL-Null:                       -1127.0\n",
       "Covariance Type:                   HC0   LLR p-value:                 3.647e-42\n",
       "=================================================================================\n",
       "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "inflate_const     1.2974      0.363      3.570      0.000       0.585       2.010\n",
       "inflate_x1       -0.5643      0.147     -3.843      0.000      -0.852      -0.277\n",
       "const             1.5979      0.283      5.647      0.000       1.043       2.152\n",
       "x1               -1.0428      0.309     -3.374      0.001      -1.649      -0.437\n",
       "x2                0.8340      0.388      2.148      0.032       0.073       1.595\n",
       "=================================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = m.predict(exog=exog_vars, exog_infl=exog_vars_infl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAFyCAYAAABbdsanAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5ykVX3n8U9fZnqm6e6hZ6ZHLsNFGfqIgjBcjCPqIMyGGMKuSnaJrBpeGJM10WgUMQq6aCDieolxN4mu0WRZYzL7wo2aEC/BUVYRFGFEMXKGQUFmQLtnppnupmeqr/tHVY9NU119redU9fN5v16+nDpPPVW/U011f5/znOc8DRMTE0iSJCmdxtQFSJIk5Z2BTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhJrTl3AYvX2Dixq3Y7Ozlb6+oaWqpy6kdd+g33PY9/z2m+w7/Y9X+qh311d7Q3l2nM/Qtbc3JS6hCTy2m+w73mU136Dfc+rvPa9nvud+0AmSZKUmoFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmScq0wMkZP3xCFkbFkNdT9rZMkSZIWYmx8nO07drNzVy8H+gus7Whhc3cXl1+4iabGbMesDGSSJCmXtu/YzW3f23Pk8f7+wpHHV2zrzrQWT1lKkqTcKYyMsXNXb9ltO3fty/z0pYFMkiTlzsHBAgf6C2W39Q0c5uBg+W3VYiCTJEm5s6athbUdLWW3dbavYk1b+W3VYiCTJEm507Kiic3dXWW3be5eT8uKpkzrcVK/JEnKpcsv3AQU54z1DRyms30Vm7vXH2nPkoFMkiTlUlNjI1ds6+ayradwcLDAmraWzEfGJhnIJElSrrWsaGJDZ2vSGpxDJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpRY1VbqDyEcA9wAnBljPG/atv8MfAZojzEOltq2Aa8EeoCJGON7q1WbJElSLanmCNmLgC8ADVMbQwinAc+Z1tYKfBz4oxjj9cDzQggXVbE2SZKkmlG1QBZjvAUYmNpWCl7XANNHv7YAj8QYC6XHdwCXVKs2SZKkWpL1zcVvBN4XYxwOIUxt38BTw1t/qW1WnZ2tNDcv7s7sXV3ti9q/XuW132Df8yiv/Qb7nld57Xu99juzQBZCOAHoBC6fEsbeGkL4F4rzxqZ+gh2ltln19Q0tqq6urnZ6ewdmf+Iyk9d+g33PY9/z2m+w7/Y9X+qh3zMFxswCWYzxUeDKycchhPcDH4kxDpZOZZ4UQmgpnbY8H/jLrGqTJElKqWpzyEIIW4HXAMeGEK4LIawutXeFEK4rPe2aEMLxMcYh4A3Ax0IINwA/iDF+rVq1SZIk1ZKqjZDFGG8Hbi/T3ktxOYwbprX/K/Cv1apHkiSpVrkwrCRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCmx5mq+eAjhGOAG4MwY43mltncAxwCPA+cC74kxPlDa9mpgMzAGPBRj/EQ165MkSaoF1R4hexHwBaBhSlsb8NYY438DPgd8ECCEsBG4Grg6xngN8DshhFOrXJ8kSVJyVQ1kMcZbgIFpbe+OMU5Mef/B0r8vBu6Zsu1O4GXVrE+SJKkWVPWUZSUhhJXAbwN/UGrawFPDW3+praLOzlaam5sWVUtXV/ui9q9Xee032Pc8ymu/wb7nVV77Xq/9ThLISmHsr4BrY4wPlZp7gE1TntYB7J7ttfr6hhZVS1dXO729A7M/cZnJa7/Bvuex73ntN9h3+54v9dDvmQJj5ldZhhBagU8AH4kx3hNCuKy06SvAOSGEyflmW4AvZV2fJElS1qp9leVW4DXAsSGE64APA58BTgeeGUIAOAr4XIxxTwjhQ8CfhRDGgL+OMT5YzfokSZJqQVUDWYzxduD2ac2vrPD8z1AMbJIkSbnhwrCSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxJqr+eIhhGOAG4AzY4znldpWAR8C9gKnAjfFGHeVtr0a2AyMAQ/FGD9RzfokSZJqQbVHyF4EfAFomNL2FuBnMcb3A38GfAoghLARuBq4OsZ4DfA7IYRTq1yfJElSclUNZDHGW4CBac2XAHeWtv8QODOE0AFcDNwTY5woPe9O4GXVrE+SJKkWVPWU5Qw28NSQ1l9qm6m9os7OVpqbmxZVUFdX+6L2r1d57TfY9zzKa7/BvudVXvter/1OEch6gKmfVkeprQfYNK1992wv1tc3tKhiurra6e2dPoi3/OW132Df89j3vPYb7Lt9z5d66PdMgTHFVZa3AlsAQghnAPfFGPuBrwDnhBAm55ttAb6UoD5JkqRMVTWQhRC2Aq8Bjg0hXBdCWA38OXBSCOE64G3A6wBijHsoXn35ZyGEDwN/HWN8sJr1SZIk1YKqnrKMMd4O3F5m0x/M8PzPAJ+pZk2SJEm1xoVhJUmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYk1p3jTEMLbgZOBfcCpwOuA1cBNwE9Kbe+KMf4iRX2SJElZynyELIRwDPBO4E0xxv8KHAW8EvhT4LYY403A54EPZV2bJElSCgsKZCGE/7qI9xwChoGO0uM24EfAJcCdpbY7So8lSZKWvYaJiYmKTwgh/C5wPbABaCj9byLG2LTQNw0hvAZ4NfB46fX+ANgPPCPG+EQIoRkYAVbEGEcrvdbo6NhEc/OCS5EkScpSQ7nGucwhewuwFXgoxjgOEEJ450KrCCGcBbwdODvGOBpC+DDwHqAHaAeeoDh61jdbGAPo6xtaaCkAdHW109s7sKjXqEd57TfY9zz2Pa/9Bvtu3/OlHvrd1dVetn0ugey+GOOD09q+tIhajgcOTAlbjwMnArcCW4BHgfNLjyVJkpa9GQNZCOE9pX8+GUL4GnAXUCi1/TrwggW+55eBXy+NjD0BnE5xFK4AfCCE0A2cAly9wNeXJEmqK5VGyC4F/pniiNWjpbaGaf8/bzHGMYpzxsp5/UJfV5IkqV5VCmTXxBi/Xm5DCOErVapHkiQpd2Zc9mIyjIUQfn9qewjhSuBV1S1LkiQpP+ayDtmzpz6IMf4t0FqVaiRJknKo0qT+nwITwNoQwm9M2dQE3F/twiRJkvKi0hyyCyhO3r8emLoy/2HvMSlJkrR0ZgxkMcZHSv+8cvq2EMLpMUZHySRJkpbArAvDhhBeW6b51cCvLn05kiRJ+TOXlfr/GPhO6d8rgDOBh6tVkCRJUt7MJZC9Mca4Y/JBCKEB+JPqlSRJkpQvsy57MTWMTdnnjOqUI0mSlD9zmUM2ufwFFJe8WAP8RTWLkiRJypO5nLK8nV8uezEB9MYYD1WvJEmSpHyZSyDrBP5djPGvq12MJElSHs3l1kknAZ+udiGSJEl5NZdA9g1g/dSGEMINValGkiQph+ZyyvJsYFcI4UdAodR2KnBd1aqSJEnKkbkEslbgP0x5vBF4ZXXKkSRJyp+5BLLHYoy3Tz4IIVxEccV+SZIkLYEZA1kI4cTSP/eHEE4AGkqPfwqMVLswSZKkvKg0QjY5KrYWuGBK+yHg76tVkCRJUt7MGMhijM8ECCH8YYzxY9mVJEmSlC9zuZelYUySJC1bhZExevqGKIyMJathLpP6JUmSlp2x8XG279jNzl29HOgvsLajhc3dXVx+4SaaGueyVOvSMZBJkqRc2r5jN7d9b8+Rx/v7C0ceX7GtO9Naso1/kiRJNaAwMsbOXb1lt+3ctS/z05cGMkmSlDsHBwsc6C+U3dY3cJiDg+W3VYuBTJIk5c6athbWdrSU3dbZvoo1beW3VYuBTJIk5U7LiiY2d3eV3ba5ez0tK5oyrcdJ/ZIkKZcuv3ATUJwz1jdwmM72VWzuXn+kPUsGMkmSlEtNjY1csa2by7aewsHBAmvaWjIfGZtkIJMkSbnWsqKJDZ2tSWtwDpkkSVJiyUbIQggBeBXFm5VvBa4HeoB3A7uBk4G3xRgHE5UoSZKUiSQjZCGEJuAjwPtijB8AXgf8FPg48IkY4/uB+4F3pKhPkiQpS6lOWZ4HNABvCiG8E7gUeAJ4KXB36Tl3AJekKU+SJCk7qU5ZngRsAV4VYzwYQvgMsA44FGOcKD2nH9gw2wt1drbS3Ly4KyK6utoXtX+9ymu/wb7nUV77DfY9r/La93rtd6pA1g88EGM8WHr8LeDFwOoQQkMplHVQnFNWUV/f0KIK6epqp7d3YFGvUY/y2m+w73nse177DfbdvudLPfR7psCY6pTld4B1pblkUBwx+xHwdYqnMwHOB25NUJskSVKmkoyQxRgPhBDeAXw0hNALdAHvAz4LvCeE8KvAicBbU9QnSZKUpWTLXsQY/xH4x2nNDwNXZV+NJElSOi4MK0mSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJEmJGcgkSZISM5BJkiQlZiCTJElKzEAmSZKUmIFMkiQpMQOZJElSYgYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCXWnOqNQwirge8AX40xXh1CWAV8CNgLnArcFGPclao+SZKkrKQcIbsB2Dnl8VuAn8UY3w/8GfCpJFVJkiRlLEkgCyG8BrgD+OmU5kuAOwFijD8EzgwhdCQoT5IkKVOZn7IMITwHOC3G+K4QwvOmbNoADEx53F9q66/0ep2drTQ3Ny2qpq6u9kXtX6/y2m+w73mU136Dfc+rvPZ9If0+PDxKX3+Bzo4WVq1MM5srxbu+AjgcQvhj4EXAyhDCW4AeYOqn2FFqq6ivb2hRxXR1tdPbOzD7E5eZvPYb7Hse+57XfoN9t+/5Mt9+j42Ps33Hbnbu6uVAf4G1HS1s7u7i8gs30dRYnZOIMwXGzANZjPHGyX+XJvK3xRg/Wvr3FuCbIYQzgPtijBVHxyRJkhZq+47d3Pa9PUce7+8vHHl8xbbuTGtJNqk/hHAZ8BLgBSGEVwF/DpwUQrgOeBvwulS1SZKk5a0wMsbOXb1lt+3ctY/CyFim9SRb9iLG+Dngc9Oa/yBFLZIkKV8ODhY40F8ou61v4DAHBwts6GzNrB4XhpUkSbmzpq2FtR0tZbd1tq9iTVv5bdViIJMkSbnTsqKJzd1dZbdt7l5Py4rFreAwX8lOWUqSJKV0+YWbgOKcsb6Bw3S2r2Jz9/oj7VkykEmSpFxqamzkim3dXLb1FA4OFljT1pL5yNgkA5kkScq1lhVNmU7gL8c5ZJIkSYkZyCRJkhIzkEmSpFwbGBrmxw8fYGBoOFkNziGTJEm5NDw6yo0338ve3kHGJ6CxAY7vauPa157NyuZsI5IjZDlXGBmjp28o81tESJKU2o0338ujPcUwBjA+AY/2DHLjzfdmXosjZDk1NjbOZ2/blekd7iVJqhUDQ8Ps7R0su21v7yADQ8O0t67MrB7/8ubUp//pR9z2vT3s7y8wwS/vcL99x+7UpUmSSjyLUT17poyMTTc+UdyeJUfIcqgwMsZd9z9edtvOXfu4bOspyRbGkyTB2Pg423fs9ixGFW3c0EZjA2VDWWNDcXuW/Knm0MHBAr1PHCq7bfIO91K1eeQvzWz7jt2exaiy9taVHN9VPnQd39WW6elKcIQsl9a0tdB19Gp6+p4eylLc4V754pG/VFlhZIydu3rLbvMsxtK69rVnz3iVZdYMZDnUsqKJF5x+LF/85k+eti3FHe6VL5NH/pMmj/wBrtjWnaosqWYcHCxwoL/8mYrJsxipb/OzXKxsbua9Vz2fgaFh9vQMsnFD9iNjkzwczamrLn0u287dyLqOVTQ2wLqOVWw7d2OSO9wrP2Y78vf0pVQ8i7G2o/yZCs9iVEd760pOO3ltsjAGjpDlVlNT7dzhXvnhkb80u5YVTWzu7nrKSPIkz2IsXwaynKuFO9wrPyaP/PeXCWUe+Uu/NHm2YueuffQNHKazfRWbu9d7FmMZM5BJyoxH/tLcNDV6FiNvDGSSMuWRvzR3nsXIDwOZpEx55C9JT2cgk5SER/6S9EsueyFJkpSYgUySJCkxA5kkSVJiBjJJkqTEDGSSJCnXCiNj9PQNJb19m1dZSpKkXBobH2f7jt3s3NXLgf4Cazta2NzdxeUXbqKpMdsxKwOZJOmIgaFh9vQMsnFDW9IbLUtZ2L5j91PuHLK/v3Dk8RXbujOtxUAmSWJ4dJQbb76Xvb2DjE9AYwMc39XGta89m5XN/qnQ8lMYGWPnrt6y23bu2sdlW0/JdNFq55BJkrjx5nt5tKcYxgDGJ+DRnkFuvPnetIVJVXJwsMCB/kLZbX0Dhzk4WH5btSQ57AkhnALcANwLbAT2xxjfF0JYC9wE/AQ4FXhXjPEXKWqUpLwYGBpmb+9g2W17ewcZGBr29KWWnTVtLaztaGF/mVDW2b6KNW0tmdaTaoRsLfAPMcYPxhjfDPxWCOEc4E+B22KMNwGfBz6UqD5Jyo09U0bGphufKG6XlpuWFU1s7u4qu21z9/rM77GbZIQsxnj3tKZG4EngEuDGUtsdwP/Ksi5JyqONG9pobKBsKGtsKG6XlqPLL9wEFOeM9Q0cprN9FZu71x9pz1LDxMQMh0UZCSG8ArggxvjmEEIBeEaM8YkQQjMwAqyIMY7OtP/o6NhEc3O2KVaSlps3f/jr/OSx/qe1P+u4Dv78bS9NUJGUncPDo/T1F+jsaGHVyqqPVTWUa0x66UwI4aXAS4G3lJp6gHbgCaAD6KsUxgD6+oYWVUNXVzu9vQOLeo16lNd+g33PY9/z2m+Ye9+vueKssldZXnPFWXX72flzz1/fF9PvZmDg4CGq/al1dbXP+P5JhBAuAV4MvBk4NoRwEnArsAV4FDi/9FiSVGUrm5t571XPdx0yKZFUV1meA2wHvgd8HTgK+AvgXcAHQgjdwCnA1Snqk6S8am9dyWknr01dhpQ7qSb13wPMNEv09VnWIi0nhZExDg4WWNPWkvkVQpKkhXP5ZWkZqKX7sUmS5s9AJi0DM92PbWx8govPO8ERM6lOOeqdHwYyqc5Vuh/b7Tv38vV797KuNGL2xv+0OePqJC2Eo975YyCT6lyl+7FNLvQ5OWLWunolLz//5OyKk7QgM416A1yxrTtVWaoiY7ZU5ybvxzYXd93/OIWRsSpXJGkxKo1679y1z+/wMmUgk+pcy4omWletmNNz9z1xiIOD5UfTJNWGSqPefQOH/Q5XQWFkjJ6+oaRh11OWUp0rjIzx5KHhOT13/dGrWdM2t9E0SWlMjnrvLxPKOttX+R1eQrU0V88RMqnOHRws0Dcwt0D2gtOP9Uotqca1rGhic3dX2W2bu9f7HV5Ck3P19vcXmOCXc/W279ideS0GMqnOVZpD1tgADQ2wrmMV287dyFWXPjfj6iQtxOUXbmLbuRtZ17GKxinf4csv3JS6tGWj1ubqecpSqnOTR9NTr8iatPWs47j4+SceWcOoqcljMKkeNDU2csW2bi7beorrkFXJXObqbehszaweA5m0DEweNe/ctY++gcN0tq9ic/d61yyS6lzLiqZMQ0GerGlrobN9JQfKTPk4uq0l87l6BjJpGfBoWlqeXKm/elpWNHHU6vKB7KjVKzL/vA1k0jLi0bS0PNTS1X/LVWFkjKHDI2W3DR0eoTAylmko86cqzVMtrFczV/VUq6RfqqWr/5arynPICpmv9+YImTRH9XTEOlOt3stSqn2zXf132dZTPH25BGptvbfa+isi1bB6OGKdHBH77L/uKlvrp//pR6lLlDQLV+rPRq2t9+YImTQHtX7EOn1ErKGh/PPuuv9xXvb8Ezy6lmpYrY3cLGe/ecGziD97gr29g4xPFNduPL6rjd+84FmZ1+IImTQHtX7EOn30bnyi/PPq5V6Wzn1TntXayM1ydss3fsKjPYNHfmeOT8CjPYPc8o2fZF6LI2TSHNTyEWul0bvpav1elvU0T0+qploauVmuau3Mh7/hpDmo5SPWg4OFskGxnFq6l2W5UbB6mKcnZaGWRm6Wq1o78+EImTRHlVbDT2l1SzONDTOfpoTiffA2d6/nqkufy4EDT2ZXXBkzjYK9/MXPqqmjVSmVwsgY9zzwi7Lb7nmgx+/CEnGlftWM5bQCdBZ9qdXV8A8VRiuGsd+55DTOOGUd7a0ra+JelpOjYJMmR8EOHR6tqfvKSakcHCzQN1h+wdK+wWG/C0vElfqV3Nj4OJ/8/A+54769dT9PJ8Wco1pbDX9NWwudbSvK/gJvAD51649rZh2ySkf+//bw/hmPVlPP05Oy1NQ4w2XSc9yuuXGlfiW3fcduvvjNnyyLeTq1OOco6ysEW1Y0MTJWfohsovS/WlmHrPKR/wjPPmlt2W2p5+nliVe4ptfTd2hR2zU3lebf7u93pX5VWa1dVbIYhZEx7o09ZbfdG3sz70uqKwQHhoYZOjw6p+emXodstiP7l/3KRn7w0D4GD/2yP22rm3nFS55Z7dJyzytca8fRbSsXtV1zs7qlmQaKB63TNZS2Z8lvWc4s1VUl1TiKHhga5scPH2Bg6OmnrMo5OFgoe3oL4EAV70M2U99TjdbtmXIl1mxSr0M225H9x265/ylhDGDw0Cg3fWZnNcsStTnanFdPDFb+HTjbds3NocJo2TAGxZB2qDC3A92l4ghZzix2Pa1qHEUPj45y4833Pm29nWtfezYrm2f+T7TS1YWNDUt/dFOp76NjExVGHqs7WrdxQ9ucn5t6HbKNG9oq/sx6Dx4uu9+enkEGhoZpb3VkoBoqj5xnP9qcd7N9p+fzndfMsv4bMhtHyJaB+YxWLXY9rWocRd94871l19u58eZ7K+5X6erC8YmlP7qp1PeUcxHaW1cy1ym+qdcha29dyfFd5f+YrFuzasb9JiiGMlVHrc2lybv21pW0rS4fBtpWN3tgskSy/hsyGwNZHRsbH+ezt+3iuk/exTs/cRfXffIuPnvbLsbGxyvud/mFm/j3L34W6zpW0dhQXKNq27kbZ11Pa7aj6IWcvhwYGmZvb/k/tHt7Byuevpy8urCczrYVSzoSVHldoF/Q1NjATNOjqn2ktf/goRmH3adqW93Mqy8OVatjrt5+xVk0Nz31w2puauD1lz6n4n7Om6meybk05aSYS5N3hZExmmf469zciBdcLJE1bS2s6yj/d2JdR/brkBnI6thMIzZ/8y8PVPzCNjU28vqXn8F7rjyXP/zN5/GGlz+Hy7aeMuspx2ocRVea/zQ+UXlUpNLVhSNjE0s6EjTb1YE9fYeSHWnFnz0xp+cNHhrlj//yjqrVMVcf/Oz3GZ32cxsdm+CTX/y3ivs5b6Z6am0uTd4dHCzwxJPlP/Mnnhx1xHKJVD5j1OU6ZJqbSiM2377/58Sf9c04t2tsfJxP/OMPuO27P+PwcDG4rVrZyAvPOJZXXXTqjMGsGufbZ5tTVGmuRKWrC4cOjy7pnKNDhfJhbNKK5oYZ1wJb6tG66cKJR8/5uQ//vD/pXKxKI6IzzR+b5AhZ9cz23/ds27W0XIcsO7V0z1BHyGbx+L5BvnzXIzy+b+7zV3bv6eN/fO77fOLzP+SRnx8E4Fv37eU9n7qLb923d0nqqjRiA5Xndm3fsZt//tZPj4QxgMPD4+y4Z2/FuWDVON/e3rqy4mtWCg6LGV2brx8/XHkUavee/sxG66Zbt2b1004BzmR8PO1crPlcETqdI2TVM9t/37Nt19JyHbLs1NI9Q2tuhCyEsA14JdADTMQY35uijsHDw7z1v99x5NTK//nGQzQ3NfCRN51P26ryIeHA4GGu/h/ffkrbdx546pyrT38p8ukvRd591Xk8c0P7guub6xHSXffvfcoVUoWRMe6J5eeBAXz3R4/NeEXVbCNgCxkhmy3oPr5vkGPXlx8la11VOeTMtn0+Nm3sqLj9uPWtPHmofCB98tDSjtZNVxgZo311c8WAPtWGztVVqWMuKo2IzmamFbW1eMetr3znidm2a2kt5syB5q7W1uWsqRGyEEIr8HHgj2KM1wPPCyFclKKWqWFs0ujYBG/97zPPwZkexir5k0/fveDaYO5HSIOHJ54y3+DgYIG+gZnnH/QfGp9xfsJsI2ALGSG7b/f+BW/f2ztUcd/Zts/HyGjlBPHYvqGKc3CqOSo122jpdGMLHaJaApWusly1svIvvrsfKL8IsBbvsX2VvyuzbdfSam9dSeMMB92NjQ1eZblElmpdzqVSU4EM2AI8EmOc/BTuAC7JuojH9w0+LYxNGh2bKDuqs3tP37zfZzGnL/cfnPuQ9b/99JehZmys8hWYlV775/ufrLjfbNvL6e2rHFQqbX/4scqf+Wzb52O2z3u2fiy8p4wAAAuTSURBVFRzdGdybt9cpZ5/cu1rz+aE0ggAFI/4T9jQxrOecVTF/YYOecqyWlwZvrYMDA0zNsPfoLGxiTkvnq3KJtflLCfF/XNr7ZTlBmBgyuP+UtuMOjtbaW5e3JBiV9dTTx1+8/7yk+Un7f75kzzvtGOf0rb9Gw/N+32/9v29vGLbs+e9H8B3ZxhmLef2Hz7Of7y4uKTADx+ZfS7IDx5+gpecd/LT2u/58gMV97tn934u2jK/iZB3P1h5hOzuB/fz1teWP7V7V9xXcd+74j7e/Ory+07/mc9mts97tn784OE+fu3FlZcVWajRfU/O6xRgYaJh3v1fan/5jos4OFjg4cf7OfnYDta0tfCKa75QcZ8H9hxcVN2p+5zSbH1/ZF/l31+P7Bvi0jr9/Orx5/7Yg70VR9wHhsd51kmz96se+74U5tPv8888ni9+8+nzxc4/8zg2Hjf3C6aWQq0Fsh5g6ifZUWqbUV/f4obSu7ra6e0deErbpmMqH6lvOuaop+1zzqa1fOWuR+b13heddfzTXmeutp5xLPc/dGDe73Pc0bMn/nM2rS1b1wue3cU3v//YjPu94Nld8+7PhWcdyxe//WjF7TO95kL3Lfczn81sn/dstSzks5mrsZEx1ravnPE2UlM1NkL7ysaq1TJfxx29iuFDw/QeGuaFz9nA//vhzF/3Fz5nw4LrXsjPfLmYS99n+/010++EWlevP/f2lY0V55DN5Ttcr31frPn2+9ItJzJ0aJidu/bRN3CYzvZVbO5ez6VbTqza5zdTYKy1U5Z3AieFECZTw/nArVkXcez6thmvWmtuaig7yXzTxs55v8+Lzjx+3vtMOufZz1jQ+6xbs3rGFaAnzdSX005eV3G/2baX8/KXnLrg7YvZd75m+7xne6+FfDZz1bKiibNDxYHkI04+pqNm559cecnpi9quhZvt99dCfr9p4SrNtTy+q61mv8P1qKmxkSu2dXPD63+FP/3dF3DD63+FK7Z1L/hWgItRU4EsxjgEvAH4WAjhBuAHMcavpajlI286v+xq4h950/kz7vOhN75wzq//7qvOW3Btk/7kd39lQe/zgTdsoXVV+R/9bH34wO9vmVf7XLztt86aV/tS7TtfM33ek+3V+Gzm6vILN7Ht3I1H7r7Q2Va89cr0eVoffNOLq17LYrzh3z93Xu1aOjN99+fze01LZ6a5lte+9uy0hS1TLSua2NDZmvTWcg0TE+muuFoKvb0Di+rAbMObj+8b5L7d+zlz07oZl1+YbveePr78nUdY0dTEr73gRE46Zg3fum8vX/3eo/zquScsamSsnHse+AW33bOHbeds5FBhdE7v09XVzgO7e7h95x56+w9z0dkb53UU/OOH93P79x9j61nHLdnoz+f/34PsuPcxLjz7uHmPbs1n38UO5U/9vMuNnFXjs5mrwsgYBwcLrGlroWVFEwNDw+zpGWTjhuJRdb2cxvjbW+/n2z/q4YXP3bAkI2P10u9qmG/fd+/p444f/pzzzzim7kfGlsPPffp3eK6WQ98Xoh763dXVXvYUnIGsDn541ZDXfoN9z2Pf89pvsO/2PV/qod8zBbKaOmUpSZKURwYySZKkxAxkkiRJiRnIJEmSEjOQSZIkJWYgkyRJSsxAJkmSlJiBTJIkKTEDmSRJUmJ1v1K/JElSvXOETJIkKTEDmSRJUmIGMkmSpMQMZJIkSYkZyCRJkhIzkEmSJCXWnLqAVEII24BXAj3ARIzxvYlLykQI4RjgBuDMGON5qevJSgjhFIr9vhfYCOyPMb4vbVXZCCE0Av8EfAdYCZwCXBVjPJS0sAyFEFZT7P9XY4xXp64nKyGEu4DDpYdjMcaLUtaTlRBCAF4FHAK2AtfHGL+btqrqCyGcDHwNeLTU1AH8IMZ4ZaqashRCeDtwMrAPOBV4XT39nstlIAshtAIfB54bYyyEED4XQrgoxvi11LVl4EXAF4CzUheSsbXAP8QYvwAQQvi3EMKtMcZ7EteVlTtjjDcAhBC+QPFg5O/SlpSpG4CdqYtI4MsxxutTF5GlEEIT8BHg0hjjeAjhZmA0cVlZGQB+L8Z4G0AI4XrgtqQVZaQ02PBOYH3p5153v+dyGciALcAjMcZC6fEdwCUUjyyWtRjjLSGEC1LXkbUY493TmhqBJ1PUkrUY4zjFQEIIoZniCGFMWlSGQgivofgdfx7QlricrJ0RQngHsBq4O8Z4a+qCMnAe0AC8qXTwvR/4ZNqSshFj3E8pgIUQWoBzcxTIh4BhiqOCT1D8rv8oaUXzlNc5ZBsoHklM6i+1KQdCCK8AvhJjfCB1LVkKIVwM/DPwzzHG76WuJwshhOcAp8UY/2/qWhL5QIzxA8CfAO8KIbwkdUEZOIniQfffxhjfD7wE+O20JSXxKuAfUheRlRhjP/B2YHsI4W+BPcDupEXNU14DWQ/QPuVxR6lNy1wI4aXAS4E/Sl1L1mKMX4kx/hrwzBDC76euJyOvAA6HEP6Y4un654cQ3pK4psxMzpuKMY4B36T43/5y1w88EGM8WHr8LeCCdOUk8x+B7amLyEoI4SyKgeyS0py5fcB7khY1T3kNZHcCJ5WGdAHOB/IwlJ9rIYRLgIuBNwPHhBC2JC4pEyGE55T6PumnwLNS1ZOlGOONMcb3xRhvoviH+bsxxo+mrisLIYRnhxBeN6XpVOChVPVk6DvAutJcMiiOmO1KWE/mStNS7owxjqSuJUPHAwdijJPzBR8HViWsZ95yOYcsxjgUQngD8LEQQi/Fq1CW/fwxgBDCVuA1wLEhhOuAD9fTVSgLFUI4h+LR4veArwNHAX9BMZwvdwXgdSGEzcAK4DTgD9OWlK0QwmUUT12tDCG8Ksb496lrykA/cEkI4TiKZwEeBT6btqTqizEeKM2b+2jp93sXkIsrqqf4PeBNqYvI2JeBXw8hfJjiHLLTgboaDW+YmJhIXYMkSVKu5fWUpSRJUs0wkEmSJCVmIJMkSUrMQCZJkpSYgUySJCmxXC57IUkhhDOBj1K8jdqVpZuwP0jxdjN9Ffa7Erhg8obNc91PkipxhExSLsUY7wP+15TH48C2+Yaqhe4nSVO5DpmkmhdC+C/Au4F/objQ7enA54DRUvvfUbz7wAXA64CdwHsprtZ9AvDZGOO/lF7r/cC5wE+AcWB1aYTs9cD1wH+OMX4jhHAycBPwM4qrvUfgfwN/RXFV8K8B/wScWGa/91NciPVE4FMxxn8tLUb9buD/lNpPB66OMX5xqT8vSfXHETJJNS/G+HGKpwUfjzG+Efh14FqKt0P6KnByjPEy4FLgYeAzwF/HGN8JvAH4mxDC0SGE3wBeBlwcY/w9ptzTNsb4ydJ7TPo74O9jjNdQvLvFmTHGB0uv/Z0Y4xtL9wctt98/lPZ7A/C/QwjPiDH+VanWzhjjK4E3Au9cwo9JUh0zkEmqJ3dA8fZnFG+DNXmz7NtK7d+meL/GFwFXhhA+Dnyw1HYCcCHwrdJpRijecPtpQgjtwAunvN9wjPE/zFZcmf36gD3Ai6c87dul/38IeMasPZaUC07ql7QcFMq0XRtj/AVACGE1MJxtSTOarHUMD4ollfjLQFI92QIQQmilOA/s69OfEGMcoDhC9aul5zYCXwJagB3Ai0ttUBxJe5opr3F+6TVWhRC2lzYfBppCCA0hhN+eZb9OYCMzjMRJ0iQn9UuqCyGEbwB3A2uAANwC/JTiKcm9wMcmJ8iHEJ4JfJjixPqjgFtijF8ubXs/cA6wG+gEzgLeAxxNcXL+dynO71pJcVL/XmAd8D9jjHeEEI6neEHBbooBr6nCfhspzmX7amn+2mStvwO8C/gt4N0xxj9f4o9LUp0xkEmqC6VAdn2M8RuJS5GkJecpS0k1L4Tw+0A38NYQwsbU9UjSUnOETJIkKTFHyCRJkhIzkEmSJCVmIJMkSUrMQCZJkpSYgUySJCkxA5kkSVJi/x+qmPf4A09/dQAAAABJRU5ErkJggg==\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 ZIP to ordinary Poisson 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
}