{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e9b462e4",
   "metadata": {},
   "source": [
    "# Generalization to Bayesian Softmax Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8510f432",
   "metadata": {},
   "source": [
    "Ref: Chap 4 of Mar18\n",
    "\n",
    "https://cfteach.github.io/brds/referencesmd.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5015e78e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import theano.tensor as tt\n",
    "import seaborn as sns\n",
    "import scipy.stats as stats\n",
    "from scipy.special import expit as logistic\n",
    "import matplotlib.pyplot as plt\n",
    "import arviz as az\n",
    "import requests\n",
    "import io "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6c2e3cbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "az.style.use('arviz-darkgrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "224771aa",
   "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>sepal_length</th>\n",
       "      <th>sepal_width</th>\n",
       "      <th>petal_length</th>\n",
       "      <th>petal_width</th>\n",
       "      <th>species</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal_length  sepal_width  petal_length  petal_width species\n",
       "0           5.1          3.5           1.4          0.2  setosa\n",
       "1           4.9          3.0           1.4          0.2  setosa\n",
       "2           4.7          3.2           1.3          0.2  setosa\n",
       "3           4.6          3.1           1.5          0.2  setosa\n",
       "4           5.0          3.6           1.4          0.2  setosa"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_url = 'https://raw.githubusercontent.com/cfteach/brds/main/datasets/iris.csv' \n",
    "\n",
    "download = requests.get(target_url).content\n",
    "iris = pd.read_csv(io.StringIO(download.decode('utf-8')))\n",
    "\n",
    "iris.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f145d8b",
   "metadata": {},
   "source": [
    "## Recipe 1: Dealing with correlated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a3e9f488",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/r2/_2532dgx683084s9v9ss0cfc0000gq/T/ipykernel_21894/3442237557.py:1: FutureWarning: The default value of numeric_only in DataFrame.corr is deprecated. In a future version, it will default to False. Select only valid columns or specify the value of numeric_only to silence this warning.\n",
      "  corr = iris[iris['species'] != 'virginica'].corr()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[Text(0, 0.5, 'sepal_length'),\n",
       " Text(0, 1.5, 'sepal_width'),\n",
       " Text(0, 2.5, 'petal_length'),\n",
       " Text(0, 3.5, 'petal_width')]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbcAAAEoCAYAAADbp799AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1EklEQVR4nO3dd5xU9dXH8c/uKkivCigiiHKwRbBroqDEaEwzamJiL7HGBGKMxhI1scWOibHEWBOe51HTLInRqIANRIOIiBxDL9IWkA7L7s7zx+/ucnfZNlvmzs5836/XvJi5bc69rnPu73fPvb+CVCqFiIhILilMOgAREZHmpuQmIiI5R8lNRERyjpKbiIjkHCU3ERHJOUpuIiKSc7ZLOgCpatWqVbo3Q0RaXLdu3Qoau275kkH1/k4V9v600dtvDkpuIiKSli2p0nqXaZuBOOqi5CYiImkpJ/s7mJTcREQkLVtSZfUu0y4DcdRFyU1ERNKilpuIiOScMiU3ERHJNVtS5UmHUC8lNxERSUv2pzYlNxERSZO6JUVEJOdsyf7cpuQmIiLpKSPRh480iJKbiIikpVwtNxERyTUlreCZ+0puIiKSlvKUuiVFRCTH6JqbiIjknC2poqRDqJeSm4iIpEUtt2ZiZv2BOcDB7v5+cy2bCWZ2DnC/u3dMOhYRkeZQlsr+gpLsj7AVMbO5ZnZF0nGIiLSkLRTV+0paq2i5iYhI9mgNLbcGJTczOwq4A9gXKAMcOM/dp5nZEcBtwMHAKuB54Cp3XxOtOw6YAWwGzoo2+YdomfJomTOAkcBgYCMwHhjl7ouaYR8xs72BO4Gjou2/BvzE3ZdE858AegL/Bq4E2gN/B37o7huiZToADwInAeuB0cAXgWJ3Pyfaz92AO83sTgB3r+yYNrMRwH3AAGAS4fjNaY79ExHJpGxomdWn3vRrZtsBzwFvAfsDhxJ+2MvMbD/gFUJC25/wwz8EeKzaZk6Pvutw4CLgQmBUbH4b4IZoG18nJJr/bdQebRt/H+ANYBpwCPBloCPwnJnF9/9IQvL+MnAq8G1Cwq1wNzAsmn5MFOuRsfknAQuBXwF9oleFtsDVwHmEY9AVeKg59k9EJNPKUoX1vpLWkJZbZ8KP8QvuPiuaNgPAzJ4Cnnb3uysWNrNLgA/MbCd3XxZNXgz82N1TwAwzGwRcDtwD4O7xZDg72sYnZtbX3Rc2fvcAuAT40N2visV4FrASOIjQigJYA1zs7mXRdz8LjABuM7OOhMR0lrv/O9rG+YRkRrQPK82sDFhb0SKM2Y7QCvRo3buAx8ysIDomIiKtRnkrKNeoN7lFP9pPAC+b2WuELr0/u/t84EBgDzM7NbZKRVfcQKAiuU2s9iM+AbjJzDq7+xozO4DQchsCdI9tox+xBNJIBwJHmdm6GuYNZGtymx4ltgqfEVqpFcttH1sWd19vZtMaGMPmisQW23YboBshyYqItBoluXKfm7ufa2ajgeOBbwK3mNmJhK7GPwD31rBag66XRdeyXgZeBc4kJMSewJuEBNBUhcA/gJqqGJfG3m+pNi9F81WTltawbZpx+yIiGVOeBd2O9WlwtaS7fwh8CNxuZi8BZwOTgX3cfWY9qx9arQvuMOCzqNV2ICGZXVNRYGFmJ6W7I3WYDHwXmOfu1RNYQ80iJL+DgdkAZtaecI1uVmy5EmgFV1pFRJqgrBWcl9eb3MxsAKEI5HlCa2x34AuEysHngYlm9hDwMLCWUPH4DXe/KLaZnYHRZvYAsB/wM+DmaN58QiXlZWb2O2Av4Kam71ql3wEXAE+b2e3A8mgfvgv81N3X1rcBd19nZo8REnsx4RridYSWV7y7dS5wpJn9idAVWdyM+yEikhVaw+O3GpJ+NwCDgGeBT4EngTHA7e4+lVBe359Qvv8h4baApdW2MYbQonkXeAR4lKgr092XE1qBJwLTCdfeLm/8LlXl7p8RSvbLgX8BHxMS3ubo1VBXELpKnwfGAlOB94FNsWWuB3YltOaWNzV2EZFs1BqqJQtSqZYt1ovu/5rm7pe16BdlmJm1BeYBd8arRZtq1apVqp4UkRbXrVu3Rj8g8g+fHlnv79QPBr2Z6AMo9YSSBjKzoYQu00lAJ+Cq6N+nk4xLRCTTSlLZnzqyP8Jqout7Z9Qy+0/ufnELfv3lgBGqH6cARzXDfXgiIq2KBisF3H14M2/yeuCuWuataebvquTuHxBu+hYRyWs5US2ZbaKnniyrd0EREWkRraFastUlNxERSVZO3cQtIiICGolbRERy0Jby7E8d2R+hiIhklXK13EREJNdkwxNI6qPkJiIiaVG1pIiI5BzdxC0iIjmnKS03M7uUMDJMH8KD7Ee5+5t1LH8acCXhAf5rCGN/XuHuS+r6nuzvOBURkaxSniqo91UTMzsVuA+4FRgKvAO8ZGb9aln+i8AfCaPR7EMYPWZvwkgzdVLLTURE0tKEm7gvB55w90eizz8ys+OBS4Cra1j+cGChu98bfZ5jZr8FflvfFym5iYhIWrY0IrmZWRvgQLZ9NvArwBG1rPY2cKuZfQN4EegBfA/4Z33fp25JERFJS3mqsN5XDXoSBq2uPpj1UqB3TSu4+wRCMhsDlBAGgS4gDHBdJyU3ERFJSzkF9b6ag5ntTeiCvInQ6juekAgfrm9ddUuKiEhatpQ3qlqyGCgDelWb3guorfLxamCSu98ZfZ5qZuuBN83smrrG01RyyzLf7XFh0iFkrWdW/D7pEESExt3n5u4lZvYf4Fjg2disY4G/1LJae0JCjKv4XGfPo5KbiIikpQndjvcAfzSzSYRikYuBnYGHAMzsKQB3Pyta/gXgETO7BHiZcG/caGCyu8+v64uU3EREJC2ljeuWxN2fNrMewHWERDUNOMHd50WL9Ku2/BNm1gm4DLgbWA28DlxV33cpuYmISFqa8vgtd38AeKCWecNrmNag+9qqU3ITEZG0aMgbERHJOaXl2X8XmZKbiIikRaMCiIhIzlFyExGRnFOqkbhFRCTXqOUmIiI5RwUlIiKSc1JquYmISK7RfW4iIpJzytQtKSIiuUYFJSIiknN0zU1ERHJOWXn2J7fs7zitg5n1N7OUmR3UAtt+wsxerGeZF83siXqWOcfM1jVrcCIiCSqnoN5X0tRyq91ISO+/kJnNBe5397taJCIRkSygbslWzN1XJx1DSzrouCFcOvpcCosKeenR13j69r9XmX/yT77OV88fQVlpGauXr+Gu8x9g2fxiAG7957XsddieTHtrBr/45q8TiF5EktQauiWblNzM7CjgDmBfoAxw4Dx3n2ZmRwC3AQcDq4DngavcfU207jhgBrAZqBhS/A/RMuXRMmcQWlCDgY3AeGCUuy9qRKwTgb+7+6+jz38CTgf6uPsSM2sfxTnC3d+Kuht7uvvXo+XbEwbYOwVYD9xXbfvjgN2AO83sTgB3L4jNHxGtMwCYFB2nOenuR3MoLCzkR/efz1VfuYnihSu5f9JtTHj+feZ/srBymZkfzOGHB1/F5o0lfP3ir3DB7Wdyy/fvBeDZu56jbfu2fO3CY5MIX0QS1hpabo2+5mZm2wHPAW8B+wOHAqOBMjPbD3iFkND2B04ChgCPVdvM6VEMhwMXARcCo2Lz2wA3RNv4OtAT+N9GhjwOGB77PAwojk07AiglJJ6a3AUcC5wMjACGAkfF5p8ELAR+RRg+vU9sXlvgauA8wr52BR5q3G40nR2yB5/NXMKSOcso3VLKuKff5ohvVb1s+eG4j9m8sQSATyZ+yo59u1fO++D1aWxYuzGjMYtI9kilCup9Ja0pLbfOhB/pF9x9VjRtBoCZPQU87e53VyxsZpcAH5jZTu6+LJq8GPixu6eAGWY2CLgcuAfA3ePJcHa0jU/MrK+7LyQ944DLoqTcH+gC/AY4Gvg/QpKb4O4l1Vc0s47A+YTW1svRtHMJyYwo1pVmVgasdfcl1TaxHfBDd/do3buAx8ysINr3jOq5S3eWL1xR+bl44UoGH7pnrct/9fwRTPrXB5kITURagZzulox+zJ8AXjaz14DXgD+7+3zgQGAPMzs1tkrF0RgIVCS3idV+3CcAN5lZZ3dfY2YHEFpuQ4DusW30I5ZYGugtQgvqYGCf6POrwMPR/OHAv2pZdyChFTmhYoK7rzOzjxr43ZsrElvks2h73YCVDdxGIkacfiSDDtydnw6/IelQRCRLZEPLrD5NuhXA3c8ldEe+AXwTcDM7LtruHwhJqeK1P7AnMKUh2zazDsDLwAbgTEJSOj6a3aYRsa4D/kNoqQ0HxgITgX5mtke0/XHpbreBSqt9rkjoidyKUbxoJTv27VH5uWff7hQvWrHNckNH7Mdp15zE9d+6nS0l1XdBRPJVrndLAuDuHwIfAreb2UvA2cBkYB93n1nP6odW65o7DPgsarUdSLjGdk1F4YWZndTEcMcRkttg4D5332Rm7wLXUvf1tlnAlii+2VEsHQiFNLNiy5UARU2MscX5ezPZZc8+9O6/E8WLVjL81C9y2+lV6mMYOKQ/ox66kGu+egufL1+TUKQiko1y+vFbZjaAUATyPLAI2B34AvBgNG2imT1E6PZbS0go33D3i2Kb2RkYbWYPAPsBPwNujubNJ1RSXmZmvwP2Am5qbLyRccBPCa3BybFp1wLja7reBpVdkI8SEvhyQrfi9WybyOYCR0aVmJvdvbiJ8baI8rJy7v/Ro9z2r2spLCrk5cfHMm/6Qs7+5al8+v4sJrzwPhfecSbtOu7AL575KQDL5hdz/Ym3A3DP+F+x6+BdaNdxB/5n/kPc84MHef+VD5PcJRHJpIxXCqSvKS23DcAg4FlCC2spMAa43d23RLcJ3Ewo3y8itHj+Vm0bY6J57xIO16PAvQDuvtzMzgZuBX4ITCUUm9R2Xawh3or+fdPdy6L34wjX9cbVs+4VQIdoHzYAv40+x11PSOazCNf3svb0ZtJLHzDppapFIk/e8HTl+6u+Uvt5xOXDrm+xuEQk+5W3goKSglQqmRQc3Rc2zd0vSySALHVs4XdawTlRMp5Z8fukQxDJGd26dWt0hhr49C31/k7NOvXaRDOgnlAiIiLpyeVrbtkmur53Ri2z/+TuF2cyHhGRXJUqTzqC+iWW3Nx9eDNv8nrCU0RqonI/EZFmkg2l/vXJmZZb9NSTZfUuKCIiTdMKKgNyJrmJiEhmpFpBtaSSm4iIpEnJTUREco26JUVEJOeoW1JERHJNQs/+SIuSm4iIpEfJTUREck2BuiVFRCTnqOUmIiI5R08oERGRnNOEZ0ua2aWEsTv7AB8Do9z9zTqWbwNcB5xJGAN0KXCXu/+mru8pbHyIIiKSl1INeNXAzE4F7iOM0zkUeAd4ycz61fFt/wccD1wIGPAdwviedVLLTURE0tKEgpLLgSfc/ZHo84/M7HjgEuDq6gub2VeAEcBAdy+OJs9tyBcpuYmISHoaUVASdS8eyLajt7wCHFHLaicC7wGXm9lZwEbgJeAad19X1/cpuYmISCb0BIoI18zilgJfrmWd3YEvAZuBk4GuwG8J195OqevLlNyyzL6Ti5IOIWsdMOaxpEPIWpNPPy/pECSPZPA+t0JCO/E0d18NYGaXAS+bWS93r54oq6woIiLScI0rKCkGyoBe1ab3ApbU8k2LgUUViS3ySfRvXUUoSm4iIpKmRiQ3dy8B/gMcW23WsYSqyZq8DexsZh1j0wZF/86rK0R1S4qISFoKGn+f2z3AH81sEiFxXUy4fvYQgJk9BeDuZ0XL/w/wC+BxM7uRcM3tPuDP7r6sri9Sy01ERNLTyPvc3P1pYBThpuwphGKRE9y9ohXWj1h3Y1QR+WWgC6Fq8hlgPFDvRWa13EREJC0FTXi2pLs/ADxQy7zhNUxz4Cvpfo+Sm4iIpEejAoiISK5pSsstU5TcREQkPUpuIiKSa5pQLZkxSm4iIpIetdxERCTX6JqbiIjkHiU3ERHJNWq5iYhI7lFBiYiI5Bq13EREJPcouYmISK5pDfe5tYpRAcwsZWZ1DinemGVbmpkNj+LpmXQsIiLNppGjAmRSxlpuZjYcGAvs6O7FmfreTDGzccA0d78s6VgaYvmUNUx/YiGp8hS7HtODgSf2rjJ/Y3EJH/5uHqUbykiVp7DTdmanoV0oWVvK5HvmsHrWBvoO784+5+2a0B5kxlH9+nPDkUdTWFDA09On8dDkSdss87U9BjHykCNIpVJ8smI5o175ZwKRimSOrrlJVkqVp/j4sQUccu0e7NBje96+2tnpoC506tuucpmZf11Cn8O7sttXdmTtwo28/+vZ7HR/Fwq3L2DQqX1Yu2AT6xZsTHAvWl5hQQG/GjaCM5/7M0vWreW5757Oq3NmMnPVyspl+nfpyiUHHsopf/lf1mzeTI927erYokiOaAXdkg1OblHLZAawGagYJfUPwFXuXm5mbYCbgNOB7sDHwHXu/rKZ9Se02gCWmxnAk+5+jpkdD1wL7EtozL4HjHL3T5q4bxVx7wLcDRwXTXon2v5/o/k3AqcANwO3ADsBrwE/qGhhmtl2wJ3AOdE2ngB2APZy9+Fm9gQwDBhmZj+MlhkQC2N/M7sV2A+YDlzo7pObY/8a4/OZG2jfqy3te7UFoM8R3Vj63uoqyQ2gdGP4Cy7dUE7bbtsDsN0ORXQf3JENSzZnNugE7N+rN/NWf86CNasBeOG/zrG778HM/2xtvX1vny/wx4+msGZzOB4rNuZ2wheB1tFyS/ea2+nROocDFwEXEkZVBXic8AN/GiFRPQm8YGb7AwuAk6Pl9gH6ACOjzx2A0cAhwHBgdbRem3R3pjoza09Iqpui2A4HFgOvRvMq9AdOBb5NGBRvKCHRVbiCkNh+ABxGOAanxeaPBCYQjkGf6LUgNv824OfAAcAKYIyZJTYg0qaVJezQY+vhbdejDZtXbamyzJ7f6cOiN1fy+iXTeO/Xs9jn3L6ZDjNxvTt0ZPHatZWfl6xbS+8OHassM6BrNwZ07cazJ3+Pv57yfY7q1z/DUYokIAevuS0GfuzuKWCGmQ0CLjez54DvA/3dfX607P1m9mXgIne/1Mwq+nKWxa+5uftf4l9gZucCawjJ7q30d6mK7wEFwLlRzJjZRcAy4OuEIcshHIdz3H11tMzvgXNj2xkJ3F4Rq5mNAo6P7cNqMysBNrj7kti+VLz9hbuPjab9KtqvXYCFTdy/FvPZ26voO6w7u3+jF6s+Xc+H98/jyLsGU1CY/YMUZlJRYQH9u3Tl+397ht4dOvL0Sd/j+P99krUlud+ylfzVGqol001uEyuSRGQCoSvyS4QkMj32gw7QFni9rg2a2cBoG4cCOxJaRYVAvzRjq8mBhO7BtdXiag8MjH2eV5HYIp8Ruicxsy5Ab6CyL8rdU2Y2CWhoNcXUatsm2n4iyW2H7m3YtKKk8vPGFSWV3Y4VFo5dwcFXh0PUbVAHyraUU7K2lLZdqi6Xy5asX0efTp0qP/fu2Ikl69dVXWbdOqYsXUxpeTkL165hzucrGdC1K1OXLc10uCKZkwUts/o0560AKeBgYEjstRdwXj3rvUhIahcREtxQoBRocrckYf+mVItpCDAIeDi2XNU+ubAvzXls4tuv+LNI7DaMLgPbs37JZjYs20x5aTmL31lFr4O6VFmmXc/tWTEtdMmtW7iJ8i3ltOmcX/VHU5cuoX+XrvTt1JntCwv5xp7Gq3NmVVnmldkzOWyXcI7TbYd2DOjanflrVte0OZGcUdCAV9LS/bU61MwKYq23wwgtkQmE/eld0f1Wg4qmQlHFBDPrAQwGLo112x3QiLhqM5nQXVrs7p83ZgNRl+MSQuJ+PYqxIPq8JLZoCbF9y2aFRQXsc15fJt06C8pT9B3eg067tuPTZxbTZff29DqoC4PP3IVpDy9gzj+WQUEBX7hkNwoKwp/s2Ms+pnRDGeWlKZa+t5qDrx24TTFKLihLpbjhjdd56lsnU1hQyLPTp/HflSv4ySFH8NGypbw6dxZvzJ/Lkf1245XTzqEsVc5t74zn802bkg5dpEXlYrfkzsBoM3uAUPn3M+Bmd//UzMYAT5jZTwlJpTuhQGS2u/8VmEdotXzNzF4ANgKrgGLgAjNbQLgOdSeh5dYcxhCKQZ4zs+uB+YSuxG8BD1VUTDbAfcCVZvYpodrxIkLRyOLYMnOBQ6LK0HXASrLYTkO7sNPQqq21Qd/tU/m+U992HH7ToBrXPfr+fVo0tmwybt4cxs2bU2XavZPeqfL5lrfGcwvjMxmWSLJysFtyDKF18i7wCPAocG8071xCteAdhFsGXgSOIiQ13H0RcAOhCnEpcL+7lxOqFL8ATAN+B/yCcLtBk7n7hiiG2cCzUVxPAt0IibWh7gL+SNi/idG0vxGqMOPLlBCS33Ka55qhiEj2aQXVkgWpVMOiaG1P4GhpZvYB8Ja7/6g5t/uTKd/Lgj+L7PT3tw5OOoSsNfn0+i5ti1TVrVu3Rl8aG3LZvfX+Tk25/yeJXnrLrwqBRjKz3Qg3gY8HtgcuILQ2L0gyLhGRJLSGm7hbVXIzs2uAa2qZ/aa7f7WFvrqc8FSWOwldudOBr7r7+y30fSIiWSunCkrcfXgLxtFQD7H1xuvqWuy5R+6+gHAvn4iIqOXWvNx9JVlehSgikvOU3EREJNfkVLekiIgIQEEDq+yTpOQmIiLpyf7cpuQmIiLpUbekiIjkHN3nJiIiuUfJTUREco26JUVEJOeoW1JERHKPbgUQEZFc05RuSTO7lDAWaB/gY2CUu7/ZgPW+BIwDZrj7vvUtn+54biIikucKyut/1cTMTiUM/nwrMBR4B3jJzOoc/9LMugFPAa81NEYlNxERSU/jByu9HHjC3R9x90+i8TAXA5fU842PEgaantDQEJXcREQkLQXlqXpf1ZlZG+BA4JVqs14Bjqjtu6JuzF7AzenEqOQmIiJpKUjV/6pBT6AIWFpt+lKgd00rmNl+wA3AGe5elk6MKijJMp9vaZd0CFmraIPOxWrTZfOhSYeQtVa3fTfpEHJOJu5zM7O2wNPAFe4+J931ldxERCQ9jbsVoBgoI3QxxvUCltSwfB9gL+BxM3s8mlYIFJhZKXCCu1fv4qykU2EREUlLY7ol3b0E+A9wbLVZxxKqJqtbBOwHDIm9HgJmRu9rWqeSWm4iIpKWJnRL3gP80cwmAW8DFwM7E5IWZvYUgLuf5e5bgGnxlc1sGbDZ3atMr4mSm4iIpKeGasiGcPenzawHcB2h23EaoXtxXrRInfe7pUPJTURE0tOEp2+5+wPAA7XMG17PujcCNzbke5TcREQkLTXdx5ZtlNxERCQtGhVARERyj5KbiIjkmoKy7M9uSm4iIpKWAo3nJiIiOSf7c5uSm4iIpEfVkiIiknvULSkiIrkmE6MCNJWSm4iIpEfdkiIikmtaQ7Vkqx7yxsxSZnZKC2z3RjOr86nTZna/mY2rZ5nhUYw9mzVAEZEklaXqfyUs8eSWpQngLmBYOiuY2Tgzu7+F4hERyRoFqVS9r6SpW7IG7r4OWJd0HC1p1YermPPHOVAOOw3fib7f7Ftl/pw/zWH19NUAlJeUs2XNFg79/aEAzP2/uayasgqAXU/clZ6HZdN5SdMdNWA3rhsxnKKCQp6ZOo2H332vyvyT9t2bnw8/kiVrw5/Inz74kGemTuOwfn255uit50QDe3Rn5PP/5NWZszIaf0t681249bdQXg6nfA0uOL3q/EVL4LrbYeXn0KUz3HEt9N4pzLvzQRg/EVLlcMRBcM2PoaAg47sgzSELkld9mpzcoq65GcBm4Kxo8h+Aq9y93MzaADcBpwPdgY+B69z9ZTPrD4yN1lluZgBPuvs5ZnY8cC2wL+GWwfeAUe7+SSNi/D/gc3e/OPp8c7Ttw919YjRtAXC1u//JzG4ETnH3faN5RcDtwPnRJp8EimLbf4LQ0htmZj+MJg+IhbC/md1KGFV2OnChu09Odz+aS6o8xewnZ7PPz/ehTfc2TL1+Kt0P7E77XdpXLjPgjK3hL35lMevnrgdg5QcrWT93PUNuGUL5lnKm3TKNrl/oynbtc+M8qbCggBu/fAxnP/NXlqxdy1/POo3XZs5i5oqVVZb7x4xP+eWrY6tMmzh/Id98cgwAXXZoy2sXnMdbc+eRK8rK4KbR8Ojd0GtH+O5FcPQXYY/+W5e58wH41nFw4vEwcTLc83u44zr4YFp4PfdYWO70y+C9KXDI0AR2RJouC7od69Nc3ZKnR9s6HLgIuBAYFc17nPDDfxohUT0JvGBm+wMLgJOj5fYhDF43MvrcARgNHAIMB1ZH67VpRHzjom1UGA4UV0wzsz2AvtFyNfkpcEG0b4cTElv8nHUkMIGwr32i14LY/NuAnwMHACuAMWaW2DnrulnraNerHTvstAOF2xXS87CerPzPylqXL55QTM/DQ+ts46KNdLbOFBQVULRDER36deDzqZ9nKPKWt3+f3sz7/HMWrF7NlvJy/vGJ8+U9Bqa9neNtEOPnzGFTaWkLRJmMqZ9Av11g152hzfZwwjHw+ltVl5k5Dw49ILw/dCi8/vbWeZtLYEsplGyB0jLo0S1zsUvzyqduycXAj909Bcwws0HA5Wb2HPB9oL+7z4+Wvd/Mvgxc5O6XmlnFr+oydy+u2KC7/yX+BWZ2LrCGkOyq/S9Vr3HAg2bWh5AkDwauB44Bfk1IcrPcfWEt648C7nD3Z6JYRgLHxWJdbWYlwAZ3XxKLueLtL9x9bDTtV1H8uwC1fV+L2rxqM226bz1HaNO9Detm1dwLu6l4E5uWbaLLPl0A6LBbBxb8dQE7n7Az5SXlrJ6+mna7tMtI3JnQq2NHFq9dW/l5ydp17L9z722WO27QnhzcdxfmrvqcW14fx+K1VY/f1wcP4rH3E2uct4hlxVu7GCG03qZW60cZPBD+/QacdQr8+01Yv6GAVatTDN03JLujTgo9Wqd/Gwb2z2j40pyyIHnVp7mS28QosVWYQOiK/BJQAEyP/dADtAVer2uDZjYw2sahwI6ElmEhjRiG3N1nmNkSQhJbDswCngZ+YWbbR9PH1RJHF0JLbEJse+Vm9i6wawNDmBp7/1n0704klNzSUTyhmB6H9KCgMDQ0u+7XlXWz1/HRLz9i+87b02nPTpXz8sXrM2fz4idOSVkZ39t/P+444TjOfHrrudiOHTpgO/bkzTm50yXZUFdeGrou//4SHLQ/9NoxRVEhzFsIs+bB2GfDcuf/FN7/MCwjrVB59t/FnYkLJSlCS2lLtekb61nvRcKP/0XAIqCUcL2qMd2SAOOBo4FlwFh3n2tmxVFsw4CrG7ndhojve8VJQGKVqm27taVkZUnl55KVJbTpVvNhXTFxBQPOHlBlWt9v9aXvt0IByqe/+5R2vXOn5bZ03Tr6dOpU+bl3p44srdYq+3zTpsr3z0ydxlXDj6wy/4TBg3jlv7MobQU/AOnYqScsWbb189Ll0Kvntsv89ubwfv0GeOUN6NwJnn0R9t8bOkSXdY88FKZ8rOTWarWCP+3m+oE9tNo1pMMILZQJhJZbb3efWe21KFq24lc2XqDRAxgM3Orur0ZFJJ1oWjIeR0huw9naShtHuJZW6/U2d19N6HY9LBZfAaF7NK4kvg/ZrOPuHdm4ZCOblm2ivLSc4onFdD+g+zbLbfhsA6XrS+m059Yf+1R5ii1rQ65eP3896xesp+t+XTMVeoubungJu3XrRt8undm+sJCv7WW8NnN2lWV27NCh8v2IPXZnVrVik2/sZbz4yYyMxJtJ+w0OLbCFi8N1s3++HgpK4lZ9vvWk/pExcNJXw/s+veC9D6G0NFx3e/9DGLhbRsOXZpRP19x2Bkab2QOEisCfATe7+6dmNgZ4wsx+CkwmVEwOB2a7+1+BeYTWzNfM7AVCi24VoeDjgqiKcRfgTkLrrbHGAQ8Cu1E1uT1C3dfbAO4DrjazT4GPgEsJXZWLY8vMBQ6JKkDXAbVXaCSsoKiA3c/enel3TCdVnqLXsF6079ue+X+eT8cBHel+YEh0xROK6XlYTwpi9dqp0hTTbgr3txe1K2LQJYMoKMqdbsmyVIpfvvo6j3/nJIoKCnj2o4/574oVjPzS4UxbspTXZs7m7AOHMGKPgZSWl7N60yau/OfLlevv0rkzvTt14t35Wd/jnLbttoPrRsEPrggJ7KQTYM8B8JtHYd/BcMwXYdKUUCFZUBBaZdePCuseNwzenQzfOjfM+9Ih2yZGaUXKsr/p1lzJbQyh1fIuIVE9CtwbzTuXUHZ/B6GFtBKYRHQLgLsvMrMbgFsItxA8Fd0KcCrwG2AaMJNQsVilyCQdsetuK9x9eTR5HOEYjKtn9buB3lF8AH+M9nmv2DJ3ESpBpwPtqHorQNbpNqQb3YZULVfrd0rVy5n9Tt728mZhm0KG3pHb9dvjZ89l/Ownqky7763KS67c9cbb3PXG29Rk0Zo1fOnBR1oyvEQNOyy84n58/tb3xw0Pr+qKiuCXV7RkZJJRWdAyq09BqolBRve5TXP3y5olojx37nvnZv9fTULeHPuFpEPIWp+e9WDSIWSt1W3fTTqErNStW7dGd7l8dc8r6/2deum/dyTapZMbd96KiEjm5FG3ZKLM7Brgmlpmv+nuX81kPCIiOS2VB8nN3Yc3QxxN9RDwTC3z6rvlQERE0tEKrrnlRMvN3VeSxdWJIiI5Rd2SIiKSc9RyExGRnFNWlnQE9VJyExGR9KjlJiIiOUfJTUREck1K3ZIiIpJzytVyExGRXKNuSRERyTnqlhQRkVyTagUD8Sq5iYhIetQtKSIiOUfdkiIikmtSqpYUEZGckw9D3oiISH5pDTdxF6RawYVBERGRdBQmHYCIiEhzU3ITEZGco+QmIiI5R8lNRERyjpKbiIjkHN0KIDUys65UO/lx95XJRJNddGxEsp+Sm1Qys92Ah4DhQJvYrAIgBRQlEFZW0LGpn5ntDOzEtol/cjIRZQ8dm8xTcpO4x4GuwPnAZ4QfbQl0bGphZkOBPwGDCck+Lq8Tv45NcpTcJO4Q4DB3n5Z0IFlIx6Z2vwcWABegxF+djk1ClNwkbg7QNukgspSOTe32Boa6+6dJB5KFdGwSompJiRsJ3GZmeyQdSBbSsandR0DvpIPIUjo2CdGzJfOcma2lalfJDoTrAJuB0viy7t45g6ElTsemdmbWPfZxCHArcB3hx3xLfNl8qyTVsckO6paUy5IOIIvp2NSumKqJvwB4pYZp+Vg0oWOTBZTc8py7P5l0DNlKx6ZORycdQBbTsckC6paUSmZWBvRx92XVpvcAlrl73p5l6tjUzsz6AQvcPVVtegGwq7vPTyay5OnYJEcFJRJX/T6cCm2BkkwGkoV0bGo3B9ixhundo3n5TMcmIeqWFMzs8uhtCrjYzNbFZhcBRwIzMh5YFtCxaZCK60fVdQQ2ZTiWbKNjkxAlNwH4UfRvAfADID6GfAkwF7g4wzFlCx2bWpjZb6K3KcJtEhtis4sIN75PyXRc2UDHJnlKboK7DwAws7HASe6+KuGQsoaOTZ32i/4tAPaiavdsCTAZuCvTQWUJHZuEqaBERJrEzB4HRrr7mqRjyTY6NslRcpNKZvZYLbNShOsDM4Gn3f2zzEWVnDqOxzbc/byWjEVE0qNuSYnbkVAgUQ5UPCB4X0LXyn+Ak4BfmdmR7j4lkQgzq3qV21GEY/NR9HlfQsXxG5kMKtuY2eu1zIqfFD2ZL8O7RF3YDWo1uPsxLRxO3lJyk7i3gXXA+e6+AcDM2gOPAB8CJwBPAXcDI5IKMlPc/RsV783samAjcK67r4+mdQAeZWuyy1czgNOAJcCkaNrBhGcq/p1wwnSpmR3v7q8lEmFmxUeOKAJOJxybd6NphwB9CEPhSAtRt6RUMrPFwDHu/km16XsDr7l7n2h8qlfdvUciQSYkOjYj3H16ten7EI5N3j4c18zuAQrdfVS16XcDKXe/wszuAw5x98OTiDEpZnYvIcGNjN/IbWajgQJ3H5lUbLlON3FLXEfCGWV1vaN5AGvIzxZ/R2DnGqb3AdpnOJZsczbwuxqmPwycG71/hDD8S745C7i/+hNKgAeAMxOIJ2/k44+U1O5vwKNmdiXwXjTtYOAO4K/R50OAfByb6i/A42b2M2BiNO0w4Ha2Hpt8VQDsA/y32vS92fpklxLC9cp8U0C4LaD6/zP71bCsNCMlN4m7GLiHcC2g4m+jFHgMuCL6/AlhVOF8cwnhWuMTwPbRtFLCNbcralknXzxJOCnak6onRVcRjhfAMKpei8oXjwF/iI5N/KToSuDxxKLKA7rmJtuICiUGRh9nVRRQiI5NTcysCPgZ8GO2Dsy5BLgPuMvdy6IHCJe7+8KEwkyEmRUSTn5GsrXLfzHh2Nzt7mW1rStNo+QmIs3GzDoD6KblbenYZJaSm1Qysx0IZ5gjgJ2oVnDk7l9IIq6kmNnzwBnuviZ6Xyt3/2aGwhKRBtA1N4l7APg28CzwDg28ETWHrWDrMViRZCDZzMy6A7dQ+0lR5yTiSoqZTQWGufsqM/uIOv4/yrcTxkxScpO4E4HvuPurSQeSDdz93JreyzYeBYYCvwc+QydFfwE2x97n+/FIhLolpZKZLSTcqOxJx5JtzOwIYJK7lyYdS7YxszXAse7+br0Li2SIbuKWuDuAy82stlGn89nrwCoze8XMrjGzI8xMPR/BMsJj26QaMzvNzGp6MIK0MLXcpJKZvUB4DuBqYDqwJT4/n4smzKwd8EXC/VrDCfdxbQEmAGPd/bbkokuWmZ0KfBc4292V5GLMbD6wCzALGFfxypeRNZKk5CaVorGnaqXrTluZ2UDgWuAMoMjdixIOKTFR0UR/wjMU57HtSVFeF02Y2R6EE6Jh0asi2Y1194sSDC2nKbmJNICZ7UT4gTo6+rcf4Qn44whn4uOTii1pZnZDXfPd/ZeZiiWbRTe7H0J4wk/enxS1NCU32YaZHUR4CseL7r4+eirH5nwupjCzcmA54WHArwLvuvvmuteSfGdmh7D1pOiLQDEwnq0nRfMSCy7H6YK4VDKzXsBzhLPLFLAnMJvwvMlNhBu889X/EAYrHQkcAIw1s3HA5Bqe+J53ogcAfJ1wUvSwu38edd2ucveVyUaXqImEk6K7gIvcfX7C8eQNVUtK3L3AUqAHsCE2/VngK4lElCXc/Qx370dIbH8DhhBGA1hpZs8lGVvSomtKM4CHCDdzd49mXUKowM1ntxJGBLgJ+KeZ/dbMTjazvBoPMQlKbhI3ArjW3VdVmz6LcI1JYA7h6fbTAQc6AMcnGlHyRgOvAL0Io5VXeJ7QHZe33P06dz8S6EZo9X8e/bvIzD5MMrZcp25JiWtHGHeruh0J3ZJ5KxrjbjjwJaAt8B/CtZO7gbeSiywrHAEcFj39Pz59PjUP8JqPOgM9CY8n6w20iT5LC1Fyk7g3gHOAa6LPqajC6yrgtaSCyhLfJhQB3Ae8paFutrF9DdP6Ee6ZzFtm9iCh/N8IXf4VJ0Tj9CSglqXkJnFXAuPN7GBC6+RuwgjLXQiVXnnL3Q9vyHJm9gBwvbsXt3BI2eQV4HLg/OhzKhre5ZfAPxKLKjt0JZwQKZllmG4FkCrMrDehEOBAwjXZycDv3H1xooG1EtFzFoe4++ykY8kUM9sZGBt93B34ANiD0FI5yt2XJxVba2Fm/wB+oP/Pmo9ablKFuy8B6rwpV+qUd8/ldPfPzGwI8H1CNWkhYYSAMe6+sa51pdJRhGve0kyU3PKcmR3Q0GXdfXJLxiKtV5TEHoteIolTcpP3CTds19fiSBGeHSiCmZ3U0GXd/a8tGYtITZTcZEDSAUir9OcGLqeTIkmEkluea8yz7fK0IlBi3F0PgJCspj9QaYwzCDelyrb+BKxJOohsZGb/0MCdkilquUlj5EVFYGOKbdz9kpaLqNVTRWDtbgXy+QHTzU7JTaR2KraRtDWm2CafR3JvKUpuIrVTsY00hoptsoCSm0gtNJCkNIaKbbKDkptIGqJHTfUjPNW9kru/kUxEIlITJTdpjLyrCIySWsVo3BXX4eIPZlX3ktTIzLYjjG5f00nRU4kElQeU3PKcKgIbbDRQBuwNvEcYoLQX8CvgJ8mF1arkXUWgmQ0GXiBcvy0g/A1tB2wBNgNKbi1EyU1UEdgww4CvufsMM0sBy939bTPbDNwE/DvZ8DJLFYENNpowsO0QYEn0bxfgQeC6pILKB0puoorAhmkHVDyRZSVhROVPgenAF5IKKkGqCGyYg4Fh7r7ezMqB7dx9cjSy+2/Jz7+djFByy3OqCGywGcBgYC4wBbjYzBYAPwQWJRdWMlQR2GAFwIbo/XJgF8CBhYQx76SFKLnJNlQRWKP7gN7R+18B/yKMX7YZODupoCTrTQP2B2YDk4CrzKwMuACYmWRguU7JTSqpIrB27j4m9n6ymfUntOTm6wHSqgiswy1Ah+j9dcA/CKOWFwOnJhVUPihIpVL1LyV5wcyeAXoQutq2qQh097wqmqiNmXUEcPd1SceSDeqrCHR3PWQ7xsy6A6vcXT++LUj95hI3DLjK3WcQWmzLo0q3qwgVgXnNzEaZ2XxgNbDazBaY2U/MLC8eJF2H0YSKwC6E60t7AQcRrk2enFhUWcDMHjOzTvFp7r4SaG9mGrW8BSm5SVxNFYGQvxWBlczsDuBG4GHg2Oj1EHA9cHtykWWFg4Gb3X09UFkRCFwJ3J1oZMk7m5pHQmgHnJXhWPKKrrlJnCoCa/cD4AfuHi+Bf93MnJDwrkwmrKygisBqoq7HgujVzcxKY7OLgK8BS5OILV8ouUmcKgLrNrWWafneA6KKwG0VE7r2U4Sej+pSwA0ZjSjPKLlJJVUE1ukpQgt2ZLXplwB/zHw4WUUVgds6mtBqe51w3TH+2LESYJ67f5ZEYPlC1ZJSI1UEVmVmDwKnAYuBidHkQ4GdgTFAZbeTu/844wFmGVUEBma2G+HkMK+PQxKU3KQKMxsFXE64bgLwGXAPMDqf/wc1s7ENXDTl7se0aDBZJqr6G+nua6tN7wD81t3PSyay7GBm+wEXAQOB89x9sZmdSGi9fZBocDlM3ZJSKaoIvBC4E5gQTT6cUBHYhzwumnD3o5OOIYudDfwcWFttekVFYN4mNzP7CvA88BJwDFsrJwcC5wAnJhJYHlBykzhVBNbDzHoSfpimuPvmpONJkioCG+Qm4HJ3f8DM4sl/HPDTZELKD0puUp0qAmsQ3Yj7GKE4IAXsCcw2s4eAJe5+Y4LhJUUVgfXbF/hnDdNXAt0zHEteyesfLNlGRUVgdaoIDDdq7wwcAGyMTX8R+HYiESXvaGAEoeV2CqHbreL1JaCfu9+SXHhZYSVbr1/HHUC4D1BaiFpuEtcWOM3MjqOGikAz+03FgnlYEfhN4NvuPiUarLTCJ8DuCcWUKHcfD2BmA1BFYG3+B7jTzL5LaMluZ2bDgLuAxxONLMcpuUncYGBy9H636N8l0Wuv2HL5+CPWDVhRw/ROhAcF5y13n2dm+5mZKgK3dR3wBDCP0MKdTugxG0O4P1BaiJKbVFJFYJ3eI7TeRkefKxL8RcA7SQSULVQRWDt33wKcbma/IHTVpoAJ7p6vT27JGCU32YYqAmt0DfCyme1D+P/m8uj9ocCRiUaWPFUE1qGme0fNLO/vHW1pKiiRSmbWycyeBZYRWiO7RNMfMrMbk4wtae7+DuGevzbALEIhxSLgsOgJ+PlMFYG10GgSyVFykzhVBNbCzPYGtrj72e6+LzCKcA3l62aWtyOUR1QRWLuKe0dvcffXo9cthIdKn59wbDlNyU3ivgmMcvcpVC0ayduKwJjHgKEAZrYr8DdCq+SHwM0JxpUNKioC+7JtReBTiUaWHXTvaAJ0cCVOFYG1i1eSngJMcvcTgDMJwwLls+uAOYSKwI6EisCxwFuoIlD3jiZEBSUSp4rA2hURhiqBcL2t4hrTLKBXIhFlCVUE1kn3jiZEyU3iVBFYu2nAJWb2IiG5XR1N34XwGKq8porAWune0YQouUkld3/HzA4HfsbWisD/ECoCP0o0uORdBfwduAJ4MnY8vkkYfTpvaTSJ2une0eRoPDepFFUElrm7R5+/Qhiy5GPgDnfP6+tuUVVkZ3dfFZvWH9jg7ssSCyxhZrYSuLDaaBKY2SnAw+7eI5nIJJ+p5SZxjxGut3msInA84YJ4Z7Z2xeWlKLmvqjZtbjLRZB1VBEpW0R+exKkiUBpDFYGSddRykzhVBEpjqCJQso6Sm8SpIlAaQxWBknWU3CROFYGSNlUESjZStaRUoYpAEckFSm4iIpJzVC0pIiI5R8lNRERyjpKbiIjkHCU3ERHJOUpuIiKSc/4fpMENJjHpVDkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "corr = iris[iris['species'] != 'virginica'].corr() \n",
    "mask = np.tri(*corr.shape).T \n",
    "g = sns.heatmap(corr.abs(), mask=mask, annot=True, cmap='viridis')\n",
    "g.set_xticklabels(g.get_yticklabels(), rotation = 90, fontsize = 14)\n",
    "g.set_yticklabels(g.get_yticklabels(), rotation = 0, fontsize = 14)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6144913",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\" style=\"background-color: white; border: 2px solid; padding: 10px\">\n",
    "    <b><i class=\"fa fa-info-circle\" aria-hidden=\"true\"></i>&nbsp; Notes</b><br>\n",
    "    <p style=\"color: black\">\n",
    "        (i) Correlated data has typically less power to restrict the model; correlated variables translate into wider combinations of coefficients that are able to explain the data.\n",
    "    </p>\n",
    "    <p  style=\"color: black\">\n",
    "       (ii) One solution when dealing with highly correlated variables is to remove one (or more than one) correlated variable.\n",
    "    </p>\n",
    "    <p  style=\"color: black\">\n",
    "        (iii) Another option is scaling all non-binary variables to have a mean of 0, and then using:\n",
    "    </p>    \n",
    "    <center>\n",
    "        $\\beta \\sim StudentT(0,\\nu,sd)$\n",
    "    </center>\n",
    "    <p  style=\"color: black\">\n",
    "    $sd$ should be chosen to weekly inform us about the expected value for the scale. The normality parameter $\\nu$ is typically chosen to be in the range (3,7). \n",
    "    <b> This prior is saying that in general we expect the coefficienct to be small, but we use wide tails because occasionally we will find some larger coefficients.\n",
    "    </b>\n",
    "    </p>\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91fa75de",
   "metadata": {},
   "source": [
    "$\\beta \\sim StudentT(0,\\nu,sd)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3907e55c",
   "metadata": {},
   "source": [
    "## Recipe 2: Dealing with unbalanced classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f65d95c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]\n"
     ]
    }
   ],
   "source": [
    "df = iris.query(\"species == ('setosa', 'versicolor')\") \n",
    "df = df[45:]   # let's select two unbalanced classes\n",
    "y_3 = pd.Categorical(df['species']).codes \n",
    "x_n = ['sepal_length', 'sepal_width'] \n",
    "x_3 = df[x_n].values\n",
    "\n",
    "print(y_3) #this is why is unbalanced"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c25dc3ba",
   "metadata": {},
   "source": [
    "Doing the usual thing: build the logistic regression..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "987cf733",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/deprecat/classic.py:215: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n",
      "  return wrapped_(*args_, **kwargs_)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [β, α]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\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:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='8000' class='' max='8000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [8000/8000 00:04&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 10 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_3: \n",
    "    α = pm.Normal('α', mu=0, sd=10) \n",
    "    β = pm.Normal('β', mu=0, sd=2, shape=len(x_n)) \n",
    "     \n",
    "    μ = α + pm.math.dot(x_3, β) \n",
    "    θ = 1 / (1 + pm.math.exp(-μ)) \n",
    "    bd = pm.Deterministic('bd', -α/β[1] - β[0]/β[1] * x_3[:,0]) \n",
    "     \n",
    "    yl = pm.Bernoulli('yl', p=θ, observed=y_3) \n",
    " \n",
    "    trace_3 = pm.sample(1000, target_accept=0.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "976f39e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/arviz/plots/hdiplot.py:157: FutureWarning: hdi currently interprets 2d data as (draw, shape) but this will change in a future release to (chain, draw) for coherence with other functions\n",
      "  hdi_data = hdi(y, hdi_prob=hdi_prob, circular=circular, multimodal=False, **hdi_kwargs)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'sepal_width')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABYP0lEQVR4nO3dd3xU15n4/8+dpi7UC2poBjjGxgaMwca4xnGcxOl27KyTTZx1kk3fkrLJpu5uNskmu2mb+Jdvks1umhPHcQGMMcXYxjbFGGNsDFxAEggJCQnUNZp67++PqxnUpQGNZqR53q+XXqBp99HV1TxzznnOOZppmgghhBBzjS3RAQghhBDxIAlOCCHEnCQJTgghxJwkCU4IIcScJAlOCCHEnCQJTgghxJzkSHQAk+ns7IzOY8jMzMTr9SYynFlBztPUybmaGjlPUyfnamqm6zzl5+dr4903q1pwdrs90SHMCnKepk7O1dTIeZo6OVdTMxPnaVYlOCGEEGKqEtJFqZQ6AdSMcdcTuq7fPsPhCCGEmIMSNQa3ChjaPi0H9gF/Tkw4Qggh5pqEJDhd19uHfq+Uug/oQRKcEEKIaZLwMTillAbcB/xe1/WBRMcjhBBibkiGaQK3ArXAL8e6MzMzM1ptY7fbycnJmcHQZic5T1Mn52pq5DxNnZyrqZmJ85QMCe6jwF5d1w+MdefQeRI5OTn09vbOVFyzlpynqZNzNTVynqZOztXUTNd5ys/PH/e+hHZRKqVKgHcyTutNCCGEuFCJHoO7F/ADf0xwHEIIIWaQYRjEe8PthCW4weKSjwB/0nW9L1FxCCGEmDl+v5/nnnuOL37xizQ0NMT1WIkcg7sJWAR8IIExCCGEmAGBQIC9e/eyfv366NhbIBCI6zETluB0XX8aGHeRTCGEELNfMBjk5ZdfZt26dXR1dVFUVER+fj5tbW1xP3YyVFEKIYSYY0KhEPv372fdunV0dHRQWFhIdXX1jMYgCU4IIcS0CYfDHDhwgHXr1tHe3k5BQcGMJ7YISXBCCCEummEYvPbaazz22GOcOXOG/Pz8hCW2CElwQgghLphhGBw+fJhHH32U5ubmpEhsEZLghBBCxMw0TY4cOcJjjz1GY2Mj8+bNo6ZmrF3QEkcSnBBCiCkzTZOjR4+yfv166uvryc3Npbq6Gk1LvqJ4SXBCCCEmZZomdXV1rF+/nuPHj5OVlZW0iS1CEpwQQohxmabJiRMn2LBhA4cPHyY7O5uqqqqkTmwRkuCEEEKMqbGxkccff5yDBw+SkZGR9C22kSTBCSGEGKapqYknnniC/fv3k5GRMWtabCNJghNCCAFAS0sLmzZt4qWXXiItLY2qqipstkRvOnPhJMEJIUSKO3PmDJs3b2bPnj04HA4qKytndWKLkAQnhBAp6uzZs2zZsoUXXngBh8PB/PnzsdvtiQ5r2kiCE0KIFNPR0cG2bdvYsWMHdrt9ziW2CElwQgiRIrq6uti+fTvbt29H0zTKyspwOOZuGpi7P5kQQggAenp6eOaZZ9i2bRumaVJaWorT6Ux0WHEnCU4IIeaovr4+duzYwZYtWwiFQimT2CIkwQkhxBzT39/P888/z5NPPkkwGKSkpASXy5XosGacJDghhJgjBgYG2LVrFxs3bsTv91NcXExaWlqiw0oYSXBCCDHL+f1+9uzZw4YNG/B6vRQXF1NcXJzosBJOEpwQQsxSgUCAvXv3sn79enp7eykuLqawsDDRYSUNSXBCCDHLBINBXn75ZdatW0dXVxdFRUXk5+cnOqykIwlOCCFmiVAoxIEDB3jsscc4d+4cBQUFVFdXJzqspCUJTgghklw4HOa1117j0Ucfpb29nfz8fElsUyAJTgghkpRhGBw6dIhHH32UlpYW8vLyJLHFQBKcEEIkGdM0OXLkCI899hiNjY2S2C6QJDghhEgSpmly7Ngx1q1bR319Pbm5ubNuF+1kIglOCCESzDRN6uvr2bBhA7quk52dLYltGkiCE0KIBDFNk5MnT7JhwwYOHTpEVlaWJLZpJAlOCCES4NSpUzz++OO89tprZGRkSGKLA0lwQggxg5qamti0aRP79+8nPT2dqqoqSWxxIglOCCFmwOnTp3nyySd56aWXSEtLo7KyEpvNluiw5jRJcEIIEUetra1s3ryZF198EafTKYltBkmCE0KIOGhra2PLli3s2rULp9NJRUWFJLYZlrAEp5QqB74LvBXIAeqBT+i6/myiYhJCiIvV1tbGI488wq5du7Db7cyfPx+73Z7osFJSQhKcUioPeAF4HrgdaAfcQFsi4hFCiIt19uxZtm7dyp49ezBNk/LycklsCZaoFtwXgRZd1z845LaGBMUihBAX7OzZs2zbto3nn38em81GdXU1oVAo0WEJEpfg3gU8qZR6ELgZOA38CviZrutmgmISQogpO3fuXDSxaZpGWVkZDocDu90uCS5JJCrBuYFPAj/EGodbDvz34H0/HfrAzMzMaDPfbreTk5Mzc1HOUnKepk7O1dTIeTrv7NmzbN68mWeeeSbaYnM4zr+VappGWlpaAiOcHTRNIysrK67XVaISnA14Sdf1Lw9+v18ptQj4FCMSnNfrjf4/JyeH3t7eGQtytpLzNHVyrqZGzhN0dHSwbds2duzYgc1mo6SkBIfDQTgcJhwORx+XlpaG3+9PYKSzg2ma9Pf3X/R1NdFO5olKcC3AoRG3HQb+LgGxCCHEuDo6Onjqqad49tlnsdls0a5IkfwS9Vt6AVAjblsMnExALEIIMcrQxDZ0jE3MHon6bf0Q2KmU+grwILAC+CzwzwmKRwghAElsc0lCfmu6ru9VSr0L+DbwNaBx8N/7ExGPEEKcO3eOp556KjrGJolt9kvYb0/X9Y3AxkQdXwghwKqKfOqpp0aV+4vZT36LQoiUNHKCdmlpqSS2OUZ+m0KIlNLW1sa2bdt44YUXsNvt0mKbw+S3KoRICWfOnGHr1q3s3r1bFkFOEZLghBBzWmtrK1u2bGHPnj04HA5ZBDmFSIITQsxJLS0tbN68mb179+J0OqXFloIkwQkh5pTTp0+zefNmXnrpJdloNMVJghNCzAnNzc3RxOZyuSSxCUlwQojZrbm5mSeffJJ9+/bhcrmorKyUxCYASXBCiFlKEpuYjCQ4IcSscvr0aTZt2iSJTUxKEpwQYlZoaWlh06ZN0TE2SWxiMpLghBBJrbW1lc2bN/Piiy/idDolsYkpkwQnhEhKbW1tbNmyhV27dkm5v7ggkuCEEEll5CLIMkFbXChJcEKIpNDZ2cm2bdvYsWMHmqbJklriokmCE0IkVE9PD9u3b2f79u2Ypinb1ohpI1eRECIh+vr62LFjB1u2bCEUClFaWorT6Ux0WGIOkQQnhJhRPp+PF154gY0bN+L3+yktLcXlciU6LDEHSYITQsyIYDDIiy++yLp16+jv76ekpIS0tLREhyXmMElwQoi4MgyDAwcO8Mgjj3Du3DmKioooKChIdFgiBUiCE0LEhWmaHDt2jIcffpjGxkYKCgqorq5OdFgihVxQglNK2YFRfQu6rnsvOiIhxKzX3NzMY489xsGDB8nJyaG6uhpN0xIdlkgxU05wSqlc4NvAe4ASYKyrVSatCJHCurq62LRpE88//zxpaWmS2ERCxdKC+3/A24BfAYeAQFwiEkLMOj6fj2effZaNGzdimqasPiKSQiwJ7jbgH3Rd/1W8ghFCzC7hcJj9+/fz0EMP0dvbS1lZmZT8i6QRS4LrB5riFYgQYvYwTZO6ujoeeughGhsbKSoqkgISkXRiSXD/BXxSKbVF13UjXgEJIZJbW1sb69at4+WXX5YCEpHUJkxwSqnvjbhpGaArpZ4GukbcZ+q6/k/TGJsQIon09fWxdetWnnrqKRwOB1VVVbJ9jUhqk7Xg3jvie2PwObeO8VgTkAQnxBwTCoV45plneOCBB/D7/ZSVlcliyGJWmPAq1XW9dqYCEUIkF9M0OXz4MH/+85/p6OggPz+f4uLiRIclxJTFMg/ug8BGXdfPjXFfAfA2Xdd/O53BCSESo6mpiUceeYTDhw+Tl5dHbW0tfr8/0WEJEZNY+hn+F1gDjEpwQO3g/ZLghJjFOjs7eeKJJ9i5cyfp6elSQCKmVTgcpqmpibq6Onp7e/nYxz4W1+PFkuAmusoLgZ6LjEUIkSADAwM8++yzPPHEEwAyUVtMm/7+fhoaGqirq6OhoYFAIIDdbqempgbDiG9B/mRVlO8E3jnkpq8ppdpHPCwduB7YO82xCSHiLBwOs2/fPh5++GH6+vpkbzZx0UzT5MyZM9TX11NXV0dLSwsA2dnZXHLJJXg8Hmpqaujq6op7Fe5kLbgS4PIh33uAshGPCQBbgG9N9aBKqW8C3xhx8xld10e+thAiDkzTRNd1/vKXv9Dc3ExRURFVVVWJDkvMUoFAgJMnT1JXV0d9fT19fX0AlJeXc9111+HxeCgpKZnx7u7Jqih/CfwSYHDu2yd1XT88TcfWgZuGfB+eptcVQkygqamJxx57jNdff53c3FxqamoSHZKYhbq6uqirq6Ouro5Tp04RDodxuVzU1tbi8Xiora0lKysroTFOeQxO1/Wbp/nYIV3XW6f5NYUQ4zh37hxPPPEEu3btIiMjQwpIREzC4TDNzc3RVtq5c1a9YUFBAVdeeSVut5vKysqkGrudbAzu67G8mK7r/xrDw91KqdOAH9gD/LOu6/WxHE8IMbne3l62b9/O1q1bsdvtVFZWygokYkq8Xm90LO3EiRP4/X7sdjtVVVUsW7YMj8dDfn5+osMc12QtuM+M+D4DyBz8fx+QPfh/7+DXVBPcHuBe4AjWON9XgZ1KqcvGmmcnhIidz+fj+eefZ+PGjYRCIUpLS3E6nYkOSyQx0zRpa2uLJrXTp08DkJWVxeLFi6MFImlpo/a7TkqTjcFFly1QSq0B/oCVjB7VdX1AKZWBtQHqvwHvn+pBdV3fNPR7pdRuoB74EPCDofdlZmZGm7x2u52cnJypHiZlyXmaurl4roLBILt37+ahhx6ir6+PsrIy0tPTL+o1NU2bNW9qiTbbzlUgEKChoYHjx49z7Ngxent7AWuqyI033siiRYsoKyub9u5sTdPIysqK699fLPPgfgJ8W9f1ByI36Lo+APxBKZUF/Ay48kKC0HW9Tyn1OrBo5H1erzf6/5ycnOjJF+OT8zR1c+lchcNhXnnlFR599FE6OjooKiqivLwc4KJXIUlLS5OVTKZoNpyrrq6uaCutsbExWiCyYMECrrvuOmpra8nOzo4+PhCY/v2tTdOkv7//ov/+JuoijSXBLQVOj3NfM7AkhtcaRimVDlwCPH2hryFEqjIMg0OHDvHII4/Q0tJCQUGB7M0mhjEMI1ogUldXFy0Qyc/PZ/ny5SxcuDDpCkSmQywJ7ijwj0qpp3Rdj348GUxO/4hV9j8lSqn/BDYAjVhjcF8DsoDfxBCPECktMpdt3bp1nDhxgry8PCn5F1Fer3fYCiJ+vx+bzUZVVRVXXHEFHo+HgoKCRIcZV7EkuM8ATwBNSqmtQBtWcroVq/DkLTG8ViXwR6AIaAd2A9foun4yhtcQIiWZpklDQwPr1q3j6NGjsumoAKzror29PVrGf/r0aUzTjBaIuN1uFixYMKvGBy9WLPPgdiilFgH/AKwCVgCtWIss/0jX9fG6L8d6rffFGqgQAk6dOsX69et5/fXXyczMlMSW4oLB4LAVRCLjWaWlpaxZswaPxxOXApHZIqZdC3VdbwG+GKdYhBDjaGlp4YknnmDfvn2kp6dTVVWVsm9aqa67u3tYgUgoFMLpdLJgwQLWrl2L2+0eViCSymRbXiGSWHt7O5s3b2bXrl04nU6ZpJ2CIgUikaR29uxZAPLy8qKTrSsrK2WX9TFMtpLJi8C9uq4fUkrtBcyJHq/r+urpDE6IVNXZ2cnWrVvZsWMHdrtdtq9JMQMDA8MKRHw+HzabjcrKSm6++WbcbjcFBQXSip/EZCn/dWBgyP8nTHBCiIvT09PD9u3b2b59O6ZpUlZWJp/MU4Bpmpw9ezY6ltbc3IxpmmRmZrJw4ULcbje1tbUpVSAyHSZbyeTDQ/5/b9yjESJF9ff3s2PHDjZv3izLaqWIYDBIY2NjNKn19Fh7RpeWlnLNNdfg8XgoLy+XVtpFmPJHQ6XULcAuXde9kz5YCDElAwMD7Ny5k40bN+L3+2XD0Tmup6cnOtl6aIFITU0Na9aswe12z7ml4xIplr6PLUBYKbUfeG7w63lZHFmI2Pn9fl588UXWr19Pf38/paWl0v00BxmGwenTp6MFIu3t7QDMmzcvOtm6qqpKuqHjJJazWgLcAFwH3Aj8HWBTSh1hMOHpuv6H6Q9RiLkjGAyyb98+1q1bR3d3N8XFxRQWFiY6LDGNBgYGOHLkCPX19TQ0NDAwMICmaVRWVnLTTTdFVxCRrsf4i2Wi9zng0cEvlFLZwM3A54CPAR/F2m1ACDFCZCHkxx57jHPnzlFYWCjrRc4Rpmly7ty5aNdjpEAkIyMDt9uNx+NhwYIFF72jg4hdTO3iwaR2LXD94NdqwAdsxGrFCSGGMAyDgwcP8uijj3LmzBny8/Mlsc0BoVAoWiBSV1cXLRApKSlh7dq1LFiwgLKyMpmzmGCxFJm8BCwDzmAls4eAzwKv6bou0weEGMI0TQ4fPsxjjz1GU1MTeXl5kthmud7e3mEFIsFgMFogEql6zMnJmRXb5aSKWFpwy4AgsAvYCbyg6/qrcYlKiFnKNE2OHj3KunXraGhoIDc3V5bVmqUMw6ClpSVaINLW1gZYBSJLly7F4/FQXV0tBSJJLJbfzDzOd0++B/gPpZQfeAHYAezQdX339IcoRPIzTZPjx4+zYcMGjh07Jiv8z1I+n48TJ05E56ZFCkQqKiq48cYb8Xg8FBYWyu91loilyMQLbBv8QinlBG4BvgR8F2uVE1lLSKQU0zSpr69nw4YNHD16lKysLElss4hpmnR0dES7HpuamjBNk/T09GiBSG1trRSIzFKxFpkUc77A5Hqsbksb1jJeUmQiUkZkT7bHH3+cI0eOkJmZKV2Rs0QoFOLUqVPRpNbd3Q1AcXExV199dXQFESkQmf1iKTLRgYVAGHgZeBr4V6z5b53xCU+I5CKJbXbq7e2NjqWdPHmSYDCIw+GgpqaG1atX4/F4yM3NTXSYYprF0oL7I1YrTZbrEikn0hW5ceNGjhw5QkZGhiS2JGaaJi0tLdFWWqRAJDc3l6VLl+J2u6murpb1Pue4WMbgvjnVxyqlbMBx4O26rr9+AXEJkRQixSMbN27k6NGjktiSmN/vp6Ghgfr6eurr6/F6vcMKRNxuN0VFRfK7SyHxqm/VgAWALK4nZqVIuf/jjz9OXV2ddEUmoaEFIvX19TQ1NWEYBunp6dTW1kYLRDIyMhIdqkgQmcAhxBCGYXDkyBE2bNjAiRMnyM7OlsSWREKhEE1NTdGux66uLsAqEFm1ahUej4f58+dLgYgAJMEJAVhrRR46dIj169fT1NQk89iSSF9fX7RA5MSJE9ECkerqalatWoXb7WbevHmJDlMkIUlwIqWFw2H279/P+vXraW1tZd68eZLYEsw0TVpbW6OttDNnzgCQk5PDZZddhtvtpqamRgpExKQkwYmUFAwGeeWVV3jyySdpaWkhLy+PmpqaRIeVsvx+PydOnIgWiPT396NpGvPnz+eGG27A7XZTXFwsHzxETCTBiZQSCATYt28fGzZsoKuri9LSUlkEOUGGFoicOnUKwzBIS0uLFoi43W4pEBEXJV4JzgROArKktkgKPp+PvXv3snHjRnp6eigoKKC6ulpWfp9B4XCYpqYmjh8/Tn19PZ2d1voQRUVFXHXVVXg8HioqKqRAREybuCQ4XdcNoDYery1ELLxeL7t372bTpk309/dTVFQkLbYZ1NfXR0NDQ7RAJBAIYLfbqa6uZuXKlXg8HikQEXEzYYJTSu3Fao1Nia7rqy86IiGmQU9PDzt37mTLli34/X6KioooLCxMdFhzXqRAJFL12NraCkB2djZLliyJbjHjcrkSHKlIBZO14F4nhgQnRKJ1dHSwY8cOnn76acLhMMXFxaSlyXoD8eT3+zl58mR0PK2/vx+A+fPnc/311+N2uykpKZECETHjJkxwuq7fO0NxCHFRWltb2b59Ozt37kTTNEpKSqSMPI46OzujZfxDC0QWLFjAwoULqa2tJTMzM9FhihQnVZRi1jJNkxMnTrB161YOHDiAw+GgvLwcu122JZxu4XA42kqrq6uLFogUFhZGx9IqKirk3IukEut+cAuADwCLgVE7AOq6ftf0hCXE+MLhMEeOHGHTpk3U19eTlpYm1Xdx0N/fH52XduLECfx+P3a7naqqKq688ko8Hg95eXmJDlOIccWyH9xKYAfQiJXgXgXmYS2q3IS1e4AQcePz+di/fz9PPPEEHR0dZGVlyTqR08g0Tc6cORMtEGlpaQEgKyuLJUuWUFtbS01NjRSIiFkjlhbc94GHgPuAIHCfrusvK6Wuxdor7ntxiE8IOjo62LVrF9u3b2dgYICCggKqqqoSHdacEAgEhhWI9PX1AVBeXs51112H2+2mtLSU9PR0mS8oZp1YEtxy4D8AY/D7dABd13cqpf4F+C7w5LRGJ1JWZHztmWeeYd++fWiaRlFREcXFxYkObdbr6uqKJrTGxkbC4TAul2vYFjNZWVmJDlOIixZLgjOBgK7rplKqDagBdg7edwpYNN3BidTj9/t57bXX2LJlC01NTaSlpUnhyEUKh8M0NzdHk9q5c+cAKCgoYMWKFXg8HiorK+UcizknlgR3CPAATwO7gH9QSr0EBIAvAnUXEoBS6svAt4Gf6br+6Qt5DTH7nTlzhj179vDss8/i8/nIzc2V8bWL4PV6aWho4Pjx49ECEZvNRnV1NcuWLcPj8ZCfn5/oMIWIq1gS3C+wWm0A/wxsAY4Mft8P3BnrwZVS1wAfwypYESkmEAhw5MgRtm/fzrFjx7DZbNINeYFM06StrS1aIHL69GnAKhBZvHgxHo+HmpoamfQuUsqUE5yu678b8v/DSqklwLVYY3G7dV1vi+XASql5wB+AvwG+EctzxexlmianT59m7969PPfcc/h8PrKysqisrJTWWowiBSKRpBYpECkrK2Pt2rV4PB5KS0vlvIqUdcETvXVd78NqxV2oXwB/0XX9aaWUJLg5rre3l4MHD/L000/T3NyM3W6nsLBQWmsx6urqiia0SIGI0+kcViCSnZ2d6DCFSAqxTvQuAf4eWA2UAy3AHuAnuq6fieF1PgosxJo0PqHMzMzo4LfdbicnJyeWkFNSspynUCiErus899xz7Nu3D8MwyMvLw+PxJE2rQtO0pO62MwyDU6dOcezYMY4dO8bZs2cBq0DkqquuYtGiRVRXV8e9QCTZz1MykXM1NZqmkZWVFdf3qlgmeq8FngBCwFasopMS4OPAZ5RSb9F1/YUpvI7CKiq5Ttf14GSP93q90f/n5OTQ29s71ZBTViLPk2matLS0sG/fPp577jn6+/tJS0ujsLAw+iYcCAQSEttYknE/uEiBSF1dHQ0NDdECkcrKSm6++WY8Hg8FBQXRx4dCIUKhUFxjSsbzlKzkXE2NaZr09/df9HvVRMVSsbTgfgrsA96u63p/5EalVDbwOPDfwJVTeJ01QBHwupXrALADNyilPg5k6bouV8cs09/fH+2CbGxsxG63U1RUNOyNWIzNNE3a29ujZfynT5/GNE0yMzNZtGgRHo+HBQsWSKtAiBjFkuAuAe4cmtzAGotTSv0n1ionU/EY8NKI2/4XOIbVskuej/diQoZh0NDQwM6dO9m7dy/hcFjK+6coGAwOW0Ek8im2tLSUNWvW4PF4KCsrk/MoxEWIdR5c2Tj3lXN+ysCEdF3vArqG3qaU6gc6dF0/GEM8IkG6u7vZv38/Tz31FB0dHbhcLkpKSnA4ZHOKiXR3dw8rEAmFQjidThYsWMDatWupra1NirFTIeaKWN6RPgP8TinVBzym67pfKZUGvBv4EvDBeAQoZsbuPSZ/fNDk3Dm45mq4530aBQXnWw/hcJijR4/ywgsv8MorrwBW37esCTk+wzA4ffp0dIuZSIFIXl4ey5Ytw+12U1VVJR8MhIiTWP6y1gGZwAMAg4kuUo/sAx4dMqaGruslU31hXddviiEOMc3+/JDBL/4HfD7r++Zm2LLV5De/hlDoLPv27eOZZ56ht7eX9PR0WTprAgMDA8MKRHw+HzabjYqKCm666aZogYh0PQoRf7EkuJ9hrUcp5pCBAZNf/Ap8Q8p6/AEffa2H+fgnnqWosA5N0ygsLJS9v8ZgmiZnz56NjqU1NzdHC0Q8Hk90bpoUiAgx82JZyeSbcYxDJEh9A9gdYPoMgoETePtfxOvdg2kGCQezWb5MVhgZKRgM0tjYGB1P6+npAaCkpIRrrrkmWiAiG7AKkVgxd/4rpfKBpUAVsEnX9U6lVDrWTgPGxM8WycYIt3Gu/RW6unYQDnWi2ZzY7YVompP8POKe3EwT2s+CaUBxMSRrTujp6Ym20k6ePBktEKmpqWHNmjW43W4pEBEiycQy0duBVcb/KSADq7tyFdAJPIxV+i9Lbs0CPT09HDp0iB07dnDy5ElMA2y2Amyu8wUjDgcsXRrf5Hb2LGx9yiQUAg1Ag5tu0KisjOthp8QwDFpaWqIFIu3t7QDMmzePK664ArfbTXV1tRSICJHEYvnr/Hfgo8CnsbbMqR9y3zqsFU0kwSUpr9eLruvs3LmTw4cPY5omubm5VFZW8o4Sje1Pm7S1n29BrV4F5eNNCpkGoRA8ucVk5KIm258xuePdGonYb9Pn8w0rEBkYGEDTNCorK7nxxhtZuHChFIgIMYvEkuA+CHxJ1/X/VUqNLKGrA9zTF5aYDj6fj+PHj7Nnzx4OHDhAOBwmMzOT+fPnDxsfSkuDt7xZo7/fqqTMy4N4F0k2nrK6J0cyTTheZ7LsivgnkcgKIocPH6auri5aIJKRkYHb7cbtdlNbW0t6enrcYxFCTL9YElwe429q6sJabkskmM/no76+nh07drB//37C4TBpaWmUlpZOWtqflcWMtZz8fjDGGLENh89PV4iHUCg0rECku7sbsApErr76ajweD+Xl5VIgIsQcEEuCOwi8E9g2xn1vAV6elohEzLxeL8ePH2fv3r28+uqraJqGzWajuLg4aceIysthrJ4+pxMqKqa39dbb2zusQCQYDOJwOKipqWHt2rVUV1eTm5s7rccUQpxnmibhcJhwOEwwGCQcDs/IgtSxvPt9C3hYKZWBte6kCSxXSr0b+FvgHXGIT4yjq6uLo0ePsnfvXo4cOYJhGKSnp1NcXExWVlbSr2aeNw88bqirt8bjwCpsKS6GivkX99qGYdDa2hotEGlrs/bizc3NZenSpXg8HqqqqnA6nbLyuxBTZBgG4XA4untFKBQa9r1hGNhstugYtaZpmKYZ/XK5XGRlZUW/li9fHvf9IDVzrIGQcSil7gK+B1QPubkZ+Jyu63+e5tgA6OzsjAaYytvlRJZ90nWdvXv3curUKcDaLy8vL29Y9+NsetM+cRL0oyZGGBYu1PC4L2yqgM/n48SJE9GWWqRApKKiIjrhurCwcFSByGw6V4kk52nqkvVcRVpRYyWnyPeapkV7gCLPiXzZbDbS09OHJans7Gyys7PJzc0lKyuL9PR00tLShv0b+Ro5RDJd7+f5+fnjdvnElOAilFKLsba86QJODy6gHBepnOB6eno4ceIEBw8e5MCBA/T3Wxs55ObmkpOTM241X7L+gU0n0zTp6OiIttKam5ujrVi32x1dQWSyApFUOFfTQc7T1MXrXJmmOSoxjWxBRRJU5L0h8v4eSVJpaWlkZWWRmZk5LEHl5OSQnZ0dTUZDk1RGRgZpaWm4XK5prSCeiQQXyzy4TwA5uq5/T9f1o0qpTGALUK6UegV4p67rTRcdbQrzer00NjZy7NgxXnnlFVpbW9E0DYfDQX5+/oQb+6WCUCjEqVOnogUiXV1dABQXF7N69WrcbveoClEhEs00zWj3XiQhjfV/YFhyGtnFp2kaGRkZZGRkkJOTM2aSGtl6GvpvWlpayv1txLqbwE+GfP8T4DTweeCfgO8CH5i+0Oa+vr4+Tp06RV1dHQcPHqSpyfp8oGka8+bNo7JSlsnq7e2lvr6e+vp6Tpw4ES0Qqa6uZtWqVXg8HikQEXExtMU0MhkN/YokjUhycjgcBIPBYa0np9NJenp6NCllZWWRkZERbU1lZ2dHW0rjfU13CyoVxJLgqgEdQClVDKwFbtF1/RmlVABrx28xjsiivM3NzRw7doxDhw7R3t4evWBzc3OpqKhI+QvYNE1aW1s5fvw49fX1nDlzBrDOz2WXXYbH46G6uhqn05ngSEWyiySosb7C4XC0VTSyxQREu/vS09OjrabMzMxhX5HkFElAkaKlgoICgsEgLpcrmphk943EiCXB+bHmuwHcDHiB5wa/78CaJycG9fb20tLSQlNTE7quU1dXh8/nwzRNHA6HJLQh/H7/sAIRr9cbLRC54YYb8Hg8FBUVyblKYZECiWAwOOoLiBZGjExQGRkZZGdnk5eXF+3KGzrm5HK5hnXhpaenD0tMF3LNpVqtQDKLJcG9CHxKKdUEfBZ4Utf18OB9bqzuypRjmibd3d2cOXOG1tZWjh07Rn19Pd3d3dE/jszMTHJzcykqKkpwtMljaIFIU1NTtECktrY2WiCSkZGR6DBFHA1NWoFAgEAgEJ0jNbJQIlIgEfk7ysvLY968eeTn50e79zIyMqLdgJFiiVQbcxLDxZLgPgdsAF4DTgF/M+S+u4EXpjGupNTf309HRwft7e00NzdTX19PY2Mjfr8fTdMwDCPadSHjZ8OFw+HoeOPQApGioiJWrVqF2+2moqJC3pCSjNbfjv30K5iBXuy5FYTLrgDn5B88TNMkGAzi9/vx+/0EAoFRicswDFwuF3l5eZSVlZGfn09hYSH5+fnDKv0i/5duaRGrWPaDOwR4lFKFQIeu60PnF3weaJ3u4BIhFArR1dVFZ2cnHR0dNDc309TURHNzM319fdFuEJvNFp2DJn94Y+vr64tWPEYKROx2OzU1NdGkNm/evESHKcZh62jAcWwLmGEwTex9bdjOvE7oirsxnJkEAgH8fj8+ny/6IW9o8srOzqaoqIiioiKKi4spKiqKdg1GvmQjWBFPMa/jpOv6uTFue216wom/QCBAb28vvb299PT00NnZGe1ebGtro6ura9gfqcPhiA4w5+XlSatsApECkchYWmur9ZknJyeHSy+9FI/HQ01NjXwgmBVMtLrteP1BfCGTgZBJ2Aij2UIYB7YTKltOTk4OxcXFlJWVMX/+fPLz85k3bx65ubnk5ubK71kkXHIuVBgHra2t/PznP6e9vX3ULP3IoHJ6eroUfsTI7/dz8uTJaFKLTEavqKjg+uuvx+PxUFxcLOc0SZmmic/nY2BgAJ/PFy3OwN+Po9NHQTrU5tmpzLUzP8dGXoaNvCIbrk/9QFpfIumlTILr7u7m7Nmzc35sbGAAenoMMjKshYvjoaOjI9r1eOrUKQzDIC0tbViBSGZmZnwOLmJmmiZ+vz+axEKh0LAPeQUFBdEx0MhYWH6GRsWfX8JhBke9Xri4iIFkTm6midZRj2aGMQoXgibjuqkqZRIcMGwh0LkmFIJnd5g0NYPdHsAwYNkVTMu+auFwmKampmiBSGdnJwCFhYVcddVVeDweKRBJMMMw8Pl80a+hE5BN0yQ/P5/q6mrmz58f7U4sKCggLy9v/B0naq7BbNyNZpxPcqYjg+DKD8/Ej3RBbO066es/g9Z/FjQw03Lwve1HGPNXJDo0kQApleDmsud3WsktHLa+AA68Cjk54K6N/fX6+/uHFYgEAgHsdjvV1dWsXLkSt9tNXl7etP4MYmKRlpjX6x3WnRiZsFxUVITH42H+/PmUlpZGl3ebMIlNwPeW75G+7lPY2w6h2Z2YoQDBFe8ntOTtcfjppkFwgIw/fxD8PUQ+1mnBATIe/gj9H9kGGam91F0qkgQ3BwSDcPLk+cQWEQrBq6+auGsnb8WZpsmZM2eirbRIgUh2djZLlizB7XZTU1ODy+Wa5JXExTIMg4GBgWgii7TEDMOIdidWVlZSWlpKQUEBBQUF5ObmTv9qGRl5+N73B7TOBrLCvfRnVSV1knAcfwqMEKOudjOM8/DjBK/860SEJRJIEtwcEAiMvXkowMAEu2OPVyAyf/58rr/+etxuNyUlJXO2WzfRTNOMJrKBgYFhXbylpaUsXLiQmpoaSkpKooksER8wzPxatJwcSPLVOTTvWTBGjxlqIT9a35kERCQSbU4nOJ/PZOMmk+eeBwyTzk6T+Re5mWYyysy0NguNbBwaoWlQXjb8ts7OzmhCO3XqFOFwGJfLFS0QcbvdUiASB4FAAK/XS39/f3RjSMMwKC4u5tJLL2XBggWUlZVRVFREQUHBtO7ErnU34Xz5t9jOHsUoW0Zwxfsxs0um7fWTRbhiJWgOYHiSM52ZhKtWT/4CRhjH0SdxHFoPdifBy+8kXHvj+J8exQWxtR/Fsf+3mP1ncFZeTfDyuyA9Pgumz9kENzBg8tGPm7S0gt8PAT90nDOxO0wuUXPrgtU0WHMNPPf8+SRns1lJb9kVBidPNkWTWkdHBwAFBQVceeWV0QIRWQx2eozVvWiaJunp6VRXV0e7FyMTn+Ndam9reZWMv3wYwgE0I4TZvB/nq3/Ee8+fMfMXxPXYM80ou5xwzRrsJ3ehhQYAMB3phEuWEF5w3cRPNg3S130Se9NetKD1XPvJnQSXvofAG74a79BThr1uO+kbPwfhAJgGrsa9OPf/Hu8HHoHMgmk/3pxNcBs2nk9uAKZpjVG9uBc87viV0CdK7QKNzAx49aBJV5cXh70eI1zPb37TEC0QqaqqYsWKFbjd7pTfW246BINB+vv7h7XKTNOkpKSE5cuXR7t4S0tLJ9ygNp7Stn0dLeiNfq8ZAUx/kLRnv4fvXffPeDzx5nv7j3G89jDOg38BI0zo0ncSXPZXk04VsJ/chb3ppWhyA9BCAzgP/sVq8eZfQKWWGM4Ik77lq2ih8+MmWtgP3g5ce39J4MZ/mvZDztkE9+yO88ltKJsGZ89CefnMxxQvpmnS1tZmFYecrqOlpQWArKwsLrnkkugKIlIgcmEMw8Dr9eL1evH7/RO2yoqLi5PnPIf82M4eG3Wzhom9cXcCApoBNgehZXcTWnZ3TE+zN+yAIR8Eokwr+YUkwV00reskhEa/KWtGEMfx7ZLgYjFeBbthQrK8/1yMQCAwrECkr68PgPLycm688UZqamooLS2VApEYmKZJIBCgv7+fgYGBYfeVlZWxZMkSamtrKS4upqSkhNzc3OQ+vzY72BxWd9AIpkvGWYcyM/LA5hxdpGKzQ1qSb6g70IW9eR/hwoWQXzNzx/X1YOs7gzGvckoLcOPKBiM05l2mjMHF5o53a+x50cQ3pIpQwyrIKCxMWFgXpaurK5rQGhsbhxWIuN1u3G43WVlZpKWl4R+r+SqiTNPE6/XS1dVFb29vtFWWlZUVPZ8VFRXRxYKns+hjxtgchC65HceRjWhDkpzpSLe67URU6NJ34HrxF6MTnKYRWviGxAQ1BemPfgJ7wzPR782c+Xg/8DBk5MXvoOEArqf+Defh9daHAjNMYNVHCF7zyQkLcszsEsJll2M/fQDNPJ/oTEcGwSs/GJdQZ+Ff7dRcuULjb+41+dWvz4+3ZWbCbbcm8SfuEcLhMM3NzdGkdu6ctc51fn4+K1aswOPxUFlZKQUik4h0Mfb19REMBqPJrLi4mGXLlkUnRpeUlJCdnZ3crbIY+d/wVbTeVuynXwabC8IBQu6bCa7+WKJDSypmboU1sf3JLw2O15mYNqc1TulMztau85nvYm94Zvi8v97TZDxwFwP3bYnbcV07vo/zyOPWh6bBD06uvf+DmV1C6PL3Tvhc/9t+RPojH8XWdRLN5sAM+Qkuex+hS94Wl1i1yA64yaqzszMa4IXslNvTa/L6Iejs0Nn0xM+oqKiY9hink9frpaGhgbq6OhoaGqJjPlVVVXg8Hjwez6QFIqncgovMLevt7SUYDEZX+igvL2fhwoW43W7KysooKSkhLS0tZXZf1joasHWdxChchDkv9r+BVDlPBH3YT+/DtDkxKq60unhjNFPnKuvHy60ijRFMoP9jOyC7ePoPGg6S9bPVwwpFIox51Xjv2zz5a5gmtvYjZBq99Oe4MbMubiPo/Pz8cT+RztkWXERujsaaq0HXtaT8ZG6aJu3t7dEVRE6ftjZGz8rKYvHixdECEVm5fbTI0lV9fX3RidKRltmqVatYuHAh5eXl0WSWysyCWsIFUigxKWc64Zq1iY5iasYYW42w9Z7GiEeCC/nGHUfTBkbtpDY2TcMoWYKWk4MZ5w8CCUlwSqlPAX8LLBi86XXgW7qub0xEPNPN74cjuklbm1XssmSJRnbW+fuDweCwApHIp72ysjLWrl0bbWUkY0JOpGAwSG9vL/39/dg0DbP7NAWhMywvz8N989spv/JNlJWXk5ExhQHvCzXQRdrT38He+AJm+jwCaz5NWL1las/19+E8+Bfsp/ZgzKsiuPyeGZmLpvW24nzlAWxnjxIuX0bwirvjMudoJFv7ERyv/BFbfzsh983WGpbO9Lgfd9YJDuA4tA5Hw7MY2aWElt+DUbR40qeZ2SXQd2b00mRoGMVL4hEpuLIxs4rQeofvb22Ctdt7kklUC64J+CfgGGADPgQ8ppRaqev6qwmKaVr09cH6x01CQQiFofm0lezWXttDV6eV0E6ePEk4HMbpdFJbW8t1111HbW0t2dnZiQ4/aRiGQX9/P729vRiGAYDL5cLj8aCUwlP/R6q765ln92PSAS2/IHjyHAH3V+IXlPccWb98gzVpGjC950jf+I8ET71I4I3fmOS5HWT+4Q60gS60kA9Tc+B87SF87/xpXFsMtjOvk/HnD4ERQAsHsTfuxvXyb62J3nlVcTuu4/AG0rZ+3TpXpoG9cTfO/b9j4J4/Je2YVkIE+sl44C5sPS1ooQFMzYbz0Dp8t32HsHrzhE/1X/c50p/8IiZEk5wJhN03gSNOpeKahv/mr5L+xOch5LP+DjQbONIJ3PCF+BzzIiQkwem6vm7ETV9RSn0CWAPM6gS3d5+J32+9QYeCpwkE6gj46/nLQ2cByMvLY/ny5Xg8HqqqqqRAhPObbvb29uL3+9E0qzu5srKSVatW4Xa7KS8vp7CwEJvNhq15HxmH96PZrfEHDSA4gPO1hwhe8T7MQk9c4kzf8tVocoseF3C++icC1/3DhMsNufb8HK3/LNpg945mhiAUIm3zV/B+9Om4LQeVtvXraMH+6Pda2I9pBEnb8T187/jvuByToI+0bd8cPqE3NICt+xTOA38meNW98TnuLOTc/3ts3c3RsTTNNCDkI33r1+lf+Aawj5+o7M0vATY0jOhtGmA7e9Ra2SJO11R44S0M3PErXLt/jq3rJOHSpQTWfBKzcGFcjncxEj4Gp5SyA+8FsoGdCQ7nogwMDHDsWD0D/fUEAg2Ypg+w4XRWkpF5E3fesZCSkvh3DSW7yLiZ1+uNjpsVFBSwYsUKFi1aREVFBaWlpeNOmHY0PAvBMVaRNk0cJ58nGKcEZz/14hjdQYMxHdlIaPn4pfeO409Fk9tQmq8bracZc17lNEU5RMiPrf3I6GOaBvYTL0z/8QbZ2g6NuXKIFvLhOPakJLghHMe2jFkoAga2dh2j7PLxn1v/zLDkFqH1n0Xra8PMKZ3GSEdEV7ES3x2/jNvrT5eEJTil1OXALiAd6APerev6ayMfl5mZGW3l2O12cnJyLuh4WVlZOByOaS02iKwgcuzYMY4fP05TU9Pg3lyZuNIW4nK5cboWYLOlYbfB/Plp2O3xH1fTNC1piirC4TC9vb309vYSqdjNzMxk6dKlLFmyhMrKSioqKmLqnjVzisHuHDXIrtkcpOUWkR7DNRLLNWU60sZc7UID0gsrrRX3x3tueg6MMZ6umQZZ+SVoWRd2XU/ENDLGneituTKn/nP7e7F1N5GdXYpmn3yNOzO/GMzRb7wA9qyCqR+3/xwEB2BexYyNR5umCT0tYHegXeCC1DFdU5l5Y96umQaZ+aUTX1Np2dDfNvq5GGQVFKOlT/81NZ0u5v18qhLZgtOB5cA84E7gN0qpm3RdPzj0QV7v+TeUiym/7e/vJxQKXXT5fDAYpLGxMVr1GImntLSUNWvWEA67OdFYRjh8/g/SbofqKgiFAqNW/I+HRE4T8Pv99PT0RBcattlsLFiwgGuuuYaamhrmz5/PvHnzhr1hmaYZ0+9Vc7+RzGf/a1RrysSkv/L6mLZ1ieWaci77K1y77x92XBPA5qC/Yu2Ex3Usu4e0p/5t+ARXNMJll+MzXHHbiiZt8W04jm4eMdE7jcDS9xKc7JghP2nbvoFD30TY5sC02fHd8IVJ5zqRUUlmZiFat3f4ubK78C29i/Akx9X62kh74nPYWw6AZsfMyMf/5u9ObUeAi2BrPUj6E58fLKAwCRctxnf7D2Ieq4zlmrJf/j7ST700/PcDGDnlDKSVTHhdOJf9Fa7n/mtYV7BpcxCuXI0vaINgck/rmK7pFBNNm0pYgtN1PQAcH/x2n1JqFfAPwH2Jimk8PT090YTW2NhIKBTC6XSyYMEC1q5dS21tbfSTiGFA+FmTU01WYjMMyM+HtdfOzYpIv99Pd3d3dOxsZOusrKxs2lcBMXPK8b15cFKu7fwYpu8dP4nbthsAYffNsPv/YzCtnb+9wGNt3zABM6t41PMAzMz4Lqvjv+UbaN3N2NsOWa05I0i4Zi3BNZ+Y9LlpW79+PjkOjj2mPf1tzKwSwu4bx3+ipmFk5GPvPjX8diOMmTnJnCfTJP2he7F1NaKZYSCI1ttC+qMfx/uh9fHpygUY6CTjL/eiBc6PV9raDpHx4AfwfmSb1WMQB2ZmIRjh0benTz6UEVz2V9jOvI5D32TFZ4Yx8qrxveV78Qh1Vkr4GNwQNiAp+tUMw+D06dPRpHb27PkCkWXLluF2u6mqqhrzjdtmgzfcrNHTA52dkJMDBXNo2C0cDtPd3U1/fz+appGRkcGyZcu47LLLqKmpoaioaEa6k8KL30R/7XXYT+0BzWF9unfE9/JxvvTrUbdpgL3rBFrXqQk/6bte/MXgG/bQ55o46p/B7++FtDh11biy8L3vD9jadbTOkxjFi6a2Mr6/F8fRJ4e1LMAaR3Pt+TkDEyQ4resU9rP6GDtrGzhf+jX+t/1g3OfaTu/H1ndm1LnCCOE88CCBGz43eewXwHlo3aj5XZppQLAfe8MOwgtvictxXS/9elR3rgbYz7yG1tuCmTPBqvA2O/43f4fAmk9ja3sdM2c+Rullsn/dEImaB/ddYCNwCsgB7gFuAm5PRDxgFYgMXUEk0sVWUVHBTTfdhMfjoaCgYMpv3rm51tdsF6lw7OzsxDAMNE1DKcWKFSuora2lrKxs2E7UM8qZabWqZoit+xTaGK0wbE60vjMTJriR84bOP9eBNtCJGa8EN8goVlCspvx4zdsx/vhd3zg/y9D7bU5geDe5homtu2nC59p6W2CsmV1GEK27cdK4L5TW3YQ2xkr3hIPYelsZ3caaxuOOdU3ZnYOFIpNve2LOqyB8AavTpIJEteDKgN8P/tuNNTXgLbquT2Gdl+lhmiZnz56NTrZubm7GNE0yMjKiS2ItWLCA9PTUm5hqGAbd3d3RHQoKCgp4wxvewJIlS1iwYEHSFLDMtHDVaqslNHJB3nAAo2jRxM+dvxzt6ObRb2aabUpvYjPNzJ0/ZiWkqdkIz79ywucaRYvH3sEAO+Gqqyd8brjs8jFXyjAdGZM+92IYFSsxX3902N55AGh2wuVTmMBshLAffwrHsa2YuUXY1DsxSiafbB2uWo2to250ha0RxCiITzVw1EAXzoMPY2t9DaP4EkJXvDfuXeYzLVHz4O5NxHFDoVA0odXV1dHT0wNASUkJ11xzDR6PJ7EtkgQKBAJ0dnZG129cuHAhK1euZNGiRZSUlMiqKkBw5b04Dj4C/t5oF5rpyCBw5V9D+rwJnxu47E4cR58cNSk3WHNt3MZ3Lordif/6fyTt2e+f3x1bs4Ejg8CaT0/83PR5hAsV9rbXhv2sECZ42bsmfKqZV01o8W04jm09f1ybEzMjn9Cl77yYn2hCoYW34Nx9vzX2N5icTUc64cqVE5bqA2CESH/4I9hbX7MSpGYnY/+D+G/6EqErJt6XLnjVfTgPr8f09w2/plbdB2nxW/hB6zpFxgN3oYV81sID9c/g2vdrvO97ICnns12oOb/YMsDmzZv56U9/yu7du6MriNTU1ODxeHC73XEvVZ1pU62i9Hq9dHZ2YpomaWlpLF++PLoTdVZW1qTPnwtivaa03hacu36G48QLmBl5BK/6G2sl9Ek+AGT87g5s7YdGV35qNvr/7rVJi1QSxV73NK49P8fe30awfDmBaz+LOdmalr4+su5fNUaVK4QW3or/HT+Z+PlGGMerD+J65QEIegktupXA1R+HjDjvQu/vw7X3lziObASbk+DldxK88q8nnGwN4Dj8+ODO6cP3EDTtafT/7Y5JC5+0nmZcO3+GvXEnZmahdU2pt8Z1LC39sU9gb9hhjTNG4kUjXLES392/i9txh5rGKspxT1RKJLg3vvGNtLa2kpOTw+WXX051dfXs3N9risZLcKZp0tfXR3d3d3Ry9erVq1m6dOm4RTNz3Yyt/P7DpaMLJ7De9Afu+j1G5cq4x3AxYjlPjlcftFYyGeM+05lJ/2f2TW9wCZa+7jM46raNut10ZeN7y38Q9iTffnJZP142qoAIIh+4Xh1WnRwvM5HgUuIdbevWrei6zv3338/8+fMTHc6MMk2Tnp6e6ETriooKbrnlFpYsWTK9Czr7+6y19DKL5n4Vl2mi9bdjurLANcWWrmaDMRIcgJkxcfdmVMCLraMeI28BpM/guqXBAczefjDSp/TGZ7ommJw8SWtoNjJdWZhoYxSLmFPb6ToR7K6xdyOwOcYce52tUiLBRdY2TBWmadLd3U1PTw+aplFTU8Ptt9/OkiVLKJzu7cwHOkl/8kvYT+4CTcPMLsV327cxKq+a3uMkCXv9s6Rt+wbaQBdgEPLcgv9N35o00YUWvhHH0U2jJombzqzJxzwMg7T1n8ZR/3T0pnD5cqsr6QL2K5uycADX9n+3Sug1jSxHujWuNMlYWHjxm2HTFzHN8KifN3j5XfGLN0GCV7wXx7HN1lYyQ5g2F+HKVQmKamLBpe/BeeBPwyeY211x7xqdaSmR4FLB0O5Hh8NBeXk5t99+O5dddtmkG6RexEHJePgj2M4ePb+IcPcpMh75GN4ProvravWJYGs7RPrjfz9s5QhH3Xa0DX+H745fTfhc/5u/i731VehpHvKCDnzv/b9Jj+va+jUc9U8PSxb2lldI/8t9+O76TYw/xdS5nvpXnEc2nl8IOGQtomxmFhFeMMEOCDYbA+/8GRmPfRJzyFqJRulSgmv/Lm7xJopRsZLANZ/EteunVnm/pmFqdmutxnh+ALkIgbX/gK1dj64WAwZG0WL8N8dxN44ESM6zL6bM5/Nx7tw5DMOgtLSUO++8k2uuuWZGpjfY2g5h62wYs8TZ+cofCNz0pbjHMJOcL/169PqX4QD25pfQupsn3inb4cL7kW3YTryAo247RsECQsveP6XiEuehdaPGszTA3vSiVVIfjzdRfx/OwxvGnui9+34GJkpwgOG+kf7P7sP58m/R+toILXknRvkk1YizWHD1Rwld9m7sp/aQnldMf/HySYtTEsqZju+9/4et/Qi2s8cw8mvn5CRxSXCzUCgU4uzZswSDQXJycrjttttYsWIF8+fPR9O0GSuc0HqaBz/9jbjdCGHrPBH34880W+fJYVVn5+9wovWenjjBDTIWrCUwSXIYZZyxOwB8PXHZvFTznht/onfP6am9iCOd4OqPTXNkycvMKiJ0ye3WAskz8Pc3HYziSzCKL0l0GHEjCW6WiIyr9fb24nA4uOqqq7j66qvxeDwJ21POKLkUwsFRt1vzh5Jz7AEA08De8Bz2um2Y2YXYFt0+6URtgHDFynEmegcxpjJ3KOTHoW/CfupFjHmVhC6/09qVeTLODGtV/ZE0G6TnTf78C2Dmlo/5ad7UbITLl8XlmEJMN0lwSW5oF2RNTQ3vec97WLp0KZmZid8V2ZxXOTgpd0t0XMrUHJiubIKX35ng6MZhGqSv+zT2U3vOT8rd+xv8N32Z0BUTF0AEV34Y5+uPYgbC0Zac6cgguOx9k8/R8vWQ+ce70frOoAUHMO0uXC/9DwN3/Apj/ooJn+pf9VHSdv5kdMHGJe+I3/w5u4vAtZ/F9fyPzk+4RrN2br52koneQiQJSXBJKBwOc+7cOfx+P9nZ2dx2222sXLlyesv6p4n/tm9jlFyK85UHINhPyH0TwWs/O+nKHolir9t+PrkBmGG0UJi0p79NaPFtE8Zt5pTi/cDDuF74kVU1mp5H4KoPE1p6x6THdb34C7Tu5mjrL7JCf/oTX8B739YJxz7sZ49hjboNL0O3d52Y9LgXI3jlBzGySwcnercTLruCwHV/P6dWuhBzmyS4JBGpguzq6sJms7F8+XLWrl3LwoULk3sCts1OcOWHCK78UKIjmRKH/uTo9QYB7A7sjbsJL75twueb8yrxv/U/Yz/u0SdHd21ijXVpvacxc8cfv3M0PDtqjpUG2FpfhVAAHPErZggvvo2BxbeRk5ODb5aMKwkRkcTvnKkhFArR3t5OKBSipKSE973vfSxbtozcRGxFEApYO1anz5tz1VQRpjN9nEm5TH27HdMEfw840qf+nPEeZxpgn+Q17E4YnRutMbgkXeIr4SL7uk11Ir6YkyTBJUBkB+vInLWrr76aa6+9lpqamsR0QQZ9uJ7+d5yH14NpYmYX43/jvxBecN3MxxJnoaV34Dz8OIRHL2UWrr520ufbT+60lqHqbQVNI6Teiv+Wr4Fz4jHR4BV343r+h8N3X9ZsGMVLMLMm3gQ0eOm7cR54YPikXJuTkOcNSTvPKlG03hbSnvwy9mZrOTCj5DJ8b/nO1PbAE3OO/HXMoEh5fyAQoLy8nLe//e1cccUVCV/YOH3TF7A3PBd9A9V6TpO+/jMM3P17a27MHGLmlI8quzcBI7No0lX9be066es+NXyit74JbaAT37t/PuFzg8vvwda8D0fDDkCztslJn4dvgs0/IwJrP4vtzEHsba9bwWoaxrxK/G/8xqTPTSnhIBl/ugetry1aBGRrfZWMP95j7cotrbmUIwluBvT399PZ2YmmaaxevZrrr78+ca21EbS+NmtV8ZHznUJ+nHt/hf9tP0xMYHHiePVBrM3jz9MAW3+btS/WBHt/OV/6nzEmevuxN+5G6zlt7aE2HpsD/9t/TPDsUWytr2FmlxGuvmZqi9o6M/Dd9VtsZw5ia9cx8mswKq6as93IF8pe/wyar3fYXEUNE8J+HEc2TlolK+YeSXBxYhgGHR0deL1eCgoKuPPOO7nyyisTM7Y2Aa3n9JgLr2qYc3Oy9rk6NGOMRWY1G7buUxMmOFtHw9gTve2uyRPcIKNosbUhaKw0DaPs8sn3Jkthtu6mMbueteAAtq747QYukpckuGkWCARob2/HMAwuvfRS3vCGN7B48eKpT8b29+E4vAHbuWMYJZcSuuStk47vXAyjoHbsydqaY05O6DXmX4l54rlh3YzWHSHCk+zAHJ6/Alv7kdFLk4X9GAXuaY5UxMooWWJ9WBvx+zGdmYTLliYoKpFIkuCmSV9fHx0dHbhcLm6++WbWrl1LaWlpTK+hdZ0i4493owV91tYzzkxcO/+bgXv+jJkT22tNWfo8gsvvsVYWHzqh15lOcNVH4nPMBApe+g5cz31/1M7aRkYB5iRJKnjVh3G+/hhmsH/IRO90gkvviMtyWSI24aqrMQo91mozkV25bU7M7BLCnlsSHJ1IBElwF8E0TTo6Oujv76ewsJB77rmHK6+88oJXGUnb9g00X3f0zVMLeiHkx/XMd/C//UfTGPlwgRu+gDGvCte+/0Ub6CJceRX+6z+HOa8ybsdMFOf+34FpDFsVRANsvS3Q0wK55eM+18wpx/v+B0l77gfYT+3BTMsluPJegsvviXvcYgo0jYH3/h+unT/FcXg9mmkQXPxmAmv/btICIjE3SYK7AOFwmLa2NoLBIIsXL+ZNb3oTSqmLWxPSCFsrbIwY49HMMI6GZxg9sjCNNI3Q8r8itPyv4nmUpOA8vH7MnabBSn7BG7844fPN/Fp87/jv6Q9MTA9nJoEbv0hgkt+jSA2S4GLg9/tpb28H4Oqrr+bmm2+moqJieqohNW1w1+cxihi0GH5NpoEZ8FqTkaXKbhTT5hzWPTnMVCdti+QXDlh9z3Fc5UUkP0lwU+D1ejl37hwul4vbbruN6667bvo3EdVs1q7Px7cNK2IwbU6Cl9w++fONMM7d9+Pa9xsI+cjMLsF/05cJL7p1euOc5YIr/pq0p7819n0rPzzD0YjppvW1kbb1a9hPvgCmtQOE/03fmnOb74qpkQQ3gZ6eHrq6usjNzeXuu+/mqquuiusq/v5bvo7t7DFsva1ghkCzY+QvIHDDFyZ9ruu5H1irXQxWB9p6W0jf9EV8afcTrl4Tt5hnm/CiW+Hpb2NiDC8yyS6DtJxEhiYulhGyJnr3tqINTua3N79Exh/fh/cjW+NajSySkyS4EUzTpLOzk76+PsrKyrjvvvu44oorcDpnYJA6I5+BD63H3rgbW2cD4cJFGJWrJu9qDPqGJbcILeTDtfOnDEiCi3K8+iDYHcMmtmuAzd+N7czrGFJOPmvZG3agDXRFkxtgjWmHBnDom6a064OYWyTBDTJNk7NnzzIwMIDb7ebee+/lkksuwTbTi9lqNsI11xKumXxdxOhTBjoYZ1QJTSa4DmM7d2z0qi1gTfTuOikJbhazdZ4cewfy4AC2joYERCQSLeUTnGEYtLe34/f7WbJkCW9961txu91JsYzWVJlZRWOuKm8CRrGa+YCSmFG2DLNhrIne4QtbYWSu83XjOLIR03cOe9GlhN03zsgCz1pPM44jT6AF+gi5b8IoXz5pT4ZRvNiaDjBiWyLTmWlNAhcpJ2UTnGEYnDlzhmAwyIoVK7jtttuoqqqaVYktyu4icPXHce26PzpZG7B2X177d4mLKwkFl94RXVMyOlnbnka4cjVG0aIER5dcbG2HyPjzh6yVQUI+0p2ZGPk1DNz9+7iOZzmOPEHaln8GwwAjhPPl3xFa9Cb8b/7OhEkuXL0GI696cDk2K8mZNgdmZgGhhVJslYpSLsFF5rCFQiGuuuoq3vSmN1FRMf5mk7NF8Kr7MDPyce35OTbvOcJFCv8NX5C1C0fKyGPg/Q/hevb7OE48h+bMIHjZHQTWfCrRkSUX0yRt4+ch0Bft/NaCXmzn6nDu/R+C134mPsf195G25StooSEzP0MDOI5vJXTirYRrbxj/uZqNgbt+h+v5H+I88jiYBqGFt+K/4fMyXSBFpVSCa2trIxwOc80113DrrbdSVlaW6JCmj6YRWnoHoaV3kJOTw4DsvjwuM7cC/9t/hB/IyckhIOdqFK2vFVvP6VEju1o4gPPwhrglOHvjrjF3WNCCXhxHHp84wQGkZRO45WsEbvlaXOITs0tKJbja2lo+/elPz4kWmxBxpdlhrF3Po/fFyTjbB5losrmriFnK7He/aNEivva1r0lyE2IKzOwSjIJaK7EMvd2RTnDpe+J23HD1GmsVnpEc6YQufVfcjivmppRJcDabLa6TtIWYa3xv+yFmRj6mMwvsTkxnBuHy5QRXfih+B3Vm4HvbjzAd6ZjODEy7C9ORRnDZ+whXrY7fccWcJG1+IcSYzPxavB/djqNuO+nBLgbyFcb8FXFf4zRcez39H33aWrYu6CW04HrMgtq4HlPMTZLghBDjc6QRUm9By8nBmMlinIw8QpffOXPHE3NSQhKcUurLwHsABfiB3cCXdV0/mIh45oygF8fRzZi+s9gLFOEF11k7FAghRApKVAvuJuB+YC/WGlP/CmxTSl2q63pHgmKa1bRzdWQ8+AFrGaqg15qUW1DLwF2/lUVmhRApKSEJTtf124Z+r5T6a6AbWAtsSERMs136E5+3dgMfLO3Wgl5sZ4/jfPGXBGU1EyFECkqW/qscrFg6Ex3IbKR5z2HrqI8mt+jtYT/OQ+sSFJUQQiRWshSZ/Bh4Bdg18o7MzEzsdmvyp91uJydH9uwaybT5x73PZrPJOZuAXFNTI+dp6uRcTc1MnKeEJzil1A+A64DrdF0Pj7zf6/VG/5+Tk0OvLKs0hjQyChdiazs8rBVn2tMIXPIOgnLOxiXX1NTIeZo6OVdTM13nKT8/f9z7EtpFqZT6IfBXwBt0Xa9PZCyzne+t/4mZkYfpzATNZm0RUqwIrv5ookMTQoiESFgLTin1Y+Bu4GZd148kKo65wiyoxfuRp3Ac20p6oANf/iJr2SOZJiCESFGJmgf3M+CvgXcBnUqpyLL+fbqu9yUipjnBmUHo0neg5eQQli4SIUSKS9TH+09iVU4+BbQM+fp8guIRQggxxyRqHtws3DZbCCHEbCIDNEIIIeYkSXBCCCHmJElwQggh5iRJcEIIIeYkzRxre3ghhBBilpMWnBBCiDlJEpwQQog5SRKcEEKIOUkSnBBCiDkp4dvljKSU+jLwbeBnuq5/epzHLAAaxrjrLbquPxnH8BJKKfVN4Bsjbj6j63rZGA+PPOdy4KfAaqAD+H/Av+m6Pqeri2I9V6l6TQEopcqB7wJvxVpCrx74hK7rz07wnFS9rmI6V6l4XSmlTgA1Y9z1hK7rt4/znGrgZ8AbgAHgAeDzuq4HLiaWpEpwSqlrgI8Br07xKW8GDgz5vmPag0o+OnDTkO9H7aEXoZTKBbYCO4BVwCXA/wL9wH/FL8SkMeVzNURKXVNKqTzgBeB54HagHXADbRM8JyWvqws5V0Ok0nW1CrAP+b4c2Af8eawHK6XswEbgHHA9UAj8BtCAz1xMIEmT4JRS84A/AH/D6E/e4zmn63pr/KJKSqEYfub3A5nAh3RdHwAOKqUuAf5RKfWDuf5pm9jOVUSqXVNfBFp0Xf/gkNvGanEMlarX1YWcq4iUua50XW8f+r1S6j6gh3ESHPAm4DKgRtf1U4PP+SLwK6XUV3Rd77nQWJImwQG/AP6i6/rTSqmpJrhHlFLpwDHgh7qu/yV+4SUNt1LqNOAH9gD/PMFmsWuA5wbfhCI2A/8GLGDqf5yzVSznKiLVrql3AU8qpR4EbgZOA7/CGiIYL1Gl6nX1LmI/VxGpdl0BoJTSgPuA34+4XoZaAxyOJLdBm4E0YCXw9IUePymKTJRSHwUWAl+d4lP6sLbWuQurL/wp4EGl1AfiE2HS2APci9Xd8VGgDNiplCoc5/FlwJkRt50Zct9cFuu5StVryo21fVU9cBvwY6wxpk9N8JxUva4u5Fyl6nUVcStQC/xygseMdT2dxRpSuKjrKeEtOKWUwioquU7X9eBUnqPr+lmG9/W/pJQqwupC+P30R5kcdF3fNPR7pdRurD+2DwE/SEhQSSrWc5Wq1xTWh9yXdF3/8uD3+5VSi7DetH+auLCSUsznKoWvq4iPAnt1XT8w6SPjIBlacGuAIuB1pVRIKRUCbgQ+Ofh92hRfZw+wKF5BJqPB3c9fZ/yfuxUoHXFb6ZD7UsYUztVYUuGaagEOjbjtMFA9wXNS9bq6kHM1llS4rlBKlQDvZOLWG4x9PRVhFapc1PWUDAnuMeByYPmQr5eAPw3+f6plosuxLsCUMdinfwnj/9y7gOsHHxdxK9bYwYn4RpdcpnCuxrI8xsfPRi8AasRti4GTEzwnVa+rCzlXY1nO3L+uwBoi8AN/nORxu4AlSqnKIbfdOvjcfRcTQMK7KHVd7wK6ht6mlOoHOnRdPzj4/XeA1bqu3zL4/YeAILAfMIC3Y3UT/NOMBZ4ASqn/BDYAjUAJ8DUgC6ukdtR5wppL8g3g/5RS38L6Y/wS8C9zuNINiP1cpeo1BfwQa2zyK8CDwArgs8A/Rx4g11VUzOcqVa+rweKSjwB/Guw9GXrfp4FP67p+yeBNW7B6V36rlPoc1jSB7wO/vJgKSkiOFtxUlAOeEbd9Faultxd4H/A3uq7/cKYDm2GVWJ+GdOARrE841+i6HvkEOew86brejfVJaD7WufoZ1nhAKozXxXSuBqXcNaXr+l6s6sC7gIPAv2N9GLh/yMPkuuLCztWglLuusOafLmLs7skihrSEdV0PY80r9GK1kh8EHsYqzrkosl2OEEKIOWm2tOCEEEKImEiCE0IIMSdJghNCCDEnSYITQggxJ0mCE0IIMSdJghNCCDEnSYITIsGUUv+nlHophsefGJzInhBKqRKl1DcHN/McevtNSilTKbU0QaEJMYwkOCFErEqwVjJZkOA4hJiQJDghhBBzUsLXohQiEZRSl2EtL7Uaa2PFRuCnuq7/bPD+d2Itw7QUa63U3wJfiWzppJT6JvBprNXS/xu4FDiCtcbe80OO80HgY4P3a8ArwBd0XZ9yl+QUf57rgW8Bq4ABrOXJ/lHX9d7B++8F/he4AmtJrWuBU1ibwD4y5HU04F+BvwXSgb9grRX4R6x9vQBeG/z3aWu3K9B1XRsSTpFS6iHgLUAb8J+6rg9dzkqIGSEtOJGqNmBtqPgB4B1YSSoHQCl1F1aCeHHwvn/BSlLfGfEamVh7ev0ceC9WItyklBq6SeMCrOT4XuAerKTynFLKPV0/iFJqLbANa2uRO4G/x9pc83/HePgDwHrg3Vi7S/9pxCruf4+1ePDPB19rAPjekPtbgPcP/v9TWNtdrRlxjF8CBwaP8QzwM6XU6gv52YS4GNKCEylncMPJWuCduq5HWiNPDd6nYa1k/ltd1z855Dl+rDfq7+i6fm7w5gysVt0Dg495Gqsl+PdYq+uj6/q/DnkNG7AVq9X4AayW0nT4LrBT1/W7hxyrGXhKKbU0sivHoB/quv7rwcfsw9pJ+W3Az5VSdqyNOH+u6/rXBx+/RSlVC1QN/jx+pdSrg/cd0nV99xjx/FHX9W8NHuMZrBX034P1gUGIGSMJTqSiDqyW1M+VUj8BntZ1vW3wvsVYG1j+WSk19O9jO1aX3VLg2SG3Pxr5j67rfUqpSAIDQCm1BGvH+muxijMiFk/HD6KUysRqQX1mRLzPY23TshJr5fuILUPiPaeUasPaeQGsJFaG1cIbaj1Wd+NUDT1GUCl1bMgxhJgx0kUpUo6u6wbwJqwuvV8DrUqp55RSK7C28gB4AitBRL4aBm+vGvJSfbquD4x4+TasLVNQSuVgvdlXAf8IXI81RnYAK1lOh3ysnY/vHxGvH3COiBdG7L2ItaFwJJZI12r7iMeM/H4yEx1DiBkjLTiRknRdPwLcoZRyYiWe/wA2Yu1zBtaY2/4xntow5P/ZSqmMEUmuhPO7Na/BarncOng8AJRS86bnpwCsZGIC38RKyiOdjuG1Wgf/LR5x+8jvhZgVJMGJlDZYFbldKfUDrAKMFqAZWKDr+libNY707sHnoZTKxkqQvxi8L2PwX3/kwUqpa7EKT/ZNU/z9SqndgBo63neBTmEluXcCm4fc/o4RjwsM/iutMpHUJMGJlKOUugL4T6ydg+uxuvn+CTig63qHUupzwO+UUrnAJqw3dDfWbs536rruHXypAeDfBxPbaawdiF3Ajwfv3w30Ab9USn0PqzX3TawEOp2+iFVQYmCV9fdijSPejlUEc3QqL6Lrelgp9X3g+0qpdqzdld8BXD74EGPw30asn/1DSqluIDjd0x6EmA4yBidSUStW9eBXsBLY/cBhBlsquq4/iNWKWQ48hDVl4JPAy5xvvQB4gQ8O3vcwVqJ8q67rLYOvcwZrekAZsA6ruvLjwPHp/GEG593dgNWV+DusKRBfxGqRnYnx5X6INR1i6M/07cH7egaP5wM+ilXA8iyw9+J+AiHiQzNNM9ExCDHrRCZ667peNNljZzul1K+wxhFrEh2LELGQLkohRNTgQsl3AzuxuiTfAnwYqwtXiFlFEpwQSWRwsrU2zt2mruvhOIfQD1yHtQxZFnASK7n9V5yPK8S0ky5KIZKIUuoEMF5X4Eld1xfMXDRCzG7SghMiubwda/HnsfjHuV0IMQZpwQkhhJiTZJqAEEKIOUkSnBBCiDlJEpwQQog5SRKcEEKIOUkSnBBCiDnp/wdCXxOC+wPMrAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "idx = np.argsort(x_3[:,0]) \n",
    "bd = trace_3['bd'].mean(0)[idx] \n",
    "\n",
    "plt.scatter(x_3[:,0], x_3[:,1], c= [f'C{x}' for x in y_3]) \n",
    "plt.plot(x_3[:,0][idx], bd, color='k')\n",
    "\n",
    "az.plot_hdi(x_3[:,0], trace_3['bd'], color='k')\n",
    "\n",
    "plt.xlabel(x_n[0]) \n",
    "plt.ylabel(x_n[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4bfe7a3",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\" style=\"background-color: white; border: 2px solid; padding: 10px\">\n",
    "    <b><i class=\"fa fa-info-circle\" aria-hidden=\"true\"></i>&nbsp; Notes</b><br>\n",
    "    <p style=\"color: black\">\n",
    "        (i) In case of an unbalanced dataset, logistic regression can run into some trouble: the boundary cannot be determined as accurately as when the dataset is more balanced. \n",
    "    </p>\n",
    "    <p style=\"color: black\">\n",
    "        (ii) The decision boundary is \"shifted\" towards the less abundant class, and the uncertainty band is larger.\n",
    "    </p>\n",
    "    <p style=\"color: black\">\n",
    "        (iii) It is always good to have a balanced dataset. If you do have unbalanced data though, you should be careful when you interpret results: check the uncertainty of the model, and run some posterior predictive checks for consistency. Another option is to input more prior information if available and/or run an alternative model.  \n",
    "    </p>    \n",
    "   \n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d981ce29",
   "metadata": {},
   "source": [
    "## Generalization to multiple classes: Softmax Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa61bd29",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\" style=\"background-color: white; border: 2px solid; padding: 10px\">\n",
    "    <b><i class=\"fa fa-info-circle\" aria-hidden=\"true\"></i>&nbsp; Notes</b><br>\n",
    "    <p style=\"color: black\">\n",
    "        In order to generalize to mutliple classes, two modifications are needed: \n",
    "    </p>  \n",
    "    <p style=\"color: black\">\n",
    "        (i) We use a softmax (see also Boltzmann distribution in physics), which is defined as:\n",
    "    </p>\n",
    "    <br></br>\n",
    "    <center>\n",
    "        $softmax_{i}(\\mu)= \\frac{exp(\\mu_{i})}{\\sum_{k}exp(\\mu_{k})}$ \n",
    "    </center> \n",
    "    <p style=\"color: black\">\n",
    "        (ii) We then replace the Bernoulli distribution with the \n",
    "    <a href=\"https://en.wikipedia.org/wiki/Categorical_distribution\">categorical distribution</a>.\n",
    "    As the Bernoulli (single coin flip) is a special case of a Binomial (n coin flips), the categorical (single roll of a die) is a special case of the <a href=\"https://en.wikipedia.org/wiki/Multinomial_distribution\">multinomial distribution</a> (n rolls of a die).\n",
    "    </p>   \n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "250279e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n",
      "[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
      " 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
      " 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n",
      " 2 2]\n",
      "Index(['sepal_length', 'sepal_width', 'petal_length', 'petal_width'], dtype='object')\n",
      "     sepal_length  sepal_width  petal_length  petal_width    species\n",
      "0             5.1          3.5           1.4          0.2     setosa\n",
      "1             4.9          3.0           1.4          0.2     setosa\n",
      "2             4.7          3.2           1.3          0.2     setosa\n",
      "3             4.6          3.1           1.5          0.2     setosa\n",
      "4             5.0          3.6           1.4          0.2     setosa\n",
      "..            ...          ...           ...          ...        ...\n",
      "145           6.7          3.0           5.2          2.3  virginica\n",
      "146           6.3          2.5           5.0          1.9  virginica\n",
      "147           6.5          3.0           5.2          2.0  virginica\n",
      "148           6.2          3.4           5.4          2.3  virginica\n",
      "149           5.9          3.0           5.1          1.8  virginica\n",
      "\n",
      "[150 rows x 5 columns]\n"
     ]
    }
   ],
   "source": [
    "iris = sns.load_dataset('iris')\n",
    "y_s = pd.Categorical(iris['species']).codes\n",
    "x_n = iris.columns[:-1]\n",
    "x_s = iris[x_n].values\n",
    "\n",
    "x_s = (x_s - x_s.mean(axis=0)) / x_s.std(axis=0)\n",
    "\n",
    "print(np.shape(x_s))\n",
    "\n",
    "print(y_s)\n",
    "print(x_n)\n",
    "print(iris)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "8d9ed80a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/deprecat/classic.py:215: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n",
      "  return wrapped_(*args_, **kwargs_)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [β, α]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\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:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='12000' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [12000/12000 00:19&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 33 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_s:\n",
    "    α = pm.Normal('α', mu=0, sd=5, shape=3)\n",
    "    β = pm.Normal('β', mu=0, sd=5, shape=(4,3))\n",
    "    μ = pm.Deterministic('μ', α + pm.math.dot(x_s, β))\n",
    "    θ = tt.nnet.softmax(μ)\n",
    "    yl = pm.Categorical('yl', p=θ, observed=y_s)\n",
    "    trace_s = pm.sample(2000, target_accept=.95)\n",
    "    idata_s = az.from_pymc3(trace_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "1c7e1c4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy is: 0.980\n"
     ]
    }
   ],
   "source": [
    "data_pred = trace_s['μ'].mean(axis=0)\n",
    "\n",
    "y_pred = [np.exp(point)/np.sum(np.exp(point), axis=0)\n",
    "          for point in data_pred]\n",
    "\n",
    "res_t = np.sum(y_s == np.argmax(y_pred, axis=1)) / len(y_s)\n",
    "print(\"accuracy is: {:1.3f}\".format(res_t))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "bd3265cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import softmax "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "32eac468",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy is: 0.980\n"
     ]
    }
   ],
   "source": [
    "y_pred2 = softmax(data_pred, axis=1)\n",
    "res_t2 = np.sum(y_s == np.argmax(y_pred2, axis=1)) / len(y_s)\n",
    "print(\"accuracy is: {:1.3f}\".format(res_t2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "57a64211",
   "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>α[0]</th>\n",
       "      <td>-1.135</td>\n",
       "      <td>3.746</td>\n",
       "      <td>-7.921</td>\n",
       "      <td>6.067</td>\n",
       "      <td>0.048</td>\n",
       "      <td>0.039</td>\n",
       "      <td>6080.0</td>\n",
       "      <td>5193.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>α[1]</th>\n",
       "      <td>5.830</td>\n",
       "      <td>3.251</td>\n",
       "      <td>-0.181</td>\n",
       "      <td>11.987</td>\n",
       "      <td>0.044</td>\n",
       "      <td>0.031</td>\n",
       "      <td>5382.0</td>\n",
       "      <td>5418.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>α[2]</th>\n",
       "      <td>-4.846</td>\n",
       "      <td>3.421</td>\n",
       "      <td>-11.267</td>\n",
       "      <td>1.436</td>\n",
       "      <td>0.046</td>\n",
       "      <td>0.035</td>\n",
       "      <td>5554.0</td>\n",
       "      <td>5354.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[0, 0]</th>\n",
       "      <td>-2.604</td>\n",
       "      <td>4.099</td>\n",
       "      <td>-10.736</td>\n",
       "      <td>4.568</td>\n",
       "      <td>0.047</td>\n",
       "      <td>0.041</td>\n",
       "      <td>7711.0</td>\n",
       "      <td>5642.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[0, 1]</th>\n",
       "      <td>1.985</td>\n",
       "      <td>3.267</td>\n",
       "      <td>-4.214</td>\n",
       "      <td>7.966</td>\n",
       "      <td>0.044</td>\n",
       "      <td>0.034</td>\n",
       "      <td>5404.0</td>\n",
       "      <td>4740.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[0, 2]</th>\n",
       "      <td>0.640</td>\n",
       "      <td>3.265</td>\n",
       "      <td>-5.464</td>\n",
       "      <td>6.770</td>\n",
       "      <td>0.045</td>\n",
       "      <td>0.036</td>\n",
       "      <td>5341.0</td>\n",
       "      <td>5130.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[1, 0]</th>\n",
       "      <td>3.186</td>\n",
       "      <td>3.414</td>\n",
       "      <td>-3.514</td>\n",
       "      <td>9.272</td>\n",
       "      <td>0.048</td>\n",
       "      <td>0.034</td>\n",
       "      <td>5167.0</td>\n",
       "      <td>5913.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[1, 1]</th>\n",
       "      <td>-1.009</td>\n",
       "      <td>3.037</td>\n",
       "      <td>-7.135</td>\n",
       "      <td>4.367</td>\n",
       "      <td>0.045</td>\n",
       "      <td>0.032</td>\n",
       "      <td>4579.0</td>\n",
       "      <td>5063.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[1, 2]</th>\n",
       "      <td>-2.408</td>\n",
       "      <td>3.057</td>\n",
       "      <td>-8.328</td>\n",
       "      <td>3.168</td>\n",
       "      <td>0.045</td>\n",
       "      <td>0.032</td>\n",
       "      <td>4636.0</td>\n",
       "      <td>5119.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[2, 0]</th>\n",
       "      <td>-6.333</td>\n",
       "      <td>4.256</td>\n",
       "      <td>-14.192</td>\n",
       "      <td>1.548</td>\n",
       "      <td>0.047</td>\n",
       "      <td>0.036</td>\n",
       "      <td>8195.0</td>\n",
       "      <td>6219.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[2, 1]</th>\n",
       "      <td>-1.395</td>\n",
       "      <td>3.540</td>\n",
       "      <td>-8.308</td>\n",
       "      <td>5.034</td>\n",
       "      <td>0.042</td>\n",
       "      <td>0.035</td>\n",
       "      <td>7265.0</td>\n",
       "      <td>6032.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[2, 2]</th>\n",
       "      <td>7.810</td>\n",
       "      <td>3.743</td>\n",
       "      <td>0.967</td>\n",
       "      <td>15.085</td>\n",
       "      <td>0.044</td>\n",
       "      <td>0.032</td>\n",
       "      <td>7166.0</td>\n",
       "      <td>5728.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[3, 0]</th>\n",
       "      <td>-5.792</td>\n",
       "      <td>4.358</td>\n",
       "      <td>-14.132</td>\n",
       "      <td>2.098</td>\n",
       "      <td>0.051</td>\n",
       "      <td>0.038</td>\n",
       "      <td>7316.0</td>\n",
       "      <td>5879.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[3, 1]</th>\n",
       "      <td>-1.078</td>\n",
       "      <td>3.558</td>\n",
       "      <td>-7.501</td>\n",
       "      <td>5.823</td>\n",
       "      <td>0.046</td>\n",
       "      <td>0.037</td>\n",
       "      <td>5998.0</td>\n",
       "      <td>5723.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>β[3, 2]</th>\n",
       "      <td>6.744</td>\n",
       "      <td>3.683</td>\n",
       "      <td>-0.055</td>\n",
       "      <td>13.745</td>\n",
       "      <td>0.047</td>\n",
       "      <td>0.034</td>\n",
       "      <td>6223.0</td>\n",
       "      <td>5727.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",
       "α[0]    -1.135  3.746  -7.921    6.067      0.048    0.039    6080.0   \n",
       "α[1]     5.830  3.251  -0.181   11.987      0.044    0.031    5382.0   \n",
       "α[2]    -4.846  3.421 -11.267    1.436      0.046    0.035    5554.0   \n",
       "β[0, 0] -2.604  4.099 -10.736    4.568      0.047    0.041    7711.0   \n",
       "β[0, 1]  1.985  3.267  -4.214    7.966      0.044    0.034    5404.0   \n",
       "β[0, 2]  0.640  3.265  -5.464    6.770      0.045    0.036    5341.0   \n",
       "β[1, 0]  3.186  3.414  -3.514    9.272      0.048    0.034    5167.0   \n",
       "β[1, 1] -1.009  3.037  -7.135    4.367      0.045    0.032    4579.0   \n",
       "β[1, 2] -2.408  3.057  -8.328    3.168      0.045    0.032    4636.0   \n",
       "β[2, 0] -6.333  4.256 -14.192    1.548      0.047    0.036    8195.0   \n",
       "β[2, 1] -1.395  3.540  -8.308    5.034      0.042    0.035    7265.0   \n",
       "β[2, 2]  7.810  3.743   0.967   15.085      0.044    0.032    7166.0   \n",
       "β[3, 0] -5.792  4.358 -14.132    2.098      0.051    0.038    7316.0   \n",
       "β[3, 1] -1.078  3.558  -7.501    5.823      0.046    0.037    5998.0   \n",
       "β[3, 2]  6.744  3.683  -0.055   13.745      0.047    0.034    6223.0   \n",
       "\n",
       "         ess_tail  r_hat  \n",
       "α[0]       5193.0    1.0  \n",
       "α[1]       5418.0    1.0  \n",
       "α[2]       5354.0    1.0  \n",
       "β[0, 0]    5642.0    1.0  \n",
       "β[0, 1]    4740.0    1.0  \n",
       "β[0, 2]    5130.0    1.0  \n",
       "β[1, 0]    5913.0    1.0  \n",
       "β[1, 1]    5063.0    1.0  \n",
       "β[1, 2]    5119.0    1.0  \n",
       "β[2, 0]    6219.0    1.0  \n",
       "β[2, 1]    6032.0    1.0  \n",
       "β[2, 2]    5728.0    1.0  \n",
       "β[3, 0]    5879.0    1.0  \n",
       "β[3, 1]    5723.0    1.0  \n",
       "β[3, 2]    5727.0    1.0  "
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(idata_s).head(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c2e3d4",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\" style=\"background-color: white; border: 2px solid; padding: 10px\">\n",
    "    <b><i class=\"fa fa-info-circle\" aria-hidden=\"true\"></i>&nbsp; Notes</b><br>\n",
    "    <p style=\"color: black\">\n",
    "        (i) 98% is the accuracy on our data; \n",
    "        <b>a true test to evaluate the performance of our model will be to check it on data not used to fit the model</b>\n",
    "    </p>    \n",
    "    <p style=\"color: black\"> \n",
    "        (ii) You can check that we obtained a wide posterior. This is a result of the fact softmax normalizes probability to 1. Therefore, when we used priors on the parameters of 4 species, in reality we can \"eliminate\" one species\" from the problem, in that one of them can be calculated from the other 3 once we know their probabilities (again, they have to sum up to 1!)\n",
    "    </p>\n",
    "    <p style=\"color: black\"> \n",
    "        (iii) Below is a suggested solution, that does fix the extra parameters to some value, e.g., zero \n",
    "    </p>\n",
    "    \n",
    "    \n",
    "   \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "3573ce58",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/deprecat/classic.py:215: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.\n",
      "  return wrapped_(*args_, **kwargs_)\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [β, α]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\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:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='12000' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [12000/12000 00:10&lt;00:00 Sampling 4 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 24 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_sf:\n",
    "    α = pm.Normal('α', mu=0, sd=2, shape=2)\n",
    "    β = pm.Normal('β', mu=0, sd=2, shape=(4,2))\n",
    "    α_f = tt.concatenate([[0] ,α])\n",
    "    β_f = tt.concatenate([np.zeros((4,1)) , β], axis=1)\n",
    "    μ = pm.Deterministic('μ', α_f + pm.math.dot(x_s, β_f))\n",
    "    θ = tt.nnet.softmax(μ)\n",
    "    yl = pm.Categorical('yl', p=θ, observed=y_s)\n",
    "    trace_sf = pm.sample(2000, target_accept=.92)\n",
    "    idata_sf = az.from_pymc3(trace_sf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "a7df5514",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy is: 0.973\n"
     ]
    }
   ],
   "source": [
    "data_pred_sf = trace_sf['μ'].mean(axis=0)\n",
    "\n",
    "y_pred_sf = softmax(data_pred_sf, axis=1)\n",
    "res_sf = np.sum(y_s == np.argmax(y_pred_sf, axis=1)) / len(y_s)\n",
    "print(\"accuracy is: {:1.3f}\".format(res_sf))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "fdaf3b4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#az.summary(idata_sf) --- it will complain as one value of mu is 0 by construction "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf11cde",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmpd_df = az.compare({'model_s':idata_s, 'model_sf': idata_sf}, method='BB-pseudo-BMA', ic='waic')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9852437a",
   "metadata": {},
   "source": [
    "## Final remarks: Robust Logistic Regression (extra, for the curious...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a5eca01",
   "metadata": {},
   "source": [
    "Let's take the dataset for the species setosa and versicolor only. \n",
    "Let's complicate the problem by assuming the presence of unusual seros and/or ones in our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "46f49670",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'sepal_length')"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAf2klEQVR4nO3df7RdZXng8e/1EMpNvdoEYgjEjEibJxrNkoU6JSKpq6bL4liNmGJbC5hZdQClzWKCLUqm/sRWIoaiLhhapTC1pbFmRjsyYvFHKUEptJaSylM7RSy/ApiMjQ2FGO/8sfeFk5tz7j3n3nNy733P97PWWcl597v3fp67zznP2Xu/Z++h0dFRJEkqzTNmOgBJkvrBAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIh0x0wEA7Nmzp+VvFebPn8++ffsOdzgzypzLN2j5gjkPgpnKd8GCBUPtps3qPbhGozHTIRx25ly+QcsXzHkQzMZ8Z3WBkyRpqixwkqQiWeAkSUWywEmSimSBkyQVqaOfCUTEacAm4GTgOOCtmXntJPO8GPgY8HJgN3A18P7M9PYFkqS+63QP7pnA3cBvAI9P1jkingV8CdgFvKye7yLgwqmFObFt27axatUqjj76aFatWsW2bdt6uvxNmzaxaNEiFi5cyKJFi9i0aVPP4xmb58gjj+xonn7n3G/N8Z944omTxn84t0G//qbdbuOVK1eycOHCpx4rV66c0f7r1q07qP+6desm7A9Pb7d58+Z1tN1Wr1590DpWr17d0/7dvo66zbl5+UcddVTPX6fdxj+VHLqNaVZ/do2Ojnb1WL58+Q+WL19+ziR9zlu+fPm/Ll++fLip7ZLly5c/sHz58qHx/Xfv3j3a6rF///6W7c2Pq6++enR4eHgUeOoxPDw8evXVV086byePDRs2HLTssceGDRt6Fk+38/Q7534/uo1/Nm6Dfue8ZMmSljkvWbJkRvqvWbOmZf81a9b07L2zYsWKlv1XrFjRk/7dxtNtzv1+nXa7/Knk0O/Pon68zyaqRUPd3vA0In4AvGOiQ5QRcR1wdGa+tqntZcDtwPMz897m/u2uZDIyMsLevXsnjGfVqlXcf//9h7QvXbqUu+66a8J5O7Fo0SIOHDhwSHuj0eDRRx/tSTzdztPvnPut2/hn4zboVrfLX7hwYdtl7d69e9b3h+63W79j6nc8/X6ddrt86D6Hfn8W9eN9NtGVTPpV4G4C7s/MDU1ty4D7gNWZeVtz/yeeeGK01a/gG41Gyw3a7Mgjj6RVDkNDQzz55JMTJ9OBefPmtZ22f//+nsTT7Tz9zrnfuo1/Nm6DbvU759nWfzbGNNv69/s1cThimg2fXUcccUTbAjcrrkXZ7vplnezBHX/88S2/ERx//PGTztuJdkW20Wi0XP5U4ul2nn7n3G/dxj8bt0G3ern8udK/2+3W75j6HU+/X6e9jB9a59Dvz6J+vM8WLFjQdlq/fibwMLB4XNvipmk9s3nzZoaHhw9qGx4eZvPmzT1Z/tlnn91V+1Ti6Xaefufcb93GPxu3Qbe6Xf6SJUtmVfuaNWu6aofut9uKFSv62t5tPN3m3O/XabfLh+5z6Pdn0eH+7OpXgbsNeGVEHNXUthZ4EPhOL1e0fv16tm7dytKlSxkaGmLp0qVs3bqV9evX92T5W7ZsYcOGDU9dSLTRaLBhwwa2bNnSs3i6naffOffb+PiXLVs2YfyzcRt0q9vl79y585Bis2TJEnbu3Dkj/bdv337Ih+KaNWvYvn17y/7Q/XbbsWPHIcVpxYoV7Nixoyf9u42n25z7/TrtdvlTyaHfn0WH+7Oro3NwEfFM4CfrpzuA3wE+B+zOzO9GxIeAl2fmz9b9nw0k8FXgA8By4FrgvZn5kfHLn84gk9KYc/kGLV8w50EwU/n24nY5LwX+tn4MA++t//++evoS4MSxzpn5fao9tuOAO4CPAx8BLu8ydkmSpqSjQSaZ+VWgbZXMzHNatP09cNpUA5MkaTq8FqUkqUgWOElSkSxwkqQiWeAkSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHCSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIlngJElFssBJkopkgZMkFckCJ0kqkgVOklQkC5wkqUgWOElSkSxwkqQiWeAkSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSEd02jEizgcuApYAO4GNmXnLBP1/GXgnsBz4V+AvgE2Z+fC0IpYkqQMd7cFFxJnAFcClwEnADuDGiFjWpv8rgOuBPwRWAm8AXgj80fRDliRpcp3uwV0IXJuZ19TPL4iI1wDnARe36H8KcH9mfrR+fm9EXAlcOa1oJUnq0KQFLiKOBE4GtoybdBOwus1stwKXRsTrgD8HjgbeDHxh6qFKktS5TvbgjgEawK5x7buAV7eaITNvi4g3Ux2SHK7X8yXg7Fb958+fT6PROKS90WgwMjLSQYjlMOfyDVq+YM6DYDbm2/Egk25ExAupDke+H/gi1cCUy4CrgbPG99+3b1/L5YyMjLB3795+hDhrmXP5Bi1fMOdBMFP5LliwoO20TgrcY8ABYPG49sVAuxGRFwO3Z+Zl9fO7IuLfgFsi4l2ZeX8H65UkacomHUWZmU8CdwJrx01aSzWaspX5VEWx2dhzf3snSeq7Tg9RXg5cHxG3Uw0gORc4DrgKICKuA8jMscOPnweuiYjzePoQ5VbgbzLzuz2LXpKkNjram8rMG4CNwCXAN4FTgdMz8766y7L6Mdb/WqqfFrwDuBv4DPCPwOt7E7YkSRMbGh0dnekY2LNnT8sgBu0kLZjzIBi0fMGcB8EMDjIZajfN82GSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIlngJElFssBJkopkgZMkFckCJ0kqkgVOklQkC5wkqUgWOElSkSxwkqQiWeAkSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHCSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIlngJElFOqLTjhFxPnARsATYCWzMzFsm6H8kcAnwq8BxwC5gS2b+3rQiliSpAx3twUXEmcAVwKXAScAO4MaIWDbBbH8CvAZ4GxDAeuCuaUUrSVKHOt2DuxC4NjOvqZ9fEBGvAc4DLh7fOSJ+DvhZ4MTMfKxu/s40Y5UkqWOTFrj6UOPJwJZxk24CVreZ7Q3AXwMXRsRZwOPAjcC7MvMHU45WkqQOdbIHdwzQoDqH1mwX8Oo28zwfOBV4AjgD+AngSqpzcW8a33n+/Pk0Go1DFtJoNBgZGekgxHKYc/kGLV8w50EwG/PteJBJl54BjAK/nJnfB4iIdwBfjIjFmXlQsdy3b1/LhYyMjLB3794+hTg7mXP5Bi1fMOdBMFP5LliwoO20TgaZPAYcABaPa18MPNxmnoeAB8aKW+1b9b8TDUyRJKknJi1wmfkkcCewdtyktVSjKVu5FTguIp7Z1La8/ve+boOUJKlbnR6ivBy4PiJupype51KdT7sKICKuA8jMs+r+nwY2A5+KiPdQnYO7AvhMZj7Sq+AlSWqno9/BZeYNwEaqH25/k2oAyemZObY3toymQ4/1SMlXA8+mGk35p8DXgA09iluSpAkNjY6OznQM7Nmzp2UQg3aSFsx5EAxavmDOg2AGB5kMtZvmtSglSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHCSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIlngJElFssBJkopkgZMkFckCJ0kqkgVOklQkC5wkqUgWOElSkSxwkqQiWeAkSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHCSpCJZ4CRJRTqi044RcT5wEbAE2AlszMxbOpjvVOCrwD2Z+aIpxilJUlc62oOLiDOBK4BLgZOAHcCNEbFskvkWANcBN08zTkmSutLpIcoLgWsz85rM/FZmXgA8BJw3yXx/APwhcNs0YpQkqWuTFriIOBI4Gbhp3KSbgNUTzHc+sBj4wHQClCRpKjo5B3cM0AB2jWvfBby61QwR8WLgt4GfzswDETHhCubPn0+j0TikvdFoMDIy0kGI5TDn8g1avmDOg2A25tvxIJNORcSPATcAmzLz3k7m2bdvX8v2kZER9u7d28PoZj9zLt+g5QvmPAhmKt8FCxa0ndZJgXsMOEB1uLHZYuDhFv2XAC8APhURn6rbngEMRcQPgdMzc/zhTkmSemrSc3CZ+SRwJ7B23KS1VKMpx3sAeDHwkqbHVcA/1f9vNY8kST3V6SHKy4HrI+J24FbgXOA4qsJFRFwHkJlnZeZ+4O7mmSPiEeCJzDyoXZKkfumowGXmDRFxNHAJ1SHIu6kONd5Xd5nw93CSJB1uQ6OjozMdA3v27GkZxKCdpAVzHgSDli+Y8yCYwUEmQ+2meS1KSVKRLHCSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIlngJElFssBJkopkgZMkFckCJ0kqkgVOklQkC5wkqUgWOElSkSxwkqQiWeAkSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHCSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpGO6LRjRJwPXAQsAXYCGzPzljZ93wicC5wEHAX8A/DBzPzctCOWJKkDHe3BRcSZwBXApVRFawdwY0QsazPLGuDLwGvr/l8AtkfEK6cdsSRJHeh0D+5C4NrMvKZ+fkFEvAY4D7h4fOfM/I1xTe+NiNcCbwBa7vVJktRLk+7BRcSRwMnATeMm3QSs7mJdI8CeLvpLkjRlnRyiPAZoALvGte8Cju1kJRHxdmApcH1X0UmSNEUdDzKZqog4A7gMODMz72vVZ/78+TQajUPaG40GIyMjfY5wdjHn8g1avmDOg2A25ttJgXsMOAAsHte+GHh4ohkj4k3AdcBZmfn5dv327dvXsn1kZIS9e/d2EGI5zLl8g5YvmPMgmKl8FyxY0HbapIcoM/NJ4E5g7bhJa6lGU7YUEb9IdUjynMz8TEeRSpLUI50eorwcuD4ibgdupfqN23HAVQARcR1AZp5VP38zVXHbBPxlRIydq3syM3f3LnxJklrr6HdwmXkDsBG4BPgmcCpwetM5tWX1Y8y5VMVzK/BQ0+OzPYhZkqRJDY2Ojs50DOzZs6dlEIN2DBvMeRAMWr5gzoNgBs/BDbWb5rUoJUlFssBJkopkgZMkFckCJ0kqkgVOklQkC5wkqUgWOElSkSxwkqQiWeAkSUWywEmSimSBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHCSpCJZ4CRJRbLASZKKZIGTJBXJAidJKpIFTpJUJAucJKlIFjhJUpEscJKkIlngJElFssBJkopkgZMkFckCJ0kqkgVOklQkC5wkqUgWOElSkSxwkqQiWeAkSUWywEmSinREpx0j4nzgImAJsBPYmJm3TNB/DXA5sBJ4EPhwZl41vXAlSepMR3twEXEmcAVwKXASsAO4MSKWtel/AvCFut9JwIeAKyPijF4EPdtt27aNVatWcfTRR7Nq1Sq2bds26Tzr1q1j4cKFzJs3j4ULF7Ju3bqermP16tUsXLjwqcfq1asn7L9p0yYWLVrEwoULWbRoEZs2bZqw/8qVKw9a/sqVKzuOZ968eZPGs3jx4oOWv3jx4p72h6e3wdij19tg7G80to0n+xt1uw261e3yp/K6lmbS0Ojo6KSdIuIbwF2Z+WtNbd8GPpOZF7fo/7vAGzPzp5rafh9YmZmnjO+/Z8+elkGMjIywd+/ejhKZLbZt28bGjRt5/PHHn2obHh5m69atrF+/vuU869at42tf+9oh7WvWrGH79u3TXsfq1au55557DmlfsWIFO3bsOKR906ZNfPKTnzykfcOGDWzZsuWQ9pUrV/LQQw8d0r5kyRJ27tw57XgWL17M/v37D2mfN28eu3btmnZ/6P826PZv1O026Fa3y5/K67rZXHwvT9eg5TxT+S5YsGCo3bRJC1xEHAnsA34pM7c1tX8ceFFmrmkxz18Cf5+Zb29qWw98GpifmQd9+pRU4FatWsX9999/SPvSpUu56667Ws6zcOHCtsvbvXv3tNfR7fIXLVrEgQMHDmlvNBo8+uij017+bOs/lXlm2zboVrfLn8rrutlcfC9P16DlPBsLXCfn4I4BGsD4r767gFe3medY4C9a9D+iXt5BX2Xnz59Po9E4ZCGNRoORkZEOQpw9HnjggbbtU8ml1Ty9XEer/q0++Mbae7H8udS/3TxzaRu0W043y59uvnPxvTxdg5bzbMy340Em/bRv376W7XPxG9Dxxx/f8pvu8ccfP6VcWs3Ty3W06t9oNNp+u+/F8udS/3bzzKVt0Eq3y59uvnPxvTxdg5bzDO7BtZ3WySCTx4ADwPiz9IuBh9vM83Cb/j+sl1eszZs3Mzw8fFDb8PAwmzdvbjvPmjWHHOWdsL3bdaxYsaKr9rPPPrur9iVLlnTV3m088+bN62s79H8bdPs36nYbdKvb5U/ldS3NtEkLXGY+CdwJrB03aS3VKMlWbmvT/47x599Ks379erZu3crSpUsZGhpi6dKlk56I3759+yEfpO0GN0xlHTt27DikeLQb0AGwZcsWNmzY8NRh40ajMeHghp07dx7yQd1u8MRU4tm1a9chxWmiASPd9of+b4Nu/0bdboNudbv8qbyupZnW6SjKM4HrgfOBW4Fzgf9MNSryvoi4DiAzz6r7nwDcDVwDXA28AvgE1UCVPxu//JIGmUyXOZdv0PIFcx4Es3GQSUe/g8vMG4CNwCXAN4FTgdMz8766y7L6Mdb/XuB04LS6/7uBX29V3CRJ6oeO9uD6zT24p5lz+QYtXzDnQTBn9+AkSZprLHCSpCJZ4CRJRbLASZKKNCsGmUiS1GvuwUmSimSBkyQVyQInSSqSBU6SVKRZcbucMRHxNuCXgJOAZwMnZOZ3JpnnHOBTLSYNZ+a/9zrGXppKvvV8ZwDvB04E/i/w7sxsfVXgWSYifgzYQpX3MHAzcH5mHnovlqfneQ/w2+Oad2Xmsf2Kczoi4nzgImAJsBPYmJm3TNB/DXA5sBJ4EPhwZl51OGLtlW5yjoifAb7SYtILMvPQW73PMhFxGrAJOBk4DnhrZl47yTwvBj4GvBzYTXWN3vdn5pwY5ddtzhHxPODeFpN+PjP/Tz9ibGW27cHNB24C3tPlfPuo3lhPPWZ7cat1nW9EnALcAPwR8JL6320R8R/7EF8/bAXOoCpwrwSeBfx5RBx6x9uDJQdv4xf3McYpqy9MfgVwKdUXlx3AjRGxrE3/E4Av1P1OAj4EXFl/iZkTus25yUoO3qbf7mecPfRMqovJ/wbw+GSdI+JZwJeobvr8snq+i4AL+xhjr3WVc5PXcPA2/nLvQ2tvVu3BZeZWgIh4aZezjmZmu3vTzVpTzHcj8JXM/GD9/IMR8aq6/Zd6GV+vRcSzqe5C8dbM/FLd9qvAfVR3h//iBLP/cI5s4wuBazPzmvr5BRHxGuA84OIW/c8FHszMC+rn36q/rGwC5srFybvNecwjmTnn7g+ZmV+g+lJCRFzbwSy/QvVl9uzMfBy4OyJWABdGxOVzYS9uCjmP+d5Mvm9n2x7cVA1HxH0RcX9E/HlEnDTTAfXRKVR7fc2+CKyegVi6dTIwj6b4M/NfgG8xefzPj4gHI+LeiPiTiHh+H+Ockog4kirH8dvnJtrn1257vjQi2t+hdZaYYs5j7oiIhyLi5vpLWqlOAW6pi9uYL1Id6nvejER0+Hw2Ih6JiFsj4k2He+UlFLgENgCvp9qD+Xfg1oj4qRmNqn+OpTrU0WxX3T7bHUt1d/jx39oni/8bwDlUhzt+re67IyKO7kOM03EM0KC77dNuex5RL2+2m0rOD1Ht3Z0BvJHqPXxzRLyyX0HOsHbbeGxaiX5AdRTiF6lunXYzcENEvOVwBtH3Q5QR8QGq+8FN5FWZ+dWpLD8zb6O6g/jY+nZQ3YPuAuDXp7LM6eh3vrNRpzlPdfmZeeO49X0d+GfgbKrBGZpDMjOpitqY2+pBCRcBbQfjaO6oDz1/pKnpjog4Bngn8D8OVxyH4xzcViZP6Lu9WllmHoiIO4CZ2oPbSn/zfRhYPK5tcd0+U7bSWc4/TfVt/xjg0aZpi+nigy0zfxARO5m5bdzOY1R7qN1sn3bb84ccuqc7G00l51a+Aby5V0HNMu228di0QfEN4K2Hc4V9L3B1JT9sb9SIGAJWAX93uNbZ7DDkexuwFrisqW0t1ci1GdFpzhFxJ7CfKt5P121LgRfQRfwRcRSwgtZDzWdMZj5Z57gW2NY0aS3tB4zcBqwb17YWuCMz9/c+yt6aYs6tvITq0GWJbgN+NyKOahrdvZbqJyHfmbGoDr+XcJi38awaRRkRx1Idk15eN70wIn4C+G5m7q773AzcnpkX189/G/g61RDjZ1EdllxFdYx/VptKvlTDsf8yIn4L+J9UH46vAk49jKFPSWZ+PyL+APhwRDwCfI/qEONdwF+M9YuIe4CPZebH6udbgM9T7QU+B9gM/Djwh4c3g45cDlwfEbcDt1KNkjwOuAogIq4DyMyz6v5XAe+IiK1Uv416BdX5xlk9InacrnKOiI1UH+w7gSOBtwBvoDonN+tFxDOBn6yfPgNYFhEvAXZn5ncj4kPAyzPzZ+s+n6b6Hee19eH85cBvAe+dCyMoofucI+Jsqi+zfwv8CHgd8HbgNw9n3LNtkMm5VH+QP6qf/+/6+S809TmR6vcUY34C+O9UI/FuAo4HTsvM2/sdbA90nW9m7qA6lHMOVWE4CzgzM79xGOLthY3Adqrf8t1KdTL6dZl5oKlPcPAAi6XAH1Odt/ks8ATw05l53+EIuBuZeQNVjpdQnQs+FTi9KdZl9WOs/71UJ+FPq/u/G/j1zJwrPxHoOmeqonYZ1ev3lrr/azPzs4cp5Ol6KdX79G+pLlbw3vr/76unL6F63wLVFzuqPbbjgDuAj1Odn5pL54+7yrl2CVW+f031mbUhMz96WKKtebscSVKRZtsenCRJPWGBkyQVyQInSSqSBU6SVCQLnCSpSBY4SVKRLHDSDIiIa+tLynXa/zv1D95nREQ8JyLeU18zsrn9ZyJiNCJeNEOhSW1Z4CR14jlUV+N43gzHIXXMAidJKtKsuhal1E8RsZLqEkkvB36M6tqWH8vMj9fTX091ncsXAf8PuA5499hFjyPiPcA7qO49eCXwQuAe4B2Z+VdN6zkLeFs9fYjq8lUXZWbHhyQ7zOeVwAeAlwGPU13G7MLM3FtPPwf4FNW1WS+nugHpvwDvar4sVn2B8vcB/wU4CvgM1WXv/hg4oe729/W/X4kIADJzqCmcYyJiG/DzwCPAlsz8RC/zlbrlHpwGyeepbu3yFqrrfV4JjABExC9SFYjb62nvpSpSHxq3jPlUtwa6ClhPVQhvrC+cPeZ5VMVxPfDLVEXlll7ehTwiXkF1geqHgTdRXQvydKqCNt6ngc9RXZj728Cf1HdxGLMReBdVTm+iKpYfbpr+EPAr9f/fTnWH6lPGreMaqjt4rAO+Cnw8Il4+ldykXnEPTgOhvtniCcDrM3Nsb+TmetoQ1cV/r8vM85vmeYLqg/pDmfm9unmYaq9u7HY/X6HaE9xIdYV4MvN9Tct4BvAlqr3Gt/D0xWmn63eAHZl5ZtO6HqC6M/aLMvPupr4fzcxP1n3upLqb9H8CroqIBtVNKK/KzP9W978pIk4Anlvn80RE3FVP+4fM/HqLeP44Mz9Qr+OrVFePfyPVFwZpRljgNCh2U+1JXRURvwd8JTMfqactp7ra/Z9GRPN74stUh+xeBHytqX372H/qm6+OFTAAIuIFwKVUhwSf0zTfcnogIuZT7UFdMC7ev6K6RcnJQHOBu6kp3u/Vtyoa24N7LtUtmz43bjWfozrc2KnmdeyPiG83rUOaER6i1EDIzB8BP0d1SO+TwMMRcUtEnMTTt+b5AlWBGHvcW7c/t2lRP8jMx8ct/hHqWxpFxAjVh/1zgQuBV1KdI/s7qmLZCwuo7oz+iXHxPgHMGxcvVIdRmz3ZFMvYodVHx/UZ/3wyE61DmhHuwWlgZOY9wBkRMY+q8Pwu1T341tZd3kZ1j6vx7m36/zMjYnhckXsOT9+p+BSqPZe19foAiIhn9yYLoComo8B7qIryeA92sayH638XjWsf/1yacyxwGjj1qMgvR8TlVAMwHgIeAJ6Xmdd0sIh19XxjdzpeS3XTXajO0UG1N0XdZzXVwJM7exT/v0XE14FoPt83Rf9CVeReD3yxqf0XxvV7sv7XvTLNGRY4DYSIWAVsobqT+D9THeb7TeDvMnN3RPxX4PqIeBZwI9UH+vOBNwBvysx99aIeBz5YF7YHgU1Ud6i+op7+daq7lF8TER+m2pt7D1UB7aV3Ug0o+RHVsP69VOcRX0s1COYfO1lIZh6IiMuAyyLiUaq7rP8C8OK6y4/qf79LlfvZEfF9YH+vf/Yg9Zrn4DQoHqYaPfhuqgL2CeBb1HsqmXkD1V7MS4BtVD8ZOB/4G57eewHYB5xVT/szqkJ5emY+VC9nF9XPA44F/hfV6MpzgX/qZTL17+5OozqUeD3VTyDeSbVHtqvLxX2U6ucQzTldWk/713p9/w78GtUAlq8Bfz29DKT+GxodHZ3pGKQ5YeyH3pl5zGR957qI+H2q84j/YaZjkabKQ5TSgKsvlHwmsIPqkOTPA2+lOoQrzVkWOGmG1T+2HmozeTQzD/Q5hH8DTqW6DNmPA/dRFbeP9Hm9Ul95iFKaYRHxHaDdocD7MvN5hy8aqRzuwUkz73VUF39u5Yk27ZIm4R6cJKlI/kxAklQkC5wkqUgWOElSkSxwkqQiWeAkSUX6/xFQyvjUgPJoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "iris = sns.load_dataset(\"iris\") \n",
    "df = iris.query(\"species == ('setosa', 'versicolor')\") \n",
    "y_0 = pd.Categorical(df['species']).codes \n",
    "x_n = 'sepal_length'  \n",
    "x_0 = df[x_n].values \n",
    "y_0 = np.concatenate((y_0, np.ones(6, dtype=int))) \n",
    "x_0 = np.concatenate((x_0, [4.2, 4.5, 4.0, 4.3, 4.2, 4.4])) \n",
    "x_c = x_0 - x_0.mean() \n",
    "plt.plot(x_c, y_0, 'o', color='k');\n",
    "plt.xlabel(x_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e439e03",
   "metadata": {},
   "source": [
    "We have some versicolor (category 1) with some unusually short sepal_length... \n",
    "\n",
    "We can fix this with a **mixture model**. We say that the outpu variable comes with $\\pi$ probability of random guessing (0.5 chance for category 1 to be indeed 1), and with 1-$\\pi$ probability from a logistic regression model:\n",
    "\n",
    "$p= \\pi \\ 0.5 + (1-\\pi) \\ logistic(\\alpha+X\\beta)$\n",
    "\n",
    "Notice that when $\\pi=1$, we get $p=0.5$ (random guess), whereas when $\\pi=0$ we get the logistic regression. \n",
    "\n",
    "This model can be implemented with a slight modification of what we saw in mod2_part2. \n",
    "\n",
    "N.B. $\\pi$ is a new variable in our model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "daec8bbe",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [π, β, α]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\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:not([value]), progress:not([value])::-webkit-progress-bar {\n",
       "        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n",
       "    }\n",
       "    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "        background: #F44336;\n",
       "    }\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "      <progress value='16000' class='' max='16000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [16000/16000 00:02&lt;00:00 Sampling 4 chains, 1 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "/Users/cfanelli/Desktop/teaching/BRDS/jupynb_env_new/lib/python3.9/site-packages/scipy/stats/_continuous_distns.py:624: RuntimeWarning: overflow encountered in _beta_ppf\n",
      "  return _boost._beta_ppf(q, a, b)\n",
      "Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 8 seconds.\n",
      "There was 1 divergence after tuning. Increase `target_accept` or reparameterize.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_rlg:\n",
    "    α = pm.Normal('α', mu=0, sd=10)\n",
    "    β = pm.Normal('β', mu=0, sd=10)\n",
    "    \n",
    "    μ = α + x_c *  β  \n",
    "    θ = pm.Deterministic('θ', pm.math.sigmoid(μ))\n",
    "    bd = pm.Deterministic('bd', -α/β)\n",
    "    \n",
    "    π = pm.Beta('π', 1., 1.) \n",
    "    p = π * 0.5 + (1 - π) * θ \n",
    "    \n",
    "    yl = pm.Bernoulli('yl', p=p, observed=y_0)\n",
    "\n",
    "    trace_rlg = pm.sample(2000, target_accept=0.95,  tune = 2000, return_inferencedata=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "7b37f4b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([<matplotlib.axis.XTick at 0x2b1dc96a0>,\n",
       "  <matplotlib.axis.XTick at 0x2b1dc9610>,\n",
       "  <matplotlib.axis.XTick at 0x2a727b760>,\n",
       "  <matplotlib.axis.XTick at 0x296a8e1f0>,\n",
       "  <matplotlib.axis.XTick at 0x2a7265820>,\n",
       "  <matplotlib.axis.XTick at 0x2a4abac40>,\n",
       "  <matplotlib.axis.XTick at 0x2adcbd160>,\n",
       "  <matplotlib.axis.XTick at 0x296a44fd0>,\n",
       "  <matplotlib.axis.XTick at 0x2aad9d1c0>],\n",
       " [Text(-2.0, 0, '3.4'),\n",
       "  Text(-1.5, 0, '3.9'),\n",
       "  Text(-1.0, 0, '4.4'),\n",
       "  Text(-0.5, 0, '4.9'),\n",
       "  Text(0.0, 0, '5.4'),\n",
       "  Text(0.5, 0, '5.9'),\n",
       "  Text(1.0, 0, '6.4'),\n",
       "  Text(1.5, 0, '6.9'),\n",
       "  Text(2.0, 0, '7.4')])"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABN7UlEQVR4nO3deWAkZZn48W9V9ZFzMpnJTOZmDpiXa7hmOEUduVQQERHxYPFWRFldV91F9LfourK7Ki6usii7K4IXiuKxiKDIfcPAwFzvcMx9JZlJZpLupLur6v39UZ2kk3Qn3Ukn6e48H4zp6nqr+klNup+8b72HZYxBCCGEqDT2ZAcghBBCjAdJcEIIISqSJDghhBAVSRKcEEKIiiQJTgghREWSBCeEEKIihSY7gHwkEgkTj8cnO4xRq6mpQeKfHOUcO5R3/OUcO5R3/OUcO0BjY6NVjPOURQ3OcZzJDmFMJP7JU86xQ3nHX86xQ3nHX86xF1NZJDghhBCiUJLghBBCVCRJcEIIISqSJDghhBAVSRKcEEKIiiQJTgghREWSBDdVJbrAS052FEIIMW7KYqC3KCJjiNz/VcLrfg1OhO6Lf4C/YNVkRyWEEEUnNbgpxjq4k/CG32L5LlYqTvSBb0x2SGIieKnJjkCICSc1uKkmFIX0Iu7GsjHRusmNR4yvRBfVv/oAdstG/OZj6L70VojUTnZUQkwIqcFNMaZuNomzrsWvmYnftJzEef8y2SGNC7t1E1V3XYn5zd9ixfdPdjiTJrzuTuz9r2BhsPe/THj9b8d0PqurBeeV+7E69xQnQFFcyRiRB67H/OpK7JaNkx3NpJMa3BTkrrgUd8Wlkx3G+PFSVP/yA5A4BHaI6IHt9Fx2+2RHNTmccP9jAyZzu0BW+zZqfnoJYIHxib/vDszMw8ceYwVzXn2A8Eu/wpt3IqlVHwZ7fOeIjN57LaHXHgAvSfUrDxL7+INTusYuCU5UnmQXpLqxAHwXu31rzqJ2q8Zu2Yi38BTMtHkTFWFRWO3biP75K2A5JM77OqZh/pAy3twTwHeDVmk/hTf3+FG/XujVv4KbwPJdjGUTevnPpMopwRmDvXsNVqobb9FpYI/vx5/duomquz+H5fbgbH8SnAiplR8Y19d0Wjdi9faO9l2sWBtGEpwQFaS6EW/JmTjbn8LCkFr5oazFnB1PU3XXJ8CywbKJX/G70Se5ZIyquz6Bs3cd7tLVJC741rh/gNbcfjG43QBU/+QS4p96ckgZZ+sjAFgE91xDWx8lNUsNKRd64edEHwiaqxNvuoZDx1xER3cbKS8BlkVtuJ4ZMxYTscPgu+BE8ZuWj98PNw4ij36H8PNBTd6bv4qeS24Z19ezD2wJrhVgud3Ye18c19cDSJ74N0QfuQHLtvEbl2IaFoz7a5YySXCTzRhCG/+A3bGN1FEXYhoXj+l0VlcL1qHd+M1HgxMpToxlqOft38Pe/QI1jbNJ1WR/kzv6j1huDwAmVI2z4yncYy4e1euF1/wYZ++LWF6K0NaH8Tbfi3vkBaOOf0S+C266lgqQOJi92KwjwYkGidAOBdtZRP/6z/j4bHRcnnzsWnZ0Po2V5Rb9wqNO4fiODpYccQHRw8/Oei5733qif/pHDD72OV/Dn78yS2Ae0fu+TOjVB/AWrKTnghuCDlDjKPz87X3/3s62x4KxoOPYycoPRfsSnAGwsjdPWvtfoer3f4uVOETyjf+Ae9SFo35N98TL8eedSI2J0910/Lg3iZa6oiQ4pdQbgM8DK4F5wIe01reOcMwK4HvAKcAB4AfAP2utTTFiKhfhp35A5Okfgpsg/PztxD58L1Q3jupczvYnqfrtJ8Gy8RsW0P2+X477h0bJsmz8+Sdh1ddDZ2fWIv68EzEbfo/ldgMGf/ZRuU/Xvg1CUUz9nOwFfA9Mb/dU+j7Yxo0dgnA1JhXU4ExVQ9Zi3tLVJM69Due1h/CWnY23+HVZy+23PO6KJNhje9RhMatmPpY1cM1J3/i0hg7x+2gUq/M5jtoEq+a/kea6hQPKVv32k9ixVgCqf/MJYp9+BgadK6TvIfTyvVipbpytjxFe+3NSKz842quRHy/z38TQ1514nFi+C+FqSHUDVs77n1X3fBG7fQsWEL3vy7hL3ghV00b9un7zMcP+3k8lxarB1QHrgNvSX8NSSk0D/gw8DJwMHAn8CIgB3y5STGUhtPWR9AcsGGOw97866oHX4Wf+u+8vVPvgLpzdzwf3GkRW7lFvB8DetQZPvTVn7SZy/9cIr/8NGENi9T/iHv/eIWVSJ15O6OU/Y+9/BW/ucbjL3zqusQPE3/9rqv7vs2CH6HnbjTnL+Y1LsJr34zcelnX/1nbN7+fMItS+nTnGxsw4AtcauqCybdlMM9DQ3YVXO4vN+9eyvuUZFjUs47RF53HYdIVt2QN7raZiYPyhtZdkF7iJ4LGXhJ7sNdBi8puWYbfqYCNSGySfceQtfgN+4xLsA6+CU0Xq5I9mLWcl4/01cQyWnxrn1Dt1FCXBaa3/CPwRQCl1ax6HvB+oAT6gte4G1imljgQ+p5S6YSrV4tzlb8Zu3Qi+h7HD+Fnuj+TLb1yM2f4klvHAT+Hnqm3kIxkjet9XcFo3Yc74OKh3jP5cpcqycI++CI6+KHeZVJzwi78MrikQffw/syY4qhvp/sDvg5rcBDULRZ74HvbBnYAh/MwtJM/96pAy9s5nqf7VB8F4YDl0v+en+BkdTbZ3vMyv1v0XdUvPptY3uICJ1md9PSt+gPC6XwHBB0fTMe/Er11Aa2wPv3zpv5he1cSpC8/hBAvqRnoHW6H+Gi+mqNfMGINvfAw+vu/h42OMIXbBNwk/9E1MqovEaZ/CS7T3l8cHYzAYjDEYfDpNLbGurvRzfWfHYKBzH5FHvwPJLpKnfAwzZ0XwfEYMAFzwDYi1Yqqmg52Cdk3mmQDs0z5E9P5/xvgu7tHnk0q0QKJlTNegJlFDLBYftkzIDrGw4fAhNfVKMln34E4HHkknt173Av8MLAa2TEZQkyF10hX40xdhd+zAPeI8yPHhkg9TNT343vtEuGbU54o89h+EXr0/6JF139ewG9WwTXgVy4lAuAqSMYxl49c2DV++SB/Udttm7D1rMWo1RGYNLWB8QvoerPS/dnjDb7MmuPBzt4Lxgk4mxiP83K0k3vYdANpie7hz/Q+oizRQE65LN6UNE1PHtiCB934wt2/F1M6ioWomDcwknurivpd/wV9n1bPk0H6UH2b6zKOIdLdQFaohGqrGMy7t3a0cir3GoahPh0kStyxibY8Se6aDHjdOwu0e8L3H7SbpJfD8FK7v4vopPBN8d32373nPd/GNNyDR5PSXv45cJl+PPTf2c/R2dNxxa/A1QY6evYobL/hDxSa5yUpwc4Cdg57bl7FvQIKzbZv6+tF/8E82x3GGj//40d9UzmT2rQ3+UgcIRant3o01d9noztWzv38yZsumxsSDdv0yM+K1z4O5/Gfwp69ghWtx3v7Nov8uXnPNNezbt69ve0mkgy80PwtAz70W/7rnVPa4ddx3330AnHfeeQBcP6+KmaFuDLCnO8J1n/lM3zmam5u5/vrrMdH+LuIWEK6uI1JfT3cqxt1rb6MmWsP06pmYHc/ClsfSAZyBtfDkodfBgsz7Vg6GULT/Hm80GmVazXRaQxEe2/MMd6c6STg7SPzhzSTcbhJuN27mvcnMFsLWB4MvMaE2tjyHF0rQWJPlj6gKUBa9KH3fp7OMb5jW19dPSPwh9TaiO54DDMYOE69fMuobzfZJH6H6lYfB+FhNy+iauaIsb1qPdO1D+o/Yu9bgLn8z/oKhH+oATF8O77mjf7vI12Hnzp0sWrSob/tMZw8Ry8eyIIXN6c0pnvDnEokEvWLnzp0LwMOcyDt4AoPFI9bKvucBtm/fTmdnJ+G6OURIDxMAktXNJA8d4o+bf0rLoT3Mrp1PIpEgsuWx/j+OtjxGYvZxQ+J0sHD6zmThWTaxWDutsT20xnbTGtvN/vhePDPOHWzyZGFjWRYWVsZ3e8B28D87uAfWWyb9vG3ZfU2TQfn+M1tuAisZCzZD0YyWl8xS/Q8yjx5QxrKwYsEfk70xmLo5Qzrl4CawutvT9zOtoINRusWml8n4/yD23mey1GiNwbYdzjv8MkJeVcl9vjY2jq6j3WCTleD2As2DnmvO2CdGwT36Ivxp87Hbt+ItXQ05etblw28+htjHH8SKtVK78BgYoT2/HIU2/h/RP/8/LLeb8Lpf0/3eX4zpHmixbPdnk7IdIngYbHaabM2ihneGnsC2gt6A7w49wnWpxUNK+XOPh1B1MEwgVI0/9zhePbCel/Y+SXPdwv6zhSJYvT0ynew9b/2Zh5Pa9QzbSLLPdtl38AUOtT5U0M9mWzZVoRqqnWqq4wcImeDj3mpcjB2qIRKqoiHaSENVE9OrZjJj5xqm73mJOmMTnXUU4dVfpjZcR1WohpATIWSHcOxw8N0KYVtO0NFlDE1udstGnNceJLrkZLqac3f4srpaIBXDTF88NCEVoOoX78PZ/Xzwp0MoSuwDjw65vWDveJrq316JlerG2CFSx15C8pzrcp5zov6oLnWTleCeAP5NKVWlte5JP3cusBvYOkkxVQR/wariLX8TrcNE67AqdCyNvevZvh6sYGG3bCiJBLfZLOBX7htYZu/mVWcZO8zQ5iMb03f/DcDBz3oub9lZJN4UTN/kHn42sYWr+NMz/0JD1Uxsq3+cm6ltxurYGjyuG/y3Z9Bp4uXO13g+miThpd+yyew9H+siDcyqnces+rlErCqqQrVUh2uJ2FESbjdJP4HTuonIwYMs9GwW+w6z553AtHO+QW1k2oDkVPvoCiwT/P6ZnZuIzTh26NAXY4g88A1Cm+7Gm78yGGSfZXiM1b6Nqnu/BKk4iTddm/V9YrVvofoX7w96eD4TxXnLv+EdcW7Wn9PUzc76fKHsrn39dTpjsHoOYQYlOH/ByaRWXEp43W/wZywj+brPDDlPIayO7YRe/jP+zKV4S980pnOVsmKNg6sDeufssYFFSqkTgANa6+1KqeuBU7TWvSNDfwb8E3CrUurrwHLgH4GvTqUelGJyecvfjNnwO4L2KBtv4SmTHVKfjWYRG71FRENRIDFkv49Nu6mlkaCZrMXkHjflrrgEd8UlADy99Q90u13Mrh04+N0+uL3/cce2Afv2x/fx9M77aYsPnWDZtmxmVDcHCS39VRMOBk9Ho1F6enroSh4inurEtVMcPuNYljcdz4KeXzF72zacdJr2UtATzdLiMLi50x/68eBs/APhF34S1ARf/QvmsRtJvvGLQ8pV/eFq7LZg4unquz5B7KonB87VCTi716abKn1IdeNsfThngiuW5KlXBrPIWDbeglOz/oGBZZFcfQ3J1deM+fWsWBs1P7kE3B6ww8HQl+PePebzlqJi1eBWAQ9kbH81/fVj4IPAXKCvt4PW+qBS6lzg+8CzQDvB+LcbihSPECPyFp1O92U/wW7ZgLfodMy0oXM5lrLvuRdxir0Jg8XTfvYxfJnau1t5avtfmFGTpYbmhLF6x6WlZ8BJuD28sPdRNretHVC2NjINNfN4ZtXOZ2ZNM06OKcm6Egdp62phTt0Czjn8EpY2HkUkVBW8xJFJ7A13g58CO4x75Pn5/dBeAqga8FRox1MDtp2dz2Q91Irt76/1eomgE9WgBOfNO5HgHnYIywlPSO3GXXEp3oKTsRKd+M3HjKm5Mx92ywYgPRDddwm98hdJcMPRWj8I5PxX0Vp/MMtzLwFvKMbri8lj73yW6GM34tfOInH2V0Y9C8tk8ZuPCT5UylCSMI/6K/Iu//j2e3HsECF76Iwalpvse2y8BK/uX8dzex4m4fYPHbAth2Nmn8yxzadkPUcvz3dpi++hqX4Olx57JUtnHD3knpi36DQSb/46If0n3MNeh6uyJzhTNR16OtIBhLPO8OEuexOh9XfR28HCXbo667mSZ/4d0b/+M2CRWvGurLPsm8bD6L7sZzj6bqJHvBFvzsSsdm8aF0/Y4O5gQgODsRxwIrjLpIlSiKESXVTf9XGsVDe2HcJKdtHzzh9OdlQii7b4Xtbve5pZtdlrqSZSi5XsYr/l83jEp2XHvQP2z6tfzMkLzmJadPg/YHrcOO3drZy68BzefMylJLtzryTuHvX2vtlkcol/8G6i916L5XaTOOerWWs33uHnkDzz7wht+C3eotNJnfqJ7K+34l1B5ysvkbu2nowR/ePnsQ/tgnV3Yr3n52OeH7bUmLrZxN//K0Kb78OfuQxvWfY5RSuBJDgxalZPR9BtmaC5Y/C9G1E6ntx+LyE7PKBjSabkkRey5pW72Oh3DKhJ1ITrOXn+m/Ka8SKWPEQ81cXFR38UNesEoqEqkuROcHmpmUHi4v8asVjqlI+ROuVjI5YzIwzUD736AHbn7mDKOzdBeM3tJM/+St7hlgvTuCTnHwKVRBKcGDUzbT7egpNxdq0B45E87arJDklk0dHdxobWNcyqyb0U0NqDG9ngd/Rt25bN0bNWcWzzaYTzWCS1K3GQpJ/gfcd/hnnTFhch6slhamb0b4QimNrKHAA9VUiCE6NnWfS842bslo2YqgbM9IUjHyMm3At7HsPGzll760x0sL6lv2PGnLpFnLLgbBqqZmQtP1g81UWP1837jv9b5tZnn9C5XHiLTie56qOEN/wW+7BTSK3KvpagKA+S4MTY2A7+nGMnOwqRQ3cqxprdD9NYnbsm8tzuh/DTs5jMrJnDOcvelfdA6ZSX5FDPAS477lNln9wAsCxSp19F6vSrginZZLB0Wcv+J50QoiLotrW4fipnr8c9ndvYcfCVvu2T578p7+RmjKEtvoezl13C4saRhykIMdEkwQlRoXzj88zO+6nP0fPRNx7P7Oofvrq08Whm1ea+TzfY/vhejpq9kpPmy2gfUZokwQlRoXYf2kp7d1vfzCKD6ba1HOwJFicN2WFOnPf6vM8dT3URCVVxzrJ35by3J8Rkk99MISrU2j2P52ya7HHjvLj38b7tFc2n5UyEgxljONizn7cufx+1kfJbQklMHZLghKhAvu2yqW0N06tmZt3/wp7HSHrB1Fz1kekcNeukvM/d0dPG0sajWDajPGeAEVOHJDghKlCi+iCe72WdJ/JAfB8v73+xb3vV/NU555MczDceCa+H1UvfUbGrQIvKIQlOiAoUr2/N2uRojBnQsWRe/WLmT1ua93nbu1s5ZvYqZteV18TUYmqSBCdEhUlacZJVXdRFhi4/s7VD0xLbBQQrWa+avzrvmphvfFJeklMWVO7chaKySIITosJ0hlqxjDUkcaW8FGt296/AfeSsE2nIcY8um4M9+zl85gqpvYmyIQlOiApiMLRHdmJ7Q++prW95mniqC4CqUA3HzTm9oHMnvR5OWXhWUeIUYiJIghOigiTtOAk7huUPTHCdiYMD5ps8ce6ZRJxo3ueNp7qYXt1U0P06ISabJDghKkin0wKANWj94TWZ801WN7NsRmHzh3YmOjhlwTkyqFuUFfltFaJCBM2TuwiZgTWzPZ3b2H7w5b7tkxecVVAXf893cWwH1XR80WIVYiJIghOiQiStOEk7jmP6Zy8xxvDcrv6OJYXONwlB55Ijm06kOlxbtFiFmAiS4ISoEF2h/YAZ0DwZT3XS3tMKgGOFCppvspfrpzh2zqnFClOICSMJTogKcTA8tHmyLb6n7/Gs2nl5zzfZK+H2UBWuZYF0LhFlSBKcEBUgZfXQ7XTimMiA51tj/QmuqXZuwec9lDjACXPPyHsqLyFKiSQ4ISpAzGkHhvaebIvt7nvcVFN4gjPG54iZx40tOCEmiSQ4ISrAwfBebDOwluX5Lvu7W/q2C01wPW6c2mgDs+sWFCVGISaaJDghypyHSyzURnjQ/bf27ta+sW91kQaqwzUFnbcz0cFxzafJ2DdRtuQ3V4gy1+0cxBiDNejt3Bof2/03g+GIphVjjk+IySIJTogy1xlqzTpwuy2jg8msApsnE24PNeE6ZtXKxMqifEmCE6KMGQyHwvsI+VVD9rUNqMEVNri7M9nO0bNXSfOkKGvy2ytEGUvacVwrgcPADia+7dKVPAiAbTk0Vs0q6Ly+73H4zMLmqxSi1BR1cItS6irgC8BcYD3wWa31I8OUfx/wRWA5cAj4C/B5rfXeYsYlRKXqHR4wmBvu7ns8s6YZx3byPqfrpwjZEebWHzbm+ISYTEWrwSmlLgNuBL4BnAg8DtyjlFqUo/zrgNuBHwPHAO8AjgZ+WqyYhKh0h0J7B8w92cuN9PQ9LnR4QGeigyOaVhCyh55XiHJSzBrc54Bbtda3pLevVkq9BfgkcE2W8qcDO7XW30lvb1FK/Sfwn0WMSYiK5eMSD7UT8Yd2/3cj/TW4WQX2oEz6CZbLygGiAhQlwSmlIsBK4FuDdt0HnJHjsMeAbyilLgT+D5gJvAf4YzFiEqLSdTuHMAwdHmAwA5ooC6nBGWOwDLKwqagIxarBNQEOsG/Q8/uAc7IdoLV+Qin1HoImyep0LH8GPjC4rG3b1NfXFynUiec4jsQ/Scoh9nA4TDSafXVty7L69tl2kMh6t9tDh3Bsh5A18G2csGJgGwBqwnU01jXlvf5bV+IQC2cuZc7MwnpdZlMO13445Rx/OcdeTJM2g6pS6miC5sh/Bu4l6JjyTeAHwBWZZX3fp7Ozc8JjLJb6+nqJf5KUQ+ypVIpEIpF1XzQa7dvn+z5A3/b+0C4sL4Rr3AHHxEIH+x431cwlmUzmHUt7rI3jZp9RlGtWDtd+OOUcfznHDtDY2FiU8xQrwbUBHtA86PlmIFePyGuAp7XW30xvv6iUigGPKKW+pLXeWaTYhKg4rpUgYXcR9Yf+lZ6wu/oej2YGk0UNh48pNiFKRVF6UWqtk8BzwLmDdp1L0JsymxqCpJipd1vG5wkxjLhzkGDtgKFNjwmn/y/3QmYwCYYHhGX2ElExitlEeQNwu1LqaYIOJFcC84CbAZRStwForXubH/8A3KKU+iT9TZT/AazRWm8vYlxCVJwup21I5xIIJl5O2cEQAQubGTWDG1WGOWfiIMsajy5ozJwQpaxoNSWt9R3AZ4EvAy8AZwLna623pYssSn/1lr+VYGjBp4F1wJ3AZuCiYsUkRCUyGDrDrYT8oR1TEk5/82Rj9ayCxrIlvB6WzZTJlUXlKGonE631TcBNOfatzvKcjHsTokApqwfXShA1dUP2Zd5/K2T8mzFBr8v50xaPOT4hSoXc6xKizHQ7BwGT/f5bZgeTAu6/Jb0epkUbaaiaWYwQhSgJkuCEKDNdoTYsht4nM5gBTZSF9KDsSh7k8JnH5j1eTohyIAlOiDLTFWojnOX+W8rqwbeCjsiW51AfmZ73OT3fY3HjkcUKUYiSIAlOiDLi4+FaSewst88za2+hVFXetbHg/pthTt3CYoUpREmQBCdEGemtoWWTef8tlKzO+5w9bjeNNbOoizaMKTYhSo0kOCHKiGelso5/g9EnuK7kQY6YcdyYYxOi1EiCE6KMeFYq6/03H4+kHQ82DDipqrzPaYzHwulHFCtEIUqGJDghyoSPh8HPfv/NjtE7aiBsqrFNfrORGGPAgua6BcUMVYiSIAlOiDIx7P23jA4mUW/oAPBcetxuGqtnURuRpVVE5ZEEJ0SZ8KwUZBncDQPvv1X5+Se4WPIgy2YcO9bQhChJkuCEKAMGk7ODicEMSHCF1OB8fBbK8jiiQkmCE6IMuFYCg5+1/uZaSTw7BYBlHMImvx6UxhiMjH8TFUwSnBBloMc5lH40/PyTUb826xyV2SS9HhqiM2T8m6hYkuCEKAMxp52c998yOphUFdA8GUseYolMzyUqmCQ4IcpAMMHyyB1MogV0MHFNikXTl485NiFKlSQ4IUqcR4qEHcua4Ax+MAYurZAOJmDJ+DdR0STBCVHiepzO9KNs99/iYAWLlYb8KA75reCd8pJUhaqZXtVUrDCFKDmS4IQocXGnHUz2faNtnowlD3FYw3JZ/01UNElwQpS4rtB+HCJZ9yX6aneFdTBJeD0smXHUmGMTopRJghOihPl4xJ2DhEz2pseBNbj8p9uyLIvZdfPHHJ8QpUwSnBAlrMfqAkzWGUxcK4lrJwGwjEXEz2+At+d72NjMrJlTzFCFKDmS4IQoYd32QXLdgMusvUX8upzrxA05Z6qLudMWE7Lz65AiRLmSBCdECeu0W7FzNk9mDA/wa/M+Z9ztYmmj3H8TlU8SnBAlymDocvYTMtk7mKTsnr7HEb8m7/NawNxph401PCFKniQ4IUpUyurGI4VN9sVLXas/wYX9/FbwNsZgjGFW7byixChEKZMEJ0SJ6nE6MTmm5zKYATW4fBNcwuumoWomNeFCZjwRojxJghOiRMWcduwcCc6zUhjLB8A2Tt4zmMSTnSxuVEWLUYhSJglOiBIVC+0nZKJZ96VG0TwJkPJTLJp+xJhjE6IcSIITogT1TrCcq2aW2TwZMvknOMuyaJL7b2KKCBXzZEqpq4AvAHOB9cBntdaPDFM+AnwZ+BtgHrAP+JbW+rvFjEuIctPjdGFBziVyRtPBxPNdQnaIxqpZxQhRiJJXtBqcUuoy4EbgG8CJwOPAPUqpRcMc9gvgLcDHAQVcCrxYrJiEKFfd9qFc8ysDjKqDSTzVxfxpS3Ds7L0yhag0xazBfQ64VWt9S3r7aqXUW4BPAtcMLqyUOg84G1imtW5LP721iPEIUbZioTYcE861iPfABGfym6KrO9XFYhngLaaQoiS4dFPjSuBbg3bdB5yR47B3AM8An1NKXQF0A/cAX9Jad+U4RoiKZzDEQ+2EcswtaTCDOplk74gyhGUxp25hMUIUoiwUqwbXBDgE99Ay7QPOyXHMUuBMIAFcAkwH/pPgXty7ihSXEGUnacfx8bFz3EFwrWTfIqe2CWHn8TY2xoAxNNXOLWqsQpSyonYyKZBNMIvs+7TWBwGUUp8G7lVKNWut+5KlbdvU1+e/FEipcRxH4p8k5RB7OBwmGu2vhXU7+7Ftm5AVAgtCTvA27e1w4odSfWUjpqZvfybHdgacsycVZ9a0uTTPmLgEVw7XfjjlHH85x15MxUpwbYAHNA96vhnYm+OYPcCu3uSWtjH9fREZtUHf9+ns7KRc1dfXS/yTpBxiT6VSJBKJvu326D6MA67vEnJCuJ4LBE2TAAnT34If8qJ9+zN5vjfwnPH9HNW8ckKvRTlc++GUc/zlHDtAY2NjUc5TlF6UWusk8Bxw7qBd5xL0pszmMWCeUipzzqDl6e/bihGXEOUoFjpAyM8+wTJAyu5PXOEcA8GHHGOSLGqQAd5iailmE+UNwO1KqacJkteVBPfTbgZQSt0GoLW+Il3+Z8BXgB8ppa4juAd3I3Cn1rqliHEJUTY8UiTt2LCrc49mFhMLi1ly/01MMUUbB6e1vgP4LMHA7RcIOpCcr7XurY0tSn/1lu8i6IDSQNCb8pfAQ8CHixWTEOWmx+nESv+Xy8AxcCMPEfB8D9uymVE9uygxClEu8q7BKaVOBL4HnARsB67TWv88s4zW+ibgpmzHa61XZ3lOA+cVEK8QFa3b7hx2gDeYAbOY5JqrcsA5U13MrV+MY09mnzIhJl5eNTil1Dzgz8DLwPnA7cBt6YHcQogi6RvgnYPB9A3+dvxwzrXiMnW7MQ6bvnzEckJUmnz/pPt7YBfwIa21AR5QSs0F/gX403gFJ8RUYjDEnY5hJ082+H2P853BxBjD3HpZwVtMPfneg7sI+HU6ufX6NXCSUmpB8cMSYupJ2nF8y8s5wBvAt7y+x3mv4o2RDiZiShqxBqeUqiaYdWSzUiqz/Ob092OAneMQmxBTSsIeeYa6ATW4PBJc0ktQE6mjLtIwptiEKEf51OAaCVr9fw6kMr52pPfPHJ/QhJhaYk77sL0ngb5VvCG/deDiqS4WTTscyxr+vEJUokK6Vf098HDGdh3wAIzQ6UsIkZdYaD+hESZONhlvt3xqcAm3m0WNasyxCVGO8klw7QRJbK/W+tneJ5VSvdOSHxiPwISYSjxcEnaMqF83TCmTbqJ0wOQ3i4llWcyWFbzFFDViE6XWuht4FTh80K7ePwvXFTsoIaaahBPcfxuuiTLz/lvIRLFGePsaYwBDU82cosQoRLnJtxflb4F3KaUyy78XeFZrvavoUQkxxXTbB0cs4xfYwaTHjTOzeg6RUH69LYWoNPneg7sB+BDwU6XUD4HXAx8gGPQthBijWOgAthn+7TiaDiaHzzx2zLEJUa7yqsFprfcQzBu5iGDV7cuBy7XW941jbEJMCcEA73ZCJvcKAkG5wmpwnp9iQcOyMccnRLnKuxel1voF4HXjF4oQU5MXSuJZ7ogzk2TW4MJ51OCwLFnBW0xpRVtNQAgxOm64O69yhdTgDD6O5dBYNWtMsQlRziTBCTHJUlVdMMIAbx+vfwycsUZcRcB3XOZOW4xjjzwZsxCVShKcEJMsUX1oxPtvA9aAM9GRZzyxPVlBQEx5kuCEmERJL4Eb7h52iRxg4BpweU2ybJhbv2jkYkJUMElwQkyi/fG9wPADvGFwDS7PDiY10sFETG2S4ISYRC1duwfML5lLKqMGN1IHEx8X23NkBQEx5UmCE2ISbWvfhG1G7ggyoAY3QoJzrRSRnnpZQUBMeZLghJgkxhi2H3wF2xt5OGpmghtpFhPfcgl3TxtzfEKUO0lwQkySzmQH3W4Ma4QanI+Lb7l92yP1uASLcGr4QeNCTAWS4ISYJG2xPXmVy6y9WdgjrDgQrCAQSkqCE0ISnBCTZPehrdjWyG/BzA4mlhm+vI9L1K/N676eEJVOEpwQk2RL+yaqQ7UjlsuswdkjvGVdO0mN2zjm2ISoBJLghJgErp9iX9cOqsPDreAdGNxEORyDR603Y8zxCVEJJMEJMQnaYnswxuTVROkW0EQJFlF/5KQpxFQgCU6ISdAay2+At8HkXYMz+FhYRPyaosQoRLmTBCfEJNjWsZmIM/KUW8EQAa9ve7gelK6VotqdPuJ9OiGmCnknCDHBjDFs69hMbaR+xLJJq3+tuKD2ljvB+aTk/psQGfJe0TsfSqmrgC8Ac4H1wGe11o/kcdyZwIPAJq31scWMSYhS05U8SCzZSW3tyLONZA4RsEe4/2YsQ40n808K0atoNTil1GXAjcA3gBOBx4F7lFLDrtmhlGoEbgPuL1YsQpSy1thuLMvKa67IpD24BjccIx1MhMhQzCbKzwG3aq1v0Vpv1FpfDewBPjnCcf8D/Bh4ooixCFGydh3aklfvSYDUkCbK7HxcQiYy4krfQkwlRUlwSqkIsBK4b9Cu+4AzhjnuKqAZ+Hox4hCiHGxp30h1KL+aVr6zmLhWkhpvxojrygkxlRSrBtcEOMC+Qc/vA+ZkO0AptQL4J+ByrbWXrYwQlSblJdnXtZPq8MgzmBhMlk4m2fmWR607sygxClEpitrJJF9KqShwB/B5rfWWkcrbtk19/cg9zkqV4zgS/yQptdh3H9xKJBymuqp/MmTHdgg5Q9+KLkmM5QNgGyfd/d/KWtaxHBqcmUTtoIkyHA5P+s9date+UOUcfznHXkzFSnBtgEfQ3JipGdibpfxc4CjgR0qpH6WfswFLKeUC52ut+5o7fd+ns7OzSKFOvPr6eol/kpRa7K/u0yRTSRKJRN9znu/heu6Qst12V9/jkF+FwQfMkLIGg2d70BMmQXDeVCo16T93qV37QpVz/OUcO0BjY3HmUy1KE6XWOgk8B5w7aNe5BL0pB9sFrABOyPi6GXgl/TjbMUKUva3tm4g6+S1l42au4j3MIqeelaLKq8OZnAYZIUpWMd8RNwC3K6WeBh4DrgTmESQulFK3AWitr9Bap4B1mQcrpVqAhNZ6wPNCVApjDDsOvkxNOL+mo8wOJmG/CohnLedZSRq8rLe6hZjSipbgtNZ3KKVmAl8maIJcR9DUuC1dZNjxcEJUuo6eNnpScaZF85ttJHMOyiDBZWcw1HjTxxqeEBWnqG0aWuubgJty7Fs9wrHXAdcVMx4hSklrbDfDTbU1WCrPJkqAKk86FAgxmMxFKcQE2d7xMo6d39+UBjNgmZxQjhqcj4dtHMImv/t6QkwlkuCEmCBbOzblff/Ns5IYK1hOxzahnB1IPCtJjdcoA7yFyEISnBAToDsVoz3eSlUov5rW0A4m2XmWKwO8hchBEpwQE6AltivvCZYh/w4mANX+yKsSCDEVSYITYgLsOrSFYncw6V0RvMqTFQSEyEYSnBAT4LUDG6gJ55+I8mmi9HGJ+rU4hMccnxCVSBKcEOMs5SXZ27md6gISXOYsJqEcNTjXTsj9NyGGIQlOiHHWGtsDmLzXgDP4pKz+uSpz1eAMhlqvOHP2CVGJJMEJMc72dG7DNybv8gk7DukhAiE/io2Ts6wM8BYiN0lwQoyzLe0bqA7V5F0+YffPAl9tsveQ9PFwTEgGeAsxDElwQowjz/fY3vEKtZH8u/L3OIf6HucaAuDKAG8hRiQJTohx1Bbfg+e7BU3R1eP0rwOXK8H5lkud21SUGIWoVJLghBhHezt39I1Xy0fK6sa3ggVN7RGaIKs8GeAtxHAkwQkxjl5r35D39FwAPU7//bcqrz5rE2SQMC2q/NpihChExZIEJ8Q48Y3Ptg5d4P23gQkuG89KUu3VY8sK3kIMSxKcEOPkQHwfKTdJyM5/ppGejB6UVX6OBEdK7r8JkQdJcEKMkz2d2wq8/5bAs5MAWMYmkqMJ0liygrcQ+ZAEJ8Q4eeXA+oLuv2WOf4v6dTmGABjA5KzdCSH6SYITYhx4vse29rHcf8t+nMEQ9msImeiYYxSi0kmCE2Ic7I/vJeUXeP8tY4B3rg4mBkOdTLAsRF4kwQkxDnZ3bqOA6SfxSPWvAWcson6ulQcMdZ4kOCHyIQlOiHHwyv4XqR7l+LeoX4s9zFtTJlgWIj+S4IQoMtdPsa1jc9HHv4HBwpYJloXIkyQ4IYqspWsXvu/nPf8k5Df+zeDjmLBMsCxEniTBCVFkOw6+QiE5yFg+STvWtx0dpoOJY/LvtCLEVCcJTogi29z2IrXh/O+TueHuvoQY9qtxhpmCyza5Fz8VQgwkCU6IIupOxdjbtZ3qcK5ekEOlovG+x7nuv/l4WFjDru4thBhIEpwQRbSncxsYg23l/9ZyI919j3MlONdKSPOkEAWSBCdEEb16YH1BnUs83wuaKNNydTDxLVcSnBAFkgQnRJEYY9jctpb6yPS8jznQvQ/sYER4yI8OMwWXhW1keRwhClHUd4xS6irgC8BcYD3wWa31IznKvhO4EjgRqAI2AP+itf59MWMSYqLs795HLNlJXV1D3se0xHb1PR7u/ptjQnL/TYgCFa0Gp5S6DLgR+AZB0nocuEcptSjHIW8E/gpckC7/R+AupdTrixWTEBNpZ8erBR/T0rWz73E0R/OkayWolfknhShYMWtwnwNu1Vrfkt6+Win1FuCTwDWDC2utPzPoqa8qpS4A3gFkrfUJUcrWtzxLTQG9J40xtMR2923nrMFZLvXurDHHJ8RUU5QanFIqAqwE7hu06z7gjAJOVQ+0FyMmISZSdyrG7s4tBU3P1dGzn6QXTLBsmxBhU5WjpEW1n3+zpxAiUKwaXBPgAPsGPb8POCefEyilPgUsAG4fvM+2berry3eCWcdxJP5JMlGx72p5hXA4RHVV/vNEHujY2/e42p9G2BnaS9KzXCJWFfXhRmw7+Hs0Gh15LbhwODzp/2bl/HsD5R1/OcdeTCXRLUspdQnwTeAyrfW2wft936ezs3PogWWivr5e4p8kExX789sex3gWiUQi72N2d/T/qkfdOlzPHVLGDfVQ3zObZCKJ7/sAeb1GKpWa9H+zcv69gfKOv5xjB2hsbCzKeYrVyaQN8IDmQc83A3uHFu+nlHoXQa3tCq31H4oUjxATxvVTvNz2EtOi+b8pg/tv/R1Mck+w7FHnyf03IUajKAlOa50EngPOHbTrXILelFkppd5NkNw+qLW+sxixCDHRdh3agmtSBa3eHUseIp7qCjZ8m4hfO6SMwYAF1V7+9/WEEP2K2UR5A3C7Uupp4DGCMW7zgJsBlFK3AWitr0hvv4cguX0eeFgpNSd9nqTW+kAR4xJiXG1uW4tjFTZGLXP8WyhZlXUJHB+XqF87TOcTIcRwijYOTmt9B/BZ4MvAC8CZwPkZ99QWpb96XUmQYP8D2JPx9ZtixSTEePN8jw0tzzEtOqOg4/ZlNE+GkzVZy7hWggZv7pjiE2IqK2onE631TcBNOfatHm5biHK0t2s7Sbeb6VWFDcRu6cqswVVDts6XlqHeaxpjhEJMXTIXpRBjsLFlDbZdWPNkjxvnUCJohbctO0hwgxh8ZPybEGMjCU6IUfJ8l/UtT9MQHX3tbWZ1M1aWt2EwPdeMYRc/FUIMTxKcEKO089BrJLwewk6koOMyO5jMrluQtYxnpZjmDh51I4QohCQ4IUZpQ8uzhKzC12gbkOBq5w/ZbzCARY1XnMGuQkxVkuCEGIUeN86GlucK7lyS8pIciPfPaDcrS4LzcQmbKiJ+9t6VQoj8SIITYhReO7ARz3cLWr0boC2+J11Dg+lVTURDQ8e4uVaChuTcrGPjhBD5kwQnxCis2fVQQUvj9Nre8XLf49l1Q2tvAEaGBwhRFJLghChQW3wvuzu3URcprAv//vg+Xt7/Yt/2gmnLhpTx8bGxqZLpuYQYM0lwQhToxT1P4NgOlpV/E6Ixhqd33t/XPDm3/jDm1S8eUi5ld1Ofmo1NYWPrhBBDSYITogDdqRgv7HmU6VWFNSG+cuAl2uJ7ALAth1Pmn5U1QRo8Gtw5Q54XQhROEpwQBdBta3F9t6CVAxJuN8/vfqRv++jZq5hWNXTuymD2Epsat7B5LYUQ2UmCEyJPrp/i8W33FDyx8po9j5DwegCojUxjRfOpWculrAR1bpPMXiJEkUiCEyJPuvUFupKHqA7nPz6tNbabV/a/1Ld98vyzctb+fMulISXNk0IUiyQ4IfLg+ike2Xo39dHpeR/jG5+nd97ft71g2lIWNgztOQm9zZNQ5xU2cFwIkZskOCHysG7vUxxMHCho7NvmtrUc6G4BwLFCrJp/Vs6yvc2TIVPYvJZCiNwkwQkxgniqiwe3/J7Gqll5HxP0tnysb/vY5lOpj+YeN+dbLo2p7AO/hRCjIwlOiBE8svVuXD+VdVqtXNbsfpiUnwCgPtrIMbNX5SxrMFhY1ErvSSGKShKcEMPY2r6J53c/wsya/Dt/7OvayWvtG/q2T5l/1rBzVnqhBNNSzTgUvjKBECI3SXBC5BBLdnK3/gnTojOwrfzeKr7xeGrnX/q2FzUsZ960xcMeYyyf6dI8KUTRSYITIgvf+Nz78s/pcePURurzPm5T6/Mc7NkPQMgOs2r+6mHLJ90ebD8ka78JMQ4kwQmRxZPb72Nz24vMrM6/aTKe7GTt3sf7to+bc/qIyfFgYj81h2Zhy1tRiKKTKRPElBaPG/buhQULIBIJ5oZcv+9ZHt76B2bXLsCyLPbtgyeeMoTD8PrXWUzLMtF/wu3hiR334fopABqqZnLUrJOGfW1jDL7xqY7PgIbB+8DzwkRk1IAQoyYJTkxZ27cbPn6VwfNgxgz4nx/A3sR67ta3MaN6Do4dwvfhT/cZ/GAcNn/6s+Hdl/RPkmyMYWvHJp7d9SA9brzv+VMXnI1tDb8iwKFEO4umL2drasOA55PJKh576v3E4400NLRy2qqfEQqliveDCzFFSIITU9av7zLEusAA+/fDnQ+8SEvD/9JQNbNvSEAq1Z/cAGJd/Y8PJdp5asdf2Nu1fcB5j5x1Es11C0d8/R43xqr5q3mYgQlu244TiMenY4xN56EZ7NpzNIctXDvqn1OIqUoSnCg7XV2GvfvgsEUQDue/Jttgs2dDJAKJpIHZj7DevZPDq2ZRFeqfazIUsiC9hhtAOAKe77K+5Rle2vcUvvH69tWE6zh5/lksbDh8xNfucePUROpZ3KiG7LNtF2OCn8s3NrYttTchRkMSnMhp127D7t1w2qlm5MIT5NXXDFddHdSqmpqCZsWamtEluXe/y2L77hjPtP6GWUc9hVowj7Az8KaX48CcObB3b7A9//Ad/J/+C4cS7X1lLCzUrBM5Yc7rhhyfy8Ge/Zyz7F1ZJ162rP4qozFgZ2wLIfInCU5k9cyzhn+81uA4ML2hk1v/Z/SJpJh+eachHg8++Nva4Mmn4azVQ8sZY9iyFZqafKZl6chojGFn1yYaXvczTk910VSzKOtYN9dNJzcnjj3vYXaGN0Cif//M6mZOXXguM2ua8/4ZuuJJjBfh6Obss5u4bhTLMkFys31Sbv4zqIzFa68Ff8gsXTr5/85CFIMkOJHVnb8xJNIf5JZlWPuixemnTW5MAHOaIRKGRDJIck05Jt+/6mrDS+sAOrn8/XDlx4LkZYxhb9d2Htt2D5tb17HhhUYOtMxl2TJYldHp0Tc+B7pb2HNoB/binVh1O7GcZN/+sB3hhLlnsrzp+LwHgQM8/4LhxVdbMLvO5aiuKt5+4dAyixa8xNbtJ5BKVVEVjTFv7sa8zz9aN93s8+u7gscXX2T49FUybEGUP0lwIqsjlsGzz0EiAZ4H8+dNdkSB97/XoqXFsH4DvP1tcNyKobWNlhY/ndwCP/0ZfPhDPWzt0Dyz86/sPPQaEbuK1i0L2b3dwvdhwwaP6hktmOod7OvaSWtsNyk/SGj2oC78h01fzqr5bypoZQEA34e16xKYSBh35xu46QdkTXDB/Tcb2wKDDWb8a1S/vDOorQL86tfw6avG/SWFGHdFTXBKqauALwBzgfXAZ7XWjwxT/o3ADcAxwG7g37XWNxczJjE6H/yAhecbNm+Gy99fw6JFPaM+V1eX4Z++anj5VbjoQvjIh7LXDpJJw1NPw7RpcPxx2T/UIxGLf/jC8B/4vR/UhGJY017Dmvk8//nkSyTcbsDCxqIj2caW2BaY24UdPQA1u1nT4UJH7vOankaW163mtMVLR/6hs7AscGpbSb56Ibh11OcYA75j5wqSySrAIZGA3XuOYvFhz4/qNfPVNBP2tfQ/FqISFC3BKaUuA24ErgIeTX+/Ryl1tNZ6e5byS4A/Av8LXA6cCdyklGrVWv+6WHGJ0QmFLD7xsSCR1NeH6ewcfYL731sNa16AVAp+fgesWmmGJDDfN3z6M8F9M2PgPZeluPxvXFwvRcpP4vopUl4K10/SGU/Stj9JXWOMpB8nnuqiO9XFwZ4DtHe3sK9zL5Hzd0GkA0I92KEe/vJKd1+NrE8t2LW5464J1zG7dgGvvrQA07UQktPRtsVpx47uOnS7XZx5aj0v7zyD6BFw7TXZE3W0KoZj+3i+g2UZotFYznN2ds7Esgx1dQdGF1Tad75t8R/fDe7BfeZquQcnKkMxa3CfA27VWt+S3r5aKfUW4JPANVnKXwns1lpfnd7eqJQ6Ffg8MKUS3D2bf8aDW36P77vDljPk15vRmP5yA4/JeN4M3N9/jBnyvO3YeK7b/xwGTLBlMOljg+8+Psb4fc8bY9gf8bEvMEQsH2yPr77gE1rv4RsPz/jBd8+j+0QPTvLBcvmF8fnFbXn9uFnZTQO3U3l0RKwJTWNO/UKa6xbQXLeAukgDnZ3wSkbu8Ic5z8aNht17YekSWLJ4YJIwxnCw5wDvPfUjqLcNk1WBRQte5FBnE21tS1gw/xXmNG/OWm79xjexbccJACxb8hTqiMezlsvHwgUW3/53SWyishQlwSmlIsBK4FuDdt0HnJHjsNPT+zPdC3xAKRXWWk+JwT8tXbu44dG/zzt5laUQWNOg9+PzkAtkyeXWON8Rti2b6lAdNeE6aiJ1uN317NhSi0nWY/XM5ey3NjB90P22+vqB/y6hHCvavPIqPLsmaB7dvRtqqqE5o2Nle3cLSxuPZHnT8SPGaVmGFUffD0A0Gu3r7JPJGNiybSW908m+uuXUMSW4RMJwz5+Cx299C0Sj2ZPdq68ZXlgLxx0LRxwhCVGUtmJ9pDQBDrBv0PP7gHNyHDMH+Mug5/alY2oC9vQ+ads29bluWJQBx3Fyxh+tXsLsuvns69o5wVGVPsdysO0QjuXg2A625dB5yMb4DvgOxo9QWxVh/pwqQk6EiBMh7EQI2xHCTpSwEyYaqibiVBG2I1iWRdJP4HpJwObxJzz8WLr26XSxdUeEU2Y1YFkDP7gve7fPXx9IEY5YnHd2OOuH/8GDKVw3GPRtDHR2OSxaFLy9elJxwpEIl5z0UaZVD53IcsGCBezZs2fAc51dS+mKLWVa/cvU1gT7Vq9eDdBXNhTqwnXr0o8PDjjHggULCnrPfO4LXax90cMCHnzY4X9+OLQDzSbt8YmrujAmuJ94y801HLci90fIcL/35aCc4y/n2IupLHpR+r5PZ2fnZIcxavX19cPG/6N3PsqG1ufwzcjtaBb5/dU88DO6f6OnG75yHfQkIOxYnHACfOgKe0A5y7L4/n/5rFsfbEfCFp/4qM2KY+2M81pYlkXvf1gWNnbwXPoxloVt2VhYvLTOYscOi5NXOsyf62DbTpDALBvbCpLX+vU2t/y3ReO0MF/8fIimmUN/1m9/x+eu3/Vvf/MGWHnS0E4rGzcZbrrZMH16iM982qOpqf9cnu+R9Hq491tx3J5urHAnVt029s54hh0HXiNkh5heNQvHDuaSrKmGt50PQfNtMmuN6rDDYMPG4BLaFsyd45FIeHi+S0tsF+84+iM4bjTr78GXvvSlAdvPrTF88RpDKgUHOuD737VQy4dei09e7bNuXfCvduIJTXz7328csL+Q98xza3y89KQsa573sh770MMG1w1qqZYFDz0cH9IUm2mk3/tSV87xl3PsAI2NxVk+qlgJrg3wgMGjXZuBvTmO2ZujvJs+35QRCVVxwtzXTchrbd1mMC0G0wNJYF8YVswZmiBm2z6mNXicAI6bB0fOHt3YqD//xefb3wo+GO+8DX5yq0XT7IEfjN3dhmuvMXT3gG3DP33V8P3vDv3wvPpTFl1dhhfXwSUXZ09uqZThs58zxOLgOC4trfCD7/efy7Edqu1aTE81xIK0ZTqOomH6m7niyl08v/tR1rU8BdjMrJ7N3r0hHn/CEArB6jdYZHvvNc2Eiy+yONAOs2dBVVVw3601tpvTF70Z1XRC3tfruTX9YxB9H9a+CGr5wDK+b1i3LqgtGmDNGDtZHnsMbNwUPFZDZw/rK+M4wb9jJBJsC1HKijKaU2udBJ4Dzh2061wg142BJ3KUf3aq3H+bDAvmw8IFQa2kKgqXvjN7uZ6MTpM1NdDeMfr7LQ8/Epyvt/t+7wdpplgc0i18+D7sHdzYnbZzJzz+JBw6BP93d7DczWDd3cFAcAjG8O3dM6QIAKtWDtx+8zkWc+oX8Vb1Pj5+8j+xct4b2B9v4f4ndtEVS9HRAX99MPe90s4u6OgIfhbf+Ozt2s4xzadw5mHnD2n2HM6qlRbRaJBMbBtOzHLbzrYtli6BUCj4GpwAC/Wtf7O48uPB1w05Opsct8Li36+3eP974fqvW6xaKffgRGkrZhPlDcDtSqmngccIeknOA24GUErdBqC1viJd/mbg00qp/wB+ALwO+CDw3iLGJAYJhSxu/n7wF/+MGWRt+gK49BKL554z2A40Trc4/rjs53Ndw4svQUMDLMsxxdOSJfDAQ8Hj7m5YstjAoKbWmTPgzDPgiaeCBPfhD2Z/vcypulpzTNU1bZrFG15vePyJoNwVl2c/1xf/3mLbNsO27XDqKXDh2/pjaqiawVnLLua4Wau5+38fgjkPg+WS8GYAQ3tB7toF9z8QLL3zwro4p72xldVHruaspe/sa+rM10knWtzwTXhpHbzxDbUsXBDPWu6737G441fBdGqXXTq2ZFNdbfHud41cbuVJFitPksQmykPREpzW+g6l1EzgywQDvdcB52utt6WLLBpUfotS6nzgOwRDCXYDfytj4MZfNDrytFurVlrcdmvwwX3aqfV4XteQMr5v+MznDJtfDpLS1VcZ3nHR0EaBHRn9ZyIR2L7DYtGigWUsy+Jr18FrW6C2FuY0Z/8QbW6GcBiSI0zV9dX/Z/HyKzBrVi2N07MniIYGi5/8ePgP66ZpjfzNG9/O7XechdX0LCde9CAtXTvBsqgO1RJxqnBsh1e2u3ihOFTF8bwq5nZ+mHOWnVRQzS3T8ccFf1TU1zvkupXS0GDx8Y9KshEil6J2MtFa3wTclGPf6izPPQQMv+yxmDTz5lrMmxtMspztQ3bPXtik6btf9PM74B0XDS3XnF6WJpkMOifMmJH99SzLYtkIk4S8/70W+/YFU3VdmGOqrt5zLT9i+ASRr498yObSS+oJhd5EVfVqWrp2sq1jM9s6NtPe3Uo8FcdNRfAPLsW0H4vZfzy1R1ePOrkJIYqjLHpRitLUOD24TwTBfaAlS7KX++AVFhs2GF7dAu+4EI4+avQf/JGIxT9+ceITx7Rpva8Z3KebU7+IUxf2j4D52TafH/5fcJ8xGpEZ+YUoBTJluBi1mhqL795gccbpcP5b4cs5pp568GHDhk3Q1QW//DXsa6m8Qe2XXGzxhjNh5ky48EJ4w5ljP2fmjDRCiMJJDU6MyZFHBj3rhvPYY/29Mi0LNm0Kmi0rSTRq8bXrilNri8UMf/d5w6ZNcPLJXVz/dUMkIjVCIQolNTgx7l53RjAuzHGCjiFHHjnZEZW2u35nePkV8A2sfdHjvsHz/Qgh8iI1ODHu3vJmm4YGwyuvwuvPhObZUhsZzuC+KdJXRYjRkQQnJsTpp5XGiuDl4OKLLB59zLBhI6w8McR553iTHZIQZUkSnBAlpqbG4r++17sWX21ZzykoxGSSe3BCCCEqkiQ4IYQQFUkSnBBCiIokCU4IIURFkgQnhBCiIkmCE0IIUZEkwQkhhKhIkuCEEEJUJEtmLBdCCFGJpAYnhBCiIkmCE0IIUZEkwQkhhKhIkuCEEEJUpEldTUAp9SngE8Di9FPrga9rre/O49gjgDWApbWuG7cgh4+h4PiVUu8GvgQsB1qB72mtvznOoY5IKXUN8A3g+1rrT+dRftKv/6B48oq/VK6/Uuo64J8GPb1Paz0nj2Mn9dqPJvZSue4Z8cwF/hU4H6gHXgM+qbV+KI9jJ/v6Fxx7qVx/pdRW4LAsu/6otb5ghGObgLXAPGCW1rptpNeb7BrcTuAfgJOAVcBfgd8qpY4b7iClVAT4BfDwuEc4vILiV0q9FfgZ8EPgWOAq4O+UUiMmlPGklDoN+DjwYp7lS+X6A/nHX4LXXwNzM75WjHRACV37vGMvteuulJoOPAZYwAXAUcDVQEsex07q9R9N7CV2/U9m4O/NSYABfpnHsT8CXijkxSa1Bqe1/t2gp65VSn0SOJ3hP6z+Lb3/IeCN4xTeiEYR/98Af9Ba35Tefk0pdT3wD0qp72utJ3zMhlKqAfgp8GGG/lWeS0lcfyg4/lK7/q7Wem+Bx5TKtS8k9lK77l8E9mitr8h4bkuex0729R9N7CVz/bXWrZnbSqmPAIcYIcEppT4D1AD/QlBzzUvJLHiqlHKAS4E64PFhyl0AvA04EXjXxEQ3sjzjjwI9g57rBhYQVNu3jld8w/ghcKfW+gGl1IgJrgSvfyHxl9r1X6qU2g0kgKeAL2mtX8tVuMSufSGxl9p1fwfwJ6XUHcCbgN3AfxM0b+f8sC+R6/8OCo+91K4/AEopC/gI8BOtdfcw5U4kaCk7GTiikNeY7CZKlFIrlFJdBG+Um4GLtdYv5Sg7D7gFuFxr3TWBYeZUSPzAvcBFSqnzlFK2Umo58PfpfXMnINwBlFIfAw4Hvpxn+ZK6/oXGT2ld/6eADwJvAT4GzAEeV0rNzFa4xK59QbFTWtcdYClBM91rwJuBGwnuaX0q1wEldP0Ljp3Su/69zgWWEFzXrJRStQRNwldrrXcV+gKlUIPTwAlAA8FfRT9WSq3WWq/LUvZ24L+01k9NYHwjKST+W4BlwO+AMEHV/EbgOsCfiGB7KaUUQaeMM7XWqTwPK5nrP8r4S+b6a63vydxWSj1J8KH1AeCGLIeUzLUfRewlc93TbOBZrfU16e3n0x1HPgV8L8cxpXL9RxN7qV3/Xh8DntFarx2mzHeBR7XWvx7NC5TcVF1Kqb8A27TWH8myzwBexlMWwT+4B1yltf7hxESZ23DxZ5RxCP7qbQXOBv4IzB7cPj2elFIfJLhpm3k9HYIbvj5Qq7VODDqmZK7/aOLPOHbSr3+OuB4ANmmtP5llX8lc+2yGiz2jTElcd6XUNuDPWuuPZjz3N8DNWuvaHMeUxPUfTewZ5Uri+qdjmU3QSe9TWuvhanBbgYUE72sYeN3/TWt97XCvUwo1uMFsgjbjbAb31LoIuBY4BSi4+jpOhosfAK21RzpepdR7gScm4Zfst8Czg577EfAyQc0omeWYUrr+v6Xw+IGSuf4DKKWqgCOBB3IUKaVrP0AesQMldd0fA9Sg55YD24Y5plSu/2hiB0rq+kPQxJ0Afj5CufOASMb2ycD/AqsJ3uvDmuxxcP8K3A3sIBjP8T6CwC9I778eOEVrfTbA4GY/pdQqwM/RHDjuCo0/PY7jUuBBgiT4ofT2hPfG0lp3AB2ZzymlYsCB3utZytd/NPGX0vVXSn0L+AOwHZgNfAWoBX6cLfZSuvaFxl5K1z3tOwT3DK8F7iDoNPK3BOPEgJK+/gXHXmrXP9255KPALwbfz0wPXfi01vpIAK315kH7m9IPN5XDOLg5wE8I7mPdT5Cd35rRxj+XoO24VI0m/iuAZwj+EjsGWK21fnpiwi1YqV//kZTy9V9A8NerBn5D8NfsaVrr3r/ES/najyb2UrnuaK2fIeiN+G5gHUHX868AN2UUK8nrP4bYS+b6E1QCjiB755ImhtZQR63k7sEJIYQQxTDZNTghhBBiXEiCE0IIUZEkwQkhhKhIkuCEEEJUJElwQgghKpIkOCGEEBVJEpwQE0wpdatSavAsLMOV35oeXD0plFKzlVLXKaUWD3p+tVLKKKWOnaTQhBiWJDghxEhmE6y1t3iS4xCiIJLghBBCVKRSnGxZiKJTSh0DfJtgctwowTyK39Nafz+9/yKCKY+OJZjj8jbg2t6leJRS1wGfJphk9z+Bo4FNBPPmPZrxOlcAH0/vt4AXgC9orfNukszz53k98HWC6eG6CabM+pzWujO9/4MEk08fR7CEzRkEc6Z+SWv9m4zzWMDXgE8AVcCdwH0EU3EtSRfrXd/wgWCVItBaWxnhNCmlfgW8FWgBvpWxerQQk0ZqcGKq+APBEhuXA28nSFL1AEqpdxMkiKfT+75KkKSuH3SOGoK5R28mmKy2A7hHKTUno8xiguR4KcHk2zuAR5RSS4v1gyilXgf8BdhLsAbhZ4HzCRLaYD8Dfg9cTDD7+i+UUgsy9n+WYKLem9Pn6gb+PWP/HuD96cefAk5Pf2W6BVibfo0Hge8rpU4Zzc8mRDFJDU5UvPQM5EuAizJWW78/vc8CvgncprW+KuOYBMEH9fVa6/3pp6sJanU/S5d5gKAm+FngHwG01l/LOIcN/Jmg1ng5QU2pGP4VeFxrfVnGa+0C7ldKHTtohvvvaK3/N13mOWAf8Dbg5vT6YF8kWEvs/6XL36eUWkKwBhda64RS6sX0vg1a6yezxPNzrfXX06/xIHAh8E6CPxiEmDSS4MRUcICgJnWzUuq7wANa65b0vuXAIuCXSqnM98NfCZrsjgUeynj+rt4HWusupVRvAgNAKXUUwXp0ZxB0zui1vBg/iFKqhqAGdfWgeB8FUsBKglnme92XEe9+pVQLwWoAECSxOQQ1vEy/J2huzFfma6SUUi9nvIYQk0aaKEXF01r7BAsn7iVYLHGvUuoRpdSJBMtzQLC6cSrja0v6+YUZp+rSWncPOn0LwfIkKKXqCT7sFwKfA15PcI9sLUGyLIZGgpXLbxoUbwIID4oXBq2ZR7AQbG8svU2rgxe9LHQRzOFeQ4hJIzU4MSVorTcBlyilwgSJ598IFqs9N13k48DzWQ7dkvG4TilVPSjJzSa4TwVBzWoBcG769QBQSjUU56cAgmRigOsIkvJguws4197091mDnh+8LURZkgQnppR0r8i/KqVuIOiAsQfYBSzWWmdbgHGwi9PHoZSqI0iQP0zvq05/T/QWVkqdQdDx5LkixR9TSj0JqMz7faO0gyDJXQTcm/H82weVS6a/S61MlBVJcKLiKaWOA74F3AG8RtDM9w/AWq31AaXU3wO3K6WmAfcQfKAvJVg5+V1a63j6VN3Av6QT227g80AEuDG9/0mgC7hFKfXvBLW56wgSaDF9kaBDiU/Qrb+T4D7iBQSdYDbncxKttaeU+ibwTaVUK8Fqz28HVqSL+Onv2wl+9g8opQ4CqWIPexBiPMg9ODEV7CXoPXgtQQK7CdhIuqaitb6DoBZzAvArgiEDVwFr6K+9AMSBK9L7fk2QKM/XWu9Jn2cfwfCAOcDvCHpXXgm8UswfJj3u7g0ETYm3EwyB+CJBjWxfgaf7DsFwiMyf6RvpfYfSr9cDfIygA8tDwDNj+wmEmBiWMWayYxCi5PUO9NZaN41Uttwppf6b4D7iYZMdixBjIU2UQkxh6YmSLwMeJ2iSfCvwIYImXCHKmiQ4ISZRerC1lWO30Vp74xxCDDiTYBqyWmAbQXL79ji/rhDjTpoohZhESqmtQK6mwG1a68UTF40QlUVqcEJMrgsJJn/OJpHjeSFEHqQGJ4QQoiLJMAEhhBAVSRKcEEKIiiQJTgghREWSBCeEEKIiSYITQghRkf4/WdxK1Xos+aAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta = trace_rlg.posterior['θ'].mean(axis=0).mean(axis=0)\n",
    "idx = np.argsort(x_c)\n",
    "\n",
    "np.random.seed(123)\n",
    "\n",
    "\n",
    "plt.plot(x_c[idx], theta[idx], color='C2', lw=3)\n",
    "\n",
    "plt.vlines(trace_rlg.posterior['bd'].mean(), 0, 1, color='k')\n",
    "\n",
    "bd_hpd = az.hdi(trace_rlg.posterior['bd'])\n",
    "\n",
    "\n",
    "plt.fill_betweenx([0, 1], bd_hpd.bd[0].values, bd_hpd.bd[1].values, color='k', alpha=0.5)\n",
    "\n",
    "\n",
    "plt.scatter(x_c, np.random.normal(y_0, 0.02),\n",
    "            marker='.', color=[f'C{x}' for x in y_0])\n",
    "\n",
    "\n",
    "az.plot_hdi(x_c, trace_rlg.posterior['θ'], color='C2')  #green band \n",
    "\n",
    "\n",
    "plt.xlabel(x_n)\n",
    "plt.ylabel('θ', rotation=0)\n",
    "# use original scale for xticks\n",
    "locs, _ = plt.xticks()\n",
    "plt.xticks(locs, np.round(locs + x_0.mean(), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28aaf11",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "jupynb_env_new",
   "language": "python",
   "name": "jupynb_env_new"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}