14 Overdispersion, Zero-inflation, and offsets in the GLM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "93cb91e9-a035-40aa-a1fd-43f2b2a0c349",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import arviz as az\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import scipy.stats as st\n",
    "from toolz import partition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ce30c5f-1930-4abc-9eb1-b6d4a8fd9560",
   "metadata": {},
   "source": [
    "### 14.1 Overdispersion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1c6819d-68ae-44de-b5fb-cec6a9ba0d92",
   "metadata": {},
   "source": [
    "#### 14.1.2. Data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c5b1aa17-a4e0-47e5-88b7-88d958e1a372",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGQCAYAAAC+mlmCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfIklEQVR4nO3de7hddX3n8ffHBETkJpKm3CLeBmvtCM4pFUGLaFurorS1HR3F4tjJdKqVTLWMqVar1VJ9pprM2NpGUKhQWwWbVsc66igCKpSEi4LQ0lKUCIaAooC0qHznj72ih0NIdnL2OWv9znm/nmc/WXutddb67p39+332uuy1UlVIkqT2PKjvAiRJ0q4xxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEYZ4otUkqcm+YftTD8sSSVZOp91Sdp1tuvFxxBfIJKsTvKxGeOue4BxL6yqC6vq8Gnjb0jyzFnW8JQkn05yR5JvJflIksdPm35cknuT3Nk9NiX5YJKfnM16pYVqCO16J2o9M8lbdjBPkvx2V+/dSb6a5A+TPHjGcu7p+pE7klyV5LQk+879q2iPIb5wXAAck2QJQJIfBXYDnjRj3GO6eScqydHAJ4C/AQ4CHglcCXwuyaOmzXpTVe0F7A08GbgWuDDJMyZdk7QA9Nqu58D/AlYCL2XUB/w8cDzwwRnzvb2q9gaWAS9j1Fd8LslD57HWJhjiC8eljBr3Ed3zpwGfAf5hxrh/rqqbuq3iTQBJ3g+sAD7SbSGfOm25L+6+Ld+a5HXbWf/bgT+vqrVVdUdVfaOqXg9cDPzezJlrZFNVvQE4HXjbrr1saUGb13ad5MFJ1iS5qXus2bqVnOTkJBdNL67bNf+YJCuBFwOnduv6yMwXkuSxwG8AL66qL1TV96rqauCXgGclOX7m31TVv1bVpcDzgIczCnRNY4gvEFV1D3AJowZN9++FwEUzxt3v23pVnQR8FTihqvaqqrdPm3wscDjwDOANSX5s5t8n2RN4CvChbZT2QeBndlD+hxltWfgtW5qmh3b9OkZbvUcATwSOAl4/Rp3rgHMYbUHvVVUnbGO2ZwCbqurvZ/ztjYy+7D9gP1FVdwCfBJ66o1oWG0N8YfksP2zYT2XU2C+cMe6zO7nMN1XV3VV1JaPd40/cxjz7M/os3byNaTcDB+xgHTcBAfbbydqkxWA+2/WLgTdX1S1VtQV4E3DSbIqf5gC23UfA+P3E/hOqZcEwxBeWC4BjkzwMWFZV1wGfB57SjXsCO3/c7OvThr8D7LWNeb4J3AscuI1pBwK37mAdBwMF3L6TtUmLwXy264OAr0yb9pVu3CTcyrb7CBi/n/jGhGpZMAzxheULwL6MThz5HEBVfZvRN9iVjE4q+5cH+Ntdvp1dVd3VrfuXtzH5V4D/t4NF/AJwWbccSfc1n+36JuAR056v6MYB3AXsuXVCd0Ldzqzr08ChSY6aPjLJoYx24T9gP5FkL+CZjPZAaBpDfAGpqruBDcBvcd8P+0XduO19W98MPGo703fktcCvJnlVkr2TPKz7ucnRjHbJ3Uf3U5ODk7wR+DXgd2axbmnBmud2/QHg9UmWJTkAeANwdjftSuDHkxyRZA/uf8LqdtdVVf8I/ClwTpInJ1mS5MeB84BPVdWnZv5Nd6LdfwDWM9rj976deC2LgiG+8HwW+BFGDXyrC7tx22vspzFqvLcnec3OrrSqLgJ+DvhFRse3vgIcCRzb7f7b6qAkdwJ3Mjrz9ieA46rqEzu7TmkRma92/RZGXxi+CHwJuKwbtzWE3wx8CrhuRi0AZwCP79a1/gGW/0pGv0Y5m1Ef8HHgfEZnqE93apI7GO0+/3NgI/AU99bdX6p2eS+qJEnqkVvikiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEYZ4pIkNWpp3wXsrAMOOKAOO+ywvsuQ5tXGjRtvraplfdfRN9u/FqPttf/mQvywww5jw4YNfZchzaskX9nxXAuf7V+L0fbav7vTJUlqlCEuSVKjDHFJkhpliEuS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJapQhLqk3SfZLcm6Sa5Nck+TovmuSWtLcZVclLShrgY9X1QuS7A7s2XdBUksMcUm9SLIP8DTgZICquge4p8+apNYY4hq8JH2X8ANV1XcJC8mjgC3A+5I8EdgInFJVd02fKclKYCXAihUr5r1I9W8ofcAQ27/HxDV4VTXrxySXo4lZCjwJeHdVHQncBbx25kxVta6qpqpqatmyRX831kXJ9v/ADHFJfdkEbKqqS7rn5zIKdUljMsQl9aKqvg7cmOTwbtQzgC/3WJLUHI+JS+rTbwLndGemXw+8rOd6pKYY4pJ6U1VXAFN91yG1yt3pkiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQofycuSZozp512Gps3b+61hn333ZdVq1b1WgPA8uXLWb169USXaYhLkubM5s2bOfHEE3utoe/1b7V+/fqJL9Pd6ZIkNcoQlySpUYa4JEmNMsQlSWqUIS5JUqMMcUmSGmWIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHFJkhpliEuS1KjeQzzJf09ydZKrknwgyR591yRJUgt6DfEkBwOvAqaq6gnAEuCFfdYkSVIret8SB5YCD0myFNgTuKnneiRJakKvIV5VXwP+J/BV4GbgW1X1iZnzJVmZZEOSDVu2bJnvMiVJGqS+d6c/DHg+8EjgIOChSV4yc76qWldVU1U1tWzZsvkuU5KkQep7d/ozgX+pqi1V9V3gw8BTeq5JkqQm9B3iXwWenGTPJAGeAVzTc02SJDWh72PilwDnApcBX+rqWddnTZIktWJp3wVU1RuBN/ZdhyRJrel7d7okSdpFhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEb1/jtxSYtXkhuAO4DvA9+rqql+K5LaYohL6tvTq+rWvouQWuTudEmSGuWWuKQ+FfCJJAX8WVV574QFZHRfK1i7dm3PlQzHmjVrJro8Q1xSn46pqpuS/AjwySTXVtUF02dIshJYCbBixYo+apQGyxCX1Juquqn795Ykfw0cBVwwY551dHc3nJqaqnkvUrusahj/XUkGU8ukeUxcUi+SPDTJ3luHgZ8Fruq3KqktbolL6sty4K+746ZLgb+oqo/3W5LUFkNcUi+q6nrgiX3XIbXM3emSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHFJkhpliEuS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEYZ4pIkNcoQlySpUYa4JEmNMsQlSWqUIS5JUqMMcUmSGmWIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHFJkhpliEuS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEYZ4pIkNar3EE+yX5Jzk1yb5JokR/ddk6T5k2RJksuTfLTvWqTWLO27AGAt8PGqekGS3YE9+y5I0rw6BbgG2KfvQqTW9LolnmQf4GnAGQBVdU9V3d5nTZLmT5JDgOcAp/ddi9SivnenPwrYAryv2512epKHzpwpycokG5Js2LJly/xXKWmurAFOBe59oBls/0oyq8cklrF1OUPTd4gvBZ4EvLuqjgTuAl47c6aqWldVU1U1tWzZsvmuUdIcSPJc4Jaq2ri9+Wz/qqpBPIao7xDfBGyqqku65+cyCnVJC98xwPOS3AD8JXB8krP7LUlqS68hXlVfB25Mcng36hnAl3ssSdI8qarVVXVIVR0GvBD4dFW9pOeypKYM4ez03wTO6c5Mvx54Wc/1SJLUhN5DvKquAKb6rkNSf6rqfOD8nsuQmtP3MXFJkrSLDHFJkhpliEuS1Kjej4lr4TvttNPYvHlzrzXsu+++rFq1qtcali9fzurVq3utQdLCYohrzm3evJkTTzyx1xr6Xj/A+vXr+y5B0gLj7nRJkhpliEuS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJatRYIZ7kGw8w/pbJliNpiOwDpGEad0t8t5kjkuwGLJlsOZIGyj5AGqDtXrEtyYVAAXskuWDG5EOAz89VYZL6Zx8gDduOLrt6OhDgJ4Ezpo0vYDPw6TmqS9Iw2AdIA7bdEK+qswCSXFxV185PSZKGwj5AGraxboBSVdcm+VngCGCvGdPeMAd1SRoQ+wBpmMYK8STvAn4F+AzwnWmTai6KkjQs9gHSMI17K9IXAUdU1Y1zWYykwbIPkAZo3J+Y3QbcPod1SBo2+wBpgMbdEv8j4JwkpzE6I/UHqur6iVclaWjsA6QBGjfE3939+9wZ4wsv9iAtBvYB0gCNe3a611iXFjH7AGmYbJiSJDVq3J+Ybb304v1U1dMmWpGkwbEPkIZp3GPip894/qPAy4GzJ1uOpIGyD5AGaNxj4mfNHJfkPOB9wJsnXZSkYbEPkIZpNsfEvwb8+0kVIqk59gFSz8Y9Jv6fZ4zaE/hF4OKJVyRpcOwDpGEa95j4STOe38XoPsLvnGw5kgbKPkAaoHGPiT99rguRNFz2AdIwjbslTpLHMroJwsGMjoV9oKqum6vCJA2LfYA0PGOd2JbkBGAj8DjgG8DhwIYkz5vD2iQNhH2ANEzjbon/AfD8qvrM1hFJjgPeBfzt5MuSNDD2AdIAjfsTs0OAC2eMu6gbL2nhsw+QBmjcEL8CePWMcb/VjZe08F2BfYA0OOPuTv9vwEeSnALcCBzK6CcmHg+TFgf7AGmAxv2J2bVJfgx4MnAQcBNwSVV9dy6LkzQM9gHSMI17xbYjgNuq6qJp4w5Nsn9VXTlXxUkaBvsAaZjGPSZ+NrDbjHG7A++fbDmSBmrifUCSPZL8fZIrk1yd5E2zqlBahMY9Jr6iqq6fPqKq/jnJYZMvSdIAzUUf8G/A8VV1Z5LdgIuS/F1VeT12aUzjbolvSvKk6SO65zdNviRJAzTxPqBG7uye7tY9atdLlBafcbfE3wn8TZK3A/8MPBp4DfDWuSpMC0MSANauXdtzJcOwZs2avkvYVXPSByRZwuhKcI8B/riqLpltodJiMu7Z6e9JcjvwckY/LbkReHVVnTuHtUkaiLnqA6rq+8ARSfYD/jrJE6rqqunzJFkJrARYsWLFbFY3K1u/kPatyp0V+qGxb4BSVR8CPjSHtWgBqipWrVrFiSee2HcpvVu/fn3fJczKXPYBVXV7kvOBZwFXzZi2DlgHMDU11VuCTSI8kxjCmqhxj4lL0kQlWdZtgZPkIcAzgWt7LUpqzNhb4pI0YQcCZ3XHxR8EfLCqPtpzTVJTDHFJvaiqLwJH9l2H1DJ3p0uS1KixQjzJg5O8Ncn1Sb7VjfvZJK+c2/IkDYF9gDRM426JvxN4AvBifngxhqsZ3dlI0sJnHyAN0LjHxH8BeExV3ZXkXoCq+lqSg+euNEkDYh8gDdC4W+L3MCPwkywDbpt4RZKGyD5AGqBxQ/xDjH4K8kiAJAcC7wL+cq4KkzQo9gHSAI0b4r8D3AB8CdgPuI7RjQ+8daC0ONgHSAO0w2Pi3YUYXg/8j6pa1e1Cu7W8dqC0KNgHSMO1wy3x7gYFrwC+2z3fMunGm2RJksuTeLUmaWDmow+QtGvG3Z1+FvDrc1jHKcA1c7h8SbMz132ApF0w7k/MjgJ+M8mpjG5B+INv4VX1tNkUkOQQ4DmM7kv8W7NZlqQ5M2d9gKRdN26Iv6d7zIU1wKnA3g80wxDuJzyUewmD9xNWL+ayD5C0i8YK8ao6ay5WnuS5wC1VtTHJcdtZf+/3E/ZewlrM5qoPkDQ7Y9/FLMlyRrvUDgB+sFlaVe+dxfqPAZ6X5NnAHsA+Sc6uqpfMYpmS5sAc9QGSZmGsEE9yInA2o9+G/jijayY/AbgI2OUGXFWrgdXdOo4DXmOAS8MzV32ApNkZ9+z0twAvq6ojgbu6f1cCG+esMklDYh8gDdC4Ib6iqj40Y9xZwEsnVUhVnV9Vz53U8iRN1Jz3AZJ23rghfkt3PAzghiRHA48GlsxNWZIGxj5AGqBxQ/w9wLHd8DuBzwBXAn8yF0VJGhz7AGmAxv2J2dumDf95kvOBh1aVV1mTFgH7AGmYxv6J2XRV9dVJFyKpHfYB0jBsN8STXMi0yytui5dclBYu+wBp2Ha0JX76tOEAfwz8xtyVI2lg7AOkAdtuiM+81GKSd3j5RWnxsA+Qhm3cs9MlSdLAGOKSJDVqRye2HT9z/iRP5743P/j0XBQmqX/2AdKw7ejEtjNmPL+N+97soIBHTbQiSUNiHyAN2I5ObHvkfBUiaXjsA6Rh85i4JEmNMsQlSWqUIS5JUqMMcUmSGmWIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHFJkhpliEvqRZJDk3wmyTVJrk5ySt81Sa3Z0a1IJWmufA94dVVdlmRvYGOST1bVl/suTGqFIS6pF1V1M3BzN3xHkmuAg4GJhniSSS5u1oZQT1X1XYImxN3pknqX5DDgSOCSbUxbmWRDkg1btmyZ99qkIXNLXFKvkuwFnAesqqpvz5xeVeuAdQBTU1M7vQlZVaxatYoTTzxxtqUuCOvXr++7BE2QW+KSepNkN0YBfk5VfbjveqTWGOKSepHRweEzgGuq6h191yO1yBCX1JdjgJOA45Nc0T2e3XdRUks8Ji6pF1V1EdD/qdpSw9wSlySpUYa4JEmNWjS700877TQ2b97caw377rsvq1at6rWG5cuXs3r16l5rkCRNxqIJ8c2bN/f+O9G+1w/+RlSSFhJ3p0uS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEYZ4pIkNcoQlySpUYa4JEmNMsQlSWqUIS5JUqMWzf3E1Z/ly5f3fh/zM888k5NPPrnXGpYvX97r+iUtPIa45tzq1av7LoG1a9eyZs2avsuQpIlyd7okSY0yxCVJapQhLklSo3oN8SSHJvlMkmuSXJ3klD7rkSSpJX2f2PY94NVVdVmSvYGNST5ZVV/uuS5Jkgav1y3xqrq5qi7rhu8ArgEO7rMmSZJaMZhj4kkOA44ELtnGtJVJNiTZsGXLlnmvTZKkIRpEiCfZCzgPWFVV3545varWVdVUVU0tW7Zs/guUJGmAeg/xJLsxCvBzqurDfdcjSVIr+j47PcAZwDVV9Y4+a5EkqTV9b4kfA5wEHJ/kiu7x7J5rkiSpCb3+xKyqLgLSZw2SJLWq7y1xSZK0iwxxSZIaZYhLktQoQ1ySpEYZ4pIkNcoQlySpUYa4JEmNMsQl9SbJe5PckuSqvmuRWmSIS+rTmcCz+i5CapUhLqk3VXUB8I2+65BaZYhLktQoQ1zSoCVZmWRDkg1btmzpuxxpUAxxSYNWVeuqaqqqppYtW9Z3OdKgGOKSJDXKEJfUmyQfAL4AHJ5kU5KX912T1JJe7ycuaXGrqhf1XYPUMrfEJUlqlCEuSVKjDHFJkhpliEuS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSoxbFZVeTALB27dqeKxmGNWvW9F2CJGkCFkWIq21bv4QNYTlVNYFKJGkyFkWID6XjTTKYWlrieyZJ2+YxcUmSGmWIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHFJkhpliEuS1ChDXJKkRhnikiQ1yhCXJKlRhrgkSY0yxCVJapQhLklSowxxSZIaZYhLktQoQ1ySpEYt7bsASZpry5cvZ/369X2XwZlnnsnJJ5/caw3Lly/vdf2arFRV3zXslKmpqdqwYUPfZeySJLT2fmsYkmysqqm+6+hby+0f7AO0a7bX/t2dLklSowxxSZIaZYhLktQoQ1ySpEYZ4pIkNcoQlySpUYa4JEmN6j3EkzwryT8k+ackr+27Hknzx/YvzU6vIZ5kCfDHwM8DjwdelOTxfdYkaX7Y/qXZ63tL/Cjgn6rq+qq6B/hL4Pk91yRpftj+pVnq+9rpBwM3Tnu+CfipmTMlWQmsBFixYsX8VHb/GgazHC/bqAWimfbf1TGI5dj+NV3fIb6tT/P9PqFVtQ5YB6NrJ891Udtiw5Emrpn239XR16qlB9T37vRNwKHTnh8C3NRTLZLml+1fmqW+Q/xS4LFJHplkd+CFwN/2XJOk+WH7l2ap193pVfW9JK8E/i+wBHhvVV3dZ02S5oftX5q9vo+JU1UfAz7Wdx2S5p/tX5qdvnenS5KkXWSIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHFJkhqV1u7Mk2QL8JW+69hFBwC39l3EItX6e/+IqlrWdxF9a7z9Q/ufw1a1/r4/YPtvLsRblmRDVU31Xcdi5HuvIfBz2I+F/L67O12SpEYZ4pIkNcoQn1/r+i5gEfO91xD4OezHgn3fPSYuSVKj3BKXJKlRhnhPktyQ5IAJLevOSSxnoUtyfpL7naGa5OQk7+qjJi1e9gHza6G2f0N8JyVZ2ncNemBJlvRdgxY2+4DhWozt3xCfIcnvJrk2ySeTfCDJa7pvcH+Q5LPAKUlOSHJJksuTfCrJ8u5vfzrJFd3j8iR7JzkwyQXduKuSPHUb61yfZGOSq5OsnDb+ziRvTXJlkounreeRSb6Q5NIkvz9vb84AbOu96t6nNye5BDg6yRu69+aqJOuSZNoiXpLk8920o7ax/GVJzuv+/tIkx8zXa9Mw2AcMl+1/G6rKR/cApoArgIcAewPXAa8Bzgf+ZNp8D+OHJwX+GvBH3fBHgGO64b2ApcCrgdd145YAe3fDNwAHdMP7d/8+BLgKeHj3vIATuuG3A6/vhv8WeGk3/Argzr7fu3n8P7rfe9W9T78yc55u+P3T3sPzgfd0w08DruqGTwbe1Q3/BXBsN7wCuKbv1+xjXj9f9gEDftj+7/9wt9B9HQv8TVXdDZDkI9Om/dW04UOAv0pyILA78C/d+M8B70hyDvDhqtqU5FLgvUl2A9ZX1RXbWO+rkvxCN3wo8FjgNuAe4KPd+I3Az3TDxwC/1A2/H3jbrrzYRm3rvfo+cN60eZ6e5FRgT2B/4GpGnSvABwCq6oIk+yTZb8bynwk8ftqX932S7F1Vd0z8lWiI7AOGzfY/g7vT7yvbmXbXtOH/zeib208A/xXYA6Cq/pDRt/KHABcneVxVXcDoW9/XgPcneel9Vpgcx+iDc3RVPRG4fOvygO9W95WQ0Qd1+peuRffbwO28V/9aVd/v5tkD+BPgBd3/z3v44fsJ93/fZj5/ULf8I7rHwUNuwJo4+4CBsv1vmyF+XxcBJyTZI8lewHMeYL59GTVIgF/dOjLJo6vqS1X1NmAD8LgkjwBuqar3AGcAT9rGsr5ZVd9J8jjgyWPU+Tnghd3wi8d5YQvEOO/V1gZ7a/d/+IIZ0/8jQJJjgW9V1bdmTP8E8MqtT5IcMYnC1Qz7gOGy/W+DIT5NVV3K6FjTlcCHGTXCmf/JAL8HfCjJhdz3zjiruhMmrgTuBv4OOA64IsnljHZ/rZ2xrI8DS5N8Efh94OIxSj0FeEW3m27f8V7dgrDD96qqbmf07ftLwHrg0hmzfDPJ54E/BV6+jXW8CphK8sUkXwZ+fWLVa/DsAwbN9r8NXrFthiR7VdWdSfYELgBWVtVlfdclaX7YB6glnth2f+uSPJ7RbpmzbLzSomMfoGa4JS5JUqM8Ji5JUqMMcUmSGmWIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjTLEJUlqlCEuSVKjDHHtUJIbkjyz7zok9cM+YLgM8cbZuKTFzT5gcTPEJUlqlCG+ACV5WJKPJtmS5Jvd8CHTpp+f5PeTfC7JHUk+keSAadNPSvKVJLcled0O1nVckk0zxv1gyyDJUUk2JPl2ks1J3jFtvicn+XyS25NcmeS4Sb0H0mJmH7B4GOIL04OA9wGPAFYAdwPvmjHPfwJeBvwIsDvwGoDuFozvBk4CDgIeDhzCrlsLrK2qfYBHAx/s1nMw8H+AtwD7d+s/L8myWaxL0oh9wCJhiC9AVXVbVZ1XVd+pqjuAtwI/PWO291XVP1bV3Ywa1RHd+BcAH62qC6rq34DfBe6dRTnfBR6T5ICqurOqLu7GvwT4WFV9rKrurapPAhuAZ89iXZKwD1hMDPEFKMmeSf6s2x32beACYL8kS6bN9vVpw98B9uqGDwJu3Dqhqu4Cbpu27DunPVaMUc7LgX8HXJvk0iTP7cY/Avjlbjfa7UluB44FDty5VytpJvuAxWNp3wVoTrwaOBz4qar6epIjgMuBjPG3NwM/tvVJkj0Z7U4DoKr2mj5zkuXAntOeLwGWTZv/OuBFSR4E/CJwbpKHM+ok3l9V/2WnX52kHbEPWCTcEl8Ydkuyx9YH8DBGx8BuT7I/8MadWNa5wHOTHJtkd+DNbP9z8o/AHkmek2Q34PXAg7dOTPKSJMuq6l7g9m7094GzgROS/FySJV3tx00/+UbS2OwDFilDfGH4GKMGu/WxH/AQ4FbgYuDj4y6oqq4GXgH8BaNv5N8ENm1n/m8BvwGcDnwNuGvG/M8Crk5yJ6MTXF5YVf9aVTcCzwd+B9jC6Fv5b+NnUtoV9gGLVKqq7xokSdIu8BuPJEmNMsQlSWqUIS5JUqMMcUmSGmWIS5LUKENckqRGGeKSJDXKEJckqVGGuCRJjfr/70VFXFRsw6oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(456)\n",
    "n_site = 10\n",
    "x = pd.factorize(np.repeat((\"grassland\", \"arable\"), n_site))[0]\n",
    "eps = st.norm.rvs(0, 0.5, size = 2 * n_site)\n",
    "λ_OD = np.exp(0.69 + 0.92 * x) + eps\n",
    "λ_Poisson = np.exp(0.69 + 0.92 * x)\n",
    "\n",
    "C_OD = st.poisson.rvs(λ_OD, size = 2 * n_site)\n",
    "C_Poisson = st.poisson.rvs(λ_Poisson, size = 2 * n_site)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 6), gridspec_kw={'wspace': 0.35})\n",
    "ax[0].boxplot(list(partition(10, C_OD)), \n",
    "                    widths = 0.7, \n",
    "                    boxprops = dict(facecolor = 'gray', alpha = 0.5),\n",
    "                    medianprops = dict(linestyle='-', linewidth=2, color='black'),\n",
    "                    patch_artist=True, \n",
    "                    labels = ['grassland', 'arable'])\n",
    "ax[0].set_title('With OD', fontsize = 12, pad = 10)\n",
    "ax[0].set_xlabel('Land-use', fontsize = 12, labelpad = 10)\n",
    "ax[0].set_ylabel('Hare count', fontsize = 12, labelpad = 10)\n",
    "\n",
    "ax[1].boxplot(list(partition(10, C_Poisson)), \n",
    "                    widths = 0.7, \n",
    "                    boxprops = dict(facecolor = 'gray', alpha = 0.5),\n",
    "                    medianprops = dict(linestyle='-', linewidth=2, color='black'),\n",
    "                    patch_artist=True, \n",
    "                    labels = ['grassland', 'arable'])\n",
    "ax[1].set_title('Without OD', fontsize = 12, pad = 10)\n",
    "ax[1].set_xlabel('Land-use', fontsize = 12, labelpad = 10)\n",
    "ax[1].set_ylabel('Hare count', fontsize = 12, labelpad = 10);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4024698d-a5a2-41c7-9169-3253a88ecf0b",
   "metadata": {},
   "source": [
    "#### 14.1.3 Analysis using PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5e93355c-c659-4cb6-840a-61b7352b945b",
   "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: [eps, sigma, beta, alpha]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='14000' class='' max='14000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [14000/14000 00:33<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 3_000 tune and 4_000 draw iterations (6_000 + 8_000 draws total) took 48 seconds.\n",
      "The number of effective samples is smaller than 10% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_14_1_4:\n",
    "    alpha = pm.Normal('alpha', 0, tau = 0.001)\n",
    "    beta = pm.Normal('beta', 0, tau = 0.001)\n",
    "    sigma = pm.Uniform('sigma', 0, 10)\n",
    "    \n",
    "    eps = pm.Normal('eps', 0, sigma = sigma, shape = len(x))\n",
    "    λ = pm.math.exp(alpha + beta * x + eps)\n",
    "    C_OD_pred = pm.Poisson('C_OD_pred', λ, observed = C_OD)\n",
    "    trace_14_1_4 = pm.sample(4000, tune = 3000, target_accept = 0.95, return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "5a02c617-aa1d-44b5-8be5-03c45c51e6cd",
   "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.341</td>\n",
       "      <td>0.314</td>\n",
       "      <td>-0.233</td>\n",
       "      <td>0.950</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.004</td>\n",
       "      <td>3142.0</td>\n",
       "      <td>3102.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>1.073</td>\n",
       "      <td>0.376</td>\n",
       "      <td>0.348</td>\n",
       "      <td>1.766</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.005</td>\n",
       "      <td>3622.0</td>\n",
       "      <td>3535.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.439</td>\n",
       "      <td>0.226</td>\n",
       "      <td>0.058</td>\n",
       "      <td>0.830</td>\n",
       "      <td>0.008</td>\n",
       "      <td>0.006</td>\n",
       "      <td>580.0</td>\n",
       "      <td>447.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "alpha  0.341  0.314  -0.233    0.950      0.006    0.004    3142.0    3102.0   \n",
       "beta   1.073  0.376   0.348    1.766      0.006    0.005    3622.0    3535.0   \n",
       "sigma  0.439  0.226   0.058    0.830      0.008    0.006     580.0     447.0   \n",
       "\n",
       "       r_hat  \n",
       "alpha    1.0  \n",
       "beta     1.0  \n",
       "sigma    1.0  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_14_1_4, var_names = ['~eps'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48106e7f-851b-4413-9c3a-78661edba64a",
   "metadata": {},
   "source": [
    "### 14.2. Zero-inflation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74960a6a-62ac-451f-9fd6-31a7d21dd5ac",
   "metadata": {},
   "source": [
    "#### 14.2.2 Data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "766e768f-4e4c-4b41-aa1f-b5116c962046",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(12345)\n",
    "psi = 0.8\n",
    "n_site = 20\n",
    "x = pd.factorize(np.repeat((\"grassland\", \"arable\"), n_site))[0]\n",
    "w = st.bernoulli.rvs(p = psi, size = 2 * n_site)\n",
    "λ = np.exp(0.69 + 0.92 * x)\n",
    "C = st.poisson.rvs(w * λ, size = 2 * n_site)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d47ddb96-6833-4391-9b88-8113b0bd43cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.figure(figsize = (8, 6))\n",
    "#plt.hist(C, color= 'gray', edgecolor='black', alpha = 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89dcf857-a1a9-477d-89e8-c9a8ea4ffee9",
   "metadata": {},
   "source": [
    "#### 14.2.4 Analysis using PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8182a860-4c94-4809-ab3d-1490eee2ccb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "CompoundStep\n",
      ">NUTS: [psi, beta, alpha]\n",
      ">BinaryGibbsMetropolis: [w]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='12000' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [12000/12000 00:31<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 3_000 tune and 3_000 draw iterations (6_000 + 6_000 draws total) took 47 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as ZIP:\n",
    "    alpha = pm.Normal('alpha', 0, tau = 0.001)\n",
    "    beta = pm.Normal('beta', 0, tau = 0.001)\n",
    "    psi = pm.Uniform('psi', 0, 1)\n",
    "    R_lpsi = pm.Deterministic('R_lpsi', pm.math.logit(1 - psi))\n",
    "    λ = pm.math.exp(alpha + beta * x)\n",
    "    w = pm.Bernoulli('w', psi, initval = w, shape = len(x))\n",
    "\n",
    "    C_pred = pm.Poisson('C_like', w * λ, observed = C)\n",
    "    trace_ZIP = pm.sample(3000, tune = 3000, return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "76809327-df3e-4d71-a572-d0bf27451bb7",
   "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.814</td>\n",
       "      <td>0.196</td>\n",
       "      <td>0.453</td>\n",
       "      <td>1.179</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1887.0</td>\n",
       "      <td>2796.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>0.819</td>\n",
       "      <td>0.227</td>\n",
       "      <td>0.417</td>\n",
       "      <td>1.263</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>1971.0</td>\n",
       "      <td>2554.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>psi</th>\n",
       "      <td>0.761</td>\n",
       "      <td>0.073</td>\n",
       "      <td>0.619</td>\n",
       "      <td>0.892</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>2211.0</td>\n",
       "      <td>2452.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>R_lpsi</th>\n",
       "      <td>-1.206</td>\n",
       "      <td>0.428</td>\n",
       "      <td>-2.005</td>\n",
       "      <td>-0.429</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.007</td>\n",
       "      <td>2211.0</td>\n",
       "      <td>2452.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "alpha   0.814  0.196   0.453    1.179      0.005    0.003    1887.0    2796.0   \n",
       "beta    0.819  0.227   0.417    1.263      0.005    0.004    1971.0    2554.0   \n",
       "psi     0.761  0.073   0.619    0.892      0.002    0.001    2211.0    2452.0   \n",
       "R_lpsi -1.206  0.428  -2.005   -0.429      0.009    0.007    2211.0    2452.0   \n",
       "\n",
       "        r_hat  \n",
       "alpha     1.0  \n",
       "beta      1.0  \n",
       "psi       1.0  \n",
       "R_lpsi    1.0  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_ZIP, var_names = ['alpha', 'beta', 'psi', 'R_lpsi'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b75f342e-eabe-4265-a1c3-1fd40ae65913",
   "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: [psi, beta, alpha]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='6000' class='' max='6000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [6000/6000 00:17<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 33 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as ZIP1:\n",
    "    alpha = pm.Normal('alpha', 0, tau = 0.001)\n",
    "    beta = pm.Normal('beta', 0, tau = 0.001)\n",
    "    psi = pm.Beta('psi', 1, 1)\n",
    "    R_lpsi = pm.Deterministic('R_lpsi', pm.math.logit(1 - psi))\n",
    "    λ = pm.math.exp(alpha + beta * x)\n",
    "    \n",
    "    C_pred = pm.ZeroInflatedPoisson('C_like', psi, λ, observed = C)\n",
    "    trace_ZIP1 = pm.sample(2000, target_accept = 0.99, return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9129526c-6671-471a-9dec-e23cad7e1492",
   "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.821</td>\n",
       "      <td>0.195</td>\n",
       "      <td>0.452</td>\n",
       "      <td>1.184</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.004</td>\n",
       "      <td>1160.0</td>\n",
       "      <td>1602.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>0.812</td>\n",
       "      <td>0.226</td>\n",
       "      <td>0.402</td>\n",
       "      <td>1.254</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.005</td>\n",
       "      <td>1194.0</td>\n",
       "      <td>1742.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>psi</th>\n",
       "      <td>0.760</td>\n",
       "      <td>0.073</td>\n",
       "      <td>0.620</td>\n",
       "      <td>0.890</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1742.0</td>\n",
       "      <td>1746.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>R_lpsi</th>\n",
       "      <td>-1.196</td>\n",
       "      <td>0.423</td>\n",
       "      <td>-1.943</td>\n",
       "      <td>-0.376</td>\n",
       "      <td>0.010</td>\n",
       "      <td>0.008</td>\n",
       "      <td>1742.0</td>\n",
       "      <td>1746.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "alpha   0.821  0.195   0.452    1.184      0.006    0.004    1160.0    1602.0   \n",
       "beta    0.812  0.226   0.402    1.254      0.007    0.005    1194.0    1742.0   \n",
       "psi     0.760  0.073   0.620    0.890      0.002    0.001    1742.0    1746.0   \n",
       "R_lpsi -1.196  0.423  -1.943   -0.376      0.010    0.008    1742.0    1746.0   \n",
       "\n",
       "        r_hat  \n",
       "alpha     1.0  \n",
       "beta      1.0  \n",
       "psi       1.0  \n",
       "R_lpsi    1.0  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_ZIP1, var_names = ['alpha', 'beta', 'psi', 'R_lpsi'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd583047-0488-438f-9e96-11ecfd290855",
   "metadata": {},
   "source": [
    "### 14.3 Offsets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48699adf-e8d2-4306-ab5e-382f9158828c",
   "metadata": {},
   "source": [
    "#### 14.3.2 Data Generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5d41eab5-8d62-40aa-8751-1076accfb607",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123456)\n",
    "n_site = 10\n",
    "A = st.uniform.rvs(2, 3, size = 2 * n_site)       # Areas range in size from 2 to 5 km2\n",
    "x = pd.factorize(np.repeat((\"grassland\", \"arable\"), n_site))[0]\n",
    "linear_predictor = np.log(A) + 0.69 + 0.92 * x\n",
    "λ = np.exp(linear_predictor)\n",
    "C = st.poisson.rvs(λ, size = 2 * n_site)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6ebe6f89-7edf-4802-bf32-59d485f36f10",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.figure(figsize = (8, 6))\n",
    "#plt.hist(C, color= 'gray', edgecolor='black', alpha = 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c1e6277-b2cd-4e4a-a2d5-49ee00f4dc49",
   "metadata": {},
   "source": [
    "#### 14.3.4 Analysis using PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "956ae0a1-4563-4d6e-94f2-cc0f5a076026",
   "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, alpha]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:14<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 3_000 tune and 5_000 draw iterations (6_000 + 10_000 draws total) took 29 seconds.\n",
      "The number of effective samples is smaller than 25% for some parameters.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as Offset:\n",
    "    alpha = pm.Normal('alpha', 0, tau = 0.001)\n",
    "    beta = pm.Normal('beta', 0, tau = 0.001)\n",
    "    λ = pm.Deterministic('λ', pm.math.exp(1 * np.log(A) + alpha + beta * x))\n",
    "    C_pred = pm.Poisson('C_pred', λ, observed = C)\n",
    "    \n",
    "    trace_Offset = pm.sample(5000, tune = 3000, return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "f722228c-6620-4392-bef5-bdb6a7feabdc",
   "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.701</td>\n",
       "      <td>0.123</td>\n",
       "      <td>0.45</td>\n",
       "      <td>0.915</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2078.0</td>\n",
       "      <td>2812.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>0.964</td>\n",
       "      <td>0.145</td>\n",
       "      <td>0.70</td>\n",
       "      <td>1.241</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2031.0</td>\n",
       "      <td>2828.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "alpha  0.701  0.123    0.45    0.915      0.003    0.002    2078.0    2812.0   \n",
       "beta   0.964  0.145    0.70    1.241      0.003    0.002    2031.0    2828.0   \n",
       "\n",
       "       r_hat  \n",
       "alpha    1.0  \n",
       "beta     1.0  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_Offset, var_names = ['alpha', 'beta'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32557bb7-34a0-4ea3-9f90-d9b1adbcc994",
   "metadata": {},
   "source": [
    "**EXERCISES**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed38a45b-0bc6-46da-a5db-f8b4d0cef418",
   "metadata": {},
   "source": [
    "**1.** *Estimating a coefficient for an offset covariate*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "69f440a4-ecad-4b00-9e38-3bd868ceec6e",
   "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, alpha, omega]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:25<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 3_000 tune and 5_000 draw iterations (6_000 + 10_000 draws total) took 41 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as Offset1:\n",
    "    omega = pm.Normal('omega', 0, tau = 0.001)\n",
    "    alpha = pm.Normal('alpha', 0, tau = 0.001)\n",
    "    beta = pm.Normal('beta', 0, tau = 0.001)\n",
    "    λ = pm.Deterministic('λ', pm.math.exp(omega * np.log(A) + alpha + beta * x))\n",
    "    C_pred = pm.Poisson('C_pred', λ, observed = C)\n",
    "    \n",
    "    trace_Offset1 = pm.sample(5000, tune = 3000, return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "f40dbec6-0bff-4370-a0c3-d21bf713216c",
   "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>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>omega</th>\n",
       "      <td>1.037</td>\n",
       "      <td>0.275</td>\n",
       "      <td>0.527</td>\n",
       "      <td>1.570</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2909.0</td>\n",
       "      <td>3011.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.643</td>\n",
       "      <td>0.375</td>\n",
       "      <td>-0.025</td>\n",
       "      <td>1.387</td>\n",
       "      <td>0.007</td>\n",
       "      <td>0.005</td>\n",
       "      <td>2803.0</td>\n",
       "      <td>2790.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>0.973</td>\n",
       "      <td>0.144</td>\n",
       "      <td>0.703</td>\n",
       "      <td>1.243</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4835.0</td>\n",
       "      <td>4190.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "omega  1.037  0.275   0.527    1.570      0.005    0.004    2909.0    3011.0   \n",
       "alpha  0.643  0.375  -0.025    1.387      0.007    0.005    2803.0    2790.0   \n",
       "beta   0.973  0.144   0.703    1.243      0.002    0.001    4835.0    4190.0   \n",
       "\n",
       "       r_hat  \n",
       "omega    1.0  \n",
       "alpha    1.0  \n",
       "beta     1.0  "
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_Offset1, var_names = ['omega', 'alpha', 'beta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "51b2694d-3a28-4b8e-9855-5e04f1c7baa3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAAF1CAYAAADssDCjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAd7ElEQVR4nO3deXRV1f338feGBNGCaHFChYIgOJAouLTyK0FlEJSHCkbUKj62Dg9YqGJxrBEcUNEOBFw4oDKpoHVEpUVFQcBqlSW1UZBWbR8FFBCEUiGQxPP8QczDEMhFMt2d92utu8g959yd77nHk3669z7nhCRJkCRJilm9mi5AkiSpqhl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkpSyEkB1C+EsIYW3pv9mly5MQwqAQwtIQwsoQQt8QwrkhhM9DCMtDCL22aqNlCOHPIYQ1IYQPQgjdtlp3VAjhnRDC+hDC5NLtfl667qzS7deHEJaEEPpW9/5LSl8GHkkpCSE0AF4AJgEHlv77QulygE5AG+DXwINA99L31wH3lrZRH3gJmA4cBFwBTAsh/LC0jceB54AfAq+VtvGddUAu0KS0zSlbfU6Sdin4LC1JqQghdAEeSpKk3VbL/gFcBrwBHJskyaIQwt7ABuCYJEkWb/V+H+B44MEkSbK3auMlYDLwLvA+0DRJkuLSdYuBu5MkmVROPQuBXydJMrsq9ldSXDJqugBJaeNQ4LPtlv1foFnpzysBkiTZGEIAWLXd+x8AzYGjQwhrt2ojE3gZOARY8V3YKbX0ux9CCJ2Be4Cj2NI73YgtPUGSVCEDj6RULQdabLesBfDFbrSxDFiYJMlJ268IIbQEDg4hZGwVeg7fapNHgVuAqUmSFJX28ITd+N2S6jDn8EhK1dtAwxDCZSGEjBDCZUDD0uWp+iuwdwjh0hBCZghhrxDCqSGEQ5Ik+TfwMTCsdN0AtswB+k5jYDVQEkL4GdC+UvZKUp1g4JGUkiRJNgNnsWXOzurSf88qXZ5qG8XA/wJ+ypaeoaXAtfz/v0UXAucAa4AebJkbtKl03a+A8aW/Owd4a8/2SFJd4qRlSbXWd5OikySZW9O1SEpv9vBIqjVCCJ1CCD8qHTK7gi3DWO/UdF2S0p+TliXVJs2BPwL7A0uA3CRJCmu2JEkxcEhLkiRFzyEtSZIUPQOPJEmKXkVzeBzvkiRJ6WKnNyO1h0eSJEXPwCNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Aw8kiQpegYeSZIUvYyaLkBS7VdQUJDSdllZWVVciSR9P/bwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKi5314pDos1fvrSFK6s4dHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ978Mjqdqlev+frKysKq5EUl1hD48kSYqePTySKo13bpZUW9nDI0mSomcPjxQhe1okaVv28EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0nLUuqtbxBoaTKYg+PJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl63odHUtrbnYeles8eqW4y8EhpwiegS9L355CWJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDyS0taCBQu44IILOOGEE7j33ntruhxJtZj34ZGUtpo1a0ZeXh6vvfZayp9J9X5G3qBQiouBR9IuZWdnM2TIEObPn8+6deu46667ePzxx1m0aBEHH3ww+fn5NGzYkPXr13P33Xfz6aefUlhYSLdu3Rg8eDAA99xzDwsXLqSoqIjDDjuM22+/nX333Zd3332X0aNH07p1az744AP23Xdf8vPz2X///VOq7bDDDuOwww5jzpw5lJSUVOXXICnNOaQlqUIHH3wwU6ZMoV+/fgwcOJDLLruM5557joyMDF555RUAfvvb33LKKacwdepU/vjHP1JQUMC8efMAuPzyy5k2bRpPP/002dnZTJgwoaztf/zjH1xyySU899xztGzZkqeffhqAefPm0b9//3Jfb7/9dvV/CZLSmj08kirUo0cPAI466igOP/xwWrVqVfZ+6dKlALzxxhssWrSI8ePHA7Bhwwb+9a9/kZOTw9y5c3nyyScpLCxk06ZNNG/evKztI488sqy99u3b88EHHwCQk5NDTk5Ote2jpLgZeCRVqEGDBgDUr1+/7GeAevXqlQ0lJUnC/fffz4EHHrjNZ5ctW8aYMWN44oknOOigg5g1axZTp07doe3t25s3bx5jx44tt55hw4Zx8sknV87OSaoTDDySKsUpp5zChAkTuPbaa6lXrx4rVqwghMA333xDgwYN2H///SkqKuLZZ59NqT17eCRVJufwSKoU119/PRs3buScc86hX79+DBs2jPXr19O2bVtycnLo27cvgwYNom3btpX2OxcvXkz37t159NFHmTp1Kt27dy8bEpOkrYUkSXa1fpcrJVWfVC+nVuXwsnQpLYWdrbCHR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A4+kSpUkCUVFRSRJUtOlSFIZA4+kSlVcXExeXh7FxcU1XYoklTHwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A48kSYpeRk0XINV1BQUFNV2CJEXPHh5JkhQ9A48kSYqegUeSJEXPwCNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToZdR0AZJUGxUUFKS0XVZWVhVXIqky2MMjSZKiZw+PVEVS7SGQJFU9A48k7QGHvqT04JCWlCZKSkoYMGAA2dnZ9OvXr2z59OnTyc7O3uaVn5+/y3bGjRvH6aefzgknnECvXr2YPXs2ALNmzaJ79+7k5OTwxBNPlH3mpZde4qKLLiJJkj3ah+nTp9O3b1+OO+44srOz+fDDD/e4vYr2fWff284sWbKEjh07kp2dzaRJkwDYuHEjQ4YM4aSTTuLiiy/mq6++AqCoqIh+/foxf/78PdoPSVXPHh4pTUyePJl//vOfO10/cOBAjjjiCABat2690+0eeOABHnzwQTp37sygQYP48ssvKSoqAmDUqFEcf/zxNGnShHvuuYe+ffuyefNm8vPzue+++wgh7HbdS5Ys4YgjjiAzM5ONGzdy8sknU1xczGeffZbS5zds2MCXX35Ztm/l2dW+V/S9ba2oqIibbrqJzMxMiouLy5b/6U9/4r333uOqq65i0qRJTJs2jV/96ldMnDiR1q1b07lz55Tal1Rz7OGR0sDHH3/Mfffdx9VXX73TbTp27EjXrl0544wzaNu2bbnbbN68mcmTJ3PooYcyZswYevfuzS9/+UtOP/10YEtPxlFHHcXxxx9PcXExmzdvZvTo0fTq1WunbZZn48aNTJ48mdzcXM4991w2bdoEwPnnn88NN9zAQQcdtMvPl5SU8Oabb3LjjTdy6qmn8uyzz+5y+53teyrf29YeeOABvvnmG/r377/N8g0bNtCoUSM6derEAQccwIYNG1i6dClTp07luuuuS6ltSTXLHh6plisuLiYvL4/evXuTk5Oz0+0GDRoEwJFHHsmIESPKnTPy+eefU1hYSHFxMd26dWPdunUcc8wx3HnnnbRq1Yr+/fszduxYAHr16sWnn37KX/7yF5577rkK6/zmm2+YNWsWL774Iu+88w5ffPEFvXv35owzzqBRo0Yp7evixYt56aWX+POf/8zGjRs57bTT+P3vf0+nTp12+bny9j3V7+07H374IRMnTmT8+PG8884726zr0aMHDz/8MH379mWfffYhLy+PkSNHcvnll1cY3iTVDgYeqZabNm0ay5Yt47bbbmPFihXAlhD0xRdf0KxZM5o3b851111H8+bN+dvf/sYjjzzCDTfcwIwZM3Zo67uellWrVpGXl8eaNWsYN24ct912GxMnTmTo0KH06dOHzZs306ZNG84//3yuv/56Xn31VR5++GEALr/8cn7605/u0PaoUaOYPn06Xbt2pUuXLowePZrMzMyU93P58uWcd9557LfffgwbNoxevXqx11577fIzu9r3ir637d1yyy10796dZs2asX79egD+85//sG7dOg455BBmzJjBJ598QosWLXj77bf5+uuv6dmzJ0OGDGHJkiW0a9eOkSNHst9++6W8z5KqT6hgEuKezVCU6rDKuiz97rvv5vHHH99heePGjXnzzTd3WJ6Tk8O6det499132WuvvSgpKSFJEjIyMli3bh1dunShadOmvP7662zatIkTTzyRAw88kNdee22bdiZMmMDChQu5++676dy5MyNGjODbb7/l9ttv580332TvvffeZvv33nuPxx57jLlz55KZmcmAAQM488wzadWq1Q41XnLJJSxYsIBp06Zx7LHHAlBYWMikSZOYMWMGK1asICcnh549e9KlSxcaNmyY0ne19b7n5+fv8ntLkoSSkhLq1atHvXr1+MlPflIWdLZ24YUXcv3115e9X79+Pf369SM/P5/XXnuNWbNm8eSTT9K/f3969uzJlVdeWW5tXqUlVYudTjS0h0eq5fr168eJJ54IwJo1a7jtttto1qwZN910EwAjR46kUaNGtGrVivfff79smOq73pHevXuzfPly3njjDfbff3969OjBK6+8wvjx41m7di0AJ5100ja/c9myZUyePJlp06YBkCQJ8+fPJ0mSstf2OnbsSMeOHVmzZg2DBw/mr3/9K+PHj6dt27ZMmTKFvffem0WLFrF48WJWrVoFwOzZs/noo4/Izc2lYcOGDBo0iEGDBvH+++8zY8YM7rjjDoYPH86wYcN2mFdT0b5X9L298MIL3HzzzZx//vn85je/4a677iqbqDxz5kxmzpzJWWedtcOVXWPGjKFr1660b9+eWbNmsXr1ap555hm+/vprvv32290/wJKqhYFHquXatm1bNhF32bJlAOyzzz506dIFgDZt2vDkk08ydepUGjZsSLdu3bjmmmt22t5NN91EkiQ88sgjNGzYkD59+mzTgwFw1113cfHFF3PooYcCcM011/Dggw8SQuD6669nn3322Wn7jRs3pkWLFowcOZKVK1duM7Q2Z84cHnjggbL348ePByA3N3ebNo477jiOO+44rr32WubPn192Fdn2drXvFX1v29t6+UcffQRsueJr60nQBQUFzJkzp2xO0wUXXMDChQu59957OfrooxkwYMBOvxdJNcshLamK1NU7LRcVFZVN6t2dOTyxc0hLqhY7HdLysnRJkhQ9A48kSYqegUeSJEXPwCNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A48kSYqegUeSJEXPwCNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoZNV2AlG4KCgpqugRJ0m6yh0eSJEXPwCNJkqJn4JEkSdFzDo8kVYNU535lZWVVcSVS3WQPjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A48kSYpeRk0XINUWBQUFNV2CJKmK2MMjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A48kSYqegUeSJEXPwCNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Aw8kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A48kSYpeRk0XIEnafQUFBSlvm5WVVYWVSOnBwCNJtcjuBBlJqXNIS5IkRc8eHkkpeeaZZ5g6dSrFxcW0a9eOW2+9lb333ptly5bRp08fWrduDUCSJLRo0QKAOXPmMGbMGOrXr8+tt97KscceC8Dw4cMZMGAAbdu2Lfd39erVi/Hjx5e1A5Cdnc17771HRkYGvXr1omHDhmRkZFBYWEi7du0YOHBgWXt5eXl06NCB3NzcqvxKJKURe3gkVeiTTz7hoYceYuLEiUyfPp0mTZowZcqUsvU//OEPeeqpp3jqqaeYNm0amZmZANx///3k5+czbNgwHn74YQAWLFjAvvvuu9Owk6qxY8fy9NNP8+KLL5KTk8MvfvELPv/88z1qU1K8DDySKvTJJ59w9NFHs++++wLwP//zP8ycObPCz2VmZrJp0yYKCwtp0KABRUVFPPTQQ1xxxRWVVlsIgb59+9K5c2eefPLJSmtXUlwc0pJUobZt23LPPfewYsUKDjjgAF599VW+/PLLsvVr167lvPPOA6B3795ly6+++mpuueUWMjMzGT58OFOmTCE3N5cf/OAHFf7OK6+8sqynKBXt27fnrbfe2o29klSXGHgkVahly5b8+te/ZujQoWRkZPDjH/+Y+vXrA3DggQfy8ssv07RpU1atWsWQIUNo0KABACeccAJTp04FYPny5XzwwQdceOGF3H777axZs4auXbvSp0+fcn/n2LFjd5jDsytJklTGrkqKlENaklJy5plnMm3aNB599FHatm1Lq1atAGjQoAFNmzYFtoSfM844g7Vr1+7w+T/84Q8MHTqUGTNm0LhxY373u99x3333sXHjxkqp78MPP6RNmzaV0pak+Bh4JKVk9erVAKxfv56JEycyYMCAsuVFRUUAbNy4kXnz5tG4ceNtPjtnzhxatWrFj370IwoLCwkhEEKgpKSEkpKSPaorSRKef/555s+fz7nnnrtHbUmKl0NaklJy4403smrVKjZv3sw555xDz549AVi4cCHjxo2jfv36lJSUcOqpp7Js2bKyzxUWFvLYY48xbtw4YMscn6uuuoqzzz6bs846i0aNGn2veq688krq169PYWEhRx11FBMnTuTwww/f8x2VFKVQwbi3g+KqM7zDbeUoKioiLy+PkSNH7takY1UdHy2hOiTsbIVDWpIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9LzTsqLnDQUlSfbwSJKk6Bl4JElS9BzSkqTIpTqs6zO3FDMDj9KWc3MkSalySEuSJEXPwCNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6PnwUEkS4FPVFTd7eCRJUvTs4VGtkur/w5QkaXfYwyNJkqJn4JEkSdEz8EiSpOgZeCRJUvQMPJIkKXoGHkmSFD0DjyRJip6BR5IkRc/AI0mSomfgkSRJ0TPwSJKk6Bl4JElS9Hx4qCRpt6T6kN+srKwqrkRKnT08kiQpegYeSZIUPQOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToeR8eVYtU79shSVJVsIdHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6Bh5JkhQ9A48kSYqegUeSJEXPOy1LkqpEqndYz8rKquJKJHt4JElSHWAPj/aIz8iSJKUDe3gkSVL0DDySJCl6Bh5JkhQ95/BIkmqUV3OpOtjDI0mSomfg0fe2aNEifvazn9GnTx8uvfRSVq5cWe52l1xyCb1796Z///7079+fF154oZorlVTbvfLKK5xzzjl06NCBZ555Zpfbvvzyy3To0IH27dtz7LHHsnz58mqqUunMIS19bwMHDuSSSy6hR48eTJ48mfz8fO68885yt7355ps5+eSTq7lCSemiTZs2jBo1igkTJuxyu88++4zBgwcza9YsWrZsybp162jQoEE1Val0ZuDR97JixQqWLFlC9+7dAcjNzaVbt241XJWkdHXEEUcAUK/ergce7r//fi699FJatmwJQJMmTaq6NEXCIS19L0uXLqV58+aEEABo1KgRmZmZrF27ttztR40axdlnn01eXh6rV6+uxkolxeSjjz5izZo1nHbaaXTo0IHhw4eTJElNl6U0YODR97KzPzDfBaCt3XnnnTz//PM89dRTNGvWjJtuuqmqy5MUqeLiYt566y1efPFF5s+fz9y5c5k8eXJNl6U04JCWylXRZaL//e9/+fe//02SJIQQ+OabbyguLi63e/mQQw4BoH79+lx00UWcfvrpVVKzpPi1aNGCLl260KhRIwDOPvtsFixYwM9//vOaLUy1nj08+l6aNm1Ky5Ytef311wF49tln6dq16w7bFRcXbzOENXPmTNq1a1dtdUqKS//+/Xn99dcpLi6mpKSEWbNmkZ2dXdNlKQ0YePS95eXl8dBDD9GnTx/mzJnD0KFDAVi5ciX9+/cHYPPmzQwePJjc3Fxyc3OZO3cud9xxRw1WLak2mj17Nt27d+eVV15h9OjRdO/ena+++gqAESNGMHv2bABOPfVUOnXqRFZWFscddxwtW7bk0ksvrcnSlSZCBZO9nAlWR/kUdH1fRUVF5OXlMXLkSDIzM2u6HEXEOy0rBTtOJC1lD48kSYqek5YlSWnBZ25pTxh46hiHqiRJdZFDWpIkKXoGHkmSFD0DjyRJip6XpdcxzuFRVUuShOLiYjIyMsp91IhU1Zy0XKft9I+Ok5YlVaoQgvffkVTrGHgkSXWWl7rXHc7hkSRJ0bOHR5IUFecqqjz28EiSpOgZeCRJUvQc0opEunfhrlu3rqZLkFSNmjRpUtMlqI7xPjy1XLoHmVRlZ2fXdAmSqtHf//73mi6hSng1V43zPjySJFU1L3OvvezhqeXqSg+PQ1pS3eKQVmoMRrvNHp7apK6EmN3hHz9JUlUy8KTALkpJktKbgacS2XMjSVLtVKcDjwFFkqS6YZeTlm+99daZwAHVV07KDgWW13QRqlQe07h4POPjMY1LrMfzqxEjRvQqd02SJGn3uuWWW5KarsGXx9SXx7MuvTymcb3q4vH00RKSJCl66Rp4bq3pAlTpPKZx8XjGx2Malzp3PCu68aAkSVLaS9ceHkmSpJQZeCRJUvRqdeAJIfQKISwJIXwcQrihnPUhhDC2dP3fQwgda6JOpSaF43lqCGFdCOFvpa/hNVGnUhNCmBBCWBlC+GAn6z0/00wKx9RzNI2EEJqHEGaHEBaHED4MIVxVzjZ15jyttYEnhFAfGAecARwD/CyEcMx2m50BHFn6+j/A/dVapFKW4vEEmJckyfGlr9uqtUjtrklA+fe72MLzM/1MYtfHFDxH00kxMCxJkqOBk4HBdfl/R2tt4AFOAj5OkuTTJEk2A08AZ223zVnAlGSLt4H9QgjNqrtQpSSV46k0kiTJXGDNLjbx/EwzKRxTpZEkSb5IkuS90p/XA4uBw7bbrM6cp7U58BwGfL7V+6XseKBS2Ua1Q6rHqlMI4f0Qwp9DCMdWT2mqIp6fcfIcTUMhhJZAB+Cv262qM+dpbX6WVihn2fbX0KeyjWqHVI7Ve8CPkiT5bwjhTOB5tnSzKj15fsbHczQNhRAaAc8AQ5Mk+c/2q8v5SJTnaW3u4VkKNN/q/eHs+NyPVLZR7VDhsUqS5D9Jkvy39Oc/AZkhhNr4LDelxvMzMp6j6SeEkMmWsPN4kiTPlrNJnTlPa3PgeRc4MoTQKoTQADgfeGG7bV4A/nfpLPOTgXVJknxR3YUqJRUezxDCISGEUPrzSWz573N1tVeqyuL5GRnP0fRSeqweARYnSfKHnWxWZ87TWjuklSRJcQhhCPAyUB+YkCTJhyGEQaXrHwD+BJwJfAxsAH5RU/Vq11I8nucAV4QQioGNwPmJtwKvtUII04BTgQNCCEuBEUAmeH6mqxSOqedoevkJcBFQEEL4W+my3wAtoO6dpz5aQpIkRa82D2lJkiRVCgOPJEmKnoFHkiRFz8AjSZKiZ+CRJEnRM/BIkqToGXgkSVL0DDySJCl6/w8cWwMfmeLdAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_posterior(trace_Offset1, \n",
    "                  kind = 'hist', \n",
    "                  hdi_prob = 0.95, \n",
    "                  ref_val = 1,\n",
    "                  ref_val_color = 'black',\n",
    "                  bins = 50, \n",
    "                  var_names = ['omega'], \n",
    "                  figsize = (10, 6), \n",
    "                  textsize=10,\n",
    "                  color = 'gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f77ae54-2648-4d1c-b509-f30566801cfd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}