16 Poisson mixed-effects model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "00051566-8f3b-425c-b091-b3dc323f09e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import arviz as az\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") \n",
    "import scipy.stats as st\n",
    "from patsy import dmatrix, dmatrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "361af0a2-c6cf-4fbd-adc1-75b3e6908fdb",
   "metadata": {},
   "source": [
    "### 16.2 Data generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a609e921-29cd-4386-97fe-810c8714fc7e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAGACAYAAACncLuXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAecElEQVR4nO3df/RcdX3n8eeLH6IICpSQhl8CGqn8kKgptUe3i0oLelqJtuzCVsu2ttQtbk2PW0vYrj/apdo9q3b3WO2i0tJTkdIigWWtlaLUaq0YMPLTlAgIAUyiFkFQUvC9f8z96vD1m+Q7yUzmO/N5Ps6ZM/d+7r0z7893JnnN/dw7d1JVSJKk6bfbuAuQJEm7hqEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1Yqyhn+TJSa5L8qUktyR5e9d+QJKrk9ze3e/ft82qJOuTrEtyyviqlyRpsmSc39NPEuCpVfXtJHsCnwHeCLwa+GZVvTPJucD+VfU7SY4BPgKcCBwM/B3w7Kp6fExdkCRpYox1T796vt3N7tndCjgNuKhrvwhY0U2fBlxSVY9W1Z3AenofACRJ0nbsMe4CkuwOXA88C/jjqvp8ksVVdT9AVd2f5KBu9UOAf+rbfEPXtlUHHnhgHXHEEcMvXJKkBer666//elUtmt0+9tDvhuaXJdkPuDzJcdtYPXM9xA+tlJwNnA1w+OGHs2bNmmGUKknSREjy1bnaF8zZ+1X1AHAtcCqwMckSgO5+U7faBuCwvs0OBe6b47EuqKrlVbV80aIf+qAjSVKTxn32/qJuD58kTwFOBr4MXAmc1a12FnBFN30lcEaSvZIcCSwFrtulRUuSNKHGPby/BLioO66/G3BpVV2V5HPApUleB9wNnA5QVbckuRS4FXgMOMcz9yVJmp+xfmVvV1i+fHl5TF+S1JIk11fV8tntC+aYviRJGi1DX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktSIcV+GdyK94x3vYOPGjeMuY6sWL17MqlWrxl2GJGmBMfR3wMaNG1mxYsW4y9iq1atXj7sESdIC5PC+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqxFhDP8lhST6V5LYktyR5Y9f+tiT3Jlnb3V7Rt82qJOuTrEtyyviqlyRpsuwx5ud/DHhTVd2QZF/g+iRXd8veU1X/s3/lJMcAZwDHAgcDf5fk2VX1+C6tWpKkCTTWPf2qur+qbuimHwJuAw7ZxianAZdU1aNVdSewHjhx9JVKkjT5Fswx/SRHAM8DPt81vSHJjUkuTLJ/13YIcE/fZhvY9ocESZLUWRChn2Qf4DJgZVU9CLwfeCawDLgfeNfMqnNsXnM83tlJ1iRZs3nz5tEULUnShBl76CfZk17gf7iqPgpQVRur6vGq+h7wAX4whL8BOKxv80OB+2Y/ZlVdUFXLq2r5okWLRtsBSZImxLjP3g/wIeC2qnp3X/uSvtVeBdzcTV8JnJFkryRHAkuB63ZVvZIkTbJxn73/IuC1wE1J1nZt5wFnJllGb+j+LuDXAarqliSXArfSO/P/HM/clyRpfsYa+lX1GeY+Tv+xbWxzPnD+yIqSJGlKjf2YviRJ2jUMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNWKsoZ/ksCSfSnJbkluSvLFrPyDJ1Ulu7+7379tmVZL1SdYlOWV81UuSNFnGvaf/GPCmqnoO8ELgnCTHAOcC11TVUuCabp5u2RnAscCpwPuS7D6WyiVJmjBjDf2qur+qbuimHwJuAw4BTgMu6la7CFjRTZ8GXFJVj1bVncB64MRdWrQkSRNq3Hv635fkCOB5wOeBxVV1P/Q+GAAHdasdAtzTt9mGrk2SJG3Hggj9JPsAlwErq+rBba06R1vN8XhnJ1mTZM3mzZuHVaYkSRNt7KGfZE96gf/hqvpo17wxyZJu+RJgU9e+ATisb/NDgftmP2ZVXVBVy6tq+aJFi0ZXvCRJE2TcZ+8H+BBwW1W9u2/RlcBZ3fRZwBV97Wck2SvJkcBS4LpdVa8kSZNsjzE//4uA1wI3JVnbtZ0HvBO4NMnrgLuB0wGq6pYklwK30jvz/5yqenyXVy1J0gQaa+hX1WeY+zg9wMu2ss35wPkjK0qSpCk19mP6kiRp1zD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUiHmHfpJXJtljlMVIkqTRGWRP//eB+5O8N8lPjKogSZI0GvMO/ao6ATgZ+A5wWZJ1SX43yRGjKk6SJA3PQMf0q+pLVfXbwGHAOcDpwFeSfDrJLybxHAFJkhaogY/RJ3km8Jru9j3gLcDdwBuAnwdePcwCJUnScMw79JOcA7wWeBZwKfDaqvqnvuWXAZuGXqEkSRqKQfb0Xw68C7iiqrbMXlhVjyRxL1+SpAVqkND/BeDxqvrXmYYkewK7VdWjAFX1iSHXJ0mShmSQE+8+AbxgVtsLgL8dXjmSJGlUBgn95wKfn9V2HXDC8MqRJEmjMkjoPwAsntW2GHh4aNVIkqSRGST0LwMuTnJckr2THA/8Ob0z+SVJ0gI3SOj/V+A2ekP6DwH/BKwDzhtBXZIkacjmffZ+VX0XOCfJG4ADga9XVY2sMkmSNFQDXZEvydOBo4F9unkAquqTQ69MkiQN1SBX5PuPwB8D3wYe6VtUwFHDLUuSJA3bIHv65wO/UFV/M6piJEnS6AxyIt8e9C7QI0mSJtAgof+HwO/687mSJE2mQYb3fwv4UeDNSb7Rv6CqDh9qVZIkaegGCf3XjKwKSZI0coN8T//vR1mIJEkarXkfn0+yV5Lzk9yR5Ftd2890F+uRJEkL3CAn5b0HOA74RXrfzQe4BfhPwy5KkiQN3yDH9F8FPKuqHk7yPYCqujfJIaMpTZIkDdMge/pbmPUhIcki4Btzr759SS5MsinJzX1tb0tyb5K13e0VfctWJVmfZF2SU3b0eSVJatEgof9XwEVJjgRIsgR4L3DJTjz/nwGnztH+nqpa1t0+1j3fMcAZwLHdNu9LsvtOPLckSU0ZJPTPA+4CbgL2A24H7gPevqNPXlWfBr45z9VPAy6pqker6k5gPXDijj63JEmtmXfoV9WWqlpZVfsAi4F9q+q3qmrLCOp6Q5Ibu+H//bu2Q4B7+tbZ0LX9kCRnJ1mTZM3mzZtHUJ4kSZNnkK/sHTVzA/YFjuybH6b3A88ElgH3A++aKWGOdWuONqrqgqpaXlXLFy1aNOTyJEmaTIOcvb+eXsj2h+9M6A7t2HpVbZyZTvIB4KpudgNwWN+qh9I7vCBJkuZhkOH93apq9+5+N+Bg4ALgtcMsqDtBcMargJkz+68EzuguEnQksBS4bpjPLUnSNBtkT/8JquprSVYC/wxcvCOPkeQjwEnAgUk2AG8FTkqyjN4owl3Ar3fPd0uSS4FbgceAc6rq8R2tX5Kk1uxw6HeOBvbe0Y2r6sw5mj+0jfXPB87f0eeTJKll8w79JP/AE0+c25ved+Z/b9hFSZKk4RtkT/+Ds+YfBr5UVbcPsR5JkjQig/y07kWjLESSJI3WIMP78xrGr6q37Hg5kiRpVAYZ3l8K/DzwBeCrwOH0LoN7GfDdbp05L5YjSZLGb5DQD3BmVV32/Ybk1cDpVfXLQ69MkiQN1SA/uPNyYPWstiuAV/zwqpIkaaEZJPTXA+fMavsN4CvDK0eSJI3KIMP7vwpcnuTNwL30fuHuMeDVoyhMkiQN1yBf2ftikqXAC+ldd/9+4HNV9a+jKk6SJA3PIMP7T1BVnwaelOSpQ6xHkiSNyLxDP8nx9H5c5wP84Pr4/xa4cAR1SZKkIRtkT//9wFuq6seAmSH9vwdePPSqJEnS0A0S+scCf9FNF0BVPQw8ZdhFSZKk4Rsk9O8CXtDfkOREel/lkyRJC9wgX9n7b8D/S/In9E7gWwW8Hvi1kVQmSZKGat57+lV1Fb2r8i2idyz/GcCrq+oTI6pNkiQN0bz29JPsTu/M/WOq6jdGW5IkSRqFee3pV9XjwOPAk0dbjiRJGpVBjun/EXBpkj8ANtD3M7pVdceQ65IkSUO23dBP8qNV9TXgvV3TyfR+ZndGAbuPoDZJkjRE8xne/2eAqtqtqnYDrpyZ7m4GviRJE2A+oZ9Z8/92FIVIkqTRmk/o16z52R8CJEnSBJjPiXx7JHkJPwj73WfNU1WfHEVxkiRpeOYT+pt44i/pfWPWfAFHDbMoSZI0fNsN/ao6YhfUIUmSRmyQH9yRJEkTzNCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRFjDf0kFybZlOTmvrYDklyd5Pbufv++ZauSrE+yLskp46lakqTJNO49/T8DTp3Vdi5wTVUtBa7p5klyDHAGcGy3zfuS7L7rSpUkabKNNfSr6tPAN2c1nwZc1E1fBKzoa7+kqh6tqjuB9cCJu6JOSZKmwbj39OeyuKruB+juD+raDwHu6VtvQ9cmSZLmYSGG/tZkjraac8Xk7CRrkqzZvHnziMuSJGky7DHuAuawMcmSqro/yRJgU9e+ATisb71DgfvmeoCqugC4AGD58uVzfjCYZmvXrmXlypXjLmObFi9ezKpVq8ZdhiQ1ZSGG/pXAWcA7u/sr+tovTvJu4GBgKXDdWCpc4LZs2cKKFSvGXcY2rV69etwlSFJzxhr6ST4CnAQcmGQD8FZ6YX9pktcBdwOnA1TVLUkuBW4FHgPOqarHx1K4JEkTaKyhX1VnbmXRy7ay/vnA+aOrSJKk6TVJJ/JJkqSdYOhLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRe4y7gK1JchfwEPA48FhVLU9yAPCXwBHAXcC/q6p/GVeNkiRNkoW+p/+SqlpWVcu7+XOBa6pqKXBNNy9JkuZhoYf+bKcBF3XTFwErxleKJEmTZSGHfgGfSHJ9krO7tsVVdT9Ad3/Q2KqTJGnCLNhj+sCLquq+JAcBVyf58nw37D4knA1w+OGHj6o+SZImyoLd06+q+7r7TcDlwInAxiRLALr7TVvZ9oKqWl5VyxctWrSrSpYkaUFbkKGf5KlJ9p2ZBn4GuBm4EjirW+0s4IrxVChJ0uRZqMP7i4HLk0Cvxour6uNJvgBcmuR1wN3A6WOsUZKkibIgQ7+q7gBOmKP9G8DLdn1FGra1a9eycuXKcZexVYsXL2bVqlXjLkOShmpBhr6m35YtW1ixYsW4y9iq1atXj7sESRq6BXlMX5IkDZ+hL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIjDH1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIasce4C5AWorVr17Jy5cpxl7FNixcvZtWqVeMuQ9IEMfSlOWzZsoUVK1aMu4xtWr169bhLkDRhHN6XJKkRhr4kSY0w9CVJaoShL0lSIzyRT5pQfsNA0qAMfWlC+Q0DSYNyeF+SpEYY+pIkNcLQlySpEYa+JEmNMPQlSWqEoS9JUiP8yp6kkVno1xLwOgJqjaEvaWQW+rUEvI6AWuPwviRJjZjI0E9yapJ1SdYnOXfc9UiSNAkmbng/ye7AHwM/DWwAvpDkyqq6dbyVSZo0C/2cA1j45x284x3vYOPGjeMuY5tuv/12li5dOu4ytmpXvsYTF/rAicD6qroDIMklwGmAoS9pIAv9nANY+OcdbNy4ccH/Dc8777wFXeOufI0ncXj/EOCevvkNXZskSdqGVNW4axhIktOBU6rqV7v51wInVtV/7lvnbODsbvZoYN2QyzgQ+PqQH3OhaqmvYH+nWUt9Bfs7zebT12dU1aLZjZM4vL8BOKxv/lDgvv4VquoC4IJRFZBkTVUtH9XjLyQt9RXs7zRrqa9gf6fZzvR1Eof3vwAsTXJkkicBZwBXjrkmSZIWvInb06+qx5K8AfhbYHfgwqq6ZcxlSZK04E1c6ANU1ceAj42xhJEdOliAWuor2N9p1lJfwf5Osx3u68SdyCdJknbMJB7TlyRJO8DQH8C0X/43yYVJNiW5ua/tgCRXJ7m9u99/nDUOS5LDknwqyW1Jbknyxq59Wvv75CTXJflS19+3d+1T2V/oXb0zyReTXNXNT3Nf70pyU5K1SdZ0bdPc3/2S/HWSL3f/hn9yWvub5OjudZ25PZhk5Y7219Cfp77L/74cOAY4M8kx461q6P4MOHVW27nANVW1FLimm58GjwFvqqrnAC8Ezulez2nt76PAS6vqBGAZcGqSFzK9/QV4I3Bb3/w09xXgJVW1rO+rXNPc3/8FfLyqfgw4gd7rPJX9rap13eu6DHgB8AhwOTva36ryNo8b8JPA3/bNrwJWjbuuEfTzCODmvvl1wJJuegmwbtw1jqjfV9D7PYep7y+wN3AD8BPT2l961++4BngpcFXXNpV97fpzF3DgrLap7C/wNOBOunPSpr2/s/r4M8Bnd6a/7unPX6uX/11cVfcDdPcHjbmeoUtyBPA84PNMcX+74e61wCbg6qqa5v7+EfBm4Ht9bdPaV4ACPpHk+u6KpDC9/T0K2Az8aXf45oNJnsr09rffGcBHuukd6q+hP3+Zo82vPky4JPsAlwErq+rBcdczSlX1ePWGCA8FTkxy3JhLGokkPwtsqqrrx13LLvSiqno+vcOP5yT5qXEXNEJ7AM8H3l9VzwMeZkqG8reluxjdK4G/2pnHMfTnb7uX/51SG5MsAejuN425nqFJsie9wP9wVX20a57a/s6oqgeAa+mdvzGN/X0R8MokdwGXAC9N8hdMZ18BqKr7uvtN9I73nsj09ncDsKEbqQL4a3ofAqa1vzNeDtxQVTO/Y7xD/TX056/Vy/9eCZzVTZ9F79j3xEsS4EPAbVX17r5F09rfRUn266afApwMfJkp7G9VraqqQ6vqCHr/Tj9ZVa9hCvsKkOSpSfadmaZ33PdmprS/VfU14J4kR3dNL6P30+pT2d8+Z/KDoX3Ywf56cZ4BJHkFvWOFM5f/PX+8FQ1Xko8AJ9H7BaeNwFuB1cClwOHA3cDpVfXNMZU4NEleDPwDcBM/OO57Hr3j+tPY3+cCF9F77+4GXFpVv5fkR5jC/s5IchLwX6rqZ6e1r0mOord3D72h74ur6vxp7S9AkmXAB4EnAXcAv0z3vmY6+7s3vXPKjqqqb3VtO/T6GvqSJDXC4X1Jkhph6EuS1AhDX5KkRhj6kiQ1wtCXJKkRhr4kSY0w9CVJaoShL0lSIwx9SZIaYehLktQIQ1+SpEYY+pIkNcLQl3ZQklu6X3Eb9fMcneSLSR5K8pujfr7uOd+S5L1DeJx7kjx/GDXNesznDemxBupnkruSnDzgc3z/fbIj20vDZOired1/xN9J8u0kG5P8aZJ9trddVR1bVdfughLfDFxbVftW1f+evTDJ/kmqq/+RJF9N8rqdfM5j6P3s8A5Lsh9wMPDlnayl/zH3B5YAtw3pIXe6n9uzC98n0nYZ+lLPz1XVPsDzgR8HfndHHyjJHkOrqucZwC3bWL4M+HpV7VNVewOrgP+T5MCdeM5jgRt3YnuA44G7q+qRnXyc2Y+5vqq+O6THG0Y/5zSC94G00wx9qU9V3Qv8DXAcQJJzk3ylG1q/NcmrZtbtH6rtpn8nyY3Aw0n26Obv7bZdl+Rlcz1nkuckuTbJA91Q8Cv7ln0SeAnw3m5P/tlzPMQy4Ia++b8Hdgf273ucX+vq/1aSv0lyUN+y3ZKsSnJ3kvuSnAE8C7h5wD/fbMcD65K8O8m/JLk9yYu75zwoyZXdyMqDSf5vkqfN+rucmeT6ruavdEPkz52pK8neSS5O8tGZkZlh9XM7r92yJDd2z/GXSZ7ct91c74M5h/ST/FiSO7s6SHJwksuSbO7af3PW+vN6P0nbVFXevDV9A+4CTu6mD6O3V/373fzp9IaodwP+PfAwsGSO7e4C1nbbPwU4GrgHOLhbfgTwzDmee09gPXAe8CTgpcBDwNF961wL/Oo26v9z4A+66f26+TVAurbzgOvpBdyTgA8CF/Rt/zbgs/RGFJ7eTd8xhL/rnwD/Aqzo+vlW4LZu2bOAnwb2Ag4APgf8dt+2b+peh+d3f/vju7/hBd3jHEnvg85bh93Pbb123et8XfeeOIDeYYbXz3ovff99sJX3ycldv+4GfrZr362r/S1d7UcBdwCnbK8mb94GuY29AG/exn3r/iP+NvAA8FXgfTP/Yc+x7lrgtL7t+v8z/5W+9Z4FbOr+g99zG8/9b4CvAbv1tX0EeFvf/LVsO/RvpPdh5EGggI8DB3bLDgIeAZ7dt/5PAl/sphd1fX9m3/LzgNVD+Lt+FnhL3/whwPeAPeZY9/f4wQeXRV1fTphjvc/R+zBx58zrMOx+buu1617n1/TN/w/gT2Yt/5U5tul/n7wd2AC8pG+dn6B3KKR/u1XAnw7yfvLmbXs3h/elnhVVtV9VPaOqfqOqvgOQ5JeSrO2G3h+gN+y/tWPl98xMVNV6YCW9vctNSS5JcvAc2xwM3FNV3+tr+yq9gNyuJHsBzwGeW1VPA34BeCHwr90qL6O353hdXx8+Dnyrb/ltVfWVvoddzKyT27rDD7WV22e2Ut5xwF/3zR8IfKuqHktyepLPJtnU1XQu8M/deicDN1XVl2bVkO4xX0UvaK/oWzyUfsK8Xruv9U0/Asw+6fMetu31wD9W1af62p4BHDxTe1f/eV2Ng7yfpG0y9KWtSPIM4APAG4Afqar96B3/zVY2qSfMVF1cVS+m9x96AX84xzb3AYcl6f+3eDhw7zzLPA54lN5QMFV1Gb1h45/vlh8AXN59oJm5Pb2qTuqWH0hvDxKAJHvSG45/wsltVXVSVWUrtxfPLqr72z0N2NzX/GrgqiQvpfe3WEnvQ89MDWv7an5gjr4e2d2fDLwpyfK+ZUPpZ19/5/PabU1tZ/nrgcOTvKev7R7gzln171tVrxhSTRJg6Evb8lR6/7luBkjyy3Qn+G1Pet+tf2m3J/5d4DvA43Os+nl6Q/NvTrJnd7LazwGXzLPG5wE3V1V/0HwMmDkZ8AbgJem+K5/kaUlO6/aaAdYBL07y7CRPB95P70PHzn6N7XjgMeA/dCfQvQL4dXrD+CfQC7kv0zvZ8EJ6w/O3dtt+savphPQsTfIceifx3VhVNwFnA5cnWTLsfg7w2u2oh4BTgZ9K8s6u7Trgwe5kvack2T3JcUl+fBfVpEYY+tJWVNWtwLvoHUfeSC/IPjvPzfcC3gl8nd5w8EH0hmtnP8cWegH98m7d9wG/VFXz/W77Mn54b/XjwE8neXJVfY5e0F6W5Nv0gvXUmQ8JVXU1vQ8Ya4Av0PuA813g9nk+/9YcT++EwhfRO5nvbfSOwd8OfJjeiX1fA67qnuvW7m9BVf0j8N+7ZQ8Bl9M7OfL4mb5W1Wp6J/WtHkE/5/Xa7YyqeoDeiYwvT/L7VfU4vQ97y+idr/B1eiciPn1X1aQ25Ik7CJIkaVq5py9JUiMMfUmSGmHoS5LUCENfkqRGGPqSJDXC0JckqRGGviRJjTD0JUlqhKEvSVIj/j9H+oOip/qtDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(123456)\n",
    "n_groups = 16\n",
    "n_years = 30\n",
    "n = n_groups * n_years\n",
    "pop = np.repeat(range(1, n_groups + 1), n_years)\n",
    "original_year = np.tile(range(1, n_years + 1), n_groups)\n",
    "year = (original_year - 1) / 29\n",
    "data_ds = pd.DataFrame({'pop': pd.Series(pop, dtype = 'category'), 'year': year})\n",
    "Xmat = dmatrix(\"~ pop * year - 1 - year\", data_ds, return_type = 'dataframe')\n",
    "\n",
    "intercept_mean = 3                       # Choose values for the hyperparams\n",
    "intercept_sd = 1\n",
    "slope_mean = -2\n",
    "slope_sd = 0.6\n",
    "intercept_effects = st.norm.rvs(intercept_mean, intercept_sd, size = n_groups)\n",
    "slope_effects = st.norm.rvs(slope_mean, slope_sd, size = n_groups)\n",
    "all_effects = np.concatenate((intercept_effects, slope_effects))\n",
    "\n",
    "lin_pred = Xmat.to_numpy() @ all_effects\n",
    "λ = np.exp(lin_pred)\n",
    "C = np.random.poisson(λ)\n",
    "\n",
    "data = pd.DataFrame({'C': C, 'pop': pop, 'or_year': original_year})\n",
    "\n",
    "plt.figure(figsize = (8, 6))\n",
    "plt.hist(C, color= 'gray', edgecolor='black', alpha = 0.5)\n",
    "plt.xlabel('Pairs of $Red-backed$ shrikes', fontsize = 12, labelpad = 10)\n",
    "plt.ylabel('Frequency', fontsize = 12, labelpad = 10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "424979ec-2e59-47df-a0ee-109f3a288796",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8YAAAMvCAYAAAAQ2smLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAB4/klEQVR4nO3de7xcd1no/8/TNLGkDZLaBEPTEqlF9KgFEqkekLO1xAv3niOIKFZEq+cnCopHC945XqqiWO9WAcsBrD1opFZBSz3bY730kA1UhCIhNaWhoSk0QNJQkqbP74+ZXXbTPbNnr1lrzXdmfd6v137tmbX2zHpm9vqumWc93+93RWYiSZIkSVJXnTLpACRJkiRJmiQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMVYuI2BcRn4mIIxFxZ0S8ISLOaHibWyLi2oi4IyIyIrY1uT1pWkyoPT4jIm6MiE9GxMci4g8jYkOT25SmxYTa5NdFxPv6bfITEbErIs5ucpvSNJhEezxp+2/of2/94ra2qdGYGKtOz8rMM4AnAl8F/GTD27sfeAfw3xrejjSN2m6Pnw/8PPAo4EuBrcCvNrxNaZq03SY/AHxjZj6CXrvcA/xew9uUpkXb7RGAiHgKcF4b29LqmRirdpn5UeDtwJcDRMSzI+L9/bPW8xHxpYt/2z9r98qI+EBEHOqfRTttxO3cmZm/C7yrkRcizYAW2+NbMvMdmXk0Mw8Bfwg8uYnXJE2zlj8j71iy6ARghUpaoq322H/8qcBvAS+t+3WoHibGql1EnAM8HXhPRDwW+BPg5cAm4K+Bv4yIdUse8u3AN9I7g/ZYlpy16x+YntJS6NLMmWB7fCrw/rFfgDRj2myTEXFuRHwS+Azwo8Cv1PpipCnX8mfkDwP/NzP/tdYXodpEZk46Bs2AiNgHnAXcB3wK+CvgFfQ+iL8iM5/f/7tTgNuBb8/M+f7jLs/M3++vfzrwW5k5cjeT/hm448AXZea+ul6TNK0m2R77j9sJXANcmJkfquVFSVOsgDZ5JvC9wN9n5r/U8qKkKTWJ9thPwP8PsD0zPxURCZyfmR+u+/WpulMnHYBmynMz851LF0TEo4DbFu9n5v0RcTuwdAKQ25fcvo3eWChJ45lIe4yIrwbeAnyLSbH0IBP7jMzMuyPiKuDmiDg7M+9b7XNIM6bt9vgbwKsz81PVwlUb7Eqtpt0BPHrxTkQEcA7w0SV/c86S2+f2HyOpfo22x4h4AnAt8N2ZecN4oUqd0OZn5KnAZuDhFR8vzbom2+NFwK9G76oNH+sv++eIeOEY8apmJsZq2jXAMyLioohYS6+rymeBf1ryNz8QEVv7Xb1eBfzpqE/en/Tg8/p3P281kyBIHdRYe4yIL6c3S/wPZuZf1hy3NKuabJP/NSK+JCJOiYhNwK8D78nMu2t+DdKsaPI762OBC4DH938AngXsqiFu1cTEWI3KzH8HvoPeLHwfp3cQeFZmHlvyZ28B/ha4tf/z84sroneNua8dsonPAEf6tz/Yvy9pGQ23x1fQm6zkdf2/OxIRTr4lDdFwmzyb3smqw8D76F3i8OK6X4M0K5psj5l5MDM/tvjTX/zxzPR7a0GcfEsT1Z/I4HtOHuchqX22R6kstkmpHLbH2WfFWJIkSZLUaSbGkiRJkqROsyu1JEmSJKnTrBhLkiRJkjrNxFiSJEmS1GmnTjqAJpx11lm5bdu2SYchdcbCwsLHM3PTcutsj1L7BrVJ26PUPj8jpbIMapMzmRhv27aN3bt3TzoMqTMi4rZB62yPUvsGtUnbo9Q+PyOlsgxqk3alliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdOKSowj4hER8daI+GBE3BIRXxMRZ0bE9RGxp/9746TjlCRJkiTNjqISY+AK4B2Z+TjgAuAW4DLghsw8H7ihf1+SJEmSpFoUkxhHxMOBpwKvA8jMY5n5SeA5wFX9P7sKeO4k4pMkSZIkzaZiEmPgMcBdwBsi4j0R8UcRcTrwyMw8AND/vXmSQUqSJEmSZktJifGpwBOB38vMJwD3sIpu0xFxaUTsjojdd911V1MxShqB7VEqh+1RKottUipTSYnxfmB/Zt7Uv/9WeonynRGxBaD/++ByD87MKzNzR2bu2LRpUysBS1qe7VEqh+1RKottUipTMYlxZn4MuD0ivqS/6CLgA8C1wCX9ZZcAb5tAeJIkSZKkGXXqpAM4yQ8Cb46IdcCtwIvpJe/XRMRLgI8Az5tgfJIkSZKkGVNUYpyZ7wV2LLPqopZDkSRJkiR1RDFdqSVJkiRJmgQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqNBNjSZIkSVKnmRhLkiRJkjrNxFiSJEmS1GkmxpIkSZKkTjMxliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeq0ohLjiFgTEe+JiOv698+MiOsjYk//98ZJxyhJkiRJmi1FJcbAy4Bblty/DLghM88HbujflyRJkiSpNsUkxhGxFXgG8EdLFj8HuKp/+yrguS2HJUmSJEmaccUkxsBvAD8G3L9k2SMz8wBA//fmCcQlSZIkSZphRSTGEfFM4GBmLozxHJdGxO6I2H3XXXfVGJ2k1bI9SuWwPUplsU1KZSoiMQaeDDw7IvYBVwNfHxFvAu6MiC0A/d8HBz1BZl6ZmTsyc8emTZvaiFnSALZHqRy2R6kstkmpTEUkxpn5yszcmpnbgBcAf5eZ3wFcC1zS/7NLgLdNKERJkiRJ0owqIjEe4nJgZ0TsAXb270uSJEmSVJtTJx3AyTJzHpjv3/4EcNEk45EkSZIkzbbSK8aSJEmSJDXKxFiSJEmS1GkmxpIkSZKkTjMxliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE47ddIBTMLNe+5lfuEoBw+dYPPGNcxtX88F55826bAkSZIkSRPQucT45j33cv1N93Dx3Aa2bVnLvgPH2TV/GMDkWJIkSZI6qHNdqecXjnLx3AbO27qONWuC87au4+K5DcwvHJ10aJIkSZKkCehcYnzw0Am2bVn7oGXbtqzl4KETE4pIkiRJkjRJnUuMN29cw74Dxx+0bN+B42zeuGZCEUmSJEmSJqlzifHc9vXsmj/M3v3HOHEi2bv/GLvmDzO3ff2kQ5MkSZIkTUDnJt9anGDruhuPPDAr9c4LT39guTNWS5IkSVK3dC4xhl5yvFyy64zVkiRJktQ9netKPYwzVkuSJElS95gYL+GM1ZIkSZLUPSbGSzhjtSRJkiR1j4nxEs5YLUmSJEnd08nJtwZZacbqKpzlWpIkSZLKZmJ8kkEzVlfhLNeSJEmSVD4T45osVxleOss18MAs19fdeMTEWJIkSZIKYWJcg0GV4Tvuus9ZriVJkiSpcCbGNRhUGX7tW+5m34HjDyyHcme5diy0JEmSpK4yMa7BoOsfxymwa/7wQyrJOy88fUKRLs+x0JIkSZK6zMS4BovXPz65MvzFW9cxt319rbNcQ/3VXcdCS5IkSeoyE+MaLF7/eLnKcJ2zXEMz1d1BFW/HQkuSJEnqAhPjGjRx/eNBmqjuDqp4lzgWWpIkSZLqZmJck7orw4M0Ud0dVvGWJEmSpFlXTGIcEecAbwS+ELgfuDIzr4iIM4E/BbYB+4DnZ+ahScU5aU1Ud9useEuSJElSaYpJjIH7gFdk5rsjYgOwEBHXA98F3JCZl0fEZcBlwI9PIsASLmnUVHW3rYq3JEmSJJWmmMQ4Mw8AB/q3D0fELcDZwHOAuf6fXQXMM4HEuJRLGlndlSRJkqR6FZMYLxUR24AnADcBj+wnzWTmgYjYPImYSrqkkdVdSZIkSarPKZMO4GQRcQbwZ8DLM/PTq3jcpRGxOyJ233XXXbXH5SWNpNE13R4ljc72KJXFNimVqajEOCLW0kuK35yZf95ffGdEbOmv3wIcXO6xmXllZu7IzB2bNm2qPbbFSa+WKvGSRjfvuZcrrr6bn/i9u7ji6ru5ec+9kw5JHdR0e5Q0OtujVBbbpFSmYhLjiAjgdcAtmfnrS1ZdC1zSv30J8La2Y4PPTXq1d/8xTpxI9u4/xq75w8xtXz+JcJa1OA76mU85g1dfehbPfMoZXH/TPSbHkiRJkjRESWOMnwy8CHhfRLy3v+xVwOXANRHxEuAjwPMmEdw0THpV0jhoSZIkSZoWxSTGmXkjEANWX9RmLIOUPumV46AlSZIkafWKSYw1vsVx0IsVY5jsOOgSrvssSZIkSSsxMZ4hi+OgT77W8s4LT289llKu+yxJkiRJKzExnqC6K6pNjYOuEqfjnTUqexZIkiRp0kyMJ6Spimrd46Crxul4Z43CngWSJEkqgYnxhExLRbVqnE2Nd7a6OFuaagfuJ5IkSVoNE+MJmZaKatU4mxjvbHVx9jTRDtxPJEmStFomxhOyUkW1lIpX1cpvE+Odp6XKrtGN07NgUBtZaT+p2rZKaZOSJEmqn4nxhAyrqJZU8Rqn8lv3eOdpqbJrdFX3r2FtZNh+UrVtldQmJUmSVD8T4wkZVlG94uq7i6mMrlT5bbOKVtp1mjW+qvvXsKrwsP2kajW5y70VrJRLkqQuMDGeoEEV1dIqo4PibLuKVtJ1mlWfKvvXsDby/KdtGLifXPPOw5WqyaW1ybZYKZckSV1hYlygaamMtl1Fa+o6zSpT1arwsP1kfuFopWrytLTJunW5Ui5JkrrFxLhATVVG6+4SOYkqWt3jllWuqlVhGLyfDGtbw6rJK21vVnW1Ui5JkrrHxLhATVRGm+gS2dUqmtpRtSo8TNVqcld7K9jGJUlSV5gYF6ruymgTXSId81s/Jzr6nJX2r6ptpEo1eZztTTPbuCRJ6goT445ooktkV6toTXGiowdre/9yf34o3xNJktQVJsYd0VSXyC5W0RbVXd11oqOHanv/moX9ue79chbeE0mSpJWYGHeEXSLr1UR114mONC57HUiSJFVjYtwRdomsZlD1rYnqrhMdTa8mxoYPe84298u2Oc5ekrrDY75KYmLcIXaJXJ1h1bcmqrtW9adTE1XaYc8JtLpftsmKtyR1h8d8lcbEWJ1XpfrWRHXXqv50aqJKO+w5gUr7Zdtn5atsbxIVb6sVkjQZs9DLSbPFxFidVrUq/PynbWikumtVf/o0UaVd6TlXu18+5uy1rZ6Vr1oFaLvibbVCkiZn2ns5afaYGKs4bY7XrFoVXqm6axWqO1bqPVBlrPBKz7na/bLts/Irba/q63YmeEkaTdXjZZvfX6ZlbhW/03WHibGK0vZ4zXGqwoOqu1ahumXY2PCqY4VXGm++2v3ymncebvWs/LB2New9qfpeOhO8JH1O1eNl299fpmFuFb/TdYuJsSaiiVl1264KD2IVqluG7SdXXH13pbHCL3vBmQOfc9Fq9sumzspXqfwOax/DXvew93JYFXqYaalWlMTKiVS+qt9DmvgONsw0zK3id7puMTFW65qY7bntqvAwVqG6Z9B+UnWs8LDnXGndcpo4K1+18rtS9brKe1n1jP40VCtKYuVEmg5Vv4c08R1slOS45OOH3+m6xcRYrWtitue2q8LDNFGFqjJOdaV1al7VscJ1a2I/r1r5nV84Wul1V61CD3uNk6hWTHObXGm29Gl9XYumYUymuqeJeSqqPm6QWa6q2rOoW0yM1bomZntuuyo8TN1VqKrjVIetm/YPqmkxzljhutW9n690Fn3Q9qq2j3Gq0MO0Wa2Y9orroP/5h/cf4/h9ObWvC6ZnTKa6pYl5Kqr29hlmlquq9izqFhNjta6JCm7bVeFh6p6xuuo1bYetW2lMphWQeoyy75U8tmqYqmfRq7bHYY+rWoUeR9vXaS6hTQ76n+f9ZR1r2vzfTOJ/WsK+MIppibNkVT//q87XMMr8Fstp6ooCJexD9izqFhNjtW6ls29VKjhNPOc46pyxepxxqlXGZIKV5jrVOVa4JOOcRa/6uuuuQlfV9nWaS6lKDnqf162NYo41bf9v2v6flrIvrGRa4ixdE/NUVO3tM0wTVxQoaR+yZ1F3mBirdU2cfZuGmQ2hWnVhnHGqVcZkwvDqj2bPtM8mWncvjUV1z55ftarS9vi9QXEMep+HVezbPtY0dQ3tQZoak1l1X2i70tRmnKX0LmgqluU0MU/FOGNmV3tsWKlC3XZvjCZ6r9T9nNPes2jamRhrIpo4+zYN1bcq1YVxxqlWHZM5q2OF9FCzMptonb00Vnpc1Sph1apKm+P3Vnq/Br3PpRxrqs5e3sTY9ybibGJ29qrajLOJ1zYNlcwm5qmous9WPTaU0htj7/5j3PrR47X2XmmiR8y09yyadibGUouqnKkdd5xqlTGZzsDYHbM8myg0M3a0auW3alWlzVlRq7xfw17XOMeaumfjrTqT+rDtVe2tUDXOcWZnr7ua1ESci89bZwW67qp2m8fMcT//qz7ncprqKVP1cav9/7z2LXfzwy88s9beK030iGmqZ1FJ88qUXNk2MZZaVPVM7bDKXJV1Jc2WrMma5dlEoZlqxbCZ7puoqlSdrb+Kqu9X3ceaJmbjrXoN7ar/0yau9T1sXxj2uCaqSVX32UFxDpvdvGoFGgbvK01U7ZtQ9fO/6nMO0kRPmaqPq/L/OXLv/bX3XhlnDHhbr7u0eWVKr2xPRWIcEd8EXAGsAf4oMy+fcEhSJaWMy2ziLLSmU5vVyEloolrRxHi6qturW937Q9VjzbD3EarNxlt19vImeh1UjbNqdb6JKmfVfXZQnMNmNx+nAl33c876MXOQqq+76vGr7h42Z5x2Suvjtet+v6r2LCppXpnSe6kVnxhHxBrgd4CdwH7gXRFxbWZ+YLKRSdWUMi6zibPQmj6zfo3GpsaO1j2erpSZ9ZvYH6oca5qYjbfta7Su9LiqcVZ53ErV8irdLOuOc9js5lUr5Yu363zONntwlKSkKxFU6a3w9V/V/njtYeuGtbm6exaNM9dD1S7Yg9atdFyse3urVXxiDDwJ+HBm3goQEVcDzwFMjCVpTKX0YmhKE9WKYdquqtStlDjarsSME0vdj2tinx1WTR63m2VdcTZRKV+8XedzltJG2lbS667aW+HmPfe2Pl57UBxVuhPX3UtjpTZS9dgwbN2w19DE9la7f0ZmruoBbYuIbwG+KTO/p3//RcCFmfnSQY/ZsWNH7t69u60Qpc6LiIXM3LHcOtujumTQB3vbXyAHtclpaY/D3kdY/ktQU+9x1f9pKfvCSrHMLxzlmU8540FfVPfuP/ZAN8tB6xa7g7cRY9X3GartK3X/7/yMrFdJbauKK66+u1K7aruNVD02DFs3t319q9sb9H4OapPTUDGOZZY9JJuPiEuBSwHOPffcpmOSNITtUV1VUlVl0TS2x5LmQWi710EThsVSyuX7mnyfS/vfTWObLElJbauKcSY5hPbayDjHhkHrmjoW1XWcmobEeD9wzpL7W4E7Tv6jzLwSuBJ6Z9/aCU3ScmyP6rLSxuhPa3ssaR6EusdITsKgWJrotl53jOM8rsT/3bS2yZKU1LZWa5wJ3NpsI+McG4ata3t7q3HKqh/RvncB50fEF0XEOuAFwLUTjkmSJGnqLU4gtHf/MU6cSPbuP8au+cPMbV8/dJ2kaqalXVU9NlR9fW1vbznFV4wz876IeCnwN/Qu1/T6zHz/hMOSJEmaeiV1W5e6YFq6go97bGh7WEQd72fxk29V4UQGUrucWEQqy7RPviXNEj8jpbIMapPT0JVakiRJkqTGmBhLkiRJkjptJrtSR8RdwG39u2cBH59gOItKiQOMZRBjeahR43h0Zm5abkWh7RGMZTmlxAHGMshYbfKk9ria52taKXGAsQxSSiylxAF+RrallFhKiQOMZZDxPiNnMTFeKiJ2DxrX0cU4wFgGMZbm4yjldYGxlBwHGMsgs9omS4kDjGWQUmIpJQ6Y3fYIxlJyHGAsg4wbi12pJUmSJEmdZmIsSZIkSeq0LiTGV046gL5S4gBjGcRYHqruOEp5XWAsyyklDjCWQWa1TZYSBxjLIKXEUkocMLvtEYxlOaXEAcYyyFixzPwYY0mSJEmShulCxViSJEmSpIFMjCVJkiRJnWZiLEmSJEnqNBNjSZIkSVKnmRhLkiRJkjrNxFiSJEmS1GkmxpIkSZKkTjMxliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYqxaRMS+iPhMRByJiDsj4g0RcUbD25yLiPv721z8uaTJbUrTYBLtsb/dTRHxloj4ZEQciog3N71NaRpM6DPyVSd9Pn6m/5l5VpPblUo3wc/IH4yI/4iIT0fE7oh4StPb1OqYGKtOz8rMM4AnAl8F/GQL27wjM89Y8nNVC9uUpsEk2uOfAx8DHg1sBl7TwjaladFqm8zMX1z6+Qj8MjCfmR9vcrvSlGi1PUbEhcDlwLcAnw+8DtgVEWua3K5Wx8RYtcvMjwJvB74cICKeHRHv71eR5iPiSxf/tn/W7pUR8YF+hekNEXHapGKXZk1b7TEivgE4B/gfmfmpzDyeme9p4jVJ02wSn5EREcCLAE8eS0u02B63Ae/PzIXMTOCNwFn0TiKrECbGql1EnAM8HXhPRDwW+BPg5cAm4K+Bv4yIdUse8u3ANwLnAY9lyVm7/oFpWFeTzf1uMP8REa+NiNPrfTXSdGuxPX418O/AVRHxiYh4V0T8l7pfjzTtWv6MXPS1wCOBP6vjNUizosX2+HZgTURc2K8SfzfwXnq9rFQIE2PV6S8i4pPAjcDfA78IfCvwV5l5fWYep9e18mHAf17yuN/OzNsz827gF4BvW1yRmY/IzBsHbO+DwOOBLcDXA9uBX6/1FUnTq+32uBX4BuD/AF8I/BrwNsczSg9ou00udQnw1sw8Us9LkaZe2+3xML0TUzcCnwV+Bri0Xz1WIU6ddACaKc/NzHcuXRARjwJuW7yfmfdHxO3A2Uv+7PYlt28DHjXKxjLzY3zuTNt/RMSPAX8FfF+F2KVZ02p7BD4D7MvM1/XvXx0RPwE8GXjbaoOXZlDbbXJxGw8Dngc8Z9URS7Or7fb4PfSqxP8J+DC9E8nXRcQTMvOOCvGrAVaM1bQ76E3EAzwwzukc4KNL/uacJbfP7T+migSi4mOlLmiyPf4rvTYoaXRtfEb+V+BuYL5aiFJnNNkeLwD+MjM/lJn3Z+Y7gAM8uBqtCTMxVtOuAZ4RERdFxFrgFfS6kPzTkr/5gYjYGhFnAq8C/nSUJ47e5ZrOjZ5z6M32Z2VKGqyx9gjsAjZGxCURsSYivoXeWfZ/rDF+adY02SYXXQK80S6b0oqabI/v6j/3Y/rfW3fSG6P8bzXGrzGZGKtRmfnvwHcAvwV8HHgWvSnyjy35s7cAfwvc2v/5+cUV0bvG3NcOePonAv8M3EPvoPVvwA/V/RqkWdFke+yPt3o28KPAp4DLgOd4aRhpsIY/I4mIs+nNwfHG+qOXZkvD7fGNwNX0em58GvhN4Psy84M1vwyNITyBqEmKiH3A95w8zkNS+2yPUllsk1I5bI+zz4qxJEmSJKnTTIwlSZIkSZ1mV2pJkiRJUqdZMZYkSZIkdZqJsSRJkiSp006ddABNOOuss3Lbtm2TDkPqjIWFhY9n5qbl1tkepfYNapO2R6l9fkZKZRnUJmcyMd62bRu7d++edBhSZ0TEbYPW2R6l9g1qk7ZHqX1+RkplGdQm7UotSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqNBNjSZIkSVKnmRhLkiRJkjrNxFiSJEmS1GkmxpIkSZKkTjt10gGU5uY99zK/cJSDh06weeMa5rav54LzT5t0WJIkSZKkhpgYL3Hznnu5/qZ7uHhuA9u2rGXfgePsmj8MYHIsSZIkSTPKrtRLzC8c5eK5DZy3dR1r1gTnbV3HxXMbmF84OunQJEmSJEkNMTFe4uChE2zbsvZBy7ZtWcvBQycmFJEkSZIkqWkmxkts3riGfQeOP2jZvgPH2bxxzYQikiRJkiQ1rajEOCIeERFvjYgPRsQtEfE1EXFmRFwfEXv6vzc2tf257evZNX+YvfuPceJEsnf/MXbNH2Zu+/qmNilJkiRJmrDSJt+6AnhHZn5LRKwD1gOvAm7IzMsj4jLgMuDHm9j44gRb19145IFZqXdeeLoTb0mSJEnSDCsmMY6IhwNPBb4LIDOPAcci4jnAXP/PrgLmaSgxhl5ybCIsSZIkSd1RTGIMPAa4C3hDRFwALAAvAx6ZmQcAMvNARGwed0NVr1XsNY4lSZIkafaUlBifCjwR+MHMvCkirqDXbXokEXEpcCnAueeeO/Dvql6r2GscS6MbtT1Kap7tUSqLbVIqU0mTb+0H9mfmTf37b6WXKN8ZEVsA+r8PLvfgzLwyM3dk5o5NmzYN3EjVaxV7jWNpdKO2R0nNsz1KZbFNSmUqJjHOzI8Bt0fEl/QXXQR8ALgWuKS/7BLgbeNsp+q1ir3GsSRJkiTNppK6UgP8IPDm/ozUtwIvppe8XxMRLwE+AjxvnA0sXqv4vK3rHlg2yrWKqz5OkiRJklS2ohLjzHwvsGOZVRfVtY3FaxWfPFZ454WnN/I4SZIkSVLZikqM21D1WsUrPc4ZqyVJkiRpOnUuMYbq1yoe9DhnrJYkSZKk6VXM5FvTzBmrJUmSJGl6mRjXwBmrJUmSJGl6mRjXYHHG6qWcsVqSJEmSpoOJcQ0WZ6zeu/8YJ04ke/cfY9f8Yea2r590aJIkSZKkFXRy8q26VZ3pWpIkSZI0eSbGNak607UkSZIkabLsSi1JkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqNBNjSZIkSVKnmRhLkiRJkjrNxFiSJEmS1GkmxpIkSZKkTjMxliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6rSiEuOIWBMR74mI6/r3z4yI6yNiT//3xknHKEmSJEmaLUUlxsDLgFuW3L8MuCEzzwdu6N+XJEmSJKk2xSTGEbEVeAbwR0sWPwe4qn/7KuC5LYclSZIkSZpxxSTGwG8APwbcv2TZIzPzAED/9+YJxCVJkiRJmmFFJMYR8UzgYGYujPEcl0bE7ojYfdddd9UYnaTVsj1K5bA9SmWxTUplKiIxBp4MPDsi9gFXA18fEW8C7oyILQD93wcHPUFmXpmZOzJzx6ZNm9qIWdIAtkepHLZHqSy2SalMRSTGmfnKzNyamduAFwB/l5nfAVwLXNL/s0uAt00oREmSJEnSjCoiMR7icmBnROwBdvbvS5IkSZJUm1MnHcDJMnMemO/f/gRw0STjkSRJkiTNttIrxpIkSZIkNcrEWJIkSZLUacV1pe6Sm/fcy/zCUQ4eOsHmjWuY276eC84/bdJhSZIkSVKnmBhPyM177uX6m+7h4rkNbNuyln0HjrNr/jCAybEkSZIktcjEuGGDqsLzC0e5eG4D521dB8B5W9dx8dwGrrvxiImxJEmSJLXIxLhBw6rCBw+dYNuWtQ/6+21b1nLw0IlJhCpJkiRJneXkWw1aWhVesyYeqArPLxxl88Y17Dtw/EF/v+/AcTZvXDOhaCVJkiSpm0yMGzSsKjy3fT275g+zd/8xTpxI9u4/xq75w8xtXz+haCVJkiSpm+xK3aDFqvDiOGL4XFV4cRzxdTceeWD88c4LT3d8sSRJkiS1zMS4QYtV4ZPHGO+88HSgN/t0lUS46mWehj3OS0dJkiRJ6ioT4wY1URWuepmnYY8DvHSUJEmSpM4yMW5Y1arwIFUv8zTsccDQ57SaLEmSJGmWmRhPmaqXeVrpcYPWVa1QS5IkSdK0cFbqKVP1Mk/DHjds3bBLTkmSJEnSLDAxnjJVL/M07HHD1lWtUEuSJEnStLAr9ZRZaUKvQeOBR5kIbLl18wtHB15yatj2JEmSJGlamBhPoUETeq00HnjYRGCD1g275JTjjyVJkiTNAhPjGVJ1xuphhlWar7j67tq3J0mSJEltMzGeIU2NBx5UTXb8sSRJkqRZUDkxjogvBvZn5r01xqO+KmN3F2eXHjQeuG5tb0+SJEmSmjDSrNQR8YsRcUn/dkTE9cCHgAMRcWGTAXbR4tjdZz7lDF596Vk88ylncP1N93DznuHnIKrOWF1V29uTJEmSpCaMWjH+duBb+7e/GXg88NX95ZcDX1d7ZB1WdazwKDNP16nt7UmSJElSE0ZNjB8J7O/ffjpwTWb+v4i4G9jdSGQdNs7Y3WEzTzeh7e1JkiRJUt1GTYw/ATyaXnL8DcArlzw+Goir02Zl7K7XOJYkSZI0DUZNjP8MeEtEfAg4E3hHf/njgQ83EFenDbt28LTwGseSJEmSpsWoifGPALcB5wI/lpn39JdvAX6vicC6bBbG7jZxTWVJkiRJasKoifGjgNdm5v0nLf8N4JxaIxIw/WN3vcaxJEmSpGkx0uWagP8Azlpm+Zn9ddKDLI6TXmoax0lLkiRJmn2jJsYB5DLLzwCGX1xXneQ1jiVJkiRNi6FdqSPiN/s3E/iliDi6ZPUa4EnAe5sJTdNsFsZJS5IkSeqGlcYYf0X/dwBfChxbsu4Y8G7gNQ3EpRkw7eOkJUmSJHXD0MQ4M78OICLeALwsMz/dSlSSJEmSJLVkpFmpM/PFTQciSZIkSdIkjJQYR8RpwMuAi4DNnDRpV2Z+5biBRMQ5wBuBLwTuB67MzCsi4kzgT4FtwD7g+Zl5aNztSZIkSZIEo1/H+HeBi4H/DfwTy89QPa77gFdk5rsjYgOwEBHXA98F3JCZl0fEZcBlwI83sH0V4uY99zK/cPSBSbvmtq93rLIkFWzYcbvtY7qfIZKkKkZNjJ8LPC8z39lUIJl5ADjQv304Im4BzgaeA8z1/+wqYB4T45l18557uf6me7h4bgPbtqxl34Hj7Jo/DOAXG0kq0LDjNtDqMd3PEElSVaMmxkeB25sMZKmI2AY8AbgJeGQ/aSYzD0TE5rbiUHMGndGfXzjKxXMbOG/rOgDO27qOi+c2cN2NR/xSI0ktWG31d9hxG2j1mO5niCSpqlET418BfiQi/ntm3t9kQBFxBvBnwMsz89MRMerjLgUuBTj33HObC1BjG3ZG/+ChE2zbsvZBf79ty1oOHjoxiVBVke1RKsdq2mOV6u8dd9039Ljd5jHdzxBNAz8jpTKdsvKfALAT+FbgPyLi7RFx7dKfuoKJiLX0kuI3Z+af9xffGRFb+uu3AAeXe2xmXpmZOzJzx6ZNm+oKSQ1YekZ/zZp44Iz+/MJRNm9cw74Dxx/09/sOHGfzxjUTilZV2B6lcqymPQ47Pg9ad+x4Djxut31M9zNE08DPSKlMoybGHwd2AX8HfAz4xEk/Y4teafh1wC2Z+etLVl0LXNK/fQnwtjq2p8kZdkZ/bvt6ds0fZu/+Y5w4kezdf4xd84eZ275+QtFKUncMOz4PWhenMPC43fYx3c8QSVJVJV3H+MnAi4D3RcR7+8teBVwOXBMRLwE+AjyvhVjUoMUz+otjwOBzZ/QXx4Bdd+ORB8aw7bzw9InNbipJXTLs+Lx4++R1X7x1HXPb1w88bsPgY3rdVvoMkSRpkFHHGDcuM28EBg0ovqjNWNSsxTP6J49T23nh6UDvi81yX2KcbVSSmrXS8XnQukHHbRh8TG9K29uTJM2GkRLjiHgfQ65dnJlfWVtEmnlVz+ivNNtoE9fRtEI9e/yfSoONcnyusxpre5QklWLUivFbT7q/Fng8ve7Pv1NnQOqGKmf0h419a+I6mlaoZ4//U2llbVV/bY+SpJKMOsb455ZbHhH/A3h0rRFJAwwb+zbOdTS9pvLs8X8qlc/2KEkqybhjjP8c2A28tIZYpKGGjX275p2HK11H02sqzx7/p9J0sD1KkkoybmL8VOBoHYFIKxk29m1+4eiqZ1JdqdK80uysjo0r0zj/U0ntGac9evyVJNVt1Mm3rj15EbAFeAKwbDdrqQmDxrdVnUl1WKX5+U/bMPBxjo0r17Aq1LD/qaR2rXTcHsTjrySpCaNWjD9x0v37gfcDr8rMv603JGn1qs6kOqzSPOw5r7j6bsfGFWqc62RLak9TVyiQJKmKUSffenHTgUjjqjKTatVrKpc0Ns4uhQ9W9X8qqX11X6FAkqSqVjXGOCIeA3wZvWsa35KZtzYSldSSqhWLUsaq2qXwoawKS7OtlOOvJGm2jDrG+OHA64D/Rq8bdX9x/Bnwksw8PPDBUuGqVCyqjo2r26x3KaxaDW+7Ktx21d5eAuqycY6/th1J0iCjVoyvAL4S+Drgn/rLngz8PvAbwEtqj0wqWClVyVnuUjgt1fC245yW90VqStXjr21HkjTMqInxs4HnZuY/LFk2HxGXArswMVYHNVGVXG01Y5a7FK5UDS+l8tN21X7WewlIo6hy/B2n7ZRyvJEkNWfUxPhhPHRmaoC7AT8ZpBpUqWaU0qW7CcOq4SVVftqu2s9yLwGpSVXbTknHG0lSc0ZNjP8R+J8R8aLMPAoQEafTu4bxPw19pKSRVKlmNNWlu83qyKBtDauGN1U1rfK6x6naD9telfdFmiV1H4eqth17aWga2ctBWr1RE+MfBt4BfDQi/pXerNQXAPcA39hQbFKnVK1m1N2lu83qyLBtDauGX/POw7VXTau+7qpV+2HbAyq9L9KsaOI4VLXt2EtD08ZeDlI1o17H+N8i4nzgO4DHAQG8CXhzZn6mwfikziilEthmdWTYtl72gjOB5avh8wtHa3+vqr7uqlX7YdsDKr0v0qxo4jg0icvzWbXTJIwzR0fVfbaJ55TaNvJ1jPsJ8B82GIvUaaVUAtusjqy0rUHV8Cbeq3Fed5Wq/Urbq/K+SLOiqeNQm5fns2qnSak6RwcM7q00bJ9t4jmlSRj1Osa/ANyemb9/0vLvB87OzJ9qIjipS4ZVM9o829pm5brqtpoYW93U6646VtgKlbpsEj1oBrWdJnqF2CY1qrrnvqjaW2nYd5FxnrPN92UatqXJGrVi/CLgecssXwBeCZgYSzVYrprRdtWhzcr1ONuqu2raxOuuOoYasEKlTmu7B81KbaeJXiHSSpqY+2KlOTqqVJrH6QFVRSlzofi5OntGTYw3A3cts/wTwCPrC0fSydquOjQ10/WktzVuLFXOGFcdQ72o7nHLbY4pmwWz/vpK1vaxoYlrHM/K2OSSYumaJua+WGmOjiqV5qZ6QNX9vpS+LU3eqInxR4CvBW49aflTgf21RiTpQSZRdWhzDGtJ42UHxVL1jHHVMdQrrVvt9j68/xjH78vWxpSV8v8cx6y/vmnQ5rGhiWscz8LY5JJi6aIm5r6o2ltpWKX5+U/bUHsPqGFKmgtFs2XUxPgPgNdGxDrg7/rLLgJ+CfjlJgKTplmdszOWMlt1l1U9Y9z2/27Q9vL+wWO8oNr4r3FmPZ0GVgm6pYlrHFedQb6ktmU7mKwmPkNG6Y2x2kpz1eesahrmQgFn6p5Go16u6dci4izgN4HFPeMYcEVm/kpTwUnTqO7ZGUuZrbrLqp4xbvt/N2h769ZG7eO/qs56Oi0f/FYJuqWpaxzXPTa57bZlO5ispj5DqvRWWimWuntADTMNc6E4U/d0Ws3lml4ZET8PfBm96xh/IDOPNBaZNKXqnp1xnLF2npGsR0mzZ1fZXtUxZcNUnfV0Evtf3bO6NrE9TVbVttrEflJS27LH0mRN0zwcVVU5Xk7DXCizPlP3rBo5MQbIzHuAdzUUizQTmpidscrZ1lmo2pWipNmzq26v7vFf48x62qYmZnVtYnuavCpttYn9pKS2ZY+lyZuGeTiqGud4WfpcKLM8U/csW1ViLGllbc/OOEhJ49SmXUln7asYZ/xXlWu7rlShHqbuGbKbmNV1GNtdtzSxnwwbmzxO2xqm7ms4N8G281DT/p60fbxsc8zvLM/UPctMjKWaNXF92ipKGqc2C0o6a19FlfFfVa/t2na1dZzraw4z7eND1Y4mqkl1t61hmriGc91sOw81C+9Jm8fLtsf8tv1d0PkA6mFiLNWs7dkZBylpnJqmU1PV1qrV3SqPK2VmcNudlmpq3oIqFbEmruFcVd3Hhlk2C8eTpo6Xy+0PbY/5bfu74ErHlC62kSpMjKUGVKnO1a2kcWqaTk1UW6tWd6s+bqXra9bNdqdRNDFvQdWKWBPXcK67omdPjIeahQphE8fLQfvDHXfd1+qYXyhnpu6utpEqRk6MI+IrgO8DzgO+OzMPRMRzgdsy8z0NxSepoqbGgKo7mqi2Vq3uVn1c22MkbXcaRRP7ZdWKWBPXcK7yOpo4Nszyl/5ZmDG8iePloP3htW+5u4j5X5oy7L284uq7O9lGqhgpMY6IbwCuBd4OfD3wsP6q84DvAp7bQGySxtTmODXNnib2k6rV3WHVg5WqwsPOzDfRvayJdmc3uNlTd8Wo6iy4VXtVrLS91XbrbuLYsFIcw5Te5mblc7zu4+Wg/ShOqf/qDE2pOknYoPeytN4FJbetUSvG/xP4kcz83Yg4vGT5PPCK2qOS1KiSZhpVuZrYT6pWd4dVD6rG2Xb3smmJU9Op6iy4VffLYdur0q173dpqMQ47NjQxoV8pbW7WP8fr3i+/eOs65ravn/j8LytpYpKwknoXlN62IjNX/qOII8CXZ+a+fmJ8QWbeGhFfBNySmZN/JUvs2LEjd+/ePekwpM6IiIXM3LHcOtujlhr0objSl5Cqjxvmiqvv5plPOeNBXxb27j/2oEvmDIunrTPeVeIc1CZtj7NrWBuB5b9Qj9N+hm1vfuHowH0WWHbdH1/3KT7/jFNqPTYMi+NlLzhzYDse59iwHD8j29PEZ8VK26vzs2DYvgfLt52V9suV3pOqFeoqr32lttXWZ+ugNjlqxfgQcDaw76TlTwT2jxfayiLim4ArgDXAH2Xm5U1vU5JUv6pVgCaqI6VMOtRUnOqWtmfBHba9lbo3L7fu/oSdF55e67FhWBxNXeZNk9VmJb2Jz4KqQyKGGfaeVK1QD1s37LWXPpHeqInxW4BfjYjnAwmcGhH/BXgN8IamggOIiDXA7wA76SXh74qIazPzA01uV5LUjKpjK+sek1nKpEOLBp0pL6kbnMrW9hURBj3nON266zw2NDGhH5Q9RlLtXf2jic+Cqm1nJYPek6qT9g1bN6wKPc5EenVXr5czamL8k8AfA7cBAXyg//stwC+sequr8yTgw5l5K0BEXA08px+DJEmV1D25yzjVpGFnymdlkh11x0r7bFv7c9XLAQ2b7KuEqpbK0MRnQdttZ5wKdZXKb9U22cTY6+WMmhifyMxvj4ifotd9+hTgPZm5JyIeDnx6VVtdnbOB25fc3w9c2OD2JEkd0MSkQ1UNO1O+OHashIlhpFG03a27ShxVJ/Tz0jda1MRnQdttZ5wKdZXK77DPs2Ftsmplu6nE+HXAi/tV21sXF0bE5wN/S7OJaiyz7CEzhkXEpcClAOeee26D4Uhaie1R06JKl7u2L2NVNc5FtkdNQtvdulcbx0rtuMlL39gmZ0NTvXnabDvjVKirXkKtSpusMm9Blcr9qInxV0fEr2bm/1hcEBGPAK4HPrnqra7OfuCcJfe3Anec/EeZeSVwJfRm+Gs4JklD2B41y9q+jNW4bI/SQ02yx4htcjbMwiWzxq1Qr7byWzWWlZ6zrs/PURPjbwRujIhPZOblS5LiTwHPXPVWV+ddwPn9S0N9FHgB8MKGtylJ0kBtn7WXVL9SeoxoerXZ+6EpVSvUVXtjVImlrbHXIyXGmfmRiPhG4P9GxHF6yekngWdl5mdXvdVVyMz7IuKlwN/Qu1zT6zPz/U1uU5KkNs1C5UHqAtuqNFwTbaStsdejVozJzFsi4unADcA/A89uOilesu2/Bv66jW1JkjQJs1B5kLrAtioN1+bl4erc3sDEOCLexzKTXAH3AdvoXU8YgMz8yrEjkSRJkiRpAoZVjN/aWhSSJEmSJE1IZM7eZHgRcRdwW//uWcDHJxjOolLiAGMZxFgeatQ4Hp2Zm5ZbUWh7BGNZTilxgLEMMlabPKk9rub5mlZKHGAsg5QSSylxgJ+RbSklllLiAGMZZLzPyFlMjJeKiN2ZucM4PsdYlmcszcdRyusCYyk5DjCWQWa1TZYSBxjLIKXEUkocMLvtEYyl5DjAWAYZN5ZhY4w/DTwmMz8eEYdZfrwxAJn58KoBSJIkSZI0ScPGGP8gcLh/+6UtxCJJkiRJUusGJsaZeRVARJwK3AXclJmfaCuwGl056QD6SokDjGUQY3mouuMo5XWBsSynlDjAWAaZ1TZZShxgLIOUEkspccDstkcwluWUEgcYyyBjxTLSGOOIuBd4XGbuG2djkiRJkiSV5pQR/+5m4IubDESSJEmSpEkYtWL8zcDlwM8AC8A9S9dn5t2NRCdJkiRJUsNGTYzvX3J36QMCyMxcU3dgkiRJkiS1Ydis1Et9XaNRSJIkSZI0ISNVjCVJkiRJmlWjVoyJiPXA44HNnDRpV2b+eb1hSZIkSZLUjpES44h4GvAnwBcsszoBxxhLkiRJkqbSqJdrugL4K2BrZp5y0o9JsSRJkiRpao06K/U9wFdm5t7mQ5IkSZIkqT2jVoz/EfiSJgORJEmSJGkSBo4xjognLrn7+8BrIuJRwPuA40v/NjPf3Ux4kiRJkiQ1a1jFeDfwrv7vtwKPA64E/rm/bPeSv1HHRcS+iPhMRByJiDsj4g0RcUbD24yI+ImI+EhEfDoiro6Ihze5TalUE2qDWyLi2oi4IyIyIradtP7zIuL1/fb5sYj4kSbjkUpRaHt8fkT8U0QcjYj5JmORSlJoe3xNROyJiMMR8cGI+M4m49FohiXGXwQ8pv972M9jGo5R0+NZmXkG8ETgq4CfbHh73wm8CHgy8CjgYcBvNbxNqWRtt8H7gXcA/23A+p8FzgceDXwd8GMR8U0NxySVorT2eDfwG8DlDcchlai09ngP8Czg84FLgCsi4j83HJNWMDAxzszbRv1pM2CVLzM/Crwd+HKAiHh2RLw/Ij4ZEfMR8aWLf9s/i/fKiPhARBzqn8U7bcRNPQt4XWbenplHgF8GvrV/zW2ps9pqg5l5Z2b+LoN7Dn0n8D8z81Bm3gL8IfBd47w2adqU0h4z852ZeQ1wx/ivSppOBbXHn8nMD2bm/Zl5E/APwNeM/QI1lpEm3+p3v/mGJfd/OiL2R8TfRMSW5sLTNIqIc4CnA++JiMfSuwb2y4FNwF8DfxkR65Y85NuBbwTOAx7LkrN4/QPVUwZtqv+z9P7n0atQSZ3VYhscFsNGej05bl6y+GbgP632uaRpVkJ7lNRTYnuMiIfRq2K/f9zn0nhGnZX6ZxdvRG9SrlcBvwmsBX6t/rA0pf4iIj4J3Aj8PfCLwLcCf5WZ12fmceA19Lo8L+0u8tv9qu/dwC8A37a4IjMfkZk3Dtje24HviYhtEfH5wI/3l1sxVle13QaHWRy/9aklyz4FbKjwXNI0Kqk9Sl1Xcnv8fXonjv+mhufSGAbOSn2SRwP/3r99MfAXmfkrEfG3+E/U5zw3M9+5dEH0ZjJ/oLt9Zt4fEbcDZy/5s9uX3L6NXpVpFK8HzgHm6e3Lv0ave/X+VUcuzYa22+AwR/q/Hw7cu+T24RqeW5oGJbVHqeuKbI8R8av0unV/XWZmnc+t1Ru1YnwvnzvLfxGwuGN59l8ruYPeiRWgN5M0vWT2o0v+5pwlt89lxPFP/XEZP5OZ2zJzK70uKB896bmlrmusDQ6TmYeAA8AFSxZfgF3F1G0TaY+SljXR9hgRPwd8M/ANmfnpup5X1Y2aGP8D8GsR8VPADnp98KHX1/72gY+S4BrgGRFxUUSsBV4BfBb4pyV/8wMRsTUizqTXTf9PR3niiDgzIs6Lni8Dfh14dWbeX/NrkKZZY20QoD8Ryef1737eSROTvBH4yYjYGBGPA74X+OPqL0WaehNrjxGxpn//VOCUiDitH4PUVZNsj68EXgjszMxPjPk6VJNRE+OXAseAbwG+PzMXz5Z8M3al1hCZ+e/Ad9C7jNLH6XV1flZmHlvyZ28B/ha4tf/z84sronfNua8d8PRn0TtJcw+98cavz8wra38R0hRruA0CfIbPdZv+YP/+op8B9tLrfvb3wK9m5jvGekHSFJtwe3xR//7vAV/bv/2H47weaZpNuD3+Ir0K9J7+8xyJiFeN+ZI0prA7uyYpIvYB33PyuA9J7bANSuWwPUrlsD12z6gVY0mSJEmSZpKJsSRJkiSp0+xKLUmSJEnqNCvGkiRJkqROO3XSATThrLPOym3btk06DKkzFhYWPp6Zm5ZbZ3uU2jeoTdoepfb5GSmVZVCbHJgYR8RPj/rkmfnqqoE1Ydu2bezevXvSYUidERG3DVpne5TaN6hN2h6l9vkZKZVlUJscVjF+3kn3Hw2sBxavYfwo4CiwDygqMZYkSZIkaVQDE+PM/IrF2xHxYuA7gUsy8yP9ZecCbwDe3HSQkiRJkiQ1ZdQxxj8NPHcxKQbIzI9ExCuAtwGvbyI4SZIkSQ918557mV84ysFDJ9i8cQ1z29dzwfmnTTosaWqNmhg/EnjYMstPA86qLxxJkiRJw9y8516uv+keLp7bwLYta9l34Di75g8DmBxLFY16uabrgT+MiK+OiDX9n68G/qC/TpIkSVIL5heOcvHcBs7buo41a4Lztq7j4rkNzC8cnXRo0tQaNTH+HuB24J+Ae/s//wh8FPjeZkKTJEmSdLKDh06wbcvaBy3btmUtBw+dmFBE0vQbqSt1Zt4FPD0iHgs8Dgjglsz8UJPBSZIkSXqwzRvXsO/Acc7buu6BZfsOHGfzxjUTjEqabqNWjAHoJ8I3AX9pUixJkiS1b277enbNH2bv/mOcOJHs3X+MXfOHmdu+ftKhSVNrpIpxRKwFfgH47/Qm4XoscGtE/DJwW2b+bnMhSpIkSVq0OMHWdTceeWBW6p0Xnu7EW9IYRp2V+meAZwHfAbxlyfL/B/w4YGIsSZIkteSC808zEZZqNGpi/G3Ad2fm30fE/UuW/xu96rEkSZIkSVNp1DHGjwJuW2b5qYyeXEuSJEmSVJxRE+P3A09dZvnzgYX6wpEkSZIkqV2jVnt/DnhTRJwDrAGeFxGPA14IPKOp4CRJkiRJatpIFePM/Et61eFvAO6nNxnX+cCzMvOdzYUnSZIkSVKzRr1c03mZ+TfA3yyz7qLMvKH2yCRJkqQOu3nPvcwvHH3gkkxz29c7E7XUkFG7Uv9tRDw5Mz+2dGFEPA3YBWyoPTJJkiSpo27ecy/X33QPF89tYNuWtew7cJxd84cBTI6lBow6+dY7gOsj4hGLCyJiJ/A24EcbiEuSJEnqrPmFo1w8t4Hztq5jzZrgvK3ruHhuA/MLRycdmjSTRk2MX0rvmsV/FREP6yfFfwH8cGb+QVPBSZIkSV108NAJtm1Z+6Bl27as5eChExOKSJpto06+lcCLgE8D/4de9+mXZ+aVdQYTEY+IiLdGxAcj4paI+JqIODMiro+IPf3fG+vcpiRJklSazRvXsO/A8Qct23fgOJs3rplQRNJsG5gYR8QTl/4AXwn8IvAo4I3AwpJ1dbkCeEdmPg64ALgFuAy4ITPPB27o35ckSZJm1tz29eyaP8ze/cc4cSLZu/8Yu+YPM7d9/aRDk2bSsMm3dgMJxJJli/e/H/i+/u2kd23jsUTEw4GnAt8FkJnHgGMR8Rxgrv9nVwHzwI+Puz1JkiSpVIsTbF1345EHZqXeeeHpTrwlNWRYYvxFrUXR8xjgLuANEXEBsAC8DHhkZh4AyMwDEbF5uQdHxKXApQDnnntuOxFLWpbtUSqH7VEqy2ra5AXnn2YiLLVkYFfqzLxt1J+aYjkVeCLwe5n5BOAeVtFtOjOvzMwdmblj06ZNNYUkqQrbo1QO26NUFtukVKaRJt+KiF+IiO9fZvn3R8T/rCmW/cD+zLypf/+t9BLlOyNiS397W4CDNW1PkiRJkqSRL9f0IuA9yyxfAL6zjkAy82PA7RHxJf1FFwEfAK4FLukvu4TetZMlSZIkSarFsDHGS22mN/73ZJ8AHllfOPwg8OaIWAfcCryYXvJ+TUS8BPgI8LwatydJkiRJ6rhRE+OPAF9LL1ld6qn0ukDXIjPfC+xYZtVFdW1DkiRJkqSlRk2M/wB4bb+S+3f9ZRcBvwT8chOBSZIkSZLUhpES48z8tYg4C/hNYF1/8THgisz8laaCkyRJkiSpaaNWjMnMV0bEzwNfBgTwgcw80lhkkiRJkiS1YOTEGCAz7wHe1VAskiRJkiS1buTEOCK+Dvg24Fw+150agMz8+prjkiRJkiSpFSNdxzgivgt4O7ABmKN36aaNwBPpXWtYkiRJkqSpNFJiDPwo8NLM/DbgOPDKzHwC8CbAccaSJEmSpKk1amL8GOCd/dufBc7o3/5t4LtqjkmSJEmSpNaMmhh/gl43aoCPAl/ev/0FwMPqDkqSJEmSpLaMOvnWPwDfALwPuAb4zYjYCVwEXN9QbJIkSZIkNW7UxPilwGn9278E3Ac8mV6S/PMNxCVJkiRJUitGSowz8+4lt+8HfrmxiCRJkiRJatFqrmN8GvBC4Mv6iz4A/ElmfqaJwCRJkiRJasOo1zF+InAr8GvAk/o/rwFu7a+TJEmSJGkqjTor9ZXAjcDWzHxqZj4VOAf4v/11kiRJkiRNpVG7Uv8n4Dsz857FBZl5T0S8GtjdSGSSJEmSJLVg1IrxB4FHLbN8C/Ch+sKRJEmSJKldAyvGEXHmkrs/Se/axa8G/qW/7Kv7yy9rLjxJkiRJkpo1rCv1x4Fccj+AtyxZFv3fbwPW1B+aJEmSpDrdvOde5heOcvDQCTZvXMPc9vVccP5pkw5LmrhhifHXtRaFJEmSpEbdvOderr/pHi6e28C2LWvZd+A4u+YPA5gcq/MGJsaZ+fdtBiJJkiSpOfMLR7l4bgPnbV0HwHlb13Hx3Aauu/GIibE6b9TJtx4QEe+LiHOaCEaSJElSMw4eOsG2LWsftGzblrUcPHRiQhFJ5Vh1YgxsA9au9EeSJEmSyrF54xr2HTj+oGX7Dhxn80anC5JGvY6xJEmSpCkwaIKtue3r2TV/+CFjjHdeePqkQ5Ymrkpi/A/AZ+oORJIkSdJ4Rplg67objzyQNO+88HTHF0tUSIwz8+lNBCJJkiRpPCtNsLX4I+nBBibGEfGdoz5JZr6xnnAkSZIkVeUEW1I1wyrGv3PS/XX0Jt26v3//FOA48FnAxFiSJEmasMUJthYrxuAEW9IoBs5KnZkbFn+AFwD/CnwtcFr/52uB9wIvbCFOSZIkSStYnGBr7/5jnDiR7N1/jF3zh5nbvn7SoUlFG3WM8WuA787Mf16y7B8j4uXAHwPX1RyXJEmSpFVygi2pmlET423APcssPwqcW1s0kiRJksbiBFvS6g3sSn2Sm4DfjIizFxf0b78W+JcmApMkSZIkqQ2jJsYvAb4A2BcR+yJiH7AP2Ax8b13BRMSaiHhPRFzXv39mRFwfEXv6vzfWtS1JkiRJkmDExDgz9wJfCTwD+HV6leKnA1+RmR+uMZ6XAbcsuX8ZcENmng/c0L8vSZIkSVJtRh1jTGYm8Lf9n9pFxFZ6ifcvAD/SX/wcYK5/+ypgHvjxJrYvSZIkSeqmkRPjiPj/gB8Avgj48sy8NSIuA27NzGtqiOU3gB8DNixZ9sjMPACQmQciYnMN25EkSZJ0kpv33Mv8wtEHZrOe277eSbzUGSN1pe5flukngSuBWLLqo8BLxw0iIp4JHMzMhTGe49KI2B0Ru++6665xQ5I0BtujVA7bo1SWUtvkzXvu5fqb7uGZTzmDV196Fs98yhlcf9M93Lzn3kmHJrVi1Mm3vh/43sy8ArhvyfJ3A/+phjieDDy7P6nX1cDXR8SbgDsjYgtA//fBQU+QmVdm5o7M3LFp06YaQpJUle1RKoftUSpLqW1yfuEoF89t4Lyt61izJjhv6zountvA/MLRSYcmtWLUxPjRwL8ts/w48LBxg8jMV2bm1szcBrwA+LvM/A7gWuCS/p9dArxt3G1JkiRJerCDh06wbcvaBy3btmUtBw+dmFBEUrtGTYxvBZ64zPKnAx+oL5yHuBzYGRF7gJ39+5IkSZJqtHnjGvYdOP6gZfsOHGfzxjUTikhq16iTb70G+O2IWE9vjPHXRMSL6E2W9d11BpSZ8/RmnyYzPwFcVOfzS5IkSXqwue3r2TV/mIvnNrBty1r2HTjOrvnD7Lzw9EmHJrVipMQ4M98QEacCvwisB/4XvYm3figz/7TB+CRJkiQ1bHH26etuPPLArNQ7LzzdWanVGau5jvEfAn8YEWcBp2TmQYCIOCczb28qQEmSJEnNu+D80wYmwlUv5eQloDQtRk6MF2XmxwEi4guBn6LXlXrsCbgkSZIklWfxUk4nd7MGhia5VR8nTcLQybci4hER8eaIuCsi7oiIH4qen6E3IdeTqHmMsSRJkqRyVL2Uk5eA0jRZqWL8i8BTgauAbwJeS2926NOBb87Mv282PEmSJEmTVPVSTl4CStNkpcs1PQN4cWb+KPBsejNS783MrzcpliRJkmZf1Us5eQkoTZOVEuNH0b9OcWbeCtwL/GHTQUmSJEkqw+KlnPbuP8aJE8ne/cfYNX+Yue3rG3mcNAkrdaU+BVh6mucE4KAASZIkqSOqXsrJS0BpmqyUGAfwpoj4bP/+afQu2fSg5Dgzn91EcJIkSZImb9ilnJp4nNS2lRLjq066/6amApEkSZIkaRKGJsaZ+eK2ApEkSZIkaRJWmnxLkiRJkqSZZmIsSZIkSeq0lcYYS5IkSdJAN++5l/mFow/MPD23fb0TbmnqmBhLkiRJquTmPfdy/U33cPHcBrZtWcu+A8fZNX8YwORYU8Wu1JIkSZIqmV84ysVzGzhv6zrWrAnO27qOi+c2ML9wdOUHSwUxMZYkSZJUycFDJ9i2Ze2Dlm3bspaDh05MKCKpGhNjSZIkSZVs3riGfQeOP2jZvgPH2bxxzYQikqoxMZYkSZJUydz29eyaP8ze/cc4cSLZu/8Yu+YPM7d9/aRDk1bFybckSZIkVbI4wdZ1Nx55YFbqnReePtLEW85mrZKYGEuSJEmq7ILzT1t1Quts1iqNXaklSZIktcrZrFUaE2NJkiRJrXI2a5XGxFiSJElSq5zNWqUxMZYkSZLUKmezVmmcfEuSJElSq8aZzVpqgomxJEmSpNZVmc0avMyTmmFiLEmSJGkqeJknNcUxxpIkSZKmgpd5UlNMjCVJkiRNBS/zpKbYlVqSJElSUQaNI168zNN5W9c98Lde5kl1MDGWJEmSVIxh44gXL/N08rqdF54+4ag17UyMJUmSJBVj6Thi4IFxxNfdeISXveBMoPzLPDlz9vQpJjGOiHOANwJfCNwPXJmZV0TEmcCfAtuAfcDzM/PQpOKUJEmS1JyVxhFXvcxTW5w5ezqVNPnWfcArMvNLga8GfiAivgy4DLghM88HbujflyRJkjSDFscRLzVN44idOXs6FZMYZ+aBzHx3//Zh4BbgbOA5wFX9P7sKeO5EApQkSZLUuMVxxHv3H+PEiWTv/mPsmj/M3Pb1kw5tJM6cPZ2K6Uq9VERsA54A3AQ8MjMPQC95jojNk4xNkiRJUnMWuxuXMI64ylhhZ86eTsUlxhFxBvBnwMsz89MRMerjLgUuBTj33HObC1DSimyPUjlsj1JZbJOjKWEccdWxws6cPZ2K6UoNEBFr6SXFb87MP+8vvjMitvTXbwEOLvfYzLwyM3dk5o5Nmza1E7CkZdkepXLYHqWy2CanR9Wxwhecfxo7Lzyd6248wk9f+XGuu/FIkTNn68GKqRhHrzT8OuCWzPz1JauuBS4BLu//ftsEwpMkSZLUIeOMFS6h4q3VKSYxBp4MvAh4X0S8t7/sVfQS4msi4iXAR4DnTSY8SZIkSV2x0lhhr1U8W4pJjDPzRmDQgOKL2oxFkiRJUrcNGyvstYpnTzGJsSRJkiSVYtjs2FdcffcD44+BB8YfX3fjERPjKWViLEmSJEnLGDRW2GsVz56iZqWWJEmSpNItjj9eymsVTzcrxpIkSZK0CuNcq9hJu8pkYixJkiRJqzBs/PEwTtpVLhNjSZI0FqsfkkpR5XhU9RhW5VrF8wtHnbSrUCbGkiSpMqsfkkpR5XjU9jHMSbvKZWIsSZIqs/ohqRTDjkeL60+uCrd9DFuctGtxe+CkXaUwMZYkSZVZ/ZBUikHHow/vP8bx+3LZqnDbx7BxJu1Ss0yMV8ExVJIkPZjVD0mlGHQ8yvsZWBVu+xhWddIuNc/EeESOoZIk6aGsfkgqxaDj0bq1MbAq/PynbWj9GFZl0i41z8R4RI6hkiTpoax+SCrFoOPR/MLRgVXh0o5hw3qo2nu1WSbGI3IMlSRJy7P6IakUg45Hw6rCpRzDhvVQBey92rBOJsZVzrZUHX/QxJkdzxZJkppitULSrCmtKjzISrNq23u1WZ1LjKuOFa4yhqqJccmOdZYkNcVqhaRZVUpVeJiVeqjae7VZnUuMq44VrnKmqYlxyY51liQ1ZdarFVa8JZVspR6qVWfP9tg3ms4lxuOMFV7tmaYmxiU71lmS1JRZrlbY40pS6VbqoVpl9myPfaPrXGLc5rXKmtiW14uUJDWliWpFKWOW7XElqXSj9FBd7Tjp0o59Jc+/1LnEuM3rLTaxLa8XKUlqSt3VipLGLNvjStI0GNZDtco46ZKOfaXPv9S5xLjNWema2Na0zKo3DsdBSNJk1F2tKGnMsj2uJHVRSce+larXVXKAOivinUuMod1Z6ZrY1jTMqleV4yAkabLqrFaUNGbZHleSuqikY9+wz4SqOUCdFfFOJsYqV2njICRJ1TU1w+owgyoOK1XDSxkLLUlVVD32tWnYZ0LVHKDOiriJsYpS0jgISdJ4mphhdZiVKg6DKt4ljYWWpNWqeuxr27DPhGveebhSDlBnRdzEWEUpaRyEJGk8TcywOkzVikNJY6ElabWmpcflsM+E+YWjlXKAOiviJsYCqncTq7t72ThnfezqJknlqXuGVRh8vK/a66iksdCStFpN9bhs4rv1oOP+SjnAsFjqqoibGKvyYPcmJsqqetbHSbskqRuGHe+r9jqaxFhoSapLEz0u2/5uPSwHaCsWE+Oa1H1Gpc0KbtWp05vqtjHsrE8TsbRdLS+lsl1KHJK0GsOO91V7HbU9FnoSPOZLs6uJHpdNXFppJYNygLa6ipsY16DusxhtV3CrTp3e9kRZTcTS9ntdSmW7lDgkabWGHe+r9jpqeyx02zzmS7OtiR6XTVxaqaq2cg4T45NM+sLS4zxfE9OcD3vOlbpt1H0maZxYqjxnE9XySZx9qxJHKaxwaCn3h24Z9P9e6XhfdaxZ1bHQ03CZp0n0qpLUrirHvqrfrdv+PtvW5LwmxkuUcGHpcZ6v6uOqTp3+/KdtGPi4Js4kDXt9w2Kp+pxNVKhLOfs2DZfFssKhpdwfumXY/7vOS3M0GSeUc5mntntVSZoOVb9bD8sPmjhutHXcNzFeooQLS4/zfFUruFWnTh/2uCuuvrv2iuSw17dSF5IqlYcmquVVtwfUWpGYhstijXM20grH7JmWXg6qx7D/98tecCZQRtfm0i7zVLXKXuX1TWIelTZNQ4zSuKp+tx6WHzRx3KjzkkzDmBgv0UTFtYomJg+peuHvlWIZ9LgmKpJVY6laeWiiWl5lex/ef4zj92WtFYmSKi6DVK2uQzlVGtVnGno5qD4r/b/rujTHuEq6zFMTVfYm2t00VKGnIUapDlW/W1f9/jyONo77JsZLVD2jOs6A99VWcKvGUbWCWzWWJsYfV42lauWhiWp5le3l/YOrDlCtIlG1wt6EJqrrTVSam5iFvJSKREkzs1eteDXxXjaxr0zD/lCCaejVAuNd5qnufWGcKnub7a7t3h9V3udS5gORmtZEzjHs+3PppiIxjohvAq4A1gB/lJmXN7Gdcapoqz2LUbWCWzWOcc76VollnOp13bFUrTw0VS1f7fbWrY1GKhJVKux1f/A3Uc1fvL3cuqqV5mHrhr0n01DZLmlmdhj8njR1TKk7zibWdfEL9zT0aoHql3lq4thQ9bOn6nG4lHlZhqn6PpcyH4jUhrpzjmk5fi+n+MQ4ItYAvwPsBPYD74qIazPzA3Vvq63+69D+GdO2z743Ub2uqu2eAHVvb6Uzb9N8Nr+Jav7i7TorzcPWDaselFTZrrsy0sTjhr0nw/aHlY4pdVeMqu4rTexjs6rNz+NxjBLnavdZqLYvNDGOuIl2NyzOuo97Vdtx1R5Lpe2f0iRMy/F7OcUnxsCTgA9n5q0AEXE18Byg9sQY2hu31PZ4uUmcvWmiel1Fmz0Bmtpe3RWJUs7mN1XNr7vSPGzdOLOXt1XZ3rv/GLd+9HitlZEmHrfS/6DKMaWJitFKcba5j03DF42q2vo8HtewOKt+DlbZF5oaR1x3uxsU52POXlt7b4uq73PV2Xgl9UzL8ftk05AYnw3cvuT+fuDCCcVSm5IquG2b9dde9/aaqEgMi6XN/0+T1fw6K83D1o0ze3lble3XvuVufviFZ9ZaGWnicSv9DwZpYjz6OHG2uY9N4xcPVR+b3MRs3VWPw1Xb3aA4m+ilUfV9ntXxk5KGm4bEOJZZlg/5o4hLgUsBzj333KZjGltJFdy2deG11729JioSg9Tx/xm1PTZVza+70jxsXdXZy6s+5+Lt1aw7cu/9tVdGmnjcSv+DQZoYjz5OnG3uY6OYts/HLmiq10ubvZzGaXfLxdlEj55x2nGT4ydtk1KZpiEx3g+cs+T+VuCOk/8oM68ErgTYsWPHQxLn0pRUwW1bl197m9oeW73UqO2xxGp+nVXoUirbZ5x2Su2VkSYeN8r/YLX/16rvZRP7StV141aopu3zsQuaOjY0FctqH1clziZ6abR9vBmVbVIqU2SW3R4j4lTgQ8BFwEeBdwEvzMz3D3rMjh07cvfu3S1FKJVp0BivJhLPiFjIzB3LrZv19tjE+zzsOWH58XTD1j3m7LXLjjFeKcaqr63NfW8lVd/Lkk7UVXk/B7XJWW+Ps66ktjVMlTjrPu6V9J50+TNSKtGgNll8xTgz74uIlwJ/Q+9yTa8flhRL6rEy344m3ucmqpU377m3tcpISfteExWjtpX0fmqypmVfqBJnk700JGkUxVeMq/Dsm9Quz4ZLZbFiLJXDz0ipLIPa5CmTCEaSJEmSpFKYGEuSJEmSOm0mu1JHxF3Abf27ZwEfn2A4i0qJA4xlEGN5qFHjeHRmblpuRaHtEYxlOaXEAcYyyFht8qT2uJrna1opcYCxDFJKLKXEAX5GtqWUWEqJA4xlkPE+I2cxMV4qInYPGtfRxTjAWAYxlubjKOV1gbGUHAcYyyCz2iZLiQOMZZBSYiklDpjd9gjGUnIcYCyDjBuLXaklSZIkSZ1mYixJkiRJ6rQuJMZXTjqAvlLiAGMZxFgequ44SnldYCzLKSUOMJZBZrVNlhIHGMsgpcRSShwwu+0RjGU5pcQBxjLIWLHM/BhjSZIkSZKG6ULFWJIkSZKkgUyMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqNBNjVRIR+yLiMxFxJCLujIg3RMQZDW9zS0RcGxF3RERGxLaT1v9KRNweEZ+OiNsi4ieajEcqSYltcsnfnRkRd0XEjU3GI5WixPYYEX8cEcf6MS3+rGkyJqkEJbbH/t88LSLeHRH39L+/Pr/JmLQyE2ON41mZeQbwROCrgJ9seHv3A+8A/tuA9a8DHpeZDwf+M/DCiPivDccklaS0Nrnol4FbGo5FKk2J7fFXMvOMJT8nGo5JKkVR7TEivgx4C/ATwOcDjwcWGo5JKzAx1tgy86PA24EvB4iIZ0fE+yPikxExHxFfuvi3/bN2r4yID0TEof5Zu9NG3M6dmfm7wLsGrP/3zLxnyaL7gS+u/MKkKVVKm+w//9f043jDWC9KmlIltUep6wpqjz8J/EFmvj0z78vMT2Tm3jFfnsZkYqyxRcQ5wNOB90TEY4E/AV4ObAL+GvjLiFi35CHfDnwjcB7wWJactesfmJ4yRiyXRcQRYD9wOr2zcVKnlNIm+900fwd4KZBVnkOadqW0x77/LyLujoiFiFipp4c0cwpqj1/df473RcSBiHhTRJxZ8blUExNjjeMvIuKTwI3A3wO/CHwr8FeZeX1mHgdeAzyMXtfmRb+dmbdn5t3ALwDftrgiMx+RmZXHIWbm5cAGel1l/hfwqarPJU2h0trkDwE3Zabdw9RFpbXH3wTOBzYDPwX8cUQ8ueJzSdOmtPa4FXgRva7W5/e3+1sVn0s1OXXSAWiqPTcz37l0QUQ8Crht8X5m3h8RtwNnL/mz25fcvg14VJ1BZWbSOxP4jcDPAT9S5/NLBSumTfa3+0PA9nGfS5pSxbTH/rbeveTuX0fEm4H/CvxjHc8vFa6o9gh8BnhDZn6oH8svAu8c/hA1zcRYdbsD+IrFOxERwDnAR5f8zTlLbp/bf0wTTqXX9UXqskm1yScBW4AP9DbJw4CHRcTHgLOd9EcdVdJnZALR0HNL02CS7fFfcYhRcexKrbpdAzwjIi6KiLXAK4DPAv+05G9+ICK29sdSvAr401GfvD/pwef1737e4iQIEXFKRHxfRGyMnicBPwDcUMNrkqbZRNokvclNttGbafPxwE8D7wEeb1KsDptUeyQiviUizuh/Xn4D8B3AtWO+HmmaTaw90puQ8sUR8ZiIWA/8OHDdGK9FNTAxVq0y89/pfdj+FvBx4Fn0psg/tuTP3gL8LXBr/+fnF1dE7xpzXztkE58BjvRvf7B/f9HFwF7gMPCmfgyO11CnTapNZuZnM/Njiz/0xvsf79+WOmnCn5Evo1cJ+yTwq8D3Zub8GC9HmmqTbI+Z+XrgjcBN9Lpof5be8CNNUPSGY0rtiIh9wPecPM5D0mTYJqVy2B6lctgeu8eKsSRJkiSp00yMJUmSJEmdZldqSZIkSVKnWTGWJEmSJHXaTF7H+Kyzzspt27ZNOgypMxYWFj6emZuWW2d7lNo3qE3aHqX2+RkplWVQm5zJxHjbtm3s3r170mFInRERtw1aZ3uU2jeoTdoepfb5GSmVZVCbtCu1JEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqtKIS44h4RES8NSI+GBG3RMTXRMSZEXF9ROzp/9446TglSZIkSbOjqMQYuAJ4R2Y+DrgAuAW4DLghM88HbujflyRJkiSpFsUkxhHxcOCpwOsAMvNYZn4SeA5wVf/PrgKeO4n4JEmSJEmzqZjEGHgMcBfwhoh4T0T8UUScDjwyMw8A9H9vXu7BEXFpROyOiN133XVXe1FLegjbo1QO26NUFtukVKaSEuNTgScCv5eZTwDuYRXdpjPzyszckZk7Nm3a1FSMkkZge5TKYXuUymKblMpUUmK8H9ifmTf177+VXqJ8Z0RsAej/Pjih+CRJkiRJM6iYxDgzPwbcHhFf0l90EfAB4Frgkv6yS4C3TSA8SZIkSdKMOnXSAZzkB4E3R8Q64FbgxfSS92si4iXAR4DnTTA+SZIkSdKMKSoxzsz3AjuWWXVRy6FIkiRJkjqimK7UkiRJkiRNgomxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqNBNjSZIkSVKnmRhLkiRJkjrNxFiSJEmS1GkmxpIkSZKkTjMxliRJkiR12qmTDmCa3LznXuYXjnLw0Ak2b1zD3Pb1XHD+aZMOS5IkSZI0BhPjEd28516uv+keLp7bwLYta9l34Di75g8DmBxLkiRJ0hSzK/WI5heOcvHcBs7buo41a4Lztq7j4rkNzC8cnXRokiRJkqQxmBiP6OChE2zbsvZBy7ZtWcvBQycmFJEkSZIkqQ4mxiPavHEN+w4cf9CyfQeOs3njmglFJEmSJEmqg4nxiOa2r2fX/GH27j/GiRPJ3v3H2DV/mLnt6ycdmiRJkiRpDE6+NaLFCbauu/HIA7NS77zwdCfekiRJkqQpZ2K8Checf5qJsCRJkiTNmKK6UkfEmoh4T0Rc179/ZkRcHxF7+r83TjpGSZIkSdJsKSoxBl4G3LLk/mXADZl5PnBD/74kSZIkSbUpJjGOiK3AM4A/WrL4OcBV/dtXAc9tOSxJkiRJ0owrJjEGfgP4MeD+JcsemZkHAPq/Nw96cERcGhG7I2L3XXfd1WigkoazPUrlsD1KZbFNSmUqIjGOiGcCBzNzoepzZOaVmbkjM3ds2rSpxugkrZbtUSqH7VEqi21SKlMps1I/GXh2RDwdOA14eES8CbgzIrZk5oGI2AIcnGiUkiRJkqSZU0TFODNfmZlbM3Mb8ALg7zLzO4BrgUv6f3YJ8LYJhShJkiRJmlFFJMZDXA7sjIg9wM7+fUmSJEmSalNKV+oHZOY8MN+//QngoknGI0mSJEmabaVXjCVJkiRJapSJsSRJkiSp04rrSt0lN++5l/mFoxw8dILNG9cwt309F5x/2qTDkiRJkqROMTGekJv33Mv1N93DxXMb2LZlLfsOHGfX/GEAk2NJkiRJapFdqSdkfuEoF89t4Lyt61izJjhv6zountvA/MLRSYcmSZIkSZ1iYjwhBw+dYNuWtQ9atm3LWg4eOjGhiCRJkiSpmzrZlbqEsb2bN65h34HjnLd13QPL9h04zuaNa1qNQ5IkSZK6rnOJcSlje+e2r2fX/OGHxLHzwtNbi0GSJEmS1MHEeOnYXuCBsb3X3XiEC84/rbVq8uJzXnfjkQe2tfPC0514S5IkSZJa1rnEeNjY3raryRecf5qJsCRJkiRNWOcm31oc27vU4theZ4qWJEmSpO7pXGK8OLZ37/5jnDiR7N1/jF3zh5nbvt6ZoiVJkiSpgzrXlXrY2N75haPOFC1JkqTilXCVFWmWdC4xhsFje50pWpIkSaUr5Sor0izpZGI8iDNFS5IkqXQrXWVF0uqZGJ/EmaIlSZJUMufFkepnYixJkiQVaNA44sWrrLQ5L45jmjXrTIwlSZKkwgwbR9z2vDiOaVYXmBhLkiRJhRk2jvhlLzgTaG9eHMc0qwtMjGtSd/cSu6tIkiR110rjiNucF8cxzeoCE+Ma1N29xO4qkiRJ3TaJccTTEIvUlFMmHcAsWNq9ZM2aeKB7yfzC0SKeT5IkSdNlcRzx3v3HOHEi2bv/GLvmDzO3fX2nY5GaYsW4BnV3L7G7iiRJUrct9hJsaxzxtMQiNcXEuAZ1dy+xu4okSZLaHEe8kiZicU4dlcSu1DWou3uJ3VUkSZI0yxbn1HnmU87g1ZeexTOfcgbX33QPN++5d9KhqaOsGNeg7u4ldleRJElS29qs4HoJKJXGxLgmdXcvKanrjCRJkmZb21dFcU4dlcau1JIkSVLHtX1VlMU5dZZyTh1NkomxJEmS1HFtV3CdU0elsSt1w5xtT5IkSaUY9N10nKuiVPm+65w6Ko2JcYPaHqshSZIkDTLsu+liBffkdTsvPL3yc46SHPudWKUwMW5QU7PtWYWWJEnSag37bvqyF5wJrL6C6/ddzQoT4wY1MVbDKrQkSZKqWOm7aZUKrt93NSuKSYwj4hzgjcAXAvcDV2bmFRFxJvCnwDZgH/D8zDw0qThXY5yxGoNM4ppvnrGTJEmaHk2MIx6ktO+7fm9VVcUkxsB9wCsy890RsQFYiIjrge8CbsjMyyPiMuAy4McnGOfIqo7VGKbtGQM9YydJkjQ9mhhHPExJ33f93qpxFJMYZ+YB4ED/9uGIuAU4G3gOMNf/s6uAeaYkMW5itr0mzsoNM4kKtSRJkqppYhzxMCV93/V7q8ZRTGK8VERsA54A3AQ8sp80k5kHImLzgMdcClwKcO6557YU6crqnm2vibNyw7RdodZsKLU9Sl1ke5TK0nSbbGIc8UpK+b7r91aNo7jEOCLOAP4MeHlmfjoiRnpcZl4JXAmwY8eObC7CyWr7mm8rnbEbNo7DMR7d1ZX2KE0D26NUlqbbZNu9C5tQ9fvutLx2vyOXqajEOCLW0kuK35yZf95ffGdEbOlXi7cABycXYRnavObbsDN2w8ZxAI7xkCRJalnbvQubUuX77jS8dsdBl6uYxDh6peHXAbdk5q8vWXUtcAlwef/32yYQXmcNO2N3xdV3DxzHATjGQ5IkqWVt9y4syTS8dsdBl6uYxBh4MvAi4H0R8d7+slfRS4iviYiXAB8BnjeZ8KZDE10zBp2xW2kcx7B1diGRJElqRpu9Cydh2PfI0l+746DLVUxinJk3AoMGFF/UZizTqu2uGSuN4xi0zi4kkiRJqmLav0dOyzjoLiomMdb4VuqaUXeVdqVxHIPW2YVEo7JngSRJWqqk75FVvqdMwzjorjIxniHDumY0cXZtlHEcy6275p2H7UKiFU37GWFJklS/UroiV/2eMg3joLvKxHiGDOua0dTZtWHjOAatG6cLiRXE7ijpjLAkSSpDKV2Rx/meUvo46K46ZdIBqD6LXTP27j/GiRPJ3v3H2DV/mLnt64s5u7ZSnMMsnpl75lPO4NWXnsUzn3IG1990DzfvubelyNWmkvZZSZJUhqrfI+vm95TZY8V4hgzrmjG/cHTo2bU2K7FVu5BYQeyWUs4IS5KkcpTSFXml7yn2cpw+JsYzZlDXjGED/ScxlrNKFxLPzHWLk1NIkqTllNAVubTv1hqfiXFHDDu7dsXVd09FJdYKYreUckZYkiTpZLPw3VoPZmLcIYPOrk1LJdYKYveUcEZYkiRpOdP+3VoPZmKs4iqxg8ZkWEHUNHKMkSRJ3eIVWKaTibGKqsSuNCbDCqKmiWOMJEnqnqrfrf3eMFkmxiqqElvazNOetWveLL/Hpe3PkiR1VZevwDLL37XqZGIsoJyxnCWNyfCsXfNm/T0uaX+WJKmrunwFlln/rlUnE2MVpanxzlXOlK101s6zb+Ob9YpqaeP3JUlqWwnfl8b5Ttdm/E18b5j171p1MjFWUZoY71z1TNmws3aefavHrFdUSxq/L0lS20r5vlT1Ox3QavxNfG+Y9e9adTIxVlHGGe886Ixe1TNlw87aefatHpOYtXEaxhhJkjQLSvm+VPU7HdBq/Ct9b6jyHcbea6MzMVZxqozJGHa2r+qZsmFn7a5552HPvtWg7Vkbp2WMkSRJs6CUauU43+najn/Q94aq32HsvTY6E2PNhGFn+6qeKRt21m5+4ahn32rQ9qyNpZy5lmZNCWMIJZWnlGrlON/pSogfqn+H6ULvtbo+g0yMNROGnZF8/tM2VD5TNuisnWff6tPmrI2lnLmWZkkpYwgllaek70tVv9OVEv8432FmufdanZ9BJsaaCcPOSDZxpqwLZ99KVvUMdClnrqVZYk8MSYNMw/elUWIsIX6/wyyvzs8gE2PNhJXO9g07U1a1+0UTz6nRrPT/HvT+N3Xm2v+3umK5fX2cKoZtR5p901CtHBZjKfGXVH0vSZ29AU2MNROqnpFsogug3QqbN+z/Pcr7X+eZX//f6opB+/q6tdXG4Nl2JGl001B9n4Q6K+kmxpoZVc7oNdEF0G6F7Rj0/17p/a/7zG9p/28rcBrXai9998fXfapSFaO0tiNJpWuit+K0f2+os5JuYqxOa2IyJid4mqy23/+S/t9W4DSuKpe+uz9h54Wnr7qKUVLbkaRpNk2XsaxbnZV0E2N1WhMTGYzznIPO2k372bw2NfH+N7W9qlZb0bMCp1FVvfRdlZ4YK7WdYe2x6jpJmkWTuIxlSRXqunoDmhir05qYyKDqcw46a7d3/zFu/ejxqT6b16a6338Y/j63PRlGlYqeFTiNqqlL3y1nWNsZtp8DldZ5vJQ0q9q+jOWsVqhNjNVpJV3KadBZu9e+5W5++IVnDjybZ3Xkwep+/1c6a9r2ZBhVK3olcZ8tV5uXvhv2fFdcfffA/RyotG6c6of7rKSSNXUZy7p7qK30uEkfh02M1XlNTMNf5TkHnbU7cu/9A8/mlX7mbVLqfP9Hqba2eSmHNit6TXCfLds4l76rYtDzrdQeq64bpGqF2n1WUgmq9l6r2nOniQp1CcdhE2OpZYPOeA06a3fGaacMPJs3ibEhs6qkscJV42zqclR1jiGapn120meuJ6GUfWil9lh1XZX9EoZXobu4n0gqS9Vjd9WeO1UrzcMeN85xuC4mxlKLhp0NG3TW7uu/avDZvGveebjVsSGzrKSxwuOMaa6zotfEGKK2xzNVVcKZ60kpYR9aaT+vsm6c/bLkCockQfVjd5WeO8N6qFX5rjvKd9o25lAxMZZaNOxs2MtecCaw/Fm7m/fcu+zy+YWjlaqc44zxmFUrnW2t+z2p+j9oc0xz1RibGAfddqW5hDPXs6DJsfurXTdO9aOJCkdJM7pK0smq9lAbdqwd9l13pe+0bfTqMzGWWrRSVWLQWbtBy6tWOauO8Zj1L1+D3ucm3pNx/gdtjWmuGmMT46DbrjSPM8ZVn9PU2P0q68bZL+uucMzqjK6SZkfVHmpVv+uO01OoLibGUovqHsdatcpZdYxHV794NfGeTMP/oGqM44yDrrLPDntc1ar3OGNch+lalbCkmdLHHZ9fZ4Wjid4YTVWh29z32m4f09qupDZU7aFW9bg/bk+hOpgYSy1qYhxrlSrnOGM8uqiJ6wNPw/+gaowrVd/q3merVq+rbg+qnbnuYpWw7bH748RSpQpddT8Z1n6a2J+r7l/Q3jjpttvHNLcrqS1VeqiNc9yv2lOoLlORGEfENwFXAGuAP8rMyyccklRJKeNDxxnj0UVNVL2G7QtN/A+qVISqxlh1P6+6z1YdO1p1e4tWW/VuqkpYsjaPeZOIpYlKc9XeGE3sX9DMjNzLrWu7il769VSlaVXScX+1ik+MI2IN8DvATmA/8K6IuDYzPzDZyKRqShgfOiyOkio8pWjqPWnrfzBORahqjFX286r7bNWxoytV5qucuW6i2tdEj4U2tXXMG0UTsdRdaa7aG6OJKvTi7dU+Jww+rgxad8dd97VaRa/6uGGvrZT9XJq0ko77q1F8Ygw8CfhwZt4KEBFXA88BTIylIZoc49E1bb8ndW+viRmWm3hPqu6zVceONlGZb6LaV9I4XdWjid4YTVShF2/XOSP3oHWvfcvdrVbRqz5u2Gvr8uekNAumITE+G7h9yf39wIUn/1FEXApcCnDuuee2E5lUsKbGeIxiFttj22c/69xeUzMs1/2eVN1nq1avm+gJ0ET1uupM3otmsT3Ogrp7YzQ1J0AT1xxdbl2cUm1bTbSrpq+napuUyjQNiXEssywfsiDzSuBKgB07djxkvdQ1k6z82h7L0tQMy3Wrus+2/bhhmqhejxun7XG6NLE/j7t/1X3N0eXWffHWdcxtX99aFX2c3iTjHjNtk1KZIrPs9hgRXwP8bGZ+Y//+KwEy85cGPWbHjh25e/fuliKUFBELmbljuXW2x8kbNF5usSI0aJ3dAldv2Hs97P2s+rhBBrVJ22M31b1/rfScMPi4MmzdoFjabld1HzP9jJTKMqhNTkPF+F3A+RHxRcBHgRcAL5xsSJI0PapWhLR6JVWvpUUlzcg9yrq64m/ycbZVafYUXzEGiIinA79B73JNr8/MXxj29559k9rl2XCpLFaMpXL4GSmVZZorxmTmXwN/Pek4JEmSJEmz55RJByBJkiRJ0iRNRVfq1YqIu4Db+nfPAj4+wXAWlRIHGMsgxvJQo8bx6MzctNyKQtsjGMtySokDjGWQsdrkSe1xNc/XtFLiAGMZpJRYSokD/IxsSymxlBIHGMsg431GzmJivFRE7B40rqOLcYCxDGIszcdRyusCYyk5DjCWQWa1TZYSBxjLIKXEUkocMLvtEYyl5DjAWAYZNxa7UkuSJEmSOs3EWJIkSZLUaV1IjK+cdAB9pcQBxjKIsTxU3XGU8rrAWJZTShxgLIPMapssJQ4wlkFKiaWUOGB22yMYy3JKiQOMZZCxYpn5McaSJEmSJA3ThYqxJEmSJEkDmRhLkiRJkjptZhPjiPimiPj3iPhwRFw24Vj2RcT7IuK9EbG75W2/PiIORsS/LVl2ZkRcHxF7+r83TjCWn42Ij/bfm/dGxNNbiOOciPg/EXFLRLw/Il7WX976+zIklkm8L6dFxP+LiJv7sfxcf3kt74tt0vY4JJYi2qTtcTL8jBwai5+RtslW2R6HxuJn5Ky3x8ycuR9gDbAXeAywDrgZ+LIJxrMPOGtC234q8ETg35Ys+xXgsv7ty4BfnmAsPwv8aMvvyRbgif3bG4APAV82ifdlSCyTeF8COKN/ey1wE/DVdbwvtskHtmt7XD6WItqk7bHd//uSePyMHByLn5G2yVZ/bI9DY/Ezcsbb46xWjJ8EfDgzb83MY8DVwHMmHNNEZOb/Be4+afFzgKv6t68CnjvBWFqXmQcy893924eBW4CzmcD7MiSW1mXPkf7dtf2fpJ73xTaJ7XGQUtqk7bF7bJPLxlFEe1whltbZJptne1xeKW2yC+1xVhPjs4Hbl9zfz4T+cX0J/G1ELETEpROMY9EjM/MA9HZyYPOE43lpRPxrv9tKK11kFkXENuAJ9M40TfR9OSkWmMD7EhFrIuK9wEHg+sys632xTQ5me1yilDZpe2xVSe0RbJMPKKU9LhML2CabYnsczs/Ih8YBM9IeZzUxjmWWTfK6VE/OzCcC3wz8QEQ8dYKxlOb3gPOAxwMHgF9ra8MRcQbwZ8DLM/PTbW13xFgm8r5k5onMfDywFXhSRHx5TU9tm5wOE2uPUE6btD22zvY4mJ+Ry8dim2yO7XEwPyOXj2Nm2uOsJsb7gXOW3N8K3DGhWMjMO/q/DwK76HWZmaQ7I2ILQP/3wUkFkpl39nfs+4E/pKX3JiLW0mvUb87MP+8vnsj7slwsk3pfFmXmJ4F54Juo532xTQ7W+fYI5bRJ22P7CmuPYJsspj0OisU22Rzb42B+Rs5+e5zVxPhdwPkR8UURsQ54AXDtJAKJiNMjYsPibeAbgH8b/qjGXQtc0r99CfC2SQWyuPP2XUwL701EBPA64JbM/PUlq1p/XwbFMqH3ZVNEPKJ/+2HA04APUs/7YpscrNPtsb/dItqk7bF9BbZH6HibLKU9DovFNtkM2+NwfkZ2oD1myzOrtfUDPJ3ebGl7gZ+YYByPoTe74M3A+9uOBfgTet0ajtM7I/kS4AuAG4A9/d9nTjCW/wW8D/jX/s68pYU4nkKvi9K/Au/t/zx9Eu/LkFgm8b58JfCe/jb/Dfjp/vJa3hfbpO1xSCxFtEnbY/s/k2yP/e3bJh8aRxHtcYVYbJPNvEbb4/BY/Iyc8fYY/SeRJEmSJKmTZrUrtSRJkiRJIzExliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp00yMJUmSJEmdZmIsSZIkSeo0E2NJkiRJUqeZGEuSJEmSOs3EWJIkSZLUaSbGkiRJkqROMzGWJEmSJHWaibEkSZIkqdNMjCVJkiRJnWZiLEmSJEnqNBNjSZIkSVKnmRhLkiRJkjrNxFiSJEmS1GkmxpIkSZKkTjMxliRJkiR1momxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJExIRb4qI90bEupOWXxQRxyPiP08qNkmSusTEWJKkyXkp8AXAzywuiIiHA68HfjUz/6nOjZ2cgEuSpB4TY0mSJiQzPwm8GPixiHhSf/FrgUPAH0TE1RFxqP/zVxFx/uJjI+K8iHhbRHwsIu6JiHdHxDOXPn9E7IuIn42I10fEJ4E3t/PKJEmaLibGkiRNUGa+E/g94I0R8S3AtwOXANcD9wL/Bfga4ADwzohY33/oGcDbgZ3ABcCfAX8eEY87aRM/AnwQ2AG8qtlXI0nSdIrMnHQMkiR1WkQ8DHgPcD5wGfAJ4JXAY7P/QR0Ra4CDwH/PzGsGPM+/ANdl5s/37+8D3peZz2r8RUiSNMWsGEuSNGGZ+RngNcBngV8DtgNfBByOiCMRcQT4FLAROA8gIk6PiF+JiA/0u1ofoVcVPvekp9/d1uuQJGlanTrpACRJEgD3Afdn5v0RcQrwXuAFy/zd3f3frwG+CfhRYA9wFHgjcPIEW/c0Eq0kSTPExFiSpPK8G/g24OP9CbqW8xTgjZn5ZwARcRq9avKHWolQkqQZYldqSZLK82bgTuBtEfFfIuKLIuKpEfFrS2am/hBwcUQ8MSK+AngTcNqkApYkaZqZGEuSVJjMPAo8FbgV+N/0ZpW+it4Y40P9P/sRepNx/QO92an/pX9bkiStkrNSS5IkSZI6zYqxJEmSJKnTTIwlSZIkSZ1mYixJkiRJ6jQTY0mSJElSp5kYS5IkSZI6zcRYkiRJktRpJsaSJEmSpE4zMZYkSZIkdZqJsSRJkiSp0/5/r7Vdob7cfysAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x864 with 16 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(4, 4, figsize=(14, 12), sharey=True, sharex=True)\n",
    "fig.subplots_adjust(left=0.075, right=0.975, bottom=0.09, top=0.925, wspace=0.03)\n",
    "\n",
    "axes_flat = axes.ravel()\n",
    "\n",
    "for i, pop_ in enumerate(range(1, n_groups + 1)):\n",
    "    ax = axes_flat[i]\n",
    "    idx_ = data.index[data[\"pop\"] == pop_].tolist()\n",
    "    C_ = data.loc[idx_, \"C\"].values\n",
    "    year_ = data.loc[idx_, \"or_year\"].values\n",
    "\n",
    "    # Plot observed data points\n",
    "    ax.scatter(year_, C_, color=\"white\", ec=\"royalblue\", alpha=0.7)\n",
    "\n",
    "    # Add a title\n",
    "    ax.set_title(f\"Pop: {pop_}\", fontsize=12)\n",
    "\n",
    "fig.text(0.5, 0.02, \"Year\", fontsize=14)\n",
    "fig.text(0.03, 0.5, \"Red-backed shrike counts\", rotation=90, fontsize=14, va=\"center\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90eb2e90-eee6-48ee-b29d-89d510c78305",
   "metadata": {},
   "source": [
    "### 16.3 Random-coefficients model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf4167ed-1325-40ce-8116-23c8b471443f",
   "metadata": {},
   "source": [
    "#### 16.3.2 Analysis using PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0c837ca4-383d-41e6-ab6d-b0ec4b40ffc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = pd.factorize(pop)[0] \n",
    "groups = len(np.unique(idx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b9b03c2b-1f95-4d97-adea-74aa2d82004a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [beta, alpha, sigma_beta, mu_beta, sigma_int, mu_int]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [4000/4000 00:25<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 41 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_16_3_2:\n",
    "    \n",
    "    mu_int = pm.Normal('mu_int', 0, tau = 0.001)\n",
    "    sigma_int = pm.Uniform('sigma_int', 0, 10)\n",
    "    mu_beta = pm.Normal('mu_beta', 0, tau = 0.001)\n",
    "    sigma_beta = pm.Uniform('sigma_beta', 0, 10)\n",
    "    \n",
    "    alpha = pm.Normal('alpha', mu_int, sigma = sigma_int, shape = groups)\n",
    "    beta = pm.Normal('beta', mu_beta, sigma_beta, shape = groups)\n",
    "    λ = pm.Deterministic('λ', pm.math.exp(alpha[idx] + beta[idx] * year))\n",
    "    C_pred = pm.Poisson('C_pred', λ, observed = C)\n",
    "    \n",
    "    trace_16_3_2 = pm.sample(target_accept = 0.97, return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "59468c88-afb7-458f-b387-c8020f24227c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>mu_int</th>\n",
       "      <td>2.642</td>\n",
       "      <td>0.274</td>\n",
       "      <td>2.100</td>\n",
       "      <td>3.138</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2625.0</td>\n",
       "      <td>1215.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma_int</th>\n",
       "      <td>1.056</td>\n",
       "      <td>0.227</td>\n",
       "      <td>0.687</td>\n",
       "      <td>1.463</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.004</td>\n",
       "      <td>2634.0</td>\n",
       "      <td>1057.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mu_beta</th>\n",
       "      <td>-2.184</td>\n",
       "      <td>0.161</td>\n",
       "      <td>-2.486</td>\n",
       "      <td>-1.890</td>\n",
       "      <td>0.004</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1714.0</td>\n",
       "      <td>1731.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma_beta</th>\n",
       "      <td>0.570</td>\n",
       "      <td>0.163</td>\n",
       "      <td>0.299</td>\n",
       "      <td>0.881</td>\n",
       "      <td>0.005</td>\n",
       "      <td>0.003</td>\n",
       "      <td>1169.0</td>\n",
       "      <td>1475.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "mu_int      2.642  0.274   2.100    3.138      0.005    0.004    2625.0   \n",
       "sigma_int   1.056  0.227   0.687    1.463      0.005    0.004    2634.0   \n",
       "mu_beta    -2.184  0.161  -2.486   -1.890      0.004    0.003    1714.0   \n",
       "sigma_beta  0.570  0.163   0.299    0.881      0.005    0.003    1169.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "mu_int        1215.0    1.0  \n",
       "sigma_int     1057.0    1.0  \n",
       "mu_beta       1731.0    1.0  \n",
       "sigma_beta    1475.0    1.0  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_16_3_2, var_names = ['mu_int', 'sigma_int', 'mu_beta', 'sigma_beta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d3e6b68-102b-4116-8fe4-27e3efe9d9a2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}