11 Chapter 11/05 H6.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.style.use('seaborn-whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize(alpha=0.05):\n",
    "\n",
    "    def _summarize(ser):\n",
    "        hpd_ = pm.hpd(ser, alpha=alpha)\n",
    "        return pd.Series({\n",
    "            'mean_': ser.mean(),\n",
    "            'lower': hpd_[0],\n",
    "            'upper': hpd_[1],\n",
    "        })\n",
    "\n",
    "    return _summarize\n",
    "\n",
    "\n",
    "def agg_prediction(trace_var, source_data, group_columns, aggregate_fc):\n",
    "    df = pd.DataFrame(trace_var).stack().reset_index()\n",
    "    merged = pd.merge(\n",
    "        df, source_data,\n",
    "        left_on='level_1', right_index=True\n",
    "    )\n",
    "\n",
    "    return (\n",
    "        merged.groupby(group_columns)[0].apply(aggregate_fc)\n",
    "        .unstack().sort_index().reset_index()\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**H6** The data in `Fish` are records of visits to a national park. \n",
    "\n",
    "- `fish_caught` : Number of fish caught during visit\n",
    "- `livebait` : Whether or not group used livebait to fish\n",
    "- `camper` : Whether or not group had a camper\n",
    "- `persons` : Number of adults in group\n",
    "- `child` : Number of children in group\n",
    "- `hours` : Number of hours group spent in park\n",
    "\n",
    "The question of interest is how many fish an average visitor takes per hour, when fishing. The problem is that not everyone tried to fish, so the `fish_caught` numbers are zero-inflated. As with the monks example in the chapter, there is a process that determines who is fishing (working) and another process that determines fish per hour (manuscripts per day), conditional on fishing (working). We want to model both. Otherwise we’ll end up with an underestimate of rate of fish extraction from the park. You will model these data using zero-inflated Poisson GLMs.\n",
    "\n",
    "Predict `fish_caught` as a function of any of the other variables you think are relevant. One thing you must do, however, is use a proper Poisson offset/exposure in the Poisson portion of the zero-inflated model: \n",
    "$\\mu_i =  \\tau_i ~\\alpha~ e^{ ~\\beta x_i }$.\n",
    "Then use the `hours` variable to construct the offset. This will adjust the model for the differing amount of time individuals spent in the park.\n"
   ]
  },
  {
   "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>fish_caught</th>\n",
       "      <th>livebait</th>\n",
       "      <th>camper</th>\n",
       "      <th>persons</th>\n",
       "      <th>child</th>\n",
       "      <th>hours</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>87</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>9.611</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>80</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0.664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>236</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2.583</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>60</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2.955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1.329</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     fish_caught  livebait  camper  persons  child  hours\n",
       "87             1         1       1        1      0  9.611\n",
       "80             3         1       1        4      1  0.664\n",
       "236            2         1       1        2      0  2.583\n",
       "60             2         1       0        2      0  2.955\n",
       "50             1         1       1        2      1  1.329"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv('../data/Fish.csv', sep=';')\n",
    "data.sample(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [logit_psi, alpha]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:10<00:00, 1180.58draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_basic:\n",
    "    base_rate = data.fish_caught.sum() / data.hours.sum()\n",
    "    alpha = pm.Exponential('alpha', lam=1/base_rate)\n",
    "\n",
    "    logit_psi = pm.Normal('logit_psi', mu=0, sd=1.5)\n",
    "    psi = pm.Deterministic('psi', pm.math.invlogit(logit_psi)) \n",
    "\n",
    "    fish = pm.Deterministic(\n",
    "        'fish',\n",
    "        data.hours.values * alpha\n",
    "    )\n",
    "\n",
    "    fish_obs = pm.ZeroInflatedPoisson(\n",
    "        'fish_obs',\n",
    "        psi = psi,\n",
    "        theta = fish,\n",
    "        observed = data.fish_caught.values\n",
    "    )\n",
    "    \n",
    "    m_basic.name = 'basic'\n",
    "    m_basic.trace = pm.sample(5000, tune=1000, chains=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEcCAYAAACBPmBXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEdFJREFUeJzt3HuUnHV5wPHvcqkQkoKHhIsWiKbyGCuK3cZGuQUvKOo5oLYVqTWk2IpwPCKechEF6cV6RUGreAmIepAjp1LFCwUlRCBQMN4QlwcEg6hcIkeExASEpH+878qy7C67G3aemc33c86e3Z2Z953fb2Znv/ObnX37Nm7ciCRJnbZF9QAkSZsnAyRJKmGAJEklDJAkqYQBkiSVMECSpBJbVQ9AmqiIeBJwDvB04D7gmMy8OSJeDXwIuL296KnASuBrwLbAmzPzxxGxL7BPZr5/lP3/M/AGYAOwNXByZl4+iXEeBewCnAWckplHR8Qq4JmZuX7I5Y5oTztxHPv8HHB+Zl48yvm7A8/NzIsmOt7Hud5FwFGZedgTuV9t3lwBqRf9E7AmMxcCbwU+3p7eDxyfmYvaj+XAQTQBOho4MiL6gLcBZ4y044g4DHgp8OLMXEQToi9ExOzJDjYz78zMoye7/QS9CNinQ9clbRJXQOpFzwK+BZCZGRHz29P7gedFxLHAtcAJwBqa1c+2wFrgcODCoSuQYd4MHJeZf2j3//OI2Dsz72lXHzu2H68Ejgf2A7YETs/MC9rV1RnAb4GHgGsiYi7NqmVhex2fak+7C1g89Moj4q3tGDe225w50iDbFckJwIM0K8HzgfcBJwIzImIF8HPgTKAPuAf4R+B5wPvb7b4NLMrMA9t9fh14NzAPOIZm9bcRePWw6z4H+PP2Nj0jM78wym0pjckVkHrRD4FXRURfRCwEnhoRWwKX0qyI9gdmAkfR/JLdmSYsn6b5ZfqjiPhURBw/wr6fAtw69ITMvGfIt5dl5guBhcDTMnNf4EDg5IjYAfgk8PrMfAlNAEbyycw8AFhFs5oDICKeBbwO2JcmbIdGRIxxO+wBvLYdy/GZ+TBNhM7LzK8Bn6F5eXIR8E2aYAJsk5n7ZeZpwDYRsUdE7ArMzswfAHsCr2zn9lPgZUPGOIvm9n0N8HLg4THGJ43JFZB60dnAfOAK4CpgZWY+HBFnZ+a9ABHxVeC1mbkBOLY97SSa1cm7aEL1nojYMzNvGrLv24DdgN8NnhARLwN+3H6b7ee9gP6IuLz9fmtgLrDzkP1dRbNSGOrBzLym/XoFzct917XfP5smKt9pv38y8Iwh1znc9Zn5EPBQRKwb4fz5wCfahm0N3DxsDgBLgTcCD9D8XQ3gbuDciFgDPBO4evDCmXl/u8L8NPCnwBdHGZv0uFwBqRctAL7TPkO/ALi1/dvOjyPiz9rLvJjmDQgARMROQGTmFcAMmmfuG4Hthu37bODdEbFVu92ewGd55Jn+hvbzjcCydnXxIuDLwC3Ar4a8JLhghLH/SUTs3X69H/CTIeclcANwYLvfz/FI+EYy0oEcN/DI4zqBN7b7Oh74+rA5QPPS3atoVobnRcT2wGnAYcCbgHU0L+EB0K6U+jPz1TQvQ35g8LaSJsoAqRfdDBwbEVcD/0bzN5uNNL8wvxIRy2ki85kh27wL+Pf2608A/wvsCvxo6I4z83zgGuDKiPguzargDZl597AxXASsiYgraEK3MTPvp3mp7/MR8R2a1cxwDwBvbfe9E03cBq/7RzSrnysj4ns0q59fjf9mAeB64JD2zRRvacdyJc1Lc4+JWWauaW+DgXb899Gs3K6mWWGuo3lZctCdwC7t35guBT7UrsKkCevzaNiSpAqugCRJJQyQJKmEAZIklTBAkqQSBkiSVKJr37+/cuVK354nSdNAf39/30ind/PbsDdpYAMDA8yfP//xL9gjnE93cz7dzfmUGjE+4EtwkqQiBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSia2qByBJ3egX9/yeI8+9jltXr+Xpc7Zj6eIF7L7jjOphTSsGSNK0seSca1mWq0c599ZJ7/fmu9ew/weXTXr70RwYczhnyfOf8P32CgMkaVwO+shybrprTfUwppVluZq5J35jkltPPqiPZ8+dZ3LJ2w+Ysv0PMkCSxqUTv5CmysDAAPPnz5/QNi89fTm3rF7Dho2wRR/MmzOTS4/rjttgMvPpRr4JQZJGsHTxAubNmcmWfX3MmzOTpYsXVA9p2nEFJEkj2H3HGV2z4pmuXAFJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIkldjkAEXEnWOctygizt/U65AkTT/+H5AkjcIDkk6tcQcoIp4NnA5sCcwG3pKZK4acfzlwI/BMoA94XXvWMyLiW8BOwEWZ+Z6IOAA4lWYFNhM4PDNv2vTpSNqcTdXBSGFqDkjqwUjH7y+Ad2Tm9RFxOLAEWDHsMisy86iIOBp4J/AVYBvgUJpw/QJ4T7uvN2TmryPincDfAv8xdEerVq1i3bp1k5hSY/369QwMDEx6+27jfLrb5jKfo756O7fd+4eCEU1P3Xow0qH22GFrzjpkt0lvP9Yx6yYSoF8B746IdcAs4L4RLnNZ+3kFcEj79U8y8wGAiHhoyL7OjIg1wFOBq4bvaO7cuRMY2mNNl4P1DXI+3W1zmc/yHp3jZO+fbj0g6XT5eZvImxDOBE7NzMXA9TQvsw3X337eB7ih/XrjCJf7DLAkM48Afj3KviSplAcknVoTWQF9EbggIn4L/JLm70DDHRERxwFrgX8A9hpjX1dExFrgLuApExiHJHWEBySdWuMOUGaeTvMmhOGn7zLk25My88Yh31/efjzqspl53EQHKkmaXvxHVElSiSfs/4Ayc9ETtS9J0vTnCkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSW2qh6Aussv7vk9R557HbeuXsvT52zH0sUL2H3HGdXDkjQNGaAutuSca1mWq4eccmtHr//mu9ew/weXTeE1PHo+B8Yczlny/Cm8PkndZFoH6KCPLOemu9ZUD0PjtCxXM/fEb1QPYxOM/ARhz51ncsnbD+jwWKTuN60DNJ0e9AMDA8yfP3/Kr+elpy/nltVr2LARtuiDeXNmculxT/zt2Kn5dMp0m4/UCb4JQY+ydPEC5s2ZyZZ9fcybM5OlixdUD0nSNDWtV0CauN13nDElKx5JGs4VkCSphAGSJJUwQJKkEgZIklTCAEmSSnRFgCLiK9VjkCR1Vle8DTszX1M9Bj3C48FJ6oQpC1BEfA7oA3YDZgJvBFYBXwa2B2YAJ2fmJRFxZ2buMlVjmQ5O+fYdXHduZ48FB1N9PDiPBSdtzqZ6BXRLZi6OiFcAHwBOAmYDLwd2AvYcbcNVq1axbt26SV/x+vXrOeB9l3DbvX+Y9D7UWdPpWHB77LA1Zx2yW+FYNs369esZGBioHsYTxvnUGesQVVMdoMvazyuAj2TmDRHxKeBLwNbAmaNtOHfu3E264oGBAZafeNAm7aObdPJYY504Htx0O3aa8+luzqc7TfWbEPrbz/sAN0TEXsCszHwlsBj42BRfvybB48FJ6oSpXgEdHBGHAFsCRwB3AKdGxN/RxO+UKb5+TYLHg5PUCVMdoI9m5sXDTvub4RfyDQiStPnpiv8DkiRtfqZsBZSZR0zVviVJvc8VkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVIJAyRJKmGAJEklDJAkqYQBkiSVMECSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJQyQJKmEAZIklTBAkqQSBkiSVMIASZJKGCBJUgkDJEkqYYAkSSUMkCSphAGSJJUwQJKkEgZIklTCAEmSShggSVKJvo0bN1aPYUQrV67szoFJkiakv7+/b6TTuzZAkqTpzZfgJEklDJAkqYQBkiSV2Kp6AJMREVsAnwCeCzwAvCkzfzbk/IOBU4E+YCVwDLAN8EVgJ+B+YHFmru7w0Ec0yfkA/BK4uf366sw8qWODHsNY84mIvYGPDrn4QuBQ4HvAecC2wK+BJZn5+06OezSTnM+1wE3AT9rTL8zMMzo26DGM4+ftHcDhwAbgvZl5YURsS+8+fkaaTx89+Phpzz8BeD1wH/CBzPx6RMymSx8/Y+nVFdChwDaZ+QLgRODDg2dExCzgg8CrMvOvgVXAbOAtwPWZuR/weeBdnR70GCYzn3nA9zNzUfvRFQ+e1qjzycwfDo4Z+C/gvzPzYuAU4Lz2/vkB8ObOD3tUk5nPXwJfGnL/dEV8WmP9vO0AvA14AXAQj8S1Vx8/o82nJx8/EbEXTUwX0sznXyNiBt39+BlVrwZoX+BigMy8BvirIee9ELge+HBEXAHc1T5T++M2wLeAl3RuuI9rMvPpB54aEcsi4psREZ0e9BjGmg8AEbEdcBrNL4dHbUNv3T/AiPPpB/ojYnlEXBARu3ZqsOMw1nzWArcB27UfG4ZvQ2/dP6PNp1cfP/OByzNzfWaup1nBPYfuvn9G1asB+lPgd0O+fzgiBl9OnA0cCJwAHAwcGxF7DtvmfmD7Do11PCYznzuA/8zMA4H30rw80i3Gms+gI4ELMvM3I2zTS/fPoOHzuRE4JTMPAP4H+NjUD3PcHm8+twM/Bb4PnDnCNr12/4w0n159/FwP7B8RsyJiR5onqNvR3ffPqHo1QPcBs4Z8v0VmPtR+fQ9wXWbemZlrgO8Cew/bZhZwb6cGOw6Tmc/3gK8CZOaVwFPa17W7wVjzGfT3wGdH2aaX7p9Bw+dzGbCs/fpC4HlTN7wJG2s+BwO7Ak8DdgcOjYjn07v3z2jz6cnHT2YOAB+nWe18HPg/4Dd09/0zql4N0FXAKwAiYiHNs4JB3weeHRGz22cNC2me/fxxG5ofyis6N9zHNZn5nAoc227zXOD2zOyW/yoeaz5ExPbAkzLz9pG2obfun9Hm81ngte3XL6Z580i3GGs+vwXWAQ+0L/HcC+xA794/o82nJx8/ETEHmJWZ+wBHAbvRvNGlm++fUfXkkRCGvEvkOTTvDFtCc+P/LDO/FhGHAf/SXvzLmfn+9g9159I8G3oQODwz7+z86B9rkvN5Ms3LBjOBh4BjMvPGzo/+scYxnwXAyZl56JBtdqa5f2bRPKM7PDPXdnzwI5jkfJ4GnN1efi3NO5nu6PjgRzCO+ZwGvJzm7yVXAsfTvLuqVx8/I81nB3rw8QNcBJxF8yaXB4GTMvO73fz4GUtPBkiS1Pt69SU4SVKPM0CSpBIGSJJUwgBJkkoYIElSCQMkSSphgCRJJf4fXeQGdZKj9gcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_basic.trace, varnames=['alpha', 'psi'], rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [alpha, beta_psi_bait, logit_psi]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:28<00:00, 418.61draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_bait:\n",
    "    base_rate = data.fish_caught.sum() / data.hours.sum()\n",
    "\n",
    "    logit_psi = pm.Normal('logit_psi', mu=0, sd=1.5)\n",
    "    beta_psi_bait = pm.Normal('beta_psi_bait', mu=0, sd=1)\n",
    "\n",
    "    alpha = pm.Exponential('alpha', lam=1/base_rate, shape=2)\n",
    "\n",
    "    psi_i = pm.math.invlogit(\n",
    "        logit_psi\n",
    "        + beta_psi_bait * data.livebait.values\n",
    "    )\n",
    "\n",
    "    fish = (\n",
    "        data.hours.values\n",
    "        * alpha[data.livebait.values]\n",
    "    )\n",
    "\n",
    "    fish_obs = pm.ZeroInflatedPoisson(\n",
    "        'fish_obs',\n",
    "        psi = psi_i,\n",
    "        theta = fish,\n",
    "        observed = data.fish_caught.values\n",
    "    )\n",
    "    \n",
    "    m_bait.name = 'livebait'\n",
    "    m_bait.trace = pm.sample(5000, tune=1000, chains=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEcCAYAAADpzeJvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8VPW57/FPpFTBKLIBpbUoB6qPtHZrDwePtSpohSot0lJb2GyU4F20bqrWokW8lCotikhrkXsUjg0VVKRFKyqyUfTAsVrFxgfkGrlGCkgYQJScP2ZR0zSBZMiaNfnN9/168WIyl7W+awbyzPN71iQFlZWViIiI1NdhSQcQEZHGSQVEREQyogIiIiIZUQEREZGMqICIiEhGVEBERCQjn0s6gEhczOxwYCrQAfgIuN7dl5vZ94H7gbLorncCbwDPAM2Aa9z9bTM7G/imu/+qlu1fDQwA9gFNgZ+7+8sZ5LwWaAs8Agx398Fmtho4xd13V7lfUXTd0Dpssxgocffnarn9BOA0d59T37wH2W834Fp379eQ25XcpA5EQnYVUOHuZwI/Bn4bXd8ZuNXdu0V/FgA9SBeQwcAVZlYA/BfwUE0bNrN+QHfgW+7ejXQhmWZmrTMN6+4b3X1wpo+vp/OBb2ZpXxIodSASsq8AzwK4u5tZp+j6zsDXzWwIsBj4GVBBuvtoBuwE+gNPVe0AqrkGuMnd90bbX2Vmp7v7lujdf6voz3eAW4FzgCbAaHd/IupuHgK2Ap8Ar5tZe9Jdw5nRPsZH120CBlbduZn9OMpYGT1mbE0ho47gZ8DHpDuxEmAkMBRobmaLgFXAWKAA2AJcDnwd+FX0uBeAbu5+XrTNPwJ3AB2B60l3X5XA96vteyrw5eg5fcjdp9XyXEojpQ5EQvYW8F0zKzCzM4HjzawJMI90R3IuUAhcS/qb5HGkC8ME0t8M/2pm483s1hq2/UVgZdUr3H1LlS9fcvezgDOB/+HuZwPnAT83s2OAccB/uPsFpL+B12Scu3cFVpPupgAws68AfYGzSRem75mZHeB5OBH4QZTlVnf/lHQRedzdnwEmkl7e6wbMJV3wAI5w93Pc/W7gCDM70cy+ALR29zeBk4HvRMf2N+DbVTIeRfr57QNcCHx6gHzSSKkDkZBNAToBC4FXgTfc/VMzm+Lu2wDMbDbwA3ffBwyJrruNdHcwjHShucvMTnb3ZVW2vQZoB2zff4WZfRt4O/rSo7+/BnQ2s5ejr5sC7YHjqmzvVdLv1Kv62N1fjy4vIr1ctiT6+lTSReHF6OuWwElV9lndO+7+CfCJme2q4fZOwO+iGtQUWF7tGAAmA5cBe0jPlQA2A4+aWQVwCvDa/ju7+46ow5sAHA1MryWbNGLqQCRkXYAXo3fITwAro9nG22b2peg+3yI9QAfAzI4FzN0XAs1Jv3OuBI6stu0pwB1m9rnocScDk/jsnfa+6O/3gPnRu/vzgT8AK4B1VZbUutSQ/fNmdnp0+RxgaZXbHHgXOC/abjGfFa6a1PQD7/bx2f9/By6LtnUr8MdqxwDppa/vku7MHjezFsDdQD/gSmAX6SUwAKJOpbO7f5/0Mt6v9z9XEg4VEAnZcmCImb0G/IL0zKKS9De8J81sAekiMbHKY4YBI6LLvwP+DHwB+GvVDbt7CfA68IqZ/Tfpd+UD3H1ztQxzgAozW0i6UFW6+w7SS2WPmdmLpLuJ6vYAP462fSzp4rR/338l3X28Ymb/j3T3sa7uTwsA7wC9o5MBrouyvEJ6aetfipG7V0TPQWmU/yPSndNrpDu8XaSX9fbbCLSNZizzgPujLkgCUqCfxisiIplQByIiIhlRARERkYyogIiISEZUQEREJCMqICIikhGdl10Hb7zxhk5VE5G81blz54KartdpvHUT1JNUWlpKp06dDn7HQEycOJENGzYwfPjwpKNkVb69zpB/x5yl462xeIA6EMkDV111FaWlpUnHEAmOZiAiIpIRdSASvG7dupFKpVi8eHHSUUSCogIiwSsqKmL9+vVJxxAJjgqIBK+oqEgzEJEYaAYiwdu7dy979+5NOoZIcNSBSPC6d++uGYhIDFRAJHhXXnkl69bV99dliMjBqIBI8AYMGKAZiEgMNAOR4KVSKXbtqulXgYvIoVAHIsHr2bOnZiAiMVABkeBdd911moGIxEAFRILXt29fzUBEYqAZiARv+/bt7NixI+kYIsFRByLB6927t2YgIjFQAZHg3XjjjXzwwQdJxxAJjgqIBK9Pnz6agYjEQDMQCd6HH37I1q1bk44hEhx1IBK8Sy65RDMQkRiogEjwbr75ZsrKypKOIRIcFRAJXq9evWKdgazdkuKKR5ewsnwnHdocyeSBXTihVfPY9ieSK1RAJHgbN26kvLycTp06Ndg2B01dzHwv/5frl2+u4NxR8wE4z9owddAZDbZPkVyjAiLB69evX8YzkB4PLmDZpoqM9jvfy2k/9E8HvM/JxxXy/E+6ZrR9kaSpgEjwhg4dytq1azN6bF2+uXcfvYAV5RXsq4TDCqBjm0Lm3aSiIOHTabwSvAsvvJBzzjkntu1PHtiFjm0KaVJQQMc2hUwe2CW2fYnkkpztQMxso7u3reW2bsC17t6vntssAu4BxgCPAY8DzYD1wCDgcuAWYKS7P5JxeMkpZWVlbNiwoUFnIFWd0Kq5Og7JS/nYgTzu7qOB4dHlc4A3gWvc/bdAcZLhpOFdeumlDB06NLbtr92SovvoBXS8bS7dRy9g7ZZUbPsSySWJdyBmdiowGmgCtAauc/dFVW5/GXgPOAUoAPpGN51kZs8CxwJz3P0uM+sK3Em6MBYC/d19WS27Phu4N7r8bHT5wQY8NMkRw4YNY82aNQ26TZ2FJZIDBQT4KnCzu79jZv1JLyUtqnafRe5+rZkNBm4HngSOAL5HuvCsBe6KtjXA3deb2e3AD4Ff1rLfo4Ht0eUdQIvaAq5evTqoX4m6e/fuvPrZUMcffzytWrXK6JivnV3Gmm17M9pvXc7COvGYpjzSu11G2z+YfHudIf+OORvHe6Cl31woIOuAO8xsF3AU8FEN93kp+nsR0Du6vNTd9wCY2SdVtjXWzCqA44FXD7Dfj6L97d/vttru2L59+zodSGNRWloa2zwgF61cuZKysjJ69OhR78cuqMPzlKtnYeXb6wz5d8xJH28uzEDGAne6+0DgHdLLVNV1jv7+JvBudLmyhvtNBAa5exHpwXhN29rvVaBndPkiYGH9YktjcfnllzNs2LDYtq+zsCRf5UIHMh14wsy2Ah+QnoNUV2RmNwE7gUuBrx1gWwvNbCewCfjiAfY7AnjUzK4CPgT6Z5hfctzdd9/d4DOQqnQWluSrxAtIdEbU6Bqur3oK723u/l6Vr1+O/vzTfd39pnrsdxNwYT3jSiPUtWvXvFoXF8mWXFjCyrb+UTfzL8zsBqAou3Ekbu7OqlWrko4hEpyCysqaRglSTVBPUtKDt2zr1q1bXv4+kHx7nSH/jjlLx1vrLDnxJSyRuN17772sXr066RgiwVEBkeCdddZZtGzZMukYIsHJxxmI5JmlS5eyfPnypGOIBEcdiATvhhtuIJVKcfHFFycdRSQoKiASvFGjRuksLJEYqIBI8Lp06UJhYWHSMUSCoxmIBO+tt97SBwlFYqAORII3ZMgQUqkUffr0STqKSFBUQCR4Y8aMYeXKlUnHEAmOCogE7/TTT+fwww9POoZIcDQDkeAtWbKEd955J+kYIsFRByLB++lPf0oqleJHP/pR0lFEgqICIsH77W9/qxmISAxUQCR4p556Kk2aNEk6hkhwNAOR4C1atIg333wz6RgiwVEHIsG7/fbbSaVS9O+v31os0pBUQCR448ePZ8WKFUnHEAmOCogEz8zYt29f0jFEgqMZiARvwYIFLFmyJOkYIsFRByLBu/POO0mlUlx22WVJRxEJigqIBG/KlCm8//77SccQCY4KiASvQ4cO7NmzJ+kYIsHRDESC98ILL7Bo0aKkY4gERwVEgjdixAjGjx+fdAyR4GgJS4I3bdo0li9fnnQMkeCogEjw2rVrR0VFRdIxRIKjJSwJ3nPPPcfChQuTjiESHBUQCd7IkSOZNGlS0jFEgqMlLAleSUkJy5YtSzqGSHBUQCR4bdu2ZevWrUnHEAmOlrAkeHPmzGH+/PlJxxAJjgqIBO+BBx6guLg46RgiwdESlgRv5syZmoGIxEAdiASvdevWtGzZMukYIsFRAZHgPfnkk8ybNy/pGCLBUQGR4I0dO5bp06cnHUMkOJqBSPBmz56NuycdQyQ46kAkeC1atOCoo45KOoZIcNSBSPBmzJjBunXr6NSpU9JRRIKiDkSCN27cOEpKSpKOIRIcdSASvLlz5/Lee+8lHUMkOOpAJHjNmzenWbNmSccQCY46EAne9OnTNQMRiYE6EAnepEmTmDVrVtIxRIKjDkSCN2/ePEpLS5OOIRIcdSASvKZNm9K0adOkY4gERx2IBK+4uJj169drBtIIrd2S4opHl7CyfCcd2hzJ5IFdOKFV86RjSUQdiASvuLiYp59+OukYUk+Dpi7m3FHzWb65gk8rK1m+uYJzR81n0NTFSUeTiDoQCd7LL7+sGUiO6/HgApZtqqjTfed7Oe2H/umfrjv5uEKe/0nXOKLJAaiAiEjiavvm3330AlaUV7CvEg4rgI5tCpl3kwpFrtASlgRv4sSJPPHEE0nHkAxMHtiFjm0KaVJQQMc2hUwe2CXpSFKFOhAJ3owZM9i5cyfDhw9POorU0wmtmqvjyGF5VUDMrAi4Bxjj7qPNbAjQ1t2HmtklwAjgaXcfmmROaVgvvPCCZiCNlM7Cym35uIT1ODDOzP4PcP3+K919JjAysVQi8k90Flbuy6sOpIojgEeBecApB7vz6tWr2bVrV+yhsmX37t159Y7897//PXv37uWyyy5LOkpWNabX+drZZazZtrdO963pLKwTj2nKI73bNapjbgjZON4DfX4qLwuIu28Fno+WtA6qffv2sebJttLS0rz6UN2SJUuoqKjgvvvuSzpKVjWm13lBLTnrexZWYzrmhpD08ebjEpbkmWeffZYJEyYkHUMyoLOwcltediAi0jjoLKzcpg5EgvfQQw8xbdq0pGOIBCevOxB3L046g8TvxRdfpKKibj8mQ0TqLh87kP5mdlP1K6PPgejzHwF65plnePjhh5OOIRKcvOpAoo6juJbbZgIzs5lHRKQxy8cORPLM/fffz5QpU5KOIRKcvOpAJD+99tpr7NixI+kYIsFRAZHgzZo1K68+nSySLVrCEhGRjKiASPBGjhzJxIkTk44hEhwtYUnw3nrrLT766KOkY4gERwVEgldSUqIZiEgMtIQlIiIZUQGR4P3iF79g3LhxSccQCY6WsCR47s727duTjiESHBUQCd706dM1AxGJgZawREQkIyogErzhw4fzm9/8JukYIsHREpYEr6ysjG3btiUdQyQ4KiASvKlTp2oGIhIDLWGJiEhGVEAkeLfddhujR49OOoZIcFRAJHhbtmzR50BEYqAZiARvwoQJmoGIxEAdiIiIZEQdiATvlltuYcuWLUydOjXpKCJBUQciwdu1axd79uxJOoZIcNSBSPAefvhhzUBEYqAOREREMqIORII3ZMgQ/v73v/PYY48lHUUkKOpAREQkI+pAJHhjxozRDEQkBupAREQkI+pAJHjXX389W7du5fHHH086ikhQ1IFI8Jo1a8bhhx+edAyR4KgDkeDdf//9moGIxEAdiIiIZEQdiATv6quvZtu2bfzhD39IOopIUFRAJHitWrWisrIy6RgiwVEBkeDdd999moGIxEAzEBERyYg6EAneoEGD2LZtG0899VTSUUSCogIiwWvXrh3NmjVLOoZIcFRAJHj33HOPZiAiMdAMREREMqIORII3YMAAtm/fzpw5c5KOIhIUFRAJnplRXl6edAyR4KiASPDuuOMOzUBEYqAZiIiIZEQdiASvX79+fPTRR8ydOzfpKCJBUQGR4J1++uls3rw56Rg5Z+2WFFc8uoSV5Tvp0OZIJg/swgmtmicdSxoRFRAJ3tChQxOfgQyaupj5nsQgf2Wd7rV8cwXnjpofc5bMnWdtmDrojKRjSDUF+imldRLUk1RaWkqnTp0yfnyPBxewbFNFAyYSkTidfFwhz/+ka6YPL6jtBnUgUm+H8A8xET/4wQ/YsWMHzz//fNJRsupgbxS6j17AivIK9lXCYQXQsU0h825qXK9tdYf65qixSfp4dRaWBO8b3/gGp512WtIxcs7kgV3o2KaQJgUFdGxTyOSBXZKOJI2MOhAJ3i233JL4DCQXndCqeaPvOCRZOdGBmFmxmV1Yz8e0NbPfRZfPNbN/b6AsY8zshIbYlohIyBptB+LuG4HB0ZeXAyXA2w2w3SGHug3JLRdffDEVFRW89NJLSUfJSTqdVzKVUwXEzJoCU4EOQBNgtLvPMLMzgIeBHcBmYDdwF+micT1wIfA/zexv7r62hu0WAd8DjgJaA/e4+ywz+yVwHunnYZa7/8rMXgaudff34jxWyZ5vfetbbNq0KekYCZ3KW7fTePfT6bxSHzlVQIBrgHJ3H2BmRwF/MbMXgUeAS9393eib/vH7H+Dub5jZc0BJTcWjiiOB7kAbYLGZzQb+E+gGbACKanvg6tWr2bVr16EdWQ7ZvXt3vWYC184uY822vTEmituXgS/z+NA/JR1EDtF8L6f9QV/H+hXNxq/m4z3xmKY80rvdIW/9QGd55VoB6QS8AODuO8zsb0BH4Ivu/m50n4VAvwy2vcDd9wGbzGwr6ULyn8BIoC3wbG0PbN++fQa7y131PfVvQQCnRSZ9umMS6nrMIZ3Om2+vc9LHmxND9CpKgXMAog7ka8AqoMzMvhLd58waHrePgx9L52i7xwFHA38Hfgj8B+llrCIzO/FQD0Byz0UXXcTVV1+ddIycpdN5JVO51oFMACaa2StAM+Bud99sZoOBKWZWAXwMrKv2uP8LjDSzVe5e29pM22g5rAUw2N33mNnfgdeBXcDzwIGWwKSR6tWrFxs3bkw6Rs7S6bySqZwoIO5eVOXLgTXc5Qygl7uXm9kI4GN3X03Ujbj7eGD8QXazwN2HVtvvPcA91e7Xre7JpTEYPHiwPgciEoOcKCB1sAl4PupAtlNzkSH6XMhXarhpRozZRETyUqMoIO4+E5hZh/sNPth9JP9ccMEF7Ny5k9deey3pKCJBaRQFRORQ9O3blw0bNiQdQyQ4KiASvKuuukozEJEY5NppvCIi0kioA5HgdevWjVQqxeLFi5OOIhIUFRAJXlFREevXr086hkhwVEAkeEVFRZqBiMRAMxAJ3t69e9m7tzH/MEiR3KQORILXvXt3zUBEYqACIsG78sorWbeu+o9PE5FDpQIiwRswYIBmICIx0AxEgpdKpYL6hWAiuUIdiASvZ8+emoGIxEAFRIJ33XXXaQYiEgMVEAle3759NQMRiYFmIBK87du3s2PHjqRjiARHHYgEr3fv3pqBiMRABUSCd+ONN/LBBx8kHUMkOCogErw+ffpoBiISA81AJHgffvghW7duTTqGSHDUgUjwLrnkEs1ARGKgAiLBu/nmmykrK0s6hkhwVEAkeL169dIMRCQGmoFI8DZu3Eh5eXnSMUSCow5EgtevXz/NQERioAIiwRs6dChr165NOoZIcFRAJHgXXnihZiAiMdAMRIJXVlbGhg0bko4hEhx1IBK8Sy+9VDMQkRiogEjwhg0bxpo1a5KOIRIcFRAJ3gUXXKAZiEgMNAOR4K1cuVKfRBeJgToQCd7ll19OKpWiR48eSUcRCYoKiATv7rvv1gxEJAYqIBK8rl27agYiEgPNQCR47s6qVauSjiESHHUgErxrrrmGVCpFz549k44iEhQVEAnevffey+rVq5OOIRIcFRAJ3llnnUXLli2TjiESHM1AJHhLly5l+fLlSccQCY46EAneDTfcQCqV4uKLL046ikhQVEAkeKNGjdJZWCIxUAGR4HXp0oXCwsKkY4gERzMQCd5bb72lDxKKxEAdiARvyJAhpFIp+vTpk3QUkaCogEjwxowZw8qVK5OOIRIcFRAJ3umnn87hhx+edAyR4KiASPCWLFnCqlWr6NSpU6I51m5JccWjS1hZvpMObY5k8sAunNCqeaKZRA6FCog0eoOmLma+lx/kXkdy61/+lJU8dbF8cwXnjpqfhT3Vf+nuPGvD1EFnxJBFQlNQWVmZdIbGICtPUo8HF7BsU0U2diUiWXLycYU8/5OusWy7tLQ0G511QW03qAPJIXH9I6suS//ockouHHP30QtYUV7Bvko4rAA6tilk3k3xvea5cMzZlo/HnCR9DkSCt2jRIt58882kYzB5YBc6timkSUEBHdsUMnlgl6QjiRySOncgZlYMlLj7c3W47wnAae4+5xCyHWwfJcBl7v5xLbdvdPe29dkW0JaYc0v23X777aRSKfr3759ojhNaNY+14xDJtriWsM4HTgFi+0bs7v0aeltmFntuyb7x48ezYsWKRPatM68kZHUeokcdyL8BR5IuPFe4+/tm9mOgP+lBcwnwMPAu0By4AdgO3El6uawQ6O/uy2rZx12kv4EfC7QEfuzur5jZVODLQDPgIXefZmargVPcfXct29oK/BloB7wNDAaOB8YBRwBfAIa5+9PRtr4KvLE/t7s/U2VzQZ1pEMI6cd3OvJL6COHsqxD+bddHYxuiL3L3kWbWE/i1mQ0D+gJnR7fPI/1NeyTpb+7PmNlgYIC7rzez24EfAr88wD5S7n6+mX0VeNzMzgbOBc4k/Y28Rx2zNgN+5u5rzOwPQC8gBTzg7i+b2VnA3cDT0f0/rZq76oZWr17Nrl276rjbf3bt7DLWbNub0WPjpU9myz+b7+W0H5o7pzpnLv5/2yce05RHereLfT8Hs3v37th/ztuBClR9C8h/R38vAkYBpwInAi9G17cETqr2mHXAWDOrIN0BvHqQfbwE4O7vmllbd99hZkOACcDRwPQ6Zl3r7muq5DVgLjDMzK4gXYya1mVD7du3r+Mu/9WCHHw3lG/v0rp160YqlWLx4sVZ33e2z7yqKt9eZ8i/Y076eOt7Ftb+/vYcYCngpJerznP3bkAx6eWifVW2PREY5O5FwHoO0A5FOgOY2anAOjP7AtDZ3b8PfId051OXwvel6LGQ7pCWAr8AHnP3S4H5NWSpmlsCMWXKFEaMGJHIvnXmlYSsvh3ImWb2Eul375dHy0MvAq+Y2eHAYtIdxzvAz83sL6Q7hoVmthPYBHzxIPv4erTNI4GrgI1AWzNbRHqZ6X53/8TMDpZ1C+nO50ukl96eNbNjgPvN7DbgA6B1tcf8I7e7l9Th+ZBGoEOHDuzZsyeRfevMKwlZTn0SPRqib3T3R5LOUk3uPEkNIOm2N9teeOEF1qxZwxVXXJF0lKzKt9cZ8u+YG9sQvUGY2ZOkz+iqajtQr097mdnFwE013PSQuz+VYTwJzIgRI0ilUnlXQETilkgBcfcG+c0+0dlSzxz0jpLXpk2bxvLly5OOIRIc/SwsCV67du2oqNAPqRRpaDrjSIL33HPPsXDhwqRjiARHBUSCN3LkSCZNmpR0DJHgaAlLgldSUsKyZTX+9BwROQQqIBK8tm3bsnXr1qRjiARHS1gSvDlz5jB/fjZ+faxIflEBkeA98MADFBcXJx1DJDhawpLgzZw5UzMQkRioA5HgtW7dmpYtWyYdQyQ4KiASvCeffJJ58+YlHUMkOCogEryxY8cyfXpdf42MiNSVZiASvNmzZ+PuSccQCY46EAleixYtOOqoo5KOIRIcdSASvBkzZrBu3bq8+j0RItmgDkSCN27cOEpK9AsmRRqaOhAJ3ty5c3nvvfeSjiESHHUgErzmzZvTrFmzpGOIBEcdiARv+vTpmoGIxEAdiARv0qRJzJo1K+kYIsFRByLBmzdvHqWlpUnHEAmOOhAJXtOmTWnatGnSMUSCow5EgldcXMz69es1AxFpYOpAJHjFxcU8/fTTSccQCU5BZWVl0hly3htvvKEnSUTyVufOnQtqul4FREREMqIlLBERyYgKiIiIZEQFREREMqLTePOAmTUDpgPHAjuAge5eXu0+s4HWwF5gl7tflPWgh8jMDgN+B5wG7AGudPf3q9x+FXAN8Akwwt3/mEjQBlSHY34IOJv06w7Q2923Zz1oDMzsfwO/cvdu1a7vBQwn/TpPcfeJCcSLxQGO+SfAlcD+/9fXeBZ+i5oKSH64DnjH3e8ys37AMOC/qt3nJOCr7t6Yz6r4HnCEu3/DzM4EHgB6A5hZW+BG4H8BRwCvmNk8d9+TWNqGUesxRzoD33b3DxNJFxMzuxW4FNhZ7fqmwINAl+i2V83sGXfflP2UDau2Y450Bi5z9zeymUlLWPnhbOC56PKzwAVVbzSz44BjgDlm9oqZfTfL+RrKP47T3V8nXSz2OwN41d33RO/A3wf+PfsRG1ytxxx1JycBE8zsVTO7PJmIsVgB9Knh+k7A++6+1d0/Bl4Bzs1qsvjUdsyQLiC3Rf9/b8tWIBWQwJjZFWa2tOofoAWwf9liR/R1VZ8n/c71e6T/gT5oZsdmLXTDOZrPjhPgUzP7XC231fQ8NEYHOuYjgd8AA4ALgcFmFkLRxN1nkV5urS7U1/lAxwxQAlwLnA+cna03gVrCCoy7TwYmV73OzJ4E9v9S8KOAbdUethF4xN0/ATab2ZuAAZtjjtvQPuKz4wQ4LDqmmm6r6XlojA50zCngIXdPAZjZS6RnJW9nN2JWhfo618rMCoAx+2dbZvYn4OtA7DM+dSD54VWgZ3T5ImBhtdsvAJ4AMLNC4FSgMf742n8cZzQPeKfKbYuBc8zsCDNrQXqpY2n2Iza4Ax3zyaRnAE2i2cDZwF+yHzGrSoGTzOzfzOzzpJevXks4U9yOBpaaWWFUTM4HsjILUQeSH8YBj5rZK8DHQH8AM/s1MNPdnzWzb5vZ68A+4PZGOnR9CuhuZouAAmCQmd1Eek38GTMbS7p4Hgb83N13J5i1oRzsmKcBr5Ne+njM3d9NMGtszKw/UOjuE6Lj/zPp13mKu69LNl08qh3z7cB80mfivejuc7ORQT/KRESNXOyRAAAAMklEQVREMqIlLBERyYgKiIiIZEQFREREMqICIiIiGVEBERGRjKiAiIhIRlRAREQkI/8faPTofeSs+JYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_bait.trace, varnames=['alpha', 'logit_psi', 'beta_psi_bait'], rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [alpha, beta_psi_c, beta_psi_bait, logit_psi]\n",
      "Sampling 2 chains: 100%|██████████| 15000/15000 [00:41<00:00, 361.52draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_bait_camper:\n",
    "    base_rate = data.fish_caught.sum() / data.hours.sum()\n",
    "\n",
    "    logit_psi = pm.Normal('logit_psi', mu=0, sd=1.5)\n",
    "    beta_psi_bait = pm.Normal('beta_psi_bait', mu=0, sd=1)\n",
    "    beta_psi_c = pm.Normal('beta_psi_c', mu=0, sd=1)\n",
    "\n",
    "    alpha = pm.Exponential('alpha', lam=1/base_rate, shape=2)\n",
    "\n",
    "    psi_i = pm.math.invlogit(\n",
    "        logit_psi\n",
    "        + beta_psi_bait * data.livebait.values\n",
    "        + beta_psi_c * data.camper.values\n",
    "    )\n",
    "\n",
    "    fish = (\n",
    "        data.hours.values\n",
    "        * alpha[data.livebait.values]\n",
    "    )\n",
    "\n",
    "    fish_obs = pm.ZeroInflatedPoisson(\n",
    "        'fish_obs',\n",
    "        psi = psi_i,\n",
    "        theta = fish,\n",
    "        observed = data.fish_caught.values\n",
    "    )\n",
    "    \n",
    "    m_bait_camper.name = 'livebait & camper'\n",
    "    m_bait_camper.trace = pm.sample(5000, tune=2500, chains=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEcCAYAAADpzeJvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYVOWV7/FvyxAF2guDrSQZkAPRZSfMjAmnPcagoBGCJEBCLjAMSiNe8RKCmkFDvMUYIkbBOCIi0ArHaeMlChMxomIPih44Rkdw2gWCDYRrawBpChBDzx+1SXravhZVtavf+n2eh6erau9637Wrmlq13rWrq6C2thYREZHWOiLuAEREpG1SAhERkZQogYiISEqUQEREJCVKICIikhIlEBERScnfxB2ASKaY2ZHAXKAn8BFwpbuvMbPvAHcBG6NdbwbeABYAHYDL3P1tM+sLfM3df9nI+JcCo4GDQHvgJ+7+cgpxXg50BR4AbnL38WZWBZzq7vvq7Fca3TapBWOWAeXu/lwj27sD/+juC1sbbzPz9gcud/eR6RxXcpMqEAnZJUCNu58BXA3cF93eB/ixu/eP/lUAA0kmkPHAODMrAH4ITG9oYDMbCQwAvu7u/Ukmknlmdnyqwbr7Vncfn+r9W+lc4GtZmksCpQpEQvZFYBGAu7uZFUe39wG+bGYTgOXAvwA1JKuPDsAeYBTw27oVQD2XARPd/UA0/vtmdpq7fxi9++8S/fsm8GPgLKAdcLe7Px5VN9OBHcAnwOtm1oNk1XBGNMfM6LZtwJi6k5vZ1VGMtdF97m0oyKgi+BfgY5KVWDkwBZgEdDSzZcD7wL1AAfAhcBHwZeCX0f1eAPq7+znRmP8O/BToBVxJsvqqBb5Tb+65wBeix3S6u89r5LGUNkoViITsLeBbZlZgZmcAnzezdsBikhXJ2UAhcDnJF8kTSSaGB0m+GP6nmc00sx83MPbngHV1b3D3D+tcfcndzwTOAP6Xu/cFzgF+YmbHATOAf3L380i+gDdkhrv3A6pIVlMAmNkXgRFAX5KJ6dtmZk08DicB341i+bG7/5lkEnnU3RcAs0gu7/UHniWZ8ACOcvez3P1W4CgzO8nMPgsc7+5vAqcA34yO7b+Ab9SJ8WiSj+9wYBDw5ybikzZKFYiEbA5QDCwFXgXecPc/m9kcd98JYGbPAN9194PAhOi2G0hWB5NJJppbzOwUd19dZ+z1QDdg16EbzOwbwNvRVY9+/j3Qx8xejq63B3oAJ9YZ71WS79Tr+tjdX48uLyO5XLYiut6bZFJ4MbreGTi5zpz1rXT3T4BPzGxvA9uLgfujHNQeWFPvGABmAxcC+0n2lQC2Aw+bWQ1wKvDaoZ3dfXdU4T0IHAPMbyQ2acNUgUjISoAXo3fIjwProt7G22b2d9E+XyfZQAfAzE4AzN2XAh1JvnOuBTrVG3sO8FMz+5vofqcAD/HXd9oHo5/vAkuid/fnAr8B1gKb6iyplTQQ+2fM7LTo8lnAqjrbHHgHOCcat4y/Jq6GNPQH7w7y1///DlwYjfVj4N/rHQMkl76+RbIye9TMjgVuBUYCFwN7SS6BARBVKn3c/Tskl/HuPPRYSTiUQCRka4AJZvYa8DOSPYtaki94T5lZBckkMavOfSYDt0eX7wd+D3wW+M+6A7t7OfA68IqZ/QfJd+Wj3X17vRgWAjVmtpRkoqp1990kl8oeMbMXSVYT9e0Hro7GPoFkcjo093+SrD5eMbP/T7L62NTyhwWAlcCw6GSAK6JYXiG5tPWpZOTuNdFjUBnF/xHJyuk1khXeXpLLeodsBbpGPZbFwF1RFSQBKdBf4xURkVSoAhERkZQogYiISEqUQEREJCVKICIikhIlEBERSYnOy26BN954Q6eqiUje6tOnT0FDt+s03pYJ8kGqrKykuLi4+R3buFmzZrFlyxZuuummuEPJuHx5TnWcWdVg8gBVIJIHLrnkEiorK+MOQyQ46oGIiEhKVIFI8Pr3708ikWD58uVxhyISFCUQCV5paSmbN2+OOwyR4CiBSPBKS0vVAxHJAPVAJHgHDhzgwIEDcYchEhxVIBK8AQMGqAcikgFKIBK8iy++mE2bWvt1GSLSHCUQCd7o0aPVAxHJAPVAJHiJRIK9exv6KnARORyqQCR4gwcPVg9EJAOUQCR4V1xxhXogIhmgBCLBGzFihHogIhmgHogEb9euXezevTvuMESCowpEgjds2DD1QEQyQAlEgnfNNdfwxz/+Me4wYrfhwwTjHl7Buuo99CzqxOwxJXTv0jHusKQNUwKR4A0fPjzjPZCxc5ezxKszOkfLrWt2jzXbazh76pIsxHJ4zrEi5o49Pe4wpBH6RsKWCfJBqv9tZwPvqWD1tpoYIxKRTDjlxEKe/1G/VO+ubySU5h3GL1hOy6fvA2nqK1AH3F3B2uoaDtbCEQXQq6iQxRPb5nOeI1/1mnG5fpw5cRaWmZWZ2aBW3qermd0fXT7bzP4hTbFMM7Pu6RhLcsO1115LaWlp3GHEbvaYEnoVFdKuoIBeRYXMHlMSd0jSxrXZCsTdtwLjo6sXAeXA22kYd8LhjiG5ZciQIXn9OZD6zfMl1/VX81zSIqcSiJm1B+YCPYF2wN3u/piZnQ78K7Ab2A7sA24hmTSuBAYBXzGz/3L3DQ2MWwp8GzgaOB64zd2fNLOfA+eQfByedPdfmtnLwOXu/m4mj1WyZ+vWrVRXV2dsKSC3GujQVBNdzXNJp5xKIMBlQLW7jzazo4E/mNmLwAPABe7+TvSi//lDd3D3N8zsOaC8oeRRRydgAFAELDezZ4B/BvoDW4DSxu5YVVUV5B/j27dvH/2mPM/6nXnyZUvP/i7uCKSFlng1PSY193w1f7ZZGFI7zpOOa88Dw7od9uxNvfHKtQRSDLwA4O67zey/gF7A59z9nWifpcDIFMaucPeDwDYz20EykfwzMAXoCixq7I49evRIYbrcV1lZScWkgXGHkXHPPfccGzZs4NJLL407lIxrqOkaUvP8kFxvLqdLrh9nTjTR66gEzgKIKpC/B94HNprZF6N9zmjgfgdp/lj6ROOeCBwD/An4PvBPJJexSs3spMM9AMk9gwYN4qyzzoo7jNioeS6ZkmsVyIPALDN7BegA3Oru281sPDDHzGqAj4H6f1r1/wFTzOx9d2+sW9o1Wg47Fhjv7vvN7E/A68Be4HmgqSUwaaM2btzIli1bcvqdXCZ179KxzVcckptyIoG4e2mdq2Ma2OV0YIi7V5vZ7cDH7l5FVI24+0xgZjPTVLj7pHrz3gbcVm+//i2PXNqCCy64IG8+ByKSTTmRQFpgG/B8VIHsouEkQ/S5kC82sOmxDMYmOW7y5MmsX78+7jBEgtMmEoi7PwE80YL9xje3j+Sf8847L68/ByKSKbnWRBdJu3Xr1rFx48a4wxAJTpuoQEQOx0UXXUQikWDgwPBPWRbJJiUQCd6tt96qHohIBiiBSPD69eunHohIBqgHIsFzd95///24wxAJjioQCd5ll11GIpFg8ODBcYciEhQlEAneHXfcQVVVVdxhiARHCUSCd+aZZ9K5c+e4wxAJjnogErxVq1axZs2auMMQCY4qEAneVVddRSKRYOjQoXGHIhIUJRAJ3tSpU3UWlkgGKIFI8EpKSigsLIw7DJHgqAciwXvrrbf0QUKRDFAFIsGbMGECiUSC4cOHxx2KSFCUQCR406ZNY926dXGHIRIcJRAJ3mmnncaRRx4ZdxgiwVEPRIK3YsUKVq5cGXcYIsFRBSLBu/7660kkEvzgBz+IOxSRoCiBSPDuu+8+9UBEMkAJRILXu3dv2rVrF3cYIsFRD0SCt2zZMt588824wxAJjioQCd6NN95IIpFg1KhRcYciEhQlEAnezJkzWbt2bdxhiARHCUSCZ2YcPHgw7jBEgqMeiASvoqKCFStWxB2GSHBUgUjwbr75ZhKJBBdeeGHcoYgERQlEgjdnzhzee++9uMMQCY4SiASvZ8+e7N+/P+4wRIKjHogE74UXXmDZsmVxhyESHCUQCd7tt9/OzJkz4w5DJDhawpLgzZs3jzVr1sQdhkhwlEAkeN26daOmpibuMESCoyUsCd5zzz3H0qVL4w5DJDhKIBK8KVOm8NBDD8UdhkhwtIQlwSsvL2f16tVxhyESHCUQCV7Xrl3ZsWNH1ubb8GGCcQ+vYF31HnoWdWL2mBK6d+mYtflFskUJRNqUsXOXs8SrU7x39r+VcM32Gs6euiTLs7bsOM+xIuaOPT3DsUjICmpra+OOoS1I+UEaeE8Fq7fpDCCRkJxyYiHP/6hfxueprKykuLg44/M0o6CxDapAMiwbv2SpypFfzoz74IMPWL16NWeeeWZW5htwdwVrq2s4WAtHFECvokIWT8zO70G+PKf5cpy5rsVnYZlZmZkNauG+3c1sSOphtWiOcjP7TBPbt7Z2rGzELdl3/PHH07lz56zNN3tMCb2KCmlXUECvokJmjynJ2twi2ZSp03jPBb6WobEBcPeR7v5xmsfKeNySfU899RSLFy/O2nzdu3Rk9pgSehZ1Yl31HsY9vIINHyayNr9ItrS4B2JmZcDfAp1ILn2Nc/f3zOxqYBTJPkE58K/AO0BH4CpgF3AzyWRVCIxy9wbPqTSzW4BTgROAzsDV7v6Kmc0FvgB0AKa7+zwzqwJOdfd9jYy1A/g90A14GxgPfB6YARwFfBaY7O5PR2N9CXjjUNzuvqDOcEE2itrqMsDhNdKlvrbYTG+rv7utlSPHmbYeyDJ3n2Jmg4E7zWwyMALoG21fTPJFewrJF/cFZjYeGO3um83sRuD7wM+bmCPh7uea2ZeAR82sL3A2cAbJF/KBLYy1A/Av7r7ezH4DDAESwK/c/WUzOxO4FXg62v/PdeOuO1BVVRV79+5t4bSfdvkzG1m/80DK98+s7J+ZJLlliVfTY9Lv4g4jBfH87p50XHseGNYtK3Pt27ePysrKrMzVmKYSWGsTyH9EP5cBU4HewEnAi9HtnYGT691nE3CvmdWQrABebWaOlwDc/R0z6+ruu81sAvAgcAwwv4WxbnD39XXiNeBZYLKZjSOZjNq3ZKAePXq0cMqGVcT/DqJBOfLuJiuyfaxxNdLz5TnVceaG1vZADtW5ZwGrACe5XHWOu/cHykguFx2sM/YsYKy7lwKbaaIcivQBMLPewCYz+yzQx92/A3yTZOXTksT3d9F9IVkhrQJ+Bjzi7hcASxqIpW7cEojHHnuMRYsWZXVONdIlH7S2AjnDzF4i+e79omh56EXgFTM7ElhOsuJYCfzEzP5AsmJYamZ7gG3A55qZ48vRmJ2AS4CtQFczW0Zymekud//EzJqL9UOSlc/fkVx6W2RmxwF3mdkNwB+B4+vd5y9xu3t5Cx4PaQNmzJhBIpFg4sSJWZuze5eOWTt1VyQuOfVBwqiJvtXdH4g7lnpy50FKo1wvj9MlkUjw7rvv8pWvfCXuUDIuX55THWdW5dYHCc3sKZJndNW1C3izleMMBRp6Wznd3X+bYngSmI4dO9KhQ4e4wxAJTiwJxN2Hp2mcBcCCZneUvDZ//nw2bdqUC+/kRIKihrEE76GHHuLJJ5+MOwyR4OhvYUnwFi9eHPu59CIhUgUiwWvfvj3t27foIz8i0gqqQCR4ZWVlbN68WT0QkTRTBSLBKysr4+mnn25+RxFpFVUgEryXX35ZPRCRDFAFIiIiKVECkeDNmjWLxx9/PO4wRIKjBCLBi+OPKYrkA/VAJHgvvPCCeiAiGaAKREREUqIEIsG7//77+bd/+7e4wxAJjpawJHgLFy6kpqYm7jBEgqMEIsFbtGiReiAiGaAlLBERSYkSiARv+vTpzJs3L+4wRIKjJSwJ3osvvqgeiEgGKIFI8BYsWKAeiEgGaAlLRERSogQiwbvrrruYM2dO3GGIBEdLWBK81157jd27d8cdhkhwlEAkeE8++aR6ICIZoCUsERFJiRKIBG/KlCnMmjUr7jBEgqMlLAneW2+9xUcffRR3GCLBUQKR4JWXl6sHIpIBWsISEZGUKIFI8H72s58xY8aMuMMQCY6WsCR47s6uXbviDkMkOEogErz58+erByKSAVrCEhGRlCiBSPBuuukmfv3rX8cdhkhwtIQlwdu4cSM7d+6MOwyR4CiBSPDmzp2rHohIBmgJS0REUqIKRIJ3ww038MEHHzT497A2fJhg3MMrWFe9h55FnZg9poTuXTrGEKVI26MEIrEbO3c5S7w6gzP0hS7QY9LvmtxrzfYazp66JINxfNo5VsTcsadndU6RdFECyVED76lg9baaLMy0LgtzSGOWeHWzia31cus5PeXEQp7/Ub+4w5AMUALJUdn4D1dZWUlxcXHG58kFjR3rgLsrWFtdw8FaOKIAehUVsnhi232xy6fnVOKXkSa6mZWZ2aAW7tvdzIZkIo46c5Sb2WcyOYfkruuuu44777yzwW2zx5TQq6iQdgUF9CoqZPaYkixHJ9J25cJZWOcCX8vkBO4+0t0/zuQckrv27t3L/v37G9zWvUtHZo8poWdRJ9ZV72HcwyvY8GEiyxGKtE2ZXMIab2bXR3OMc/f3zOxqYBRQC5QD/wpMAjqa2TJgF3AzycRWCIxy99UNDW5mtwCnAicAnYGr3f0VM5sLfAHoAEx393lmVgWc6u77GhlrHHAF0A5Y4O43p+H4pRUy2kg/ejAc3XwTHdRIF2mNTCaQZe4+xcwGA3ea2WRgBNA32r4Y+D0wheSL+wIzGw+MdvfNZnYj8H3g503MkXD3c83sS8CjZtYXOBs4g2SSGthckGZ2Askk9g/APuAXZlbo7n/pYFdVVbF3797WHX0aXP7MRtbvPJDhWXKr4ZqP0t9Iz53n9KTj2vPAsG5pH3ffvn158eHQXDjOpnpqmUwg/xH9XAZMBXoDJwEvRrd3Bk6ud59NwL1mVgN8Hni1mTleAnD3d8ysq7vvNrMJwIPAMcD8FsTZE1jl7ocyxKT6O/To0aMFw6RfRYabofnScJ0wYQJ/+tOfeOSRRxrcHlIjPV+eUx1nbshkD+RQTX4WsApw4B3gHHfvD5QBbwMH68QxCxjr7qXAZqCgmTn6AJhZb2CTmX0W6OPu3wG+SbLyaS5JrgVONbMjo7GeMLPPt/AYJQBqpIukJpMVyBlm9hLJpaSL3H29mb0IvBK9WC8nWXGsBH5iZn8gWTEsNbM9wDbgc83M8eVozE7AJcBWoGvUT/kzcJe7f2JmjQ7g7tVm9kugwsxqgYXuvukwjltyzLRp05pcBujepWObrThE4pSRBBJVEA3dPpXkclZdbwKHXuHLWzlVubs/UO+2yxuYt0dTg7h7GcmKSEREWijnP0hoZk8Bf1vv5l0kE09rxhkKTGxg03R3/22K4UkbcOWVV7Jjxw4effTRuEMRCUrOJxB3H56mcRYAC9IxlrQtHTp0IJHQZztE0i3nE4jI4brrrrtiPxVSJES58El0ERFpg1SBSPAuvfRSdu7cyW9+85u4QxEJihKIBK9Lly7U1tbGHYZIcJRAJHi/+MUv1AMRyQD1QEREJCWqQCR4Y8eOZefOnfz2t/q4j0g6KYFI8Lp160aHDh3iDkMkOEogErzbbrtNPRCRDFAPREREUqIKRII3evRodu3axcKFC+MORSQoSiASPDOjujpDX5crkseUQCR4P/3pT9UDEckA9UBERCQlqkAkeCNHjuSjjz7i2WefjTsUkaAogUjwTjvtNLZv3x53GCLBUQKR4E2aNEk9EJEMUA9ERERSogpEgvfd736X3bt38/zzz8cdikhQlEAkeF/96lfZtm1b3GGIBEcJRIJ33XXXqQcikgHqgYiISEpUgUjwhg4dSk1NDS+99FLcoYgERQlEgvf1r39dPRCRDFACkeD98Ic/VA9EJAPUAxERkZSoApHgnX/++dTU1LB06dK4QxEJihKIBG/IkCFs3bo17jBEgqMEIsEbP368eiAiGaAeiIiIpEQViATvvPPOY8+ePbz22mtxhyISFCUQCd6IESPYsmVL3GGIBEcJRIJ3ySWXqAcikgHqgYiISEpUgUjw+vfvTyKRYPny5XGHIhIUJRAJXmlpKZs3b077uBs+TDDu4RWsq95Dz6JOzB5TQvcuHdM+j0iuUgKR4JWWlqa1BzJ27nKWePX/uG3N9hrOnrrkL9fPsSLmjj09bXOK5CIlEAnegQMHOHDgQKvvN/CeClZvq0lpziVeTY9Jv2t0+yknFvL8j/qlNLZIrlACkeANGDAgpR5Icy/wA+6uYG11DQdr4YgC6FVUyOKJSgqSP3I2gZjZVnfv2si2/sDl7j6ylWOWArcB04BHgEeBDsBmYCxwEXAdMMXdH0g5eMkpF198MZs2bUr7uLPHlHyqByKST/LxNN5H3f1u4Kbo8lnAm8Bl7n4fUBZncJJ+o0ePZujQoXGHIRKc2CsQM+sN3A20A44HrnD3ZXW2vwy8C5wKFAAjok0nm9ki4ARgobvfYmb9gJtJJsZCYJS7r25k6r7AHdHlRdHle9J4aJIjEokEe/fuTdt4aqKLJMWeQIAvAde6+0ozG0VyKWlZvX2WufvlZjYeuBF4CjgK+DbJxLMBuCUaa7S7bzazG4HvAz9vZN5jgF3R5d3AsY0FWFVVldYXoFyxb9++vPiE9pgxYzh48CDz5s1r8X0uf2Yj63e2vvF+SHNN9JOOa88Dw7qlPH5j8uU51XFmT3FxcaPbciGBbAJ+amZ7gaOBjxrY56Xo5zJgWHR5lbvvBzCzT+qMda+Z1QCfB15tYt6PovkOzbuzsR179OjRogNpayorK5v85QjFtddey6ZNm1p1rBUt2DcXm+j58pzqOHNDLvRA7gVudvcxwEqSy1T19Yl+fg14J7pc28B+s4Cx7l5KsjHe0FiHvAoMji6fD+jr6gI1YsQIzj///LSPO3tMCb2KCmlXUECvokI10SXv5EIFMh943Mx2AH8k2Qepr9TMJgJ7gAuAv29irKVmtgfYBnyuiXlvBx42s0uAD4BRKcYvOW7Xrl3s3r077eN279Ix9opDJE6xJ5DojKi7G7i97im8N7j7u3Wuvxz9+x/7uvvEVsy7DRjUynClDRo2bJj+FpZIBuTCEla2jYqqmU8xs6uA0uyGI5l2zTXXMHr06LjDEAlO7BVIc9y9fxrHKqOJz3lEnwO5L13zSW4YPnx47GeyiIQoHysQyTMffPABO3bsiDsMkeDkfAUicri+973vqQcikgFKIBK8a6+9lo0bN8YdhkhwlEAkeEOGDFEPRCQD1AOR4G3dupXq6urmdxSRVlEFIsEbOXKkeiAiGaAEIsGbNGkSGzZsiDsMkeAogUjwBg0apB6ISAaoByLB27hxI1u2bIk7DJHgqAKR4F1wwQXqgYhkgBKIBG/y5MmsX78+7jBEgqMEIsE777zz1AMRyQD1QCR469at0yfRRTJAFYgE76KLLiKRSDBw4MC4QxEJihKIBO/WW29VD0QkA5RAJHj9+vVTD0QkA9QDkeC5O++//37cYYgERxWIBO+yyy4jkUgwePDguEMRCYoSiATvjjvuoKqqKu4wRIKjBCLBO/PMM+ncuXPcYYgERz0QCd6qVatYs2ZN3GGIBEcViATvqquuIpFIMHTo0LhDEQmKEogEb+rUqToLSyQDlEAkeCUlJRQWFsYdhkhw1AOR4L311lv6IKFIBqgCkeBNmDCBRCLB8OHD4w5FJChKIBK8adOmsW7durjDEAmOEogE77TTTuPII4+MOwyR4KgHIsFbsWIFK1eujDsMkeCoApHgXX/99SQSCX7wgx/EHYpIUJRAJHj33XefeiAiGaAEIsHr3bs37dq1izsMkeCoByLBW7ZsGW+++WbcYYgERxWIBO/GG28kkUgwatSouEMRCYoSiARv5syZrF27Nu4wRIKjBCLBMzMOHjwYdxjSjA0fJhj38ArWVe+hZ1EnZo8poXuXjnGHJU1QApHgVVRUsH79eoqLi+MORRoxdu5ylnj1X66v2V7D2VOXcI4VMXfs6TFGJk1RApHg3XzzzSQSCS688MK4Q5HIwHsqWL2tptn9lng1PSb97n/cdsqJhUwfdEKmQpNWUAKR4M2ZM4f33nsv7jCkjud/1O9Ttw24u4K11TUcrIUjCqBXUSGLJ356P0B/XTlH5NVpvGZWamYbzGxidH2CmU2JLn/PzN49dF3C0bNnT7p16xZ3GNKM2WNK6FVUSLuCAnoVFTJ7TEncIUkz8iqBRB4FZpjZ/wWuPHSjuz8BKHkE6IUXXmDZsmVxhyESnHxdwjoKeBhYDJwacyySYbfffjuJRIJx48bFHYo0Qk30tqmgtrY27hiyxsxKgVPdfVJLrh9SVVVVu3fv3uwGmwX79u3jqKOOijuMjNuyZQv79++nR48ecYeScW3lOb38mY2s33kgpfuedFx7pn2jqE0c5+HKheezuLi4oLFt+VqBtEqoLzyVlZV5cWprcXFx3hxrWznOigZibG0TvS0c5+HK9ePMxx6I5JnnnnuOpUuXxh2GNENN9LZHFYgEb8qUKSQSCS699NK4Q5EmdO/SsdGKQ3JTXicQdy+LOwbJvPLyclavXh13GCLBycclrFGHPgdSl5l9D5jUwP7SxnXt2pWioqK4wxAJTl5VIFHFUdbItieAJ7IZj2THwoUL2bhxY043I0XaonysQCTP/OpXv6KsrCzuMESCk1cViOSnJ554Qj0QkQxQBSLBO/744+ncuXPcYYgERwlEgvfUU0+xePHiuMMQCY4SiATv3nvvZf78+XGHIRIc9UAkeM888wzuHncYIsFRBSLBO/bYYzn66KPjDkMkOKpAJHiPPfYYmzZt0udARNJMFYgEb8aMGZSXl8cdhkhwVIFI8J599lnefffduMMQCY4qEAlex44d6dChQ9xhiARHFYgEb/78+eqBiGSAKhAJ3kMPPcSTTz4ZdxgiwVEFIsFbvHgxlZWVcYchEhxVIBK89u3b0759+7jDEAmOKhAJXllZGZs3b1YPRCTNVIFI8MrKynj66afjDkMkOAW1tbVxx5Dz3njjDT1IIpK3+vTpU9DQ7UogIiKSEi1hiYhISpRAREQkJUogIiKSEp3Gm0fMrAMwHzgB2A2Mcffqevs8AxwPHAD2uvuek4rAAAADGUlEQVT5WQ80RWZ2BHA/8I/AfuBid3+vzvZLgMuAT4Db3f3fYwn0MLXgOKcDfUk+xwDD3H1X1gNNEzP7P8Av3b1/vduHADeRfD7nuPusGMJLqyaO9UfAxcCh/6+XeQ58S5oSSH65Aljp7reY2UhgMvDDevucDHzJ3dvi2RXfBo5y96+a2RnAr4BhAGbWFbgG+N/AUcArZrbY3ffHFm3qGj3OSB/gG+7+QSzRpZGZ/Ri4ANhT7/b2wD1ASbTtVTNb4O7bsh9lejR2rJE+wIXu/kZ2o2qalrDyS1/guejyIuC8uhvN7ETgOGChmb1iZt/KcnyH6y/H5+6vk0wWh5wOvOru+6N34+8B/5D9ENOi0eOMqpOTgQfN7FUzuyieENNmLTC8gduLgffcfYe7fwy8Apyd1cjSr7FjhWQCuSH6f3lDFmNqkhJIoMxsnJmtqvsPOBY4tJSxO7pe12dIvpv9Nslf5HvM7ISsBX34juGvxwfwZzP7m0a2NXT8bUVTx9kJ+DUwGhgEjDeztpoocfcnSS6n1hfS8wk0eawA5cDlwLlA31x5c6clrEC5+2xgdt3bzOwp4NCXgx8N7Kx3t63AA+7+CbDdzN4EDNie4XDT5SP+enwAR0TH0tC2ho6/rWjqOBPAdHdPAJjZSyR7JW9nN8SMC+n5bJKZFQDTDvWxzOx3wJeB2Ht4qkDyy6vA4Ojy+cDSetvPAx4HMLNCoDfQlv6M7V+OL+oNrKyzbTlwlpkdZWbHklwCWZX9ENOiqeM8hWQ/oF3UJ+gL/CH7IWZcJXCymf2tmX2G5PLVazHHlCnHAKvMrDBKJucCOdELUQWSX2YAD5vZK8DHwCgAM7sTeMLdF5nZN8zsdeAgcGMba8T+FhhgZsuAAmCsmU0kuVa+wMzuJZk0jwB+4u77Yoz1cDR3nPOA10kuhzzi7u/EGGtamdkooNDdH4yO+fckn8857r4p3ujSq96x3ggsIXnW3Yvu/my80SXpT5mIiEhKtIQlIiIpUQIREZGUKIGIiEhKlEBERCQlSiAiIpISJRAREUmJEoiIiKTkvwEqxeIH18Vw9QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_bait_camper.trace,\n",
    "              varnames=['logit_psi', 'beta_psi_bait', 'beta_psi_c', 'alpha'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [beta_fish_p, alpha, beta_psi_p, beta_psi_bait, logit_psi]\n",
      "Sampling 2 chains: 100%|██████████| 15000/15000 [01:04<00:00, 234.04draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_bait_persons:\n",
    "    base_rate = data.fish_caught.sum() / data.hours.sum()\n",
    "\n",
    "    logit_psi = pm.Normal('logit_psi', mu=0, sd=1.5)\n",
    "    beta_psi_bait = pm.Normal('beta_psi_bait', mu=0, sd=1)\n",
    "    beta_psi_p = pm.Normal('beta_psi_p', mu=0, sd=2)\n",
    "\n",
    "    alpha = pm.Exponential('alpha', lam=1/base_rate, shape=2)\n",
    "    beta_fish_p = pm.Normal('beta_fish_p', mu=0, sd=2)\n",
    "\n",
    "    psi_i = pm.math.invlogit(\n",
    "        logit_psi\n",
    "        + beta_psi_bait * data.livebait.values\n",
    "        + beta_psi_p * data.persons.values\n",
    "    )\n",
    "\n",
    "    fish = (\n",
    "        data.hours.values\n",
    "        * alpha[data.livebait.values]\n",
    "        * pm.math.exp(beta_fish_p * data.persons.values)\n",
    "    )\n",
    "\n",
    "    fish_obs = pm.ZeroInflatedPoisson(\n",
    "        'fish_obs',\n",
    "        psi = psi_i,\n",
    "        theta = fish,\n",
    "        observed = data.fish_caught.values\n",
    "    )\n",
    "    \n",
    "    m_bait_persons.name = 'livebait & persons'\n",
    "    m_bait_persons.trace = pm.sample(5000, tune=2500, chains=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEcCAYAAADpzeJvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X14VPWZ//F3oBSBqKUQpSqYn1Rv09Krdvmla61IsMJSukqLbWVZLA8+oPiEVC1afG4tXRSF9QkRQeHnxqoosBVWRExRbGGtVrDhBsGE8BwpIGEAkeT3xzlomiYkMyY5czKf13V5kZk58z2fHIe55z73zJBVVVWFiIhIslpFHUBEROJJBURERFKiAiIiIilRARERkZSogIiISEpUQEREJCVfiDqASFMxs7bADOAU4CPgKndfa2Y/Au4FysJNbwfeAuYB7YBR7v6umZ0NfNfdf1vH+pcDQ4FKoA3wS3d/LYWcVwBdgEeB29x9tJmVAKe7+/5q2w0PrxvXgDVnAoXuvrCO27sB33T3+cnmrWe/BcAV7j64MdeV9KQORFqyy4AKdz8TuAZ4MLy+J3CTuxeE/xUB/QgKyGjgEjPLAq4DJte2sJkNBvoC33P3AoJCMsvMOqca1t23uvvoVO+fpHOB7zbTvqSFUgciLdnXgAUA7u5mlhde3xP4lpmNAZYDvwAqCLqPdsBeYAjwQvUOoIZRwFh3Pxiu/4GZneHuO8JX/53C/34A3AT0AloDk9z92bC7mQzsBD4B/mhmuQRdw5nhPqaG120DhlXfuZldE2asCu8zpbaQYUfwC+Bjgk6sEJgAjAPam9ky4ANgCpAF7ABGAt8Cfhve7xWgwN37hGv+N3Ar0B24iqD7qgJ+VGPfM4Cvhsd0srvPquNYSkypA5GW7B3gX80sy8zOBE40s9bAIoKO5BwgG7iC4EnyeILC8BjBk+FfzGyqmd1Uy9onAOurX+HuO6pdfNXdzwLOBP6Pu58N9AF+aWZfAh4B/s3dzyN4Aq/NI+7eGygh6KYAMLOvARcBZxMUph+amR3hOJwMXBhmucndDxEUkafdfR4wjeD0XgHwEkHBAzjK3Xu5+53AUWZ2spl9Bejs7m8DpwE/CH+3vwL/Ui3j0QTHdxDQHzh0hHwSU+pApCV7AsgDlgJvAG+5+yEze8LddwGY2VzgQnevBMaE191M0B2MJyg0d5jZae6+ptrapUBXYPfhK8zsX4B3w4se/vkNoKeZvRZebgPkAsdXW+8Nglfq1X3s7n8Mf15GcLpsRXi5B0FRWBxe7gicWm2fNa1090+AT8xsXy235wEPhzWoDbC2xu8AMB34GXCAYK4EsB140swqgNOBNw9v7O57wg7vMeAYYHYd2STG1IFIS5YPLA5fIT8LrA9nG++a2UnhNt8jGKADYGbHAebuS4H2BK+cq4AONdZ+ArjVzL4Q3u804HE+e6VdGf65GlgSvro/F/gdsA7YVO2UWn4t2b9oZmeEP/cCVlW7zYH3gD7hujP5rHDVprYvvKvks7//DvwsXOsm4L9r/A4QnPr6V4LO7GkzOxa4ExgMXArsIzgFBkDYqfR09x8RnMb7j8PHSloOFRBpydYCY8zsTeBugplFFcET3hwzKyIoEtOq3Wc88Kvw54eB/wG+Avyl+sLuXgj8EXjdzP5A8Kp8qLtvr5FhPlBhZksJClWVu+8hOFX2lJktJugmajoAXBOufRxBcTq8778QdB+vm9n/EnQfmxp+WABYCQwM3wxwZZjldYJTW/9QjNy9IjwGxWH+jwg6pzcJOrx9BKf1DtsKdAlnLIuAe8MuSFqQLH0br4iIpEIdiIiIpEQFREREUqICIiIiKVEBERGRlKiAiIhISvS+7AZ466239FY1EclYPXv2zKrter2Nt2Fa5EEqLi4mLy+v/g1bOB2HwLRp09iyZQu33XZb1FHSgh4Xn6q1eIA6EBEJXXbZZRQXF0cdQ2JEMxAREUmJOhARAaCgoIBEIsHy5cujjiIxoQIiIgAMHz6czZs3Rx1DYkQFRESAoIBoBiLJ0AxERAA4ePAgBw8ejDqGxIg6EBEBoG/fvpqBSFJUQEQEgEsvvZRNm5L9Z0Ukk6mAiAgAQ4cO1QxEkqIZiIgAkEgk2Levtn8yXaR26kBEBIABAwZoBiJJUQEREQCuvPJKzUAkKSogktE27Egw6sUyNj31AafkdGD6sHy6dWofdaxIXHTRRZqBSFJUQKTZjZixnCVeHnWMf7B2ewXnTFwSdYy/08dymDHi282yr927d7Nnz55m2Ze0DPo694ZJ+SD1u7+INdsqGjOLiACnHZ/Ny9f3brL19XXun0rvr3M3s5lAobsvTOI+XYDb3H20mZ0D7HL3dxshywPAJHff8HnXApr0Af556S8I9J1UxPvbK6gCWmVB95xsFo1N3/9nTWnOnDls3LiRa6+9NuooEhOxfRuvu29199HhxZHACY207pjGKh6S/qYPy6frsW1onZVFty+35+ChSrrf/BJ9JxWxYUci6njNatCgQfTt2zfqGBIjadGBHGZmbYAZwClAa4JO4Bkz+zbwELAH2A7sB+4ACoGrgP7AP5nZX2t78jez4cAPgaOBzsBd7v68mf0a6ENwHJ5399+a2WvAFe6+uil/V0m/WUhJtYKRTvOQ5pqDfPjhh+zcubPJ9yMtR1oVEGAUUO7uQ83saODPZrYYeBS42N3fC5/0Tzx8B3d/y8wWEpwCO1Ln0AHoC+QAy81sLvDvQAGwBRhe1x1LSko+1wesrphbRumudP2SuvVRB5B6LPFycsf9vvl2OK8Z99VETv5SGx4d2PVzrbF//369Kw2OeJo73QpIHvAKgLvvMbO/At2BE9z9vXCbpcDgFNYucvdKYJuZ7SQoJP8OTAC6AAvqumNubm4Ku6u24zSdM2gGEjh8HPpOKmJdeQWVVZk5D5k/fz5lZWWMHj26/o0zgP5+1C/dZiDFQC+AsAP5BvABUGZmXwu3ObOW+1VS/+/SM1z3eOAY4G/AT4B/IziNNdzMTv68v4DE1/Rh+XTPyaZ1Vhbdc7KZPiw/6kjN6vzzz6dPnz5Rx5AYSbcO5DFgmpm9DrQD7nT37WY2GnjCzCqAj4GaH5f9EzDBzD5w97p6zi7h6bBjgdHufsDM/gb8EdgHvAxoeJ7BunVqn1EdR01bt26lvLxcr7qlwdKigLj78GoXh9WyybeB89293Mx+BXzs7iWE3Yi7TwWm1rObIncfV2O/dwF31diuoOHJRVqOwYMH67uwJClpUUAaYBvwctiB7Kb2IoOZPQx8rZabnmnCbCItwrhx49iwQU24NFwsCoi7Pwc814DtNP0TSVH//v31riNJSroN0UUkImVlZWzZsiXqGBIjsehARKTpXXzxxZqBSFJUQEQEgPHjx1NaWhp1DIkRFRARAeC8887TDESSohmIiACwfv16ysrKoo4hMaIOREQAGDlyJIlEgn79+kUdRWJCBUREALjzzjs1A5GkqICICAC9e/fWDESSohmIiADg7nzwwQdRx5AYUQciIgCMGjWKRCLBgAEDoo4iMaECIiIA3HPPPZSUlEQdQ2JEBUREADjrrLPo2LFj1DEkRjQDEREAVq1axdq1a6OOITGiDkREALj66qtJJBJccMEFUUeRmFABEREAJk6cqHdhSVJUQEQEgPz8fLKzs6OOITGiGYiIAPDOO+/og4SSFHUgIgLAmDFjSCQSDBo0KOooEhMqICICwAMPPMD69eujjiExogIiIgCcccYZtG3bNuoYEiOagYgIACtWrGDlypVRx5AYUQciIgDceOONJBIJfvrTn0YdRWJCBUREAHjwwQc1A5GkqICICAA9evSgdevWUceQGNEMREQAWLZsGW+//XbUMSRG1IGICAC33HILiUSCIUOGRB1FYkIFREQAmDp1KuvWrYs6hsSICoiIAGBmVFZWRh1DYkQFREQAKCoqorS0lLy8PDbsSHDJkytYX76XU3I6MH1YPt06tY86oqQZFRCRNDZixnKWeHkz7rETt437/d9ds3Z7BedMXNKMGT7Tx3KYMeLbkexb6pdVVVUVdYY4aJEHqbi4mLy8vKhjpKTf/UWs2VYRdQyRz+W047N5+freUceoT1ZdNzS4AzGzmUChuy9swLbdgG+6+/yGrp8sMysEfubuH9dx+1Z375LMWkAXmji3NI7G/EsX50La2A4fi76TilhXXkFlFbTKgu452Swam/ZPdI1Kj4v6NdUprHOB04EmeyJ298GNvZaZNXlukXT1yiuvUFpaSofjTubgoUoqw76725fbM31YfrThJC01+BRW2IF8GehAUHgucff3zewaYAjBaZ5C4CHgPaA9cDWwG7id4EOL2cAQd19Txz7uIHgCPw7oCFzj7q+b2Qzgq0A7YLK7zzKzEuB0d99fx1o7gf8BugLvAqOBE4FHgKOArwDj3f3FcK2vA28dzu3u86otp1NYLVi6Hofmn3+kj3SYfaTr4yICn/8UVmiZu08wswHAf5jZeOAi4Ozw9kUET9oTCJ7c55nZaGCou282s1uAnwC/PsI+Eu5+rpl9HXjazM4GzgHOJHgi79fArO2AX7h7qZn9DjgfSAD3uftrZnYWcCfwYrj9oeq5qy9UUlLCvn37Grjb+Ni/f3+s/gW6K+aWUbrrYBOtru+ASidLvJzcGsP8aDTt4+LkL7Xh0YFdm3Qfn9eRimiyBeQP4Z/LgIlAD+BkYHF4fUfg1Br32QRMMbMKgg7gjXr28SqAu79nZl3cfY+ZjQEeA44BZjcw6wZ3L62W14CXgPFmdglBMWrTkIVyc3MbuMt4idsrrKImyhq349CUiouLuXbB9oyff4AeFw2R7HdhHe4pewGrACc4XdXH3QuAmQSniyqrrT0NGOHuw4HNHKEdCvUEMLMewCYz+wrQ091/BPyAoPNpSOE7KbwvBB3SKuBu4Cl3vxhYUkuW6rlFMsrChQtZunQp04fl0z0nm9ZZWXTPydb8Q+qUbAdyppm9SvDqfWR4emgx8LqZtQWWE3QcK4FfmtmfCTqGpWa2F9gGnFDPPr4VrtkBuAzYCnQxs2UEp5nudfdPzKy+rDsIOp+TCE69LTCzLwH3mtnNwEagc437fJrb3QsbcDxEWowJEyaQSCS4/PLLM7LjkOSl1edAwiH6Vnd/NOosNaTPQWpEatEDOg6BrVu3smbNGs4555yoo6QFPS4+1WhD9EZhZnMI3tFV3W4gqe+SNrMLgLG13DTZ3V9IMZ5IRurSpQs7d+6MOobESCQFxN0HNdI684B59W4oIvWaP38+ZWVletUtDaaBsYgAcN999zFz5syoY0iM6MsURQSA5557jjVrav2Mr0it1IGICACdO3emY8eOUceQGFEBEREA5syZw6JFi6KOITGiAiIiAEyZMoXZsxv6RQ8imoGISGju3Lm4e9QxJEbUgYgIAMceeyxHH3101DEkRtSBiAgAzzzzDJs2bdLnQKTB1IGICACPPPIIhYX6CjhpOHUgIgLASy+9xOrVq6OOITGiDkREAGjfvj3t2rWLOobEiDoQEQFg9uzZmoFIUtSBiAgAjz/+OM8//3zUMSRG1IGICACLFi2iuLg46hgSI+pARASANm3a0KZNm6hjSIyoAxERAGbOnMnmzZs1A5EGUwciIkBQQF588cWoY0iMqAMREQBee+01zUAkKepAREQkJSogIgLAtGnTePbZZ6OOITGiAiIiQPBligsWLIg6hsSIZiAiAsArr7yiGYgkRR2IiIikRAVERAB4+OGH+a//+q+oY0iM6BSWiAAwf/58Kioqoo4hMaICIiIALFiwQDMQSYpOYYmISEpUQEQEgMmTJzNr1qyoY8TGhh0J+k4qovvNL9F3UhEbdiSijtTsdApLRABYvHhxRs9ARsxYzhIvr3Ht+gbdd+32Cs6ZuCSp/fWxHGaM+HZS90k3KiAiAsC8efPSfgbS7/4i1mxrGUVuiZeTO+73zbKv047P5uXrezf6uiogIhIbTfEkWJfi4uIjfrV930lFrCuvoLIKWmVB95xsFo1tvnzpoElmIGY208z6N3DbbmZ2flPkqLaPQjP7YlPuQyTu7r33Xp544omoY8TG9GH5dM/JpnVWFt1zspk+LD/qSM0uHTqQc4HTgflNtQN3H9xUa4u0FG+++SZ79uyJOobESFZVVVWjL2pmM4EvAx0IitQl7v6+mV0DDAGqgELgIeA9oD1wNbAbuJ2gM8oGhrj7mjr2cQdB4TkO6Ahc4+6vm9kM4KtAO2Cyu88ysxLgdHffX8da64E/Ad2BVcCl7l5ZbZPGP0hpoL4WPVPoOHwmk49F7UP0phOjIXpWXTc0ZQeyzN0nmNkA4D/MbDxwEXB2ePsi4H+ACQRP7vPMbDQw1N03m9ktwE+AXx9hHwl3P9fMvg48bWZnA+cAZxI86fdrYNaTgH5hkfsd8ENgzuEbS0pK2LdvX0N/79jYv39/2g9Nm4OOw2fS8VhcMbeM0l0Ho47R6JpziH7yl9rw6MCuKd33SC8omrKA/CH8cxkwEegBnAwsDq/vCJxa4z6bgClmVgGcCLxRzz5eBXD398ysi7vvMbMxwGPAMcDsBmbd4O7vV8tr1W/Mzc1t4DLxksmvNqvTcQhMmDCB7du3M2nSpKij/J2iiP7faIhev6b8IOHh3qwXwWkhJzhd1cfdC4CZwLtAZbUc04AR7j4c2MwRWqdQTwAz6wFsMrOvAD3d/UfADwg6n4YUyRPNrEv483fDnCIZ5Z133mH16tVRx4gNDdGbtgM508xeJTiVNNLdS81sMfC6mbUFlhN0HCuBX5rZnwk6hqVmthfYBpxQzz6+Fa7ZAbgM2Ap0MbNlwCHgXnf/xMyOtAbAAeBBM+sK/JEmHOiLpKvCwsK0O32Vzrp1ap9xHUdNTVJAwg6itusnEpzOqu5tPjtlVJjkrgrd/dEa111Ry35z61lnv7v/OMl9i4hktHR4G+8Rmdkcgnd0VbeboPAks84FwNhabpqcYjSRFuXuu++mvLycKVOmRB1FYiLtC4i7D2qkdeYB8+q4+YXG2IdInLk7u3fvjjqGxEjaFxARaR6zZ8/WDESSoq9zFxGRlKiAiAgAt912G//5n/8ZdQyJEZ3CEhEAysrK2LVrV9QxJEZUQEQEgBkzZmgGIknRKSwREUmJCoiIAHDzzTen3fdgSXpTARERAHbs2KHPgUhSNAMREQAee+wxzUAkKepAREQkJepARASAG264gR07djBjxoyoo0hMqAMREQD27dvHgQMHoo4hMaIOREQAeOihhzQDkaSoAxERkZSoAxERAMaMGcPf/vY3nnrqqaijSEyoAxERkZSoAxERAB544AHNQCQp6kBERCQl6kBEBICrrrqKnTt38vTTT0cdRWJCHYiIANCuXTvatm0bdQyJEXUgIgLAvffeqxmIJEUdiIiIpEQdiIgAcPnll7Nr1y5+97vfRR1FYkIFREQA6NSpE1VVVVHHkBhRARERAH7zm99oBiJJ0QxERERSog5ERAAYMWIEu3bt4oUXXog6isSECoiIANC1a1c+afVF+k4qYn35Xk7J6cD0Yfl069Q+6miSpnQKS0QAuOuuuyj96oWs3V7Boaoq1m6v4JInV0QdS9KYCoiIfGrjRwf/7vL68r0RJZE4SNtTWGa21d271HFbAXCFuw9Ocs3hwF3AA8BTwNNAO2AzMAIYCdwATHD3R1MOLxJDQ4cOpeK9LWT3H0tlFbTKglNyOkQdS9JYJnYgT7v7JOC28OdewNvAKHd/EJgZZTiRqJgZBd86jdatsgBo3SqLuwf2iDiVpLPIOxAz6wFMAloDnYEr3X1ZtdtfA1YDpwNZwEXhTaea2QLgOGC+u99hZr2B2wkKYzYwxN3X1LHrs4F7wp8XhD/f34i/mkis3HrrrZzzm5c5FJ7GOlRZxa1zV7FobO+Ik0m6iryAAF8Hfu7uK81sCMGppGU1tlnm7leY2WjgFmAOcBTwQ4LCswG4I1xrqLtvNrNbgJ8Av65jv8cAu8Of9wDH1hWwpKSEffv2pfK7pbX9+/frg2PoOFS38aODVIYfRq+sgnXlFRl7bPS4COTl5dV5WzoUkE3ArWa2Dzga+KiWbV4N/1wGDAx/XuXuBwDM7JNqa00xswrgROCNI+z3o3B/h/e7q64Nc3NzG/SLxE1xcfERHxyZQschMHjwYPb8dStH/+DGT2cg3XOyM/bY6HFRv3SYgUwBbnf3YcBKgtNUNfUM//wu8F74c21f2jMNGOHuwwkG47WtddgbwIDw5+8DS5OLLdKynHHGGfyw1zfpnpNN66wsuudkM31YftSxJI2lQwcyG3jWzHYCGwnmIDUNN7OxwF7gYuAbR1hrqZntBbYBJxxhv78CnjSzy4APgSEp5hdpEcaNG6dX3ZKUyAtI+I6oSbVcX/0tvDe7++pql18L//u7bd19bBL73Qb0TzKuiIiE0uEUVnMbEnYz/8DMrgaGN28ckfRw4YUXct1110UdQ2Ik8g6kPu5e0IhrzeQIn/MIPwfyYGPtTyROvvOd77Bt27aoY0iMpH0BEZHmccMNN+htq5KUTDyFJSIijUAdiIgAcMEFF1BRUcGrr75a/8YiqICISOh73/ueZiCSFBUQEQHguuuu0wxEkqIZiIiIpEQdiIgA8P3vf5+KigqWLtW3+kjDqICICADnn38+W7dujTqGxIgKiIgAMHr0aM1AJCmagYiISErUgYgIAOeddx579+7lzTffjDqKxIQKiIgAcNFFF7Fly5aoY0iMqICICACXXXaZZiCSFM1AREQkJepARASAgoICEokEy5cvjzqKxIQKiIgAMHz4cDZv3hx1DIkRFRARAYICohmIJEMzEBEB4ODBgxw8eDDqGBIj6kBEBIC+fftqBiJJUQEREQAuvfRSNm3aFHUMiREVEBEBYOjQoZqBSFI0AxERABKJBPv27Ys6hsSIOhARAWDAgAGagUhSVEBEBIArr7xSMxBJigqIiADBlylqBiLJ0AxERADYvXs3e/bsiTqGxIg6EBEBYODAgZqBSFJUQEQEgGuvvZaNGzdGHUNiRAVERAAYNGhQnTOQDTsSXPLkCtaX7+WUnA5MH5ZPt07tmzmhpBvNQEQEgA8//JCdO3fWelv/yX9g7fYKDlVVsXZ7BZc8uaKZ00k6UgEREQB+/OMfM2bMmFpvS3x86O8ur91eQb/7i5ojlqSxjCogZjbczDaY2djw8hgzmxD+/GMzW334skim+fnPf87w4cNrve3U47JplRX83CoruPzy9b2bL5ykpYwqIKGngUfM7P8BVx2+0t2fA1Q8JGOdf/759OnTp9bb7h7Yg9ZhBWndKou7B/ZozmiSpjKxgAAcBTwJ/DrqICLpYuvWrZSXl9d628gnV3DwUBUABw9VcevcVc0ZTdJURr4Ly913Ai+b2fCGbF9SUtIiv2Ru//79+uQxOg6HDRs2jMrKSmbNmvUPt9U2A+k94WUeHdi1ueI1Oz0uAnl5eXXelpEFJFm5ublRR2gSxcXFR3xwZAodh8Bdd93Fhg0baj0Wpx63nXXlFVRWBTOQ7jnZLBrbsmcgelzUL1NPYYlIDf3796dXr1613jZ9WD7dc7JpnZVF95xspg/Lb+Z0ko7UgYgIAGVlZWzZsqXWV93dOrVv8R2HJC+jC4i7z4w6g0i6uPjii/VdWJKUTDyFNeTw50CqM7MfA+MiyCOSFsaPH8+oUaOijiExklEdSNhxzKzjtueA55ozj0g6Oe+88/SuI0lKJnYgIlKL9evXU1ZWFnUMiZGM6kBEpG4jR44kkUjQr1+/qKNITKiAiAgAd955J6WlpVHHkBhRARERAHr37q0ZiCRFMxARAcDd+eCDD6KOITGiDkREABg1ahSJRIIBAwZEHUViQgVERAC45557KCkpiTqGxIgKiIgAcNZZZ9GxY8eoY0iMaAYiIgCsWrWKtWvXRh1DYkQdiIgAcPXVV5NIJLjggguijiIxoQIiIgBMnDhR78KSpKiAiAgA+fn5ZGdnRx1DYkQzEBEB4J133tEHCSUp6kBEBIAxY8aQSCQYNGhQ1FEkJlRARASABx54gPXr10cdQ2JEBUREADjjjDNo27Zt1DEkRjQDEREAVqxYwcqVK6OOITGiDkREALjxxhtJJBL89Kc/jTqKxIQKiIgA8OCDD2oGIklRARERAHr06EHr1q2jjiExohmIiACwbNky3n777ahjSIyoAxERAG655RYSiQRDhgyJOorEhAqIiAAwdepU1q1bF3UMiREVEBEBwMyorKyMOobEiGYgIgJAUVERK1asiDqGxIg6EBEB4PbbbyeRSPCzn/0s6igSEyogIgLAE088wfvvvx91DIkRFRARAeCUU07hwIEDUceQGFEBEREAXnnlFUpLS8nLy4s6Sr027EhwyZMrWF++l1NyOjB9WD7dOrWPOlbG0RBdRAD41a9+xdSpU6OOUa8RM5ZzzsQlrN1ewaGqKtZur6D/5D9EHSsjqQMREQBmzZrF2rVro45Rp373F7FmW0WttyU+PkTuuN9/evm047N5+frezRUtYzVaATGzmUChuy9swLbdgG+6+/wk99ERWAzsCP981d2X17JdAXCFuw9OZn2RTNa1a1cqKmp/gk4H1QtC30lFrCuvoLIKWmVB95xsFo1VwWhuUXUg5wKnA0kVEOAbwAfufmHjRxLJbAsXLmTDhg2xmIHcPbAHFz/xJyoPVdG6VRZ3D+wRdaSMlFVVVdUoC4UdyJeBDgSF6RJ3f9/MrgGGAFVAIfAQ8B7QHrga2A3cTjCPyQaGuPuaWtb/IrAcOAF4BDg5XG89MAP4JFxjCNAduA/YDhwHzHf3O+rInQs8C2wBTgIWuPsva2zWOAcpzRQXF8fiyaKp6TgECgoKSCQSLF/+D0192vnabQtJfHzo08vtv9iav97Vv1H3ocfFp7LqvKGRC8hqd59gZgOAS4HxwGPA4d5yEXAl8B3gdHcfZ2ajgRfdfbOZ3QJkufuv69hHAeGpqcOnzAiKxVeBm4BeBEWjM0Gh+iegNbDB3TvXsWYusAIwgmL2OnCVu//58DYlJSVV+/btS+WwpLX9+/dz1FFHRR0jcjoOgfLycg4cOMBJJ50UdZRaXTG3jNI4FOC+AAAELElEQVRdBxu07clfasOjA7t+rv3pcRHIy8urs4A09imsw2+FWAZMBHoQdAqLw+s7AqfWuM8mYIqZVQAnAm8kuc/pwC+AhQQF4Jbw+lXufgDAzD6pZ42/uPvfwm3/RFBMPi0gubm5SUaKB73CCug4BPLy8tL6WBRVy9UcM5B0PhbporHfxvvt8M9ewCrACU5X9XH3AmAm8C5QWW3f04AR7j4c2MwR2qU6DASWuvv3CE5F/SK8PpnWKs/M2ptZa+Cfgb8mmUEk9ubPn8+SJUuijtEg04fl0z0nm9ZZWXTPyWb6sPyoI2Wkxu5AzjSzVwmevEe6e6mZLQZeN7O2BDOMTcBK4Jdm9mdgNrDUzPYC2whmHMn4X+BJMxtPcLrqeuCYJNf4mKD4HA885+5/SfL+IrF33333kUgkGD16dNRR6tWtU3u96yoNNNoMJK7CGUihu595hM1a5EFSix7QcQh8+OGHrFmzhrPOOivqKGlBj4tPNdsMpFGY2RyCd3RVt9vdB36ONS8neIdWTTenuqZIS9K5c2fKy8ujjiExkpYFxN0HNcGajxG8I6w2R+o+RDLCnDlz2Lhxo151S4Ppu7BEBIApU6Ywe/bsqGNIjKRlByIizW/u3Lm4e9QxJEbUgYgIAMceeyxHH3101DEkRtSBiAgAzzzzDJs2bdIMRBpMHYiIAPDII49QWFgYdQyJEXUgIgLASy+9xOrVq6OOITGiDkREAGjfvj3t2rWLOobEiDoQEQFg9uzZmoFIUtSBiAgAjz/+OM8//3zUMSRG1IGICACLFi2iuLg46hgSI+pARASANm3a0KZNm6hjSIyoAxERAGbOnMnmzZs1A5EGUwciIkBQQF588cWoY0iMZPy/B9IQb731lg6SiGSsnj171vpvgqiAiIhISnQKS0REUqICIiIiKVEBERGRlOhtvBnMzI4FZgPHAF8Exrr7m9Gmio6Z/Qj4ibsPiTpLczOzVsDDwDeBA8Cl7v5+tKmiY2b/DPzW3QuizpLO1IFktrHAYnfvDQwHHoo2TnTMbDLwGzL378QPgaPc/TvAOOC+iPNExsxuAh4Hjoo6S7rL1L8sErgfmBr+/AVgf4RZorYMuDLqEBE6G1gI4O5/BP5vtHEitQ4YFHWIONAprAxhZpcA19e4eoS7rzCzLgSnssY0f7LmdYTj8IyZFUQQKV0cA+yudvmQmX3B3T+JKlBU3P15M8uNOkccqIBkCHefDkyveb2ZfQMoBG5w96JmD9bM6joOwkdA9X8QvVUmFg9Jjk5hZTAz+xrwLDDE3RdEnUci9QYwAMDMzgRWRhtH4kAdSGb7DcGgcLKZAex294HRRpKIvAD0NbNlQBYwIuI8EgP6KhMREUmJTmGJiEhKVEBERCQlKiAiIpISFRAREUmJCoiIiKREBURERFKiAiIiIin5//WiAlnL9FMUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_bait_persons.trace,\n",
    "              varnames=['logit_psi', 'beta_psi_bait', 'beta_psi_p', 'alpha', 'beta_fish_p'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [beta_fish_p, alpha_bait, beta_psi_multi, logit_psi]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:49<00:00, 242.51draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_bait_persons_2:\n",
    "    base_rate = data.fish_caught.sum() / data.hours.sum()\n",
    "\n",
    "    logit_psi = pm.Normal('logit_psi', mu=0, sd=1.5)\n",
    "    beta_psi_multi = pm.Normal('beta_psi_multi', mu=0, sd=1)\n",
    "\n",
    "    alpha_bait = pm.Exponential('alpha_bait', lam=1/base_rate, shape=2)\n",
    "    beta_fish_p = pm.Normal('beta_fish_p', mu=0, sd=2)\n",
    "\n",
    "    psi_i = pm.math.invlogit(\n",
    "        logit_psi\n",
    "        + beta_psi_multi * (data.persons.values > 1).astype('int')\n",
    "    )\n",
    "\n",
    "    fish = (\n",
    "        data.hours.values\n",
    "        * alpha_bait[data.livebait.values]\n",
    "        * pm.math.exp(beta_fish_p * data.persons.values)\n",
    "    )\n",
    "\n",
    "    fish_obs = pm.ZeroInflatedPoisson(\n",
    "        'fish_obs',\n",
    "        psi = psi_i,\n",
    "        theta = fish,\n",
    "        observed = data.fish_caught.values\n",
    "    )\n",
    "    \n",
    "    m_bait_persons_2.name = 'livebait & persons, v2'\n",
    "    m_bait_persons_2.trace = pm.sample(5000, tune=1000, chains=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEcCAYAAADQqlM0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucVXW9//HXSKTgmBJO0gWcA+nHKetYHMw8KniB0EISSziEMoiK4SXCNCS8ZkVppqQiIDAEP0MFb6SSiEooduCYFtj44TowcpGRABk2IDrz+2MtbBzntoc9s2Z/eT8fDx+z99prf9fnu7esz/58P/uSU1lZiYiISEMdlHQAIiKSXZQ4REQkLUocIiKSFiUOERFJixKHiIikRYlDRETS8omkAxBpamZ2MDAV6Ay8C1zh7ivM7DzgDqA03vUm4FXgSaANMMzd/2FmpwD/7e6/rmX8y4BBQAXQGviZu7/YiDgvBzoA9wM3uvtwMysBjnP33VX2K4y3jWrAmEXATHefW8vtnYD/dPc56cZbz3F7AJe7+4BMjistgyoOORBcCpS7+0nAVcA98fauwHXu3iP+bwHQiyhxDAeGmlkO8CPg7poGNrMBQE/gTHfvQZRAppvZkY0N1t03ufvwxt4/TWcA/91Mx5JAqOKQA8GXgGcA3N3NrCDe3hX4mpmNABYDPwXKiaqNNsBOYCDwWNVX/NUMA0a6+954/DVmdoK7b4lf7beP//s2cB1wKtAKuNPdH4mrmbuBrcD7wF/NLJ+oSjgpPsaEeNvbwOCqBzezq+IYK+P7jKspyLgC+CnwHlHlNRMYC4wC2prZImANMA7IAbYAFwNfA34d3+85oIe7nx6P+SfgBqALcAVRtVUJnFft2FOBL8aP6d3uPr2Wx1KyhCoOORC8DnzHzHLM7CTg82bWCphHVIGcBuQClxOdHI8iSggTiU6CfzezCWZ2XQ1jfw5YXXWDu2+pcvV5dz8ZOAn4D3c/BTgd+JmZHQGMB/7H3c8iOnHXZLy7dwdKiKonAMzsS0B/4BSihPRdM7M6HoejgfPjWK5z9w+IkseD7v4kMIloGa8H8DRRogM4xN1PdfdbgEPM7Ggz+yxwpLu/BhwLfDue2z+Bb1WJ8TCix7cf0Bv4oI74JEuo4pADwRSgAFgIvAy86u4fmNkUd98GYGZPAOe7ewUwIt52PVE1MIYowdxsZse6+/IqY68FOgLb920ws28B/4ivevz3K0BXM3sxvt4ayAeOqjLey0SvzKt6z93/Gl9eRLQstiS+fjxRMpgfX28HHFPlmNUtdff3gffNbFcNtxcA98W5pzWwotocACYDFwF7iPpGAJuBaWZWDhwHvLJvZ3ffEVd0E4FPATNqiU2yiCoOORB0A+bHr4gfAVbHvYt/mNkX4n3OJGqMA2BmnwHM3RcCbYleKVcCh1Ybewpwg5l9Ir7fscAD/PuVdUX8903ghfjV/BnAw8AqYH2VpbNuNcT+STM7Ib58KrCsym0OvAGcHo9bxL8TVk1q+mK6Cv59HnDgonis64A/VZsDREtc3yGqxB40s8OBW4ABwCXALqKlLgDiyqSru59HtFz3m32PlWQvJQ45EKwARpjZK8DPiXoSlUQnukfNbAFRcphU5T5jgNviy/cBfwY+C/y96sDuPhP4K/CSmf2F6FX4IHffXC2GOUC5mS0kSlCV7r6DaEnsD2Y2n6h6qG4PcFU89meIktK+Y/+dqNp4ycz+j6jaWN/whwWApUDfuMn/wziWl4iWsD6WhNy9PH4MiuP43yWqlF4hquh2ES3f7bMJ6BD3UOYBd8RVj2SxHH07roiIpEMVh4iIpEWJQ0RE0qLEISIiaVHiEBGRtChxiIhIWvR+6kZ69dVX9XY0EQle165dc6pv09txGy8rHrhJkyaxceNGbrzxxqRDyZji4mIKCgrq3zGLhDYnzadlS2M+H0saoIojeJdeeinFxcVJhyEiAVGPQ0RE0qKKI3A9evQglUqxePHipEMRkUAocQSusLCQDRs2JB2GiAREiSNwhYWF6nGISEapxxG4vXv3snfv3qTDEJGAqOIIXM+ePdXjEJGMUuII3CWXXML69en+RIOISO2UOAI3aNAg9ThEJKPU4whcKpVi166afl5aRKRxVHEE7pxzzlGPQ0QySokjcD/84Q/V4xCRjFLiCFz//v3V4xCRjFKPI3Dbt29nx44dSYchIgFRxRG4vn37qschIhmlxBG4q6++mrfeeivpMDJm3ZYUwx4vZf0f1tA571AmD+5Gp/Ztkw5L5ICixBG4fv36NVuPY8jUxbzgZc1yLIAVm8s57fYXmuVYp1seU4ec2CzHEmnplDgS0Ot3C1j+dnkzH3V1Mx8vLC94GfmjnmriozTfc3TsUbk8++PuzXY8CYsSRwKa8x9saL/H0fPOBazcXE4lcFAOdMnLZd7I7D8BhvbTpBK2FvuuKjMrMrPead6ng5ndF18+zcy+mqFY7jKzTpkYq7ldc801FBYWJh1Gxkwe3I2Oh7emVU4OXfJymTy4W9IhiRxwgqo43H0TMDy+ejEwE/hHBsYdsb9jJKVPnz7BfI5j3ZYUQ6ct4a13936YNNQYF2l+LT5xmFlrYCrQGWgF3OnuD5nZicC9wA5gM7AbuJkoWVwB9Aa+bmb/dPd1NYxbCHwXOAw4ErjV3Web2S+A04kem9nu/mszexG43N3fbMq5NoVNmzZRVlbW5MsgaoyLHDhafOIAhgFl7j7IzA4D/mZm84H7gQvd/Y34ZP/5fXdw91fNbC4ws6akUcWhQE8gD1hsZk8APwB6ABuBwtruWFJS0iRfHnj5E6Ws3dYEP7z0dFM3dsPVPI1xSOoNDEcf0Zr7+3bM6Ji7d+8OptKFA3c+tb3gzIbEUQA8B+DuO8zsn0AX4HPu/ka8z0JgQCPGXuDuFcDbZraVKIH8ABgLdACeqe2O+fn5jThcAwLKcGUwd+5c1q1bx2WXXZbRcZPQ884FrCorp6IyrMY4hNcc13xatv2dT4ttjldRDJwKEFccXwHWAKVm9qV4n5NquF8F9c+vazzuUcCngH8B3wf+h2i5qtDMjt7fCSSpd+/enHrqqUmHkRGTB3ejS17uh0lDjXGRZGRD4pgItDezl4AXgVvcfTNRE3yKmT0HnAhUX9/5X2CsmdWVVjvEy15PAcPdfQ9R8vgr8ALwLFDXUleLV1paysaNG5MOIyM6tW/LvJHdeeqizswb2V2NcZGEtNilKncvrHJ1cA27nAj0cfcyM7sNeM/dS4irD3efAEyo5zAL3H1UtePeCtxabb8eDY+8ZbnwwguD+hyHiCSvxSaOBngbeNbMyoHt1JxciD/X8aUabnqoCWNrMcaMGcPatWuTDkNEApK1icPdZwGzGrDf8Pr2CdlZZ50V1LtBRCR52dDjkP2wevVqSktLkw5DRAKStRWHNMzFF19MKpWiV69eSYciIoFQ4gjcLbfcoh6HiGSUEkfgunfvrh6HiGSUehyBc3fWrFmTdBgiEhBVHIEbNmwYqVSKc845J+lQRCQQShyB++Uvf0lJSUnSYYhIQJQ4AnfyySfTrl27pMMQkYCoxxG4ZcuWsWLFiqTDEJGAqOII3JVXXkkqleLcc89NOhQRCYQSR+Buv/12vatKRDJKiSNw3bp1Izc3N+kwRCQg6nEE7vXXX9cHAEUko1RxBG7EiBGkUin69euXdCgiEggljsDdddddrF69OukwRCQgShyBO+GEEzj44IOTDkNEAqIeR+CWLFnC0qVLkw5DRAKiiiNw1157LalUigsuuCDpUEQkEEocgbvnnnvU4xCRjFLiCNzxxx9Pq1atkg5DRAKiHkfgFi1axGuvvZZ0GCISEFUcgRs9ejSpVIqBAwcmHYqIBEKJI3ATJkxg1apVSYchIgFR4gicmVFRUZF0GCISEPU4ArdgwQKWLFmSdBgiEhBVHIG76aabSKVSXHTRRUmHIiKBUOII3JQpU1i5cmXSYYhIQJQ4Ate5c2f27NmTdBgiEhD1OAL33HPPsWjRoqTDEJGAKHEE7rbbbmPChAlJhyEiAdFSVeCmT5/OihUrkg5DRAKixBG4jh07Ul5ennQYIhIQLVUFbu7cuSxcuDDpMEQkIEocgRs7diwPPPBA0mGISEC0VBW4mTNnsnz58qTDEJGAKHEErkOHDmzdujXpMCQQ67akGDptCavLdtI571AmD+5Gp/Ztkw5LmpkSR+DmzJlDaWkpBQUFSYciGTBk6mJe8LKkwwBgxeZyTrv9hTr2aNpfnjzd8pg65MQmPYbULKeysjLpGLJVVjxwPXr0IJVKsXjx4qRDyZjuY59l7ba9SYchkqhjj8rl2R93b9R9i4uLG/piMqemjao4Ajdr1qzgehz39+0YXAWVxj/kRPW8cwGrysqpqISDcqBLXi7zRn785JUt82mo0Oazvxr8riozKzKz3g3ct5OZ9Wl8WA06xkwz+2QTjLsp/vsVMzutKY/VHI488kjatWuXdBgSiMmDu9ElL5dWOTl0yctl8uBuSYckCWiqiuMM4DhgThONj7sPaKqxY+cDm4C/NMOxmsyjjz7KW2+9pVdLkhGd2rdl3sjuHzbJT7/jRTXJD0DpJo7hZnZtfL+h7r7SzK4CBhKt+c8E7gVGAW3NbBGwHbiJqLrJBQa6e41rJ2Z2M1HC+QzQDrjK3V8ys6nAF4E2wN3uPt3MSoDj3H13LWOtBBYBxwLzgcOBEwF39wvNrAiY6e5z40pqgLsXxvf9PFAIvGdmfwMerutYLdm4ceNIpVJcffXVSYciGZQdTfKmbY5XpUZ580o3cSxy97Fmdg7wGzMbA/QHTolvnwf8GRhLdKJ90syGA4PcfYOZjQa+D/yijmOk3P0MM/sy8KCZnQKcBpxElJx6NTDWfKLKZyPwL+AbwFXAajM7oq47uvv6OLFscvfFZvaxfUpKSti1a1cDQ0nOb37zG3bv3k1xcXHSoWRMffO5/InSLG2eN9+JNjQveBn5o55q4qNk5/Nz9BGtub9vx49sa+g5obaVinQTx1/iv4uA24HjgaOJXtFDVCUcU+0+64FxZlYOfB54uZ5jPA/g7m+YWQd332FmI4CJwKeAGQ2MdYu7rwMws53u/s/48nbgkGr71vjOgbrk5+ene5fEhNbYq28+C7Jwrtn2HNXXJM+2+dRH8/modL9yZF8teCqwDHDgDeB0d+8BFAH/ACqqjD0JGBIvA22g/pN0VwAzOx5Yb2afBbq6+3nAt4kqnYYkvPreLrsb+Gx8+es13F51DlnroYce4plnnkk6DAmMmuQHtnQrjpPM7Hmik/LF7r7WzOYDL5nZwcBiogpjKfCzuD8wA1hoZjuBt4HP1XOMr8VjHgpcStSg7hD3Sz4A7nD392taPkrTA8AUM/sBUFPP5VXgdjPL6jWe8ePHk0qlGDlyZNKhSED2NcnlwNSiPgAYN8c3ufv9ScfSAC3ngatDKpXizTff5Otfr6moyk6hLRtAeHPSfFq2rPwAoJk9Cny62ubtwGtpjnMuUNNL6bvd/bFGhheUtm3b0qZNm6TDEJGAJJI43L1fhsZ5EngyE2OFasaMGaxfvz6oV0sikqysb/5K3R544AFmz56ddBgiEhB9V1Xg5s2bF9RnOEQkeao4Ate6dWtat26ddBgiEhBVHIErKipiw4YN6nGISMao4ghcUVERjz/+eNJhiEhAVHEE7sUXX1SPQ0QyShWHiIikRYkjcJMmTeKRRx5JOgwRCYgSR+D0JYcikmnqcQTuueeeU49DRDJKFYeIiKRFiSNw9913H3/84x+TDkNEAqKlqsDNmTOH8vLypMMQkYAocQTumWeeUY9DRDJKS1UiIpIWJY7A3X333UyfPj3pMEQkIFqqCtz8+fPV4xCRjFLiCNyTTz6pHoeIZJSWqkREJC1KHIG74447mDJlStJhiEhAtFQVuFdeeYUdO3YkHYaIBESJI3CzZ89Wj0NEMkpLVSIikhYljsCNHTuWSZMmJR2GiARES1WBe/3113n33XeTDkNEAqLEEbiZM2eqxyEiGaWlKhERSYsSR+B+/vOfM378+KTDEJGAaKkqcO7O9u3bkw5DRAKixBG4GTNmqMchIhmlpSoREUmLEkfgbrzxRn7/+98nHYaIBERLVYErLS1l27ZtSYchIgFR4gjc1KlT1eMQkYzSUpWIiKRFFUfgrr/+et555x1uGHs3Q6ctYXXZTjrnHcrkwd3o1L5t0uGJSBZSxRG4LVu2sH37doZOW8KKzeV8UFnJis3lDJ22JOnQRCRLqeII3MSJEykuLuY7f1jzke0rNpcnFJGIZDtVHAeIznmHclBOdPmgHDjmM7nJBiQiWatJKg4z2+TuHWq5rQdwubsP2I/xbwY2ufv9Ddi3N9DJ3Sea2WXAVHffW22sgcDNwGKgCKgElgFXALcClwCF7j63sTEn5Sc/+Qlbtmxh8h33fqzHISLSGMFXHO4+190nxldHA61q2O1Od38QuBMY4+6nAjlAX3cfA2Rdwthn165d7NmzJ+kwRCQg+1VxmNnxRCfbVsCRwA/dfVGV218E3gSOIzoR949vOsbMngE+A8xx95vNrDtwE1EyywUGuvvyOg5/npldALQFrnb3xWZ2JdAPOBR4BziPqJo4DlgBdABmAt+tZcyuwIL48jNAL+Cxhj0aLdO9995LcXHxh81x4MPm+LyR3ROOTkSy0f4uVX0ZuMbdl5rZQGAIsKjaPovc/XIzG070iv9R4BCik3crYB3RMtGXgUHuvsHMRgPfB35Rx7HXxON+GZhuZv8FtAfOcvcKM/sz8OF6jLtPNrMbgLqWyHLcvTK+vAM4vLYdS0pK2LVrVx1DtRy7d+9mVdlHm+ErNpdn7QcDd+/enbWx1ya0OWk+LVtD51NQUFDj9v1NHOuBG8xsF3AYUNNvlD4f/10E9I0vL3P3PQBm9n6VscaZWTnweeDleo79FwB3f8PMOsTJ4j3gj/EYXwBapzmfiiqXDwNq/a6O/Pz8NIdOxogRI/jXv/5FlxOGsqqsnIrKqDneJS+31v8pWrri4uKsjb02oc1J82nZ9nc++9vjGAfc5O6DgaVEy1HVdY3//jfwRny5sob9JgFD3L0Q2FDLWFWdCGBmXwHWmdlXge+6e3/gKqK5VR+jgrrn/FrcvAc4G1hYTwxZY/LgbnTJy6VVTg5d8nLVHBeRRtvfimMG8IiZbQXeIupzVFdoZiOBncCFwFfqGGuhme0E3gY+V8+x/8PMngcOBoYBK4GdZravUtlYwxgLgafN7PQqS1JVXQNMMrNPAsXArHpiaPHuuusuiouL6dS+rXoaIpIR+5U43P1OouZ49e1V34p7vbu/WeX6i/F/H9nX3Uemcdyba7npjHruN7ie25cDOruKiNShRX9y3MweBT5dbfN2d+9b0/77YaSZvRu/Jbd6DLcBvYnejZV1rrjiCrZu3cqDD35saiIijdKkicPde+zn/ftlKJS6jnEz0bu6art9DDCmqeNoKm3atCGVSiUdhogEpEVXHLL/7rjjjqDeRigiyQv+k+MiIpJZqjgCd9lll7Ft2zYefvjhpEMRkUAocQSuffv2VFbW9M5jEZHGUeII3K9+9Sv1OEQko9TjEBGRtKjiCNyQIUPYtm0bjz2W1V/yKyItiBJH4Dp27EibNm2SDkNEAqLEEbhbb71VPQ4RySj1OEREJC2qOAI3aNAgtm/fzpw5c5IORUQCocQRODOjrKws6TBEJCBKHIG74YYb1OMQkYxSj0NERNKiiiNwAwYM4N133+Xpp59OOhQRCYQSR+BOOOEENm/enHQYIhIQJY7AjRo1Sj0OEcko9ThERCQtqjgCd/7557Njxw6effbZpEMRkUAocQTum9/8Jm+//XbSYYhIQJQ4AveTn/xEPQ4RySj1OEREJC2qOAJ37rnnUl5ezvPPP590KCISCCWOwJ155pnqcYhIRilxBO5HP/qRehwiklHqcYiISFpUcQTu7LPPpry8nIULFyYdiogEQokjcH369GHTpk1JhyEiAVHiCNzw4cPV4xCRjFKPQ0RE0qKKI3BnnXUWO3fu5JVXXkk6FBEJhBJH4Pr378/GjRuTDkNEAqLEEbhLL71UPQ4RySj1OEREJC2qOALXo0cPUqkUixcvTjoUEQmEEkfgCgsL2bBhQ4P2XbclxdBpS1hdtpPOeYcyeXA3OrVv28QRiki20VJV4AoLCznvvPMatG/vu//Cis3lfFBZyYrN5QydtqSJoxORbKTEEbi9e/eyd+/eevfr9bsFpN774CPbVmwuJ3/UU/T63YKmCk9EspCWqgLXs2fPBvU4nv1xd3reuYBVZeVUVMJBOdAlL5d5I7s3U6Qiki0O+IrDzArNbJ2ZjYyvjzCzsfHl75nZm/uuZ6NLLrmE888/v0H7Th7cjS55ubTKyaFLXi6TB3dr4uhEJBup4og8CIw3s/8HnAjMBnD3WWaWCxyXZHD7Y9CgQfoch4hk1AFfcVRxCDAN+EXSgWRSKpVi165dDdpXzXERaQhVHDF33wo8a2aFDdm/pKSkwSfkJA0ePJiKigqmT59e536XP1Faa3P86CNac3/fjk0ZZlp2794dXBUV2pw0n5atofMpKCiocbsSRyPl5+cnHUKDXHPNNaxfv77W/wH2WVBQkDXN8eLi4nrnk21Cm5Pm07Lt73y0VBW4/v37c/bZZzdoXzXHRaQhVHEEbvv27ezYsaNB+3Zq37ZFVhgi0rIocVTj7kVJx5BJffv21XdViUhGaakqMnDf5ziqMrPvAaMSiCdjrr76agYNGpR0GCISkAO+4ogrjKJabpsFzGrOeDKtX79+Qb0bRESSp4ojcO+88w5bt25NOgwRCcgBX3GE7nvf+556HCKSUUocgbvmmmsoLS1NOgwRCYgSR+D69OmjHoeIZJR6HIHbtGkTZWVlSYchIgFRxRG4AQMGqMchIhmlxBG4UaNGsW7duqTDEJGAKHEErnfv3upxiEhGqccRuNLSUjZu3Jh0GCISEFUcgbvwwgvV4xCRjFLiCNyYMWNYu3Zt0mGISECUOAJ31llnqcchIhmlHkfgVq9erU+Oi0hGqeII3MUXX0wqlaJXr15JhyIigVDiCNwtt9yiHoeIZJQSR+C6d++uHoeIZJR6HIFzd9asWZN0GCISEFUcgRs2bBipVIpzzjkn6VBEJBBKHIH75S9/SUlJSdJhiEhAlDgCd/LJJ9OuXbukwxCRgKjHEbhly5axYsWKpMMQkYCo4gjclVdeSSqV4txzz006FBEJhBJH4G6//Xa9q0pEMkqJI3DdunUjNzc36TBEJCDqcQTu9ddf1wcARSSjVHEEbsSIEaRSKfr165d0KCISCCWOwN11112sXr066TBEJCBKHIE74YQTOPjgg5MOQ0QCoh5H4JYsWcLSpUuTDkNEAqKKI3DXXnstqVSKCy64IOlQRCQQShyBu+eee9TjEJGMUuII3PHHH0+rVq2SDkNEAqIeR+AWLVrEa6+9lnQYIhIQVRyBGz16NKlUioEDByYdiogEQokjcBMmTGDVqlVJhyEiAVHiCJyZUVFRkXQYLcq6LSmGTlvC6rKddM47lMmDu9GpfdukwxLJGkocgVuwYAFr166loKAg6VBahCFTF/OCl314fcXmck67/QVOtzymDjkxwchEsocSR+BuuukmUqkUF110UdKhNLtev1vA8rfLG7TvC15G/qinPrLt2KNyefbH3ZsiNJGspsQRuClTprBy5cqkw0hEbSf9nncuYFVZORWVcFAOdMnLZd5IJQiRhmrSt+OaWZGZ9W7gvp3MrE8jjtHOzP5mZvPMbJSZ1bjeYGY9zGxmuuNnu86dO9OxY8ekw2hRJg/uRpe8XFrl5NAlL5fJg7slHZJIVmlJFccZwHHAnDTv9xVgjbufn/mQst9zzz2nHoeIZFRzJI7hZnZtfKyh7r7SzK4CBgKVwEzgXmAU0NbMFgHbgZuIKqJcYKC7L68+sJl9EhgHfM7MbgGOjsdbDUwF3o/H2PchhmPM7BngM8Acd7+5poDNLB94BNgIfAF4xt1/tp+PQyJuu+02UqkUQ4cOTTqUFkHNcZH9l1NZWdlkg5tZEfCmu481s3OAS4AxwERg36LyPOCHwDeB49x9lJkNBx539w1mNhrIcfdf1HKMHsDl7j4gPt5MoAvwReA64FRgM3AkUYL6OtAKWOfuR9YyZj6wBDCiJPYScIW7/23fPiUlJZW7du1qzMPSrDZu3MiePXvIz89POpSM2b17N4cccki9+13+RClrt+1t9HGOPqI19/dtnmW+hs4pW2g+LVtD51NQUJBT0/bmqDj+Ev9dBNwOHE9UGcyPt7cDjql2n/XAODMrBz4PvJzmMScDPwXmEp34R8fbl7n7HgAze7+eMf7u7v+K9/1foiTyYeLIlhNxQUEBxcXFQS1VNXQ+C2rZpyU2xw/U5yhbaD4f1RzfVbWv/j8VWAY48AZwurv3AIqAfwAVVeKZBAxx90JgA1Bj1qtDX2Chu59JtOT003h7OuVVgZm1NbNWwDeAf6YZQ4swd+5cFi5cmHQYLYqa4yL7pzkqjpPM7Hmik/bF7r7WzOYDL5nZwcBiogpjKfAzM/sbMANYaGY7gbeBz6V5zP8DppnZGKJlqR8Dn0pzjPeIks5RwCx3/3ua928Rxo4dSyqV4rLLLks6lBajU/u2iVcYItmsSXsc2Sruccx095Pq2C0rHrhNmzaxfPlyTjvttKRDyZjQlg0gvDlpPi1bGvNJrMeREWb2KPDpapu3u3vf/RjzMv79jquqrm/smC1Nhw4d2Lp1a9JhiEhAsiZxuHu/JhhzItE7vGpSV7WRNebMmUNpaWlQr5ZEJFn6IafA/fa3v6WoqCjpMEQkIFlTcUjjzJo1i+XLP/bZSRGRRlPFEbgjjzySdu3aJR2GiAREiSNwjz76KPPmzUs6DBEJiBJH4MaNG8eMGTOSDkNEAqIeR+CeeOIJ3D3pMEQkIKo4Anf44Ydz2GGHJR2GiAREFUfgHnroIdavX6/PcYhIxqjiCNz48eOZOfOA++FDEWlCqjgC9/TTT/Pmm28mHYaIBEQVR+Datm1LmzZtkg5DRAKiiiNwM2bMUI9DRDJKFUfgHnjgAWb6kiLaAAADe0lEQVTPnp10GCISEFUcgZs3bx7FxcVJhyEiAVHFEbjWrVvTunXrpMMQkYCo4ghcUVERGzZsUI9DRDJGFUfgioqKePzxx5MOQ0QCot8cb6RXX31VD5yIBK9r164f+91xJQ4REUmLlqpERCQtShwiIpIWJQ4REUmL3o4bEDM7CLgP+E9gD3CJu6+scvulwDDgfeA2d/9TIoE2UAPmczdwCrAj3tTX3bc3e6BpMrNvAL929x7VtvcBbiR6fqa4+6QEwktbHfP5MXAJUBZvGuYt+FfFzKw1MAXIBw4m+jfyZJXbs+75acCcGvUcKXGE5bvAIe7+TTM7Cfgt0BfAzDoAVwP/BRwCvGRm89x9T2LR1q/W+cS6At9y93cSia4RzOw64EJgZ7XtrYHfAd3i2142syfd/e3mj7LhaptPrCtwkbu/2rxRNdogYIu7X2hmnwZeB56E7H1+qGNOsUY9R1qqCsspwFwAd/8rUZLY50TgZXffE78qXwl8tflDTEut84mrkWOAiWb2spldnEyIaVsF9KthewGw0t23uvt7wEvAac0aWePUNh+ITkrXm9lLZnZ9M8bUWI8AN8SXc4gqi32y9fmpa07QyOdIiSMsnwKqLtV8YGafqOW2HcDhzRVYI9U1n0OB3xO9ouoNDDezlp4IcffZwN4absrG56eu+QDMBC4HzgBOMbPvNFtgjeDu5e6+w8wOA2YBY6rcnK3PT11zgkY+R0ocYXkXqPoD4we5+/u13HYYsK25AmukuuaTAu5295S77wCeJ+qFZKtsfH5qZWY5wF3u/k78Cv0p4GsJh1UvM+sIvABMd/cHq9yUtc9PbXPan+dIPY6wvAz0AR6OewJLq9y2GPiFmR1C1CQrAJY1f4hpqWs+xwIPmdnXiF4AnQJMa/4QM6YYOCZehy4nWga5I9mQ9sungGVmVkDUEziDqEnbYpnZUcCzwJXuPr/azVn5/NQzp0Y/R0ocYXkM6Glmi4jWM4eY2UiitdknzWwcsJDoRPszd9+dYKwNUd98pgN/JVoq+YO7v5FgrI1iZgOBXHefGM/tz0TPzxR3X59sdOmrNp/RRK909wDz3f3pZKOr12igHXCDme3rC0wCDs3i56e+OTXqOdJXjoiISFrU4xARkbQocYiISFqUOEREJC1KHCIikhYlDhERSYsSh4iIpEWJQ0RE0vL/AcsOqccsGuxpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_bait_persons_2.trace,\n",
    "              varnames=['logit_psi', 'beta_psi_multi', 'alpha_bait', 'beta_fish_p'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [beta_fish_ch, beta_fish_p, alpha_bait, beta_psi_ch, beta_psi_multi, logit_psi]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:52<00:00, 229.69draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_all:\n",
    "    base_rate = data.fish_caught.sum() / data.hours.sum()\n",
    "\n",
    "    logit_psi = pm.Normal('logit_psi', mu=0, sd=1.5)\n",
    "    beta_psi_multi = pm.Normal('beta_psi_multi', mu=0, sd=1)\n",
    "    beta_psi_ch = pm.Normal('beta_psi_ch', mu=0, sd=2)\n",
    "\n",
    "    alpha_bait = pm.Exponential('alpha_bait', lam=1/base_rate, shape=2)\n",
    "    beta_fish_p = pm.Normal('beta_fish_p', mu=0, sd=2)\n",
    "    beta_fish_ch = pm.Normal('beta_fish_ch', mu=0, sd=2)\n",
    "\n",
    "    psi_i = pm.math.invlogit(\n",
    "        logit_psi\n",
    "        + beta_psi_multi * (data.persons.values > 1).astype('int')\n",
    "        + beta_psi_ch * (data.child.values > 0).astype('int')\n",
    "    )\n",
    "\n",
    "    fish = (\n",
    "        data.hours.values\n",
    "        * alpha_bait[data.livebait.values]\n",
    "        * pm.math.exp(\n",
    "            beta_fish_p * data.persons.values\n",
    "            + beta_fish_ch * data.child.values\n",
    "        )\n",
    "    )\n",
    "\n",
    "    fish_obs = pm.ZeroInflatedPoisson(\n",
    "        'fish_obs',\n",
    "        psi = psi_i,\n",
    "        theta = fish,\n",
    "        observed = data.fish_caught.values\n",
    "    )\n",
    "\n",
    "    m_all.name = 'all'\n",
    "    m_all.trace = pm.sample(5000, tune=1000, chains=2)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEcCAYAAADQqlM0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt4VPW1//F3oBTBWIsYpVowheoyLZ5qOXio12DBQ+lRqrWVw0G5eEdtEU8tWrzbU3q0tFLrDbn94PiLrZcqp0pFxIjG/qC2VrFxgUCAco0UkDBBkOT3x97UEBLIJHtnz0w+r+fhSWbP3t+98n2GWbO+a89MXm1tLSIiIk3VLukAREQkuyhxiIhIWpQ4REQkLUocIiKSFiUOERFJixKHiIik5VNJByASNzPrCEwHegIfAte6+zIzuwC4D1gT7no78CbwHNAJuMrd3zazM4DT3f2njYx/JTAcqAE6AD9y91eaEefVQDfgYeA2dx9jZhXAie6+s85+I8Nt45sw5gygxN3nNnJ/D+Ar7j4n3XgPct5i4Gp3HxrluJIZVHFIW3AFUOXu/YDrgQfC7X2Am9y9OPxXCpxLkDjGAJeZWR7wfeD+hgY2s6HAQODr7l5MkEBmmdmRzQ3W3Te4+5jmHp+mc4DTW+lckiNUcUhb8CXgBQB3dzMrCrf3AU4xs7HAIuCHQBVBtdEJ2AEMA56p+4q/nquAce6+Oxx/pZmd7O6bw1f7XcN/3wRuAs4E2gOT3P03YTVzP7AF+Bj4g5kVElQJ/cJzPBJu2wiMqHtyM7s+jLE2PGZyQ0GGFcAPgV0ElVcJMBEYD3Q2szJgJTAZyAM2A6OBU4Cfhse9BBS7e/9wzP8FbgV6AdcSVFu1wAX1zj0d+GI4p/e7+6xG5lKyhCoOaQveAv7NzPLMrB9wrJm1B+YRVCBnAfnA1QRPjkcTJIRHCZ4E/2Jmj5jZTQ2MfQywou4Gd99c5+bL7n4a0A/4grufAfQHfmRmnwUeAv7d3QcQPHE35CF3PxuoIKieADCzLwEXA2cQJKRvmZkdYB6OA74dxnKTu+8hSB6Pu/tzwBSCZbxi4HmCRAdwiLuf6e53AoeY2XFm9jngSHf/M3AC8M3wb/sr8K91YjyMYH4vBAYBew4Qn2QJVRzSFkwDioCFwOvAm+6+x8ymuftWADN7Fvi2u9cAY8NtNxNUAxMIEswdZnaCuy+tM/YqoDuwbe8GM/tX4O3wpoc/TwL6mNkr4e0OQCFwdJ3xXid4ZV7XLnf/Q/h7GcGy2OLwdm+CZDA/vN0FOL7OOet7x90/Bj42s+oG7i8CHgxzTwdgWb2/AWAqcCnwEUHfCGATMNPMqoATgTf27uzu28OK7lHgM8DsRmKTLKKKQ9qCvsD88BXxb4AVYe/ibTP7fLjP1wka4wCY2VGAuftCoDPBK+Va4NB6Y08DbjWzT4XHnQA8xievrGvCn+8BC8JX8+cAvwaWA2vrLJ31bSD2T5vZyeHvZwJL6tznwLtA/3DcGXySsBrS0AfT1fDJ84ADl4Zj3QT8b72/AYIlrn8jqMQeN7PDgTuBocDlQDXBUhcAYWXSx90vIFiu+++9cyXZS4lD2oJlwFgzewO4m6AnUUvwRPe0mZUSJIcpdY6ZANwT/v4g8Hvgc8Bf6g7s7iXAH4DXzOxVglfhw919U70Y5gBVZraQIEHVuvt2giWx/2Nm8wmqh/o+Aq4Pxz6KICntPfdfCKqN18zsjwTVxtqmTwsA7wBDwib/NWEsrxEsYe2XhNy9KpyD8jD+DwkqpTcIKrpqguW7vTYA3cIeyjzgvrDqkSyWp0/HFRGRdKjiEBGRtChxiIhIWpQ4REQkLUocIiKSFiUOERFJi66nbqY333xTl6OJSM7r06dPXv1tuhy3+WrLy8spKio6+J5tiOZkX1OmTGH9+vXcdtttSYeSUfQ42V+Gzsl+SQNUcYjE6oorrqC8vDzpMEQipR6HiIikRRWHSIyKi4tJpVIsWrQo6VBEIqPEIRKjkSNHsm7duqTDEImUEodIjEaOHKkeh+Qc9ThEYrR79252796ddBgikVLFIRKjgQMHqschOUeJQyRGl19+OWvXpvsVGSKZTYlDJEbDhw9Xj0NyjhKHSIxSqRTV1Q19vXfrW705xWUzF7Oicgc9Cw5l6oi+9OjaOemwJAspcYjEaPDgwZH1OEZNX8QCr4wgKli2qYqz7l3Q4nH6WwHTR50aQUSSTfRZVc2nz6pqQC7Oybk/L2Xpxqqkw5ADOOHofF684eykw2iRDP2/k12fVWVmM4ASd5+bxjHdgNvcfYyZnQVsdfe3I4jlF8Akd1/d0rEk+7T0CSlTnhAGTipleWUVNbXQLg96FeQzb1wyT7aZMifSPDn1Pg533+DuY8Kbo4FjIhp3rJKGNMe2bdvYvn170mEABD2NI4KeRk0t7N5Tw+rNqYSjkmyUsRXHXmbWAZgO9ATaE7zyf8LMTgV+BWwHNgE7gTuAEuBaYBDwVTP7a0NP+mY2EvgWcBhwJHCXuz9lZj8G+hPMzVPu/lMzewW42t3fi/NvldwzZMiQSHocUfY39qrYnGpxn0M9jrYp4xMHcBVQ6e7Dzeww4E9mNh94GLjE3d8Nn+yP3XuAu79pZnMJlroOVCkcCgwECoBFZvYs8B9AMbAeGNnYgRUVFezcuVOXWtaTi3Ny9bNrWLW1me/+7vcDAArH/y7CiDLHAq9swd+2ItJYjvtsBx4e0j3SMVtTJv7faWw5MRsSRxHwEoC7bzezvwK9gGPc/d1wn4XA0GaMXeruNcBGM9tCkED+A5gIdANeaOzAwsJCqqurtU5bTy6uXZe28O/JpDnJlD5HJs1JpsimOcmGHkc5cCZAWHGcBKwE1pjZl8J9+jVwXA0H//v6hOMeDXwG+DvwHeDfCZarRprZcS39A6Tt+uCDD9iyZUvSYfzD1BF96VWQT/u8PHoV5DN1RN+kQ5IslA0Vx6PAFDN7DegE3Onum8xsDDDNzKqAXUD9z3X4f8BEM1vp7o3Vf93CZa/DgTHu/pGZ/R34A1ANvAioKS7NdtFFF2XUZ1X16No5sSupJHdkbOJw95F1bo5oYJdTgfPcvdLM7gF2uXsFYfXh7o8AjxzkNKXuPr7eee8C7qq3X3HTIxf5xI033siaNWuSDkMkUhmbOJpgI/BiWHFso+Hkgpk9CHypgbueiDE2EQDOO++8jGt4irRU1iYOd38SeLIJ+4052D4icdmwYQOVlZVZ0/QUaYqsTRwi2WDo0KEZ1eMQiYISh0iMxo8fz+rVur5CcosSh0iMBg0apB6H5JxseB+HSNZas2YN69evTzoMkUip4hCJ0SWXXKIeh+QcJQ6RGE2YMIFVq1YlHYZIpJQ4RGI0YMAA9Tgk56jHIRKjFStW6J3jknNUcYjEaPTo0aRSKc4999ykQxGJjBKHSIzuvPNO9Tgk5yhxiMTo7LPPVo9Dco56HCIxcndWrlyZdBgikVLFIRKjq666ilQqxeDBg5MORSQyShwiMfqv//ovKioqkg5DJFJKHCIxOu200+jSpUvSYYhESj0OkRgtWbKEZcuWJR2GSKRUcYjE6LrrriOVSnH++ecnHYpIZJQ4RGJ077336qoqyTlKHCIx6tu3L/n5+UmHIRIp9ThEYvTWW2/pDYCSc1RxiMRo7NixpFIpLrzwwqRDEYmMEodIjH7xi1+wYsWKpMMQiZQSh0iMTj75ZDp27Jh0GDlp9eYUl81czIrKHfQsOJSpI/rSo2vnpMNqE5Q4RGK0ePFiVq5cSVFRUdKhJGbU9EUs8MoG7omuElu2qYqz7l0Q2Xh19bcCpo86NZaxs1VebW1t0jFkq9ry8vI2/YTQEM3JvoqLi1vlO8fP/XkpSzdWxXoOyR4nHJ3PizecHcVQeQ1tbHLFYWYzgBJ3n9uEfXsAX3H3OU0OL01mVgJc6u67Ih53g7t3M7OTgC7u/mpc55Lc98ADD7RKjyOiJ4lWE8ULjIGTSlleWUVNLbTLg14F+cwbl13zUFc2veiK63Lcc4DTYxobAHcfGvMT+beBL7XSuSRH9e7dm+OPPz7pMHLS1BF96XFE0NOoqYXde2pYvTmVcFRtQ7o9jjFm9oPwuMvc/X0zux4YBtQCJcCvgPFAZzMrA7YBtxMkqXxgmLsvbWhwM7sDOBE4CugCXO/ur5nZdOCLQCfgfnefZWYVwInuvrORsd4HyoATgPnA4cCpgLv7JXUrKDMbBAx195HhsccCI4FdZvYn4NcHOpdIY8rKyqioqMiaV5JxaI0eB0DF5lQsfQ71OPaXbuIoc/eJZjYY+G8zmwBcDJwR3j8P+D0wkeCJ9jkzGwMMd/d1ZnYL8B3gxwc4R8rdzzGzLwOPm9kZwFlAP4Lk1NQvby4kqHzWA38H/gW4HlhhZp890IHuvjZMLBvcfZGZ7bdPRUUFO3fu1Ju76tGc7Gvs2LHU1NRwyimnxHqeq59dw6qtu2M9R1u1wCspHP+7Vjpb9Muax322Aw8P6d6sYxt7wZNu4ng1/FkG3Av0Bo4jeEUPQZVQvy5fC0w2syrgWOD1g5zjZQB3f9fMurn7djMbCzwKfAaY3cRYN7v7agAz2+Hufw1/3wYcUm/fBhtAB1JYWEh1dXWbfiXZkGxap20Ns2bNYvny5bHPSWmWzXlUj5Nc6nNk0/+ddHsce+u1M4ElgAPvAv3dvRiYAbwN1NQZewowKlwGWsfBn6T7AJhZb2CtmX0O6OPuFwDfJKh0mpLwDna52E7gc+HvX23g/rp/g0izmBlf+MIXkg4jZ00d0ZdeBfm0z8ujV0E+U0f0TTqkNiHdiqOfmb1M8KQ82t1Xmdl84DUz6wgsIqgw3gF+FPYHZgMLzWwHsBE45iDnOCUc81DgCmAD0C3sl+wB7nP3jxtaPkrTY8A0M/sPoKGey5vAvWamdRdpttLSUlatWpU1rySzTY+unbO2wshmGfU+jrA5vsHdH046libQ+zgaoDnZV2u9jyPb6HGyvwydk5a9jyNKZvY0cES9zduAP6c5zvnAuAbuut/dn2lmeCKRmTZtGu+//37SYYhEKpHE4e6RfFSouz8HPBfFWCJx6NmzJx999FHSYYhESs1fkRi99NJLlJWVJR2GSKSUOERidM899/DII48kHYZIpPTpuCIxmjVrFsuWLUs6DJFIKXGIxKh79+5UVelTayW3aKlKJEZz585l4cKFSYchEiklDpEYTZw4kcceeyzpMEQipaUqkRiVlJSwdGmDHwYtkrWUOERi1K1bN7Zs2ZJ0GCKR0lKVSIzmzJnDggXxfBe2SFKUOERi9LOf/YwZM2YkHYZIpLRUJRKjJ598Uj0OyTmqOERidOSRR9KlS5ekwxCJlBKHSIyefvpp5s2bl3QYIpFS4hCJ0eTJk5k9u6nfdiySHdTjEInRs88+i7snHYZIpFRxiMTo8MMP57DDDks6DJFIqeIQidETTzzB2rVrM/ErQUWaTRWHSIweeughSkpKkg5DJFKqOERi9Pzzz/Pee+8lHYZIpFRxiMSoc+fOdOrUKekwRCKlikMkRrNnz1aPQ3KOKg6RGD322GM89dRTSYchEilVHJLzVm9OcdnMxayo3EHPgkOZOqIvPbp2bpVzz5s3j/Ly8lY5l0hrUeKQVjVq+iIWeGVi51+2qYqz7k3gY84fX7PPzf5WwPRRp7Z+HCIRUOJIyLk/L2Xpxqqkw4jJiqQDyHgLvJLC8b9LOoyENfw4OeHofF684exWjkXSEUviMLMZQIm7z23Cvj2Ar7j7nDhiCc9RAlzq7rvSOGaDu3eLK6Zc/Y9RXl6ecY3ggZNKWV5ZRU0ttMuDXgX5zBvXOvNfXFxMKpVi0aJFrXK+bJGJjxNpukxojp8DnB7nCdx9aDpJQ3LL3UN6075dHgDt2+Vx95DerXbuV155hZkzZ7ba+URaQ5xLVWPM7AfhOS5z9/fN7HpgGFALlAC/AsYDnc2sDNgG3E6Q0PKBYe7e4LfgmNkdwInAUUAX4Hp3f83MpgNfBDoB97v7LDOrAE50952NjHUZcA3QHnjO3W8HOprZ40APYDNwkbvvbuGctHlJ9zh276ll6JQ/JHDm/Zdl1OeQbBVn4ihz94lmNhj4bzObAFwMnBHePw/4PTCR4En9OTMbAwx393VmdgvwHeDHBzhHyt3PMbMvA4+b2RnAWUA/guR07sGCNLOjCJLXPwE7gZ+YWT5B4rrF3SvM7BXgFOAf6w0VFRXs3LmzWVfMXP3sGlZtzeUcpB5HU6jP0fTHyXGf7cDDQ7rHGEvymvt8EqfGlhPjTByvhj/LgHuB3sBxwPxwexfg+HrHrAUmm1kVcCzw+kHO8TKAu79rZt3cfbuZjQUeBT4DNOWLEHoCS9y9Orw9HsDM/u7uFeG2DcA+128WFhZSXV3drHXa0hxe283EteskexwDBgxgx44dvPHGG61yvmyRiY+TpGXTnMTZ49hbg58JLAEceBfo7+7FwAzgbaCmThxTgFHuPhJYB+Qd5Bx9AMysN7DWzD4H9HH3C4BvElQ6B0uOy4ETzaxjONaTZnYsQcUiOWDqiL70KsinfV4evQrymTqib6ud+6WXXmLatGmtdj6R1hBnxdHPzF4meAIe7e6rzGw+8Fr4JL2IoMJ4B/iRmf2JoEJYaGY7gI3AMQc5xynhmIcCVxBUBt3Cfske4D53/9jMGh3A3SvN7KdAqZnVAnPcfe2BjpHs0qNr51arMETagrza2ux8YR02xze4+8MJhVCbTaVla9Gc7OvBBx9kw4YN3HXXXUmHklH0ONlfhs5Jg6s+Gf8GQDN7Gjii3uZtwJ/THOd8YFwDd93v7s80MzyRA5ozZw5VVbn6Rk9pqzI+cbj7hRGN8xzwXBRjiTTVCy+8kHFXyoi0VCa8AVBERLKIEodIjO6//35mzZqVdBgikcr4pSqRbDZ//nz1OCTnKHGIxOi5555Tj0NyjpaqREQkLUocIjG677779M5xyTlaqhKJ0RtvvMH27duTDkMkUkocIjF66qmn1OOQnKOlKhERSYsSh0iMJk6cyJQpU5IOQyRSWqoSidFbb73Fhx9+mHQYIpFS4hCJUUlJiXocknO0VCUiImlR4hCJ0d13381DDz2UdBgikdJSlUiM3J1t27YlHYZIpJQ4RGI0e/Zs9Tgk52ipSkRE0qLEIRKj2267jV/+8pdJhyESKS1VicRozZo1bN26NekwRCKlxCESo+nTp6vHITlHS1UiIpIWJQ6RGN18881MmjQp6TBEIqXEIRKjVWs38r9/WUuvm59n4KRSVm9OJR2SSIspcYjEqPLkEXQ851r21NaybFMVl81cnHRIIi2mxCESoxWVOw54WyQbxXJVlZltcPdujdxXDFzt7kNbMP4dwAZ3f7gJ+w4Cerj7o2Z2JTDd3XfXG2sYcAewCJgB1AJLgGuBu4DLgZHuPre5MUvb9PEbM9ma2k2Xcy6jXR70LDg06ZBEWiznKw53n+vuj4Y3bwHaN7DbJHd/HJgETHD3M4E8YIi7TwCUMKRZ/vnz+bBnFwDt2+Vx95DeCUck0nItqjjMrDfBk2174EjgGncvq3P/K8B7wIkET8QXh3cdb2YvAEcBc9z9DjM7G7idIJnlA8PcfekBTn+BmX0X6Ax8z90Xmdl1wIXAocAHwAUE1cSJwDKgG1ACfKuRMfsApeHvLwDnAs80bTZE9vfhVy+ly+erANi9p5Zbn13CvHFnJxyVSMu0dKnqy8CN7v6OmQ0DRgFl9fYpc/erzWwMwSv+p4FDCJ682wOrCZaJvgwMd/d1ZnYL8B3gxwc498pw3C8Ds8zsn4GuwAB3rzGz3wN99+7s7lPN7FbgQEtkee5eG/6+HTi8sR0rKirYuXOn3txVj+ZkX8srq/a7rfnR46QhmTgnRUVFDW5vaeJYC9xqZtXAYUBD35H5cvizDBgS/r7E3T8CMLOP64w12cyqgGOB1w9y7lcB3P1dM+sWJotdwP8Nx/g80CHNv6emzu+HAY1+VkRhYSHV1dWNTmxbVV5erjmp4+PXb2VrajdHDLiSdnnQqyBf84MeJw3JpjlpaY9jMnC7u48A3iFYjqqvT/jzdODd8PfaBvabAoxy95HAukbGqutUADM7CVhtZv8EfMvdLwauJ/jb6o9Rw4H/5j+HzXuAbwALDxKDyAEN/NLRHNaxHe3z8uhVkM/UEX0PfpBIhmtpxTEb+I2ZbQH+RtDnqG+kmY0DdgCXACcdYKyFZrYD2Agcc5Bzf8HMXgY6AlcB7wM7zGxvpbK+gTEWAs+bWf86S1J13QhMMbNPA+XAkweJQeSApj78q6x6JSnSFC1KHO4+iaA5Xn973Utxb3b39+rcfiX8t8++7j4ujfPe0chd5xzkuBEHuX8poM6liMgBZPSn45rZ08AR9TZvc/chDe3fAuPM7MPwktz6MdwDDCK4GkskLddeey1btmzh8cf3e2iJZK1YE4e7F7fw+AsjCuVA57iD4Kquxu6fAEyIOw7JTZ06dSKV0udTSW7J6IpDJNvdd999GXeJpUhL5fw7x0VEJFqqOERidOWVV7J161Z+/etfJx2KSGSUOERi1LVrV2prG7ryWyR7KXGIxOgnP/mJehySc9TjEBGRtKjiEInRqFGj2Lp1K888ow9ZltyhxCESo+7du9OpU6ekwxCJlBKHSIzuuusu9Tgk56jHISIiaVHFIRKj4cOHs23bNubMmZN0KCKRUeIQiZGZUVlZmXQYIpFS4hCJ0a233qoeh+Qc9ThERCQtqjhEYjR06FA+/PBDnn/++aRDEYmMEodIjE4++WQ2bdqUdBgikVLiEInR+PHj1eOQnKMeh4iIpEUVh0iMvv3tb7N9+3ZefPHFpEMRiYwSh0iMvva1r7Fx48akwxCJlBKHSIz+8z//Uz0OyTnqcYiISFpUcYjE6Pzzz6eqqoqXX355n+2rN6e4bOZiVlTuoGfBoUwd0ZceXTsnFKVIelRxiMTo61//Ov369dtv+6D7X2XZpir21NaybFMVl81cnEB0Is2jikMkRt///vcb7HGkdu3Z5/ayTVWtFZJIi7X5isPMRprZajMbF94ea2YTw98vMrP39t4WicrxR+XTLi/4vV1ecFskW7T5xBF6HHjIzP4HuHbvRnd/ElDSkGb7xje+wZVXXrnf9ruH9KZ9mDnat8vj7iG9Wzs0kWZT4vjEIcBM4MdJByK547zzzqN///77bb/12SXs3lMLwO49tdz67JLWDk2k2dTjCLn7FuBFMxvZlP0rKirYuXOnrtGvR3Oyr/79+zc4J/V7Gss2VbWpedPjZH+ZOCdFRUUNblfiaKbCwkKqq6sbndi2qry8XHNST0NzcvxRm1heWUVNbdDj6FWQ36bmTY+T/WXTnGipSiRGAwYMYPTo0fttnzqiL70K8mmfl0evgnymjuibQHQizaOKQyRGF198MevXr99ve4+unZk37uwEIhJpOSWOetx9RtIxSO644oorMm7dWqSltFQVGLb3fRx1mdlFwPgE4hERyVhtvuIIK4wZjdz3JPBka8YjuaW4uJhUKsWiRYuSDkUkMm0+cYjEaeTIkaxbty7pMEQipcQhEqORI0eqxyE5Rz0OkRjt3r2b3bt3Jx2GSKRUcYjEaODAgepxSM5R4hCJ0eWXX87atWuTDkMkUkocIjEaPny4ehySc9TjEIlRKpWiuro66TBEIqWKQyRGgwcPVo9Dco4Sh0iMrrnmGvU4JOcocYjE6OKLL1aPQ3KOehwiMdq2bRvbt29POgyRSKniEInRkCFD1OOQnKPEIRKj733ve/ztb39LOgyRSClxiMTowgsvVI9Dco56HCIx+uCDD9iyZUvSYYhEShWHSIwuuugi9Tgk5yhxiMToxhtvZM2aNUmHIRIpJQ6RGJ133nnqcUjOUY9DJEYbNmygsrIy6TBEIqWKQyRGQ4cOVY9Dco4Sh0iMxo8fz+rVq5MOQyRSShwiMRo0aJB6HJJz1OMQidGaNWtYv3590mGkZfXmFAMnldLr5ucZOKmU1ZtTSYckGUYVh0iMLrnkkqzqcYyavogF/kkzf9mmKgbd/yp/vWtQglFJplHiEInRhAkTWLVqVdJhHNC5Py9l6caqRu9P7dpD4fjf/eP2CUfn8+INZ7dGaJKhYk0cZjYDKHH3uU3YtwfwFXefk+Y5ugDzgc3hz5fdfb+Xd2ZWDFzt7kPTGV+kJQYMGJDxPY76SWDgpFKWV1ZRUwvt8qBXQT7zxilRyCcyqcdxDnB6M447CVjp7gPdfWJDSUMkKStWrMi6d47fPaQ37dvlAdC+XR53D+mdcESSaVpjqWqMmf0gPNdl7v6+mV0PDANqgRLgV8B4oLOZlQHbgNsJEls+MMzdl9Yf2Mw+DUwGjjGzO4HjwvFWANOBj8MxhoWHHG9mLwBHAXPc/Y6GAjazQuA3wHrg88AL7v6jFs6DtEGjR48mlUpx7rnnJh1Kk9TvcezeU8vomYvV45B95NXW1sY2eLhU9Z67TzSzwcDlwATgUWBv7TsPuAb4GnCiu483szHAb919nZndAuS5+48bOUcx4RLU3qUxoBfwReAm4ExgE3AkQYL6KtAeWO3uRzYyZiGwGDCCJPYacK27/2nvPhUVFbVbtmzhkEMOac7U5KydO3dqTupYvHgxu3bt4vTTm1NMt46rn13Dqq27m7TvcZ/twMNDurf4nHqc7C8T56SoqCivoe2tUXG8Gv4sA+4FehNUBvPD7V2A4+sdsxaYbGZVwLHA62mecyrwQ2AuwRP/LeH2Je7+EYCZfXyQMf7i7n8P9/1/BEnkH4mjsLCQ6upqioqK0gwtt5WXl2tO6igqKsr4OSmtF1tr9DgyfU6SkE1z0ho9jlPDn2cCSwAH3gX6u3sxMAN4G6ipE88UYJS7jwTWAQ1mvQMYAix0968TLDn9MNyeTnlVZGadzaw98C/AX9OMQQR3Z+XKlUmHkZapI/rSqyCf9nl59CrIZ+qIvkmHJBmmNSqOfmb2MsGT9mh3X2Vm84HXzKwjsIigwngH+JGZ/QmYDSw0sx3ARuCYNM/5R2CmmU0gWJa6AfhMmmPsIkg6RwNPuvtf0jxehKuuuopUKsXgwYOTDqXJenTtrKuo5IBi7XFkq7DHUeLu/Q4+JG7VAAAId0lEQVSwW202lZatRXOyr7KyMioqKhg2bNjBd25D9DjZX4bOSWI9jkiY2dPAEfU2b3P3IS0Y80o+ueKqrpubO6ZIXaeddhpdunRJOgyRSGVN4nD3C2MY81GCK7wacqBqQ6RJlixZwooVKzLxlaRIs2VN4hDJRtdddx2pVIrzzz8/6VBEIqPEIRKje++9N+uuqhI5GCUOkRj17duX/Pz8pMMQiVQmfVaVSM556623Mv5DDkXSpYpDJEZjx44llUpx4YWRX9shkhglDpEY/eIXv2DFihVJhyESKSUOkRidfPLJdOzYMekwRCKlHodIjBYvXsw777yTdBgikVLFIRKjH/zgB6RSKb773e8mHYpIZJQ4RGL0wAMPqMchOUeJQyRGvXv3pn379kmHIRIp9ThEYlRWVsaf//znpMMQiZQqDpEY3XLLLaRSKX2suuQUJQ6RGD3yyCMsX7486TBEIqXEIRIjM6OmpibpMEQipR6HSIxKS0tZvHhx0mGIREoVh0iMbr/9dlKpFJdeemnSoYhERolDJEbTpk3j/fffTzoMkUgpcYjEqGfPnnz00UdJhyESKfU4RGL00ksvUVZWlnQYIpFS4hCJ0T333MMjjzySdBgikdJSlUiMZs2axbJly5IOQyRSShwiMerevTtVVVVJhwHA6s0pLpu5mBWVO+hZcChTR/SlR9fOSYclWUiJQyRGc+fOZfXq1RQVFSVy/lHTF7HAK/fbvmxTFWfduwCA/lbA9FGntnZoksWUOERiNHHiRFKpFFdeeWWs5zn356Us3di8ymaBV1I4/neN3n/C0fm8eMPZzQ1NcpASh0iMSkpKWLp0aeznacoT+8BJpSyvrKKmFtrlQa+CfOaNU0KQ9EV2VZWZzTCzQU3ct4eZndeMc3Qxsz+Z2TwzG29mDdbXZlZsZiXpjl/n+EIz+0NzjxfZq1u3bhQUFCQdBgB3D+lN+3Z5ALRvl8fdQ3onHJFkq6QqjnOAE4E5aR53ErDS3b8dfUgi0ZszZw5r1qzJuB7H7j21DJ0SvDZSj0PSlVdbWxvJQGY2AzgCOJQgIV3m7u+b2fXAMKAWKAF+BbwLdAauA7YBtxNUP/nAMHffr7Y3s08Di4BjgIeA48LxVgDTgY/DMYYBvYCfAZuAo4A57n7HAWKfAHwrjPsh4PfAfOAd4HPA2+5+Rd1jKioqards2cIhhxySxizlvp07d2pO6hgxYgQ1NTXMmjUr1vNc/ewaVm3dHcvYx322Aw8P6R7pmHqc7C8T56SoqCivoe1RJ4733H2imQ0GLgcmAI8CexdS5wHXAF8DTnT38WY2Bvitu68zs1uAPHf/cSPnKAaudveh4flKCJLEF4GbgDMJksWRBAnqq0B7YLW7H9nImKcADwBnhfv+JLy9CDiBILG9D/Rz9011Dq0tLy9P7JVkptKc7OuDDz5g6dKlnHbaaUmHklE9Dj1O9pehc9Jg4oj6neOvhj/LAAN6E1QG88N/XYHj6x2zFpgcJoL+QIc0zzkV2ArMJahgPg63L3H3j9w9VWdbQwxY5O573H2Xu99IUB2tcPct7l5DkIx0wbuk7cgjj6RLly5JhwHA1BF96VWQT/u8PHoV5DN1RN+kQ5IsFXWP41SCpHEmsARwgmWpb7h7rZndALwNFPNJ0poC9HL37WY2k0Yy3AEMARa6+51m9u/AD4GZBE/+TfEecI2ZtSOoOJ4nSEDRlGLSpj399NP87W9/y4hXkj26dtZVVBKJqCuOfmb2MjAWuMnd/0JQabxmZn8kqDbWEvQOhpjZUGA2sNDMXgcOI+hhpOOPwF3hea8GfpnOwe7+FkG18jrwGvA/gD7OVCIxefJkZs+enXQYIpGKrMfRBqnH0QDNyb62bduGu3PqqbpqqS49TvaXoXPS4ApQRr4B0MyeJrhCq65t7j6kBWNeSXDFVX03u/sbzR1X5EAOP/xwDjvssKTDEIlURiYOd78whjEfJbjCS6TVPPHEE6xduzYTX0mKNJu+j0MkRg899BAlJc3+EAORjJSRFYdIrnj++ed57733kg5DJFKqOERi1LlzZzp16pR0GCKRUsUhEqPZs2erxyE5RxWHSIwee+wxnnrqqaTDEImUKg6RGM2bN4/y8vKkwxCJlCoOkRh16NCBDh3S/fg1kcymikMkRjNmzGDdunXqcUhOUcUhEqMZM2bw29/+NukwRCKlz6pqpjfffFMTJyI5r0+fPvt9XpUSh4iIpEVLVSIikhYlDhERSYsSh4iIpEWX47aAmR1O8A2GnwE+DYzTd3sEzOwC4Dvu3tB3oOS88KuIHwS+QvCNkpe7+/vJRpUZzOxfgJ+6e3HSsSTNzDoA04BCoCNwj7s/l2hQTaCKo2XGAfPd/WxgJPCrZMPJDGZ2P/AT2vbj61vAIe7+NWA88LOE48kIZnYT8BhwSNKxZIjhwGZ3PxMYBDyQcDxN0pb/Y0fh58Aj4e+fAnYmGEsmKQOuSTqIhJ1B8F32uPsfgH9ONpyMsRyI/IvasthvgFvD3/OAjxOMpcm0VNVEZnYZcEO9zaPcfbGZdSNYshrb+pEl5wBz8oSZFScQUib5DLCtzu09ZvYpd8+KJ4a4uPtTZlaYdByZwt2rAMzsMOBJYEKyETWNEkcTuftUYGr97WZ2ElAC/Ke7l7Z6YAlqbE4EgA+Bul823q6tJw1pmJl1B54BHnT3x5OOpym0VNUCZvYlglJzmLu/kHQ8klFeBwYDmFk/4J1kw5FMZGZHAy8CP3T3aUnH01SqOFrmJwRNvvvNDGCbuw9JNiTJEM8AA82sjGDtelTC8UhmugXoAtxqZnt7Hd9w9+oEYzoofeSIiIikRUtVIiKSFiUOERFJixKHiIikRYlDRETSosQhIiJpUeIQEZG0KHGIiEha/j/6fhH9LT+rdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_all.trace,\n",
    "              varnames=['logit_psi', 'beta_psi_multi', 'beta_psi_ch',\n",
    "                        'alpha_bait', 'beta_fish_p', 'beta_fish_ch'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n",
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n",
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/pymc3/stats.py:218: UserWarning: For one or more samples the posterior variance of the\n",
      "        log predictive densities exceeds 0.4. This could be indication of\n",
      "        WAIC starting to fail see http://arxiv.org/abs/1507.04544 for details\n",
      "        \n",
      "  \"\"\")\n"
     ]
    },
    {
     "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>WAIC</th>\n",
       "      <th>pWAIC</th>\n",
       "      <th>dWAIC</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>var_warn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>all</th>\n",
       "      <td>2139.63</td>\n",
       "      <td>92.5</td>\n",
       "      <td>0</td>\n",
       "      <td>0.52</td>\n",
       "      <td>404.3</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>livebait &amp; persons</th>\n",
       "      <td>2143.83</td>\n",
       "      <td>76.05</td>\n",
       "      <td>4.2</td>\n",
       "      <td>0.37</td>\n",
       "      <td>407.37</td>\n",
       "      <td>45.5</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>livebait &amp; persons, v2</th>\n",
       "      <td>2150.41</td>\n",
       "      <td>75.38</td>\n",
       "      <td>10.78</td>\n",
       "      <td>0.09</td>\n",
       "      <td>407.41</td>\n",
       "      <td>47.58</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>livebait</th>\n",
       "      <td>2506.19</td>\n",
       "      <td>47.62</td>\n",
       "      <td>366.56</td>\n",
       "      <td>0.01</td>\n",
       "      <td>470.08</td>\n",
       "      <td>212.94</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>livebait &amp; camper</th>\n",
       "      <td>2506.52</td>\n",
       "      <td>50.21</td>\n",
       "      <td>366.88</td>\n",
       "      <td>0.01</td>\n",
       "      <td>468.77</td>\n",
       "      <td>211.73</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>basic</th>\n",
       "      <td>2562.35</td>\n",
       "      <td>42.24</td>\n",
       "      <td>422.72</td>\n",
       "      <td>0.01</td>\n",
       "      <td>483.78</td>\n",
       "      <td>231.85</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           WAIC  pWAIC   dWAIC weight      SE     dSE var_warn\n",
       "all                     2139.63   92.5       0   0.52   404.3       0        1\n",
       "livebait & persons      2143.83  76.05     4.2   0.37  407.37    45.5        1\n",
       "livebait & persons, v2  2150.41  75.38   10.78   0.09  407.41   47.58        1\n",
       "livebait                2506.19  47.62  366.56   0.01  470.08  212.94        1\n",
       "livebait & camper       2506.52  50.21  366.88   0.01  468.77  211.73        1\n",
       "basic                   2562.35  42.24  422.72   0.01  483.78  231.85        1"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "comparison = pm.compare({\n",
    "        m_basic: m_basic.trace,\n",
    "        m_bait: m_bait.trace,\n",
    "        m_bait_persons_2: m_bait_persons_2.trace,\n",
    "        m_all: m_all.trace,\n",
    "        m_bait_camper: m_bait_camper.trace,\n",
    "        m_bait_persons: m_bait_persons.trace,\n",
    "    },\n",
    "    method='BB-pseudo-BMA'\n",
    ")\n",
    "comparison"
   ]
  }
 ],
 "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.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}