statrethink_numpyro_w08.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1YI5Zm7VDTH5"
   },
   "source": [
    "# Rethinking Statistics course in NumPyro - Week 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XQXYsrXJDTID"
   },
   "source": [
    "Lecture 15: Multilevel Models\n",
    "\n",
    "- [Video](https://www.youtube.com/watch?v=AALYPv5xSos)\n",
    "- [Slides](https://speakerdeck.com/rmcelreath/l15-statistical-rethinking-winter-2019)\n",
    "\n",
    "Lecture 16: Multilevel Models 2\n",
    "\n",
    "- [Video](https://www.youtube.com/watch?v=ZG3Oe35R5sY)\n",
    "- [Slides](https://speakerdeck.com/rmcelreath/l16-statistical-rethinking-winter-2019)\n",
    "\n",
    "[Proposed problems](https://github.com/gbosquechacon/statrethinking_winter2019/blob/master/homework/week08.pdf) and [solutions in R](https://github.com/gbosquechacon/statrethinking_winter2019/blob/master/homework/week08_solutions.pdf) for the exercises of the week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "HNiX_SQ3DXpP"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "from jax import random\n",
    "from jax.scipy.special import expit as logistic\n",
    "import numpyro\n",
    "from numpyro.infer import NUTS, MCMC, Predictive\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.distributions.transforms import OrderedTransform\n",
    "\n",
    "import seaborn as sns\n",
    "import arviz as az"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5159,
     "status": "ok",
     "timestamp": 1609361667288,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "jOlA5w2QD9Xg",
    "outputId": "20208c68-8e97-4d89-c918-58fccbe8e55a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Last updated: Fri Apr 02 2021\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.9.1\n",
      "IPython version      : 7.21.0\n",
      "\n",
      "arviz  : 0.11.2\n",
      "pandas : 1.2.3\n",
      "jax    : 0.2.10\n",
      "numpyro: 0.6.0\n",
      "seaborn: 0.11.1\n",
      "\n",
      "Watermark: 2.2.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('whitegrid')\n",
    "rng_key = random.PRNGKey(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vJkV1epwDTIF"
   },
   "source": [
    "## Exercise 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Y7mqJNA1DTIF"
   },
   "source": [
    "> Revisit the Reed frog survival data, `reedfrogs`, and add the predation and size treatment variables to the varying intercepts model. Consider models with either predictor alone, both predictors, as well as a model including their interaction. What do you infer about the causal influence of these predictor variables? Also focus on the inferred variation across tanks (the σ across tanks). Explain why it changes as it does across models with different predictors included."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e4_mQ2iwDTIH"
   },
   "source": [
    "Let's get the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 142
    },
    "executionInfo": {
     "elapsed": 843,
     "status": "ok",
     "timestamp": 1609361674428,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "T0JmJamDDTII",
    "outputId": "99b55e95-0e0c-471b-b70f-7c0a20078831"
   },
   "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>density</th>\n",
       "      <th>pred</th>\n",
       "      <th>size</th>\n",
       "      <th>surv</th>\n",
       "      <th>propsurv</th>\n",
       "      <th>tank</th>\n",
       "      <th>volume</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>big</td>\n",
       "      <td>9</td>\n",
       "      <td>0.9</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>big</td>\n",
       "      <td>10</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>big</td>\n",
       "      <td>7</td>\n",
       "      <td>0.7</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   density  pred size  surv  propsurv  tank  volume\n",
       "0       10     0  big     9       0.9     0       1\n",
       "1       10     0  big    10       1.0     1       1\n",
       "2       10     0  big     7       0.7     2       1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = pd.read_csv('./dat/reedfrogs.csv', header=0, sep=';')\n",
    "d['tank'] = d.index\n",
    "d['pred'] = pd.factorize(d['pred'])[0]\n",
    "d['volume'] = d['size'].replace({'small':0, 'big':1})\n",
    "d.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4XMOL0jMDTIK"
   },
   "source": [
    "Now to define a series of models. The first is just the varying intercepts model from the text."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Drp3P2sPJR47"
   },
   "source": [
    "Model 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "m0r8LeVQJCdE"
   },
   "outputs": [],
   "source": [
    "def model(tank, density, surv=None):\n",
    "    # priors\n",
    "    a_bar = numpyro.sample('a_bar', dist.Normal(0, 1.5))\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(a_bar,sigma), sample_shape=tank.shape)\n",
    "    # likelihood\n",
    "    logits = a[tank]\n",
    "    numpyro.sample('surv_hat', dist.Binomial(density, logits=logits), obs=surv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 15586,
     "status": "ok",
     "timestamp": 1609364410322,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "gIVH800SLAEP",
    "outputId": "bb6262b8-cc19-474a-a417-8e5568db0dd7"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:06<00:00, 163.70it/s, 7 steps of size 4.66e-01. acc. prob=0.87]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 523.87it/s, 7 steps of size 4.55e-01. acc. prob=0.87]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 526.52it/s, 7 steps of size 4.68e-01. acc. prob=0.87]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 492.06it/s, 7 steps of size 4.53e-01. acc. prob=0.88]\n"
     ]
    }
   ],
   "source": [
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "dat = {k:v.to_numpy() for k,v in d[['tank', 'density', 'surv']].items()}\n",
    "mcmc.run(rng_key, **dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "id": "d7fhmoWIfRLQ"
   },
   "outputs": [],
   "source": [
    "samples_11 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bxoOaklCeFIG"
   },
   "source": [
    "Model 2 (pred):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "9dgzr_JceHeL"
   },
   "outputs": [],
   "source": [
    "def model(tank, density, pred, surv=None):\n",
    "    # priors\n",
    "    a_bar = numpyro.sample('a_bar', dist.Normal(0, 1.5))\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(a_bar,sigma), sample_shape=tank.shape)\n",
    "    bp = numpyro.sample('bp', dist.Normal(-0.5,1))\n",
    "    # likelihood\n",
    "    logits = a[tank] + bp*pred\n",
    "    numpyro.sample('surv_hat', dist.Binomial(density, logits=logits), obs=surv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 31572,
     "status": "ok",
     "timestamp": 1609364427006,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "j1hkuGEpexVL",
    "outputId": "c472fb8a-3e10-4b49-b947-9d558711c7f1"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:06<00:00, 150.98it/s, 15 steps of size 3.38e-01. acc. prob=0.88]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 508.80it/s, 15 steps of size 2.98e-01. acc. prob=0.90]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 476.86it/s, 15 steps of size 3.17e-01. acc. prob=0.89]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 502.68it/s, 15 steps of size 3.02e-01. acc. prob=0.90]\n"
     ]
    }
   ],
   "source": [
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "dat = {k:v.to_numpy() for k,v in d[['tank', 'density', 'pred', 'surv']].items()}\n",
    "mcmc.run(rng_key, **dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "x8yH16vZfTsH"
   },
   "outputs": [],
   "source": [
    "samples_12 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YvEqo4ase-dU"
   },
   "source": [
    "Model 3 (size):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "hVc2skh9fAyL"
   },
   "outputs": [],
   "source": [
    "def model(tank, density, volume, surv=None):\n",
    "    # priors\n",
    "    a_bar = numpyro.sample('a_bar', dist.Normal(0, 1.5))\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(a_bar,sigma), sample_shape=tank.shape)\n",
    "    s = numpyro.sample('s', dist.Normal(0,0.5), sample_shape=d.volume.unique().shape)\n",
    "    # likelihood\n",
    "    logits = a[tank] + s[volume]\n",
    "    numpyro.sample('surv_hat', dist.Binomial(density, logits=logits), obs=surv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 47599,
     "status": "ok",
     "timestamp": 1609364443698,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "XggMJ9AGfA6v",
    "outputId": "96ea18d6-fdc3-4034-9a09-177591fcd8b5"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:06<00:00, 151.05it/s, 15 steps of size 2.46e-01. acc. prob=0.89]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 508.20it/s, 15 steps of size 2.28e-01. acc. prob=0.91]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 489.07it/s, 15 steps of size 2.21e-01. acc. prob=0.91]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 502.49it/s, 15 steps of size 2.36e-01. acc. prob=0.91]\n"
     ]
    }
   ],
   "source": [
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "dat = {k:v.to_numpy() for k,v in d[['tank', 'density', 'volume', 'surv']].items()}\n",
    "mcmc.run(rng_key, **dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "RqrL2ZQJqURk"
   },
   "outputs": [],
   "source": [
    "samples_13 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sqoG1za3DTIO"
   },
   "source": [
    "Model 4 (pred + size):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "id": "2oDnyX6hqaZb"
   },
   "outputs": [],
   "source": [
    "def model(tank, density, pred, volume, surv=None):\n",
    "    # priors\n",
    "    a_bar = numpyro.sample('a_bar', dist.Normal(0, 1.5))\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(a_bar,sigma), sample_shape=tank.shape)\n",
    "    bp = numpyro.sample('bp', dist.Normal(-0.5,1))\n",
    "    s = numpyro.sample('s', dist.Normal(0,0.5), sample_shape=d.volume.unique().shape)\n",
    "    # likelihood\n",
    "    logits = a[tank] + bp*pred + s[volume]\n",
    "    numpyro.sample('surv_hat', dist.Binomial(density, logits=logits), obs=surv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 64999,
     "status": "ok",
     "timestamp": 1609364462062,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "YlICg_fIqtBF",
    "outputId": "3a36db76-87e3-47eb-a60d-9670adc4370a"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:07<00:00, 138.78it/s, 15 steps of size 1.88e-01. acc. prob=0.92]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 339.34it/s, 15 steps of size 2.36e-01. acc. prob=0.87]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 449.63it/s, 15 steps of size 2.25e-01. acc. prob=0.89]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 498.12it/s, 15 steps of size 1.79e-01. acc. prob=0.93]\n"
     ]
    }
   ],
   "source": [
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "dat = {k:v.to_numpy() for k,v in d[['tank', 'density', 'pred', 'volume', 'surv']].items()}\n",
    "mcmc.run(rng_key, **dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "id": "X7N47Ammq4gS"
   },
   "outputs": [],
   "source": [
    "samples_14 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QUF74m4UDTIP"
   },
   "source": [
    "Model 5 (pred + size + interaction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "id": "unjwwWC-NHWX"
   },
   "outputs": [],
   "source": [
    "def model(tank, density, pred, volume, surv=None):\n",
    "    # priors\n",
    "    a_bar = numpyro.sample('a_bar', dist.Normal(0, 1.5))\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(a_bar,sigma), sample_shape=tank.shape)\n",
    "    bp = numpyro.sample('bp', dist.Normal(-0.5,1), sample_shape=d.volume.unique().shape)\n",
    "    s = numpyro.sample('s', dist.Normal(0,0.5), sample_shape=d.volume.unique().shape)\n",
    "    # likelihood\n",
    "    logits = a[tank] + bp[volume]*pred + s[volume]\n",
    "    numpyro.sample('surv_hat', dist.Binomial(density, logits=logits), obs=surv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 81098,
     "status": "ok",
     "timestamp": 1609364480566,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "rv_mkoHFOd8O",
    "outputId": "3284772a-740b-4de3-d86d-d0e127550b2d"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:07<00:00, 140.05it/s, 15 steps of size 1.48e-01. acc. prob=0.93]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 489.48it/s, 15 steps of size 2.01e-01. acc. prob=0.90]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 445.81it/s, 15 steps of size 1.74e-01. acc. prob=0.91]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 482.28it/s, 15 steps of size 2.07e-01. acc. prob=0.90]\n"
     ]
    }
   ],
   "source": [
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "dat = {k:v.to_numpy() for k,v in d[['tank', 'density', 'pred', 'volume', 'surv']].items()}\n",
    "mcmc.run(rng_key, **dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "id": "EzBc3XmEOfzB"
   },
   "outputs": [],
   "source": [
    "samples_15 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XoewvGYeDTIP"
   },
   "source": [
    "I coded the interaction model (`samples_5`) using a non-centered parameterization. The interaction itself is done by creating a `bp` parameter for each size value. In this way, the effect of `pred` depends upon `volume`. First let's consider the WAIC scores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:145: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:655: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rank</th>\n",
       "      <th>loo</th>\n",
       "      <th>p_loo</th>\n",
       "      <th>d_loo</th>\n",
       "      <th>weight</th>\n",
       "      <th>se</th>\n",
       "      <th>dse</th>\n",
       "      <th>warning</th>\n",
       "      <th>loo_scale</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>model_12</th>\n",
       "      <td>0</td>\n",
       "      <td>208.715087</td>\n",
       "      <td>24.077409</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>6.175129e-01</td>\n",
       "      <td>10.328203</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model_15</th>\n",
       "      <td>1</td>\n",
       "      <td>210.272722</td>\n",
       "      <td>24.416782</td>\n",
       "      <td>1.557635</td>\n",
       "      <td>3.259616e-01</td>\n",
       "      <td>10.517428</td>\n",
       "      <td>4.553388</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model_14</th>\n",
       "      <td>2</td>\n",
       "      <td>211.286669</td>\n",
       "      <td>24.792323</td>\n",
       "      <td>2.571582</td>\n",
       "      <td>2.094098e-15</td>\n",
       "      <td>10.425331</td>\n",
       "      <td>2.216642</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model_13</th>\n",
       "      <td>3</td>\n",
       "      <td>219.070288</td>\n",
       "      <td>30.446760</td>\n",
       "      <td>10.355201</td>\n",
       "      <td>1.221035e-15</td>\n",
       "      <td>7.993552</td>\n",
       "      <td>6.939584</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model_11</th>\n",
       "      <td>4</td>\n",
       "      <td>220.604288</td>\n",
       "      <td>31.012833</td>\n",
       "      <td>11.889201</td>\n",
       "      <td>5.652543e-02</td>\n",
       "      <td>8.633994</td>\n",
       "      <td>7.292377</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          rank         loo      p_loo      d_loo        weight         se  \\\n",
       "model_12     0  208.715087  24.077409   0.000000  6.175129e-01  10.328203   \n",
       "model_15     1  210.272722  24.416782   1.557635  3.259616e-01  10.517428   \n",
       "model_14     2  211.286669  24.792323   2.571582  2.094098e-15  10.425331   \n",
       "model_13     3  219.070288  30.446760  10.355201  1.221035e-15   7.993552   \n",
       "model_11     4  220.604288  31.012833  11.889201  5.652543e-02   8.633994   \n",
       "\n",
       "               dse  warning loo_scale  \n",
       "model_12  0.000000     True  deviance  \n",
       "model_15  4.553388     True  deviance  \n",
       "model_14  2.216642     True  deviance  \n",
       "model_13  6.939584     True  deviance  \n",
       "model_11  7.292377     True  deviance  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.compare({'model_11': samples_11, \n",
    "            'model_12': samples_12, \n",
    "            'model_13': samples_13, \n",
    "            'model_14': samples_14,\n",
    "            'model_15': samples_15}\n",
    "           , scale='deviance')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Vg-9gtDZDTIQ"
   },
   "source": [
    "These models are really very similar in expected out-of-sample accuracy. The tank variation is huge. But take a look at the posterior distributions for predation and size. You'll see that predation does seem to matter, as you'd expect. Size matters a lot less. So while predation doesn't explain much of the total variation, there is plenty of evidence that it is a real effect. Remember: We don't select a model using WAIC\n",
    "(or LOO). A predictor can make little difference in total accuracy but still be a real causal effect.\n",
    "\n",
    "Let's look at all the sigma posterior distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 241
    },
    "executionInfo": {
     "elapsed": 1359,
     "status": "ok",
     "timestamp": 1609365078554,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "PN0RrIp0lxWW",
    "outputId": "f6db69b3-56ee-4d2d-c492-1fa35a4548db"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApAAAADgCAYAAAC0CFa4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAw2UlEQVR4nO3de1hU9b4/8DcwXBTIUUsUsAnYMJAkM4l22CEHfULSSglvbPWUtxTQNNHjT8twREk02eZl4wY3+7SDFEto1KfYlubOvGBHA+QYIAJeQNExYAuC4gzz+8McHPEyy+Yi8H49D88z67vWfNdnfZqGN2vNGq20Wq0WREREREQGsrZ0AURERETUsTBAEhEREZEgDJBEREREJAgDJBEREREJwgBJRERERIIwQBIRERGRIAyQRER3uX79OhISEvDHP/4RcrkcM2bMQElJyQO3r62tRVBQEDZt2vTIuY8fP47x48cjICAAI0aMwM6dO/XWl5WVITIyUrffK1eu6K1fs2YNPvjgA4OOY/jw4UhISLjvupycHEilUtTW1uot3/nx9fWFTCbD6NGjkZaWhlu3buk9f9OmTZDL5QbVQUSdEwMkEdFd5s2bh6+++gozZszApk2b8PTTT2PSpEmoqKi47/aJiYm6IPYw5eXlmDlzJtzd3bFp0yaEhobigw8+wD//+U/dNh9++CH69euHlJQUXL9+HWvXrtWtu3z5Mnbu3Im5c+f+/oN8gL/97W/YsWMHtm/fjo0bNyI4OBgbN27E3LlzodFoTLZfIup4RJYugIjoSfF///d/OHToEFasWIGoqCgAQHBwMM6dO4cNGzZgw4YNett///33OHToEOzt7R85d1paGtzc3PDnP/8ZVlZWCAkJQV1dHf7yl7/g1VdfBQAUFxdj06ZNCAoKwtmzZ5GZmal7fkpKCiIiItCvXz8jHrG+AQMGoFevXrrlkJAQeHp64oMPPsBXX32FcePGmWzfRNSx8AwkEdFvzp49C+B2aLybXC7HoUOH9MYaGhqgUCiwZMkS2NnZPXLuI0eOIDQ0FFZWVrqxV155BadPn8bly5cBAG5ubjh69CgaGxvx008/wc3NDQBw4cIF5ObmIjo6+vcc3mMZO3Ys3Nzc2l1uJ6KujQGSiOg3ffv2BQBcunRJb7y6uhqNjY2or6/Xja1ZswZ/+MMf8Oabbz5y3qamJly5cgUSiURvvH///gDaguvixYuRlZWFQYMG4cSJE4iLiwMAbNy4EZMmTULv3r0FHY9Wq4VarW7309raavAcVlZWeOmll1BUVNTus5BE1HXxEjYR0W8GDhyI5557DitWrMDq1ashkUjwzTff4IcffgAANDc3QywW4+jRo/j666+xe/dug+ZtbGwEADg6OuqN31m+sz40NBSHDh1CTU0N+vfvDzs7O5SVleHHH3/Ed999h+zsbPzjH//AU089hffffx/PP//8Q/e7bds2bNu2TVAP7qdXr15Qq9X497//jaeffvp3z0dEHR8DJBHRb+zs7LB582YsXLhQ93k/uVyOmTNnYvPmzXBwcEBzczM+/PBDvPvuu7oziI+i1WoBQO/y9d2srdsuBjk6OsLLy0u3/Mknn2DGjBmorq7GqlWrkJaWhl9++QXR0dHYt2/fQy+fjxw5EjNmzGg3/q9//QubN282qHYiovthgCQiuou3tzd2796NS5cuQa1Wo3///ti8eTOsra3h7OyMtWvXwtnZGVOmTIFardY9r7W1FWq1GiJR+7dVJycnALe/Iuhud5adnZ3vW8vJkydx8uRJJCcnIzU1FYGBgRg8eDAGDRqETz75BAUFBRgyZMgDj6VXr1544YUX2o2XlZU9uhF3uXz5Muzs7CAWiwU9j4g6LwZIIqLfNDc3Y+/evQgKCtK727m0tBTe3t4QiUTYt28fqqur2wWzlJQUpKSkoLS0tN28jo6OeOaZZ3DhwgW98TvLHh4e961n/fr1iI6OhoODA3799Vf06NEDwO0zlk5OTrh69ervOl5DtLa24n//938hk8nuG46JqGviTTRERL8RiURQKBT45ptvdGMXLlzADz/8gGHDhgEAtmzZgp07d+r9dO/eHRMmTHjoncpBQUE4cOCA3vcp7tu3Dz4+Pve9OebYsWO4cOECJkyYAADo3bu3LjC2tLSgvr5e7yt3TGXXrl2oqanB+PHjTb4vIuo4+OckEdFvbG1tMW7cOPz1r39Fr1694OTkhHXr1qFXr16YOnUqAEAqlbZ7no2NDfr06aN3VvLMmTNoaWnR3egyY8YMjBs3DvPnz8f48eNx5MgR7N69u913S96xfv16zJ07F7a2tgCAYcOGIS0tDV999RVOnz6Np556CjKZzKjHf+rUKTg7O0Or1eLatWs4duwYPvvsMwwfPhxvvPGGUfdFRB0bAyQR0V0WLVoEKysrfPzxx7h58yb+4z/+A4sXL0bPnj0FzbNixQpUV1fj+++/BwD4+vpiy5YtWLduHebOnQtXV1esXr1a9yXid/v+++/R2NiI0aNH68YGDhyIuLg4rF27Fj169MAnn3wCBweH33ew95g5c6bucbdu3eDp6YlFixZhypQpD7wBiIi6JivtndsDiYiIiIgMwM9AEhEREZEgDJBEREREJAgDJBEREREJwgBJRERERIIwQBIRERGRIPwaHwHKysqgVqtx8+ZN2NvbW7qcJwJ70Ya9aMNe6GM/2rAXbdiLNuyFvielHyKRCN7e3vddx6/xeQzFxcXw8/OzdBlPBPaiDXvRhr3Qx360YS/asBdt2At9HaEfvIRNRERERIIwQBIRERGRIAyQRERERCQIAyQRERERCcIASURERESCMEASERERkSAMkEREREQkCAMkEREREQnCAElEREREgjBAEhEREZEgDJBEREREJAgDJBEREREJwgBJRERERIIwQBIRERGRIAyQRERERCQIAyQRERERCWLxAHnkyBFIpVKDts3JyUFISIjgfRw/fhyhoaH3Xbd792786U9/EjwnERERUVdl8QBpaqWlpZg/fz5aW1vbrcvLy0N8fLwFqiIiIiLquDp1gMzKykJUVBR69+7dbt3mzZvxzjvvoH///haojIiIuqorDTdw4lwdapvVli6F6LEZFCCrqqoglUqxf/9+DB8+HHK5HElJSSgtLUVkZCRkMhmio6PR1NQE4Pal5lGjRmHgwIGIjIzEsWPHdHM1NjYiLi4Ocrkc4eHhOHXqlN6+ampqEBsbC5lMhtDQUKxbtw4tLS2PdXAHDx7EmjVrMHXq1HbrDh8+jPT0dIwYMeKx5iYiIhJKmV+NEesPImHPKcTuroYyv9rSJRE9FpGQjbdu3YqUlBSUlpZi8eLFOHDgABQKBUQiEaKjo5GdnQ1HR0ckJCRAoVAgICAAOTk5mDVrFnJzc+Hq6orly5ejoqICmZmZqK2txZIlS3Tza7VazJkzBz4+PsjOzkZdXR0UCgXUarXedoZKSUkBcDvQ3mv79u0AoBduiYieJBNTj1q6BKNrampC94P1Rp2zqq4J1fU3jDqnqTiIrLHn3WB4uzij7HID3th0CO/tKLB0WY/FTewA957djTKXKV4XHdWO2UGWLsEggi5hx8TEwNfXF2PGjIFYLMZrr72GoKAgDB48GEOGDEFFRQUyMjIwefJkREREwMPDAwsXLoSvry8yMjLQ0NCA3NxcLF26FAMGDMDQoUMRExOjmz8vLw9VVVVYtWoVvLy8EBgYiPj4eGRmZkKt5ql+IiLq2PqJu8HbxRkA4O3ijH7ibhauiOjxCDoD6e7urntsb28PV1dX3bKDgwNaWlpQXl6uFwoBQCaToaKiApWVldBoNPD19dWt8/f31z0uLy/HtWvXEBgYqBvTarW4desWLl68KKRUIqIOr6OciRCiuLgYfn5+Rp1z/XensWF/mVHnNJVL9c0ou9ygOwN5qb7Z0iU9tnGD+mNBmI9R5jLF64JMS1CAFIn0N7e2bn8C08HBod2YRqOBRqPRLWu12vvOqVarIZFIkJqa2m6Ovn37CimViIi6iAVhPkYLMqamzK/G+NSjkPTqjrNXG5E0diAi5G6WLotIMKPfhe3p6YnCwkK9scLCQnh4eMDT0xO2trYoKirSrSsuLtY99vDwQE1NDcRiMSQSCSQSCVQqFZKTk/VCJxERUUcUIXfDtwtCEP/GAKSMdmN4pA7L6AFy2rRp2LZtG5RKJSorK5GcnIySkhJMmDABTk5OGD16NBITE1FQUIC8vDzdjS4AEBwcDHd3dyxatAglJSXIz8/HsmXLYG1tDXt7e2OXSkREZHZ9nB0wSNITvboJughI9EQx+qs3PDwcKpUKGzduhEqlgp+fH9LT0+Ht7Q0AiI+Px8qVKzF9+nSIxWJMmTIFa9euBQDY2Nhgy5YtSExMRFRUFOzt7REWFvZYd2ATERERkWlYaXltWDB+2LcNe9GGvWjDXuhjP9qwF23Yizbshb6O0I9O/S/REBEREZHxdcgPYOzdu/ehl7WlUimysrLMWBERERFR19EhA2RwcDCUSuUD19vZ2ZmvGCIiIqIupkMGSEdHRzg6Olq6DCIiIqIuiZ+BJCIiIiJBGCCJiIiISBAGSCIiIiIShAGSiIiIiARhgCQiIiIiQRggiYiIiEgQBkgiIiIiEoQBkoiIiIgEYYAkIiIiIkEYIImIiIhIEAZIIiIiIhKEAZKIiIiIBGGAJCIiIiJBGCCJiIiISBAGSCIiIiIShAGSiIiIiARhgCQiIiIiQRggiYiIiEgQiwbII0eOQCqVGrRtTk4OQkJCBO/j+PHjCA0NbTf+2WefYcSIEZDL5Zg4cSJ+/vlnwXMTERERdUWd+gxkaWkp5s+fj9bWVr3xXbt2YcOGDYiLi4NSqcSQIUMwc+ZMXL582UKVEhEREXUcnTZAZmVlISoqCr1792637quvvsKkSZPw6quvQiKRYOHChXjmmWfw/fffW6BS6qyuNNzAiXN1uNJww9KlEBERGdUjA2RVVRWkUin279+P4cOHQy6XIykpCaWlpYiMjIRMJkN0dDSampoA3L7UPGrUKAwcOBCRkZE4duyYbq7GxkbExcVBLpcjPDwcp06d0ttXTU0NYmNjIZPJEBoainXr1qGlpeWxDuzgwYNYs2YNpk6d2m7de++9hwkTJrQbb2hoeKx9Ed1LmV+NEesPImHPKYxYfxDK/GpLl0RERGQ0IkM33Lp1K1JSUlBaWorFixfjwIEDUCgUEIlEiI6ORnZ2NhwdHZGQkACFQoGAgADk5ORg1qxZyM3NhaurK5YvX46KigpkZmaitrYWS5Ys0c2v1WoxZ84c+Pj4IDs7G3V1dVAoFFCr1XrbGSolJQXA7UB7L5lMprd88OBBnD17FkFBQYL301VMTD163/GmpiZ0P1hv3mIEqqprQnW9uc4CVgAAHETW2PNuMLxdnFF2uQFvbDqE93YUmKmG+3MTO8C9Z3ez7MtSr4sds/n/MBGRORh8CTsmJga+vr4YM2YMxGIxXnvtNQQFBWHw4MEYMmQIKioqkJGRgcmTJyMiIgIeHh5YuHAhfH19kZGRgYaGBuTm5mLp0qUYMGAAhg4dipiYGN38eXl5qKqqwqpVq+Dl5YXAwEDEx8cjMzMTarXaJAcPAJWVlfh//+//ISIiAi+88ILJ9kNdSz9xN3i7OAMAvF2c0U/czcIVERERGY/BZyDd3d11j+3t7eHq6qpbdnBwQEtLC8rLy/VCIXD7bF9FRQUqKyuh0Wjg6+urW+fv7697XF5ejmvXriEwMFA3ptVqcevWLVy8eFHYURno9OnTmD59Ojw9PZGQkGCSfXQWDzqzU1xcDD8/PzNXI8z6705jw/4ys+7zUn0zyi436M5AXqpvNuv+72fcoP5YEOZjln11hNcFERE9PoMDpEikv6m1dfuTlw4ODu3GNBoNNBqNblmr1d53TrVaDYlEgtTU1HZz9O3b19AyDVZUVISZM2fCx8cHf/3rX2Fvb2/0fdCTYUGYj1mC092hSZlfjfGpRyHp1R3napuQNHYgIuRuJq+BiIjIHIx6F7anpycKCwv1xgoLC+Hh4QFPT0/Y2tqiqKhIt664uFj32MPDAzU1NRCLxZBIJJBIJFCpVEhOTtYLncZw4cIFzJw5E88//zzS0tLg6Oho1PmJIuRu+HZBCOLfGIBvF4QwPBIRUadi1AA5bdo0bNu2DUqlEpWVlUhOTkZJSQkmTJgAJycnjB49GomJiSgoKEBeXp7uRhcACA4Ohru7OxYtWoSSkhLk5+dj2bJlsLa2NvrZwRUrVqB79+5ISEhAY2MjVCoVVCoVrl+/btT9UNfWx9kBgyQ90ce5/Zl5IiKijszgS9iGCA8Ph0qlwsaNG6FSqeDn54f09HR4e3sDAOLj47Fy5UpMnz4dYrEYU6ZMwdq1awEANjY22LJlCxITExEVFQV7e3uEhYU91h3YD9PY2Igff/wRAPDKK6/orYuOjsaCBQuMuj8iIiKizsZKa+zrw10AbxBow160YS/asBf62I827EUb9qINe6GvI/Sj0/5LNERERERkGka9hG0Oe/fufehlbalUiqysLDNWRERERNS1dLgAGRwcDKVS+cD1dnZ25iuGiIiIqAvqcAHS0dGRX7tDREREZEH8DCQRERERCcIASURERESCMEASERERkSAMkEREREQkCAMkEREREQnCAElEREREgjBAEhEREZEgDJBEREREJAgDJBEREREJwgBJRERERIIwQBIRERGRIAyQRERERCQIAyQRERERCcIASURERESCMEASERERkSAMkEREREQkCAMkEREREQnCAElEREREglg0QB45cgRSqdSgbXNychASEiJ4H8ePH0doaGi78bS0NPznf/4nAgICMGPGDJw9e1bw3ERERERdUac+A1laWor58+ejtbVVb/zLL79Eeno6Vq5ciV27dqFHjx6YPXt2u+2IiIiIqL1OGyCzsrIQFRWF3r17t1t3/fp1LF68GCEhIXjuuecwa9YsnD17FiqVygKVEhFRV1TbrMaJc3W40nDD0qUQCfbIAFlVVQWpVIr9+/dj+PDhkMvlSEpKQmlpKSIjIyGTyRAdHY2mpiYAty81jxo1CgMHDkRkZCSOHTumm6uxsRFxcXGQy+UIDw/HqVOn9PZVU1OD2NhYyGQyhIaGYt26dWhpaXmsAzt48CDWrFmDqVOntls3depUjB07FgDQ0NCAbdu2wdvbG88888xj7YuIiEgIZX41YndXI2HPKYxYfxDK/GpLl0QkiMjQDbdu3YqUlBSUlpZi8eLFOHDgABQKBUQiEaKjo5GdnQ1HR0ckJCRAoVAgICAAOTk5mDVrFnJzc+Hq6orly5ejoqICmZmZqK2txZIlS3Tza7VazJkzBz4+PsjOzkZdXR0UCgXUarXedoZKSUkBcDvQPsiOHTsQHx8POzs7pKenw9q6056QJaIOYmLqUUuXYFJNTU3ofrDeJHNX1TWhur5jnM1zEFljz7vB8HZxRtnlBryx6RDe21Fg6bIei5vYAe49u/+uOUz5uuiIHtWPHbODzFfMAxicmGJiYuDr64sxY8ZALBbjtddeQ1BQEAYPHowhQ4agoqICGRkZmDx5MiIiIuDh4YGFCxfC19cXGRkZaGhoQG5uLpYuXYoBAwZg6NChiImJ0c2fl5eHqqoqrFq1Cl5eXggMDER8fDwyMzOhVqtNcvDBwcHIycnB2LFjERsbiwsXLphkP0RERHfrJ+4GbxdnAIC3izP6ibtZuCIiYQw+A+nu7q57bG9vD1dXV92yg4MDWlpaUF5erhcKAUAmk6GiogKVlZXQaDTw9fXVrfP399c9Li8vx7Vr1xAYGKgb02q1uHXrFi5evCjsqAzk5uYGNzc3PP/88zh27BiUSiXeffddk+yLiMgQT8KZBVMqLi6Gn5+fSeZe/91pbNhfZpK5je1SfTPKLjfozkBeqm+2dEmPbdyg/lgQ5vO75jDl66Ij6gj9MDhAikT6m97vcq+Dg0O7MY1GA41Go1vWarX3nVOtVkMikSA1NbXdHH379jW0TIMcPnwY7u7ukEgkAAArKyt4eXmhrq7OqPshIiLzWRDm87uDjLko86sxbsthPPe0E87VNiFp7EBEyN0sXRaRwYz6oT9PT08UFhbqjRUWFsLDwwOenp6wtbVFUVGRbl1xcbHusYeHB2pqaiAWiyGRSCCRSKBSqZCcnKwXOo1h8+bNyMjI0C2r1WqUlJTAy8vLqPshIiK6nwi5G1JGuyH+jQH4dkEIwyN1OEYNkNOmTcO2bdugVCpRWVmJ5ORklJSUYMKECXBycsLo0aORmJiIgoIC5OXl6W50AW5/HtHd3R2LFi1CSUkJ8vPzsWzZMlhbW8Pe3t6YZeK//uu/sGPHDnz99dcoLy9HfHw8bt26hYiICKPuh4iI6EF6dRNhkKQn+ji3v3pH9KQz+BK2IcLDw6FSqbBx40aoVCr4+fkhPT0d3t7eAID4+HisXLkS06dPh1gsxpQpU7B27VoAgI2NDbZs2YLExERERUXB3t4eYWFhj3UH9qOMGjUK165dwyeffILLly9DLpfj73//OxwdHY2+LyIiIqLOxkpr7OvDXUBH+HCrubAXbdiLNuyFPvajDXvRhr1ow17o6wj94BcfEhEREZEgRr2EbQ579+596GVtqVSKrKwsM1ZERERE1LV0uAAZHBwMpVL5wPV2dnbmK4aIiIioC+pwAdLR0ZE3uxARERFZED8DSURERESCMEASERERkSAMkEREREQkCAMkEREREQnCAElEREREgjBAEhEREZEgDJBEREREJAgDJBEREREJwgBJRERERIIwQBIRERGRIAyQRERERCQIAyQRERERCcIASURERESCMEASERERkSAMkEREREQkCAMkEREREQnCAElEREREgjBAEhEREZEgFg2QR44cgVQqNWjbnJwchISECN7H8ePHERoa+sD1586dQ0BAAI4dOyZ4biIiIqKuqFOfgSwtLcX8+fPR2tp63/VarRbLli3DjRs3zFwZERERUcfVaQNkVlYWoqKi0Lt37wdus337dmg0GjNWRR3FlYYbOHGuDlca+McFERHRvR4ZIKuqqiCVSrF//34MHz4ccrkcSUlJKC0tRWRkJGQyGaKjo9HU1ATg9qXmUaNGYeDAgYiMjNS7NNzY2Ii4uDjI5XKEh4fj1KlTevuqqalBbGwsZDIZQkNDsW7dOrS0tDzWgR08eBBr1qzB1KlT77v+0qVL2LRpE1auXPlY81Pnpcyvxoj1B5Gw5xRGrD8IZX61pUsiIiJ6oogM3XDr1q1ISUlBaWkpFi9ejAMHDkChUEAkEiE6OhrZ2dlwdHREQkICFAoFAgICkJOTg1mzZiE3Nxeurq5Yvnw5KioqkJmZidraWixZskQ3v1arxZw5c+Dj44Ps7GzU1dVBoVBArVbrbWeolJQUALcD7f3Ex8dj6tSpkEgkgud+kkxMPWrR/Tc1NaH7wXqTzV9V14TqevOeBXQQWWPPu8HwdnFG2eUGvLHpEN7bUWDgsyuMXo+b2AHuPbsbfV5TMtbrYsfsoN9fDBERGZ3Bl7BjYmLg6+uLMWPGQCwW47XXXkNQUBAGDx6MIUOGoKKiAhkZGZg8eTIiIiLg4eGBhQsXwtfXFxkZGWhoaEBubi6WLl2KAQMGYOjQoYiJidHNn5eXh6qqKqxatQpeXl4IDAxEfHw8MjMzoVarjXrQSqUSV65cwYwZM4w6L3UO/cTd4O3iDADwdnFGP3E3C1dERET0ZDH4DKS7u7vusb29PVxdXXXLDg4OaGlpQXl5uV4oBACZTIaKigpUVlZCo9HA19dXt87f31/3uLy8HNeuXUNgYKBuTKvV4tatW7h48aKwo3qIq1evYs2aNUhLS4NIJDJ6ODU3S5+hKS4uhp+fn8nmX//daWzYX2ay+e/nUn0zyi436M5AXqpvNuv+7zVuUH8sCPOxaA1Cmfp1QURElmVwgBSJ9De1tm5/8tLBwaHdmEaj0btRRavV3ndOtVoNiUSC1NTUdnP07dvX0DIf6ccff0RdXR3eeustvfF33nkHsbGxiI6ONtq+6PdbEOZj9vCkzK/G+NSjkPTqjnO1TUgaOxARcrdHPo+hiYiIugqDA6QhPD09UVhYiBEjRujGCgsLIZPJ4OnpCVtbWxQVFWHo0KEAbv/CvcPDwwM1NTUQi8Xo0aMHgNvf4fjZZ5/h448/NlqNYWFhePHFF3XLGo0GI0eOxKpVqx7reyap84mQu+GPf+iNC7XN6N+rG/o4t//DiIiIqCszaoCcNm0alixZAm9vb91NNCUlJfjoo4/g5OSE0aNHIzExEUlJSbhx44buRhcACA4Ohru7OxYtWoSFCxeiubkZy5Ytg6+vL+zt7Y1Wo5OTE5ycnHTLdy5hu7i4QCwWG20/1LH1cXZgcCQiInoAowbI8PBwqFQqbNy4ESqVCn5+fkhPT4e3tzeA23c+r1y5EtOnT4dYLMaUKVOwdu1aAICNjQ22bNmCxMREREVFwd7eHmFhYY91BzYRERERmY6V9u4PJZJB+Fm3NuxFG/aiDXuhj/1ow160YS/asBf6OkI/Ou2/RENEREREpmHUS9jmsHfv3ode1pZKpcjKyjJjRURERERdS4cLkMHBwVAqlQ9cb2dnZ75iiIiIiLqgDhcgHR0d4ejoaOkyiIiIiLosfgaSiIiIiARhgCQiIiIiQRggiYiIiEgQBkgiIiIiEoQBkoiIiIgEYYAkIiIiIkEYIImIiIhIEAZIIiIiIhKEAZKIiIiIBGGAJCIiIiJBGCCJiIiISBAGSCIiIiIShAGSiIiIiARhgCQiIiIiQRggiYiIiEgQBkgiIiIiEoQBkoiIiIgEYYAkIiIiIkEsGiCPHDkCqVRq0LY5OTkICQkRvI/jx48jNDS03Xh4eDikUqneT3FxseD5iYiIiLoakaULMKXS0lLMnz8fNjY2euMtLS24cOECtm/fjv79++vGe/bsae4SiYiIiDqcThsgs7KysGbNGvTv3x/19fV66yoqKmBlZYUXXngBtra2limQLO5Kww1cqG1G/17d0MfZwdLlEBERdRiPvIRdVVUFqVSK/fv3Y/jw4ZDL5UhKSkJpaSkiIyMhk8kQHR2NpqYmALcvNY8aNQoDBw5EZGQkjh07ppursbERcXFxkMvlCA8Px6lTp/T2VVNTg9jYWMhkMoSGhmLdunVoaWl5rAM7ePAg1qxZg6lTp7ZbV15eDnd3d4bHLkyZX40R6w8iYc8pjFh/EMr8akuXRERE1GEYfAZy69atSElJQWlpKRYvXowDBw5AoVBAJBIhOjoa2dnZcHR0REJCAhQKBQICApCTk4NZs2YhNzcXrq6uWL58OSoqKpCZmYna2losWbJEN79Wq8WcOXPg4+OD7Oxs1NXVQaFQQK1W621nqJSUFAC3A+29zpw5AxsbG8ycORPFxcXw8PDAf//3fyMgIEDwfkxlYupRS5dgkKamJnQ/WP/Yz6+qa0J1/Q3jFWQgB5E19rwbDG8XZ5RdbsAbmw7hvR0FRpi5QvAz3MQOcO/Z3Qj7fnI8zutix+wg0xRDRERGZ/BNNDExMfD19cWYMWMgFovx2muvISgoCIMHD8aQIUNQUVGBjIwMTJ48GREREfDw8MDChQvh6+uLjIwMNDQ0IDc3F0uXLsWAAQMwdOhQxMTE6ObPy8tDVVUVVq1aBS8vLwQGBiI+Ph6ZmZlQq9VGPejy8nJcu3YNkyZNQlpaGry8vPD222+jqqrKqPuhJ1c/cTd4uzgDALxdnNFP3M3CFREREXUcBp+BdHd31z22t7eHq6urbtnBwQEtLS0oLy/XC4UAIJPJUFFRgcrKSmg0Gvj6+urW+fv76x7fCXWBgYG6Ma1Wi1u3buHixYvCjuoRkpOTcfPmTTg5OQEAFAoFfv75ZyiVSsydO9eo+3pcHeVsTHFxMfz8/B77+eu/O40N+8uMWJFhLtU3o+xyg+4M5KX6ZrPXcMe4Qf2xIMzHYvs3hd/7uiAioiebwQFSJNLf1Nq6/clLB4f2NyJoNBpoNBrdslarve+carUaEokEqamp7ebo27evoWUaxNbWVu/zj1ZWVvD09MSVK1eMuh96tAVhPhYJT8r8aoxPPQpJr+44V9uEpLEDESF3+11zMjQREVFXYdTvgfT09ERhYaHeWGFhITw8PODp6QlbW1sUFRXp1t39vYseHh6oqamBWCyGRCKBRCKBSqVCcnKyXug0hnHjxiEtLU233NraitLSUnh6ehp1P/TkipC74dsFIYh/YwC+XRDyu8MjERFRV2LUADlt2jRs27YNSqUSlZWVSE5ORklJCSZMmAAnJyeMHj0aiYmJKCgoQF5enu5GFwAIDg6Gu7s7Fi1ahJKSEuTn52PZsmWwtraGvb29MctEaGgo0tPT8cMPP6CiogIKhQL//ve/MXbsWKPuh55sfZwdMEjSk1/hQ0REJJBRvwcyPDwcKpUKGzduhEqlgp+fH9LT0+Ht7Q0AiI+Px8qVKzF9+nSIxWJMmTIFa9euBQDY2Nhgy5YtSExMRFRUFOzt7REWFvZYd2A/SkxMDFpbW7F8+XLU1tYiICAAn376KZydnY2+LyIiIqLOxkpr7OvDXQA/69aGvWjDXrRhL/SxH23YizbsRRv2Ql9H6IdF/y1sIiIiIup4Otw/Zbh3796HXtaWSqXIysoyY0VEREREXUuHC5DBwcFQKpUPXG9nZ2e+YoiIiIi6oA4XIB0dHeHo6GjpMoiIiIi6LH4GkoiIiIgE4V3YApSVlRn93+UmIiIiehKJRCLdVzHeiwGSiIiIiAThJWwiIiIiEoQBkoiIiIgEYYAkIiIiIkEYIImIiIhIEAZIIiIiIhKEAfI+Wlpa8OGHH2Lw4MF4+eWXsXXr1kc+p76+Hi+//DJycnLMUKH5COlFeXk53nrrLQQEBCA8PBx79+41Y6WmJ6QXx48fR2RkJGQyGcaMGYNDhw6ZsVLzaWlpweuvv44jR448cJvq6mpMnz4dMpkMI0eOxA8//GDGCs3HkF4cPXoUY8eOhVwuR3h4OL788kszVmhehvTj3m03bdpkhsrMz5BeXL58GbGxsZDJZAgNDcXnn39uxgrNx5BenD59GlOmTNH9f7J7924zVmh658+fR3R0NAYPHoyQkBAkJSXh5s2b9922pKQEEydOREBAACIjI3Hy5EkzV/tgDJD3sXbtWhQUFOB//ud/sGLFCmzZsgVff/31Q5/z0Ucf4erVq2aq0HwM7cX169cxbdo09O3bF7t27cLkyZOxcOFCnDlzxgJVm4ahvfj1118RHR2NV199Fbt378bIkSMxZ84cVFdXW6Bq07l58ybi4uJQVlb2wG20Wi1iY2MhFouxc+dOvPnmm5g3bx4uXLhgxkpNz5BenD17FrNnz0ZYWBiUSiXmzJmDhIQEfP/992as1DwM6cfdtmzZYvC2HY0hvWhtbUVMTAxu3ryJ7OxsLFq0CKtXr8bhw4fNWKnpGdKLlpYWREdHw9fXF7t27cI777yDpUuXorCw0IyVms6d47Ozs0NWVhbWrVuHffv2Yf369e22bWpqwsyZMxEQEICcnBwMGjQIs2fPRmNjowUqb48B8h5NTU344osvsHTpUvj7++OVV17BzJkzkZmZ+cDn/PDDDzh58iR69eplxkpNT0gvlEolRCIREhMT8dxzz+Gtt97Cyy+/jPz8fAtUbnxCevHzzz8DAGbNmoVnn30W0dHRcHBw6DRvgABw5swZTJgwAefPn3/odnl5eaisrERCQgL+8Ic/YNasWZDL5di5c6eZKjU9Q3vxzTffwM/PD9HR0ZBIJBg9ejQiIiKwZ88eM1VqHob2446SkhJ8+eWX8PT0NHFl5mdoL3788UecPXsWycnJ8PLywuuvv47IyMhO8/4JGN6LM2fOoLq6GvPmzcOzzz6LcePGQSqV4tixY2aq1LROnjyJ8+fPY/Xq1fDy8sKQIUMwf/78+74PfPPNN7C1tcWSJUvg5eWF999/H87OzsjNzbVA5e0xQN6jpKQELS0tGDRokG5s0KBBKCoqgkajabd9Y2MjFAoFVq5cCVtbW3OWanJCenHs2DEMHz5crwepqakYP3682eo1JSG9EIvFaGhoQG5uLrRaLfbt24fr169DKpWau2yT+emnn/DSSy9hx44dD92usLAQzz//PJycnHRjgwYNQkFBgYkrNB9DezFy5Eh8+OGHemNWVla4du2aKcszO0P7AQAajQbvv/8+Fi1aBLFYbPrizMzQXuTl5eGll17S60FCQgLmzp1r4grNx9Be9OjRAwCwc+dOtLa2Ij8/HxUVFRgwYIA5yjQ5T09PpKWlwdHRUTf2oPeBwsJCvPjii7C2ttZt9+KLLz4xf1iILF3Ak0alUqFHjx6wt7fXjT399NO4desWfv31V/Tp00dv+48//hhDhw7F4MGDzV2qyQnpxfnz5+Hn5weFQoF9+/bhmWeewbx58zBs2DBLlG50QnoRGBiIKVOmYMGCBVi4cCE0Gg1WrVoFLy8vS5RuEpMmTTJoO5VK1e7/md69e6OmpsYUZVmEob3w8PDQW7569Sq+/vprxMbGmqIsizG0HwCQnp6Onj17IiIiwqDA2dEY2ovz58/D1dUV69evh1KphJOTE6ZOndpp/gAHDO+Fm5sb4uLi8Oc//xnr1q2DRqNBbGwsXn75ZRNXaB69evXCH//4R91ya2srMjMz9cbuUKlU7d43evfujZKSEpPXaQgGyHs0NzfDzs5Ob+zOcktLi974Tz/9hAMHDjzy85EdlZBeXL9+Henp6Zg0aRLS0tJw6NAhzJkzB1988QX8/f3NVrOpCOlFU1MTqqqqEBMTg7CwMBw+fBgfffQRvL29IZPJzFXyE6G5ubndmXk7OzvcunXLQhU9GZqamjB37lz06dNHUODqTCorK5Geno7s7GxLl2Jx169fx65duzBixAj85S9/wS+//IKEhAT07NkTr7zyiqXLMyu1Wo2zZ89i7NixGD9+PH755ResXr0afn5+GDFihKXLM7rVq1ejuLj4vh/redDvnXt/51gKA+Q97O3t2/3HubPcrVs33diNGzewbNkyfPjhh3B2djZrjeZiaC8AwMbGBj4+PoiLiwMAPP/88zhx4kSnCZBCepGeno6WlhbMnz8fwO1enDlzBlu2bEFqaqp5Cn5C2Nvbt/vAd0tLCxwcHCxUkeU1NDRg9uzZqKqqwrZt29q9froCrVaLDz74ADExMXB3d7d0ORZnY2ODp556CitXroSNjQ38/f1RUlKC7du3d7kAuWvXLvz888/Izc2FtbU1/P39UVNTg40bN3aqAKnVapGYmIjt27djw4YN8Pb2brfNg37vPCnvnwyQ93BxccG1a9fQ0tKiS/4qlQp2dna6z2YAtz8Ie+7cOSxevFg31tzcjOXLl6OgoAAJCQlmr93YDO0FAPTp0wfPPvus3piHh0enuQtbSC+Kiorg6+urNzZgwABkZWWZrd4nhYuLS7vLLVevXsUzzzxjoYosq7a2FjNmzMDVq1fx2Weftft/pqu4ePEiTpw4gV9++QUbNmwAcPuP8qKiIhQWFuJvf/ubhSs0rz59+qC1tRU2Nja6MQ8PDxw9etSCVVlGUVERpFKp7nN/wO33z/T0dAtWZVytra344IMPsGfPHqxfv/6BfyS4uLhApVLpjT1J75+8ieYefn5+sLW11fuQ6okTJzBgwACIRG15e+DAgfj222+hVCp1P08//TTmzZunO/PU0RnaCwCQy+X45Zdf9MbOnDkDNzc3s9RqakJ60adPn3bBuby8vEuGhYCAAJSUlKCpqUk3duLEiS53KR9o+/qOuro6fP75553yrmNDubi44Ntvv8WuXbt0759+fn6IiopCYmKipcszO7lcjtOnT+t9tKMzvX8K0adPH5SXl+uNdbb3z6SkJOzZswebNm166FnVgIAA5OfnQ6vVArh91vLnn39+Yt4/GSDv0a1bN0RERGDFihU4efIk9u/fj7///e946623ANw+63Tjxg04ODhAIpHo/VhbW6N3797o3bu3hY/COAztBQBMnDgRlZWV+Pjjj3H+/Hl8+umnOHr0KCZOnGjJQzAaIb2YMGECDh8+jK1bt+LChQv48ssvkZOTg7ffftuSh2A2tbW1uH79OgBgyJAhcHV1xZIlS1BWVoa0tDQUFhZ2qpsDHubuXnz66ac4deoUVq9ejW7dukGlUkGlUqG+vt6yRZrRnX6IRKJ275/29vbo0aMHXFxcLF2mWdz92hg1ahREIhGWLVuGyspK7Nq1Czk5OV3m87F392L06NG4ePEiPvroI5w/fx7ffvsttm7d2mnePwsKCvCPf/wD8+bNg7+/v+594M6Zxrt/l7z66qtoamrCypUrcebMGaxevRrXr1/HqFGjLHkIbbTUTlNTk3bx4sVamUymffnll7Xp6em6dT4+Ptrs7Oz7Pm/o0KEPXNdRCelFfn6+duzYsVp/f3/tyJEjtfv27bNEySYjpBf/+te/tG+++aZWJpNpX3/9de0///lPS5RsFj4+PtrDhw/rlocNG6bduHGjbvns2bPayZMna/39/bWjRo3S/vjjj5Yo0ywe1os333xT6+Pj0+4nKirKUuWa3KNeG3eLiop64LrO4FG9KC8v17799ttaf39/7bBhw7RffPGFJco0i0f1oqCgQPunP/1JK5PJtCNGjNB+/vnnlijTJJKSku77PuDj46O9detWu98lhYWF2oiICK2/v7927Nix2qKiIgtWr89Kq/3t3CgRERERkQF4CZuIiIiIBGGAJCIiIiJBGCCJiIiISBAGSCIiIiIShAGSiIiIiARhgCQiIiIiQRggiYiIiEgQBkgiIiIiEoQBkoiIiIgE+f/oVZ8+pnDqiQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 648x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_forest({\n",
    "    'model_11': samples_11.posterior['sigma'],\n",
    "    'model_12': samples_12.posterior['sigma'],\n",
    "    'model_13': samples_13.posterior['sigma'],\n",
    "    'model_14': samples_14.posterior['sigma'],\n",
    "    'model_15': samples_15.posterior['sigma']\n",
    "}, combined=True, figsize=(9, 3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lG7JJ--0DTIR"
   },
   "source": [
    "The two models that omit predation, `model_11` and `model_13`, have larger values of `sigma`. This is because predation explains some of the variation among tanks. So when you add it to the model, the variation in the tank intercepts gets smaller."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aMrsnV4WDTIR"
   },
   "source": [
    "## Exercise 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EvFqYEslDTIS"
   },
   "source": [
    "> In 1980, a typical Bengali woman could have 5 or more children in her lifetime. By the year 2000, a typical Bengali woman had only 2 or 3. You're going to look at a historical set of data, when contraception was widely available but many families chose not to use it. These data reside in data `bangladesh` and come from the 1988 Bangladesh Fertility Survey. Each row is one of 1934 women. There are six variables, but you can focus on two of them for this practice problem:\n",
    "\n",
    "> 1. `district`: ID number of administrative district each woman resided in\n",
    "2. `use.contraception`: An indicator (0/1) of whether the woman was using contraception"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "executionInfo": {
     "elapsed": 735,
     "status": "ok",
     "timestamp": 1609493492049,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "2gJ9qonZDTIS",
    "outputId": "6b331fc3-424e-4e7c-a0fb-685d5d04fc55"
   },
   "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>woman</th>\n",
       "      <th>district</th>\n",
       "      <th>use.contraception</th>\n",
       "      <th>living.children</th>\n",
       "      <th>age.centered</th>\n",
       "      <th>urban</th>\n",
       "      <th>did</th>\n",
       "      <th>conc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>18.4400</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-5.5599</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1.4400</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>8.4400</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>-13.5590</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   woman  district  use.contraception  living.children  age.centered  urban  \\\n",
       "0      1         1                  0                4       18.4400      1   \n",
       "1      2         1                  0                1       -5.5599      1   \n",
       "2      3         1                  0                3        1.4400      1   \n",
       "3      4         1                  0                4        8.4400      1   \n",
       "4      5         1                  0                1      -13.5590      1   \n",
       "\n",
       "   did  conc  \n",
       "0    0     0  \n",
       "1    0     0  \n",
       "2    0     0  \n",
       "3    0     0  \n",
       "4    0     0  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = pd.read_csv('./dat/bangladesh.csv', header=0, sep=';')\n",
    "d['did'] = pd.factorize(d['district'])[0]\n",
    "d['conc'] = d['use.contraception']\n",
    "d.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OeV5mwHCDTIT"
   },
   "source": [
    "Now for the ordinary fixed effect model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 16015,
     "status": "ok",
     "timestamp": 1609493507992,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "R7hkawCELhdg",
    "outputId": "6b48e36b-011d-4070-9fab-eec6f049e4b4"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:06<00:00, 164.12it/s, 7 steps of size 5.45e-01. acc. prob=0.82]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 492.30it/s, 7 steps of size 5.00e-01. acc. prob=0.85]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 498.17it/s, 7 steps of size 4.86e-01. acc. prob=0.86]\n",
      "sample: 100%|██████████| 1000/1000 [00:01<00:00, 500.56it/s, 7 steps of size 5.51e-01. acc. prob=0.82]\n"
     ]
    }
   ],
   "source": [
    "def model(did, conc=None):\n",
    "    # priors\n",
    "    a =  numpyro.sample('a', dist.Normal(0, 1.5), sample_shape=d.did.unique().shape)\n",
    "    # likelihood\n",
    "    logits = a[did]\n",
    "    numpyro.sample('conc_hat', dist.Bernoulli(logits=logits), obs=conc)\n",
    "\n",
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "mcmc.run(random.PRNGKey(0), d.did.values, d.conc.values)\n",
    "\n",
    "samples_21 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 32037,
     "status": "ok",
     "timestamp": 1609493524449,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "21s7J_FsN913",
    "outputId": "28381abc-92d0-47ea-98c1-fe51dd0a4f8e"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:06<00:00, 151.97it/s, 7 steps of size 4.27e-01. acc. prob=0.88]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 422.27it/s, 31 steps of size 4.29e-01. acc. prob=0.87]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 350.04it/s, 15 steps of size 3.71e-01. acc. prob=0.90]\n",
      "sample: 100%|██████████| 1000/1000 [00:02<00:00, 413.46it/s, 7 steps of size 4.34e-01. acc. prob=0.86]\n"
     ]
    }
   ],
   "source": [
    "def model(did, conc=None):\n",
    "    # priors\n",
    "    a_bar = numpyro.sample('a_bar', dist.Normal(0,1.5))\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(a_bar,sigma), sample_shape=d.did.unique().shape)\n",
    "    # likelihood\n",
    "    logits = a[did]\n",
    "    numpyro.sample('conc_hat', dist.Bernoulli(logits=logits), obs=conc)\n",
    "\n",
    "kernel = NUTS(model)\n",
    "mcmc = MCMC(kernel, num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "mcmc.run(rng_key, d.did.values, d.conc.values)\n",
    "\n",
    "samples_22 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ncXzdoEFDTIV"
   },
   "source": [
    "Now let's extract the samples, compute posterior mean probabilities in each district, and plot it all:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "id": "MnoBMN6KDTIW"
   },
   "outputs": [],
   "source": [
    "model_21 = pd.DataFrame(logistic(samples_21.posterior['a'].mean(axis=0).mean(axis=0)))\n",
    "model_22 = pd.DataFrame(logistic(samples_22.posterior['a'].mean(axis=0).mean(axis=0)))\n",
    "aux = pd.concat([model_21.assign(model='model_21'),\n",
    "                 model_22.assign(model='model_22')]).rename(columns={0:'prob'}).reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 544,
     "status": "ok",
     "timestamp": 1609493535698,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "vlRMRkRrDTIW",
    "outputId": "07c4378c-9991-48ec-a60a-7d566e6ae95f"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray(0.3687478, dtype=float32)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logistic(samples_22.posterior['a_bar'].mean()) # use it for the red line in the next plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0Ao9qTbIDTIX"
   },
   "source": [
    "The blue points are the fixed estimations. The _orange ones_ points are the varying effects. As you'd expect, they are shrunk towards the mean (the _red_ line). Some are shrunk more than others. The third district from the left shrunk a lot. Let's look at the sample size in each district:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 202
    },
    "executionInfo": {
     "elapsed": 1158,
     "status": "ok",
     "timestamp": 1609495030038,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "AIAT_e3xWGAD",
    "outputId": "ce4cb0f7-afe5-40b1-96c0-7a2a4a42ea50"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAC2CAYAAAAiNbCqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABjiklEQVR4nO3dd3hU1dbA4d/UZNJ7IIR0Qq8KiDRFvGJBBVGKXSx4LVeED7ugFMV27V3Ui4WiqChWEAUUg9KkhBZpISSk18n0748DkwyZMJA6Sdb7PHkg+8yc2ZMzZZ199l5L5XA4HAghhBBCCCFqpW7uDgghhBBCCOHtJGgWQgghhBDCAwmahRBCCCGE8ECCZiGEEEIIITyQoFkIIYQQQggPJGgWQgghhBDCg0YLmrdu3cr1119fo/3nn3/mqquuYvz48SxZsqSxHl4IIYQQQogGo22Mnb7zzjssX74cg8Hg0m6xWHjqqaf47LPPMBgMTJw4kREjRhAREdEY3RBCCCGEEKJBNMpIc1xcHK+88kqN9oyMDOLi4ggODkav13PWWWfx559/NkYXhBBCCCGEaDCNMtJ80UUXkZmZWaO9rKyMwMBA5+/+/v6UlZV53N+mTZtQqVQN2sfTZbfbUatl6rc3kmPjveTYeC85Nt5Ljo33kmPj3Rr6+KhUKvr06VOjvVGC5toEBARQXl7u/L28vNwliK6NwWCga9eujdm1WqWnpzfbY4tTk2PjveTYeC85Nt5Ljo33kmPj3Rr6+KSnp7ttb9LTpuTkZA4ePEhRURFms5m//vqLvn37NmUXhBBCCCGEOGNNMtL89ddfU1FRwfjx43nwwQeZPHkyDoeDq666iujo6KboghBCCCGEEHXWaEFzbGysM6Xc6NGjne0jRoxgxIgRjfWwQgghhBBCNLgmndMshBBCCCEahsViITMzk8rKyubuSrOyWCy1zkM+FV9fX2JjY9HpdKd1ewmahRBCCCFaoMzMTAIDA0lISGi2LGPewGg01qgN4onD4SA/P5/MzEwSExNP6z6SP0UIIYQQogWqrKwkPDy8TQfMdaVSqQgPDz+jUXoJmoUQQgghWigJmOvuTP92EjQLIYQQQgjhgQTNQgghhBCiUZlMplNmT0tLS2Pq1KlN2KMzJ0GzEEIIIYQQHkj2DCGEEEIIcUrLli1j9erVVFZWkpubyw033MCqVavYu3cvM2bMoKKigg8//BC9Xk9CQgJPPvkkZrOZ6dOnU1JSQlxcnHNfu3fvZs6cOQCEhIQwb9685npaZ0SCZiGEEEII4VF5eTkLFixgxYoVfPDBByxZsoS0tDQ++OADMjIy+OKLLwgICGDevHksXrwYk8lEamoqU6dOZevWraSlpQHw2GOPMW/ePFJSUli6dCnvvvsu5557bjM/O88kaBZCCCGEEB517doVgMDAQJKTk1GpVAQHB2M0GklJSSEgIACA/v37s27dOux2O8OHDwegd+/eaLVK2JmRkcETTzwBKIVJEhISmv7J1IEEzUIIIYQQwqPaUrSpVCoyMjKoqKjAz8+PDRs2OAuGbNmyhZEjR7Jz506sVisAiYmJzJ8/n5iYGDZu3Ehubm6TPYf6kKBZCCGEEELUmUaj4Z577uGGG25ArVYTFxfH9OnTAZgxYwYTJ04kKSnJWa561qxZPPDAA1itVlQqFXPnzuXYsWPN+RROi8eg2eFwsG3bNkwmk7Otf//+jdopIYQQQgjhPcaOHev8/7Bhwxg2bBigTNl47733ABg9enSN+7300ks12nr06MHChQtd2hITExk4cGBDdrnBeQya77nnHvLz82nfvj2gDMFL0CyEEEIIIdoSj0FzXl4eixYtaoq+CCGEEEII4ZU8FjdJTEwkJyenKfoihBBCCCGEV/I40rxp0ybOP/98QkNDnasm161b1+gdE0IIIYQQwlt4DJp/+OGHpuiHEEIIIYQQXsvj9Izdu3dz1VVXMWTIEK688kp27tzZFP0SQgghhBDCa3gcaZ4zZw5z586lS5cupKen88QTT3hcGGi325k1axa7d+9Gr9czZ84c4uPjndsXLFjAN998g0qlYsqUKVx44YX1fyZCCCGEEKJWX24+wrM/7CaryEhMiIH/u6gzV/bt0NzdajE8jjQDdOnSBVBy8Z0ogXgqK1euxGw2s3jxYqZNm8bTTz/t3FZSUsL//vc/Fi1axIIFC5g3b14duy6EEEIIIU7Hl5uP8NCybRwpMuIAjhQZeWjZNr7cfKRZ+2UymRgxYkSt29PS0pg6dWqt2y0WC4888giTJk1i3LhxrFq1ymX7vHnz+PTTTxukrx6DZrVazerVqyktLeXnn39Gr9d73OnGjRsZOnQoAH369GH79u3ObQaDgZiYGIxGI0ajsdaSjK2Nw+HgaJGRzMIKKi225u6OEEIIIdqQZ3/YjfGk+MNosfHsD7ubqUcNY/ny5QQHB/PJJ5/w7rvvMnv2bAAKCgq49dZb+fnnnxvssTwOG8+bN4/58+fz/PPPk5yc7OzMqZSVlREQEOD8XaPRYLVanaPU7du359JLL8Vms3HHHXd43J/JZCI9Pd3j7RpDZWVlvR9bZQhmZUYpb687QLnJxiU9orlzSBzqinwcDkcD9bTtaYhjIxqHHBvvJcfGe8mx8V7eemwsFgtGo/G0bptV5P52WUXG096HO1999RVr1qyhsrKSvLw8Jk2axC+//MK+ffu4//77qaio4OOPP0av1xMXF8djjz2GxWLh4YcfpqSkhI4dO2K32zEajezdu5f58+fjcDgICQlh1qxZmEwmbDZbrX0877zzGDZsmHMwVq1WYzQaKSgo4Pbbb2fdunWn/DtZLJbTPra1Bs0ngtzIyEiee+6509rZCQEBAZSXlzt/t9vtzoB5zZo1HDt2zDl8PnnyZPr160evXr1q3Z+Pjw9du3Y9oz40lPT09Ho9ttVm55Wf9/LSqn3Otq//zmbn0VI+ue0cooN8G6KbbVJ9j41oPHJsvJccG+8lx8Z7eeuxSU9Px2AwnNZtY0IMHHETOMeEGE57H+7o9XoqKytZsGABK1as4IMPPmDJkiWkpaXxwQcfkJGRwRdffEFAQADz5s1j+fLlmEwmunTpwtSpU9m6dSt//fUXBoOBOXPmMG/ePFJSUli6dCkff/wx5557LhqNptY+GgwGjEYjNpuNGTNmMHXqVAwGAykpKYAyvUOn09V6f51OV+PY1hZE1xo0P/DAAzz//POMGjXKOYXC4XCgUqlqzBc5Wb9+/Vi9ejWXXHIJW7ZsITU11bktODgYX19f9Ho9KpWKwMBASkpKTrm/liyr2Mjba/bXaM/ILedgfrkEzUIIIYRodP93UWceWrbNZYqGQafh/y7qXO99nwg6AwMDSU5ORqVSERwcjNFoJCUlxTn7oH///qxbtw673c7w4cMB6N27t3NgNSMjgyeeeAJQRoATEhJO6/Gzs7OZNm0akyZNYvTo0fV+PrWpNWh+/vnnAXjxxRddRoHT0tI87vTCCy/kt99+Y8KECTgcDubNm8f7779PXFwcF1xwAb///jvXXHMNarWafv36MXjw4AZ4Kt7JbHXUmEN0QoVZ5jYLIYQQovGdyJLRGNkzalufplKpyMjIoKKiAj8/PzZs2EBiYiIAW7ZsYeTIkezcuROr1QooVajnz59PTEwMGzduJDc31+Nj5+XlceeddzJz5kwGDRpU7+dyKrUGzX/99Rf79u3jgw8+4OabbwaUaRYff/wx33zzzSl3qlarefLJJ13akpOTnf+/9957uffee+vT7xYjIkBP79hgtmYWu7Rr1SpiQ+t+OUQIIYQQ4kxc2bdDk6aY02g03HPPPdxwww2o1Wri4uKYPn06ADNmzGDixIkkJSWh0+kAmDVrFg888ABWqxWVSsXcuXM5duzYKR/jzTffpKSkhNdff53XX38dgHfeeQdf34a/kq9y1LISbc+ePfz4448sW7aMsWPHKjdWqejRo4dzSL2pNOdcooZ47O1Hirn23TSKjRYAVCp4blxvLuvdHh+tpiG62SZ56xwzIcfGm8mx8V5ybLyXtx4bb+1XUzMajXWel+3ub1jb37XWkebU1FRSU1O5+uqr0el0HD58mPj4eEJCQurUqbasR4dglt89mH9yy6m02EiK9Cc+3E8CZiGEEEKI0zBr1iwyMjJqtL/zzjtN1gePKedWrVrFhx9+SEpKCvv27ePf//43V1xxRVP0rVWJD/cnPty/ubshhBBCCNHizJo1q9Zt9UmZdyY8Bs1Llixh+fLl+Pj4YDQaue666yRoFkIIIYQQbYrHioDh4eFoNMo0Al9fX5meIYQQQggh2hyPI80Oh4Mrr7ySvn37OtOCTJs2DahKSyeEEEIIIURr5jFonjJlivP/jZkwWgghhBBCNKK/l8CqJ6E4E4Jj4YLHodc1zd2rFsPj9Ixu3brx22+/8cUXX1BUVER0dDQDBgxgwIABTdE/IYQQQghRX38vga/vheLDgEP59+t7lfZmZDKZGDFiRK3b09LSmDp1aq3bLRYLjzzyCJMmTWLcuHHOqtXp6elMmjSJ66+/nsmTJ5OXl1fvvnoMmh9++GE6duzIwYMHiYiI4JFHHqn3gwohhBBCiCa06kmwnJRlwmJU2luw5cuXExwczCeffMK7777L7NmzAZg7dy6PPfYYCxcu5MILL2yQ1HQep2cUFRUxbtw4li9fTr9+/bDb7fV+UCGEEEII0YSKM8+s/TQtW7aM1atXU1lZSW5uLjfccAOrVq1i7969zJgxg4qKCj788EP0ej0JCQk8+eSTmM1mpk+fTklJCXFxcc597d69mzlz5gAQEhLCvHnzPD7+qFGjnEX3HA6HM3nFCy+8QFRUFAA2mw0fH596PU84jaAZcCaTzs7OdnZGCCGEEEK0EMGxx6dmuGmvp/LychYsWMCKFSv44IMPWLJkCWlpaXzwwQdkZGTwxRdfEBAQwLx581i8eDEmk4nU1FSmTp3K1q1bSUtLA+Cxxx5j3rx5pKSksHTpUt59913OPffcUz62v78/arWasrIy7r33Xu677z4AZ8C8adMmPvroIz7++ON6P0+PQfOjjz7Kww8/TEZGBvfeey8zZ86s94MKIYRoWwrLzaCCUD99c3dFiLbpgseVOczVp2joDEp7PZ0oOR0YGEhycjIqlYrg4GCMRiMpKSkEBAQA0L9/f9atW4fdbneODvfu3RutVglHMzIyeOKJJwBlrnJCQsJpPX52djbTpk1j0qRJLkkrvv32W9544w3efvttwsLC6v08PQbNCQkJzJw5k27durFy5UpSU1Pr/aBCCCHahpySSlam57Bg3QHUKrh1SCLnd4kiKsi3ubsmRNtyIktGI2TPUKlUtbZnZGRQUVGBn58fGzZsIDExEYAtW7YwcuRIZzpjgMTERObPn09MTAwbN24kNzfX42Pn5eVx5513MnPmTAYNGuRs/+qrr1i8eDELFy5ssBojHoPm6dOnM3z4cLp168b+/fv57rvvJD+zEEIIj8pNVl5cuYdPN1RdEn5g2TZuHpzAA6O64KuT6X5CNKle1zRpijmNRsM999zDDTfcgFqtJi4ujunTpwMwY8YMJk6cSFJSEjqdDlBKZT/wwANYrVZUKhVz587l2LFjp3yMN998k5KSEl5//XVef/11AN566y3mzp1L+/btueeeewBllPvee++t1/NRORwOx6luMH78eBYvXuz8/frrr2fhwoX1etAzlZ6e7hz6b2rN+dji1OTYeC85Nt6rKY/N7uxSRr20hpO/ZTRqFT/cN5SUqMAm6UdLIe8b7+Wtx8Zb+9XUjEYjBoOhTvd19zes7e/qcaRZpVKxf/9+EhMTOXTokGTPEG1WfpmJjNwydmWX0iHEQJd2gXQI9WvubgnhtUorLTUCZgCb3UFppbXpOySEaLFmzZrlTExRXUOkkjtdHoPmhx56iKlTp5KXl0dUVJRzgrYQbUluqYnZ3+xg+dajzraOYQb+d8tAEiP8m7FnQniviAAf/PQaKsw2l/YgXy3hAfVP/ySEUNKs1TanuDWZNWtWrduMRmOt207Fw2SLGjwWN+nduzdffvkl69atY9myZXTp0qVOHROiJdudXeISMAMcLjCyaMMhrDa5+iKEO7GhBuZc2cOlTaWCeWN7EhtSt0upQogqvr6+5Ofnn3HwJ5SAOT8/H1/f01+U7HGkedGiRbz//vtYrVYcDgdarZYff/yxXh0VoqVZ/08BoMzFDPfXU1JpodJi59vtR7l9WJKMmgnhhlajZlSPdqREBbBuXx5qlYohKRGkRPmjVrf+kTEhGltsbCyZmZmnlWWiNbNYLM7FhGfC19eX2NjTz1PtMWj++OOPWbhwIW+88QajRo3iww8/PONOCdHSxYX5MW1wGBe3KyWwaAeV/h3YbuvI0n3gq5UMAELUxk+vpVdsCL1iQ5q7K0K0OjqdzpnCrS1rqgWRHoPmqKgooqKiKC8vZ+DAgbz66qsed2q325k1axa7d+9Gr9czZ84c4uPjndt//fVXXnvtNRwOB927d2fmzJltYj6OaLn+Fa/CZ+8b+G38ytnWMTieodcswd/3tAprCiGEEKIF8zinOTAwkJUrV6JSqVi0aBFFRUUed7py5UrMZjOLFy9m2rRpPP30085tZWVlPPvss7z55pssXbqUDh06UFhYWK8nIURjCy3di9+er1za1MUHCdi+EGySBUAIIYRo7TzmaS4rK+PQoUOEh4fz/vvvc/755zNw4MBT7vSpp56iV69eXHrppQAMHTqUtWvXArB27Vq++OILdDodhw8f5uqrr2bMmDGn3N+WLVvw8WmeOaOVlZVnNElcNJ2mOjZ6vZ6EjIVofv9vzY2hCeRduZjcClmEUZ28b7xXcxwbg4+eUG0lDlQUWn2oNJmb9PFbCnnfeC9vPjYaHwPlKgPHyiwE6jWE+9ixG0uau1tNqjGOT53yNN97770sWLAAgAcffPC0HqisrMxZZxyUijBWqxWtVkthYSFpaWl8+eWX+Pn5ce2119KnT59Tzsnx8fGR4iaihiY9NsXJ7tuDOxLRLpYInwD329soed94ryY/NkWHIe0t2LgAVGpC+t8G/ScrJXyFC3nfeC9vPTYlRgsL/zjIcz/+5cyJPiAxjBeu6U1sG6oj0NDHJz093W27x+kZQUFBrFy5koyMDPbv38/+/fs9PlhAQADl5eXO3+12O1qtEp+HhITQs2dPIiMj8ff35+yzz661c0J4jbhzwDekZvt5D4IEzEK4ZyyG7x6A9a+AuRxMpbDuBfjxcTCVNXfvhGjx9h0r49kfdrsUEdqwv4AvNx+RNHSNwGPQnJ+fz4cffsisWbN4/PHHmTlzpsed9uvXjzVr1gDK1IrU1FTntu7du7Nnzx4KCgqwWq1s3bqVlJSUejwFIZpARCe4+TvodiX4BkP7PnD9lxDTr5k7JoQXKz4Eu1fUbN+5DIoONX1/hGhlfs/Ic9u++K/D5JfJNKiG5nF6xi233ML555/v/P3bb7/1uNMLL7yQ3377jQkTJuBwOJg3bx7vv/8+cXFxXHDBBUybNo1bb70VgFGjRrkE1UJ4rehucOUbYCwEvR8YQuu1O4fDwd5jZazdm8fhgnKGp0bRPSaIqCDvnDcnxBmzVLhvdzjAUrcKXm1dWaWVUpOFAB8tgb5nnpdWtC5BBvevgUAfHTqtZCVraLUGzatXr2bTpk2sWLGCzZs3A8o0i1WrVnHJJZeccqdqtZonn3zSpS05uWpO6KWXXupcJChEi6L3U34awM6sEsa9uR6jRSkx/MHvB7mkRztmX9lDiqWI1iGgnXJyaTwpQ5J/JARENfrDOxwOMguNZBYa0aghNtSPmBZaidBmd7Ajq5hnv9/N1iNFdG8fxP9d1IWescHoNB4vGotWamBiOHqNGvNJlWnvHpFCsEHfTL1qvWoNmrt06UJRURE+Pj7ORXoqlUqCXSEaQFmllfk/7HIGzCd8uz2bmwcnStAsWoeQOBj7Dnw6AezHUzNq9EpbSMdGf/i/DhYy+cM/KTEqj90uyJf3b+5P1/ZBjf7YDW1PTilXv7kek1UJjtb/U8A1b63ny7sG06NDsMtt88tM7Dxawnfbswn31zOqRzs6RQWgl0JMrU5KVAALJw/goWXb+CevnCBfLdP+lcqg5PDm7lqrVGvQ3L59e8aMGcMVV1yBWi1nsUI0pJJKCxv2F7jdtvdYKf0Tw5q4R0I0ApUKks6DKb9B7i7l98guEFZLNpoGlFlYwe3/+8sZMANkl1Ry3+LNfHrrOYS1sBPTr7YccQbMJ1jtDj5OO8jcK3s6y5KXVFp4dfU+3v/tgPN2r/+Swf9uGcDglIhG76fFaudQQQVZRUZ89RriwvyIrueUs/25ZazZm8eWw4UMSYlkYFJYm8oMcSoatYqBcUF8f30MjtJs8A1GHx6GSkaZG4XHOc3vvPMO77zzjkv+u3Xr1jVqp4Ro7Qx6DUkR/uw8WlpjWwf5MhCtiUYHUV2UnyaUVWSksMJSo313dhnZJaYWFzTvO+Y+28g/ueVYbHZ81Moo8uH8CpeAGZSpHY98sY2lU84lMrDxnrfJamPlzhymLt7qnC6QEuXPOzecTWJE3bIM7TtWxvi31pNfrixq+2JzFilR/rx/0wA6hslnJaYy2P45+u/+D6wmpa37VfCv2RDcoXn71gp5HEJesWIFa9euZd26dc4fIUT9hPrpefiSmjklU6L8SY2SFHZC1Jf2FPN8tRrvXCB1qqu6l/aMcdt+eZ8YfHRV0y6OlZrc3u5AfgUlxsbNpnAov4J7F23BbLOjVatQqWDfsXJe+HFPjalop8Nmd7D0r8POgPmEfcfK+fOA+yt1bU7+Xvj63qqAGWDH55D+dfP1qRXzONIcGxvrtVVwWqPcUhP788rJKakkJtiXhAh/md/aSp0VH8bSOwbx2up9HCkycnnvGK7oE0P7FrpQSYj6yioy8k9uGYUVFhLC/YgL9ye4luwAnnQMNZASFVBjhHZElyhiguv+HjuUX8FvGXn8susY/eJDGdk1muR6nugezC9nT04Z+SV6LJlFJET4E3BSZowBiaGcnRDKXweqFlX2iAlieGqky+1C/Nz/vSIDfPDTe/zKr5cD+eVckOTPXX00hJsysWr92WWJ4uHVeeSVms54ZLjCbGXdPvcp1X7bl8fYflIgh39+dd++4S3oOQ78G39KTlvi8R1ksVgYPXq0My2cSqXi+eefb/SOtUVZRUamLt5CWrW5rv/qFs3sK3oQHSwnLq2NQa+hf2IYb8SehcliI9igQ6XyzhEw4Z3KzVYO51dwpMhIgI+WhAj/es8fbS77c8u58f00DhVUpaL793nJ3DEsiWC/M5+fGRnoy1vXncUjX27jj38KUKngom7teOTSrgT41i14PFRQzg3vb+BAnpJK74edObzxawZL7xhEp+jAOu0z/WgJk975w2UqyUMXd+G6gfH4V+tnh1A/XpvUj705pezJKSM50p/UdoG0P+kEID7cnwu7RdHe18oFHaHSpmLBNgsTz0lq9BPy9n525iZtI/KbR5xtCQFRxF/+EXVZGmXQaegVG8yOrJoloXt3DKlHT1sRTS3vDY1eWUMgGpTHT47bbrutKfohgDV7cl0CZoAfd+Ywtl8so4LbNVOvRGMz6DQYdLKqXZyZskoLn/55mHnfpjurgaVE+fPejf2JD/dv3s6dIbPNxrvr/nEJmEFZwHZht2j6xtVtUVNyVABvX3822SWVqFUq2gf74u9T99HWDfsLnAHzCUUVFpb8dZgHL+6KRn1mQUpppYW5K9JrzL1+6rtdDEuNrJHlIzrIl+ggX4Z0ch1dri7MX89rF/qj+eFBND/8ClpfRvS/A0fC7WfUt7ropMvDZ+2jro1lx0j9+1lUnd8/4/1pNWpuGJTAd9uPMrpLEJ1DVfyVbWf94QqGNsGixhYhcQio1OBwXSTKufeCn2TQaGgez/26devG6tWreffdd1m5cqUUImkklRYbn2/KdLvtqy1Hmrg3Qghvtz+vgrkr0l3K5+47Vs67a/djtp35/NHmVFBm5tttR91u23K4qF77DjLoSI0OJCUqoF4BM8AfGe7n0f6ekU+5yep226kUlltqnX5wqKCWwjCeFGei/3QcmgPHL9tbK9Gtfwn9n2+D7cz7eCZ8Sg6Bm9LNmoyfUFe4f56edI7wYcP1Qcwsn8N1myfxrOZV1twQSWJE/U4MzVYbB/PLOZBXToW5cf8ujSqiM0z4VMl9DsoI89DpkDKyefvVSnkMmh9++GFiYmKYOnUqHTp04MEHH2yKfrU5WrWq1qT7skK4FSvNhkPrYe9PcCxdqqSJ07btSBEACeF+TOwTzvmdQtGoVXy15Qj5pS2rfK5eqyY8wAcfrZqre4Vx37lh9IlVRlkjvGhNR9/4ELft/eJC8dOf+dUinUZFmL/7UfTAugb4BfuhxM1Ay4a3oMT9wEyD8allbrchVMmichKrTUlPdyCvvNaTDvWx7egXXor24DqoyEe3ZwU+/7sYcnfXuZuZhRXMXL6DkS/8yvnP/8LURVvYn1de5/01K60PdB4Ft/0Ct/8Cd66H4Q9AYHRz96xV8viuLCws5Prrrwega9eu/PDDD43eqbZIq1Fz07kJLN+a5XKirlWruKKP+1XTooUr2A+LroNj25XfVWq45DnoNb72Lx8hjgv10/HO5RH0Kf+dyANfYwyM58jVN/HCNn2LqxAX5u/DvNHJRJXuIubv+egPHuKm+FHsO3cMER28pxDJ4OQIooN8yCmpylTgp9dw7TlxNbJ1nMhXfKigAp1GRXy4f40BkPYhBu6/MJVHv9zu0h4baiChriOpVvfZM7BWNvpIM2HJEBQDJVmu7cMegCDX9GdHi40sWLef/60/iNlm58Ku0Tx4cReSIqt99lnN8PtrYD/pyom5HPb+UKc0huUmK09/t4tv/q66svHDzhwyi4x8eMuAGidpuaUmDuSXk2sOwHakmLhwP4K8sXx5SKzyIxqVx6DZZDKRm5tLZGQkeXl52O12T3cRddS1fSCf3zEQY+5+9JYyTLpAgtulkFrHBSbCi9nM8PsrVQEzKHPSVtwPsf2hfa/m65toEQZHW/Bbcz/aXOU1ZGAjKbu/5JmJXxHgJhdvTnElR4qMVKiCySoyel056bM0GWiWX+X8PST/Nc765xtUN30DeMdJZEKEP5/edg4r/j7KL3ty6dMxhGvOjq3xGW212Vm9+xh3fbIJi00ZBQn10/HRrQPpHuNave+Snu3Ra9S8t+4fzDY7fTuGcveIlLofn5A4ZVTXdlKO6sThjV+6PKQjXP8lrHoSdn8LhhAY/iB0H+uyKM1stfH66n0s/OOQs+3HnTlk5JXxya3nVC1mtVZCwT73j1XHkeasIiMr3EwF2pFVQmah0SVoziys4D+fbmHjoaqMJbcOSeTf5ycT5u89V0AaS25pJSarnTA/PX71nNrUWnj8K/znP/9hwoQJBAYGUlZWxuzZs5uiX22Sr7WMfjnLYNUTYC4D3xC4+GmIugw0Eji3KmW5sG2J+23Zf0vQ3IgO5pfzT245JquNpMgA4sP98GmB5YWDSvZBrusIJQ47/mtmQ8elSsBy3K7sEm77318cPr7QLjrIh3duOJtesSF4BWMRmlVP1GhWFR2E7O1KIFhX5flQegRUGgiOBd9gz/c5haQAK3d3KeOuDsWo/DSoAkygch0NP1xo5J5PNzsDZoDCCgszPvubjyYPJLTalIwwfz3XJFZylTYLe9kxNPHnoPavWZTltIUmwFXvweeTqwLnoBi4eD74NsGofWRnGPu28ndXayG45pXSI0VGFv15uEZ7xrFyDhdUVAXNPoHQ7Qo4urXm46SOqlP3zDa7u2nXAJiq5ZJ2OBx8s/WoS8AM8O66/YzsFs05Sa03aC6uMPPLnlye+X432SWVnJcayYxRnencznuu+jQXj0Hz4MGDWbhwIb6+vmRmZtKrl3yZ14mpFIoOg80EgTHu5xsd2Qjf/V/V75VF8MUUuOUHiDunyboqmoBao3whmGpWBEQvJ0iNZduRYq599w9naWW1Cv47vg8X92iHvqUFzgX73Tarju1ULl8fD5rzykzc8+lmZ8AMkFNi4o6FG/ni34Np5w3pLC0VUPCP+22FB+q+35ydsOw2yDl+cpEwFEa/DOFJddtfWS6snotq4/s4x00Tz4MrXlVGWY87UliByWqnU5Q/Y1J9KLeoWLKjjB1ZJRwrNbkEzWRtgQ8vQ2MqxfkKHDodBv+nZpBbUQiF+6E8F/wilADZ/6QMCVo9pF4MU36HooOg9YWwBAjuSJPR+ys/tbDYHC4nFNVVWqpdzVaplFHqLZ9AfrUR54ShEHt2nboWHehLUoQ//5w0hznIV+syul9UYWHpRvdzwFfvOsY5Sa5/96PFRjYfKuKHHdmkRAbwr+7RpEYHtsg0ouv25fGfRVucv6/adYxNhwr58q7BLS4zT0PzGDQ//vjjxMfHM3nyZF5//XWWL1/Oo48+6uluorrCg0owvOf4fPDQRBi3ADr0q7qNzQJpb7q//5ZPJWhubQLbwdD/gxVTXdt9AiG6e/P0qZUrLDfz4Od/Exdq4NbhPvjpVHy5z8qMz7bSIya43gUqmlxkZ/ftsQPApyrYyi6uZG9OzRLMR4srySoy1itoLqwwU1Rhxkerqd90D98QSBgG6V/V3NauZ932WXIUFk10DboPrIVv7oPxH9Vt1PXoVth4Uuq0/b/A3h+h/2Rnk69Ow+uXRTGgfDURuz4BnYFrz7uLxflJ6LXVgihjEXz/QM2T57XPQbfLoX3vqrbSbFg5E7Yuqmrrerkyghx00miuVg+RqcqPF4oM8KF7TFCN/Mu+OjUdQk96HYUlKlM+srcp1e/a9YTIrhDUvk6PHRHow0sT+nD9gg0UHU/156tT8/q1/ehQ7TWs16oJ89dhNPty59kBxPjD5mMOPthcTFSQ6yhzTkkl05Zs5feMfLRqFVa7g9d+2cfSKYPo2SHEbT9OZOxo7IIzZyq/zMTzP+6p0V5YYSH9aKkEzZ5usHPnTp588kkAHn30Ua699tpG71SrYjHCz3OrAmZQRgo+Ggu3/wqh8Uqbw177Ag5LHVMPCe/WdTRUFsO655Uvzfa94bIXISKluXvWKuWVmZjYVceltp8JXf8GWMoZ0ukKdo69i6PFFS0vaI5IheQLIGNVVZvWFy54HHyrrlaoVRBk0DLj3BD6BJYADtIrgnn6t5I6FZwA5dL1lsNFPPbldrZnlRDmr+f/LurMxT3aEVKHQiTo/eC8B5QAtLK4qr3L6NpPDjwpOuR+lHr/r0oWCd9uZ75Pd0E9wOaPoNc1ykkv0CXYQq9f5qA/8IvzJjEr7+L2YY+jD+5fdT9jIRz6o/b+Vw+aMze6BswA6cuVwLnX1Wf+XJpRqL+eZ8b14vr3NlBwvES2TqPi5Ql96Xhy0AzKKH5Iw42U94wN4Yd7B6Mqywbs2P2iiQwJRF0tz7a/j5YnRsURmrWW9n/cAyVHGN6uL1eNmYUjxjVP9t5jZcT4Wvjm6kBCK49g0gWxuTyMV3/ex3+v6eMyHzi/zMSfBwpZ8Ns/OBxw07mJDEwK85osMWabnbwyEwPjg7itlw4/lZmdJT68/mcJeWW1xChtyGmd4hQWFhIaGkpJSQm2Fpb/s9mVZMGOz2q2GwuVD/QTQbPWB/rfqnygn6zPpEbtYlvmcDjYk1PGd9uPsienlIu6t6N/QljTLJIKiITB90KPscoJk3+4JKNvRIF6uLzya4I2vups89u9jL55O8gbu7QZe1ZHgdFw+atw5E/Y+bVystV1tDIKV01MiIHlV4eR8NOtziCyR2B7Blz5Hv6hdUtnmZFbxqR30jAenwNaUG7moWXbCDbouKRn3UYAie4Ot66Cg79B/j+QNAyie9Z98drJxR6qq+uC9sBanpt/JKirMioElGdCtYDZebM/noc+VynTJUD53A+IhrKcmvv0DXX9fesn7h978/+UUWmtdwRdp6t7TDDL7xpMeXkpOGzo/YKJDTHUyELSKIqPEL3xPeXqrs0CvSfB0KnKdJdqOlt2o/lxivN3bfZmEr8Zj+O21UDVSZeqsohHQn8i9OuXnG0JIYkkDn+X4kqLM2g2W20s/OMgL67c67zdnwcKuWNYElNHpuJbh7SFDS3MT8/LozvQJ+czQla/AdZKzg1L4oJLn6OyXf3WA7QGHoPmu+66i6uuuoqQkBBKSkqYOXNmU/Sr9XA4av/wPrk97hwYeCdseFO5n1qjXMKvPtogGtSOrBKufnO988v/223ZnB0fymuT+jVN6XK1purESTSqKHsu6q3v1GjX5O8mzHwEaIHHITgGgq9QFkvVIsRRQvDqe1xHXUuPEv/TbXDrSuDMU1puPFjkfM9U9/KqvZybHF630WaAiE7KT0MI7ug+IG3fp+Z0htPVdTSsfR7sJ6VuO/ce0FX7vKg+Wl6duRws1ebSBsXAyFnw5Z2ut4tIrTnvuraTB/8oUHnXJf7TYiwmNn8D/P4qmEvh7MmgP7/ux+Z0WSrgl6eVk40TNn0AubtgwidVc8RNZWjWPlfz/tZKVPt/heiqoLm3by4Bf73kcjN10X66HFiIo/PTzrbDBUZeW10zG8i76/Yz7qzYOpdib0g+Og1DrOvRbvhvVWPBPyR+OwnLbb8CIc3VNa/g8Z12/vnnM2zYMAoLCwkPD2+Rk9qbVVA7SL0Edq9wbdcH1AyWAqJgxKPQ73qoKFBGL0LilEuXosFVmK28tHJvjS//vw4WsiunpGmC5hbiSKGRkkoLAT5aOoQYXC5jthRqm7nWKVA6WysuKlN8BFVues32kiNQnFmnIKWgXPk7qlQQbNBRYbJhttkpqrBgtnlJWtKQWKVS2qKJVYFzWBKMeQv8wuq2z4gucP1XylqEvD1K7uFRT0NMH9fbBcW4L20cEl/zalLnS+CahfDLU1CRDz3GwYBbax6Xvtcr86lPTv0w4HbQNP8I5RlxOGDHF/DNf6rajmyE7uPgsuddMr80uKJM2PpxzfbDfygLJ08EzTYTVOS630eJa8o6Q66b7B6A356vYMQDgHLlsqTS4nYBpM3uoNhYj4wpDak0G+3vL9Zst1aiO7Ydos88N3ZrclqnpxqNhogIqfNeJ/oA+NdsZZrG0c1Km1+48iEZmljz9j4BshCsiZRWWtl0qJD+cYHc3ktPgNrMnlIfXv2zlJ1ZpQxPbeScpi1AhdnK6l3HePyrHeSXmwkyaHn44q5c0rM9QQbXBP8Oh4PDhRVkFVViV4eQW2oi0k2+4GYTEAnteikp/arT6Jo2s0BTO9XEZVXdgq0BieFM6R/CmDgjQaV7MRui2W7twN7KUMK9KX9t7Flw28/KyYFKowxC1KdSmlYHiUPg5m+V0WSdv/sFaaHx8K858MPDVW0aHVz+srIIuDpDiDK9ImEIBbk5hMV2Ao2br+aobjDpM/h2mnLVIDgWLn4G2vWo+/NpLkWHlEWNJ9vxmTJlrTGDZquxZrGUE6qvHzKEQe+JSs7pk3VyLVGtCYyseRtQBr60VVddQv30dIry5+FBBhLUOaiAg/Yo5q43Eu4lc5px2JUrIu7U1t6GNMo1HbvdzqxZs9i9ezd6vZ45c+YQHx9f4za33347F1xwARMnTmyMbniP8GS47jPlg8JmgqDYBl3UIOrGX6/l8fMjGVa6gtBfXwGLkUGhCZx38TNkBbWwRWF1UFppIbPQSKXFRlSgb81V68Du7FLu+mSz8/cSo5UHl20jNtTAkE5VXxR2u4P1/+Rz50cbKalULl13ig7grevOcq3w1Zz8wuHyV+B/VyjpHEEZDbzyjRpzGUGZt/vTzhz+OlDI8NQIhqVGtsyV48Gx0PEcZSStuojUOn8OdQs200P1KT4rFjrb4oM7UjlpGRpvuwoRHKv8NCT/SOWnNjoD9L0BOg5U1qn4BEL8YOVvXhu/MHIqcghzFzCDMv2j00i45UcwlSgDMnXMINHszGVV78GTGfMb97EDopUR/6KDru36ANfRfZVKGfVP/wayNlW1979NOYGprl1vJf/3ydNyzn8E/KsGHDuG+bHsMi2Bn1/lzJiSoPen39hP8PWWYkP+UdDnOlj/Ss1tMX2bvj9e5rSC5gMHDnDw4EE6d+5MdHS0xykaK1euxGw2s3jxYrZs2cLTTz/NG2+84XKbF198kZKSklr20Ar5R7i8eVoyu91BZmEFx0pN+Pto6BDiV2PUsSUI8NVyse/f6Fc+U9VYeID4b6+j3eRfgHqMSHm5I4UVzFq+k7QD+fjrtVisdl6a2JdzksJdgp7PN2VyXnIgd/TSEmwrolwbxIfp8P7vBxiQGObMbXyooILbF/7FufGBjE3VUWZ28M7WSp78ZievTuxHgK+XzLmM6QO3/6LkfLWalBPasCSX0SCAfcfKuPrN3yk8npJqZXoOHcMMfDR5YMsLnA2hcOVrsPweOPi70ta+jzJNoY6L7AxFe+Hvha6NxYfx3fgO/GueMiLb1vkGKrmE65hPuFaB0fUbLfcGhlBlUWXpSZX5VColqG1Mge2U4isLx1SNLKs1MPYdJZiuLjQeJn4KBf9gLs5GH5GoXCE+eSQ8IgVu/AZ+fEzJABPYDi6YBUnDXW6mKTtK4De3u6YYNJcT+PWtcNtq7xhM02hhwG2QuQEOpx1v08GlLyifl22cyuGorTaO4qOPPuKnn36iuLiYK6+8kkOHDvH444+fcqdPPfUUvXr14tJLLwVg6NChrF271rn9+++/Jz09Ha1WS0REhMeR5i1btuDj0zyXLiorK/H1lbmtJ2j0PvxTpmftriMM6aDiWIWDnSU+3DygPY7ygibtS32PTYSfisjlkyA/o8Y2y8UvsD90aL2yxWi1SqBotVo93LJp6X0MfLXXyDkhxSSYd6Mz5lEY1ofF+325sn8ntMdHevz8/Nh4uITzSpcT9ufzyiVNlYqSPrfzc9hEesdFUFmh5P/NdgThbzxKj31v4bfnC9AHkN/3br6yD+HsLonojE372qgPvY8Pi9ON/O+PQ1zdI4QuYSo25dhZnl7EY5d0YVCUzeuOqScqlYpQg5oAezEOm40KXSgFRgcePv7d0ul0JBxehvbXuTU3BrYjf9wXHKvwstHmFqKtfN9otVriTenoP7/BZY627dz7OZoykVJT42bp0mm1tNNXoC89CHYbluB4jlkCqHSzuPUET8dGpVIRbFAToDZjdWjIN+uwWFznKSfo8jF8dKnb+xuvXc4Bq/dMCYz0V+FvzkNtLsPq345cqx9Gs/dmT2uM907Xrl1rtHkc/lmxYgUff/wxN954IzfddBNXXXWVxwcqKysjIKDqkqxGo8FqtaLVatmzZw/ffPMNL7/8Mq+99tppddzHx8dt55tCenp6sz326TBZbGjVKjRNkaYH2JdTRsfSHTxmegHfn3+CgHbk9J/Bjnxfhvfu0qSXZut9bEqzwew+B7bObiY1tW6FAYxmG0cKy7EWHwWVCn14RzqGGbym4tyRwgpGt9tCp2/HK5d5gUjgnrPuIss+xeVvGmlchV/1kXiHg6DNb/GvKwfhF191qS70yAEiV92JJv94UnxjIeG/z+aqAVOpCJpO+4SWMzJWbrJS+OcWvhvrS+ym2egOpjOuXX/uuPo/fHKwnBsG17HYhpdIT0+na0pn6vX1XF5LWevAGMKjYghvzDmprZi3f980KEtHZXR1x1fKVI3uY9C060lsXRdp1omSmUoPuFlh5OJMj02ou8acHbXe3hAQQtf23nnsfQBvv77W0O+d9HQ3i6c5jaDZ4XCgUqmcUzL0es+phAICAigvr5owbrfbnaNuX375JTk5Odx4440cOXIEnU5Hhw4dGDZs2Gk9EW9jLT2GrfgoKq0P+vB4ZS5bEzhWVIa98ADWgsOo9H7oIhIJjuiAj65xAzP/isOk/DQJyvOUhpIjRK/6D1z4KnllCUQHtaBREv8o6Hcj/Pp0zW1xA+u0S4vNTn5WBuE7PiV050JQa8nrdRu5na+kQ5x3XNryx0jUxnnOgPmEkI2voe4+BqgKiAw73ecvNmxZAN0vdqbZCjcdqQqYqwne/CYBA28EvGRe82nw1WmY0ddG5KKrnAuGfP75gdTDa7lrwnfN3Dsv0eFs5RK7sdC1/fxHGncRVzWl+dlKsKUzYAjv2DT5fUXD0fkqc2Tb0jzZ4I7Q6V9KBcnq4ocoi1SF1/MYNF922WVce+21ZGVlcdtttzFy5EhPd6Ffv36sXr2aSy65hC1btriM2M2YMcP5/1deeYWIiIgWGzCbD29Cu/zfaHPTQaXG0uMa7Oc9gk944774i0pK0ez5jqgf7wFrJQDWqJ6UXf4ePrF1rJ51moJK91YFzNVE//UcZV0uoC45XxtLeWkhjuKjgANVYAz+wSed+6vV0OdaOPQ77F9zvE0Lo+bXOVdsSWEeEWsfwzfje2dbxG9PUl6WSX7gLMJDmz85fJCjFPWhtW63GcoOAVUVy1Ra9yfJKq1emX94XK0p2yxGNHZznfvaHDQOK2HbFtRcYW+pIOzwj5Dcp1n65VUiUuCmFbB6nlKRMCQORj6pLHxrZGaLGevBNPx/mI46dxf4R1I+fCaVqZcQECLFgYQX8w2CS56DNc8cr+7oUBYbnvewchIqvJ7HoPm6665j0KBB7Nmzh6SkJDp39hyUXXjhhfz2229MmDABh8PBvHnzeP/994mLi+OCCy44817u3g13npT8/Zpr4N//hooKuOSSmve56SblJy8Pxo2ruf3OO2H8eDh8GK6/vub2adNg9Gj0+/fXfGzAfNdk9Psfh39y4XslcNXxHvbAb3GEJaF66ik491z4/Xd4+OEa9+fFF6FPH1i5EubMqbn9rbegc2f4+mt4/nmXTQGmcrSD9kAQsN0Cf5nR8geBr52HNSwFrUYDn30GERHwwQfKz8m+/Rb8/OD112HJkprbf/lF+fe55+Cbb5zNhtIcKCmHa49frPnVBPutoErHf+lVSlWq8HD4/HNl+0MPwfr1rvuOjYWPPlL+f999sGWL6/bUVHj7beX/t98Oe04awezTR/n7AVx3HWRmum4fNAjjf25B9+ND6J/+Cioc2A2h2EKT0OgNcMEF8Nhjym0n3QFlJWBNUYIkrR58cmHA8STz551X829zitdecGUZ2nY7oI8eKuywRAkm/VWv4tNuHfj4n/Zrj9274Y47am5/9FEYOVL5u913n+s2hx3rtNtxxBtQbctA88bnqE66+qGe+4iySnzLIVjjmrdY+/1s+LBH1WvvqXVw9KQ0Q2MMcM6d8PmXcGKBr6UCsiqU+YnXGMBPDVvMsNMHfp6snIycUMfXHgAGA3x3fLR39mxYtcp1e0O89l55AU3BXvjaCPmueXY1f34I5x0/8a/ltcdTTyn/v+oqyD8pE0D1197FF4PxpJONyy6D6dOV/5/haw84vc+9Xr0a7rVnt4E9UclCsvRpmDev0T73AFi4EKuqEL+HL4UNJ6ZWleP/2s3YIrvBdz83yuce0DSvvalTlf/X8XOv0V97N0+EwlyYcIvrexoa7Du3Tp970OivPe3jj0PXrrB4cdXnXnWn+5376QpYvA1sxwdmVu2Cl29q/tdePb9zm/1zb+DAhn3tuTvGnEbQvGvXLoxGI+3bt2fevHlMmTKFQYMGnfI+arWaJ590zW2YnFzz0vQ999zj6eG9lqP0qFKA5CTqsmxsgTE05iQJlcV4PGm+6+VITUUetuD4Rk10r9K7n9nk8AtHpfGOFfOWygp039yFNjPN2aY2FuKw78EW2a3msdHolJ/qv9eRylHLQgmHA9Wpyvo2BIcDe0U+2q/uhI5qOGyFo1bsUd1R+wZV3c4/EkY+AVtudr2/zheV9qTpRTp/ZdS94J/jCwHVMPwB6HAWbKx2iVHrC+EpkLe32n39ICKx5pert/MJgG5jgHU1t4V5mvnYxqg1yk8TcuxbXbMiH6AuycRYVoxBago0PIcDe14GqoVjUWVux5Fjg9AEVD6BLlecxGlSqZXPTNHieMyeMWHCBB577DFeeeUVpkyZwrPPPsvHH7upptOImnNxRG2PXbHje/yWjq95B5WaitvT8Gtft0Vkp6Ny5/f4LnHz2PoATHf8jk94I5YDNhbDH6+7zgMOioEbvqqZg9RqUqpm/bMG7GZIHA6RnaGWwPtM1XZsyg9uwv/9893ep/zm1fjH92uQx3fHkrsP3dtDXZPkAwREY5u8Ck1o46UUMuXsweedoc4pO07hKZhv+BZ9cLXFeMYiOLBWubxengvdx8LAKTVL9wLY7UqO8cpCJRdpcJz74gvmCijcD7m7MaND3767+/21BAUH4JNrIG93VVv8ECVVVXCHZutWQ2jpi82M3z6KYYObHLIRqRivW4EhpB5LHMuOKZ9b/hFNtj6lOm89NpbDm9C9P9J1ypJGh+WWVeg69G6+jjUhbz02QtEYCwHrlD1Dr9fTqVMnLBYLffr0QX2q6lJtiCM0QZmKcFJZXlPKxTiC2rm/UwNRR6SAT1CNhVzWgXfhE1rPJP6mUjCVKYt53H1pGIJh0F3Q+WIlIPYLO14kwc087oyfYdEk11Kyl78CvSY2bh5XN6NQTrbGLVWqC0vAevmraJdNrkqlpNFhveJ1tPUJmO02JXCtyFeOfWi88vqrfpOizJoBM0D+PmwlOVA9aDaEQNfREH8uWM3gF1H7MVGrISwBSDh1H/V+SjXL6O5kpKfTtaUGzKA83+uWQc525XUe3QOiuyq5ZUXz6nQhuAmajT0m1j1griiAfavg59lKye3OF8N5Dykn+W2d3Ybqr/dqzvG3WWDrYmgjQbMQcBpBs0qlYsaMGQwbNoxvv/0Wnc47LsE3N3V4IqVX/o/A5ZOV6kaALbIb1hEz8fcP8nDv+tFHpWC97gs0y+9GlZsOGh22syajPuvGul8qtZogazP8PAdydymjasOmuy/R6hukFImI6VP7/ooOw1d3uQbMAN9OVwK18JS69fM0qINilFXKxYddNwTHom7sUUKNFm2XS+GOtdgyN4JGh6ZDP7T1eb6mMkhfDt/NUE5qNDoYfB8MvMOlKpldW8sIvkaHTVNL1hs/L184VXZMmRpSnKlUdQtLqnNBjjMWEqv8dB7VNI93MpsFLEalUpkMVjiporpROfBefNNedrZZY8+B7mPqvtN9K2HZbVW/7/gCDq1Xqu+FNuKVuxbAYjGjKz7gdpumsGaOe9HArGZlMESmwngFj0Hzf//7X7Zt28awYcNIS0vjhRdeaIp+eT2Djw/lKSMouWk19qJMVDpfCIknOLJpLt1qO54NN32jjIpofdAEx9WoauZkLoeSI8pIQWA796t0s7fB+xdXBbk7v4CMlXDrqrqNtlTkKaOiJ7OalPzIjRg0G8JjMY19H5/F46v64BeOaewHGMIbuJyuOzpfaNcTTbsGyuebuwu+rLYY1WaBNc8qo7rVAgVrcEesMWehzdrocveKXjcrJxEtTfERpYpdRrVFLykjYfTLLX6KxCnZbUo+17S3IGcbpFwIvSfUOaNLa+MbHIll2DQqe1yFo+QIKkMYqrBEDCF1vMJXmqOMMNdoz1bee208aFbrfClOGUPwgZpz/Es6XUlI03epbbBalPf/+tcgfy90GQ09xkpVvmbmMWhes0ZJxfXVV18BkJ2dzZVXXtmonWop/H31EJOq/DRLB06jNHfhQVj1hDJy4rBD7AC4/GWIqjZXx1IJv79cc1TYVKqMwNQhaDZr/NBrdG6nQ5i0gTR2fUef+P5U3rwKR9EhZRFeaDy+kS10EdfO5e7b17+mBFQ+Sg7koPD2lFz2Ftr1L+OfvgS0vhT3nYKjz7WE+Ht7ano3Dqx1DZhBeT0eWg893ayQbi1ytsO7I8F2PFXf0a2w5SO46buWO0e8gen8Q9D59wH61H9nNrMSILvj7sS/jdGoVTiSzsfcvj/6o386200dh+LoeE4z9qyVO7oZ3h9VNS3m6FbY8rGyfqiNn8g1J49Bc0aGcvnF4XCQnp5OSEiIBM0thbkcx8onUO34vKotcwN8Mh5u+V5ZwAfKpR83paQByK1ZsOJ0ZBFJaJ87CN74qkt7ReqVHFNFeZod2yB8IxOhpQbK1dgdNtxenHfYgKp1vCqViqD2KRwbMZfKgffhUKsgoD0RQU2/oKnebNbjeUzd2PqpMsLexFkbmoTVBL+/WhUwn1CaDQd/k6C5MfiF40i9GFX6VzW3RcnCLwBDVAL5l76DNWcX2pJDWIMT0UZ3JiLSy65gWYxQeEBZ/6EPUKZzBbXAdQimMmWR9snzyAv3Q/bfEjQ3I49B87Rp05z/dzgc3OEuf6LwTsVHUO38omZ70UHshQeUub8APkHYu4xGnbO9xk3tnS50H7B5YHLoWKS5jCvOT6Ld9rfBbiG36038qhlIil1S7ZwJU6fLMKx/tUZ7Rb/b8fMJdGlTqVREhwZBaOPOq290ak3tFbJC4lpnwAzKVKrsre63Hd0KuMkxKupH70f54BkEHF6vzKE/rmTA/TgMHWn+ckTNz0eroX1sIgWhMVRYbATotIT6e64O3KTM5fD3Elhxf9VV04jOMPHTljelwVwGue7LOHNsl7KIWzQLj/GQ2Wx2/mRlZZF5clJr4bXMFkvNKRfHWczVRrLUaopSrsQe4joqa+o4jMIQNwsBT0NMiC9px7SM+T2BFzq+xCsJr3H15p4s+NtEbKhfnfbZVv2YG0ru8Ker8nqq1JT2uY0fKrpQZjpFphBv5nBASZYyImRxU01QpYKzblbymVan1iilz1uyymK3Od4B8AnCljDc7SZrx3MbsVNu2Froa+sMFZabmbyihB/O+YiDFy0g67zn2Dl6OXMLR7CjlsPUqhQfgQPrlKlPeXuVubS1CPP3ITbEz/sCZlAWDH9zn+t3Xt5uWPffGlmuvJ5vMMQNdr/tVAvwRaPzONI8atQoVCoVDocDX19fJk+e3BT9al4Ox/G8tMV08Fcpv7fAVatF2nDC2p+F9qjrwjB8gijzjXHOK64wW5n6UzHX9HuLntqD+JX8Q2lod34tjCD/bxPT6nB1K9BXxxOXd+eZ73fx6vqj2B0wPDWSmaO7ERHQ2DOaW5ed+Q4+/Kcb/7ngS8IooUIdwIfpcPRIBf8a0Ny9q4OyY7Dtc1j7DFSWQNcr4PyHldLM1WTqE2D0ImLWP4E6dwf2qB5kDZqJQ5+Al10UPj3leco87d9eUqZEnT0ZOl/isqjRgpqcTpOI3fGZS2BtiepJdmCPxn/edrsywrX5Y2VuddfLlRRvrfhycIXZxqZDRaTtd6DX+OGrC6CkUsmIdHbnWsrDtxa5u+CjcVWZhjR6uOo9SL24cdOCNoZjtYzMblsKwx9UMuG0FDoDDJ0G+35U1hadEDdYWQDubcrzobJIKWjlaTrMiSJZLTCmgtMImn/++eem6If3MJXCrm/g+4fAWEiQbwhcNFf5YvcN9Hh3b2LUBLLnnKfotvpW5SQAwCeIQ/96lwp1NCeSjdnsDgrKLdz1bTF++iDC/AdwrMSE2VbAjYPq/pw7hvnxzLheTL0wFYcDooN8CfBtYdXhvMCoHu14a80/3Oi8yFMEwKsT++Lv08L+nnY7bPkEVs6satvxOWRthJu+dQkgl27J5YPfVdzc7wW69VWzI8/O+8uLuW1oLveMCGn6vteHzQqbFsKqWVVt305X5ilf9pKS/xwoM1m5b5WR6ed/SlLxBgILdlDQbjCb6Mym7RZmNvYU/Zxt8N6/qvJ97/8VIrvCtUshpEWeqngU4KulT2wIfx4sxGyzY7ZVjVTGh7fABbSny1QK3z3gmprTZobPb4E717e8bC36APftvkEtczpXu+4weSXs/QmO7YTUURB7dtVaJG9gs8KRP+GbaXBsh5LH/qK5ygJ135OmCJbnQuZG2PiBktP/7MnQrmfN23m5FvaN2wSyt8MXU6p+ryxS8g2HJkDCkObqVZ10CDHw7vYItvV8l14BRWgcVg5Yw/l+nw+PdalaHBboq+O6gXE8sGwbFWYbFeaq0ZWLe9SvUItBryUpspYPM3FaUqMDeXZcL2Yt30G52YZOo+L2oUmck+zlOZbdKTkCa5+v2V54AAr2O4Nmi83OnwcKKDZaePE31wwGG/YXYLM70Khb0EhF0SFYM79m+44vYMhUMCgFIvz1GhIi/Jjw+RFiQ3sQHXQ2+3eXU1CexwvXNHKaPYsR1jxfs0BObroyn7qVBs3BBh0PX9qV8W/94RIwj+gSSVJkKw6aS7Phn19qttssyvuxpQXNkZ2VaQ2Vxa7tw2a0zMWAAFFdlJ+GZCyG0qPKSG9QB2f2pTrJ2w0fXl61cLn0KHx2C1z/JSRXq8prKoN1L0L1tTnbP4crXoM+17qOOtttyhWQ3d8qn5udL4WYvhBYrThXM5Kg+ST2je+7nejt2PAOqrhBLeqMVatRM+6sWH7ZrefuX4wYzVomDgxn2r861JiTNqxzJJf2bM+KbUcB5TV8z/kpdGnfhGeBNovyha0PaLGXbhqDv4+Wsf1iGZAYRkG5mSBfHbGhBnx0Lee16GQz16hk6WQpd/5Xp1EzrFMkv2fUTPk1LDWyZQXMoCzscTd3G5SS5sfptRomD0lixd/ZZBYaySxU7hMbauDs+LDG7WNlCRz5y/22o1uh62WN+/jNqGeHYL66ezDfbjvKgfxyLu3Znr4dQ1v3VDK1Tlkn4a6KqLYFPu/wZLhhOXxzv3LlSh8AQ++Hblc0d8+8R95e5e9zQEklTOpFMGo+hNXxEtben2pm+gFY+wJ0HAD64yedRQfhj9dq3u6HR5TByNCEqrYjm+DDS6vmoW/6H3QfC5c86znFbhOQoLk6h0M5I3K3yVyOyuFwu82bRQT4MO6sjozoEoXNDhEBelRuAtL2wQbmjenJlOHJFBnNRAX60DHMDz99E7xEbFbl8lPaW8q/nS+GHuMkvVY1GrWK+HD/ln+52D8C4gcr0xKqU6lrZMu4qEc07/++n5ySqkU87YJ8GdnVO0YczohfmFLJsFp2BgDUWhyB7aj+juzSLpAv7jqXpX9lsu1IMRd0ieKi7u2IC2/kBbS+QdDhbOVqwMnat+5SyVqNmq7tg+jalIMEzS2oA/S/DdafVJI8sH3dg6jmFtMHrvtMKa6l8VGqiLagga5GVZ4LS29S1iqcsOcHMBth/EIwhJz5Pktz3LdX5B0Ppo9/X5XnKfHVySqLlJP1E4xF8MPDNRdu7lgG59zZMoPme++9l5dfftnzDVsguwPyOo0naveKGttyUycShpoWtjTCKczf88hBsJ+Onn7NkGDpxDzKE2esWZuUs8ubvnE9AxUtn28wXPwM/G+0awaJS5+HENfFZokRASy6/Rx+3Z3LH/vzGZQUwbDUSBIj3Jw42KxQdEApt601KK+bprycV5qtBMQ6XwiOU/6tplgXhen854j65iaX1f0F5z6KVRNN9cLgKpWKLu2CeOSSrlhs9qa7oqAzwLDpyuKj6qPiUd1bfdDcJml1MOjfyuvxr3fBasIRPwTVpc+1zAqiJ/iFKT+tRFGFGaPFRqifHt/6fBYUHXYNmE84sEb53KxL0Nz5Ikh7vWZ7n2tdKw/7RyhXj08OnH1DXOc0m0qUOdLuFB5QRq+b2RkHzbNnuyk32kqo1Sp2qjvh03cKwVvecmbNKOk1mXRNF4a3tEvCLYHFqMxxPfkST/FhyPxTgubWqF0PuPVn5VKhuQwiUpWRLX3NkdTEiAASIwK4afApRr6sFtj7A3w+uepSc2QXmPBxo5ZrB5TPiMNpsOx25RKkWgN9r4fhD7gs2DlaXMnMtBAev+xronPXobFWkNtuOG/v0jOuyEaUm+94tVqFT1OPkkX3gFtX4di6GFX23zi6XY4q+YJWO5+5rctXh7M9/m40YWPQYKVYE0Y3TSy1ZEgXTajcbGXjgULmf7+LzEIjw1MjuXtECqnRdVyc76Y6r5P9FNtOJbqHcrXiz3eq2joOhG6Xu94uJB7HOXehOqnegOOiuaiqD5bo/JTP7Ly9NR8r0DvmpXsMms1mM0uWLGH//v106tSJq6++uin61Wyi2nVg5o7RXH/xaAKshZRpQ/h0j5qb2sW6ndYg6slUpiy+dCd7B/Rs3a+3tsoRmkC2uh1Wu4OIAD2G+kwDKvwHlt4I9mp5hXN3wconYcybboPxBpOfAQvHgKVC+d1uU1aHB7ZXFiCplRUSWrWKDYfLuPQQJIQPRKdRs++3MhyOSiYO9Z7PFZPdwU5TDN+YxmMKGEt0sR+XWsORiVKtj93u4IvNR5izonqqtgK6tj/Ghzf3JypIilA1p80HC7lhwQbn78u3ZrFuXx6f33mu+6ttngTHQmC7miXjw5OVqTp14R8BIx6FvtcpefcNocr+AqJcbmbV+nEg9Vb0QWfRIWMxVq0fWZ0mQlQvEqvHVf4ROC56CtXH41zu74jphyqikQdATpPHb6oHHniAmJgYBg0axMaNG3n44YeZP9/NKvBWIjU6kBuGdeep79LZfkRHt/ZqHrqkK53renYnTs03WJnjWri/5raOZzd9f0SjO1ZSyVdbsnh19T5KKy1c1L0d0y/qTHJds6wUHnQNmE/Y9TWUPdG48zPz9lQFzNWtf00ZcT6eDSQmxMBF3drx/Y5sDuRX3T4+3M+riv1sP1LM1W+ux17tKupHaYdYcseglj+fXrjIKjby8qqaI3rpR0s4mF8hQXMzKqm08PxPe2q0F5Sb2Xq4qI5BcwcY/xF8ck3V1LjAdjDugxpB7hkxhIChzymLrmQWGhm9YDdqlR8DEu/HXOlgw7ICOoSks3RKKJGBVdNH9xp6obv8c2J3vIOu4ii5SWM4HDWCOFUYzT+j+TSC5ry8PP773/8CMHLkSK677rpG71Rz0mrU9IsP5b0b+1NmslJelE+nhNYzP8rraPUw6C4lwKmeKiimH7Tv02zdEo3Dbnfw2aZMnvl+t7Ptu+3ZbDtSzNIpg2gfbDjFvWuhqeVjTOtTs6JgQ3O3cvxEe7W5y34+Wh65tCsatfJ87Q7onxDKU2N7Eu0lwUmFycpLK/e6BMwAOSUmthwukqC5lTFZ7ZRUuq/6WGqq4+V60SCMZhtZRe6z7RwprEfBndj+cNuvUHyoavF1cOMXfTlabMRosQGwenees/1AfgXHSiudQXNRhZlpX+xhd7aFkZ3vIyoMft1qYn9eJp9N6egV2WxqDZrNx8ssx8bG8vfff9OrVy927dpFQkJCU/WtWQUZdAQZdKQfLfV8Y1E/0d1g8k/KSt7sv5U0OB3PaZI3s2haWcVGXl+dUaM9s9DIwbzyugXNYUnKZUFjoWv7gDvqftnxdEV2VuYx222u7X2urbEQsWOYH89e3Zv7L+yMzeGgXZAvQQbvWVpcbray75j77EEZue7bRcsV5qene0wQO7JcU0CqVNAhxHuufrRFof46LuwWzUd/HKqx7ayEUDf3OH252iiyNEHKcdYYaIps/wad+1BTq1a5LG4sM1nZnlWCwwHf7nStYX+40MjZCY3Zy9NTa9BcvXx2Wloaer0es9mMj4/nSN9utzNr1ix2796NXq9nzpw5xMdXTfb+4IMPWLFCyVAxfPhw7r777gZ4KqJFi+ys/IhWzWKzU2ZyP7pVYbG5bfcoNAFu+Aq+uls56dLooN9NMOD22kehG0pYMlz9IXxxB5iVPNOOuEGoBt+rlCQ+iZ9eS3KUdxb7CTboGNIpgiV/ZdbY1i+ufl/UwvuE+uuZN6YnE9/5gwpz1Xvv0Uu60jGsDievosHoNRpuHpzITztzXFJujunbgZR6fH6kHy3hzo82OqeIdW4XwKsT+9Gpkaefdgwz0KdjCFsOF7m0j+/fkQ4hVa81f72W1OgAdmfXPEmPCfaOK3K1fqOcXD47Pz+f0NBQ1GrPlztXrlyJ2Wxm8eLFbNmyhaeffpo33ngDgMOHD7N8+XKWLl2KWq1m4sSJjBw5ki5dGrjqjRDC64T7+3BOUhh//OM6iqBWUb+5ve17K1WoynKUKT/BHZumQINWz+6QIRSP+pogUw5WrYFinw6kqNvR0rJJ67UabhuaxI87cyiqqLo8f25yOF3btaH8xW1Ir9hgvr57CDuyiimptNIrNpikSP+myc8vTik5MoDPppzLzqwSMosq6NkhmOTIAMLrOEUhp6SS2xf+xeGCqukdu7PLuG/xFj6aPLBGwbOGFB7gw8sT+vD6Lxks23QEnUbFrUOTmDCgo8tIc6i/nscv685176W5ZKfrnxBKopdU5/T4zkhLS+Phhx8mMDCQkpISZs+ezeDBg095n40bNzJ06FAA+vTpw/btVdkR2rVrx7vvvotGo/yhrFbraY1eC+FUfATKc+noh5IEvSVWr2qjggw6Zo3uzrXvppFfrkwBU6lg3tiexNV3dMs/XPlpQgfyyhn35gZKnaPnFcBeHrxYwx3Dklpcxp1O0YEsu/Nc0vYXsDenlHOTI+jZIZhoLxnlEQ1LpVKRHBXgtVc/2rqOYX50DGuYqTJZRUaXgPmEHVklZBUZGzVoBogL9+fJK7pz94gU1CoV0UG+biu7nh0fyhd3nsuHvx8gs8jIVf1iGdopgqhA7/gM8hg0v/jii3zyySdER0eTk5PD3Xff7TFoLisrIyCg6k2o0WiwWq1otVp0Oh1hYWE4HA6eeeYZunXrRmLiqVe3m0wm0tPTT3mbxlJZWdlsjy1c+eg0xFr+Qf/tf6A0mwCNDtvZt1PU7TqOVbSs4KS1O9X7RqPRsPCGXmSVWqmw2IgL8cHPXsH+fW5yc3q5A+bAagFzlTd/zeC8BH8cFYVu7tW8Tucz7exQLQMjDJjNhRRkFVCQ1USda+Pk+8Z7tfRjU6kJqXWb0WgkPb1p3+TFR2vf5gPce04YDpUam7mSkuyDlGTXfntouuPjMWjWaDRERysXGqOjo09rVDggIIDy8nLn73a7Ha226qFMJhMPP/ww/v7+zJw50+P+fHx86Nq1q8fbNYb09PRme2xxkpwd8Na1VenFbBY0aa8RHpFMeP/Jzds34eJ03jfdmqgvjSnjb/dfNBUmG4FBwXSIb9fEPfJMPtO8lxwb79XSj82xkkriwvw4VOCaIrNHTBDJ7cMI9fe+z6oz0dDHp7YA3OME5YCAABYuXMiuXbtYuHAhwcGeyyz369ePNWvWALBlyxZSU1Od2xwOB//+97/p3LkzTz75pHOahhAeZf7lPh/v768oJYyFaGIpkQG4m4FxeZ/2RAY27uVOIYSrSouVPTmlrN51jA3788kpqWzuLnmNqCBf3r7+LJKq5Xju2j6QFyf0rdfUjNzSSn7bl8c7a/7hxx3ZZBa6yVvfingcaX722Wd5/fXX+e9//0tycjLz5s3zuNMLL7yQ3377jQkTJuBwOJg3bx7vv/8+cXFx2O12NmzYgNlsZu3atQDcf//99O3bt/7PRrRuNlMt7a45cYVoKokR/rxwTW9mfPY3FpuycqVr+0DuGdEJvVYGBIRoKhVmK19tzuKRL7c5c42nRPnz7g39SahLMZBWqEv7IJbcMYisIuPx1IIGwuqR+zinpJKHlv3Nz7tynW1xYQb+d8vAVvs39xg0z5o1i+eff/6MdqpWq3nyySdd2pKTk53/37Zt2xntTwhAqWnvTv/JENDSchWI1sBHp+GSnu3p2SGEo8VGDDoNceF+XrNoRYi24kBeBQ994Rpb7DtWzuu/7GPOlT3kJPa4iEAfIgIbZvH8jqxil4AZ4FCBkWWbMrlvZCpqNwv9WjqP0zPMZjO7du3CZDJhNpudRU+EaHLhnWDMm6CtFpB0vgR6XYPba+RCNAEfrYaUqACGdork7IQwCZiFaAZ7ctwXIvtqSxZ5ZRK3NIY1u3Pdtn+3PZviytZZVdLjSPOBAwf497//7fxdpVKxatWqRu2UaJtySioxWewE+2kJNriZY6X3gx7jlFKgpdkY7VoM7buAIaTJ+yqEEMJ7+Pu4H0kO8tWhkUGVRpEY6T5VYVy4n0v+5dbEY9D89ddfY7PZKCgoIDw8/LSKmwhxJkqNFtbszWXuinSyiivpGxfCzNHd6dUhuOblHY0OwlMgPIUD6el0lYBZHJdVZCSnpBJfnYYOIQavKlEthGhcKVGBBPlqKal0XSx+94gUyTPeSAanRBDoo3VJu6lSwb/PS8bQSoNmjxHwTz/9xIUXXsjtt9/ORRddxG+//dYU/RJtyKZDhdz1yWayipWVzpsPFTHh7fXsy61ZSlOIk9ntDtL253PFa78x5vXfufiltdy7aDOHC1r3Km4hRJXECH8+vm0gvWOVDF/+eg0zLurMJT1bdio1b5YSFcCSKYMY1aMdIX46zooPZdFt59AjxnOWtZbK40jza6+9xtKlSwkPDycvL48pU6Z4LG4ixOkqq7Tw8s81i1pUWuz88U8+qdGBzdAr0ZIcLKjgpgV/YrTYnG2/7M7l9V/28cTl3WUBkBBtRM8OIXxwywAKykzotRpiQgxuq86JhtO1fRD/vaY3xUYLfnptq7/C53GkOSQkhPBwpTRtRESES6U/IerLZLWTX8sijZxiybEpPDuUX+4SMJ/w2cZMydMqRBsT6qcnOSqQjmF+EjA3EYNeS7vgtjElzuNIs7+/P5MnT6Z///7s2LGDyspKXnjhBUDJryxEfYT66bmybwdeXFlztHlYamQz9Ei0Jo7m7oAQQohWw2PQPHLkSOf/T5TTFqKhqNUqxvaL5aedOezIKnG233huAp2i5aqG8Cwu3B+DTlNjtHls31jaBckCICGEEA3DY9A8ZsyYpuiHaMPiwvxYcFN/DuSVU1BhpkOIgfgwP4L9pAyx8Cw+zI/3b+7P3Z9scuZjHdIpgrtHpMh8ZiGEEA3GY9AsRFOIDvIlWkYFRR2o1SrOSQrnq7sGk1NiwlenJjbUr03MrxNCCNF0JGgWQrQKHUL96BDq19zdEEII0UpJ0CxanAqzlcIKCxq/1psLUgghhBDeRYJm0aLszi7hhR/38OveXDqEGHhglJpByeEE+sqleCGEEEI0HqmJLVqM/XnljH/7D37YmUOlxU5Gbjm3L9zIn/sLmrtrQgghhGjlJGgWLcb2I8UUVVhqtD/30x6KKtwXSBFCCCGEaAgSNIsWI7vY6LY9p7gSk8XexL0RQgghRFsiQbNoMfrEhbptv6h7NKEBMqdZCCGEEI1HgmbRYiRH+nP9OfEube2Dfbl1aBJ6jRSxEEIIIUTjkewZosUI8/fh/n+lMrZfB3ZllxLqq6ZHxzBiJTevEEIIIRpZo4w02+12Hn/8ccaPH8/111/PwYMHXbYvWbKEsWPHcs0117B69erG6IJopUL99PSNC2XigDgS9GUSMAshhBCiSTTKSPPKlSsxm80sXryYLVu28PTTT/PGG28AkJuby8KFC/n8888xmUxMmjSJwYMHo9frG6MrohVzOBzN3QUhhBBCtBGNMtK8ceNGhg4dCkCfPn3Yvn27c9vff/9N37590ev1BAYGEhcXx65duxqjG0IIIYQQQjSIRhlpLisrIyAgwPm7RqPBarWi1WopKysjMDDQuc3f35+ysrJT7s9kMpGent4YXT0tzfnY4tTk2HgvOTbeS46N95Jj473k2Hi3hjw+JpPJbXujBM0BAQGUl5c7f7fb7Wi1WrfbysvLXYJod/r06dMY3RRCCCGEEOK0NMr0jH79+rFmzRoAtmzZQmpqqnNbr1692LhxIyaTidLSUjIyMly2CyGEEEII4W1UjkZYTWW325k1axZ79uzB4XAwb9481qxZQ1xcHBdccAFLlixh8eLFOBwO7rjjDi666KKG7oIQQgghhBANplGCZiGEEEIIIVoTqQgohBBCCCGEBxI0CyGEEEII4YEEzUIIIYQQQnjQKCnnWroTCxl3796NXq9nzpw5xMfHN3e32rytW7fy3HPPsXDhQg4ePMiDDz6ISqWiU6dOzJw5E7VazgGbg8Vi4eGHH+bIkSOYzWbuvPNOUlJS5Ph4AZvNxqOPPsr+/ftRqVQ88cQT+Pj4yLHxIvn5+YwdO5YFCxag1Wrl2HiJMWPGOOtNxMbGMn78eObOnYtGo2HIkCHcfffdzdzDtuutt97i559/xmKxMHHiRAYMGNBk7xt5N7pRvQz4tGnTePrpp5u7S23eO++8w6OPPupMOP7UU09x33338cknn+BwOFi1alUz97DtWr58OSEhIXzyySe8++67zJ49W46Pl1i9ejUAixYt4r777uO///2vHBsvYrFYePzxx/H19QXkc81bmEwmHA4HCxcuZOHChTz11FPMnDmT559/nk8//ZStW7eyc+fO5u5mm5SWlsbmzZv59NNPWbhwIdnZ2U36vpGg2Y1TlQEXzSMuLo5XXnnF+fuOHTsYMGAAAMOGDeP3339vrq61eaNGjeI///kPAA6HA41GI8fHS4wcOZLZs2cDkJWVRVBQkBwbLzJ//nwmTJhAVFQUIJ9r3mLXrl0YjUZuueUWbrjhBv7880/MZjNxcXGoVCqGDBkix6aZrFu3jtTUVO666y6mTJnCeeed16TvGwma3aitDLhoPhdddJGzqiQowZlKpQKUUuylpaXN1bU2z9/fn4CAAMrKyrj33nu577775Ph4Ea1WywMPPMDs2bMZPXq0HBsvsWzZMsLCwpwDNCCfa97C19eXyZMn89577/HEE0/w0EMPYTAYnNvl2DSfwsJCtm/fzksvvcQTTzzB9OnTm/R9I3Oa3ThVGXDhHarPVyovLycoKKgZeyOOHj3KXXfdxaRJkxg9ejTPPvusc5scn+Y3f/58pk+fzjXXXOOc4gRybJrT559/jkqlYv369aSnp/PAAw9QUFDg3C7HpvkkJiYSHx+PSqUiMTGRwMBAioqKnNvl2DSfkJAQkpKS0Ov1JCUl4ePjQ3Z2tnN7Yx8bGWl241RlwIV36NatG2lpaQCsWbOGs88+u5l71Hbl5eVxyy238H//93+MGzcOkOPjLb788kveeustAAwGAyqVih49esix8QIff/wxH330EQsXLqRr167Mnz+fYcOGybHxAp999plzLVNOTg5GoxE/Pz8OHTqEw+Fg3bp1cmyayVlnncXatWtxOBzOYzNo0KAme99IRUA33JUBT05Obu5utXmZmZncf//9LFmyhP379/PYY49hsVhISkpizpw5aDSa5u5imzRnzhy+++47kpKSnG2PPPIIc+bMkePTzCoqKnjooYfIy8vDarVy2223kZycLO8dL3P99dcza9Ys1Gq1HBsvYDabeeihh8jKykKlUjF9+nTUajXz5s3DZrMxZMgQpk6d2tzdbLOeeeYZ0tLScDgcTJ06ldjY2CZ730jQLIQQQgghhAcyPUMIIYQQQggPJGgWQgghhBDCAwmahRBCCCGE8ECCZiGEEEIIITyQoFkIIYQQQggPJGgWQogWxGQyMWLECObOnUtWVlatt1m6dGmN9tzcXGbNmlXrvv/880927drVUF0VQohWRYJmIYRogR555BFiYmLcbsvNzXUbNEdGRp4yaP788885duxYQ3VRCCFaFakNLYQQXq68vJzp06dTUlJCXFwcUFUQo6ioiPnz56PVajEYDLz00ku8+eab7Nu3j1dffRWHw8HmzZupqKhg7ty5PPTQQyxZsoTVq1c7t3fv3p3x48ezdu1aduzYQUpKSq0BuRBCtFUy0iyEEF5u0aJFpKam8vHHHzNhwgSXbStXruTiiy/mo48+YuLEiZSUlDBlyhRSUlK4++67AUhKSmLRokX4+PgAYLVamT17Nm+//TbLli0jLi6OsLAwhg4dyv/93/9JwCyEEG5I0CyEEF7uwIED9OzZE4DevXuj1VZdJJwyZQrHjh3jxhtv5Pvvv3fZdkJiYqLL74WFhQQFBREeHg7AbbfdJoGyEEJ4IEGzEEJ4ueTkZLZs2QLAzp07sVqtzm3Lly9nzJgxLFy4kE6dOrFkyRLUajV2u915G7Xa9aM+PDyckpISioqKAJgzZw5///03KpUKh8PR6M9HCCFaIpnTLIQQXm7ixInMmDGDiRMnkpSUhE6nc27r1asXjz76KAaDAbVazZNPPkl4eDgWi4Vnn30WX1/fGvtTq9XMnDmTO+64A7VaTbdu3ejZsyc7d+7kueeeIzY2luTk5KZ8ikII4fVUDhlWEEIIIYQQ4pRkeoYQQgghhBAeSNAshBBCCCGEBxI0CyGEEEII4YEEzUIIIYQQQnggQbMQQgghhBAeSNAshBBCCCGEBxI0CyGEEEII4cH/A2icTeO4EJ7aAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x180 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = sns.scatterplot(data=aux, x='index', y='prob', hue='model', s=50)\n",
    "g.set(ylim=(0,1), xlabel='district', ylabel='prob. use contraception')\n",
    "g.figure.set_size_inches(12,2.5)\n",
    "g.axhline(0.3679, ls='--', c='r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "PBbjzGMkDTIY"
   },
   "source": [
    "District 3 has only 2 women sampled. So it shrinks a lot. There are couple of other districts, like 49 and 54, that also have very few women sampled. But their fixed estimates aren't as extreme, so they don't shrink as much as district 3 does.\n",
    "\n",
    "All of this is explained by partial pooling, of course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 100
    },
    "executionInfo": {
     "elapsed": 857,
     "status": "ok",
     "timestamp": 1609494963912,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "jFApQAMBDTIY",
    "outputId": "448c0b4e-89aa-40fb-ff0f-e0d538ad89c9"
   },
   "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>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>10</th>\n",
       "      <th>...</th>\n",
       "      <th>51</th>\n",
       "      <th>52</th>\n",
       "      <th>53</th>\n",
       "      <th>55</th>\n",
       "      <th>56</th>\n",
       "      <th>57</th>\n",
       "      <th>58</th>\n",
       "      <th>59</th>\n",
       "      <th>60</th>\n",
       "      <th>61</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>district</th>\n",
       "      <td>117</td>\n",
       "      <td>20</td>\n",
       "      <td>2</td>\n",
       "      <td>30</td>\n",
       "      <td>39</td>\n",
       "      <td>65</td>\n",
       "      <td>18</td>\n",
       "      <td>37</td>\n",
       "      <td>23</td>\n",
       "      <td>13</td>\n",
       "      <td>...</td>\n",
       "      <td>37</td>\n",
       "      <td>61</td>\n",
       "      <td>19</td>\n",
       "      <td>6</td>\n",
       "      <td>45</td>\n",
       "      <td>27</td>\n",
       "      <td>33</td>\n",
       "      <td>10</td>\n",
       "      <td>32</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 60 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           1   2   3   4   5   6   7   8   9   10  ...  51  52  53  55  56  \\\n",
       "district  117  20   2  30  39  65  18  37  23  13  ...  37  61  19   6  45   \n",
       "\n",
       "          57  58  59  60  61  \n",
       "district  27  33  10  32  42  \n",
       "\n",
       "[1 rows x 60 columns]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(d.district.value_counts().sort_index()).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Uur4Nc9ZDTIY"
   },
   "source": [
    "## Exercise 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "umXdqjrgDTIY"
   },
   "source": [
    "Return to the Trolley data, (`Trolley`), from Chapter 12. Define and fit a varying intercepts model for these data. By this I mean to add an intercept parameter for the individual to the linear model. Cluster the varying intercepts on individual participants, as indicated by the unique values in the id variable. Include action, intention, and contact as before. Compare the varying intercepts model and a model that ignores individuals, using both WAIC/LOO and posterior predictions. What is the impact of individual variation in these data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 142
    },
    "executionInfo": {
     "elapsed": 636,
     "status": "ok",
     "timestamp": 1609496047750,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "kkSxAc03DTIZ",
    "outputId": "5b2ff347-c218-412e-d885-10ff39c474b2"
   },
   "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>case</th>\n",
       "      <th>response</th>\n",
       "      <th>order</th>\n",
       "      <th>id</th>\n",
       "      <th>age</th>\n",
       "      <th>male</th>\n",
       "      <th>edu</th>\n",
       "      <th>action</th>\n",
       "      <th>intention</th>\n",
       "      <th>contact</th>\n",
       "      <th>story</th>\n",
       "      <th>action2</th>\n",
       "      <th>edu_new</th>\n",
       "      <th>age_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>cfaqu</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>14</td>\n",
       "      <td>0</td>\n",
       "      <td>Middle School</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>aqu</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-1.650275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>cfbur</td>\n",
       "      <td>3</td>\n",
       "      <td>31</td>\n",
       "      <td>0</td>\n",
       "      <td>14</td>\n",
       "      <td>0</td>\n",
       "      <td>Middle School</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>bur</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-1.650275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>cfrub</td>\n",
       "      <td>4</td>\n",
       "      <td>16</td>\n",
       "      <td>0</td>\n",
       "      <td>14</td>\n",
       "      <td>0</td>\n",
       "      <td>Middle School</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>rub</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>-1.650275</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    case  response  order  id  age  male            edu  action  intention  \\\n",
       "0  cfaqu         4      2   0   14     0  Middle School       0          0   \n",
       "1  cfbur         3     31   0   14     0  Middle School       0          0   \n",
       "2  cfrub         4     16   0   14     0  Middle School       0          0   \n",
       "\n",
       "   contact story  action2  edu_new   age_std  \n",
       "0        1   aqu        1        1 -1.650275  \n",
       "1        1   bur        1        1 -1.650275  \n",
       "2        1   rub        1        1 -1.650275  "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = pd.read_csv('./dat/Trolley.csv', header=0, sep=';')\n",
    "elvl = d['edu'].unique()\n",
    "idx = [7 , 0 , 6 , 4 , 2 , 1, 3, 5]\n",
    "cat = pd.Categorical(d.edu, categories=list(elvl[idx]), ordered=True)\n",
    "d['edu_new'] = pd.factorize(cat, sort=True)[0].astype('int')\n",
    "d['age_std'] = (d['age'] - d['age'].mean())/d['age'].std()\n",
    "d['id'] = pd.factorize(d['id'])[0]\n",
    "d.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sExJeEvPDTIa"
   },
   "source": [
    "First, let's load the data and re-run the old model from Chapter 12:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 785799,
     "status": "ok",
     "timestamp": 1609497392955,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "nj73M5KXbaHq",
    "outputId": "24f34729-a40f-4c5e-ce18-34dfe76b0d51"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [00:29<00:00, 34.45it/s, 31 steps of size 2.09e-01. acc. prob=0.95]\n",
      "sample: 100%|██████████| 1000/1000 [17:08<00:00,  1.03s/it, 1023 steps of size 5.34e-06. acc. prob=0.83]\n",
      "sample: 100%|██████████| 1000/1000 [00:17<00:00, 56.43it/s, 15 steps of size 2.82e-01. acc. prob=0.93]\n",
      "sample: 100%|██████████| 1000/1000 [00:18<00:00, 53.22it/s, 23 steps of size 2.20e-01. acc. prob=0.93]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "                  mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "          bA     -0.29      0.33     -0.45     -0.54      0.27      2.04      7.41\n",
      "          bC     -0.54      0.34     -0.37     -1.12     -0.27      2.05      6.02\n",
      "          bI     -0.48      0.32     -0.32     -1.03     -0.23      2.04      6.94\n",
      "         bIA     -0.79      0.62     -0.47     -1.85     -0.35      2.02      9.64\n",
      "         bIC     -0.50      1.28     -1.19     -1.35      1.72      2.01     15.92\n",
      "cutpoints[0]     -1.98      1.14     -2.62     -2.69     -0.01      2.00     26.74\n",
      "cutpoints[1]     -0.61      2.31     -1.92     -1.99      3.40      2.00     59.25\n",
      "cutpoints[2]      1.26      4.51     -1.33     -1.40      9.07      2.00    121.92\n",
      "cutpoints[3]      2.86      5.49     -0.29     -0.36     12.36      2.00    156.60\n",
      "cutpoints[4]      3.41      5.29      0.38      0.31     12.58      2.00    151.02\n",
      "cutpoints[5]      4.25      5.16      1.29      1.21     13.19      2.00    136.07\n",
      "\n",
      "Number of divergences: 0\n"
     ]
    }
   ],
   "source": [
    "def model(action, contact, intention, response=None):\n",
    "    # priors\n",
    "    cutpoints = numpyro.sample('cutpoints', \n",
    "                               dist.TransformedDistribution( \n",
    "                                   dist.Normal(0, 1.5).expand([6]), \n",
    "                                   OrderedTransform()\n",
    "                                   ),\n",
    "                               )\n",
    "    norm = lambda label: numpyro.sample(label, dist.Normal(0,0.5))\n",
    "    bA = norm('bA')\n",
    "    bC = norm('bC')\n",
    "    bI = norm('bI')\n",
    "    bIA = norm('bIA')\n",
    "    bIC = norm('bIC')\n",
    "    # likelihood\n",
    "    BI = bI + bIA*action + bIC*contact\n",
    "    phi = bA*action + bC*contact + BI*intention\n",
    "    response_hat = numpyro.sample('response_hat', dist.OrderedLogistic(phi, cutpoints), obs=response-1)\n",
    "\n",
    "# posterior sampling\n",
    "mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "mcmc.run(rng_key, d.action.values, d.contact.values, d.intention.values, d.response.values)\n",
    "mcmc.print_summary()\n",
    "samples_31 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7LD89ncWDTIb"
   },
   "source": [
    "Now to run the varying intercept model, we need to build a valid individual ID variable. The IDs in the data are long tags, so we can coerce them to integers in many ways. What is important is that the index values go from 1 (_zero in python_) to the number of individuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 188644,
     "status": "ok",
     "timestamp": 1609497862752,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "-kOUxmCuesXa",
    "outputId": "2196fb8a-0503-403c-8ae0-8ef4bd6a4950"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|██████████| 1000/1000 [01:02<00:00, 16.01it/s, 31 steps of size 1.60e-01. acc. prob=0.93]\n",
      "sample: 100%|██████████| 1000/1000 [00:49<00:00, 20.18it/s, 31 steps of size 1.54e-01. acc. prob=0.92] \n",
      "sample: 100%|██████████| 1000/1000 [01:00<00:00, 16.62it/s, 31 steps of size 1.55e-01. acc. prob=0.93]\n",
      "sample: 100%|██████████| 1000/1000 [00:46<00:00, 21.42it/s, 31 steps of size 1.49e-01. acc. prob=0.91]\n"
     ]
    }
   ],
   "source": [
    "def model(action, contact, intention, id, response=None):\n",
    "    # priors\n",
    "    cutpoints = numpyro.sample('cutpoints', \n",
    "                               dist.TransformedDistribution( \n",
    "                                   dist.Normal(0, 1.5).expand([6]), \n",
    "                                   OrderedTransform()\n",
    "                                   ),\n",
    "                               )\n",
    "    norm = lambda label: numpyro.sample(label, dist.Normal(0,0.5))\n",
    "    bA = norm('bA')\n",
    "    bC = norm('bC')\n",
    "    bI = norm('bI')\n",
    "    bIA = norm('bIA')\n",
    "    bIC = norm('bIC')\n",
    "    sigma = numpyro.sample('sigma', dist.Exponential(1))\n",
    "    a = numpyro.sample('a', dist.Normal(0, sigma), sample_shape=d.id.unique().shape)\n",
    "    # likelihood\n",
    "    BI = bI + bIA*action + bIC*contact\n",
    "    phi = a[id] + bA*action + bC*contact + BI*intention\n",
    "    response_hat = numpyro.sample('response_hat', dist.OrderedLogistic(phi, cutpoints), obs=response-1)\n",
    "\n",
    "# posterior sampling\n",
    "mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=500, num_chains=4, chain_method='sequential')\n",
    "mcmc.run(rng_key, d.action.values, d.contact.values, d.intention.values, d.id.values, d.response.values)\n",
    "samples_32 = az.from_numpyro(mcmc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rKfTh509DTIb"
   },
   "source": [
    "We can begin by comparing the posterior distributions. The original coefficients are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 204
    },
    "executionInfo": {
     "elapsed": 526,
     "status": "ok",
     "timestamp": 1609497879494,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "as16uFbXDTIb",
    "outputId": "03bc3c07-da58-4444-fd18-56fe68fc14d7"
   },
   "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_5.5%</th>\n",
       "      <th>hdi_94.5%</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>bA</th>\n",
       "      <td>-0.29</td>\n",
       "      <td>0.33</td>\n",
       "      <td>-0.53</td>\n",
       "      <td>0.28</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.12</td>\n",
       "      <td>7.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>1.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bC</th>\n",
       "      <td>-0.54</td>\n",
       "      <td>0.34</td>\n",
       "      <td>-1.12</td>\n",
       "      <td>-0.27</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.13</td>\n",
       "      <td>7.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>1.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bI</th>\n",
       "      <td>-0.48</td>\n",
       "      <td>0.32</td>\n",
       "      <td>-1.03</td>\n",
       "      <td>-0.24</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.12</td>\n",
       "      <td>7.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>1.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bIA</th>\n",
       "      <td>-0.79</td>\n",
       "      <td>0.62</td>\n",
       "      <td>-1.85</td>\n",
       "      <td>-0.35</td>\n",
       "      <td>0.30</td>\n",
       "      <td>0.23</td>\n",
       "      <td>7.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>1.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bIC</th>\n",
       "      <td>-0.50</td>\n",
       "      <td>1.28</td>\n",
       "      <td>-1.34</td>\n",
       "      <td>1.72</td>\n",
       "      <td>0.63</td>\n",
       "      <td>0.48</td>\n",
       "      <td>7.0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>1.55</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     mean    sd  hdi_5.5%  hdi_94.5%  mcse_mean  mcse_sd  ess_bulk  ess_tail  \\\n",
       "bA  -0.29  0.33     -0.53       0.28       0.16     0.12       7.0      23.0   \n",
       "bC  -0.54  0.34     -1.12      -0.27       0.17     0.13       7.0      20.0   \n",
       "bI  -0.48  0.32     -1.03      -0.24       0.16     0.12       7.0      10.0   \n",
       "bIA -0.79  0.62     -1.85      -0.35       0.30     0.23       7.0      24.0   \n",
       "bIC -0.50  1.28     -1.34       1.72       0.63     0.48       7.0      20.0   \n",
       "\n",
       "     r_hat  \n",
       "bA    1.54  \n",
       "bC    1.53  \n",
       "bI    1.54  \n",
       "bIA   1.54  \n",
       "bIC   1.55  "
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(samples_31, var_names=['bA', 'bC', 'bI', 'bIA', 'bIC'], hdi_prob=.89).round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "p_19DQmWDTIc"
   },
   "source": [
    "And the new ones, having added the individual IDs, are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 235
    },
    "executionInfo": {
     "elapsed": 659,
     "status": "ok",
     "timestamp": 1609497893241,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "Uz64hzbdDTIc",
    "outputId": "8d5ecc98-9227-4600-da73-3a6aa19e4fd8"
   },
   "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_5.5%</th>\n",
       "      <th>hdi_94.5%</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>bA</th>\n",
       "      <td>-0.65</td>\n",
       "      <td>0.06</td>\n",
       "      <td>-0.74</td>\n",
       "      <td>-0.56</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2090.0</td>\n",
       "      <td>1678.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bC</th>\n",
       "      <td>-0.45</td>\n",
       "      <td>0.07</td>\n",
       "      <td>-0.57</td>\n",
       "      <td>-0.35</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2334.0</td>\n",
       "      <td>1674.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bI</th>\n",
       "      <td>-0.39</td>\n",
       "      <td>0.06</td>\n",
       "      <td>-0.49</td>\n",
       "      <td>-0.30</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1694.0</td>\n",
       "      <td>1597.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bIA</th>\n",
       "      <td>-0.56</td>\n",
       "      <td>0.08</td>\n",
       "      <td>-0.69</td>\n",
       "      <td>-0.43</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2010.0</td>\n",
       "      <td>1711.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>bIC</th>\n",
       "      <td>-1.67</td>\n",
       "      <td>0.10</td>\n",
       "      <td>-1.83</td>\n",
       "      <td>-1.50</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2145.0</td>\n",
       "      <td>1470.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>1.92</td>\n",
       "      <td>0.08</td>\n",
       "      <td>1.78</td>\n",
       "      <td>2.04</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3166.0</td>\n",
       "      <td>1528.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mean    sd  hdi_5.5%  hdi_94.5%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "bA    -0.65  0.06     -0.74      -0.56        0.0      0.0    2090.0   \n",
       "bC    -0.45  0.07     -0.57      -0.35        0.0      0.0    2334.0   \n",
       "bI    -0.39  0.06     -0.49      -0.30        0.0      0.0    1694.0   \n",
       "bIA   -0.56  0.08     -0.69      -0.43        0.0      0.0    2010.0   \n",
       "bIC   -1.67  0.10     -1.83      -1.50        0.0      0.0    2145.0   \n",
       "sigma  1.92  0.08      1.78       2.04        0.0      0.0    3166.0   \n",
       "\n",
       "       ess_tail  r_hat  \n",
       "bA       1678.0    1.0  \n",
       "bC       1674.0    1.0  \n",
       "bI       1597.0    1.0  \n",
       "bIA      1711.0    1.0  \n",
       "bIC      1470.0    1.0  \n",
       "sigma    1528.0    1.0  "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(samples_32, var_names=['bA', 'bC', 'bI', 'bIA', 'bIC', 'sigma'], hdi_prob=.89).round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WOIHeaWMDTIc"
   },
   "source": [
    "Everything has gotten more negative. This is because there is a lot of individual variation in average rating, look at the distribution for sigma. That is on the $logit$ scale, so that's a lot of variation on the probability scale. That variation in average rating was hiding some of the effect of the treatments. We get more precision by conditioning on individual.\n",
    "\n",
    "The WAIC comparison can also help show how much variation comes from individual differences in average rating:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 233
    },
    "executionInfo": {
     "elapsed": 5571,
     "status": "ok",
     "timestamp": 1609498189373,
     "user": {
      "displayName": "Andrés Suárez",
      "photoUrl": "https://lh5.googleusercontent.com/-s0kzcIwylzA/AAAAAAAAAAI/AAAAAAAAQXA/v8Sc6WgQy7c/s64/photo.jpg",
      "userId": "06409440331868776168"
     },
     "user_tz": -60
    },
    "id": "DOt1m8HXjZJc",
    "outputId": "8ec8c748-6a0d-41db-a3b1-7c1209893213"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:145: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:1405: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
      "See http://arxiv.org/abs/1507.04544 for details\n",
      "  warnings.warn(\n",
      "/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/arviz/stats/stats.py:1405: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
      "See http://arxiv.org/abs/1507.04544 for details\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rank</th>\n",
       "      <th>waic</th>\n",
       "      <th>p_waic</th>\n",
       "      <th>d_waic</th>\n",
       "      <th>weight</th>\n",
       "      <th>se</th>\n",
       "      <th>dse</th>\n",
       "      <th>warning</th>\n",
       "      <th>waic_scale</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>model_32</th>\n",
       "      <td>0</td>\n",
       "      <td>31056.632797</td>\n",
       "      <td>356.438751</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>179.443929</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>model_31</th>\n",
       "      <td>1</td>\n",
       "      <td>348197.337167</td>\n",
       "      <td>154668.515625</td>\n",
       "      <td>317140.70437</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>2493.716070</td>\n",
       "      <td>2512.3124</td>\n",
       "      <td>True</td>\n",
       "      <td>deviance</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          rank           waic         p_waic        d_waic        weight  \\\n",
       "model_32     0   31056.632797     356.438751       0.00000  1.000000e+00   \n",
       "model_31     1  348197.337167  154668.515625  317140.70437  1.110223e-16   \n",
       "\n",
       "                   se        dse  warning waic_scale  \n",
       "model_32   179.443929     0.0000     True   deviance  \n",
       "model_31  2493.716070  2512.3124     True   deviance  "
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.compare({'model_31': samples_31, 'model_32': samples_32}, ic='waic', scale='deviance')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "95Zju_FSDTId"
   },
   "source": [
    "The WAIC difference is massive. This is consistent with individual variation in average rating being a major effect in this sample. This is all quite typical of likert-scale data, in my experience. Individuals anchor on different points and this adds noise. When we have repeat samples from the same individual, we can condition away some of that noise and get more precise estimates of the treatment effects."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "rethink_stats_numpyro_w08.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}