Chapter_6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "14c5cae5-ac7e-47f7-aa5b-ff92e3446c15",
   "metadata": {},
   "source": [
    "# Chapter 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c1c2d631-3c49-4b7d-9209-a68bf6ff022c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import arviz as az\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc3d8c63-b73a-42c2-93f3-429e2bb2889f",
   "metadata": {},
   "source": [
    "## Selection/distortion effect"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bf061b91-9d64-4ce9-b403-6413df6f1c83",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(2)\n",
    "n = 200\n",
    "p = 0.1\n",
    "newsworthiness = np.random.normal(size=n)\n",
    "trusthworthiness = np.random.normal(size=n)\n",
    "score = newsworthiness + trusthworthiness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "00a1e30a-8a65-4d0f-8e9a-f8ef0eeafee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = np.quantile(score,1-p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c3af28b-292d-4bae-a758-0663026ae1cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected = score >= threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0fa56566-0ea0-48fa-a233-8dd27adf31dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        , -0.73834648],\n",
       "       [-0.73834648,  1.        ]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.corrcoef(newsworthiness[selected], trusthworthiness[selected])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a82fa3cf-110a-455a-993a-a50446202d54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'trustworthiness')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuO0lEQVR4nO2dfZxcZZXnf6eri+4MCURIT7pJkVTrB6FflLeMoxDBMOOSUTcTHZnZDCOLZtcZ1251ZRAy5Q7DjlkEtpwZo86uM3HdVioz4w4o8mIa7RZoFJgmBNIkhFEYNIC8iLyFhHQ6Z/+oe5tblaq67291f9/P5/l09e373HvqVvU5z3POec4jqgpCCCHZoyNuAQghhMQDDQAhhGQUGgBCCMkoNACEEJJRaAAIISSjdMYtgBuWLFmixWIxbjEIISRV3Hfffc+pak/98VQZgGKxiOnp6bjFIISQVCEijzc6ThcQIYRkFBoAQgjJKDQAhBCSUWgACCEko9AAEEJIRklVFhAhYfLt+5/Atdv24MkX9uOExQtw6fknY93py+IWi5DQoAEgBFXlv/H6ndg/OwcAeOKF/dh4/U4AoBEgbQtdQIQAuHbbnnnlb7J/dg7XbtsTk0SEhA8NACEAnnxhv6vjhLQDNACEADhh8QJXxwlpB2gACAFw6fknY0E+V3NsQT6HS88/OSaJCAmf2ILAItIN4A4AXYYc/09Vr4hLHpJtzEAvs4BIlogzC+g1AOep6isikgcwJSK3qurdMcpEMsy605dR4ZNMEZsB0Opu9K8Yv+aNxh3qCSEkImKNAYhITkR2AHgGwG2qek+Dcz4qItMiMv3ss89GLiMhhLQrsRoAVZ1T1dMAFAC8TUSGG5zzVVVdqaore3qO2M+AEEKIRxKRBaSqLwCYBLAmZlEIISQzxGYARKRHRBYbrxcAeDeAh+OShxBCskacWUB9AP6viORQNUT/pKo3xSgPIYRkijizgB4EcHpc9yeEkKyTiBgAIYSQ6KEBIISQjML9AEikcNMVQpIDDQCJDG66QkiyoAuIRAY3XWnMJf/0AG7Z+ZTrfj9//lV8Z8cTkd2PtB80ACQyuOlKsOz91X58Z8eTcYtBUgxdQCQyTli8AE80UPbtuOnKqwcP4ePXbcdTLx7AYVWMnncSiscfjb+8eRdePXgIb/i1o1C+4FT8+jHdNf127n2x4Tn/9tw+lL69E7985SByHYKvXHgGrv7ew/jpM6/gd/7mTvzeGcvw4bP7cfX3Hsbdj/4SBw8dxofesQIX/uYKqCquuPEhTP3rc+hb3I18juM+UoUGgETGpeefXBMDANp305Xb9zyLpcd04/98+G0AgJcOzOLir92Lv7toJY5f2IXvPvAkrt22B9decOp8n3++by82Xr8TB+cOY9niBThj+Rvmz/nkP+7Ax859E9YM9+LA7BxUgcvWnIK/u/NRfO3i3wAAVO75GRZ1deLGkVV47dAcPvi3P8Y5J/XgoSdfxKPP7sNtnz4Xz73yGn77C7fj91eeGMtzIcmCBoBERpY2XTm5dxE+d/NuXHXrbvzWKUtx7II8Hnn6FfzRlnsBAIcPK379mK758+997HlsvfdnODh3GEA1QP6NHz+ONy9dhFdeO4SnXzyANcO9AIDuup3LTO7812fx8C9exi0zvwAAvHxgFo89tw/3PPY81p56AnIdgqXHdOOsNx0f5lsnKYIGgERKVjZdeWPPQtz8iVWY3PMM/uf4Hpz1puNx0tKFuOG/nN3w/Bt3PInXDh2uOaYAXnntkON7qgJ/sXYI5765tmru5J5nXMtPsgGdgYSEwNMvHUB3Pof3n17AH5/zRuz4+Qt4ft9B3Pf4rwAAs3OH8cjTL8+f//yrBxte54kX9mNhVyd6j+3GtoeqI/vXDs1h/8E5LOzqrDEQ57y5B9+8+3HMGrOIR599Ba8ePITf7D8O333wScwdVjzz0gH8+Ke/DOttk5TBGQAhIfDwL17GVbfshoggnxN8bt0wch2CK2/chZcOzGLusOIjq/rx5qWLAADH/dpRDY3A4gV5AMBf/cFp+LPrd+KvbnsEnTnBV/7wTJzStwg5Eaz56zvwwTML+MjZ/dj7q1fxvi9OQaE47uij8NWLVuL8oV786Ke/xLu/cDtOWLwAZ6x4Q6TPgiQXqe7MmA5Wrlyp09PTcYtBSODUL5IDqgHyqz7wlky4zEi4iMh9qrqy/jhnAIQkgCwFyElyoAEgJCFkJUBOkgODwMQXW7duxfDwMHK5HIaHh7F169a4RSKEOIQzAOKZrVu3olQqYcuWLVi1ahWmpqawYcMGAMD69etjlo4QYgeDwMQzw8PD2Lx5M1avXj1/bHJyEqOjo5iZmYlRMkKIlWZBYBoA4plcLocDBw4gn8/PH5udnUV3dzfm5uZa9CSEREkzA8AYAPHMwMAApqamao5NTU1hYGAgJokIIW6gASCeKZVK2LBhAyYnJzE7O4vJyUls2LABpVIpbtEIIQ5gEJh4xgz0jo6OYvfu3RgYGMCmTZsCDwBzG0lCwoEzAOKL9evXY2ZmBnNzc5iZmbFV/m7TRs0Vsk+8sB+K17eR/Pb97nfCIoTUQgNAIsNMG928eTMOHDiAzZs3o1QqtTQC3EaSkPCgASCRsWnTJmzZsgWrV69GPp/H6tWrsWXLFmzatKlpH24jSUh40ACQyNi9ezdWrVpVc2zVqlXYvXt30z7Ntotsx20kCYma2AyAiJwoIpMisktEHhKRT8YlC4kGL2mjl55/MhbU7YDVrttIEhI1cc4ADgG4RFUHAbwdwMdFZDBGeUjIeEkbXXf6Mlz1gbdg2eIFEADLFi9giWRCAiK2NFBVfQrAU8brl0VkN4BlAHbFJRMJF69po36rZDKNlJDGJCIGICJFAKcDuKfB3z4qItMiMv3ss89GLlvQZL16ptu0Ub8wjZSQ5sRuAERkIYB/BvApVX2p/u+q+lVVXamqK3t6eo68QIrwkgZJ/ME0UkKaE6sBEJE8qsr/OlW9Pk5ZosBLGiTxB9NICWlOnFlAAmALgN2q+oW45IgSL2mQSSGtrqskpZF++/4ncPbnJ9B/+c04+/MTdEOR2IlzBnA2gA8BOE9EdhjtPTHKEzpprZ6ZZtdVUtJIGYsgSSQ2A6CqU6oqqvpWVT3NaLfEJU8UpLV6ZppdV0lJI2UsgiQRVgONkKiqZwZNml1XQDI2W2csgiQRGoCIWb9+feIVfj2m68q69WMaXFdJ4oTFC/BEA2XPkhYkTmJPAyXJJ62uqySRlFgEIVZoAIgt69evx6ZNmzA6Ooru7m6Mjo6G6rpKa8ZRK5ISiyCkBlVNTTvzzDOVuKdSqejQ0JB2dHTo0NCQViqVuEVqSqVS0f7+fp2YmNCDBw/qxMSE9vf3J1pm0h7csH2vnnXVD7R42U161lU/0Bu2741bpMAAMK0NdCpnAG1O2lI405xxRNJLVtN0pWoc0sHKlSt1eno6bjFSxfDwMDZv3lwTwJ2cnMTo6ChmZmZilKwxuVwOBw4cQD6fnz82OzuL7u5uzM3NtehJiHfO/vxEwyD9ssULcNfl58UgUbCIyH2qurL+OGcAKcWpnzxtKZxpXSxH0k1W03RpAFKIG7dO2hSqNePoG9/4Bvr7+3HeeefhxRdfTKzbKg2wDEVrklQyJEpoAFKIGz952lI4zYyjiy66CBdddBEAYGxsDGNjY5HGLtopEymr/m03ZDVNlzGAFOLWT75161Zs2rRpfvVxqVRK/GK0OGMX5gxry5YtWLVqFaamprBhwwZ88JKrcdfLx6duY5l2928HRTtvHNQsBkADEDNelHPaArteiDMY3Oj5/o/rxvHV+/cBnUfNH1uQz6Uil7//8pvR6L9cADz2+fdGLQ6JAQaBE4jXFM20uXW8EGfsolHg/Kaf5WqUP1BbzC3JPvas+reJPTQAMeI15z3qlblx0MzIrV69OnTffCPj8+SLBxqe++QL+xPvY8+qf5vYQxdQjDDnvTX17rHVq1fj5ptvPsI3H7TxaxQD+PD1e4Gjjzvi3GXGKDrpPvZ29m8TexgDSCBZ8OUHSZTPq974vOdjV+C7v1hUU9PfjAH813/cQR87STSeYwAicrSIdBiv3ywia429fIlPsuDLD5IoF7WtX78eMzMzmJubw8zMDK75+AVNi7nRx07SipMYwB0AukVkGYBxVLdx/HqYQmWFsH357ZTLDoQTGHbzjNadvgx3XX4eHvv8e3HX5efNu1DoYyeppVGFOGsDsN34OQrgM8brHXb9wmisBuqcdqyqGfR7CvJ67VxJkqQfNKkG6sQA3A/gHQDuBjBkHNtp1y+MRgPgnKGhIZ2YmKg5NjExoUNDQzFJFAxBlrZu12dESD3NDIBtEFhEzgVwCYC7VPVqEXkjgE+p6ifCmJG0ot2CwGHCDCN7+IxIVvAcBFbV21V1raH8OwA8F4fyJ+6IYyFV2mIOAwMDuPLKK2tkvvLKKxNbKI+QwGk0LbA2ABUAxwA4GsAuAHsBXGrXL4xGF5Bzoo4BhHW/MHczGxkZ0c7OTi2Xy7pv3z4tl8va2dmpIyMjgd2DkCQAHzGAHcbPCwGUAeQBPGjXL4xGA3AkrRRklFtBhuFPD9uIDQ0NaalUqnlG5u+EtBN+DMBDhtL/FoBzjWMP2PULo3kxAGnaD9ctScr06ejo0IMHD9YcO3jwoHZ0dHi+ZthB2jBkJiSJ+DEAnwDwBIBbUF3cuALAnXb9nDQAXwPwDIAZJ+e7NQBJUpBh4FVBhmEUw1DWYStoU2br8ygWi1ooFAK5PiFJwbMBaNgJ6PTSr8F1zgFwRlgGoN3T/LwoyDB99UFfN+zPr1KpaE9Pj/b29ur4+LiOj49rX1+f9vT0tM0gIWq4HiKZ+JkBLAWwBcCtxu+DADbY9XPaABTDMgDtPsX3oiDDVKpBzyyimMEVCgUtFos1MrfTICFKbti+V0/57K264rKb5tspn72VRiAB+DEAtwL4fdPvD6ATAS4EC9MAtPsMwIuCTJtRDNKoNLpWUM+jnWNNTjnrqh/UKH+znXXVD+IWLfP4MQD/Yvy833Jsh10/p83OAAD4KIBpANPLly939abbPQag6l7xpNEoBqFcm30XCoWC7+eRhe+ZE4oNlP+Ky27S4mU3xS1a5vFjAH4I4Hi8XhPo7QBut+vntIU5A1BN9sgsDtnSpqyCkreZ4SsUCr6vH6RRTbMPnTOA5OLHAJwB4C4ALxo/HwHwVrt+TlvYBiCpxKmIk2YUW8kTlHJt5erx+zyCciM58aEn2UAwBpBcfGUBGX7/IQDDAPJO+ji87lYATwGYNVYYtwwut5MBSKMrph6/irNSqWihUFAAWiwWdWxs7AhDGJRyDfN5B3VtuxF0GhRskg1UlvFrAM4C8IcALjKbk35Bt3YyAH4UWxJG8H5nMGb/YrGol19+uRaLRQWghUKhZjVuUMo1zBlXUNe286HTxUK84scF9A0APwLwFQCbjfZFu35htHYyAH4WcQVdE9+LMfGrmM3+IqLFYlEnJiZ0fHxci8WiFotFFZF5+YJ6v2EaziCubafgGWQlXvFjAHbD2Ds47tZOBsCrYnOreO1qBXlVrn5dM2b/rq4uLZfLNf3L5bJ2dXU5eg9xErRcdi4ezgCIV/wYgG8B6LM7L4rWTgZA1ZsCcaN47RS8n1F8q77N3pf1eFdXl5ZKJRWReRnNGUB/f//8DKD+/STFEITlUmrlQw8rBkC/ffvjxwBMAvgVgG0AbjSbXb8wWrsZAC+4Udp25/qNQzRSgCMjI46Ol0ol7ezs1CVLlujGjRubxgBMRkZGNJ/Pq4jo4OCglkql0MtbtzI2cQXxg1bWaQgsE//4MQDnNmp2/cJoWTYApkISEc3n81oqlWxHnnYK3q8Sa6Qkm12zq6vriOOmEWiVBWTeJ5/Pa7lcrnnPYZVudjK6T9uK6macduU2upUyQKDF4OJqWTUA9QqpVCrNj4ZbuULsFHwYboxmihGA5zx882/W/hMTEzo4OBiKwnViGNshjfeG7XsbKn8GltsP1wYAwJTx82UAL1naywBeatYvzJZVAxBmxlDQfnU3MwCnCrOjo0MHBwdr+h88eHDeAPql/hmIiI6NjdUcGxsbqzE2rVxgSYlT2NEsqMwZQPvBGUCKSeKagVaBXjexgUql4mgGYPr8zf7lclnz+XwoFUdzuZwuXry45lhfX98R+wTUy93qPfqV0c1n6DRO0CytdMVlNzEG0Gb4XQiWA3ACgOVmc9Iv6JZVA5A0d4PdzMJJFpCpMO1WAlvvVyqV5t0++Xw+kL17Gz3bnp4ezeVyNe+vt7fX1gAEUViuHrduOjdB3WYzgFP/YptneUky8RMEHgXwHKpbQ+40GvcEjpAkFHBrlMJppT62YDdita4EHh8fr3lPjZRmWDOZRrOrjo4OBdDQBdQqGA9Ax8bGaq4V9baYbtYKMAMoO/gxAD8BcLzdeVG0rBoA1Xhz4OsNkJm7b5XBGtB1YqxMxWZVwKZiizKbppGCNVcjW6mvHDo4OKjlcrnmvTXr52cG4Nb953a1MNcAZAO/6wAC2QLSb8uyAYiTeiU5NDSk5XK5RrGZis7piNVUbNbzTcUWpXurkcHq6enRnp6eI4yY1cVjym+VdWxsTAEEOlMLcwZAsoOXLKBPG20LgCkAGy3HPt2sX5gtyQYgSatUg6Z+FFqpVOa3UaxXdE5HrKZisypg60pgP8/P7WdhPb9QKGihUFAR0a6urppUW+t7M+W3vjdzltCq9Ibb70iYMQCSHbwYgCtatD9v1i/MllQDkAQffZg0GoWWSiXt6uqqUWaVSuUIpanaeMRqfWZjY2PzK4GPO+44LRQKrg2pl4Vyja7htHSGeW65XJ5PUW11Hz/fkbCygEh28OMCusDJsShaUg1A0rJ0GuFnhuJ0PYGZqWMqxmKx2LJkQ1BplFb5GvnmnX4WbhfOOV2Q5+TahISJHwOw3cmxKFpSDUDSywIEMfps5BKxUj86tmYMOTU2XpWktV8r37yd8XPyOXo1pEn/joSBm5kIZy3h4sUF9DtG7f+nAXzR0r4O4N5m/cJsSTUAVn+2qRyKxeIReeNxEeZKYpMgFJxZ6M2qXJ1co5VvvlKpaG9vrxaLRdv34OQ5eTUATt1ofu6RJNzEIhi3CB8vBuBUAP8RwOPGT7N9AMAbmvULsyXVAFQqFe3p6dG+vj4dHx/X8fFx7e3t1Z6enkT883pVzm4Mh18Xh5+Cb61888ViUfv6+mo+h2ZyOVng5mcm1agaan2solHVVKdupiThJhuJmUvh48kFZKwArrQ6J8qWVAOgqlooFOYzY8x/1KT4eL0qZ697D1iDuoVCwZHS8lPuoZVv3u3irFaj7yCMnN1iOmvNJLeB5iThZj0CdzoLHz8xgDsBHGV3XhQtyQYgyT5eryNXtwqvUrHf5L0Z5vOzKsnBwcGGG8M0u3cjxR1k8DXIz9hJ1dRmqaZJGFTYwRlAsvBjAMYA/AuA/9bO6wD8+l2TnuURRQ66qvfnENYMKsgUXTfvze55N7uWdQbQKKCdlEGFHYwBJAs/BqDhWgC7fmG0sAxAEEqiHdYCNFJabg2Hl1GymxiKV0MWRFDV6WfsJm22/hxrDMBPSmsSYBZQcvBdDhrAQgALnZ4fRgvLAAQ1ek9z9kZQBszLs3SaReVXRq+fT6OVwq2u4fQZNJMniEVthFjxMwMYBnC/kQ30OID7AAzZ9QujhWUAkuy/j4ogjaBbJe00/dOPjF7k8hrTCPL7lOZBBUkOfgzAjwCstvz+LgA/susXRkv6DCDNRKm06kfUHR0dunTpUhURXbJkieZyufmFZ9a+fmT0EtC2K1cdReDZlIVGgPjBjwF4wMmxKFqSYwBpJyojWJ8umsvl5lfprl27Vnt7e3XhwoXa0dFxRBkJN4up6nFrPMx7NStXbZbEbvSdCfL71OhaPT09nuolkezixwDcYGQAFY32WQA32PVz0gCsAbDH2HPgcrvzk5wFlHaiMoJWf39/f78C0KuuukqXLFmiABSArlixQgGoqn0tnkaLqfz45U3sylXb7XEc1PepXu5KpaJ9fX2OVjYTYuLHALzBKAGx3Wh/HcRKYGOR2U8BvBHAUQAeADDYqk+S1wG0A0EESe36mYrVTPsEoMuXL1cRmc/9Hx8fnzcArWrxNFtM1agks1sDV2+oJiZqy1WLSCRxo/qZy9DQkI6Pj9fcJ2vuSuIePwbgTXbneGkA3gFgm+X3jQA2tupDA5A8vCjWUqmkAHR8fFyXLVumixYt0lwup8uWLVMANRuwOxmlW2m1KYsbQ2W3sjkql1n9fTo6OnR8fLzmPllLWCDu8WMAbjdG6v8A4OMA3mLXx0kD8EEAf2/5/UMAvtSqTzsYgHZzNXkJrubzeV26dKmOj49rqVTSjo4OPfbYY+ddQMcee6yjjJtG9w5yW0ZrFlB9WYtW5SeC/Fzr71MsFrW3t7fm+u00A+B6gHDwtQ7AcNGcDaAE4GcAnnfSz+aajgwAgI8CmAYwvXz58nCfUsikLdjsxFh5ycwRkfnRfrFY1I997GPzfn83G8I0ep4IcGN2u8/LT76+29mINWuq0XaVSf0OuYErgsPDzwxgleGeucVICf0KgPV2/RxcN3MuoDDdBkHPLJwaqzAXfjmV0/q+rfv2NpPH6bNyWp7Ca4pp1Avakg5rAoWHHwNwCMA9ANYhwKJwADoBPAqg3xIEbrnALIkGwG41p/V4WAvOwphZuFnN6vTeUaxwdTJqd1rOwYxTWM8bGxs74vNyW8gujIFAOxgFVgUNDz8GYDGA9wK4GsAEgO8D+Eu7fk4agPcAeMSIMZTszk+aAXBSz8V63Mno1AthKBS3paDtlE9UPnM7eZo9q/rMoUKhoEuWLKmZAZRKJS0Wiw1nAG5KWds9W7fKPG2uxWZwBhAefmMAAwD+BMB1AB4DcLuTfkG3pBmAZsqkWY54oVAI5R81jJlF0EYlrtXW9cq0Ufpmo8whALpo0SLt7e2tKVAHoOFswc1mNq2ehRdl3i4r2RkDaEwQgXE/M4BHDf//RiMeENveAEkzAE5quluPm1sUBj1VD8ulEKSxiqPeUqP3YLqerDTKHMrn8/PZNqbhMFctN/rc3Gxn2erZevks26mWFbOAagnKKPoxAOc0OHa2Xb8wWtIMgNsZQFgjsrBcAE5dO04MWhyj1GalI/L5vG3mkJn6aa4D6Ovr06VLl84fq3++bt5fq/RSL8q8XWYA5EiCcov5MQDbnRyLoiXNALiNAYTpkw0rCNjqum4DwFE/k2bK1Iw7tMocKhaL2tPTo0NDQ/PpqlaXTqPMoiD2CvCizL0823YIGmeBoALjrg2AkaZ5CYCfw7ITGIC/QJsVg/ODmyygtBG0sgrzWVldNV1dXfM/G5WKcJLJ1NPTM59vb66+tb73RqNyJ+/D7pl5NZRu1xW0Q9A4C8Q2AwBwLqq7fz2F2t3APg3gpGb9wmxJNAD1pFXxN5LbTlkF4XsOQhmZ1zAzccrlshaLRV23bp3jYnGN3r95zJwBWPt5dbE4eWZBfYeaXYcuo/SQhBjACsvrDgDH2PUJq0VpALz8E6Z1ZNVMbruCZ0EokiCvYb2W+bvTctGtCPJzNWUMYhGcV5nbKWicBeLOAqoAOAbA0QB2AdgL4FK7fmG0qAyA13/4OEZWQYwWvQazg1CMXvcQbpTaab2W9fcgFFuQo3Kn+x/7odV3kTOA7OHHAOwwfl4IoAwgD+BBu35htKgMgNd/kKhHVkGNTFsFS51sbu5HMbodEbdK7Ww0A0iiYnNaYqIZTp55q+9iWmeqxDt+DMBDhtL/FoBzjWNtHQT2qshNl4P1n7PZYqAgCGokZ7cwKcyYhjkiti646uvrazoibpXaWR8DqN9RLCn4GSg4Vd5Ogs1pjFURb/gxAJ8A8ISxGEwArABwp12/MFrSZwAjIyPa2dmp5XJZ9+3bp+VyWTs7O3VkZCQUOYOaccQ9Imw0Im7mu7dL7bRmASVVsfkx3E77xv2ZkmThqxTEEZ2ATi/9/LY0xADSOANQjXdEWK/UK5XKvEEIYqVs0vCjnIOu0USygZ8ZwJ83anb9wmhJzwLyO7XPStZRPfVKfWhoSMvlco1St7qkWu3UlRa8Kud2MIAkevwYgEssrQTgxwC+ZtcvjJb0dQBe/zn9KHK7lbppGAHWv38z+GyVt75apllKoVgsOto9rF1oF6NPoiUwFxCALgA/dNsviJZ0A5Ck9FE7WZJmHKzyOFnBG8dIeGRkRLu6uhSAdnV1hRbbsSNpnx1JPkEagDcA+InbfkG0pBsA1ehdR82wy+wJexTpR0k5ka/RMzPLOoehGMMO8FOpkzDx4wLaCeBBoz0E4BkAI3b9wmhpMABeCGM028qoOLlf2ArcyTVa3b/+PVQqFe3r65t3CXmJDbS6Z1dXl5bL5Zrzy+WydnV1eZK//ly6dUiY+CoFYWnL4soA0jY2AGEogFZK3smOVH7kicI9Uy9jsVjU3t7emkqs4+PjWiwWHRdTa1bZdWioWg9oYGCg5jr79u1TAI6v1UwGBnZJ2HgyAAByAB5udU6UrV0NgGq0m7rbKRy/CimqFdHWZ2bW9LfKbt7Tiex2C8xERAFoPp+fd/s0mwG4fX6szUPCxs8M4DsAltudF0VrZwMQBs2Mit0I1a9CimNEa97TKrt5TyeyN3rP5sbu/f3985VF165dq52dnXrhhRc2jQG4fX6cAZCw8WMA7gDwMoAfALjRbHb9wmg0AMHRasbhVyHF4dM271ksFnV8fLzmnl5nAOaqYvO4NQtIRJoGgN0+vyTEABiEbm/8GIB7jb0BzPYuAPfY9Quj0QBEQxRB3DDwsz6gWZG5+pLYExMT83v/urlWknfoSoIBIuHixwA02hKyrauBZhWrEioUClooFCJX4G4yZ5rV/nFynUbn1B8bGRnRfD6v5XK5RjE6Ke+RphE1XVDtj2sDAOBjRgroPksa6IMAHgPwzWb9wmw0AOER9yjQzf3Nc0ulkhaLRS2XyzU7gnnN+GnUzzQCIqKDg4OJrTDqBwah2x8vBuBYAEUAW+tSQY9r1ifsRgMQHnGPAt3c3zzX2sf6uxd/v52PPi2jeS/E/dmT8Am0GmhcjQagliAVU9yjQDf3t+705WUHsLjfa9KIe/ZHwocGoM0I+p827lFgkmcAWaDdZzlZJ1EGAMAFRlmJwwBWOu1HA/A6QSuxMEeBToOySYwBENIOJM0ADAA4GcAPaQC8EYYbI4xRoFvF7jULyKzQ6XQnsFZZRIS0G4kyAPM3pwHwTFrcGHHUBfKT+8+ZAGlHUmsAAHwUwDSA6eXLl4f0eNJHUMrL6ajb6+yg1UwlqBmHVyPTqKKoWUU0DbMB+u2JUyI3AAC+D2CmQftdyzmcAfjAqwKwuj/MYmetjIgfY9NMORcKhcBG317dYdZ+5nscHx9XEUn8bICzF+KG1M4ArI0GwD9WxTE4ODgfQDUVR6ORsx83TjNFVSgUAnMN+ZkBmKt6zfIR1lW+SXSpmaTFBUiSAQ0AUdVaxWGOgK2Ko9HI2W/AuVHA1mzWmYvXILbX0bB1ly8R0WuuuaamwqcfecJ2zXAtA3FDogwAgPcD2AvgNQBPA9jmpB8NgH+sisM0BlbFEfQMwKQ+dXPp0qXa29tbk7rpN43VrdI1ZwBmATkAumTJEi0UCp7eo/V9hu2a4QyAuCFRBsBrowHwj1VxmMqqXC7r4OBgKDGA+vuaPyuVivb29s6Xb3a6c1eQdHR06NjYWM0agmuuuUYBeK75E5ViZgyAuIEGgKjqkYrD3PXKLg/er1ujUfkG6ybuACJXXkNDQ1osFmsMopkF1NXVFWhA2ny+Qa+xYBYQcQINAJknDsXht3xDGFQqFQWg4+PjNaPosbExz750a2DZfL7r1q3TfD7P0TqJDRoAEit25RvMzdejHs0WCgUtFos19/VjkKyB5X379mm5XFYAum7duprz6K8nUUIDQGKnWfmFkZERT26poGQKuqhe/QwAgA4ODtacx4wdEiU0AGSepPmOvQSmgyTsstqNtpCMewaQtO8ACRcaAKKqycwe8ZKamlQaZQGZM5pWzzxKhZzE7wAJFxoAoqrJzB/3sjgtqTRTrq1iHFEr5CR+B0i40AC0GWEUZ4sLL+UpkkT9Z+E2oB21Qk7id4CECw1AGxFGcba4FazbAnVJIYjRe9QKOanfARIeNABtRBjF2ZKkYNMUoAxCmUatkNPwHSDBQgPQRgRVnC0NCraepMkexOg9DoWctOdIwoUGoI3I6hQ+iSPXoD4LKmQSJjQAbUQSFWEUJNHwZfWzIOmCBqDNiGLE6PceQcuY1OwVjt5J0qEBIK5wO7JtlAoZ9Mg4iTOAMKFhIUFBA0Bc4UbZNjIWZiqnk/5OyZK7JUvvlYQPDQBxhRt3SyNj0dHREUoBNL+LrtJC1mY7JFxoAIgr3CigZgXQRMRRf6+08yg5qfEOkk5oAIgr3ChXrwXQ/NLOo+R2fm8kemgAiGucBiHtCqDV1/8Pygi08yi5nWc3JHpoAEioNDMWYSqydh8lMwuIBAUNAImFMJV0VKNkKmKSdmgASCyE7aYJWznTFUPaARoAEgtpctM0MiZpkp+QZtAAkFhIywi6mZwi0raBZpIdaABIbKTBh95spN/V1cUZAEk9zQyAVP8WLSJyLYB/D+AggJ8C+LCqvmDXb+XKlTo9PR2ydCSL5HI5HDhwAPl8fv7Y7Owsurq6UCwWsWXLFqxatQpTU1PYsGEDNm3ahPXr18coMSHOEZH7VHVl/fGOOIQBcBuAYVV9K4BHAGyMSQ5CAAADAwOYmpqqOTY1NYXBwUFs2rQJo6Oj6O7uxujoaEPlv3XrVgwPDyOXy+HEE0/EiSeeiFwuh+HhYWzdujXKt0KIcxpNC6JsAN4P4Don59IFRMLCT6zC2ndsbEz7+vq0t7dXx8bGEhvzINkCSY0BAPgugD9q8fePApgGML18+fKQHg8h3mMV1viB+doaJ2DMgMRNMwMQWgxARL4PoLfBn0qq+h3jnBKAlQA+oA4EYQyAJBFr/MB8DQDd3d2Ym5vD7Ozs/GtC4qBZDKAzrBuq6m/bCHQxgPcB+C0nyp+QpGLGD1avXl0TSxgYGABQjSWYrwlJEqEZgFaIyBoAnwFwrqq+GocMhARFqVTChg0bsGXLFlx22WW48MILoaq45pprMDk5OZ81REjSiMUAAPgSgC4At4kIANytqn8SkyyE+MLMCBodHcXu3btxwgknAAAuvvhiDAwMMGWUJJZY1gF4hTEAQghxT9LWARBCCIkZGgBCCMkoNACEEJJRaAAIISSj0AAQYmCt58MaPiQLxJUGSkii2Lp1K0ql0hFVPwEwhZO0LUwDJQTA8PAwNm/ejNWrV88fm5ycxOjoKGZmZmKUjBD/NEsDpQEgBM33A2ANH9IOcB0AIS1oth8Aa/iQdoYGgBC8Xs9ncnISs7Oz+OxnP4vzzz8fu3btYkCYtC0MAhOC2no+u3btQmdnJz7zmc/giiuuYECYtC2MARBSBwPCpN1gEJgQhzAgTNoNBoEJcQgDwiQr0AAQUkd9QNjc1KVUKsUtGiGBwiAwIXXUb/DCTV1Iu8IYACGEtDmMARBCCKmBBoAQQjIKDQAhhGQUGgBCCMkoNACEEJJRUpUFJCLPAng8bjlcsgTAc3EL4RLKHA2UORooM7BCVXvqD6bKAKQREZlulH6VZChzNFDmaKDMzaELiBBCMgoNACGEZBQagPD5atwCeIAyRwNljgbK3ATGAAghJKNwBkAIIRmFBoAQQjIKDUAEiMi1IvKwiDwoIjeIyOK4ZbJDRC4QkYdE5LCIJDaFTkTWiMgeEfmJiFwetzxOEJGvicgzIpKK/SVF5EQRmRSRXcZ34pNxy2SHiHSLyL0i8oAh85Vxy+QUEcmJyP0iclPY96IBiIbbAAyr6lsBPAJgY8zyOGEGwAcA3BG3IM0QkRyALwP4HQCDANaLyGC8Ujni6wDWxC2ECw4BuERVBwG8HcDHU/CcXwNwnqqeCuA0AGtE5O3xiuSYTwLYHcWNaAAiQFXHVfWQ8evdAApxyuMEVd2tqnvilsOGtwH4iao+qqoHAfwDgN+NWSZbVPUOAM/HLYdTVPUpVd1uvH4ZVeW0LF6pWqNVXjF+zRst8RkvIlIA8F4Afx/F/WgAoucjAG6NW4g2YRmAn1t+34uEK6a0IyJFAKcDuCdmUWwxXCk7ADwD4DZVTbzMAP4awGcAHI7iZtwSMiBE5PsAehv8qaSq3zHOKaE6nb4uStma4URmQkxEZCGAfwbwKVV9KW557FDVOQCnGTG3G0RkWFUTG3cRkfcBeEZV7xORd0VxTxqAgFDV3271dxG5GMD7APyWJmTxhZ3MKeAJACdafi8Yx0jAiEgeVeV/napeH7c8blDVF0RkEtW4S2INAICzAawVkfcA6AZwjIh8U1X/KKwb0gUUASKyBtVp3VpVfTVuedqIfwFwkoj0i8hRAP4DgBtjlqntEBEBsAXAblX9QtzyOEFEesxsOxFZAODdAB6OVSgbVHWjqhZUtYjqd3kiTOUP0ABExZcALAJwm4jsEJH/FbdAdojI+0VkL4B3ALhZRLbFLVM9RmB9BMA2VAOT/6SqD8UrlT0ishXAjwGcLCJ7RWRD3DLZcDaADwE4z/j+7jBGqUmmD8CkiDyI6kDhNlUNPa0ybbAUBCGEZBTOAAghJKPQABBCSEahASCEkIxCA0AIIRmFBoAQQjIKDQAhISEif2Z5XWxW/VNE/ruIpH1RHkkhTAMlJGCMhVMC4CVVXWgcKwK4SVWH45SNECucAZDUYIyid4vI3xk13sdFZIGIvElEvici94nInSJyilEI7DGpslhE5kTkHOM6d4jISSJyrmVh0/0iskhEviwia43zbhCRrxmvPyIim4zXnxaRGaN9yiLbHhEZQ7XcwBYAC4xrm7WfcvWyG32/LiIfNF7/m4hcKSLbRWSniJxiHD/a2EfgXkPW3zWODxnHdkh1v4mTjHNvNmrhz4jIH0T0EZGUQQNA0sZJAL6sqkMAXgDwe6huoD2qqmcC+FMAXzEKge1BdZ+AVQC2A3iniHQBOFFV/9U49+OqehqAdwLYD+BO4zVQrSxq1r1/J4A7RORMAB8G8Juo1sb/zyJyukW2r6jqkKp+GMB+VT1NVS9sIXsjnlPVMwD8rSEjAJRQLQ3wNgCrAVwrIkcD+BMAf2O8h5WoVkRdA+BJVT3VmHF8z9mjJVmDBoCkjcdUdYfx+j4ARQBnAfiWUfr3f6NaBgCoKvNzjHYVqobgN1AtDQAAdwH4goh8AsBio7TEnagaikEAuwA8LSJ9qJbE+JFxjRtUdZ9Rb/56vG4wHlfVu13K3ojrG5zz7wBcbrzHH6JaLGw5qiUl/kxELgOwQlX3A9gJ4N0icrWIvFNVX2whE8kwNAAkbbxmeT0H4DgALxgjbbMNGH+/A1Xl/DYAtwBYDOBdqCp5qOrnAfwnAAsA3CUip6jqE8Z5a4z+dwL4fQCvGJuhtGKfS9mbVeN9rcE5AuD3LO9xubFpTwXAWlRnL7eIyHmq+giAM1A1BJ8TkT+3kYtkFBoAknZeAvCYiFwAVAOwInKq8bd7UZ0dHFbVAwB2APhjGNtcisibVHWnql6N6qzgFKPf3QA+hdcNwJ8aP2H8XCciv2a4YN5v+Vs9s0YZ5SDYBmDUCDDDdDuJyBsBPKqqXwTwHQBvFZETALyqqt8EcC2qxoCQI6ABIO3AhQA2iMgDAB6CsS2kqr6G6o5hplvmTlSrsu40fv+UESR9EMAsXt+p7U4Anar6E1RjB8fh9VnDdlT39L0X1V2x/l5V728i11cBPGgJAvvhL1Hd1vBBEXnI+B2ozk5mDNfQMIAxAG8BcK9x7AoAnwvg/qQNYRooIYRkFM4ACCEko9AAEEJIRqEBIISQjEIDQAghGYUGgBBCMgoNACGEZBQaAEIIySj/HxXFAuaoLxL5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Figure 6.1\n",
    "plt.scatter(newsworthiness[~selected], trusthworthiness[~selected], lw=1, edgecolor=\"k\", color=(0, 0, 0, 0))\n",
    "plt.scatter(newsworthiness[selected], trusthworthiness[selected], color=\"C0\")\n",
    "plt.text(0.8, 2.5, \"selected\", color=\"C0\")\n",
    "\n",
    "plt.xlabel(\"newsworthiness\")\n",
    "plt.ylabel(\"trustworthiness\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da1b3140-531c-4759-bb8b-192ac9a315b8",
   "metadata": {},
   "source": [
    "# Multicolinearity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3988f7b8-814e-44e7-b5f9-fb1a3c42064f",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "height = np.random.normal(10,2, size=n)\n",
    "leg_proportion = np.random.uniform(0.4, 0.5, size=n)\n",
    "left_leg = leg_proportion * height + np.random.normal(0, 0.02, size=n)\n",
    "right_leg = leg_proportion * height + np.random.normal(0, 0.02, size=n)\n",
    "heights = pd.DataFrame({'height': height,\n",
    "                        'left_leg':left_leg,\n",
    "                        'right_leg':right_leg})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fda3cf20-4e52-49dd-94a2-a1532694c42a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, beta_right, beta_left, 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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 02:20<00:00 Sampling 4 chains, 31 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 155 seconds.\n",
      "There were 12 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 19 divergences after tuning. Increase `target_accept` or reparameterize.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m1:\n",
    "    alpha = pm.Normal('alpha', 10, 100)\n",
    "    beta_left = pm.Normal('beta_left', 2, 10)\n",
    "    beta_right = pm.Normal('beta_right', 2, 10)\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    \n",
    "    mu = pm.Deterministic('mu', alpha + beta_left * heights['left_leg']+ beta_right * heights['right_leg'])\n",
    "    sdr = pm.Normal('height', mu, sd=sigma, observed=heights['height'])\n",
    "    m1_trace = pm.sample(return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4ebec61f-8973-4415-b5a2-09452672902c",
   "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.954</td>\n",
       "      <td>0.291</td>\n",
       "      <td>0.407</td>\n",
       "      <td>1.487</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2177.0</td>\n",
       "      <td>1996.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>0.597</td>\n",
       "      <td>0.042</td>\n",
       "      <td>0.518</td>\n",
       "      <td>0.673</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1890.0</td>\n",
       "      <td>1775.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_left</th>\n",
       "      <td>3.756</td>\n",
       "      <td>2.138</td>\n",
       "      <td>-0.143</td>\n",
       "      <td>7.874</td>\n",
       "      <td>0.053</td>\n",
       "      <td>0.038</td>\n",
       "      <td>1599.0</td>\n",
       "      <td>1795.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_right</th>\n",
       "      <td>-1.741</td>\n",
       "      <td>2.145</td>\n",
       "      <td>-5.735</td>\n",
       "      <td>2.298</td>\n",
       "      <td>0.054</td>\n",
       "      <td>0.038</td>\n",
       "      <td>1593.0</td>\n",
       "      <td>1793.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  \\\n",
       "alpha       0.954  0.291   0.407    1.487      0.006    0.004    2177.0   \n",
       "sigma       0.597  0.042   0.518    0.673      0.001    0.001    1890.0   \n",
       "beta_left   3.756  2.138  -0.143    7.874      0.053    0.038    1599.0   \n",
       "beta_right -1.741  2.145  -5.735    2.298      0.054    0.038    1593.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "alpha         1996.0    1.0  \n",
       "sigma         1775.0    1.0  \n",
       "beta_left     1795.0    1.0  \n",
       "beta_right    1793.0    1.0  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m1_trace, var_names=['alpha','sigma', 'beta_left', 'beta_right'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e4e3adb9-bf91-4e20-9c88-6b5d5dc8b4cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAHHCAYAAAABJ3dMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhtElEQVR4nO3ce5hdVX3/8fcHIiZILIZIxIiAEkREKshFlHLpA5iqtCgg3n4VqhTFimKlLaiFijfUGlDQom3FVn6KiEKDIiKFIhopULwAUtIfECBiQgyRhCQEZP3+2HvoYZjbmczMWRPer+c5z+Tsvfbe373P5Hxmrb3OSSkFSZJqtlGvC5AkaTiGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpU0xpIckOSaJGuSLE/yr0lmDbPNPyQpSb46wmNslOSkJHcmWZvkZ0kOG6DdXya5J8mSJB9PslG/9XslWZlk2xEe99S2zikDrNu+XXdUx7Kj2mV9jwfbmr+d5PVJMsB+SpKPjKQePXkYVtIYSvIHwPeBFcBhwHuAfYErkjx1kG1eAbwFeKCLQ50GnAqcBfwR8BPggiSv6tjvHwKfAD4MvBd4F/CnHes3Br4AfKyUcmcXxx6NI4C9gVcBHwIeAr4GXJ5k2jgfWxuAJ/x1JGm9nAIsAg4tpTwCkOSXwHXA24DPdzZO8hTgHOCjwLEjOUCSLYH3A58opXy6XXxlku1pwum77bI/Ai4vpXyx3W7fdtm57fp3AtOAvn2Mp5+WUv6n4/m/JrkAuAD4JPDuCahBk5g9K2lsvYwmIB7pW1BKuR74DfDaAdqfCGxMd4HxSmAToP+Q4VeBFyfZrn2+CbCmY/1qYCpAOyx5GnBcKeXhLo49ZkopFwIXA8ck2bQXNWjyMKyksfU7YN0Ayx8Cdu5c0PaEPkj3gfGidn//02/5ze3Pndqf1wIHJtmtPdYRNMOF0ITjd0opV3Zx3E4bJ5nS+aAJ3W59F3gqsPso69CThMOA0tj6b5re1WOSbANsBfQPpC8A3xpFYMwAVpQnfrHn8o71AOcDhwI3tM+vBD6bZD/gNcCOXR6309r12LbTXe3PrcZof9pAGVbS2DoT+Go7m+2zNMHxReDR9gFAkrcAewAvGK9CSim/A16f5NnAlFLKXe09srOBD5ZSliR5D80kkM2AbwEnlFLWDL7Xx7yMphfZ6TnAt7sss282oN+orSEZVtIYKqWcl2RHmgkQH6B5Ez6fZrhrZ4AkmwGfAU4HHkqyebv5RsBT2ucPDjE0eD+weZL061319aiWdzYupfyq4+l7aXpFX0hyEM19q32BxcBlwMk0s/WGc0Pnfbn2vFaMYLv+tm5/3juKbfUk4j0raYyVUj4EzAR2AbYqpbwRmANc0zaZCTwT+BhN8PQ9tgZe3/771UMc4maa+zzP77e8717VLQNtlOQ5NPfI3llKeRSYSzMZ5KellPuAL7fLJtKracLzhuEa6snNnpU0DkopDwK/AEgyl+b+0Nva1b8GDhhgs6+323wUuGmI3X+P5v7Xm4G/61j+FuCmUsodg2x3BnBeKeW6jmVP6/j3ZvzvsNy4az/E/MfAmaWU1RN1XE1OhpU0hpLsSvNZpv9qF+1DMz39k6WUHwOUUtYCVw2w7VpgSSnlqn7LHwG+Ukp5W7v90iSfAU5KsrI91pHAH9K8+Q9U1yuBP+Dx98h+ALwnyXHAr2g+63Ru1yc9Mi9JMpNmOv1zaSZ4HAFcDpw0TsfUBsSwksbWOppvafgrmqG6XwLvKKV8eT32uTFPnBb+AWAVzeSIZ9HMQnx9KeWS/hu335xxFnBiKWVF3/JSyqVJTqa5T7UpcBEwXl9zdEH7cy2wlCZg3wB8c4BZjdITxN8TSVLtnGAhSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5T18fJ3Llzy7Jly3pdhiRNKjfccMNlpZQnfJOKU9fHjxdWkro34LeoOAwoSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmq3gYTVklKksO7aL9/u83M8axLkrT+NpiwkiRtuAwr6Ulg6cq13LDofpauXNvrUqRRmTRhlWRukh8muT/J8iSXJXnhIG23bYf43pTkmiRrk9ya5OABmv9+kmuTrE5yfZLdOvazRZKvJbknyZokNyc5etxOUhoHF924mIPnXc2H59/MwfOu5qIbF/e6JKlrKaX0uoYRSXJY+8+fA9OADwK7ATuVUtYlKcARpZRvJtkWuANYDLyv3eZdwNuB7Uspi5PsD1wJXAf8NXAvcCbw3HafJcls4E3A5cADwIHAWcAflVKuGKbkyXFhNSkdec6CIdffc/9qFq9oelFTp2zE/Hfvw5xZ01m4ZCWHfO4a1j7yKLM3n8pznrHpoPs4/9i9x7RmaYQy0MJJ07MqpVzYPhaWUn4OHA1sB+w5xGZfKKV8o5RyK/Ae4G7gnf3afKiUcmXb5sPAjsDs9piLSymfKqX8tJRyeynli8C3gDeO8elJ42arzacxZ9Z0AObMms5Wm0/rcUVS96b0uoCRSvJ84DRgL+CZNEG7EU1PaDCP/flZSnk0ybXATv3a/Lzj379qf24J3JNkY+BvgCNpAuypwCbAVaM+EWkMDNfrmXf5bZx5xUIA7l2xhoVLVj7Ws7p3xRoADn/p1pxw0A7jXqs0FiZNWAGXAPcAx9IM7z0C3EITHuvj4Y5/9w3d9fU43w/8JU2v7BfAKuBjNGEmVeuEg3Z4LIguunExR5yzgG1mbMqi5av5xGG7cOius3tcodSdSRFWSbagGZ47rpRyZbtsN4av/2XAv7ftQzNk+M0uDr0PML+U8q8d+9gBWNFN/VIvHbrrbF6+/RbcvXwNW8+YxpbTp/a6JKlrkyKsgPuBZcAxSe6mGZL7FE3vaijvTHIbTa/oOGAb4AtdHPc24Mgk+7THfzfNfbIbuytf6q0tp081pDSpTYoJFqWUR2nuG+0C3AScDXwIeGiYTf+GZjbgz4C5wGtLKfd0ceiPAP8JXApcDTwInNdV8ZKk9TZppq53o2Pq+h6llOt7VMaGd2ElafxN7qnrkqQnL8NKklS9DXIYsBJeWEnqnsOAkqTJybCSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVz7CSJFXPsJIkVc+wkiRVb9KFVZL9k5QkM3tdiyRpYqSU0usaupJkE2AGsKTUXXzNtUlSrTLgwrrf7yc1L6wmzNKVa7l7+Rq2njGNLadP7XU50voYMKyqHQZMsm+SnyRZleS3Sf4zyc4DDQMm+bMkdyVZnWR+kuOSlI71pya5Kclbk9yZ5MEkX06ySdv27iS/SfKZJBt1bPeWJNclWZlkaZILksye6GshDeWiGxdz8Lyr+fD8mzl43tVcdOPiXpckjbkpvS5gIEmmABcD/wS8GXgKsBvwuwHa7g38I3AS8G1gP+BjA+x2W+BPgNcAs4ELga2Ae4GDgR2BbwA/atcBbAKcAtwKzAROB74G7LveJymNgaUr13Lyt3/Bxe96BXNmTWfhkpUccc4CXr79FvawtEGpMqyApwObA/NLKf+vXXYrQJJZ/doeD3y/lHJ6+/y2JHsAx/RrtzFwdCnlt8BNSb5HE2yzSynrgF8m+RFwAG1YlVL+uWP725O8s233nFLKPWNxotL6uHv5GmZvPo05s6YDMGfWdLaZsSl3L19jWGmDUuUwYCllOXAucFmS7yR5X5LnDtJ8R+A/+y27doB2d7VB1WcJcFsbVJ3Ltux7kmS3JBcnWZRkJXB9u2qwWqQJtfWMady36iEWLlkJwMIlK1m0fDVbz5jW48qksVVrz4pSytFJzgDmAn8MfDTJocBDo9zlw/0PMciyjQGSPA24DPgB8H+ApTRDgT+kGR6Uem7L6VM59ZAXccQ5C9hmxqYsWr6aUw95kb0qbXCqDSuAUsrPgJ8Bpye5FHgr8MV+zW4F9ui3bM8xOPyONOF0cinlDoAkrxuD/Upj6tBdZ/Py7bdwNqA2aFWGVZLtgGOBfwMWA88DdgG+MEDzzwLXJDkRuIhm8sNrx6CMu2h6cX+R5GzghcBpY7BfacxtOX2qIaUNWpX3rIDVwA7ABcBtwFeA82hm4z1OKWUBzWSK44GfA4e27dauTwGllPtoenKHArfQzAp83/rsU5I0Ohvkh4KTzAMOLKW8uIdlbHgXVpLG34AfCq5yGLBb7RDg5cAq4EDgHcDJPS1KkjRmNoieVZLzgf2B3wPuAM4BzuzxdwdO/gsrSRPP7wacYF5YSere5PpuQEmS+hhWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqGVaSpOoZVpKk6hlWkqTqjWlYJbkqyVljuc+xkKQkObzLbf48yV1JHk1y6jiVJkkagap6Vkn2b4NlZo/reAZwNvApYDbw6VqDWJKeDKoKq4psA0wBLiml3FtKWdXrgiSNztKVa7lh0f0sXbm216VoPYxHWE1JcmaS+9vHp5JsBJBkkySnJ7knyeok1yV5ZbtuW+DKdh/3tT2sc9t1c5P8sN3f8iSXJXnhaAtMMjvJ1ztq/E6SOe26o4Ab26a3d9SxH/Cu9nlp65VUsYtuXMzB867mw/Nv5uB5V3PRjYt7XZJGaco47PPNwLnA3sAuwJeAe4HPAF8Gng+8CbgHeBUwP8kewE3AYcCFwIuA5cCadp9PA84Afg5MAz7YbrdTKWVdN8Ul2ZQmFH9ME0DrgPcDP2gD8Py23u8BewJ3t3XsANwKnNzu6r5ujisN5chzFvS6hHFzz/2rWbyiN72aqVM2Yv6792HOrOksXLKSQz53De89/6cTdvzZm0/lOc/YdMKOV4Pzj917XPY7HmF1L3B8KaUAtybZAXhfkouBNwLbllLuatueleRA4NhSynFJlrfLl5ZSlvXtsJRyYecBkhwNPEATJtd0Wd8bgABHtzWS5FhgKfCaUso3kvymbXtfKeXXbZt1wOq+55Lqt9Xm05gzazoAc2ZNZ6vNp3HHsgd7XJVGYzzC6id9IdBaAJwG7EMTErck6Wz/VODfh9phkue3+9gLeCbN8OVGwHNHUd9Lge2Alf3q2JSm1ydNuPH6a7QG8y6/jTOvWNiTY9+7Yg0Ll6x8rGd174o1w280hg5/6daccNAOE3rMDdV4hNVQCrAH8HC/5cP9Bl1CM2x4LLAYeAS4BdhkFDVsBPyUpofV3/IBlklaDycctEPP3rAvunExR5yzgG1mbMqi5av5xGG7cOius3tSi9bPeITVXknS0bt6GfArmh5WgGeVUq4cZNu++08b9y1IsgWwI3Bc33ZJdluP2v+LZjhyWSllRRfbreusS1L9Dt11Ni/ffgvuXr6GrWdMY8vpU3tdkkZpPGYDPhs4I8kL2g/ingjMK6XcBpwHnJvk8CTPS7J7kvcneV277SKa3terkzwzyWbA/cAy4Jgk2yfZD/gHmt7VaJwHLAEuTrJfku2S7Jvk7/tmBA7iTmDPJNsmmdk3w1FS3bacPpWXbvMMg2qSG4833PNoeiDX0swE/CdgXrvuaJoZgZ+kmVl3CbAvTUhRSlkMnAJ8lCZQziqlPAocSTOz8CaaD+t+CHhoNMWVUla3x7wduKCt4yvAM2iCcTCfpuld3UIzE3A098skSaOQx8+F0BjywkpS9zLQQoeyJEnVm/RhleTkJKsGeVza6/okSetv0g8DJpkBzBhk9Zr2PlgvTO4LK0m9MeAw4KQPq4p5YSWpe96zkiRNToaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpeoaVJKl6hpUkqXqGlSSpesOGVZKrkpw1EcWsryT7JylJZnaxzalJbhrPuiRJ62fce1ajCZD18GNgK+A3Y7nTJOcmuWQs96neWbpyLTcsup+lK9f2uhRJIzSl1wWMlSRPKaWsA37d61pUr4tuXMyp829mmxmbsmj5ak495EUcuuvsXpclaRgppQzdILkKuBV4CPjTdvE/An9dSnk0ySbAacCbgRnAzcAHSymXJdkWuKPfLr9SSjkqyVzgA8DOQAGuA95bSvnlsEX/737fBBwD7A2cCNwEXAk8s5SyrG37Z8CpwEzgCuBS4OxSStr1pwKHAx8BPgps2bZ7eyllWbv+lH4lHFBKuWqYMoe+sEM48pwFo9100rnn/tUsXjFxPZypUzZi/rv3Yc6s6SxcspJDPncNax95dMKO32n25lN5zjM27cmxJ9L5x+7d6xI0uWSghSMdBnxz23Zv4Fjgz4H3tuu+DOxHExw7A18B5if5feBu4LC23Ytohuje0z5/GnAGsCewP/DbdrtNRlgTwMeBzwM7ARf1X5lkb5pgPRt4CfBvwN8NsJ9tgSOB1wIHA7vSBBfAp4FvAD9o69+KZrhRk9BWm09jzqzpAMyZNZ2tNp/W44okjcRIhwHvBY4vTTfs1iQ7AO9LcjHwRmDbUspdbduzkhwIHFtKOS7J8nb50r7eDkAp5cLOAyQ5GniAJryuGWFdnyulfLNjH9v3W3888P1Syunt89uS7EHTG+s0BTiqlPLbdj9fBI5u61yVZA3wUCllQoYYn0x/ic67/DbOvGLhhB3v3hVrWLhk5WM9q3tXrJmwY/d3+Eu35oSDdujZ8aXJZKRh9ZPy+PHCBTRDf/vQdNluSR7Xc3sq8O9D7TDJ89t97AU8k6bnthHw3BHWBHD9MOt3BOb3W3YtTwyrRX1B1foVzXCgxtkJB+0woW/YF924mCPOWfDYPatPHLaL96ykSWAsJlgUYA/g4X7Lh/uT9RLgHpphxcXAI8AtQDfDgA920XYo/Wsv+Bm0DdKhu87m5dtvwd3L17D1jGlsOX1qr0uSNAIjDau9kqSjd/Uymt7HApqe1bNKKVcOsu269ufGfQuSbEHT6zmub7sku3VRz0jdShOknfYcxX7W0VG/Jrctp081pKRJZqS9h2cDZyR5QZLDaWbezSul3AacB5yb5PAkz0uye5L3J3ldu+0imp7Kq5M8M8lmwP3AMuCYJNsn2Q/4B5re1Vj6LHBwkhOTzEnyNppJFN26E9i5Pf+ZSZ4yplVKkoY00rA6j6ZncS3wJeCfgHntuqNpZgR+kqYncwmwL01IUUpZTDP1+6PAEuCsUsqjNLPvdqGZbn428CGa6fFjppSygOb+1PHAz4FDgdOBbudKfwn4Jc09svuAV4xdlZKk4Qz7OasNTZJ5wIGllBeP86GeXBdWksbGgJ+z2mC+wWIwSU4ELgdWAQcC7wBO7mlRkqSuVDnjLcnJSVYN8ri0y93tDlxGM9z4HuAkmg8jS5ImiSqHAZPMoPnqpoGsae+D1a6+CytJ9RtwGLDKsNpAeGElqXvr9d2AkiT1jGElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqnmElSaqeYSVJqp5hJUmqXkopva5hg5Tke8DMHh1+JrCsR8eeCJ7f5Ob5TW7jfX7LSilz+y80rDZASa4vpeze6zrGi+c3uXl+k1uvzs9hQElS9QwrSVL1DKsN0xd7XcA48/wmN89vcuvJ+XnPSpJUPXtWkqTqGVaSpOoZVhuAJFclKf0eXx/BdocluSXJQ+3P105Evd1IMiPJ55LcmmRNkruTfCHJFsNsd9QA16QkmTpRtQ9R23FJ7kiyNskNSf5gmPb7te3WJrk9yTsmqtZuJDkpyXVJHkhyX5L5SXYeZpttB3mdnvA5m15LcuoAdf56mG1enOQ/2t/dxUn+NkkmquZuJLlzkNfiO0NsM1D7cfn9nDIeO1VPfBk4ueP5mqEaJ9kbOB84BfgW8DrggiSvKKVcO25Vdu/ZwGzgr4Bb2n9/HvgacPAw264Gnt+5oJSydhxqHLEkRwJnAscB17Q/L02yUynlrgHabwd8F/hn4C3APsDnk9xXSrlw4iofkf1pXpvrgAAfBn7QntvyYbadC/ys4/lw7Xvlv2nOs8/vBmuY5OnA5cDVwB7AjjT/Tx8E/n78Shy1PYCNO55vBdwAfGOY7Y4BLul4/tsxrqtRSvExyR/AVcBZXW5zPnB5v2U/AL7W6/MZQe2vAh4Fnj5Em6OAVb2udYC6rgW+1G/ZQuDjg7Q/HVjYb9k/Agt6fS4jONfNaN7MDxmizbZAAXbvdb0jOJ9TgZu6aP9O4AFgWseyDwKLaSe31fwAPgCs6Kx/gDYFOHwi6nEYcMPxhiTLktyc5NNJpg/Tfm/g+/2WXQa8fHzKG1NPBx6i6TkNZVqSRUnuSXJJkl0noLZBJdkEeClPvO7fZ/DrPtjrtHuSp4xthWNuOs2thvtH0PZbSZYm+VGSw8e5rvXxvCS/aodxv57keUO03Rv4YSmlc5TjMprRgm3Hs8j11Q5Vvg34ar/6B3Jm+95zXZJ3JBmXXDGsNgz/F3gzcABwGnAYMNwQ0bOAJf2WLWmXVyvJ5jTn+KVSyiNDNP1v4M+APwHeCKwFfpRkzrgXObiZNMMs3Vz3wV6nKfTuuydH6kzgp8CCIdqsAt4PvJ6mx3wFcH6St4x7dd27lqbHPpdm6OtZwI+HuH862GvXt65mBwHbAV8apt3fAkcCBwJfpxnePHnILUbJe1aVSvIRmm74UA4opVxVSun8kN4vktwOXJtkt1LKf41flaPXzfl1bLMZMJ9mGOWvhtqwlLKAjjfJJD+meeN8N3D8qIrWiCX5DM39tX1KKYPe1ymlLOPx92+uTzKT5vX96vhW2Z1SyqWdz5P8BLgdeCvwmZ4UNX6OAa4rpfxsqEallNM6nv40ycY0/68/MtYFGVb1OoPh/7M+4YZ863qaewVzgMHC6tfArH7LZrXLJ8IZdHF+bVB9t336mtLlRIlSyu+SXE9zTXplGc3r0s11H+x1eoRKv9k7yTzgDTR/bNw+il1cCxw9tlWNvVLKqiQ3M/jv1GCvXd+6KiXZkmZE4l2j2Pxa4OlJZpVS+vcq14thVan2L87Rvhm9mGa46d4h2iyg6ep/qmPZQcCPR3nMrnRzfu39t0tpZpjNLaWs6vZ47Rj8Ljx+xtmEKqWsS3IDzXW+oGPVQQw+bLsA6P+RgoOA60spD499lesnyZk0w0IHlFJuHeVuXsLQv7tVaD8GsSNw5SBNFgCnJ5na8cfVQcCvgDvHv8JRO4rmnvDXRrHtS2iG3FeMXTmtXs848bF+D5qp2X8L7E5z0/ZVwC9pelQbd7S7go4ZZzQ39B8B/obmP9xJwMPAXr0+p37nN53mP33fX7DP6nhsMsT5nQK8Enhe+x/on9vz27PH53MksA54O/BCmvs6q4Bt2vX/AvxLR/vtaKY6n9G2f3u7/WG9fm0GOLezaWa//WG/12mzjjYfB67oeP5W4E3tub2A5v7VOuCEXp/PAOf3aWC/9jXZi2a69gMdr13/c/s9mh7U14GdaT4e8gDwl70+lyHOMcBt9Jux2q77C+DWjueH0AwX7ty+D72dZtr6meNSW68vjo/1fAFha+A/gN/Q/DX0P+0b4Ix+7e4Ezu237HDg1vbN4ZfA63p9PgOc3/4002MHeuw/2PkB84BF7TVZSjMLa+9en09b23FtvQ/RfI5l3451VwFX9Wu/H80fHw8BdwDv6PU5DHJeg71Op3a0ORe4s+P5W2k+P/dg+0Z+PfCWXp/LIOf3dZpe0Tqa+6YXAjsNdm7tshfTfM5qLU1v8RQqnrZOM0mrMMAfdTRT90vH87nAjcDK9vX7BfAeYMp41OYX2UqSqufUdUlS9QwrSVL1DCtJUvUMK0lS9QwrSVL1DCtJUvUMK0lS9QwrSVL1DCtJUvX+PwBETL050mftAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x532.8 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_forest(\n",
    "    [\n",
    "        m1_trace\n",
    "    ],\n",
    "    var_names=['alpha','sigma', 'beta_left', 'beta_right'],\n",
    "    combined=True,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6413f51f-cbf9-4447-8597-2f2760217fa7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Joint Posterior distribution of beta_left and beta_right')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEXCAYAAABLZvh6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAp80lEQVR4nO3de7xUdb3/8dfH7UY3eNmal2IEMTVSQpmOqRzO+UXeUPEy+jPJ8pzsgnVOF1HDwEteAqEopUd1TqeytCCFDHcWmZfUOnEOGslWwkuIIrq9gOH2AjvZbD7nj7UGh+3MntuaWXN5Px+PebCZWfNdn5lZsz7zvazv19wdERFpTjvEHYCIiMRHSUBEpIkpCYiINDElARGRJqYkICLSxJQERESaWFMnATN7w8zeHXcclWRmK81sfIXKHmFmbmY7hv+/w8w+HlHZ/2xmT2T8f42ZHRdF2WF5FXtfBtinmdmPzewVM3swy+PnmdkfqxlTlhi2+0wL2H6Gmb1sZi9WOrYs+875fhX7OqqpmGMv6uM+m4ZMAmZ2v5l9Ot927r6Luz9VYJluZgcN8Ph5ZtYXJpbXzKzTzE4pJu5i91kIdx/l7veXU0YR+zrJ3W/Kt10hr8vd/9vdR0YRl5ndaGYz+pVftfclwz8BxwP7ufuRURYcRwIxs+HAxcCh7v7OWkhiUah0Aonq2DOz8Wb2XLnlNGQSiNH/uvsuQDtwA7DQzPaII5ByD+A4f0HV4q+3iOwPrHH3jXEHEpHhwN/cfV3cgdSDmj2u3b3hbsD9wKfDvycDTwIbgNuBoRnbOXBQ+PeNwHeBxcDrwAPAgeFjfwi33Qi8AUzKss/zgD9m/H9I+JwjgN2BnwDrgWeAy4Edwu0OAn4PvAq8DCwYaJ/AKUAn0A38D3BYxj7XAF8GHgHeBHYM7zsufHwnYC7wfHibC+wUPjYeeC58/ovAT7O8xhbgG2GcTwGfC2PcMcv7XvDryrbv9H39Xtt04FHgFeDHwM7Z3vvMzxY4H+gFNof7+1VGecW8LxcD64AXgE8McOwNJTjONhAcd5PD+z8F/B3oC+O4OscxtAT4Tvi+PQ4cm/H47gQ/Ll4AuoAZ4WdySL+yu8PtJwLLgdeAZ4GrCvjujOj3meba53FAD7A13OeCbDFkKf8TwGME37GngM9kPDbgew28I3xvXwMeBL7a/3PP8jrODz/TF4AvZTy+AzANWA38DVgI7Bk+tjZ87hvhbSxwIHBvuO3LwHygvYD3cw0DfyfbgJsIjunHgEt4+3H/pfD5r4bv884E55fM9/8NMs5tRZ0v4zpRV/JGeDICjgk/sPcTfNG/Dfyh/4ki/PvG8AM+Mvyg5gO3ZNs2xz7PSx+Q4fMvCA/0dAL4JbBreHD+FfhUuO3NwGXhQbkz8E+59gkkwy/HUQRfxI+HB0n6hLWGIEEMA9oy7ksfcNcAS4F9gL0JkshXM76AW4Cvhe9VW5bX+FmCE9MwYE/gPnIngWJe19v2TfYk8JeMfS8BZvR/7wf4bGdk+XIW875cA7QCJwObgD1yHAd/AP4jfM1jCBL/MbnizHIMbQEuDPc1ieCLnz453Qb8F8EJYB+CE+FnBngPxgOjw8/gMOAlIJXnuzOi32c60D77f0YDvr5wm4kEJ1QDPhi+l+8v5L0GbiE4WQ8B3keQlPIlgZvD7UeHn0X6M78g/Mz3Izjm/gu4Odt7EN53EEFT3k7hMfIHYG4B56I1DPydnE3wY2mPMJZHePtx/yDBj4s9CRLFZ7O9/yWfL8stoBZvvJUEbgC+nnH/LgS/CkeE/+9/ovhhxrYnA49n/L+QJLCF4Bf6y+EBdhzByXozQbtpetvPAPeHf/8E+D5BO3H/MvufLP+T8OSUcd8TwAczDphPZjkI0wfcauDkjMcmEDRPpA+ozYS/rnO8xnvTB2D4/xPInQSKeV1v23f/Azx8HZn7PhlYnfHel5ME8r0vPWx/QlgHHJ3ldQ0j+CW8a8Z9s4Abc8WZ5Rh6HrCM+x4E/gXYl+CXZFvGY+cA9xVSdrjNXOD6PNuMSH+mBeyz/2eUN4Ys++sALsj3XhN8j3qB92Y8dm2u/WW8jsztvw7cEP79GNvXst4Vlr8jWZJAlvJTwPICXt8aBv5OPgVMyHjs07z9uD+332v4Xrb3v9Rbo/cJDCVofgHA3d8g+LWfyLF95giHTQRJoxhL3b3d3fdy96Pd/R5gL4JfNc9kbPdMRgyXEPwqejAcNfDJAcrfH7jYzLrTN4ITz9CMbZ4d4PnbvR/h35nPXe/uf8/z/Mzyn8m1IcW9rkL2TZZ9D821YZHyvS9/c/ctGf/PdWwMBTa4++v9ysp1vGXT5eE3vF8s+xMcRy9kfPb/RfDrPCszO8rM7jOz9Wb2KkFNbq8iYil6n/mY2UlmttTMNoTlndwvplzv9d4EJ+hCj7+0XMfM/sBtGa/rMYIEvm+OuPc1s1vMrMvMXgPmUfh7me87mfl4tm3LPS8NqNGTwPMEHzYAZjaEoF2xq4oxvEzwC2P/jPuGp2Nw9xfdfbK7DyWoIfzHACNnngVmhokmfRvs7jdnbOM5ngv93o8wjucLfC4E7arD+j0/qyJfVyH7Jsu+07FvBAanHzCzdxZZdr73pVDPA3ua2a79yirmeEuYmWWJ5VmCX+V7ZXz2u7n7qHC7bK/xZwRt6MPcfXfgewSJuVD59tnfgO+zme0E/IKgX2lfd28HflNgTOsJatoFHX8Zch0zzwIn9fsu7ezuXTlex7Xh/aPdfTfg3ALjJkd5aS8QNANli7eccgvW6EngZuATZjYmPACvBR5w9zUllPUSUPQ1Be7eR9COOdPMdjWz/YGLCH5JYGYfNrP0QfAKwQe7Ncc+fwB8NvyFZ2Y2xMwm9jvpDORm4HIz29vM9gK+ko6jQAuBL5rZfuGop2m5NizydRXqc+G+9yTob1gQ3v8wMCr8nHcGrur3vHz7K/d9AcDdnyXoT5hlZjub2WEEHcLFlLUPwXvcamYfJuj0/Y27vwDcBXzTzHYzsx3M7EAz+2D4vJeA/cxsUEZZuxLUTP5uZkcCHy3y9eTbZ3/ZYsg0iKBNfT2wxcxOImhSLCSWPmARcJWZDTazQwn6xPK5Itx+FEGndPqY+R7Bd3J/gPCzPz18bD3BsZp5zOxK0Pn6qpklgKmFxF2AhcB0M9sjLPfzRTz3JeAdZrZ7OQE0chLwsDnmCoJfHy8QdEh9pMTyrgJuCquPZxf53C8Q/Fp9CvgjwS+0H4WPfQB4wMzeIPjVdoG/de3Cdvt092UEo52+Q3BifZKgHbZQM4BlBJ1PK4CHwvsK9QPgToKT7kMEX8pcCn5dRez/ZwQnpacI2vFnALj7Xwk6E+8BVhG8x5luAA4N99eRpdxy35dM5xC0KT9P0Kl6ZXgcFuoB4GCCGuRM4Cx3/1v42L8SnEjTI6RuJWjLhqC/ZiXwopm9HN7378A1ZvY6QWJbWMLrGWif/WWLYZuwmeyLYRyvECSl24uI5fMETSEvEvTz/LiA5/ye4HvyO+Ab7n5XeP+3wn3fFb4/SwkGXODumwje+yXhMXM0cDXBAJNXCUYQDnTsF+MaghFRTxMcv7cS1L7ycvfHCX7APBXGWVLzqG3f/NgYzOwh4Bp374g7FhGRQpnZvwEfcfdcta3INVxNIKz2HUIwPlpEpGaZ2bvMbFzY1DaS4BqJ26oZQ0MlATP7GkFzwZfdvZCRAyJNx8w+ZsH0Jv1vK+OOrd6Y2fAc7+UbFkyrkc8gghFXrxM0p/2S4DqTqmnI5iARESlMQ9UERESkOLU5odEA9tprLx8xYkTcYYiI1JU///nPL7v73v3vr7skMGLECJYtWxZ3GCIidcXMsvaTqjlIRKSJKQmIiDQxJQERkSamJCAi0sSUBEREmljdjQ4qRcfyLubc+QTPd/cwtL2NqRNGkkoWM8W7iEhjavgk0LG8i+mLVtDT2wdAV3cP0xetAFAiEJGm1/BJYM6dT2xLAGk9vX3MufOJ7ZLA5R0ruPmBZ+lzp8WMc44axozU6GqHKyJSVQ2fBJ7v7sl6f1d3D8lr7uKVTb1ve6zPnXlL1wIoEYhIQ2v4juGh7W05H8uWADL97IG1UYcjIlJTGj4JTJ0wkrbWlu3uK3Rh0K0O42bfS8fyai5JLCJSPQ2fBFLJBLPOHE2ivQ0DEu1tRa3OnO5IViIQkUbU8H0CECSCzE7gcbPvpStHX0E22TqSRUQaQcPXBLKZOmEkrS2FNgoFcnUwi4jUs6ZMAqlkgjlnHc4eg1sLfo6j/gERaTx1t7zkEUcc4VGvJ9D/grKBtO5g7LLzjnRv6tXVxyJSN8zsz+5+RP/7m6JPIJ/0STw9tUT74Fbcobvn7UNIe7f6tqGluvpYROqdkkCof+cxwAHTFucdSaROYxGpZ03ZJ1CogS40y9TV3cOIaYsZc/Vd6jMQkbqiJDCAbBeaDaS7p5epP39YiUBE6oY6hvPInIa62HdqyKAWZp4xWk1FIhK7mu4YNrMWYBnQ5e6nxB1Ppsy+gmIvMtu4uY+Lf/7wtnJERGpNrTQHXQA8FncQ+RTbPATQt9WZsqBT1xiISE2KPQmY2X7AROCHcceST3oeomIuMkvr6u7hwgWdjJi2WAlBRGpG7EkAmAtcAmzNtYGZnW9my8xs2fr166sWWDapZILlXzmBuZPG0N5WXDJI9yloUjoRqRWxdgyb2SnAye7+72Y2HvhSvj6BancMF+ryjhXbFqIpVKK9jSXTjqlQRCIib8nVMRx3TWAccJqZrQFuAY4xs3nxhlSaGanRDBlUXH9BV3ePmoZEJFaxJgF3n+7u+7n7COAjwL3ufm6cMZVj5hnFL0WppiERiVPcNYGGkkomaGst/i3t6e3j4oUPc4A6jUWkymomCbj7/bV2jUApZp15GK07FLdWAQSL2zuqGYhIddVMEmgUqWSCOR8+fNtylu1trQwusnaQnpRORKTSNG1ElVzesYL5S9cWNfWEgdYsEJFI1OrooKYxIzWa6yeNIRHOTNpi+ZuM0s1DFy3oVPOQiFRETcwd1CxKXfB+KzBlQee2MkREoqKaQIyKnYtoyoJODpi2mMs7VlQwKhFpJqoJxChzWctCZyd1YN7StduuTt5jcCtXnjpKNQQRKYlqAjFLJRMsmXbMtr6CYr2yqZept2ohGxEpjZJAjZg6YWRJ1xcA9PY5ly56JOKIRKQZKAnUiPT1BcXOTJq2qXerpqkWkaLpOoEaVcqspJkO3mcId180PrqARKSu6TqBOjMjNZq5k8bQ2lJaE9GqdRs5aubdEUclIo1GSaCGpZIJ5px1eMmdxi+9vpmP/eB/I45KRBqJmoPqyMd+8L8sWb2hpOfOnTRGw0hFmpiagxrA/MljGXfgniU9d8qCTkZ95bfqNBaR7agmUIc6lncx9eed9OZclTm/hCamE2kqqgk0kFQywaprJ5ZcK4BgYrqpP9dFZiLNTkmgjs2fPJZzjx5e8vN7tzoXaoZSkaam5qAGcegVd7CpnPYh1EQk0sjUHNTgrj3zsLLL0NKWIs0n1iRgZsPM7D4ze9TMVprZBXHGU89SyURZTUNpWtpSpLnEXRPYAlzs7ocCRwOfM7NDY46pbqWvMt5jcGnzD6V1dfcwYtpiktfcpVqBSIOLdT0Bd38BeCH8+3UzewxIAI/GGVc9y1y9rGN5Fxcu7KTUbp9XNvVy0cLObeWKSOOJuyawjZmNAJLAA1keO9/MlpnZsvXr11c9tnqVSia4/uwxZZWx1YMLzVQjEGlMNTE6yMx2AX4PzHT3RQNtq9FBxetY3sX0RY/QU+boIdBKZiL1KtfooNiTgJm1Ar8G7nT36/JtryRQnnLmH8p07tHDmZEaHUFEIlINNTlE1MwMuAF4rJAEIOWbP3kscyeNobQJqt8yb+ladRyLNIC4+wTGAf8CHGNmneHt5JhjanipZIKnZ5c37QQEHcdTFnRy/HX3RxOYiFRdrEnA3f/o7ubuh7n7mPD2mzhjaiblTjuRtmrdRiUCkToVd01AYpa+tqBcq9Zt5JAr7lDzkEidURIQUskEa2ZPZO6kMZS4miUAPb1bmbKgU6uZidQRJQHZJpVMsHrWRA7eZ0hZ5SxZvUHNQyJ1QklA3ubui8aXnQhWrdvIiGmLGTFtMZd3rIgoMhGJmpKAZHX3ReMj6SuAYDipmohEapOSgOSU7ivYbaeWsstasnoDB1+6WB3HIjVGSUDyeuTqEzn36OFYmVeY9W4N5iE6aubd0QQmImWLfdqIYmnaiPiNmLa47DJ2bjEen6nrAkWqpSanjZD6FEVfwd/7nIOml59MRKQ8sa4nIPUpPYPopYseKWtd4y3+Vq1CE9KJxEPNQVK2qGYmBS12L1Ipag6Sipk/eWzZk9GldXX3MGVBp64tEKkSJQGJRFST0aXNW7pWw0lFqkBJQCIzIzWaNbMnRpYMtKylSOUpCUjk0jOTDm4t//BSIhCpLCUBqYhUMsGjXz0pkuGkWrhGpHKUBKSiUskEifa2sstJT0inOYhEoqUkIBU3dcJIWstZqCDDktUbNDOpSISUBKTiUskEc846nB13iCYRgGYmFYmKkoBURSqZ4MlrT47segJ4q1ag/gKR0sWeBMzsRDN7wsyeNLNpcccjlTV/8ljWzJ4Y6VDSVes2ah4ikRLFmgTMrAX4LnAScChwjpkdGmdMUj3poaQRjCTdNg+RagUixYm7JnAk8KS7P+Xum4FbgNNjjkmqKJVMsOraYJH7ctcrgLdGEY36ym91fYFIAeJOAgng2Yz/Pxfetx0zO9/MlpnZsvXr11ctOKmeVDLB9WePoa21/FXMADZu7mPKgk51HovkEXcSKIi7f9/dj3D3I/bee++4w5EKSSUTzDpzdCTXFaQtWb1Bw0lFBhD3egJdwLCM/+8X3idNKpVMbDeN9AHTFlPuZOfzlq7VWgUiORRcEzCznxZyX5H+BBxsZgeY2SDgI8DtZZYpDeTp2RMjmXpCVxuLZFdMc9CozP+EI3v+oZydu/sW4PPAncBjwEJ3X1lOmdJ4UskEa2ZPLLvtUtcViLxd3u+VmU03s9eBw8zstfD2OrAO+GW5Abj7b9z9Pe5+oLvPLLc8aVzXTRpDawRXHadHEKmvQKSI5SXNbJa7T69wPHlpecnm1rG8izl3PkFXd09kZWpJS2kGuZaXLGqNYTNLAPuT0aHs7n+IJMICKQlIWpRrG4MWu5fGlisJFDw6yMxmE3TcPgr0hXc7UNUkIJI2f/JYIOj0jcK8pWs5Yv89VSOQplJMX9sZwEh3P9ndTw1vp1UqMJFCRTF6KG3Kgk5dbSxNpZgk8BTQWqlAREqVSiaYO2kMbVFMQoSuNpbmkrdPwMy+TdDskwAOB34HvJl+3N2/WMkA+1OfgAykY3kXX/7FI7y5ZWsk5bUYfPPsMWoikrpXTp9A+oz7Z3Qhl9S4zCuOj7/uflat21hWeX0eNBGlyxZpNEWNDqoFqglIMaJIBGkaSir1LFdNoJhpI1aY2SP9bv9tZteb2TuiDVckGndfNJ5zjx5OFAtbdnX3qK9AGk4xPWl3AIuBj4W3XxE0Fb0I3Bh5ZCIRmZEazdMRrmS2ZPUGTT0hDaOYK4Yfcvf3Z7vPzFa4e1WuslFzkJSjY3kXU3/eSW8E/cb77jqIBy47vvyCRKqg7IvFgBYzO9LdHwwL/ACQXgFkSwQxilRcuuO4Y3nXtg7fUr30+uZtF6qpv0DqVTHNQZ8GbjCzp81sDXADMNnMhgCzKhGcSKWkZyYdd+CekZSX7i/QRWZSb4oeHWRmuwO4+6sViSgPNQdJ1KKoFWRaM3tiZGWJRKXk5iAzO9fd55nZRf3uB8Ddr4ssSpEYpJtwokoE6SYiTUgn9aCQ5qAh4b+75riJ1L1081BUI4ggmJBOw0ml1hXUHBSuIvZFd7++8iENTM1BUg1RNxG1mHHOUcNUM5DYlHWxmLv3AedEHpVIjUpPSheVPnfmLV2r1cyk5hQzOmiJmX3HzP7ZzN6fvlUsMpGYpZKJSJuHIGgiEqklxVwsdl+Wu93djylpx2ZzgFOBzcBq4BPu3p3veWoOkmqrxJKWcydpZlKprkiWl8yzg4+7+01FbH8CcK+7bzGzrwG4+5fzPU9JQOJ0eceKyH7ND2oxvn7W4UoGUhVlTyBXgAuK2djd73L39JXGS4H9IoxFpCJmpEazZvZEhgxqyb9xHpv7nCkLOtVPILGKMgmUM1HjJwkmqMtesNn5ZrbMzJatX7++jN2IRGPmGaNp2SGKuUmDfoIR0xZrOKnEIsok8LZ2JTO7x8z+kuV2esY2lxHMPTQ/Z8Hu33f3I9z9iL333jvCkEVKk0om+OaHD4+kRpC2ZPUGjpp5d2TliRQiyj6B5e6eLPI55wGfAY51902FPEd9AlKLoly8BtRxLNGrRp/AkmI2NrMTgUuA0wpNACK1Kr14TVSmLOjUmgVSFUXVBMxsIjAK2Dl9n7tfU9KOzZ4EdgL+Ft611N0/m+95qglIrYvyauMdDZ6cpQnppHxlrydgZt8DBgMfAn4InAU8WGpA7n5Qqc8VqWVRTki3xYMJ6TQzqVRKMc1B/+ju/wq84u5XA2OB91QmLJH6lp52YseIRhBp9JBUSjFJIH255CYzGwr0Au+KPiSRxpBKJnjy2pMjW7hmyeoNSgYSuWKSwK/NrB2YAzwErAFurkBMIg1l/uSxzJ00hj0Gt0ZS3pLVGzho+uJIyhIpZu6gndz9zfTfBJ3Df0/fVy3qGJZ6FvUU1e1trVx12igNJ5W8ohgiuq0O6u5vhstLql4qUoR0X0GLRdNX0N3Ty5QFnYy5+i6tbywlKWR5yXcCCaDNzJK8NT3EbgSjhUSkCOlf7dMXraCnty+SMrt7epm+aMV25YsUopAhohOA8wgmeMtcT/g14NIKxCTS8NIn6iinqO7p7ePihZ1KAlKUYvoE/r+7/6LC8eSlPgFpVFH2F+y76yAeuOz4SMqSxhBFn8ASM7vBzO4ICzzUzD4VWYQiTS7Klcxeen0zI6Yt1jTVklcxSeDHwJ3A0PD/fwWmRB2QSDObkRrN3EljaG+LZjjpvKVrOeSKO9RpLDkVkwT2cveFwFaAcEGYaHq1RGSbVDJB55UnhKOIyi+vp3crUxZ06iIzyaqYJLDRzN5BuG6AmR0NvFqRqESEVDLB6lkTI2siSl9xrNlJJVMxSeAi4Hbg3Wa2BPgJ8IWKRCUi26SXtIxq+olV6zYqEcg2xSSBR4HbgD8BLwE/IOgXEJEqmD95LGtmT2THCJqIVq3bqOYhAYpLAj8B3gtcC3ybYAbRn1YiKBHJ7cmImojSzUPjZt+rjuMmVsx1Ao+6+6H57qs0XScg8pajZt7NS69vjqSsg/cZwt0XjY+kLKk9UVwn8FDYGZwu8ChAZ2ORGD1w2fGce/RwopiJaNW6jZqdtAnlrQmY2QqCEUGtwEhgbfj//YHHVRMQqQ1RXXGsJS0bUznLS55SgXhEJGJRLWu5xeH46+5X01CTyJsE3P2ZagQiIuVLJ4KLFnQGV3WWaNW6jUGn8YF7Mn/y2GiCk5pUTJ9ARZjZxWbmZrZX3LGINIJUMsFTs6MbQXTYlb+NICqpVbEmATMbBpxA0M8gIhGakRodSSJ47c0+DSVtYHHXBK4HLiGcikJEopWekK6ttfyveld3D1MWdHLg9MVKBg0ktiRgZqcDXe7+cAHbnm9my8xs2fr166sQnUjjSCUTPPbVk5g7aUwkQ0n7POh81jTVjaHgi8VKKtzsHuCdWR66jGBVshPc/VUzWwMc4e4v5ytTQ0RFynN5xwrmLY2mBXbIoBZmnjFaq5nVgXKGiJbM3Y/LEcxo4ADgYQsW3N6P4GK0I939xUrGJNLsZqRGA0SSCDZu7mPqrUFlXomgPlW0JlBwEKoJiFRdx/KusoeSZtrB4KNHDd+WZKS2RDFthIg0kPRQ0rmTxjA4go7jrR7ULjQ7aX2piSTg7iMKqQWISPRSyQSPhh3HUaxktmT1BnUa15GaSAIiEr/0SmZzJ40pu6x5S9dqKGmdUBIQke2kkgnWzJ7IwfsMKascDSWtDzXRMVwMdQyLVNfx193PqnUbyyoj3co0tL2NqRNGaiRRDNQxLCIlufui8WU3EXl46+ruYeqtD6uZqIYoCYhIXqlkIpJ5iAB6+5zLblMTUa1QEhCRgqTnIUq0t5Vd1sbNfRpKWiPUJyAiJTlw+mL6Ijh97DG4lStPHaV+ggpTn4CIROqbZ49hhwiuK3hlUy8XLexUP0FMVBMQkZJ1LO9izp1P0NXdU3ZZGkFUWaoJiEjkUskES6YdE0mnsUYQxUNJQETKll7FrMWiWLEgGEF09a9WRlKWDExJQEQiMSM1mtWzTo5s8ZpXNvXqauMqUBIQkUilkgk+FtE1BfOWrmXEtMVKBhWkjmERqYiO5V1c/auVvLKpN9JyE+o4LkksK4uJSPNKJRPbTtQdy7u4+OcP07e1/B+dXd09TF+0Yts+pDxqDhKRikslE3zzw4ezx+DWSMrr6e1Tx3FElAREpCpSyQTLv3ICa2ZPjGRI6SubejWMNAJKAiJSdel5iMpd1nLKgs5oAmpiSgIiEov0spbl1gpGTFvMuNn3qlZQoliTgJl9wcweN7OVZvb1OGMRkXikawXlXGjW1d3DlAWdHHLFHUoGRYptdJCZfQg4HTjc3d80s33iikVE4pUe5TP11ofpLWNq0p7erduaiDRyqDBxDhH9N2C2u78J4O7rYoxFRGKWPmlHcW2BEkHh4mwOeg/wz2b2gJn93sw+kGtDMzvfzJaZ2bL169dXMUQRqab+I4jKmX5C01MXpqJXDJvZPcA7szx0GTATuA/4IvABYAHwbs8TkK4YFmkeHcu7uOr2lXT3lFczGDKohZlnjG7qmkGuK4ZjmzbCzH4LfM3d7wv/vxo42t0H/KmvJCDSfKKaguLco4czIzU6oqjqSy2uJ9ABfAjAzN4DDAJejjEeEalRmc1E4w7cs+Ry5i1dq7WN+4kzCfwIeLeZ/QW4Bfh4vqYgEZH5k8eWdW3BktUbNCtpBs0iKiJ1qWN5Fxcu6KTUM1jrDjDnw2Oapp+gFpuDRERKlkomuH7SmJJPYr1bg6GkB0xv7vUKlAREpG6lkgmumzSGtjLmIHJ/a/GaZhxSqiQgInUtlUzw2FdPimSa6ikLOhlz9V1NlQyUBESkIVx56ijaWlvKLqe7p5eptz7cNIlASUBEGkIqmWDWmaNJtLdhUFbNoLfPuXTRI9EFV8O0vKSINIzMJS3TLu9Ywbyla4sua1PvVjqWdzX86CHVBESkoc1IjWbN7InMnTSm6OdedlvjjxpSTUBEmkIqmSh66omNm/s4YNribdcitLe1ctVpoxqqdqCagIg0jStPHUVrS3Fzk2ZejNbd09twI4hUExCRppH+BT/nzid4vruHlh2MLVuLv+a4u6eX6YtWbFdmvVISEJGm0r/zuNSO457ePubc+UTdJwE1B4lIU0uvcTxkUPHXGDzf3VOBiKpLSUBEml4qmWDlNScWPYJo5zKmq6gVag4SEQmlkgmWPbOh4Oahnt6tjJi2GCPoQE60tzF1wsi6aiKq/zQmIhKhdPNQor2t4Oeku5a7unu4cEFnXc1KqpqAiEg/mZ3Hh15xB5t6txb8XIdtNYl6WMpSNQERkQFce+ZhJT1v/tK1dXEtgVYWExHJo9RhpGktZpxz1LBYawZaWUxEpETpfoL2ttJmJu1zZ97StRwwrfZWMYutT8DMxgDfA3YGtgD/7u4PxhWPiMhAMvsJLu9Ywc8eWEuxFxvXYn9BbM1BZnYXcL2732FmJwOXuPv4fM9Tc5CI1JrLO1Ywf+naghe938HgqVkTKxpTf7XYHOTAbuHfuwPPxxiLiEjJZqRGc304rLSQ6em2OjXTaRxnTeAQ4E7ACJLRP7r7M/mep5qAiNS6QjqSE+1tLJl2TJUiyl0TqGgSMLN7gHdmeegy4Fjg9+7+CzM7Gzjf3Y/LUc75wPkAw4cP/4dnnsmbK0REYlXMiKI9Brdy5amVXacgliQwEDN7FWh3dzczA151993yPU81ARGpJ4dccQc9BVxs1tpizDnr8IolglrsE3ge+GD49zHAqhhjERGpiFlnHkZba/4ZSnv7nIsXPlz1voI4p42YDHzLzHYE/k7Y3CMi0kj6L2QzUNtLn3vVF6vRFcMiIlU0bva9dOVZh6ASnca12BwkItJ0pk4YmXebai5WoyQgIlJFqWSCPQYPPP3E0IxprDuWdzFu9r0cMG0x42bfG3mfgZKAiEiVXXnqqJydxW2tLdtqCx3Lu5i+aAVdYV9CV3cP0xetiDQRKAmIiFRZKplg1pmjty1c02LBdcaJ9jZmnTl6u87knt6+7Z6bXuA+KlpURkQkBpkT0uWSqwM5X8dyMVQTEBGpUekaQqH3l0JJQESkRvXlGMKf6/5SKAmIiNSoXIvd57q/FEoCIiI1auqEkW8bRZQ5eigK6hgWEalR/aecGNrextQJIyOdUkJJQESkhhUyiqgcag4SEWliSgIiIk1MSUBEpIkpCYiINDElARGRJlZ3i8qY2Xqg1lea3wt4Oe4giqB4K0vxVk49xQrxxru/u+/d/866SwL1wMyWZVvBp1Yp3spSvJVTT7FCbcar5iARkSamJCAi0sSUBCrj+3EHUCTFW1mKt3LqKVaowXjVJyAi0sRUExARaWJKAiIiTUxJoELMbI6ZPW5mj5jZbWbWHndM2ZjZiWb2hJk9aWbT4o5nIGY2zMzuM7NHzWylmV0Qd0z5mFmLmS03s1/HHUs+ZtZuZreGx+1jZjY27pgGYmYXhsfBX8zsZjPbOe6YMpnZj8xsnZn9JeO+Pc3sbjNbFf67R5wxgpJAJd0NvM/dDwP+CkyPOZ63MbMW4LvAScChwDlmdmi8UQ1oC3Cxux8KHA18rsbjBbgAeCzuIAr0LeC37v5e4HBqOG4zSwBfBI5w9/cBLcBH4o3qbW4ETux33zTgd+5+MPC78P+xUhKoEHe/y923hP9dCuwXZzw5HAk86e5Puftm4Bbg9JhjysndX3D3h8K/Xyc4SVVuovUymdl+wETgh3HHko+Z7Q78P+AGAHff7O7dsQaV345Am5ntCAwGno85nu24+x+ADf3uPh24Kfz7JiBVzZiyURKojk8Cd8QdRBYJ4NmM/z9HDZ9UM5nZCCAJPBBzKAOZC1wCbI05jkIcAKwHfhw2X/3QzIbEHVQu7t4FfANYC7wAvOrud8UbVUH2dfcXwr9fBPaNMxhQEiiLmd0Ttkf2v52esc1lBM0Y8+OLtLGY2S7AL4Ap7v5a3PFkY2anAOvc/c9xx1KgHYH3A//p7klgIzXQVJFL2JZ+OkHyGgoMMbNz442qOB6Mz499jL6WlyyDux830ONmdh5wCnCs1+YFGV3AsIz/7xfeV7PMrJUgAcx390VxxzOAccBpZnYysDOwm5nNc/daPVE9Bzzn7uma1a3UcBIAjgOedvf1AGa2CPhHYF6sUeX3kpm9y91fMLN3AeviDkg1gQoxsxMJmgJOc/dNcceTw5+Ag83sADMbRNCxdnvMMeVkZkbQZv2Yu18XdzwDcffp7r6fu48geF/vreEEgLu/CDxrZiPDu44FHo0xpHzWAkeb2eDwuDiWGu7IznA78PHw748Dv4wxFkA1gUr6DrATcHdwjLLU3T8bb0jbc/ctZvZ54E6C0RU/cveVMYc1kHHAvwArzKwzvO9Sd/9NfCE1lC8A88MfBE8Bn4g5npzc/QEzuxV4iKC5dTk1NiWDmd0MjAf2MrPngCuB2cBCM/sUwZT4Z8cXYUDTRoiINDE1B4mINDElARGRJqYkICLSxJQERESamJKAiEgTUxIQEWliSgLSlMxsROYUvwVsf56ZDS1xXzea2Vl5tnmvmXWG8/YcaGYfLWVfIsVSEhApzHkEc9RUSgq4NZy3ZxigJCBVoSuGpZntaGbzCSZOWwn8K3AIcB2wC/Aywcl/HHAEwdW0PcBYYCpwKtAG/A/wmULmhzKzf8hSfhKYAvSZ2bFhmYeEV0Xf5O7XR/JqRbLQFcPSlMKpqJ8G/sndl5jZjwjmnjkDON3d15vZJGCCu3/SzO4HvuTuy8Ln7+nuG8K/fwosdPdf5djXjcCvCeaJ+X2O8q8C3nD3b5jZ+HBfp1To5Ytso5qANLNn3X1J+Pc84FLgfbw131MLwVz12XzIzC4hWMxkT4KaRNYkkGFkEeWLVIWSgDSz/tXg14GV7j7g2rrhWrb/QbC04bPhr/hC1re1QsoXqSZ1DEszG56xmPpHCZYB3Tt9n5m1mtmo8PHXgV3Dv9Mn/JfDBW4GHPmT4YkBys+UuS+RilISkGb2BMFi9Y8BewDfJjihf83MHgY6CRYqgWDR8O+FnbVvAj8A/kIwDfefCtlZuI5zrvIzPULQSfywmV1YygsTKZQ6hkVEmphqAiIiTUwdwyIRMbPvElxTkOlb7v7jOOIRKYSag0REmpiag0REmpiSgIhIE1MSEBFpYkoCIiJN7P8AqjaW0Aovqe4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(m1_trace.posterior['beta_left'].values,m1_trace.posterior['beta_right'].values)\n",
    "plt.xlabel('beta_left')\n",
    "plt.ylabel('beta_right')\n",
    "plt.title('Joint Posterior distribution of beta_left and beta_right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d4cab9f0-87a3-4a31-9f8e-adc48ae63d01",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Sum of beta_left and beta_right')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEXCAYAAABPkyhHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcUUlEQVR4nO3deZhddZ3n8ffHJAIqi5J6BEIg3YjsqyWCyzRKM4IiTCsOaCtLw6TtBxroxraB5xmMjNOt0+LC0jI0S8IyCoJL2MZmBBoEWQoMSwj4RJBOaJQiQBIWg4HP/HF+BTe3blXdC3XuTVKf1/Pcp+4553fP+dYvlfu553fOPUe2iYiIie1NvS4gIiJ6L2EQEREJg4iISBhERAQJg4iIIGEQEREkDKLLJG0jaZ6k5ZKOa7H8JklH96K2hhpmSbqkzbbvlHRz+X1Or7u2Ftsfsb86+T26SdKHJD3cZtu9JS2uu6ZIGKy1JH1Q0m2Slkp6WtKtkt7b67qALwE32l7f9hnjueIeBclM4ClgA9snrg5hNh7qDBLbt9jeZjzWJWm2pK+Ox7omuoTBWkjSBsDVwJnAO4BpwFeAFb2sq9gSmN/rIsbRlsCDzrc32yJpcq9riNYSBmundwPY/p7tl22/aPtfbd8Hwz/1SZohyUP/Ucun26+WPYvnJF0laWNJl0paJukuSTNG2rikAyXNl/RsWdd2Zf4NwIeBs8p63z3CKraSdGfZ1k8kvaNh3XuWup6VdK+kvcv8/wl8qGHdZ5X535G0qKzrbkkf6rQzR9nmbOBw4Etlm7e2qqHF+n4g6bdlr+1mSTs0LJst6WxJ15ShpzskbdWwfF9JD5XXngVojPLXlXRZWdc9knZpWNdmkq6UNCjp0aFhO0n7AacAh5Tf494y/0hJC8q6HpH0l2303d6SFkv6e0m/BS5sHvqRtLukX5b1/qDU+9Wm9Zwo6UlJT0g6ssybCfx5Q/9fNVY9MQrbeaxlD2ADYAkwB9gfeHvT8lnAJQ3TMwADk8v0TcBCYCtgQ+BB4FfAnwKTgYuAC0fY9ruB54F9gSlUw0ILgTc3rPvoUWq/CXgc2BF4K3DlUK1UezhLgI9RfZDZt0z3jbRu4HPAxqXuE4HfAuuO0X+zOtjmbOCrTfWP+PuVNn8BrA+sA3wbmNewbHZZ/x6l5kuB75dlU4HlwMGlb/8GWDnS9srv8YeG9l8EHi3P3wTcDZwKvBn4Y+AR4KOt/kbKvI+XvwkBfwK8AOw+xu+6d6nx6+X3Xa/MW1yWvxl4DDi+1PVJ4KWhPm14/Wll+cfKdt/eqv/zeP2P7BmshWwvAz5I9Qb/L8CgpLmS3tnBai60/WvbS4HrgF/b/n+2VwI/AHYb4XWHANfYvt72H4BvUL0BvL+DbV9s+wHbzwP/HfivkiZRvbFfa/ta26/Yvh4YoHqDaMn2JbaX2F5p+3SqN6ROxqs73uZYbF9ge7ntFVRvurtI2rChyY9s31n6+lJg1zL/Y8B821eUvv02VbiN5u6G9t8E1gX2BN5LFWin2X7J9iNUfyuHjlL3NeVvwrb/DfhXqj2hsbwCfNn2CtsvNi3bkyr0zrD9B9s/BO5savMH4LSy/FrgOTr7N4w2JAzWUrYX2D7C9uZUn7I3o3rzaNfvGp6/2GL6bSO8bjOqT3pDdbwCLKL6hN2uRQ3PH6P6RDiVanz+02W45llJz1KF3qYjrUjSF8vQxtLSfsOyrnZ1vM3RSJok6WuSfi1pGfCbsqixpsY3+Bd4ra83o6FvbJtV+6qVxvavAIvLerYENmv6vU4BRvzAIGl/SberOiHhWapwaqcvB23/foRlmwGPl99lWM3FkhKMQxr7JMZJDuZMALYfKuPbQ2O8zwNvaWiyyThu7j+AnYYmJAmYTjX0067pDc+3oPpk+BTVm8TFtv/bCK9b5SBuOT7wJWAfqk/Ur0h6hrHH2RuNtc1Ra2jhs8BBVENuv6EKp3ZreoKGvmno29E0tn8TsDnVv9FK4FHbW4/wuua+XIdqyO4w4Ce2/yDpx23WPVqfPAFMk6SGQJgO/LqN9Y617uhA9gzWQpK2LQfcNi/T04HPALeXJvOA/yRpizI8cfI4bv5y4OOS9pE0hWqcfgVwWwfr+Jyk7SW9hWqs+ArbLwOXAJ+Q9NHyCXvdcjBy8/K631GNfQ9Zn+pNbxCYLOlUquMpnRhrm82aa2i2PlV/LKEK5H/ooJZrgB0kfVLVwf7jGDvI39PQ/oSy7duphmKWlwO765XfbUe9dvrx74AZJUCgGttfh6ovV0raH/jPHdQ+kl8ALwPHSpos6SCq4yXtGqu/o00Jg7XTcuB9wB2Snqf6z/8A1RszZdz7MuA+qoOIV4/Xhm0/TDXOfibVp/lPAJ+w/VIHq7mY6sDgb6nGuI8r615E9an6FKo3pUXA3/Ha3/F3gIMlPSPpDOCnwP+lOvj9GPB7xh5Waf59xtpms+Yaml1Uanmc6sD87S3ajFTLU8Cnga9RhcnWwK1jvOwnVMdxngE+D3yyjL2/DBxAdTziUap/q/Oo9lSgOi4EsETSPbaXU/07XF7W9Vlgbru1j/I7vUR10Pgo4Fmqv52raf806POB7ctQ14/faD0TmVYdqouI6C1JdwDn2L6w17VMJNkziIiekvQnkjYpw0SHAztT7dFFFyUMYkKSdF35olLz45Re17amkXTKCH15XZur2Aa4l2qY6ETgYNtP1FVvtJZhooiIyJ5BRESsod8zmDp1qmfMmNHrMiIi1ih33333U7b7Wi1bI8NgxowZDAwM9LqMiIg1iqTHRlqWYaKIiEgYREREwiAiIkgYREQECYOIiCBhEBERdCkMyuVxfylp2NUxJa1T7nm6UNX9Xmd0o6aIiHhNt/YMjgcWjLDsKOAZ2+8CvkV1r9SIiOii2sOg3ATk41TXSm/lIKobtwNcAexT7uAUERFd0o09g29T3XrwlRGWT6PccKTc53QpsHFzI0kzJQ1IGhgcHKyp1IiIianWMJB0APCk7bvf6Lpsn2u733Z/X1/LS2tERMTrVPeewQeAAyX9Bvg+8BFJlzS1eZxy0+5yn9YNqW7pFxERXVJrGNg+2fbmtmcAhwI32P5cU7O5wOHl+cGlTW6yEBHRRT25aqmk04AB23Opbmh9saSFwNNUoREREV3UtTCwfRNwU3l+asP83wOf7lYdERExXL6BHBERCYOIiEgYREQECYOIiCBhEBERJAwiIoKEQUREkDCIiAgSBhERQcIgIiJIGEREBAmDiIggYRARESQMIiKChEFERJAwiIgIag4DSetKulPSvZLmS/pKizZHSBqUNK88jq6zpoiIGK7uO52tAD5i+zlJU4CfS7rO9u1N7S6zfWzNtURExAhqDYNyY/vnyuSU8sjN7iMiVjO1HzOQNEnSPOBJ4Hrbd7Ro9ilJ90m6QtL0umuKiIhV1R4Gtl+2vSuwObCHpB2bmlwFzLC9M3A9MKfVeiTNlDQgaWBwcLDWmiMiJpqunU1k+1ngRmC/pvlLbK8ok+cB7xnh9efa7rfd39fXV2utERETTd1nE/VJ2qg8Xw/YF3ioqc2mDZMHAgvqrCkiIoar+2yiTYE5kiZRBc/ltq+WdBowYHsucJykA4GVwNPAETXXFBERTVSd8LNm6e/v98DAQK/LiIhYo0i623Z/q2X5BnJERCQMIiIiYRARESQMIiKChEFERJAwiIgIEgYREUHCICIiSBhERAQJg4iIIGEQEREkDCIigoRBRESQMIiICBIGERFBwiAiIkgYREQE9d8DeV1Jd0q6V9J8SV9p0WYdSZdJWijpDkkz6qwpIiKGq3vPYAXwEdu7ALsC+0nas6nNUcAztt8FfAv4es01RUREk1rDwJXnyuSU8mi+6fJBwJzy/ApgH0mqs66IiFhV7ccMJE2SNA94Erje9h1NTaYBiwBsrwSWAhvXXVdERLym9jCw/bLtXYHNgT0k7fh61iNppqQBSQODg4PjWmNMbDNOumbU5TvN2an2Gn52w1a1b6Ndi0+6ZcRlq1OdMb66djaR7WeBG4H9mhY9DkwHkDQZ2BBY0uL159rut93f19dXc7URERNL3WcT9UnaqDxfD9gXeKip2Vzg8PL8YOAG283HFSIiokaTa17/psAcSZOogudy21dLOg0YsD0XOB+4WNJC4Gng0JprioiIJrWGge37gN1azD+14fnvgU/XWUdERIwu30COiIiEQUREJAwiIoKEQUREkDCIiAgSBhERQcIgIiJIGEREBAmDiIggYRARESQMIiKChEFERJAwiIgIEgYREUHCICIiSBhERAQJg4iIoP57IE+XdKOkByXNl3R8izZ7S1oqaV55nNpqXRERUZ+2w0DSxe3Ma7ISONH29sCewDGStm/R7hbbu5bHae3WFBER46OTPYMdGifKTe7fM9oLbD9h+57yfDmwAJjWaZEREVGvMcNA0smSlgM7S1pWHsuBJ4GftLshSTOA3YA7WizeS9K9kq6TtEOL5UiaKWlA0sDg4GC7m42IiDaMGQa2/9H2+sA/2d6gPNa3vbHtk9vZiKS3AVcCJ9he1rT4HmBL27sAZwI/HqGOc2332+7v6+trZ7MREdGmye02tH2ypGnAlo2vs33zaK+TNIUqCC61/cMW613W8PxaSf8saartp9qtLSIi3pi2w0DS14BDgQeBl8tsAyOGgSQB5wMLbH9zhDabAL+zbUl7UO2tLGm3roiIeOPaDgPgz4BtbK/o4DUfAD4P3C9pXpl3CrAFgO1zgIOBv5K0EngRONS2O9hGRES8QZ2EwSPAFKDtMLD9c0BjtDkLOKuDOiIiYpyNGQaSzqQaDnoBmCfpZzQEgu3j6isvIiK6oZ09g4Hy825gbo21REREj4wZBrbndKOQiIjonU7OJrqfario0VKqPYev2s4ZQBERa6hODiBfR3VK6f8p04cCbwF+C8wGPjGulUVERNd0EgZ/anv3hun7Jd1je3dJnxvvwiIions6uVDdpPKlMAAkvReYVCZXjmtVERHRVZ3sGRwNXFCuMyRgGXC0pLcC/1hHcRER0R2dXJvoLmAnSRuW6aUNiy8f78IiIqJ72vnS2edsXyLpb5vmAzDSNYciImLN0c6ewVvLz/XrLCQiInqnnS+d/e9yV7Nltr/VhZoiIqLL2jqbyPbLwGdqriUiInqkk7OJbpV0FnAZ8PzQzKF7HEdExJqrkzDYtfw8rWGegY+MWzUREdETnZxa+uHRlks6PBe1i4hYM3XyDeSxHD+O64qIiC4azzAYdkczSdMl3SjpQUnzJQ0LDFXOkLRQ0n2Sdm9uExER9erkmMFYWt23eCVwou17JK0P3C3petsPNrTZH9i6PN4HfLf8jIiILql1z8D2E0NnG9leDiwApjU1Owi4yJXbgY0kbTqOdUVExBjGMwxuHW2hpBnAbsAdTYumAYsaphczPDCQNFPSgKSBwcHBN1hqrI52mrPTiMvO/sINq0yffsgBnH7IAW2t92c3bPXqYzxrW7DtdizYdruO1rP4pFuGzZs1a1bLtpvcOG/YvBknXcOMk67paJutamy376B1fa1qizVbR8NEkj4O7ACsOzTP9mnl57GjvO5twJXACbaXvZ5CbZ8LnAvQ39/fakgqIiJep7b3DCSdAxwC/DXVkNCngS3beN0UqiC41PYPWzR5HJjeML15mRcREV3SyTDR+20fBjxj+yvAXsC7R3uBqkubng8sGOXqpnOBw8pZRXsCS20/0UFdERHxBnUyTPRi+fmCpM2AJcBYB3o/AHye6haZ88q8U4AtAGyfA1wLfAxYCLwAHNlBTRERMQ46CYOrJW0E/BNwD9WppOeN9gLbP6fFWUZNbQwc00EdERExzjoJg/9lewVwpaSrqQ4i/76esiIiops6OWbwi6EntleU217+YpT2ERGxhmjntpebUJ33v56k3Xht2GcD4C011hYREV3SzjDRR4EjqE75bDwjaBnVweCIiFjDtXPbyznAHEmfsn1lF2qKiIgu6+SYwa2Szpd0HYCk7SUdVVNdERHRRZ2EwYXAT4HNyvSvgBPGu6CIiOi+TsJgqu3LgVcAbK8EXq6lqoiI6KpOwuB5SRtT7lswdOmIWqqKiIiu6uRLZ39LdR2hP5Z0K9AHHFxLVRER0VWdhMGDwI+orh+0HPgx1XGDiIhYw3UyTHQRsC3wD8CZVFcsvbiOoiIiors62TPY0fb2DdM3SnpwxNYREbHG6GTP4J5y0BgASe8DBsa/pIiI6LZ2rk10P9UZRFOA2yT9e5neEnio3vIiIqIb2hkmav/O2RERsUZq59pEj3WjkIiI6J1Ojhl0TNIFkp6U9MAIy/eWtFTSvPI4tc56IiKitU7OJno9ZgNnUZ2WOpJbbGcoKiKih2rdM7B9M/B0nduIiIg3rtYwaNNeku6VdJ2kHUZqJGmmpAFJA4ODg92sLyJirdfrMLgH2NL2LlTfav7xSA1tn2u733Z/X19ft+qLiJgQehoGtpfZfq48vxaYImlqL2uKiJiIehoGkjaRpPJ8j1LPkl7WFBExEdV6NpGk7wF7A1MlLQa+TPVNZmyfQ3UJ7L+StBJ4ETjUtuusKSIihqs1DGx/ZozlZ1GdehoRET3U6wPIERGxGkgYREREwiAiIhIGERFBwiAiIkgYREQECYOIiCBhEBERJAwiIoKEQUREkDCIiAgSBhERQcIgIiJIGEREBAmDiIggYRARESQMIiKCmsNA0gWSnpT0wAjLJekMSQsl3Sdp9zrriYiI1ureM5gN7DfK8v2BrctjJvDdmuuJiIgWag0D2zcDT4/S5CDgIlduBzaStGmdNUVExHC9PmYwDVjUML24zBtG0kxJA5IGBgcHu1Lc2mrBttuNuvz0Qw4YNm/xSbe84e3OOOmaVWfM2vDVpzvN2enVn0P1Ldh2u2G1nv2FG1apr1Wt7djkxnmvPobVNmvDVWprrHGovrH6sNnik255tQ+Hfs6aNWuVNj+7YatXH0N1NdbXTo2NGms8+ws3rPIYMtR/rWprrG+orkZDtY1oqL6mOof+rZvra3T6IQe8+ui2ofrq1NyXjYb9P+mSXodB22yfa7vfdn9fX1+vy4mIWKv0OgweB6Y3TG9e5kVERBf1OgzmAoeVs4r2BJbafqLHNUVETDiT61y5pO8BewNTJS0GvgxMAbB9DnAt8DFgIfACcGSd9URERGu1hoHtz4yx3MAxddYQERFj6/UwUURErAYSBhERkTCIiIiEQUREkDCIiAgSBhERQcIgIiJIGEREBAmDiIggYRARESQMIiKChEFERJAwiIgIEgYREUHCICIiSBhERAQJg4iIoAthIGk/SQ9LWijppBbLj5A0KGleeRxdd00REbGquu+BPAk4G9gXWAzcJWmu7Qebml5m+9g6a4mIiJHVvWewB7DQ9iO2XwK+DxxU8zYjIqJDdYfBNGBRw/TiMq/ZpyTdJ+kKSdNbrUjSTEkDkgYGBwfrqDUiYsJaHQ4gXwXMsL0zcD0wp1Uj2+fa7rfd39fX19UCIyLWdnWHweNA4yf9zcu8V9leYntFmTwPeE/NNUVERJO6w+AuYGtJfyTpzcChwNzGBpI2bZg8EFhQc00REdGk1rOJbK+UdCzwU2AScIHt+ZJOAwZszwWOk3QgsBJ4GjiizpoiImK4WsMAwPa1wLVN805teH4ycHLddURExMhWhwPIERHRYwmDiIhIGERERMIgIiJIGEREBAmDiIggYRARESQMIiKChEFERJAwiIgIEgYREUHCICIiSBhERAQJg4iIIGEQEREkDCIigoRBRETQhTCQtJ+khyUtlHRSi+XrSLqsLL9D0oy6a4qIiFXVGgaSJgFnA/sD2wOfkbR9U7OjgGdsvwv4FvD1OmuKiIjh6t4z2ANYaPsR2y8B3wcOampzEDCnPL8C2EeSaq4rIiIayHZ9K5cOBvazfXSZ/jzwPtvHNrR5oLRZXKZ/Xdo81bSumcDMMrkN8HBthVemAk+N2WriSH8Mlz5ZVfpjVatjf2xpu6/VgsndruT1sn0ucG63tidpwHZ/t7a3ukt/DJc+WVX6Y1VrWn/UPUz0ODC9YXrzMq9lG0mTgQ2BJTXXFRERDeoOg7uArSX9kaQ3A4cCc5vazAUOL88PBm5wnWNXERExTK3DRLZXSjoW+CkwCbjA9nxJpwEDtucC5wMXS1oIPE0VGKuDrg1JrSHSH8OlT1aV/ljVGtUftR5AjoiINUO+gRwREQmDiIiY4GEg6QJJT5bvOrRavqGkqyTdK2m+pCO7XWM3tdEfb5f0I0n3SbpT0o7drrGbJE2XdKOkB8u///Et2kjSGeVyKvdJ2r0XtXZDm/2xraRfSFoh6Yu9qLOb2uyTPy9/G/dLuk3SLr2odSwTOgyA2cB+oyw/BnjQ9i7A3sDp5ayotdVsRu+PU4B5tncGDgO+042iemglcKLt7YE9gWNaXE5lf2Dr8pgJfLe7JXZVO/3xNHAc8I1uF9cj7fTJo8Cf2N4J+B+spgeWJ3QY2L6Z6o93xCbA+uXyGG8rbVd2o7ZeaKM/tgduKG0fAmZIemc3ausF20/Yvqc8Xw4sAKY1NTsIuMiV24GNJG3a5VK7op3+sP2k7buAP/SgxK5rs09us/1Mmbyd6vtWq50JHQZtOAvYDvgP4H7geNuv9LaknroX+CSApD2ALVlN/7DHW7ma7m7AHU2LpgGLGqYXMzww1jqj9MeE1WafHAVc15WCOpQwGN1HgXnAZsCuwFmSNuhlQT32NapPvvOAvwZ+Cbzc04q6QNLbgCuBE2wv63U9vZb+GK6dPpH0Yaow+Ptu1tauNebaRD1yJPC18o3ohZIeBbYF7uxtWb1R/siPhOrAKdVY6CM9LapmkqZQ/Se/1PYPWzRp55Ira402+mPCaadPJO0MnAfsb3u1vNxO9gxG9+/APgBlbHwb1vI3v9FI2qjhAPrRwM1r8yfDEnjnAwtsf3OEZnOBw8pZRXsCS20/0bUiu6jN/phQ2ukTSVsAPwQ+b/tX3ayvExP6G8iSvkd1ltBU4HfAl4EpALbPkbQZ1Rk2mwKi2ku4pCfFdkEb/bEX1b0nDMwHjmo4MLbWkfRB4Baq40VDx4pOAbaAV/tEVMeW9gNeAI60PdCDcmvXZn9sAgwAG5Q2zwHbr60fGtrsk/OATwGPleUrV8ermU7oMIiIiEqGiSIiImEQEREJg4iIIGEQEREkDCIigoRBRESQMIgJTtKMkS7ZPUL7I8r3T17PtmZLOniMNttKmifpl5K2kvTZ17OtiE4lDCI6cwTVtarq8l+AK2zvRnWZi4RBdEWuTRQBkyVdCuxO9c3qw6iuVvtNqkuXP0UVAh8A+oFLJb0I7AX8HfAJYD3gNuAv3cY3OSW9p8X6dwNOAF6WtE9Z53blwoBzbH9rXH7biBbyDeSY0Mplhx8FPmj7VkkXUF2T/s+Ag2wPSjoE+Kjtv5B0E/DFoUtOSHqH7afL84uBy21fNcK2ZgNXAz8B/m2E9c8CnrP9DUl7l20dUNOvH/Gq7BlEwCLbt5bnl1BdW2ZH4Prq0kNMAka6+NyHJX0JeAvwDqo9i5Zh0GCbDtYf0RUJg4jqwnuNlgPzbe812oskrQv8M9Bve1H5VL9uG9tTO+uP6KYcQI6ALcoVWaE6YHs70Dc0T9IUSTuU5cuB9cvzoTf+p8rNTUY9U6jBw6Osv1HjtiJqlTCIqN6cj5G0AHg7cCbVG/vXJd1Ldbe795e2s4FzykHdFcC/AA8APwXuamdjtl8aZf2N7qM6mHyvpL95Pb9YRLtyADkiIrJnEBEROYAcMe4knU31nYRG37F9YS/qiWhHhokiIiLDRBERkTCIiAgSBhERQcIgIiKA/w+d9T1AHalO9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(m1_trace.posterior['beta_left'].values+m1_trace.posterior['beta_right'].values)\n",
    "plt.xlabel('beta_left')\n",
    "plt.ylabel('beta_right')\n",
    "plt.title('Sum of beta_left and beta_right')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95a3884b-8967-4983-9d12-215714db52a9",
   "metadata": {},
   "source": [
    "# Multicollinear milk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5b40d61f-8b6e-4c93-8474-320587e0f41f",
   "metadata": {},
   "outputs": [],
   "source": [
    "milk = pd.read_csv('./data/milk.csv', sep=';')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e84ab67c-0235-4d55-8372-0fc5161d98ee",
   "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>clade</th>\n",
       "      <th>species</th>\n",
       "      <th>kcal.per.g</th>\n",
       "      <th>perc.fat</th>\n",
       "      <th>perc.protein</th>\n",
       "      <th>perc.lactose</th>\n",
       "      <th>mass</th>\n",
       "      <th>neocortex.perc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Strepsirrhine</td>\n",
       "      <td>Eulemur fulvus</td>\n",
       "      <td>0.49</td>\n",
       "      <td>16.60</td>\n",
       "      <td>15.42</td>\n",
       "      <td>67.98</td>\n",
       "      <td>1.95</td>\n",
       "      <td>55.16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Strepsirrhine</td>\n",
       "      <td>E macaco</td>\n",
       "      <td>0.51</td>\n",
       "      <td>19.27</td>\n",
       "      <td>16.91</td>\n",
       "      <td>63.82</td>\n",
       "      <td>2.09</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Strepsirrhine</td>\n",
       "      <td>E mongoz</td>\n",
       "      <td>0.46</td>\n",
       "      <td>14.11</td>\n",
       "      <td>16.85</td>\n",
       "      <td>69.04</td>\n",
       "      <td>2.51</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Strepsirrhine</td>\n",
       "      <td>E rubriventer</td>\n",
       "      <td>0.48</td>\n",
       "      <td>14.91</td>\n",
       "      <td>13.18</td>\n",
       "      <td>71.91</td>\n",
       "      <td>1.62</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Strepsirrhine</td>\n",
       "      <td>Lemur catta</td>\n",
       "      <td>0.60</td>\n",
       "      <td>27.28</td>\n",
       "      <td>19.50</td>\n",
       "      <td>53.22</td>\n",
       "      <td>2.19</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           clade         species  kcal.per.g  perc.fat  perc.protein  \\\n",
       "0  Strepsirrhine  Eulemur fulvus        0.49     16.60         15.42   \n",
       "1  Strepsirrhine        E macaco        0.51     19.27         16.91   \n",
       "2  Strepsirrhine        E mongoz        0.46     14.11         16.85   \n",
       "3  Strepsirrhine   E rubriventer        0.48     14.91         13.18   \n",
       "4  Strepsirrhine     Lemur catta        0.60     27.28         19.50   \n",
       "\n",
       "   perc.lactose  mass  neocortex.perc  \n",
       "0         67.98  1.95           55.16  \n",
       "1         63.82  2.09             NaN  \n",
       "2         69.04  2.51             NaN  \n",
       "3         71.91  1.62             NaN  \n",
       "4         53.22  2.19             NaN  "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "milk.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "fe5f4b34-0d52-4c94-ade5-db0b967808d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def standarize(col):\n",
    "    return (col - np.mean(col))/np.std(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "af8dcd7f-b9f0-4ff5-84f9-d59be1d7190a",
   "metadata": {},
   "outputs": [],
   "source": [
    "milk['standarized_kcal'] = standarize(milk[\"kcal.per.g\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "827f33f5-637f-4135-8f2c-e9e73b247598",
   "metadata": {},
   "outputs": [],
   "source": [
    "milk['standarized_fat'] = standarize(milk[\"perc.fat\"])\n",
    "milk['standarized_lactose'] = standarize(milk[\"perc.lactose\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "87fa453c-5e3f-43ad-ad38-81552de96f5a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, beta_fat, 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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 19 seconds.\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, beta_lactose, 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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 19 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m1:\n",
    "    alpha = pm.Normal('alpha', 0, 0.2)\n",
    "    beta = pm.Normal('beta_fat', 0, 0.5)\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    mu = pm.Deterministic('mu', alpha + beta * milk['standarized_fat'])\n",
    "    sdr = pm.Normal('standarized_kcal', mu, sd=sigma, observed=milk['standarized_kcal'])\n",
    "    m1_trace = pm.sample(return_inferencedata=True)\n",
    "    m1_prior_predictive = pm.sample_prior_predictive()\n",
    "    \n",
    "with pm.Model() as m2:\n",
    "    alpha = pm.Normal('alpha', 0, 0.2)\n",
    "    beta = pm.Normal('beta_lactose', 0, 0.5)\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    mu = pm.Deterministic('mu', alpha + beta * milk['standarized_lactose'])\n",
    "    sdr = pm.Normal('standarized_kcal', mu, sd=sigma, observed=milk['standarized_kcal'])\n",
    "    m2_trace = pm.sample(return_inferencedata=True)\n",
    "    m2_prior_predictive = pm.sample_prior_predictive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "4c3b60c6-d864-44a2-a98f-5b117bb9098a",
   "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>beta_fat</th>\n",
       "      <td>0.856</td>\n",
       "      <td>0.091</td>\n",
       "      <td>0.680</td>\n",
       "      <td>1.019</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4229.0</td>\n",
       "      <td>2370.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0.085</td>\n",
       "      <td>-0.155</td>\n",
       "      <td>0.154</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>3875.0</td>\n",
       "      <td>2751.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  \\\n",
       "beta_fat  0.856  0.091   0.680    1.019      0.001    0.001    4229.0   \n",
       "alpha     0.001  0.085  -0.155    0.154      0.001    0.001    3875.0   \n",
       "\n",
       "          ess_tail  r_hat  \n",
       "beta_fat    2370.0    1.0  \n",
       "alpha       2751.0    1.0  "
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m1_trace, var_names=['beta_fat', 'alpha'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "81a6106e-a266-41c1-831b-48615b0dc29e",
   "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>beta_lactose</th>\n",
       "      <td>-0.899</td>\n",
       "      <td>0.078</td>\n",
       "      <td>-1.036</td>\n",
       "      <td>-0.748</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4820.0</td>\n",
       "      <td>2943.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>-0.002</td>\n",
       "      <td>0.072</td>\n",
       "      <td>-0.134</td>\n",
       "      <td>0.132</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4203.0</td>\n",
       "      <td>2628.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  \\\n",
       "beta_lactose -0.899  0.078  -1.036   -0.748      0.001    0.001    4820.0   \n",
       "alpha        -0.002  0.072  -0.134    0.132      0.001    0.001    4203.0   \n",
       "\n",
       "              ess_tail  r_hat  \n",
       "beta_lactose    2943.0    1.0  \n",
       "alpha           2628.0    1.0  "
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m2_trace, var_names=['beta_lactose', 'alpha'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "144cf8d3-c5d1-4579-846b-705db9c8f150",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, beta_fat, beta_lactose, 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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:06<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 20 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m3:\n",
    "    alpha = pm.Normal('alpha', 0, 0.2)\n",
    "    beta_lactose = pm.Normal('beta_lactose', 0, 0.5)\n",
    "    beta_fat = pm.Normal('beta_fat', 0, 0.5)\n",
    "\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    mu = pm.Deterministic('mu', alpha + beta_lactose * milk.standarized_lactose + beta_fat * milk.standarized_fat)\n",
    "    sdr = pm.Normal('standarized_kcal', mu, sd=sigma, observed=milk['standarized_kcal'])\n",
    "    m3_trace = pm.sample(return_inferencedata=True)\n",
    "    m3_prior_predictive = pm.sample_prior_predictive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "071e4cd9-29c8-41db-b4f0-bf5d37185b4e",
   "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>beta_lactose</th>\n",
       "      <td>-0.669</td>\n",
       "      <td>0.202</td>\n",
       "      <td>-1.064</td>\n",
       "      <td>-0.308</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>1355.0</td>\n",
       "      <td>1572.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_fat</th>\n",
       "      <td>0.249</td>\n",
       "      <td>0.201</td>\n",
       "      <td>-0.139</td>\n",
       "      <td>0.622</td>\n",
       "      <td>0.006</td>\n",
       "      <td>0.004</td>\n",
       "      <td>1324.0</td>\n",
       "      <td>1536.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.001</td>\n",
       "      <td>0.076</td>\n",
       "      <td>-0.130</td>\n",
       "      <td>0.147</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>2511.0</td>\n",
       "      <td>2034.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  \\\n",
       "beta_lactose -0.669  0.202  -1.064   -0.308      0.005    0.004    1355.0   \n",
       "beta_fat      0.249  0.201  -0.139    0.622      0.006    0.004    1324.0   \n",
       "alpha         0.001  0.076  -0.130    0.147      0.002    0.001    2511.0   \n",
       "\n",
       "              ess_tail  r_hat  \n",
       "beta_lactose    1572.0    1.0  \n",
       "beta_fat        1536.0    1.0  \n",
       "alpha           2034.0    1.0  "
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m3_trace, var_names=['beta_lactose', 'beta_fat', 'alpha'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13078c53-2f36-43e7-9f93-f0625e147672",
   "metadata": {},
   "source": [
    "# Post treatment bias"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39248847-708e-48cb-bd60-eed513207d4f",
   "metadata": {},
   "source": [
    "In general people take care about consequences of not including variables in the model, but less often about consequences of including them. Post-treatment bias is one case of problematic inclution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e84c2f86-6f73-4c67-81f5-61347c498369",
   "metadata": {},
   "source": [
    "## Fungus example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "8ca25983-c756-467d-8f97-3d05b8f77879",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "heights_0 = np.random.normal(10,2,size=n)\n",
    "treatment = np.random.binomial(n=1, p=1/2, size=n)\n",
    "fungus = np.random.binomial(n=1, size=n, p=0.5 - treatment * 0.4)\n",
    "heights_1 = heights_0 + np.random.normal(5-3*fungus, size=n)\n",
    "heights_df = pd.DataFrame({\n",
    "'heights_0': heights_0,\n",
    "'treatment': treatment,\n",
    "'fungus': fungus,\n",
    "'heights_1': heights_1,})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21a937b7-0537-48f0-89bd-86d708a92fbd",
   "metadata": {},
   "source": [
    "## Growth only wrt initial height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "95c26e74-dc7d-4dbf-89e1-bea99574ee35",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [sigma, proportion]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 20 seconds.\n"
     ]
    }
   ],
   "source": [
    "h0 = np.asarray(heights_df.heights_0)\n",
    "with pm.Model() as m4:\n",
    "    proportion = pm.Lognormal('proportion', 0, 0.25)\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    sdr = pm.Normal('growth_proportion', mu=proportion*h0 , sigma=sigma, observed=heights_df['heights_1'])\n",
    "    m4_trace = pm.sample(return_inferencedata=True)\n",
    "    m4_prior_predictive = pm.sample_prior_predictive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "fcdf9c00-bd18-4851-9c60-6d11de258291",
   "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>proportion</th>\n",
       "      <td>1.388</td>\n",
       "      <td>0.020</td>\n",
       "      <td>1.349</td>\n",
       "      <td>1.423</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>3717.0</td>\n",
       "      <td>2617.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.975</td>\n",
       "      <td>0.142</td>\n",
       "      <td>1.724</td>\n",
       "      <td>2.254</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>3611.0</td>\n",
       "      <td>3013.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  \\\n",
       "proportion  1.388  0.020   1.349    1.423      0.000    0.000    3717.0   \n",
       "sigma       1.975  0.142   1.724    2.254      0.002    0.002    3611.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "proportion    2617.0    1.0  \n",
       "sigma         3013.0    1.0  "
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m4_trace, var_names=['proportion', 'sigma'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "370c64a4-367e-4984-b8a9-b8bc61115bf0",
   "metadata": {},
   "source": [
    "## Including Fungus and treatment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "f9bc9839-302f-44e0-8cd9-521350e93593",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [alpha, beta_fungus, beta_treatment, sigma]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.\n",
      "The acceptance probability does not match the target. It is 0.8958949075991842, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The acceptance probability does not match the target. It is 0.8847142810519663, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The acceptance probability does not match the target. It is 0.8835652701820682, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m5:\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    beta_treatment = pm.Normal('beta_treatment', 0, 0.5)\n",
    "    beta_fungus = pm.Normal('beta_fungus', 0, 0.5)\n",
    "    alpha = pm.Lognormal('alpha', 0, 0.25)\n",
    "    \n",
    "    proportion = pm.Deterministic('proportion', alpha + beta_treatment * heights_df.treatment + beta_fungus * heights_df.fungus)\n",
    "    \n",
    "    sdr = pm.Normal('growth_proportion', mu=proportion*h0 , sigma=sigma, observed=heights_df['heights_1'])\n",
    "    m5_trace = pm.sample(return_inferencedata=True)\n",
    "    m5_prior_predictive = pm.sample_prior_predictive()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "31123f7b-afa4-4c40-8535-9ce781e26cca",
   "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>beta_treatment</th>\n",
       "      <td>-0.015</td>\n",
       "      <td>0.033</td>\n",
       "      <td>-0.077</td>\n",
       "      <td>0.046</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1462.0</td>\n",
       "      <td>1944.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_fungus</th>\n",
       "      <td>-0.303</td>\n",
       "      <td>0.033</td>\n",
       "      <td>-0.365</td>\n",
       "      <td>-0.239</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1595.0</td>\n",
       "      <td>2003.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.342</td>\n",
       "      <td>0.096</td>\n",
       "      <td>1.166</td>\n",
       "      <td>1.527</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>3020.0</td>\n",
       "      <td>2406.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  \\\n",
       "beta_treatment -0.015  0.033  -0.077    0.046      0.001    0.001    1462.0   \n",
       "beta_fungus    -0.303  0.033  -0.365   -0.239      0.001    0.001    1595.0   \n",
       "sigma           1.342  0.096   1.166    1.527      0.002    0.001    3020.0   \n",
       "\n",
       "                ess_tail  r_hat  \n",
       "beta_treatment    1944.0    1.0  \n",
       "beta_fungus       2003.0    1.0  \n",
       "sigma             2406.0    1.0  "
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m5_trace, var_names=['beta_treatment', 'beta_fungus', 'sigma'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164d94ae-67c0-4dfd-89e0-204471fe4fee",
   "metadata": {},
   "source": [
    "Treatment is not associated with growth, although we know it certainly is given the way data was constructed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "464a98ce-d4b3-4006-a2ca-7b5e8b6142a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [alpha, beta_treatment, sigma]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:03<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 16 seconds.\n",
      "The acceptance probability does not match the target. It is 0.8890655723841803, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The acceptance probability does not match the target. It is 0.8915858534704867, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m6:\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    beta_treatment = pm.Normal('beta_treatment', 0, 0.5)\n",
    "    alpha = pm.Lognormal('alpha', 0, 0.25)\n",
    "    \n",
    "    proportion = pm.Deterministic('proportion', alpha + beta_treatment * heights_df.treatment)\n",
    "    \n",
    "    sdr = pm.Normal('growth_proportion', mu=proportion*h0 , sigma=sigma, observed=heights_df['heights_1'])\n",
    "    m6_trace = pm.sample(return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "219416fb-0fca-4ce7-8d8c-59fa61a6458d",
   "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>beta_treatment</th>\n",
       "      <td>0.147</td>\n",
       "      <td>0.037</td>\n",
       "      <td>0.074</td>\n",
       "      <td>0.216</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1761.0</td>\n",
       "      <td>1961.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.839</td>\n",
       "      <td>0.135</td>\n",
       "      <td>1.594</td>\n",
       "      <td>2.094</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>2740.0</td>\n",
       "      <td>2672.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>1.297</td>\n",
       "      <td>0.030</td>\n",
       "      <td>1.241</td>\n",
       "      <td>1.354</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1736.0</td>\n",
       "      <td>2052.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  \\\n",
       "beta_treatment  0.147  0.037   0.074    0.216      0.001    0.001    1761.0   \n",
       "sigma           1.839  0.135   1.594    2.094      0.003    0.002    2740.0   \n",
       "alpha           1.297  0.030   1.241    1.354      0.001    0.001    1736.0   \n",
       "\n",
       "                ess_tail  r_hat  \n",
       "beta_treatment    1961.0    1.0  \n",
       "sigma             2672.0    1.0  \n",
       "alpha             2052.0    1.0  "
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m6_trace, var_names=['beta_treatment', 'sigma', 'alpha'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "a7b67e14-b02b-4cc0-af4d-ee7947eaacdc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAA2CAYAAAAcX4Z+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAALtUlEQVR4nO2de7BWVRXAfwvwAhJcFVADUjFA5PL0ivgItTST0VQcDalJJV+IzWhTNpZT2dhDx9Q0Bx9JZRqO6WTQwzIlTS0Rr4KioqKCgQ9A8U0IuPpj7StXuvf7zmufY971mzkz986cvdZeZ59v7dda+4iq4jiOE4MuVVfAcZyPLu5gHMeJhjsYx3Gi4Q7GcZxouINxHCca7mAcx4mGOxjHcaLhDsZxnGi4g3EcJxruYBzHiYY7GMdxotEtlmAR6QXsBTQDY4FGQIA3gEeBB4EHVPX1WHVw8iEiAjQBe4ZrENAArAeeZXMbPlNZJXMiIj0x25qBccB22Hv6FvAYZuN8VX21skrmRES2ByZgNo4AegGbgDXAw0AL8KCqvlu47qKTHUVkJHA6MBV4Amugh4FXwi3bYA6nGRgD3ArMVNUHC62IkxkR6QMcD8wAugP3Yy/hc8C7QA9gGNaG+wLLgZnAzar6nyrqnBYRGYK9pycAz2Dv6UPAakCBPsBozMZm4DbMxvv0/yBDOHQOnwHOAD4NPIC14aPAm9jsZUfMttbOYxZwtaouL6wiqlrIhTXI1cCLwHnAwARl+gPnAM8DNwL9iqqPX5naUIApwEvATcABhE6oRpluwJHA7dio5sCq7ahT357ARZgjuRAYnKDMNsCZwNPAH5O82xXbOBiYh43ApgO9E5TZDfgpNqo5D2gopC4FGbQf1ov9HGjMUH5r4OLgnA6tuoE64wX0Bm4BHgcmZJRxOLACuBzoVrVN7dRvDPBkcJ79M5RvCD++VcAXqrangzp+JTiJbwBdM5QfGJzoQmC33PUpwKBDQm9weAGyJobec0rVDdWZrtBDPxA6iB45ZW2LTSd+X1QvWJCN+wTH8MUCZO0B/BuYXrVdW9TrHGApMDynHAFOCx3+6FyyclZk7+Bc9ivwIY0KTsZHMiVc2BrLvcBl9aZDKWQ2AHOA2UXJzFmfJuBlYFKBMnfFpvZTq7Yv1Of0MIX7eIEypwArSTCN7FBGDuW9sMWxyREe1kTgBXxNJvoF/ACYC3QpWG4PYBFwQsX2bYUt3p4cQfbo0MHuXLGNI8K06JMRZJ8N3JX1/ci8iyQilwPbquqXMwmoL/8SzBtPjSHfARFpBv4MjFXVFyPIH4st/o5T1ZVFy09Yh+9g06PDNOvLXlv+t7DdmkNiyE+gvxtwH/ArVb0ygvyuwD3AbFW9InX5LM9ERAYDC4BhGik+QES2BpYAR6tvYUdBRO4AblTVWRF1XAj0UdXTY+moobsvtiYxSlVXRNLRDQvDOFtV/xJDRx39U4GvAhNV9b1IOoZjTmwnVX07TdmskbynAdfFci4AqvoOcCU2t3QKJrw0o4AbIqu6HJgSYmvKZhowN5ZzAVDVjdj6VVXv6QzgkljOBUBVlwD/BI7LUjjLnHYVMLTGPW9t8f+JwBXh7+7YNuFSYD6wSw05/YG1ZNj69qtuO14M/DhjG+6PrWtsBI5JoOtmSt5xwXZCngb2TlFmE7Y923rtkrBcLyyQdFDJNjZhi7BbJby/bxvbXgplW/+vueMHTAIWpK1jlhFME7BGVZ/OUBbgJGCtqg4BLsWCndpFVVcDjwDjM+pyOuZAbHE3C89jDmd2wvv/gAXtlckALD1lfooy61R1bJtrWZJCatOGu7DNiTI5ALhNVTckuVlVX2m1DbgKuLSNrfXSBO4ARohI7zQVzOJgmrGQ46wcCVwX/r4FOCiENXdECxbK7BSEiHQHdsd2eVKjqstU9REg6bC8ijbcE2jR0P2WQAv22yiTvL/FxAQnthhL80lMlmTHUdiQqhY9RaTtPduxubcciAUpoaobReR1bOi2pgNZC4FDM9SzU9Nmh+gprHddgL2ML2Bh4cvU1rk6olYbpmUJMCgs3K/DQtn3xBLw9sdyYiao6gsZ5bdHkvd0S9ra/JyqTk5RdiFwFoCINAAjMQfwKWwXa7mqfjZlfeoxGri2YJm1WBh03pO0QBYH0wd4rc4968IwDAAROZHsPdha4AgR+dAnmH0IUWB77AV/C4tNWY8lLdYbVhfWhqq6SUTexXYihmHrIxuAj7F5FL2y9kA2E+ekvP8DNqdkLdAkIk9hDvQdoCu2PgMwNMI7/HbQWxZrsd9/YrI4mLwPaSXwCWBF2OJrZHOmdXsI1li1elvnf+mPtZVgGdCCvfBrsdFi3wrq9BKwA9Av1GkDtujfyuoCdTUUKCsJrc95FTYiawj/b8KeOxRrH9iGS+FeuQZCyt9/FgfTOqXJylwsRf5fwDHAvDrz5L7AXzVSQN9HlXBsxk1Yev492PRokaquE5Em4Hcl1mUrzJEcHfT3wc5eacYWKncCPqeqqwrU+W3KdaJ9gSWqOimsKQ7E7BuP2bhSVdNv89ZARO6nfBuXpSmQxcEsIt+ayCzgehFZCrxK/b31caSfS3d6VHUxtuPXHk9hayK9VfXNtLJFZDx2js+2wOdF5Puq2pEuQj2eVdV1oW5vAHeH65K0+hOyEPhaJNntMY6waB46zBXhmhNR56Kg996IOtoyDvv9JiZ1JG/oGW9V1aGpCmZERO4DzlXVu8rQ11kQkflY9Ok/StB1EnZOTGmjUBHZATvwrG8ZO0kiMge4QVVvjq2rjc5TsUTjE0rQ1R2bXversznwAbJsUz8BbC0iozKUTYWIDMK2UxfE1tUJuRNIs0uSh8lBX2mo6svYjln0+BsRacR2w6I76y2YB0wSkR4l6DoM2/ZPtRaa2sGo6ibgGsoJjT4FS7JKlf/gJOIa4PiwdRyNkLe2N/DbmHo64EoslD42x2PrhC+XoOt9VHUpNhU8pgR1M7DnmYqsyY4DsKCbJo2QhRt0NGLxEwer6mMxdHR2RGQucKeqXhZRxxXAelX9eiwdNXT3wRYlJ+SIPK+nowe2kH5SGdPNdvQfCXwP2EstLyqGjj2wmKqdVXV9qrJZp6cicj4WdHNUjDmuiFwLbFTV6UXLdgwR2R3bYRqvqs9FkD8R28kapaq1QhGiISJnYj38ARohIVBELgCGAMeWGDXcVn8X4G/A7araYdpNDvkN2BLFT1T1+tQCciRadcc897SsMmrIPhzreeoeVuxX7mf9TSzSN1HCXAq5jViy4VEV29cFWxs5O4Ls/QixPRXbuAsW25TreMsOZP8IyyXLdDJhXuUjsaMIDyu40VYD+1bZaJ3lwoLA/gT8hoIO6sYidO8Gfla1faE+g7Et4y8VKHN0cC6FHcOZsz7HBRuHFCjzDOzUyh0zyyigEhOCk5mW1cu1kTU5OJdDqm6wznRhX3W4HQu+2yanrEFYEOUsCj6GM2e9Wo82OCtvvYCDwjt/bNV2bVGvk4ON++SU0xX4LpZSsmsuWQUZNhI7VmEOGQ4dxhLprsfOiEl8fodfxV1hyjsTS0RN3StjYeTTsFD5c/N2NpFs3BX7iNw8MhxkHUZmrc/o4Krt6aCOR4SR1UVAzwzlh4cOYh4FfP+pSMMagPOx6NyZ2MJevTLDsEjONdipYL2qbqDOfoXe+Vng78Cx1FmbwZL5TsGOjXyICOsABdvXFTvI+lXgl9gCd70yOwM/DD/cX+Qd5ZVgY38sLGBFGIkMqHO/YFnfs7G8wBlFjT5jfDp2ADZUOxVLzW9h86djFQsvH4vlaWwXGuxqTXi4jxOfkDt0FPaijcXasAVzPBv44Kdjx2C93UzgDo14dGORiEg/7CNl07GF4Lafjn0Pyxoeg9k4APg1cJWqPllJhTMQgmFbP+O8HLPxUSyzvguWeNqaL/U21obXqeprhdWhaAfzvmA7jXw4m1/CRsxTvolNp1qAxzXhaVxONYSQ+9bvMw/CRqrrMWfTAjykqmUeGVAoYZt3KGbfOKwD7IL9CBdjNi7WlPEfHyZEpCebv7M9Aht1bsQ6/YcxG5/RCM4gmoNxHMfJ+lUBx3GcuriDcRwnGu5gHMeJhjsYx3Gi4Q7GcZxouINxHCca7mAcx4mGOxjHcaLhDsZxnGi4g3EcJxruYBzHicZ/ASRfTq11cJVHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 266.457x39.685 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import daft\n",
    "\n",
    "from causalgraphicalmodels import CausalGraphicalModel\n",
    "\n",
    "plant_dag = CausalGraphicalModel(\n",
    "    nodes=[\"H0\", \"H1\", \"F\", \"T\"], edges=[(\"H0\", \"H1\"), (\"F\", \"H1\"), (\"T\", \"F\")]\n",
    ")\n",
    "pgm = daft.PGM()\n",
    "coordinates = {\"H0\": (0, 0), \"T\": (4, 0), \"F\": (3, 0), \"H1\": (2, 0)}\n",
    "for node in plant_dag.dag.nodes:\n",
    "    pgm.add_node(node, node, *coordinates[node])\n",
    "for edge in plant_dag.dag.edges:\n",
    "    pgm.add_edge(*edge)\n",
    "pgm.render(backgo)\n",
    "plt.gca().invert_yaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47151565-62de-4456-9c30-2a80256a9af3",
   "metadata": {},
   "source": [
    "# Collider bias"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "4100e54f-f585-4a99-bd36-6c51a6a8c320",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This one just copy pasting from pymc port\n",
    "def inv_logit(x):\n",
    "    return np.exp(x) / (1 + np.exp(x))\n",
    "\n",
    "\n",
    "def sim_happiness(N_years=100, seed=1234):\n",
    "    np.random.seed(seed)\n",
    "\n",
    "    popn = pd.DataFrame(np.zeros((20 * 65, 3)), columns=[\"age\", \"happiness\", \"married\"])\n",
    "    popn.loc[:, \"age\"] = np.repeat(np.arange(65), 20)\n",
    "    popn.loc[:, \"happiness\"] = np.repeat(np.linspace(-2, 2, 20), 65)\n",
    "    popn.loc[:, \"married\"] = np.array(popn.loc[:, \"married\"].values, dtype=\"bool\")\n",
    "\n",
    "    for i in range(N_years):\n",
    "        # age population\n",
    "        popn.loc[:, \"age\"] += 1\n",
    "        # replace old folk with new folk\n",
    "        ind = popn.age == 65\n",
    "        popn.loc[ind, \"age\"] = 0\n",
    "        popn.loc[ind, \"married\"] = False\n",
    "        popn.loc[ind, \"happiness\"] = np.linspace(-2, 2, 20)\n",
    "\n",
    "        # do the work\n",
    "        elligible = (popn.married == 0) & (popn.age >= 18)\n",
    "        marry = np.random.binomial(1, inv_logit(popn.loc[elligible, \"happiness\"] - 4)) == 1\n",
    "        popn.loc[elligible, \"married\"] = marry\n",
    "\n",
    "    popn.sort_values(\"age\", inplace=True, ignore_index=True)\n",
    "\n",
    "    return popn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "00ac51ee-a078-40dc-9eee-a1ae9a6e76bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "popn = sim_happiness()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "f7fcdf7c-25c1-4a82-9cdb-5c47d031e74a",
   "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>age</th>\n",
       "      <th>happiness</th>\n",
       "      <th>married</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.947368</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>1.789474</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>1.578947</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>1.368421</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   age  happiness  married\n",
       "0    0   0.947368    False\n",
       "1    0   2.000000    False\n",
       "2    0   1.789474    False\n",
       "3    0   1.578947    False\n",
       "4    0   1.368421    False"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popn.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "3b9a2c5c-ba6c-4ffa-89c1-1070f788a751",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hy/9qct0wrx5pl34rzz153ys6sw0000gp/T/ipykernel_9512/1174518190.py:2: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  popn['age'] = popn.age - 18 / (65-18)\n"
     ]
    }
   ],
   "source": [
    "popn = popn[popn.age > 17]\n",
    "popn['age'] = popn.age - 18 / (65-18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "93ee7b77-da2c-45e1-856a-c5df1a2b960c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hy/9qct0wrx5pl34rzz153ys6sw0000gp/T/ipykernel_9512/2888753241.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [alpha_married, alpha_single, beta_age, sigma]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:06<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 20 seconds.\n"
     ]
    }
   ],
   "source": [
    "popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
    "with pm.Model() as m7:\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    beta_age = pm.Normal('beta_age', 0, 2)\n",
    "    alpha = np.array([pm.Normal('alpha_single', 0, 1), pm.Normal('alpha_married', 0, 1)])\n",
    "    alpha = pm.Deterministic(\"alpha\", alpha[0] * (1-popn.marriage_indicator) + alpha[1] * popn.marriage_indicator)\n",
    "    mu = pm.Deterministic('mu', alpha + beta_age * popn['age'])\n",
    "    sdr = pm.Normal('happiness', mu=mu , sigma=sigma, observed=popn['happiness'])\n",
    "    m7_trace = pm.sample(return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "8aca61cd-0b06-44cc-8d53-dfb357e8fbd7",
   "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_single</th>\n",
       "      <td>0.065</td>\n",
       "      <td>0.102</td>\n",
       "      <td>-0.147</td>\n",
       "      <td>0.243</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1505.0</td>\n",
       "      <td>1795.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha_married</th>\n",
       "      <td>1.479</td>\n",
       "      <td>0.127</td>\n",
       "      <td>1.227</td>\n",
       "      <td>1.708</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1545.0</td>\n",
       "      <td>1827.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_age</th>\n",
       "      <td>-0.015</td>\n",
       "      <td>0.003</td>\n",
       "      <td>-0.020</td>\n",
       "      <td>-0.010</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>1442.0</td>\n",
       "      <td>1709.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  \\\n",
       "alpha_single   0.065  0.102  -0.147    0.243      0.003    0.002    1505.0   \n",
       "alpha_married  1.479  0.127   1.227    1.708      0.003    0.002    1545.0   \n",
       "beta_age      -0.015  0.003  -0.020   -0.010      0.000    0.000    1442.0   \n",
       "\n",
       "               ess_tail  r_hat  \n",
       "alpha_single     1795.0    1.0  \n",
       "alpha_married    1827.0    1.0  \n",
       "beta_age         1709.0    1.0  "
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m7_trace, var_names = ['alpha_single', 'alpha_married', 'beta_age'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "1d48057c-c454-4fe3-86e4-b03924850757",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hy/9qct0wrx5pl34rzz153ys6sw0000gp/T/ipykernel_9512/3458000315.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [alpha, beta_age, sigma]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:05<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 19 seconds.\n",
      "The acceptance probability does not match the target. It is 0.8867005657326522, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The acceptance probability does not match the target. It is 0.8987936880096942, but should be close to 0.8. Try to increase the number of tuning steps.\n"
     ]
    }
   ],
   "source": [
    "popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
    "with pm.Model() as m8:\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    beta_age = pm.Normal('beta_age', 0, 2)\n",
    "    alpha = pm.Normal('alpha', 0, 1)\n",
    "    mu = pm.Deterministic('mu', alpha + beta_age * popn['age'])\n",
    "    sdr = pm.Normal('happiness', mu=mu , sigma=sigma, observed=popn['happiness'])\n",
    "    m8_trace = pm.sample(return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "6e70df60-dd94-4023-bef6-9652be0339d3",
   "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.006</td>\n",
       "      <td>0.125</td>\n",
       "      <td>-0.218</td>\n",
       "      <td>0.253</td>\n",
       "      <td>0.003</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1590.0</td>\n",
       "      <td>1723.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_age</th>\n",
       "      <td>-0.000</td>\n",
       "      <td>0.003</td>\n",
       "      <td>-0.006</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>1529.0</td>\n",
       "      <td>1732.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  \\\n",
       "alpha     0.006  0.125  -0.218    0.253      0.003    0.003    1590.0   \n",
       "beta_age -0.000  0.003  -0.006    0.005      0.000    0.000    1529.0   \n",
       "\n",
       "          ess_tail  r_hat  \n",
       "alpha       1723.0    1.0  \n",
       "beta_age    1732.0    1.0  "
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m8_trace, var_names = ['alpha', 'beta_age'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "9cfe7552-553c-48fc-bff4-a3b0bf6ecd2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "n=200\n",
    "g_on_p =1\n",
    "g_on_c = 0\n",
    "p_on_c = 2\n",
    "u_on_p = 2\n",
    "u_on_c = 2\n",
    "u = 2 * np.random.binomial(n=1,p=0.5, size=n) - 1 \n",
    "g = np.random.normal(size=n)\n",
    "p = np.random.normal(g_on_p * g + u_on_p * u)\n",
    "c = np.random.normal(g_on_c * g + u_on_c * u)\n",
    "triads = pd.DataFrame({'uobserved':u, 'granparents':g, 'parents':p, 'childrens':c})\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "afda1132-184b-4828-9616-72da3c70a2ef",
   "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>uobserved</th>\n",
       "      <th>granparents</th>\n",
       "      <th>parents</th>\n",
       "      <th>childrens</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>-0.432948</td>\n",
       "      <td>2.194821</td>\n",
       "      <td>1.446184</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1</td>\n",
       "      <td>-0.000371</td>\n",
       "      <td>-1.662014</td>\n",
       "      <td>-2.697246</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1.077268</td>\n",
       "      <td>2.517555</td>\n",
       "      <td>2.079002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-1</td>\n",
       "      <td>-1.247351</td>\n",
       "      <td>-3.055249</td>\n",
       "      <td>-2.345441</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1.739514</td>\n",
       "      <td>1.402570</td>\n",
       "      <td>1.153817</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   uobserved  granparents   parents  childrens\n",
       "0          1    -0.432948  2.194821   1.446184\n",
       "1         -1    -0.000371 -1.662014  -2.697246\n",
       "2          1     1.077268  2.517555   2.079002\n",
       "3         -1    -1.247351 -3.055249  -2.345441\n",
       "4          1     1.739514  1.402570   1.153817"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "triads.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4832ebcb-c1d4-4cfd-a34d-363c938d7a16",
   "metadata": {},
   "source": [
    "Since there's influence of G on C through P, we control on P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "c1f85a1e-9577-4dcb-a78a-5088b4861dff",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hy/9qct0wrx5pl34rzz153ys6sw0000gp/T/ipykernel_9512/1542437129.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [alpha, beta_grandparent, beta_parent, sigma]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:03<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.\n"
     ]
    }
   ],
   "source": [
    "popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
    "with pm.Model() as m9:\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    beta_parent = pm.Normal('beta_parent', 0, 1)\n",
    "    beta_grandparent = pm.Normal('beta_grandparent', 0, 1)\n",
    "    alpha = pm.Normal('alpha', 0, 1)\n",
    "    mu = pm.Deterministic('mu', alpha + beta_parent * triads['parents'] + beta_grandparent * triads['granparents'])\n",
    "    sdr = pm.Normal('educational_achievement', mu=mu , sigma=sigma, observed=triads['childrens'])\n",
    "    m9_trace = pm.sample(return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "479c945c-1982-404a-882c-c5d264667c18",
   "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>beta_parent</th>\n",
       "      <td>0.824</td>\n",
       "      <td>0.044</td>\n",
       "      <td>0.747</td>\n",
       "      <td>0.915</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.000</td>\n",
       "      <td>4643.0</td>\n",
       "      <td>3276.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_grandparent</th>\n",
       "      <td>-0.656</td>\n",
       "      <td>0.111</td>\n",
       "      <td>-0.858</td>\n",
       "      <td>-0.443</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>3933.0</td>\n",
       "      <td>3153.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>-0.000</td>\n",
       "      <td>0.101</td>\n",
       "      <td>-0.191</td>\n",
       "      <td>0.188</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.002</td>\n",
       "      <td>6852.0</td>\n",
       "      <td>3479.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.403</td>\n",
       "      <td>0.072</td>\n",
       "      <td>1.271</td>\n",
       "      <td>1.536</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4684.0</td>\n",
       "      <td>2804.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  \\\n",
       "beta_parent       0.824  0.044   0.747    0.915      0.001    0.000    4643.0   \n",
       "beta_grandparent -0.656  0.111  -0.858   -0.443      0.002    0.001    3933.0   \n",
       "alpha            -0.000  0.101  -0.191    0.188      0.001    0.002    6852.0   \n",
       "sigma             1.403  0.072   1.271    1.536      0.001    0.001    4684.0   \n",
       "\n",
       "                  ess_tail  r_hat  \n",
       "beta_parent         3276.0    1.0  \n",
       "beta_grandparent    3153.0    1.0  \n",
       "alpha               3479.0    1.0  \n",
       "sigma               2804.0    1.0  "
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m9_trace, var_names=['beta_parent','beta_grandparent','alpha','sigma'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "34f3cac5-280f-45c8-a161-e37fc769eb7c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hy/9qct0wrx5pl34rzz153ys6sw0000gp/T/ipykernel_9512/2136256069.py:1: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [alpha, beta_u, beta_grandparent, beta_parent, sigma]\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='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04<00:00 Sampling 4 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 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.\n"
     ]
    }
   ],
   "source": [
    "popn['marriage_indicator'] = popn.married.apply(lambda x: 1 if x else 0)\n",
    "with pm.Model() as m10:\n",
    "    sigma = pm.Exponential('sigma', 1)\n",
    "    beta_parent = pm.Normal('beta_parent', 0, 1)\n",
    "    beta_grandparent = pm.Normal('beta_grandparent', 0, 1)\n",
    "    beta_u = pm.Normal('beta_u', 0, 1)\n",
    "    alpha = pm.Normal('alpha', 0, 1)\n",
    "    mu = pm.Deterministic('mu', alpha + beta_parent * triads['parents'] + beta_grandparent * triads['granparents']+ beta_u * triads['uobserved'])\n",
    "    sdr = pm.Normal('educational_achievement', mu=mu , sigma=sigma, observed=triads['childrens'])\n",
    "    m9_trace = pm.sample(return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "fee51b70-5cc0-4379-8ab7-0f272453977c",
   "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>beta_parent</th>\n",
       "      <td>0.095</td>\n",
       "      <td>0.081</td>\n",
       "      <td>-0.060</td>\n",
       "      <td>0.247</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.001</td>\n",
       "      <td>1561.0</td>\n",
       "      <td>1973.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_grandparent</th>\n",
       "      <td>-0.015</td>\n",
       "      <td>0.109</td>\n",
       "      <td>-0.211</td>\n",
       "      <td>0.194</td>\n",
       "      <td>0.002</td>\n",
       "      <td>0.002</td>\n",
       "      <td>1969.0</td>\n",
       "      <td>1938.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_u</th>\n",
       "      <td>1.839</td>\n",
       "      <td>0.184</td>\n",
       "      <td>1.490</td>\n",
       "      <td>2.174</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1642.0</td>\n",
       "      <td>2152.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.027</td>\n",
       "      <td>0.081</td>\n",
       "      <td>-0.123</td>\n",
       "      <td>0.174</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>4112.0</td>\n",
       "      <td>2300.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.130</td>\n",
       "      <td>0.056</td>\n",
       "      <td>1.026</td>\n",
       "      <td>1.234</td>\n",
       "      <td>0.001</td>\n",
       "      <td>0.001</td>\n",
       "      <td>3601.0</td>\n",
       "      <td>2653.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  \\\n",
       "beta_parent       0.095  0.081  -0.060    0.247      0.002    0.001    1561.0   \n",
       "beta_grandparent -0.015  0.109  -0.211    0.194      0.002    0.002    1969.0   \n",
       "beta_u            1.839  0.184   1.490    2.174      0.005    0.003    1642.0   \n",
       "alpha             0.027  0.081  -0.123    0.174      0.001    0.001    4112.0   \n",
       "sigma             1.130  0.056   1.026    1.234      0.001    0.001    3601.0   \n",
       "\n",
       "                  ess_tail  r_hat  \n",
       "beta_parent         1973.0    1.0  \n",
       "beta_grandparent    1938.0    1.0  \n",
       "beta_u              2152.0    1.0  \n",
       "alpha               2300.0    1.0  \n",
       "sigma               2653.0    1.0  "
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(m9_trace, var_names=['beta_parent','beta_grandparent', 'beta_u', 'alpha','sigma'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15c64a75-d859-432d-a311-481aff010a94",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}