10 Normal Two-Way ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "11ee4e0d-8a78-480f-9316-eb0dd4979b41",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import arviz as az\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") \n",
    "import scipy.stats as st\n",
    "from patsy import dmatrix, dmatrices\n",
    "from more_itertools import ncycles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "183bb97f-8107-437d-b165-67c2259f209d",
   "metadata": {},
   "source": [
    "#### 10.2 Data Generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cbc7ff85-ccaa-422d-8af4-92cf5e4413bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(73)\n",
    "# Choose sample size\n",
    "n_pop = 5\n",
    "n_elev = 3\n",
    "nsample = 12\n",
    "n = n_pop * nsample\n",
    "\n",
    "# Create factor levels\n",
    "pop = np.repeat(np.arange(1, n_pop+1), nsample)\n",
    "elev = list(ncycles(np.repeat(np.arange(1, n_elev+1), nsample / n_elev), int(n / nsample)))\n",
    "\n",
    "# Choose effects\n",
    "baseline = 40\n",
    "pop_effects = [-10, -5, 5, 10]\n",
    "elev_effects = [5, 10]\n",
    "interaction_effects = [-2, 3, 0, 4, 4, 0, 3, -2]\n",
    "all_effects = [baseline, *pop_effects, *elev_effects, *interaction_effects]\n",
    "\n",
    "sigma = 3\n",
    "eps = st.norm.rvs(0, sigma, size = n)\n",
    "\n",
    "data_ds = pd.DataFrame({'pop': pd.Series(pop, dtype = 'category'), \n",
    "                        'elev': pd.Series(elev, dtype = 'category')})\n",
    "X = dmatrix('~ pop * elev', data_ds, return_type = 'dataframe').to_numpy()\n",
    "wing = X @ all_effects + eps\n",
    "data = pd.DataFrame({'wing': wing, \n",
    "                     'pop': pd.Series(pop, dtype = 'category'), \n",
    "                     'elev': pd.Series(elev, dtype = 'category')})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01d5b6ef-1979-452e-8f1f-0f1d352df00a",
   "metadata": {},
   "source": [
    "**FIGURE 10.2**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9a40276e-72ef-41d7-87d6-80bec431651c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAGLCAYAAABk9f4xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1NklEQVR4nO3deVTWdf7//8cFsimypUml4TLWlI7LNLlMrhHieLIUxVSiESqzaLy0TLxItjBJMQvN1PbvoGVqjOmUWaJluTBNueRSnnHMAoXMQEBRQa7fH/Pr+mSxCML7zQX32zke4b29nu8LDz7O6/V6v94Wu91uFwAAAAznYnYBAAAAzRVBDAAAwCQEMQAAAJMQxAAAAExCEAMAADAJQQxAk1RYWKiSkhKzy7jE999/b2h7jfEzAHApghiAOjl69Kgefvhh3Xrrrerdu7fuuusurVmzxrF/2bJleuKJJ+q1zcWLF2vq1KmXdWxoaKiOHz9ep3b69u2r7Ozsej02KytL06dPr1M9dXUlnwEAY7QwuwAAzqeiokIPPPCAwsLC9Nxzz8nd3V3//ve/9eijj8rHx0ehoaGaMmWKqTUWFhaa2v6vnT59WhUVFYa22dg+AwC/RY8YgForKChQTk6O7rrrLnl6esrFxUV9+vTRE088obKyMkmX9l7NmjVL8+fP1/jx49WrVy/de++92rdvn8aPH6/evXsrOjraMYQWGRmpFStWONpasWKFIiMjf1PDuXPnlJSUpJCQEPXq1UvDhg3T5s2bJUlhYWGSpPDwcMe2N998U8OGDVPfvn0VExOjkydPOq61YcMGBQcH649//KPS0tKqvffqjj148KAmTZqkAQMGqGfPnoqOjtaPP/6offv2KTExUYcOHdJtt91W7bGVeeONNzRkyBD17dtXERER2r9/v2NfVfdV2WcAoPEhiAGotauuukp9+vRRVFSUFi1apF27duns2bMKDw/XnXfeWek5mZmZmjNnjrZv364ff/xRjzzyiJ5++ml9/PHHys3N1bvvvlurGl599VUdOXJEmZmZ+uKLLxQWFqaUlBRHW5K0Zs0a3XHHHdq4caNeeuklLVmyRNu2bVOHDh0cw4Rff/21Zs+erblz52rXrl2yWCxV9iTVdKzValVwcLA+/fRTffzxxyouLtaKFSvUo0cPJScn66abbtL27durPfbXvvvuO6Wnp2vlypXatWuX+vXrp9TUVEmq9r5+/RkAaJwIYgDq5JVXXtG9996rXbt26YEHHlCfPn302GOPqaCgoNLjhw4dqt/97ndq1aqV/vCHP2jw4MHq0qWLfH191bNnT+Xm5taq/YiICC1atEgtW7bUiRMn1KpVK+Xn51d67Nq1azVp0iR17dpVHh4eeuyxx7R3714dPXpUmzZt0sCBA9W3b1+5u7tr6tSpatmyZaXXqenYV199VRERESotLVV+fr78/f2rrOlyj/Xy8lJZWZlWr16tr7/+WjExMVq5cmWN9wXAOTBHDECdeHh4aNKkSZo0aZLOnz+vL774QmlpaYqLi9PSpUt/c7yvr6/ja1dXV/n4+Di+d3FxUW1fe1tcXKzk5GTt27dPHTp0UIcOHaq8xokTJ/T888/rhRdecGyzWCw6fvy4fvzxR7Vr186x3d3dXW3btq30OjUdu2/fPj344IM6c+aMbrzxRp0+fVoBAQGVXutyj23btq1efvllvfLKK3rjjTfk6+srq9WqMWPGVHtfnTp1quKTA9CYEMQA1Nr777+vBQsWKCsrSxaLRR4eHvrzn/+shx9+2DFs9msWi+Wyru3i4uKYZyZVPeE8MTFRXbp00bJly9SiRQt9/vnn2rhxY6XHtm3bVtHR0Ro7dqxj25EjR9ShQwd9+eWXOnDggGN7eXm5Tp06Vel1rr766iqPzcvLU2xsrN5880317NlTkmSz2SoNh7U59qefflLLli316quv6vz58/rggw8UGxurAQMGVHtfAJwDQ5MAaq1///46e/asnn76aZ06dUp2u13Hjh3TqlWrNHTo0Cu6dseOHbV582aVlJTo+++/1/r16ys9rqSkRJ6ennJ1ddWJEyeUnp4uSY4Q5+bm5ngAYPTo0Xr99dd17NgxVVRUKCMjQ+PGjVNpaalGjBihHTt2aOvWrSorK9OSJUuqXHurumPPnDkjSfL09JTdbtcnn3yiDz74wFGPu7u7zpw5I7vdXuOxv5Sbm6uoqCgdOHBAHh4e8vf3l4eHh1q2bFntff36MwDQONEjBqDW/P399eabb+r555/XnXfeqbNnz+qqq67SyJEjFRMTc0XXnjx5smw2mwYNGqSgoCCNGjVKO3fu/M1xNptNCQkJWrFihQICAjR+/HgdOHBAR44c0e9//3uFhYUpKipKSUlJGjVqlAoLC/Xggw/qxx9/VOfOnbV8+XL5+vrK19dXCxcu1DPPPKMffvhBI0aMUFBQUKW1denSpcpju3Tpoocfflh//etfdfHiRXXp0kXjx4/Xrl27JEm33nqr4+/t27dXe+wv/eEPf9Djjz+uv/3tb/rpp5907bXX6vnnn1fr1q119913V3lfki75DEaPHn1FPxcADcNir+3EDAAAANQLhiYBAABMQhADAAAwCUEMAADAJAQxAAAAkxDEAAAATEIQAwAAMAlBDAAAwCQEMQAAAJMQxAAAAEzilK846tu3r6677jqzywAAAKhRbm6usrOzK91nShBbvny5tmzZorKyMk2YMEE333yzpkyZoo4dO0qSJkyYoBEjRlR5/nXXXafMzEyDqgUAAKi7sLCwKvcZHsSys7O1e/duvfXWWyotLdVrr70mSYqKilJ0dLTR5QAAAJjG8CD22Wef6YYbblBMTIxKSko0c+ZMrV27VkePHlVWVpaCgoIUFxcnb29vo0sDAAAwlOGT9QsKCrR//36lp6crOTlZM2bMUI8ePTRz5kytXLlSHTp00JIlS4wuCwAAwHCGBzE/Pz8NGDBA7u7u6ty5szw8PDRkyBB1795dkhQSEqKDBw8aXRYAAIDhDA9it9xyiz799FPZ7Xbl5+ertLRUkydP1r59+yRJO3fuVLdu3YwuCwAAwHCGzxEbOnSoPv/8c40dO1Z2u10JCQkKCAhQSkqK3Nzc1KZNG6WkpBhdFgAAgOFMWb5i5syZv9m2atUqEyoBAAAwDyvrAwAAmIQgBgAAYBKCGAAAgEkIYgAAACYhiAEAAJjElKcmAQBwJqmpqcrLy6v1eYGBgbLZbA1QEZoKghgAADWoLkxZrValp6cbWA2aEoYmAQAATEIQAwAAMAlBDAAAwCQEMQAAAJMQxAAAAExCEAMAADAJQQwAAMAkBDEAAACTEMQAAABMQhADAAAwCUEMAADAJAQxAAAAkxDEAAAATEIQAwAAMAlBDAAAwCQEMQAAAJMQxAAAAExCEAMAADAJQQwAAMAkBDEAAACTEMQAAABMQhADAAAwCUEMAADAJAQxAAAAkxDEAAAATNLCjEaXL1+uLVu2qKysTBMmTFCfPn00a9YsWSwWde3aVYmJiXJxISMCAIArk5qaqry8vDqdGxgYKJvNVs8VXcrwIJadna3du3frrbfeUmlpqV577TWlpqZq2rRp6tu3rxISEpSVlaWQkBCjSwMAAA3ErEBU3XlWq1Xp6el1um59MTyIffbZZ7rhhhsUExOjkpISzZw5U6tXr1afPn0kSYMGDdL27dsJYgAANCE1BanGEIrMYHgQKygo0PHjx7Vs2TLl5OTo4Ycflt1ul8VikSS1atVKxcXFRpcFAACcVFJ8vAqKiup0rtVqrdN5/j4+SkpJqdO5v2R4EPPz81Pnzp3l7u6uzp07y8PD45KuyjNnzsjHx8fosgAAgJMqKCpSXO9AQ9ucu7tuw6y/ZngQu+WWW/T3v/9dUVFR+uGHH1RaWqr+/fsrOztbffv21bZt29SvXz+jywIANGNX0qMi1a1Xpb56VODcDA9iQ4cO1eeff66xY8fKbrcrISFB7du3V3x8vBYuXKjOnTsrNDTU6LIAAM2YM/eowLmZsnzFzJkzf7NtxYoVJlQCAABgHlOCGAAAaHoSExNVWFhY5/PrMsTr5+dX5/YaA4IYAACoF4WFhbJOnWpom+mLFhnaXn1j+XoAAACTEMQAAABMwtAkAABwes76FCpBDAAAOD1nXX6EoUkAAACTEMQAAABMwtAkAACoN86+nITRCGIAAKDesI5Y7TA0CQAAYBKCGAAAgEkIYgAAACYhiAEAAJiEIAYAAGASghgAAIBJWL4CAAA577sK4dwIYgAAyHnfVQjJ38fH8M/T38enXq5DEAMAAE4tKSWlTudZrValp6fXczW1wxwxAAAAkxDEAAAATMLQJAAAqBd+fn6Gv/vRz8/P0PbqG0EMAADUi+Tk5Dqf2xjma5mBoUkAAACTEMQAAABMQhADAAAwCUEMAADAJAQxAAAAk/DUJACg2XPmV+TAuRHEAADNXl1fkSM132UXUD8YmgQAADAJQQwAAMAkpgxNjho1Sq1bt5YktW/fXvfee6+mTJmijh07SpImTJigESNGmFEaAACAYQwPYufPn5ckZWRkOLatWbNGUVFRio6ONrocAAAA0xgexL7++muVlpYqOjpa5eXleuyxx7R//34dPXpUWVlZCgoKUlxcnLy9vY0uDQAAwFCGBzFPT0/df//9Cg8P17fffqsHH3xQkydPVnh4uLp3766lS5dqyZIlio2NNbo0AADQxKSmpiovr+qlSaxWa5X7AgMDZbPZGqIsB8ODWKdOnRQUFCSLxaJOnTrJz89PAwcO1DXXXCNJCgkJUcoVPEaM6tX0D7IqRvxjBIDGqq7/mfO703yN/fM3PIitXbtWhw8fVlJSkvLz81VSUqKYmBglJSWpR48e2rlzp7p162Z0Wc1Gdf8gWQsHACrX2P8zh/MyPIiNHTtWNptNEyZMkMVi0dy5c+Xh4aGUlBS5ubmpTZs29IgBAIBmwfAg5u7urmefffY321etWmV0KQAAAKZiQVcAAACTEMQAAABMQhADAAAwiSmvOAIAAM3L5Syf1ByXASGIAQCABtdUg9SVYmgSAADAJPSINTFJ8fEqKCqq8/nVveqhOv4+Pkpi/TcAcGp2u135+flq166dLBaL2eU0CwSxJqagqEhxvQMNb3fu7tq/NgkA0Ljs2bNHb7zxhiZNmqTevXubXU6zwNAkAACQJGVlZV3yNxoeQQwAAEiSgoODL/kbDY+hSQAAIEnq1auXbDab2rVrZ3YpzQZBDAAASJIsFosCA42fZ9ycMTQJAABgEoIYAACASQhiAAAAJiGIAQAAmITJ+gAANDHxCUkqOl1geLs+vv5KeSrJ8HadGUEMAIAmpuh0gboNe9rwdg98+KThbTo7hiYBAABMQo8YAABNEL1TzoEgBgBAE8TQpHMgiAEAnEZqaqry8vJqfV5gYKBsNlsDVARcGYIYAMBpVBemrFar0tPTDaym8fLx9Teld8rH19/wNp0dQQwAgCbmSpaQINAai6cmAQAATEIQAwAAMAlBDAAAwCTMEQMAoJmp6elTq9Va6XaePq1/BDEAQKOSmJiowsLCOp1bVYCoiZ+fn5KTk+t0rjMiTDUeBDEAQKNSWFgo69SphraZvmiRoe0BP2OOGAAAgEnoEWuC5u6u/arTAADAeKYEsVGjRql169aSpPbt22vKlCmaNWuWLBaLunbtqsTERLm40FlXV3G9Aw1vk/AHAEDtGR7Ezp8/L0nKyMhwbJsyZYqmTZumvn37KiEhQVlZWQoJCTG6NAAAAEMZ3u309ddfq7S0VNHR0brvvvu0Z88eHThwQH369JEkDRo0SDt27DC6LAAAAMMZ3iPm6emp+++/X+Hh4fr222/14IMPym63y2KxSJJatWql4uJio8sCAAAwnOFBrFOnTgoKCpLFYlGnTp3k5+enAwcOOPafOXNGPj4+RpcFAABgOMOHJteuXatnnnlGkpSfn6+SkhLddtttys7OliRt27ZNf/rTn4wuCwAAwHCG94iNHTtWNptNEyZMkMVi0dy5c+Xv76/4+HgtXLhQnTt3VmhoqNFlAQAAGM7wIObu7q5nn332N9tXrFhhdCkAAACmYrEuAAAAk7CyPgCg0eHdj2guCGIAgEaHl36juWBoEgAAwCQ19oidO3dOH374oX766SfZ7XbH9qioqAYtDAAAoKmrMYg9/vjjOnHihG644QbH6vcAAAC4cjUGscOHD2vTpk1ycWEUEwAAoD7VmK6uuuoqlZeXG1ELAABAs1Jlj9jrr78uSWrbtq0iIyMVHBwsNzc3x37miAEAAFyZKoPY4cOHJUne3t7y9vbW0aNHDSsKAJqa1NRU5eXl1encwMBA2Wy2eq4IQGNQZRBLTU2VJG3evFl33HHHJfvWrVvXoEUBQFNTXZCyWq1KT083sBoAjUWVQWzLli0qLy/X/PnzZbfbHUtXlJeXa/HixRo1apRRNQJAvaFnCkBjUmUQO3TokHbt2qVTp07p73//+/+d0KKFJk2aZERtAFDv6JmqHwRaoH5UGcRiYmIUExOjlStXKiIiwsiaAACNwPTp01VRUVHv183Ly5PVaq10H0slobmpcR2xgoICvfDCC47vLRaLvLy81LVrVw0cOLBBiwMAmKeiokLdhj1taJsHPnzS0PYAs13Wgq67d+9WaGioXF1d9dFHH+m6667Txo0btW/fPsXExBhRJwAAQJNTYxA7deqUMjMz1bZtW0nSlClTZLVatXLlSo0ZM8bpgxjzHAAAgFlqDGKFhYWOECZJ/v7+KiwslLu7u1q0qPH0Rq+mIOVsk3f9fXw0d3fdguWVtgsA9cHPz0/pixYZ3iZghhqTVIcOHfTss89q3LhxkqS1a9fq+uuv1969e5lU2QglpaTU+VxnC51AVRITE1VYWFinc6uaRF4TPz8/JScn1+lcXKqunyO/w+CMagxic+fO1Zw5czR69Gi5urpq6NChmjNnjt5//33FxsYaUSMA1EphYaGsU6ca2qbRPTgAmoYag1hAQIAWLlz4m+0TJ05skIIAAACaixqD2O7du7Vw4UKdPn3asbq+JG3YsKFBCwMAmMvH19/w5SR8fP0NbQ8wW41BLCEhQWFhYbr55ptlsViMqKnexSckqeh0QZ3Pr8ucER9ff6U8lVTnNgHAbPwOAxpejUGsRYsWioqKMqKWBlN0uoBFCQEAQKNTYxDr2rWrvvnmG914441G1NNssH4ZAACoMYh9//33GjNmjK699lp5eHg4tjNH7Mrw4mEAAFBjEJs+fboRdQAAADQ7Na7I2qdPH3l6euq///2vevXqJTc3N/Xp08eI2gAAAJq0GoNYZmambDabXnnlFRUXF+uRRx7R6tWrjagNAACgSasxiGVkZOjtt9+Wt7e3rrrqKmVmZur//b//Z0RtAAAATVqNQczFxUXe3t6O76+55hq5uro2aFEAAADNQY1BzM/PT4cOHXIs5rp+/Xr5+vo2eGEAAABNXY1PTcbFxclqteq7777TgAED5OHhoRdffPGKGj116pTCwsL02muv6dy5c5oyZYo6duwoSZowYYJGjBhxRdcHADMkxceroKioTufW5Q0ekuTv46OklJQ6neuMalqDsarPkfUX0VjVGMS6dOmid999V99++60uXryoTp06yc3Nrc4NlpWVKSEhQZ6enpKkgwcPKioqStHR0XW+ZmOVmJiowsLCOp1b11/Kfn5+Sk5OrnI/v8SAhlNQVKS43oGGtjl3d90WhnZW/B5CU1NlEHv99dcr3b59+3ZJqvNrj+bNm6fx48frpZdekiTt379fR48eVVZWloKCghQXF3fJnDRnVlhYKOvUqYa2mb5oUbX7+SUGAEDjUWUQO3z4cL03lpmZqYCAAA0cONARxHr06KHw8HB1795dS5cu1ZIlSxQbG1vvbQMAADQ2VQax1NTUem/snXfekcVi0c6dO3Xo0CHFxsZq6dKlatu2rSQpJCREKc1orgMAAGjeanxqsj6tXLlSK1asUEZGhm666SbNmzdPjzzyiPbt2ydJ2rlzp7p162ZkSQAAAKapcbJ+Q0tKSlJKSorc3NzUpk0besQAAECzYVoQy8jIcHy9atUqs8oA0ETV9OAKADQGNQaxF1544ZLvLRaLvLy81LVrVw0cOLDBCgOAK9HYnlgGgMrUGMQOHz6s3bt3KzQ0VK6urvroo4903XXXaePGjdq3b59iYmKMqBMAAKDJqTGInTp1SpmZmY4nG6dMmSKr1aqVK1dqzJgxThPEDnz4pNklAAAAXKLGIFZYWOgIYZLk7++vwsJCubu7q0UL0+f6X7Zuw542tD2CH+riSl6RU1fN7RU5ANCY1JikOnTooGeffVbjxo2TJK1du1bXX3+99u7dKxcXQ1e/AJo8XpEDZ2e325Wfn6927drJYrGYXQ7Q6NUYxObOnas5c+Zo9OjRcnV11dChQzVnzhy9//77rIAPAL/S3IPtnj179MYbb2jSpEnq3bu32eUAjV6NQSwgIEALFy78zfaJEyc2SEEA4Myae49mVlaW42+CGFCzGoPY7t27tXDhQp0+fVp2u92xfcOGDQ1aGADA+QQHB+uNN95QcHCw2aUATqHGIJaQkKCwsDDdfPPNjPcDAKrVq1cv2Ww2tWvXzuxSAKdQYxBr0aKFoqKijKgFAODkLBaLAgONHZ4FnFmNjz127dpV33zzjRG1AAAANCs19oh9//33GjNmjK699lp5eHg4tjNHDGgYjW3yNQCg4dQYxKZPn25EHQD+f839qTsAaE6qDGJHjhxRly5d1KpVKyPraXJ4ETAAAKhKlUFs3rx5eumll/S3v/3tN/ssFotjrRhUzzp1qqHtEfwAAHAeVQaxnxdx3bJli2HFAAAANCdVBrH+/fvrlltu0ZAhQzRkyBB17NjRwLIAAACaviqD2CeffKJdu3Zp586dWrFihVxcXDR48GANGTJEffr0kZubm5F1AgAANDlVBrGAgACNGDFCI0aMkCTl5uZqx44dSktL07Fjx7R7927DigSA2vDz8zN8vqSfn58KCwsNbROA87usdcS2bNmi7du36+DBg+rWrZvGjRtnRG1As+Pv42P4chL+Pj6GtmeE5OTkKvelpqYqL69un3FgYKBsNluV+61Wa52uC6D5qjKIPffcc9qyZYvOnDmjgQMHauLEierfv/8li7oCqF9JKSl1PtdqtSo9Pb0eq2maqgtSAGC0KoPY8uXLdfvtt2vy5Mnq1auXgSUBAAA0D1UGsQ8++EBbt27Vs88+q2+//Va33XabhgwZogEDBsjb29vIGgHo8obUqhoaq2lIDQBgjiqDWMeOHRUVFaWoqCgVFRVp27Zt+uijj5SamqrOnTvr9ddfN7LOK+Lj668DHz5peJtFpwsMbRNNG0EKAJqeGifrS9Lx48f1008/6cKFC3Jzc5OLi0tD11WvUp5KqvO5VzLvhom7AACgOlUGsYyMDGVnZ+vzzz+Xr6+vBg0apLFjx6pv377y9PQ0skYAAIAmqdoFXQcNGqQZM2awqj4AXAaWHwFQW1UGsVdeecXIOgDA6dV1+RGWHgGar8uaI9aU8SQaAAAwS7MPYgQpAABgFud6/BEAAKAJafY9Yg3JrBcPAwAA52BKEDt16pTCwsL02muvqUWLFpo1a5YsFou6du2qxMREp1unrCrVvXgYAADA8MRTVlamhIQEx1pkqampmjZtmt58803Z7XZlZWUZXRIAAIApDA9i8+bN0/jx43X11VdLkg4cOKA+ffpIkgYNGqQdO3YYXRIAAIApDB2azMzMVEBAgAYOHKiXXnpJkmS322WxWCRJrVq1UnFxcY3XKS4uvmTY78EHH5Qkvfzyy45tgwcP1pAhQ/Tss8+qpKREknTNNddo8uTJ2rBhg7788kvHsY899piOHz+uVatWObbdeeeduuWWWy5p54YbbtCECRP01ltv6fDhw47tiYmJ+uKLL/TPf/7TsW38+PG69tprtXDhQse2P/7xjxo5cqReeuklnThxQpLk7e2txx9/XB9//LE++eQT7ol74p6a4T0FBAQoOTm5Sd1TU/w5cU/c05XeU2UsdrvdXu0R9SgiIkIWi0UWi0WHDh1Sx44ddfDgQR08eFCStHnzZu3YsUMJCQnVXicsLEyZmZlGlAwADY4FXYGmrbrcYmiP2MqVKx1fR0ZGKikpSWlpacrOzlbfvn21bds29evXz8iSAAAATGP644mxsbFavHix7rnnHpWVlSk0NNTskgAAAAxh2jpiGRkZjq9XrFhhVhkAAACmMb1HDAAAoLkiiAEAAJiEIAYAAGASghgAAIBJCGIAAAAmIYgBAACYhCAGAABgEoIYAACASQhiAAAAJjFtZX0AzVtqaqry8vLqdG5gYKBsNls9VwQAxiOIATBFTUHKarUqPT3doGoAwBwMTQIAAJiEIAYAAGASghgAAIBJCGIAAAAmIYgBAACYhCAGAABgEoIYAACASQhiAAAAJmFBVwANJj4hSUWnC+p8vtVqrfU5Pr7+Snkqqc5tAoCRCGIAGkzR6QJ1G/a0oW0e+PBJQ9sDgCvB0CQAAIBJCGIAAAAmIYgBAACYhCAGAABgEoIYAACASQhiAAAAJmH5CgANiuUkAKBqBDEADYp1xACgagxNAgAAmIQgBgAAYBKCGAAAgEkMnyN28eJFzZ49W0ePHpWrq6tSU1NVXFysKVOmqGPHjpKkCRMmaMSIEUaXBgAAYCjDg9jWrVslSatWrVJ2drZSU1N1++23KyoqStHR0UaXAwAAYBrDg9gdd9yhIUOGSJKOHz+uNm3aaP/+/Tp69KiysrIUFBSkuLg4eXt7G10agHrm4+tv+FOMPr7+hrZ3uVJTU5WXl1flfqvVWuW+wMBA2Wy2higLgMksdrvdbkbDsbGx+uijj7Ro0SLl5+frxhtvVPfu3bV06VIVFRUpNja2ynPDwsKUmZlpYLUAjGa1WpWenm52GQBwxarLLaZN1p83b542bdqk+Ph4DRgwQN27d5ckhYSE6ODBg2aVBQAAYBjDhybXrVun/Px8PfTQQ/Ly8pLFYtGjjz6q+Ph49ejRQzt37lS3bt2MLguolZqGmarCEBMA4JcMD2LDhg2TzWZTRESEysvLFRcXp2uuuUYpKSlyc3NTmzZtlJKSYnRZQK1UF6YYUgMAXC7Dg1jLli0r/U9q1apVRpcCAABgKhZ0BQAAMAlBDAAAwCQEMQAAAJMYPkcMAKTLe/K0qkVOefoUQFNBEANgCoIUADA0CQAAYBp6xIBKJCYmqrCwsM7nV/fewKr4+fkpOTm5zm0CAJwPQQyoRGFhoaxTpxraZvqiRYa2BwAwH0OTAAAAJiGIAQAAmIQgBgAAYBKCGAAAgEkIYgAAACbhqUkAANAolJWVKScnR+fOnTO7lDrx9PRU+/bt5ebmdtnnEMSAKrCcBAAYKycnR61bt1bHjh1lsVjMLqdW7Ha7Tp06pZycHHXq1OmyzyOIAVVgHTEAMNa5c+ecMoRJksVi0VVXXaWTJ0/W6jzmiAEAgEbDGUPYz+pSOz1iAACg0crOzta0adP0u9/9zrHN399fBQUFSkpKUpcuXa7o+ufPn9f69esVHh6uzMxM+fr6Kjg4+ErLvmwEMQAA0Kj169dPzz333CXbIiMj6+XaJ0+e1Jo1axQeHq6wsLB6uWZtEMQAAIDTKi4u1pNPPqmCggJJ0uzZs5WTk6PNmzcrNTVVkjRq1Ci9+uqr2rhxoz788EOVl5erdevWWrx4sZYtW6b//Oc/euGFF2S329WmTRtNmDBBzzzzjL744gtJ0p133qm//vWvmjVrltzd3ZWbm6sffvhBzzzzjLp163ZF9RPEgEr4+fkZPnnez8/P0PYAwFns2rXrkh6wwYMHO75etmyZ+vXrp4kTJ+rbb7+VzWbTihUrlJaWprNnz+o///mPrr/+evn7+6uwsFBvvPGGXFxcdP/99+urr77SlClTdPjwYT366KNavHixJGnr1q3KycnR6tWrVV5erokTJ6pfv36SpGuvvVZPPfWUVq9erbfffltPPfXUFd0bQQyoRHJycp3PtVqtSk9Pr8dqAKB5q2xo8pNPPpEkHT58WLt27dLGjRslSUVFRXJ1dVVoaKg+/PBD7dmzR+Hh4XJxcZGbm5see+wxtWzZUnl5eSovL6+0vSNHjuhPf/qTLBaL3Nzc1LNnTx05ckSSdNNNN0mSAgMD9eWXX17xvfHUJAAAcFqdO3fWpEmTlJGRoeeff14jR46UJI0dO1br16/X3r17ddttt+nrr7/W5s2b9fzzzys+Pl4VFRWy2+1ycXFRRUXFJdfs0qWLY1iyrKxMu3fvVlBQkKT6f6qTHjGgDlJTU5WXl1flfqvVWun2wMBA2Wy2hioLAJqkXw9NSnKsvj9lyhQ9+eSTWr16tUpKSvToo49Kkjp06CBJCg4OlouLi4KCguTl5aWwsDC5u7urbdu2+uGHH9S7d2+VlZUpLS1Nnp6ekqShQ4fqX//6l+655x6VlZVp+PDhVzwXrCoWu91ub5ArN6CwsDBlZmaaXQYAAKhHhw4dcgz9OavK7qG63MLQJAAAgEkIYgAAACZhjhgavfiEJBWdLjC0TR9ff6U8lWRomwCA5ocghkav6HSBug172tA2D3z4pKHtAQCaJ4YmAQAATEIQAwAAMAlDkwAAwClkvrtJBcUX6u16/q3dFXZ3aI3H7d27VwsWLFBGRka9tf0zw4PYxYsXNXv2bB09elSurq5KTU2V3W7XrFmzZLFY1LVrVyUmJsrFhc46AADwfwqKL+ii/9D6u17B1hqPefnll7V+/Xp5eXnVW7u/ZHgQ27r1fze9atUqZWdnO4LYtGnT1LdvXyUkJCgrK0shISFGl4ZGjMnzAAAzXH/99Vq8eLFmzpzZINc3PIjdcccdGjJkiCTp+PHjatOmjT7++GP16dNHkjRo0CBt376dIIZL8NQkAMAMoaGhysnJabDrmzL+16JFC8XGxiolJUWhoaGy2+2Ol2i2atVKxcXFZpQFAABgKNMmYs2bN0+bNm1SfHy8zp8/79h+5swZ+fj4mFUWAACAYQwPYuvWrdPy5cslSV5eXrJYLOrevbuys7MlSdu2bdOf/vQno8sCAAAwnOFzxIYNGyabzaaIiAiVl5crLi5OXbp0UXx8vBYuXKjOnTsrNLTmR0kBAEDz4t/a/bKedKzN9S5H+/bttXr16npr95cMD2ItW7ZUenr6b7avWLHC6FIAAIATuZw1v5wNi3UBAACYhCAGAABgEoIYAACASQhiAAAAJiGIAQAAmMTwpyaB2vLx9Tf8lUM+vv6GtgcAqNnGjRtVWlpab9fz8vLSX/7ylyr3l5WVKS4uTrm5ubpw4YIefvhhBQcH11v7EkEMTiDlqaQ6nWe1WitdKgUA4JxKS0vVq2fPervenr17q92/fv16+fn5KS0tTQUFBRo9ejRBDAAAwAjDhw+/ZJF5V1fXem+DIAYAAFCJVq1aSZJKSko0depUTZs2rd7bYLI+AABAFU6cOKH77rtPd999t0aOHFnv16dHDAAAoBI//vijoqOjlZCQoP79+zdIG/SIAQAAVGLZsmUqKirSiy++qMjISEVGRurcuXP12gY9YgAAwCl4eXnV+KRjba9XndmzZ2v27Nn11l5lCGIAAMApVLfml7NiaBIAAMAkBDEAAACTEMQAAABMQhADAAAwCUEMAADAJDw1CQAAnMIHmWt0vvBUvV3Pw+8qDQ8Lr3L/xYsXNXv2bB09elSurq5KTU3V9ddfX2/tSwQxAADgJM4XnlI/l4J6u96uwur3b926VZK0atUqZWdnKzU1VUuXLq239iWCGAAAQKXuuOMODRkyRJJ0/PhxtWnTpt7bIIjBqaWmpiovL6/K/VartdLtgYGBstlsDVUWAKCJaNGihWJjY/XRRx9p0aJF9X/9er8iYCDCFACgoc2bN08zZszQuHHj9N5776lly5b1dm2emgQAAKjEunXrtHz5ckn/ey+lxWKRq6trvbZBjxgAAEAlhg0bJpvNpoiICJWXlysuLk4eHh712gZBDAAAOAUPv6tqfNKxtterTsuWLZWenl5/DVaCIAYAAJxCdWt+OSvmiAEAAJiEIAYAAGASghgAAGg07Ha72SXUWV1qJ4gBAIBGwdPTU6dOnXLKMGa323Xq1Cl5enrW6jzDJ+uXlZUpLi5Oubm5unDhgh5++GEFBgZqypQp6tixoyRpwoQJGjFihNGlAQAAE7Vv3145OTk6efKk2aXUiaenp9q3b1+rcwwPYuvXr5efn5/S0tJUUFCg0aNHKyYmRlFRUYqOjja6HAAA0Ei4ubmpU6dOZpdhKMOD2PDhwxUaGur43tXVVfv379fRo0eVlZWloKAgxcXFydvb2+jSAAAADGX4HLFWrVrJ29tbJSUlmjp1qqZNm6YePXpo5syZWrlypTp06KAlS5YYXRYAAIDhTJmsf+LECd133326++67NXLkSIWEhKh79+6SpJCQEB08eNCMsgAAAAxl+NDkjz/+qOjoaCUkJKh///6SpPvvv1/x8fHq0aOHdu7cqW7dulV7jdzcXIWFhRlRLgAAwBXJzc2tcp/FbvAzonPmzNHGjRvVuXNnx7Zp06YpLS1Nbm5uatOmjVJSUpgjBgAAmjzDgxgAAAD+hwVdAQAATEIQAwAAMAlBDAAAwCQEMQAAAJMQxBqhvXv3KjIy0uwyUAtlZWV64oknNHHiRI0dO1ZZWVlml4RauHjxomw2m8aPH6+IiAh99913ZpeEWjp16pQGDx6sI0eOmF0KamnUqFGKjIxUZGSkbDab2eUYzvB1xFC9l19+WevXr5eXl5fZpaAWKnuHanBwsNll4TJt3bpVkrRq1SplZ2crNTVVS5cuNbkqXK6ysjIlJCTI09PT7FJQS+fPn5ckZWRkmFyJeegRa2Suv/56LV682OwyUEvDhw+X1Wp1fO/q6mpiNaitO+64QykpKZKk48ePq02bNiZXhNqYN2+exo8fr6uvvtrsUlBLX3/9tUpLSxUdHa377rtPe/bsMbskwxHEGpnQ0FC1aEFHpbOp7B2qcC4tWrRQbGysUlJSFBoaanY5uEyZmZkKCAjQwIEDzS4FdeDp6an7779fr776qpKTkzVjxgyVl5ebXZahCGJAPfn1O1ThfObNm6dNmzYpPj5eZ8+eNbscXIZ33nlHO3bsUGRkpA4dOqTY2FidPHnS7LJwmTp16qS77rpLFotFnTp1kp+fX7P7+dH1AtSDyt6hCuexbt065efn66GHHpKXl5csFgvDy05i5cqVjq8jIyOVlJSktm3bmlgRamPt2rU6fPiwkpKSlJ+fr5KSkmb386NHDKgHy5YtU1FRkV588UXH0z/nzp0zuyxcpmHDhungwYOKiIjQ/fffr7i4OHl4eJhdFtDkjR07VsXFxZowYYKmT5+uuXPnNrvpObxrEgAAwCT0iAEAAJiEIAYAAGASghgAAIBJCGIAAAAmIYgBAACYhCAGAABgEoIYAACASQhiAAAAJiGIAQAAmIQgBgAAYBKCGAAAgEkIYgAqlZOTo5tuukl33323489dd92ltWvXNkh7mZmZeuihh2o8bvbs2dq/f78k6cknn9SOHTvqpf3Fixfrqaeeqpdr/VpkZKQ++OCDer/uxx9/rPT0dElSVlaW5syZU+9tAGhYzesV5wBqxdPTU++++67j+/z8fN15553q3r27fv/735tS044dO3TPPfdIkp5++mlTamgsvvrqK50+fVqSFBwcrODgYJMrAlBbBDEAl61du3YKCgrSt99+q9///vdasmSJ3nvvPbm6uqpTp06Kj49X27ZtFRkZqZtvvllffPGFCgoKdPfdd2vq1KnKycnRyJEjtXv3bkn6zfc/27Nnj9LS0nThwgWdPHlSf/7znzV37lw999xz+uGHHzRjxgzNnz9fCxYsUEREhIYPH67NmzfrhRdeUEVFhVq1aiWbzaYePXpo8eLFys3N1cmTJ5Wbm6t27dopLS1NV1999W/u78iRI4qIiNDp06d10003KTExUd98840ef/xxbdmyRS4uLiotLdXtt9+u9957TwEBAZecv3TpUn344YeqqKjQddddp8TERLVr1+6SY7788kstWLBApaWlcnFx0aOPPqqhQ4dq/PjxioqKUmhoqCQpLS1NkhQTE6OkpCQdO3ZMhYWFatWqlRYsWKDi4mKtWrVKFy9eVOvWrRUUFKRNmzZp+fLlysvLU1JSknJzc2W32zVq1Cg98MADysnJ0aRJkzR48GDt3btXRUVFeuKJJxQSElJv/0YA1A5DkwAu2+7du/Xdd9+pZ8+eeuedd/Tpp59q7dq12rBhg7p27apZs2Y5jj169Kjeeust/eMf/9D777+vrVu3XnY7f//73zV16lStWbNG7733nrZs2aL9+/dr+vTpuvrqq7VgwQL17NnTcfyRI0eUmJioxYsXa/369Zo6daoeeeQRlZSUSJL+/e9/Kz09XR988IG8vLy0atWqStv97rvvtHjxYm3YsEF2u11Lly7VLbfcIl9fX3366aeSpPfee0/9+/f/TQhbt26dDh8+rDVr1ujdd9/V4MGDNXv27EuOOX36tGw2m+bPn69//OMfevHFF5WUlKTjx48rPDxcmZmZkqSLFy9q/fr1Cg8P17Zt2+Tj46O3335bmzZtUvfu3bVy5Ur17NlT48eP14gRIzR9+vRL2pkxY4b69u2rDRs26K233tL69ev13nvvSZK+//57DRgwQGvXrtXjjz+uuXPnXvbPBUD9o0cMQJXOnTunu+++W9L/woG/v7/S0tJ0zTXXaNu2bQoLC1PLli0lSffdd5+WLVumCxcuSJLuueceubm5yc3NTcOHD9dnn32mrl27Xla7zzzzjLZt26Zly5bpv//9r86fP6+zZ89WefyuXbvUr18/dejQQZIcQennuWR9+vSRt7e3JOnmm292DOf9WkhIiCNgjRkzRvPnz5ckRUREaPXq1Ro8eLDefvttzZw58zfnbt26VV999ZXGjBkjSaqoqFBpaeklx+zZs0cnT55UTEyMY5vFYtE333yjESNGaP78+Tp58qQOHjyojh07Ov506NBBGRkZOnbsmP71r3+pd+/eVX4WZ8+e1ZdffqnXXntNktS6dWuFhYVp27Zt6tmzp9zc3DR48GDHZ1FYWFjltQA0PIIYgCr9eo7YL1VUVMhisVzyfXl5ueP7Fi3+79eL3W6Xi4uLLBaL7Ha7Y3tZWVml17733nt14403auDAgfrLX/6ivXv3XnJeTbX83ObP9Xh6ejq2/7qGX3J1db3kmj/fw8iRI7Vw4ULt2rVLZ8+e1a233lppDQ888IAmTpwoSbpw4cJvAt/FixfVpUsXrVmzxrEtPz9fAQEBcnNzU2hoqP75z39q9+7dCg8PlyS9+eabWr16tSIiIjRy5Ej5+fkpJyen2s/i1/f3y5+Nm5ubXFxcHJ8FAHMxNAmgTgYOHKh33nnH0VOVkZGhW2+9Ve7u7pKk9evXq6KiQqdPn9bGjRt1++23y8fHR2VlZfrPf/4jSY7hsl8qKirSV199pRkzZmjYsGHKy8vTd999p4qKCkn/C0u/DHzS/3rAPvvsM33//feSpJ07d+rEiROXDF9eji1btuj06dO6ePGiVq9erUGDBkmSvLy8dNdddykuLk7jx4+v9Nyfh/t+Hg5NT0//Tc9Zr169dOzYMX3++eeSpEOHDik0NFT5+fmSpHHjxukf//iHvvzyS8dcsc8++0yjR49WeHi4OnXqpC1btujixYtVfhbe3t7q2bOnVq5cKUkqLi7WunXr9Oc//7lWnwUAY9AjBqBOxo4dqxMnTig8PFwVFRUKCgrSggULHPvPnTunsWPH6syZM5o4caL69+8vSXriiSf04IMPKiAgQMOHD//NdX18fDR58mSNHj1aLVu2VLt27fTHP/5Rx44dU//+/RUSEqInnnhCSUlJjnN+97vfKTExUY8++qguXrwoT09PLVu2TK1bt67VPXXp0kUPPfSQioqKdMstt2jy5MmOfWFhYVq9erVGjRpV6bnh4eHKz8/XuHHjZLFYdM011+iZZ5655JiAgAAtWrRI8+fP1/nz52W32zV//ny1b99ektS9e3e5urpq+PDh8vDwkCRFR0crISHBsWxIr169dPjwYUlSv379NGPGDKWkpKhbt26OdhYsWKCnnnpKmZmZunDhgkaOHKmwsDDl5ubW6vMA0PAs9ur6+wGgDiIjIx1PMzYFdrtdL7/8snJzc5WcnGx2OQCaEHrEAKAGwcHBuvrqq/Xiiy+aXQqAJoYeMQAAAJMwWR8AAMAkBDEAAACTEMQAAABMQhADAAAwCUEMAADAJAQxAAAAk/x/VymAE47FWAwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_style(style = 'white')\n",
    "plt.figure(figsize = (10,6))\n",
    "g = sns.boxplot(x = pop, y = wing, hue = elev, palette = 'coolwarm', fliersize = 2, linewidth = 1, saturation = 0.8)\n",
    "plt.title('Simulated data set', pad = 10, fontsize = 13)\n",
    "plt.xlabel('Population by elevation', labelpad = 10, fontsize = 12)\n",
    "plt.ylabel('Wing length', labelpad = 10, fontsize = 12)\n",
    "g.legend(title = 'Elevation', loc = 'lower right')\n",
    "plt.axhline(40, 0, 5, lw = 1, ls = '--', color = 'gray');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8ea6a4c5-f178-4c0c-87bf-54e70d841b94",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABE8AAAILCAYAAAAZu7JXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABWDklEQVR4nO3de1jUdcL+8XvkqCCRYNp6QEmlEjFK0La0EhWyYrWsJMASq6W0siJQDA9NeVitVlFbO9g+CeWqWblPT1hNbdZmPLQ/VkXUDTNdU1xDXIQAOczvjx5mJXUEAr4zw/t1XV4L3xmYm7naD8w9n4PJarVaBQAAAAAAgHPqZHQAAAAAAAAAR0Z5AgAAAAAAYAflCQAAAAAAgB2UJwAAAAAAAHZQngAAAAAAANhBeQIAAAAAAGAH5Qk6nNGjRyskJMT2b/DgwYqKitLLL7/cLo+fmZmpuLi4Jt339OnTWr9+ve3zxMREvfjii20V7bwZbrnlFmVmZrbr4wJwDYy5TZORkdHoeQoJCdEf//jHdnlsAABwYe5GBwCMMGvWLN16662SpNraWn311VeaM2eOLrnkEk2YMMHYcGd4//33tXr1ak2ePFnSTy8CPDw82jXDSy+9pKKionZ9TACuhTH3woqKipSamqrY2FjbNV9f33Z5bAAAcGHMPEGH5Ovrq+7du6t79+669NJLNXHiRF177bX68MMPjY7WiNVqbfS5v7+/fHx82u3x9+7dq40bNyo4OLjdHhOA62HMvbBvv/1WoaGhtuepe/fu6ty5c7s8NgAAuDDKE+D/uLu7y9PTU5JUX1+vV199VWPGjFFYWJgSEhK0d+9e231DQkK0YcMGjR07VuHh4XriiSdUXl4uScrNzVVISIhqa2tt9581a5ZSUlLO+bhvv/22br75ZoWGhmr48OGaN2+eamtrlZubq9mzZ+vYsWMKCQnR4cOHz5pCvnnzZo0fP15hYWG6/fbblZuba7tt9OjRysrK0uTJkzV06FBNmjRJO3futN2emJioWbNmnff5qKurU3p6ulJSUuTv79+8JxMALoAx9z+OHz+ukydPqn///i14JgEAQHugPEGHV1NTow8//FB//etfNXr0aEnSqlWrtHbtWs2ePVvvvPOOevfurfvvv9/2x7okrVixQunp6XrjjTf0zTff6Omnn272Y3/99ddasGCBHn/8cW3dulULFizQ5s2b9eGHHyo8PFzp6enq3r27vvjiC1166aWNvnbz5s165pln9OCDD+q9997TddddpwcffFBHjhyx3WflypW6//779ac//UleXl4ym8222zIzMzVnzpzzZnvttdd08cUXO9SUegDOjzH3bEVFRXJ3d9fy5cs1cuRIxcbGavPmzc3++QAAQNuhPEGH9Mwzzyg8PFzh4eEKCwtTWlqa7r33XsXGxspqtSorK0szZsxQVFSULrvsMpnNZrm7u+u9996zfY/7779fN910k4YMGaI5c+Zo69atOnnyZLNyeHt767nnntO4cePUq1cvxcTE6Morr1RRUZE8PT3VtWtXderUSd27d5ebm1ujr123bp3i4+M1YcIE9e/fX08++aQuv/xyrVu3znafCRMmaMyYMbr88suVlJSkgoIC223+/v7q2rXrOXMdOHBAr732mhYsWNCsnwcAzoUx1/6Y++2330qSLr/8cr3yyiuaNGmS5s6dqw8++KBZPx8AAGg7bBiLDmnGjBmKiYmRJHl5eTX6Q7mkpEQnT57U0KFDbff38PBQaGio9u/fb7sWHh5u+zg0NFT19fU6cOBAs3KEhobK29tbK1asUFFRkfbt26eDBw9qxIgRF/za/fv366GHHmp07aqrrrL9ES5Jffr0sX3s6+ur+vp61dXVnfWi4ExWq1Vz5szRQw89pN69ezfr5wGAc2HMPf+YK0n33HOPbrnlFtsSycsvv1wHDx7UW2+9pZtvvrkZPyEAAGgrzDxBh9StWzcFBQUpKChIPXv2bPSHrZeX1zm/pq6uTnV1dbbPz/ya+vp6SVKnTp1kMpnO+toz1+Kf6fPPP9fEiRN1/PhxjRw5UitWrNDVV1/dpJ/B29v7ghkb9hM40883RPy5I0eO6G9/+5t+//vf294p/vvf/641a9bo/vvvb1I2ADgTY659JpPprL2lgoODdezYsSZlAwAAbY/yBPiZrl27qnv37tqxY4ftWk1NjXbv3t1oM789e/bYPi4oKJCHh4eCg4Ntx1qeuVb/8OHD53ysjRs3auLEiTKbzbrzzjt12WWX6dChQ7bbz/WioEFwcHCjjJK0Y8eOX7zhYI8ePfThhx/qvffe07vvvqt3331XV1xxhSZPnqznnnvuF31vAPi5jj7mStLixYv129/+ttG1PXv2cNIZAAAOhGU7wDkkJSVp5cqV6tGjh/r166dXX31V1dXVuvXWW233Wblypfr27WtbQx8bG6uuXbtq4MCB8vb21sqVK3Xvvffq448/VmFh4TmXwPj7+ys/P1979+6Vm5ub1qxZo+PHj+v06dOSpC5duujUqVM6cOBAo+ngkjR16lTNmjVLAwcO1NChQ7V582bt2bNHCxcubNLPePLkSbm5uZ21Bt/d3V1BQUGNrnl5eemiiy5Sjx49mvS9AaA5OvKYK/10Us+9996rN954QzfeeKM+//xzvfvuu/rjH//YjGcRAAC0JcoT4Bzuu+8+lZeXa968eTp16pSuuuoqrVu3ToGBgbb7TJw4UbNnz9a///1v3XrrrUpPT5f00zp3s9msF198URs3btSYMWM0ZcoUFRcXn/U4M2bM0OzZszV58mT5+vpq5MiRio+PV2FhoSRpxIgRCg4OVmxsrN56661GXxsdHa3jx49rxYoVOn78uK644gq9/vrrGjhwYJN+xkceeUS9evXS4sWLW/o0AUCr6OhjbmRkpJ5//nmtXr1aS5cuVZ8+ffTCCy9o2LBhTX4OAQBA2zJZm7IYF0AjISEhev311/XrX//a6CgA4PIYcwEAgNHY8wQAAAAAAMAOyhMAAAAAAAA7WLYDAAAAAABgBzNPAAAAAAAA7HCa03bWrFmjTz75RDU1NYqLi9OVV16p5ORk9evXT5IUFxen8ePHGxsSAAAAAAC4HKcoT3Jzc5Wfn6+33npLlZWVWrt2rSRp6tSpSkpKavL3GT58uHr16tVWMQGgQ/j++++Vm5t7wfsx5gLAL9fUMRcA0Lacojz54osvNGjQIE2fPl3l5eVKTU3Vpk2bdODAAVksFgUFBSk9PV2+vr52v0+vXr20efPmdkoNAK7p9ttvb9L9GHMB4Jdr6pgLAGhbTrHnSWlpqQoKCrR8+XItWLBAKSkpCgsLU2pqqrKzs9WnTx+tWrXK6JgAAAAAAMAFOUV54u/vr+uvv16enp4KDg6Wl5eXbrzxRoWGhkqSxo4dq8LCQoNTAgAAAAAAV+QU5ck111yjzz//XFarVceOHVNlZaUefPBB7dy5U5K0fft2DR482OCUAAAAAADAFTnFnic33XST8vLyNGnSJFmtVs2dO1fdunWT2WyWh4eHAgMDZTabjY4JAAAAAABckFOUJ5KUmpp61rX169cbkASAM7DkVSg7p0yHimvUt6eH4mP8FBXhY3QsAAAAAE7IacoTAGgqS16F1m45qZSEAA0Z4KVdRdVallUiSRQoAAAAAJrNKfY8AYDmyM4pU0pCgMJDvOXuZlJ4iLdSEgKUnVNmdDQAAAAATojyBIDLOVRcoyEDvBpdGzLAS4eKawxKBAAAAMCZUZ4AcDl9e3poV1F1o2u7iqrVt6eHQYkAAAAAODPKEwAuJz7GT8uySpS/r0q1dVbl76vSsqwSxcf4GR0NAAAAgBNiw1gALqdhU9jMDaW203aSYv3ZLBYAAABAi1CeAHBJURE+lCUAAAAAWgXLdgAAAAAAAOygPAEAAAAAALCD8gQAAAAAAMAOyhMAAAAAAAA7KE8AAAAAAADsoDwBAAAAAACwg/IEAAAAAADADsoTAAAAAAAAOyhPAAAAAAAA7HA3OgAAAACcmyWvQtk5ZTpUXKO+PT0UH+OnqAgfo2MBANBqnKY8WbNmjT755BPV1NQoLi5OkZGRmjVrlkwmkwYOHKh58+apUycm0gAAAPt4od+6LHkVWrvlpFISAjRkgJd2FVVrWVaJJPG8AgBchlO0Dbm5ucrPz9dbb72ldevWqbi4WIsWLdLMmTP15ptvymq1ymKxGB0TAAA4uIYX+o/cdbFyVvTRI3ddrLVbTsqSV2F0NKeVnVOmlIQAhYd4y93NpPAQb6UkBCg7p8zoaAAAtBqnKE+++OILDRo0SNOnT1dycrJuvPFG7d69W5GRkZKkUaNG6csvvzQ4JQAAcHS80G99h4pr9MPJWiWZj2rM9ENKMh/VDydrdai4xuhoAAC0GqdYtlNaWqojR47oD3/4gw4fPqyHHnpIVqtVJpNJkuTj46NTp04ZnBIAADi6Q8U1GjLAq9G1IQO8eKH/CwRc5Ka1f/63UhP/s2znd+tKFHCRm9HRAABoNU4x88Tf31/XX3+9PD09FRwcLC8vr0ZlSUVFhfz8/AxMCAAAnEHfnh7aVVTd6Nquomr17elhUCLn5+YmpSY2ns2TmhggN7oTAIALcYry5JprrtHnn38uq9WqY8eOqbKyUtdee61yc3MlSdu2bdOwYcMMTgkAABxdfIyflmWVKH9flWrrrMrfV6VlWSWKj+FNmJb614m6c87m+deJOoMSAQDQ+pxi2c5NN92kvLw8TZo0SVarVXPnzlXv3r2VkZGhF154QcHBwYqOjjY6JgAAcHBRET7a/W215r/yg8p/rJdvl06KiujCqTC/QMNsnvAQb9s1ZvMAAFyNU5QnkpSamnrWtaysLAOSAAAAZ2XJq1BuQaXmPxDY6FhdS3AFBUoLNczm+flRxUmx/kZHAwCg1ThNeQIAAPBLnXnajiTbaTuZG0opT1qo4XnL3FCqQ8U16tvTQ0mx/jyfAACXQnkCwCVZ8iqUnVNm+0M+PsaPP+QBcNpOG4mK8GGMBQC4NMoTAC7HklehtVtOnjWFXBJ/3AMdHPtzAACAlnCK03YAoDnOnJbfcGxmSkKAsnPKjI4GwGCcttM2LHkVSjIf1Zjph5RkPipLXoXRkQAAaFXMPAHgcpiWD+B8OG2n9THbDwDQETDzBIDLaZiWfyam5QOQGp+2szWzj+Y/EKjcgkpmSvwCzPYDAHQElCcAXA7T8gGcDy/0Wx+z/QAAHQHLdgC4HI7NBHA+vNBvfWzCCwDoCChPALgkjs0EcC680G99DbP9fr7nSVKsv9HRAABoNZQnAACgw+CFfutjth8AoCOgPAEAAB0GL/TbBrP9AACujvIEAAB0KLzQBwAAzUV5AsAlWfIqlJ1TZntnOT7GjxdLAAAAAFqE8gSAy7HkVWjtlpNn7WkgiQIFAAAAQLN1MjoAALS27JwypSQEKDzEW+5uJoWHeCslIUDZOWVGRwMAAADghChPALicQ8U1GjLAq9G1IQO8dKi4xqBEAAAAAJwZ5QkAl9O3p4d2FVU3urarqFp9e3oYlAgAAACAM6M8AeBy4mP8tCyrRPn7qlRbZ1X+viotyypRfIyf0dEAAAAAOCGn2TB2woQJ6tq1qySpd+/eSkhIUHJysvr16ydJiouL0/jx4w1MCMBRNGwKm7mh1HbaTlKsP5vFAkAb4YQzAICrc4rypLr6p+n369ats13buHGjpk6dqqSkJKNiAQAAdHiccAYA6AicojzZu3evKisrlZSUpNraWj3xxBMqKCjQgQMHZLFYFBQUpPT0dPn6+hodFYAD4A95AGg/Z55wJsl2wlnmhlLGXACAy3CKPU+8vb01bdo0vfbaa1qwYIFSUlI0ePBgpaamKjs7W3369NGqVauMjgnAQXBUMQC0H044AwB0BE5RnvTv31+xsbEymUzq37+//P39NXLkSIWGhkqSxo4dq8LCQoNTAnAU/CEPAO2HE84AAB2BU5QnmzZt0uLFiyVJx44dU3l5uaZPn66dO3dKkrZv367BgwcbGRGAA+EPeQD2WPIqlGQ+qjHTDynJfFSWvAqjIzk1TjgDAHQETrHnyaRJkzR79mzFxcXJZDJp4cKF8vLyktlsloeHhwIDA2U2m42OCcBBNPwh//M9T5Ji/Y2OBsBg7InU+jjhDADQEThFeeLp6annn3/+rOvr1683IA0AR8cf8gDOh81N20ZUhA/PHwDApTlFeQIAzcUf8nAllrwKZeeU2crA+Bg//vtuIfZEAgAALUF5AgCAA2OZSetq2BOpYeaJxJ5IAADgwpxiw1gAADoqjt5uXWxuCgAAWoKZJwAAODCWmbQu9kQCAAAtQXkCAIADY5lJ62NPJAAA0Fws2wEAwIGxzAQAAMB4zDwBAMCBscwEAADAeMw8AQAAAAAAsIOZJ4ADsORVKDunzPaucnyMH+8qA5D00/iwatMJdfFykySdrrFq1aYTkjiqGAAAoL1QngAGs+RVaO2Wk0pJCNCQAV7aVVStZVklknhhBEB6+Z2T8vLopCfju9nGiN+tK9HL75xkjAAAAGgnLNsBDJadU6aUhACFh3jL3c2k8BBvpSQEKDunzOhoABxAeWW9UhMbjxGpiQEqr6w3OhoAAECHQXkCGOxQcY2GDPBqdG3IAC8dKq4xKBEAR1JVbT3nGFFVbTUoEQAAQMdDeQIYrG9PD+0qqm50bVdRtfr29DAoEQBH0iPA7ZxjRI8AN4MSAQAAdDyUJ4DB4mP8tCyrRPn7qlRbZ1X+viotyypRfIyf0dEAOIARoZ313Os/NBojnnv9B40I7Wx0NAAAgA6DDWMBg0VF+Gj3t9Wa/8oPKv+xXr5dOikqogsbQQKQJP39H9W6+de+ytxQajuR6+Zf++qvOyqNjgYAANBhUJ4ABrPkVSi3oFLzHwhsdNqOJbiCAgWADhXX6OX0npoW62+7Vltn1Vtb2VQaAACgvbBsBzAYp+0AsId9kQAAAIznNDNPJkyYoK5du0qSevfureTkZM2aNUsmk0kDBw7UvHnz1KkTXRCcD6ftALCnYV+klISARrPTks6YiQIAAIC25RTlSXX1T++4rVu3znYtOTlZM2fO1PDhwzV37lxZLBaNHTvWqIhAizW8qxwe4m27xrvKABo0LN87c8+TpFh/lvUBAAC0I6coT/bu3avKykolJSWptrZWTzzxhHbv3q3IyEhJ0qhRo/TXv/6V8gROiXeVAVxIVIQPZQkAAICBnKI88fb21rRp03TnnXfqu+++0wMPPCCr1SqTySRJ8vHx0alTpwxOCbQM7yoDAAAAgGNzivKkf//+CgoKkslkUv/+/eXv76/du3fbbq+oqJCfn5+BCYFfhneVAQAAAMBxOcUOq5s2bdLixYslSceOHVN5ebmuu+465ebmSpK2bdumYcOGGRkRAAAAAAC4KKeYeTJp0iTNnj1bcXFxMplMWrhwoS6++GJlZGTohRdeUHBwsKKjo42OCQAAAAAAXJBTlCeenp56/vnnz7qelZVlQBoAAAAAANCROMWyHQAAAAAAAKNQngAAAAAAANhBeQIAAAAAAGCHU+x5AseTmnlMhQdOq7LKqs7eJl3Z31O/e6SH0bGcliWvQtk5ZTpUXKO+PT0UH+PH0cUAAAAA4CAoT9BsqZnH9O33NTL/truGDPDSrqJqPff6D0rNPEaB0gKWvAqt3XJSKQkBtudzWVaJJFGgAAAAAIADYNkOmq3wwGnNmRqo8BBvubuZFB7irTlTA1V44LTR0ZxSdk6ZUhICGj2fKQkBys4pMzoaAAAAAECUJ2iByiqrhgzwanRtyAAvVVZZDUrk3A4V15zz+TxUXGNQIgAAAADAmShP0GydvU3aVVTd6Nquomp19jYZlMi59e3pcc7ns29PD4MSAQAAAADORHmCZruyv6eee/0H5e+rUm2dVfn7qvTc6z/oyv6eRkdzSvExflqWVdLo+VyWVaL4GD+jowEAAAAAxIaxaIHfPdJDqZnHlLHmOKfttIKGTWEzN5TaTttJivVns1gAAAAAcBCUJ2gRipLWFRXhQ1kCAAAAAA6KZTsAAAAAAAB2UJ4AAAAAAADYQXkCAAAAAABgB+UJAAAAAACAHZQnAAAAAAAAdnDaDuAALHkVys4psx1VHB/jx+k7AAAAAOAgnKo8KSkp0e233661a9eqqqpKycnJ6tevnyQpLi5O48ePNzYg0AKWvAqt3XJSKQkBGjLAS7uKqrUsq0SSKFAAAAAAwAE4TXlSU1OjuXPnytvbW5JUWFioqVOnKikpyeBkwC+TnVOmlIQAhYf89N92eIi3UhIClLmhlPIEAAAAAByA0+x5smTJEk2ePFmXXHKJJKmgoEB/+ctfFB8fr/T0dJWXlxucEGiZQ8U1GjLAq9G1IQO8dKi4xqBEAAAAAIAzOUV5snnzZnXr1k0jR460XQsLC1Nqaqqys7PVp08frVq1ysCEQMv17emhXUXVja7tKqpW354eBiUCAAAAAJzJKcqTt99+W19++aUSExO1Z88epaWladSoUQoNDZUkjR07VoWFhQanBFomPsZPy7JKlL+vSrV1VuXvq9KyrBLFx/gZHQ0AAAAAICfZ8yQ7O9v2cWJioubPn6+HH35YGRkZCgsL0/bt2zV48GADEwIt17CvSeaGUttpO0mx/ux3AgAAAAAOwinKk3OZP3++zGazPDw8FBgYKLPZbHQkoMWiInwoSwAAAADAQTldebJu3Trbx+vXrzcwScdmyatQdk6ZbaZEfIwfL/4BAAAAAC7J6coTGM+SV6G1W04qJSFAQwZ4aVdRtZZllUgSBQoAAAAAwOU4xYaxcCzZOWVKSQhQeIi33N1MCg/xVkpCgLJzyoyOBgAAAABAq6M8QbMdKq7RkAFeja4NGeClQ8U1BiUCAAAAAKDtUJ6g2fr29NCuoupG13YVVatvTw+DEgEAAAAA0HYoT9Bs8TF+WpZVovx9Vaqtsyp/X5WWZZUoPsbP6GgAAAAAALQ6NoxFszVsCpu5odR22k5SrD+bxQIAAAAAXBLlCVokKsKHsgQAAAAA0CGwbAcAAAAAAMAOyhMAAAAAAAA7KE8AAAAAAADsoDwBAAAAAACwg/IEAAAAAADADsoTAAAAAAAAOyhPAAAAAAAA7KA8AQAAAAAAsIPyBAAAAAAAwA7KEwAAAAAAADucqjwpKSnRDTfcoP379+vgwYOKi4vTPffco3nz5qm+vt7oeAAAAAAAwAU5TXlSU1OjuXPnytvbW5K0aNEizZw5U2+++aasVqssFovBCQEAAAAAgCtymvJkyZIlmjx5si655BJJ0u7duxUZGSlJGjVqlL788ksj4wEAAAAAABflFOXJ5s2b1a1bN40cOdJ2zWq1ymQySZJ8fHx06tQpo+IBAAAAAAAX5m50gKZ4++23ZTKZtH37du3Zs0dpaWk6ceKE7faKigr5+fkZmBAAAAAAALgqpyhPsrOzbR8nJiZq/vz5Wrp0qXJzczV8+HBt27ZNI0aMMDAhAAAAAABwVU6xbOdc0tLSlJmZqbvvvls1NTWKjo42OhIAAAAAAHBBTjHz5Ezr1q2zfZyVlWVgEgAAAAAA0BE47cwTAAAAAACA9kB5AgAAAAAAYAflCQAAAAAAgB2UJwAAAAAAAHZQngAAAAAAANhBeQIAAAAAAGAH5QkAAAAAAIAdlCcAAAAAAAB2UJ4AAAAAAADYQXkCAAAAAABgB+UJAAAAAACAHZQnAAAAAAAAdlCeAAAAAAAA2EF5AgAAAAAAYAflCQAAAAAAgB2UJwAAAAAAAHZQngAAAAAAANjhbnSApqirq9PTTz+tAwcOyM3NTYsWLdKpU6eUnJysfv36SZLi4uI0fvx4Y4MCAAAAAACX4xTlyaeffipJWr9+vXJzc7Vo0SKNHj1aU6dOVVJSksHpAAAAAACAK3OK8mTMmDG68cYbJUlHjhxRYGCgCgoKdODAAVksFgUFBSk9PV2+vr7GBgUAAAAAAC7HafY8cXd3V1pamsxms6KjoxUWFqbU1FRlZ2erT58+WrVqldERAQAAAACAC3Ka8kSSlixZoq1btyojI0PXX3+9QkNDJUljx45VYWGhwekAAAAAAIArcory5N1339WaNWskSZ07d5bJZNKMGTO0c+dOSdL27ds1ePBgIyMCAAAAAAAX5RR7nowbN06zZ89WfHy8amtrlZ6erksvvVRms1keHh4KDAyU2Ww2OiYAAAAAAHBBTlGedOnSRcuXLz/r+vr16w1IAwAAAAAAOhKnWLYDAAAAAABgFMoTAAAAAAAAOyhPAAAAAAAA7KA8AQAAAAAAsIPyBAAAAAAAwA7KEwAAAAAAADsoTwAAAAAAAOygPAEAAAAAALCD8gQAAAAAAMAOyhMAAAAAAAA7KE8AAAAAAADsoDwBAAAAAACwg/IEAAAAAADADsoTAAAAAAAAOyhPAAAAAAAA7KA8AQAAAAAAsIPyBAAAAAAAwA53owM0RV1dnZ5++mkdOHBAbm5uWrRokaxWq2bNmiWTyaSBAwdq3rx56tSJLggAAAAAALQupyhPPv30U0nS+vXrlZubaytPZs6cqeHDh2vu3LmyWCwaO3aswUkBAAAAAICrcYqpGmPGjJHZbJYkHTlyRIGBgdq9e7ciIyMlSaNGjdKXX35pZEQAAAAAAOCinKI8kSR3d3elpaXJbDYrOjpaVqtVJpNJkuTj46NTp04ZnBAAAAAAALgipylPJGnJkiXaunWrMjIyVF1dbbteUVEhPz8/A5MBAAAAAABX5RTlybvvvqs1a9ZIkjp37iyTyaTQ0FDl5uZKkrZt26Zhw4YZGREAAAAAALgop9gwdty4cZo9e7bi4+NVW1ur9PR0XXbZZcrIyNALL7yg4OBgRUdHGx0TAAAAAAC4IKcoT7p06aLly5efdT0rK8uANAAAAAAAoCNximU7AAAAAAAARqE8AQAAAAAAsIPyBAAAAAAAwA7KEwAAAAAAADsoTwAAAAAAAOygPAEAAAAAALCD8gQAAAAAAMAOyhMAAAAAAAA7KE8AAAAAAADsoDwBAAAAAACwg/IEAAAAAADADvfmfoHVapXVam10rVMnOhgAAAAAAOCamlSe7NmzRxkZGdq3b59qa2vPeTsAAAAAAIAralJ5Mnv2bPn7++v3v/+9unbt2taZAAAAAAAAHEaTypNvv/1Wf/7znxUUFNTWeQAAAAAAABxKkzYrCQ0NVVFRUVtnAQAAAAAAcDjnnXmyadMm28dDhw7VrFmzNHnyZPXp0+esDWInTZrUdgkBAAAAAAAMdN7yZPXq1Y0+79q1q95///2z7mcymdq8PKmpqVF6erq+//57nT59Wg899JB69uyp5ORk9evXT5IUFxen8ePHt2kOAAAAAADQ8Zy3PPnkk0+a9A1OnDjRamHOZ8uWLfL399fSpUtVWlqqiRMnavr06Zo6daqSkpLa/PEBAAAAAEDH1aQ9T6644opzliSHDx9WVFRUq4f6uZiYGD322GO2z93c3FRQUKC//OUvio+PV3p6usrLy9s8BwAAAAAA6HjOO/PknXfese17YrVa9dBDD8ndvfHdjx8/rksuuaRtE0ry8fGRJJWXl+vRRx/VzJkzdfr0ad15550KDQ3VSy+9pFWrViktLa3NswAAAAAAgI7lvOVJdHS0vv/+e0nS3/72N1199dW2EqOBj4+Pxo0b17YJ/8/Ro0c1ffp03XPPPbrttttUVlYmPz8/SdLYsWNlNpvbJQcAAAAAAOhYzluedOnSRTNmzJAk9erVS7fccos8PT3bLdiZfvjhByUlJWnu3Lm69tprJUnTpk1TRkaGwsLCtH37dg0ePPiC3+f777/X7bff3tZxAcClNRTrTbkfYy4A/DJNHXMBAG3LZLVarRe607vvvnve2zw9PdW9e3cNHTq0zcqVZ599Vh988IGCg4Nt12bOnKmlS5fKw8NDgYGBMpvN8vX1bZPHBwAAAAAAHVeTypMpU6bo66+/lpeXl/r37y+r1aqDBw+qsrJSvXv31smTJ9W1a1e98soruuyyy9ojNwAAAAAAQLto0mk7gwYN0g033KDPPvtMmzdv1jvvvKNt27Zp3LhxGjNmjL766ivddNNNWrhwYVvnBQAAAAAAaFdNmnkybNgw/elPfzprVsn+/ft199136+uvv9bBgwc1YcIE5efnt1lYAAAAAACA9takmSddunTRN998c9b1/fv32/Y5+fHHH+Xt7d266QAAAAAAAAx23tN2zpSUlKQ5c+Zo7969Cg0NlSQVFBQoOztb06ZNU3FxsebNm6cbbrihTcMCAAAAAAC0tyYt25GkLVu26M0339S+ffvk7u6uAQMGKDExUePHj1deXp4+/vhjPfbYY+rSpUtbZwYAAAAAAGg3TS5PAAAAAAAAOqImLduRpM8//1y7du1SbW2tft63PPbYY60eDAAAAAAAwBE0qTxZuHChsrKydPnll8vHx6fRbSaTqU2CAQAAAAAAOIImLdu5/vrr9fjjj+uOO+5oj0wAAAAAAAAOo0lHFdfV1enqq69u6ywAAAAAAAAOp0nlSUJCglauXKmKioq2zgMAAAAAAOBQmrRs55577tHOnTtVX1+viy++WB4eHo1u/8tf/tJW+QAAAAAAAAzVpPLknXfesXv7xIkTWy0QAAAAAACAI2lSedKgvr5e33//vS699FLV19fL09OzLbMBAAAAAAAYrkl7ntTW1mrJkiUaOnSooqOjdfToUT311FN68skn9eOPP7Z1RgAAAAAAAMM0qTxZvny5vvjiC61du1ZeXl6SpMTERBUWFmrx4sVtGhAAAAAAAMBITSpP3n//fc2fP18RERG2a8OGDdPChQv10UcftVk4AAAAAAAAozWpPCktLVVAQMBZ1zt37qyqqqpWDwUAAAAAAOAomlSeXHvttXrllVd05t6yp06d0gsvvKARI0a0WTgAAAAAAACjNem0nWPHjmn69Ok6fPiwysrK1K9fPx09elR9+vTRSy+9pF69erVHVgAAAAAAgHbXrKOKt2/frm+//Va1tbXq37+/rr/+enXq1KTJKwAAAAAAAE7pvM1HfX39Wf+GDx+uuLg4JSYm6vrrr7fdD3Amo0ePVkhIiO3f4MGDFRUVpZdffrldHj8zM1NxcXFNuu/p06e1fv162+eJiYl68cUX2yraOR08eFBhYWGqra1t18cF4BoYc5tm+/btuuOOOxQeHq7o6Ght3LixXR4XgGthzG2aTz/9VLfddpvCwsL0m9/8Rp999lm7PC6cm/v5brjyyitlMpma9E327NnTaoGA9jBr1izdeuutkqTa2lp99dVXmjNnji655BJNmDDB2HBneP/997V69WpNnjxZ0k+/kDw8PNrt8Y8eParf/va3qq6ubrfHBOB6GHPt++677/Tb3/5WDz/8sG6++Wbt2LFDc+bMUUBAgEaPHt3mjw/AtTDm2ldUVKRHH31UaWlpuuGGG/Thhx9q+vTp+uCDD9SnT582f3w4r/OWJ2+88UZ75gDala+vr7p37277fOLEifrv//5vffjhhw71S+Xnq+r8/f3b7bE//vhjZWRkNHqeAKAlGHPt+5//+R9dccUVSk5OliQFBQUpLy9Pf/7znylPADQbY659R48eVUJCghISEiRJ06ZN05o1a/T3v/+d8gR2nXfZTmRkZJP/Aa7A3d1dnp6ekn5ajvbqq69qzJgxCgsLU0JCgvbu3Wu7b0hIiDZs2KCxY8cqPDxcTzzxhMrLyyVJubm5CgkJabTMZdasWUpJSTnn47799tu6+eabFRoaquHDh2vevHmqra1Vbm6uZs+erWPHjikkJESHDx8+azrj5s2bNX78eIWFhen2229Xbm6u7bbRo0crKytLkydP1tChQzVp0iTt3LnTdntiYqJmzZp13ufj888/1+OPP645c+Y085kEgAtjzP2Pm2++WRkZGY2umUwmZv0BaDWMuf8xcuRIpaWlSZJqamq0ceNGnT59WldddVUzn1V0NOz2ig6vpqZGH374of7617/a3uFbtWqV1q5dq9mzZ+udd95R7969df/999t+cUjSihUrlJ6erjfeeEPffPONnn766WY/9tdff60FCxbo8ccf19atW7VgwQJt3rxZH374ocLDw5Wenq7u3bvriy++0KWXXtroazdv3qxnnnlGDz74oN577z1dd911evDBB3XkyBHbfVauXKn7779ff/rTn+Tl5SWz2Wy7LTMz024xsmDBAt11113N/pkAwB7G3LP1799foaGhts9/+OEHvf/++4qIiGj2zwgAZ2LMPb/9+/dr6NChevrpp/Xwww8z6wQXRHmCDumZZ55ReHi4wsPDFRYWprS0NN17772KjY2V1WpVVlaWZsyYoaioKF122WUym81yd3fXe++9Z/se999/v2666SYNGTJEc+bM0datW3Xy5Mlm5fD29tZzzz2ncePGqVevXoqJidGVV16poqIieXp6qmvXrurUqZO6d+8uNze3Rl+7bt06xcfHa8KECerfv7+efPJJXX755Vq3bp3tPhMmTNCYMWN0+eWXKykpSQUFBbbb/P391bVr15Y9gQDQDIy5TR9zf/zxR82YMUOXXHKJbR8AAGgOxtymjbndu3fXpk2blJGRoRUrVmjr1q3N+vnQ8Zx3zxPAlc2YMUMxMTGSJC8vr0aDdklJiU6ePKmhQ4fa7u/h4aHQ0FDt37/fdi08PNz2cWhoqOrr63XgwIFm5QgNDZW3t7dWrFihoqIi7du3TwcPHtSIESMu+LX79+/XQw891OjaVVddpW+//db2+ZkNuq+vr+rr61VXV3fWLygAaEuMuU0bc0+dOqXf/va3Onz4sN5880117ty5iT8ZAPwHY27Txlw/Pz9deeWVuvLKK/WPf/xDWVlZio6ObuJPh46ImSfokLp166agoCAFBQWpZ8+ejQZZLy+vc35NXV2d6urqbJ+f+TUNR3Z36tTpnKdUne+Y388//1wTJ07U8ePHNXLkSK1YsUJXX311k34Gb2/vC2ZsWNt6pp9vzgUAbY0x98JOnDihKVOm6J///KfeeOMN9e3bt0lfBwA/x5hr3969e/X//t//a3RtwIABKi0tbVI2dFxNKk8uv/xyXXHFFef8FxYWpqioKC1btkw1NTVtnRdoc127dlX37t21Y8cO27Wamhrt3r1b/fv3t10784jugoICeXh4KDg42HbE2pnrRg8fPnzOx9q4caMmTpwos9msO++8U5dddpkOHTpku93eceHBwcGNMkrSjh07GmUEAEfHmCudPn1aycnJKi0tVXZ2toKDg3/x9wSAc2HMlT744APNnz+/0bXdu3cz9uKCmrRsZ/78+Vq5cqUeeeQRXXXVVbJarSooKFBmZqbuuOMODRo0SKtWrZLVatVTTz3V1pmBNpeUlKSVK1eqR48e6tevn1599VVVV1fr1ltvtd1n5cqV6tu3r209Z2xsrLp27aqBAwfK29tbK1eu1L333quPP/5YhYWF6t2791mP4+/vr/z8fO3du1dubm5as2aNjh8/rtOnT0uSunTpolOnTunAgQNnbWI1depUzZo1SwMHDtTQoUO1efNm7dmzRwsXLmzSz3jy5Em5ubmx7wkAw3X0MfePf/yjdu/erVdffVWdO3fW8ePHJf00lb49j6gH0DF09DF30qRJ+q//+i+9+OKLmjhxorZt26b3339ff/rTn5rzNKIDalJ58tprr2nhwoUaNWqU7drll1+uX/3qV5o/f74+/PBD9ejRQ4888gjlCVzCfffdp/Lycs2bN0+nTp3SVVddpXXr1ikwMNB2n4kTJ2r27Nn697//rVtvvVXp6emSflpzaTab9eKLL2rjxo0aM2aMpkyZouLi4rMeZ8aMGZo9e7YmT54sX19fjRw5UvHx8SosLJQkjRgxQsHBwYqNjdVbb73V6Gujo6N1/PhxrVixQsePH9cVV1yh119/XQMHDmzSz/jII4+oV69eWrx4cUufJgBoFR19zM3JyVFtba3uu+++Rtevvvrqs3IAwC/V0cfcPn366NVXX9WiRYv0+uuvq0+fPlqxYoUGDx7c5OcQHZPJ2oSFYeHh4Vq/fr1CQkIaXd+7d68mT56sv//97zp8+LBuu+025efnt1lYwFGEhITo9ddf169//WujowCAy2PMBYD2w5gLnFuT9jy57rrrtGDBgkZr1A4dOqTnnntOI0aMUF1dnTZt2qRBgwa1WVAAAAAAAAAjNGnZzrPPPqvHH39c48aNU9euXWW1WlVRUaHrr79eZrNZn332mdavX6/Vq1e3dV4AAAAAAIB21aRlOw0OHDigf/zjH3Jzc9OAAQPUr18/SVJVVZW8vLzs7pgMAAAAAADgjJpcnlitVh07duycxxH/fHdkAAAAAAAAV9GkZTufffaZMjIybEfnNbBarTKZTI3OAW8ra9as0SeffKKamhrFxcXpyiuvVHJysm32S1xcnMaPH9/mOQAAAAAAQMfSpJknt912m4KDg/Xwww/L19f3rNt79erVJuEa5Obm6vXXX9fq1atVWVmptWvXqmfPnjp16pSSkpKa/H2GDx/e5lkBwNV9//33ys3NveD9GHMB4JdjzAWA9mNvzG3SzJN//vOfWr58uYKDg1s1WFN98cUXGjRokKZPn67y8nKlpqZq06ZNOnDggCwWi4KCgpSenn7OYudMvXr10ubNm9spNQC4pttvv71J92PMBYBfjjEXANqPvTG3SUcVDxs2TDt37my1QM1VWlqqgoICLV++XAsWLFBKSorCwsKUmpqq7Oxs9enTR6tWrTIsHwAAAAAAcF1Nmnly9dVXa8GCBfrkk0/Ut29feXh4NLr9sccea5NwDfz9/RUcHCxPT08FBwfLy8tLN954owICAiRJY8eOldlsbtMMAAAAAACgY2rSzJPt27crNDRUpaWl2rFjh77++mvbv7/97W9tnVHXXHONPv/8c9uJP5WVlXrwwQdts2G2b9+uwYMHt3kOAAAAAADQ8TRp5sm6devaOoddN910k/Ly8jRp0iRZrVbNnTtX3bp1k9lsloeHhwIDA5l5AgAAAAAA2sR5y5NNmzYpNjZWnp6e2rRp03m/gclk0h133NEm4c6Umpp61rX169e3+eMCAAAAAICO7bzlyerVqxUVFSVPT0+tXr36vN+gvcoTAAAAAAAAI5y3PPnkk0/O+TEAAAAAAEBHct7yZMSIEYqIiFBkZKSGDx+uQYMGtWcuAAAAAAAAh3De8uSJJ57Qjh07tH79ei1cuFAXXXSRrUyJjIxUSEhIe+YEAAAAAAAwxHnLk7vuukt33XWXJKm8vFz5+fnasWOHPvvsM61cuVImk0kREREaPny4EhIS2i0wAAAAAABAe2rSUcW+vr4aOXKkRo4cqRMnTuhvf/ubtm7dKovFoo8//pjyBAAAAAAAuKwLlielpaXKzc3VV199pa+++koHDx5UUFCQIiMj9eyzzyoyMrI9cgIAAAAAABjivOXJokWL9NVXX+mbb75Rnz59FBkZqRkzZigyMlKXXHJJe2YEAAAAAAAwzHnLk//6r/9Sjx49lJaWpttuu03dunVrz1wAAAAAAAAO4bzlycsvv6yvvvpKW7Zs0dKlS9W/f38NHz5cw4cPV0REhPz9/dsxJgAAAAAAgDHOW56MGjVKo0aNkiSVlZXpf//3f5Wbm6sVK1Zo//79GjRokK1MGT16dLsFBgAAAAAAaE9NOm3Hz89PY8aM0ZgxYyRJhYWFeuutt7Rhwwa98cYb2rNnT5uGBAAAAAAAMMoFy5Py8nLt2rVLO3bs0M6dO7Vz5079+9//1uDBgxUfH6+IiIj2yAkAAAAAAGCI85Yns2bN0o4dO3Tw4EF5enpq6NChioiIUGJiosLDw+Xt7d2eOQEAAAAAAAxx3vLkxIkTmjhxoiIiIjRkyBC5uzdphQ8AAAAAAIBLsXvaDgAAAAAAQEfXyegAAAAAAAAAjsxp1uKsWbNGn3zyiWpqahQXF6fIyEjNmjVLJpNJAwcO1Lx589SpE10QAAAAAABoXU7RNuTm5io/P19vvfWW1q1bp+LiYi1atEgzZ87Um2++KavVKovFYnRMAAAAAADggpyiPPniiy80aNAgTZ8+XcnJybrxxhu1e/duRUZGSpJGjRqlL7/80uCUAAAAAADAFTnFsp3S0lIdOXJEf/jDH3T48GE99NBDslqtMplMkiQfHx+dOnXK4JQAAAAAAMAVOUV54u/vr+DgYHl6eio4OFheXl4qLi623V5RUSE/Pz8DEwIAAAAAAFflFMt2rrnmGn3++eeyWq06duyYKisrde211yo3N1eStG3bNg0bNszglAAAAAAAwBU5xcyTm266SXl5eZo0aZKsVqvmzp2r3r17KyMjQy+88IKCg4MVHR1tdEwAAAAAAOCCnKI8kaTU1NSzrmVlZRmQBAAAAAAAdCROsWwHAAAAAADAKJQnAAAAAAAAdlCeAAAAAAAA2EF5AgAAAAAAYAflCQAAAAAAgB2UJwAAAAAAAHZQngAAAAAAANhBeQIAAAAAAGAH5QkAAAAAAIAdlCcAAAAAAAB2UJ4AAAAAAADYQXkCAAAAAABgB+UJAAAAAACAHZQnAAAAAAAAdlCeAAAAAAAA2EF5AgAAAAAAYAflCQAAAAAAgB3uRgdoqgkTJqhr166SpN69eyshIUHJycnq16+fJCkuLk7jx483MCEAAAAAAHBFTlGeVFdXS5LWrVtnu7Zx40ZNnTpVSUlJRsUCAAAAAAAdgFOUJ3v37lVlZaWSkpJUW1urJ554QgUFBTpw4IAsFouCgoKUnp4uX19fo6MCcBCWvApl55TpUHGN+vb0UHyMn6IifIyOBQAAAMAJOUV54u3trWnTpunOO+/Ud999pwceeEAPPvig7rzzToWGhuqll17SqlWrlJaWZnRUAA7AklehtVtOKiUhQEMGeGlXUbWWZZVIEgUKAAAAgGZzig1j+/fvr9jYWJlMJvXv31/+/v4aOXKkQkNDJUljx45VYWGhwSkBOIrsnDKlJAQoPMRb7m4mhYd4KyUhQNk5ZUZHAwAAAOCEnKI82bRpkxYvXixJOnbsmMrLyzV9+nTt3LlTkrR9+3YNHjzYyIgAHMih4hoNGeDV6NqQAV46VFxjUCIAAAAAzswplu1MmjRJs2fPVlxcnEwmkxYuXCgvLy+ZzWZ5eHgoMDBQZrPZ6JgAHETfnh7aVVSt8BBv27VdRdXq29PDwFQAAAAAnJVTlCeenp56/vnnz7q+fv16A9IAcHTxMX5allVy1p4nSbH+RkcDAAAA4IScojwBgOZo2BQ2c0Op7bSdpFh/NosFAAAA0CKUJwBcUlSED2UJAAAAgFZBeQLAJVnyKpSdU2abeRIf40eZAgAAAKBFKE8AuBxLXoXWbjl51p4nkihQAAAAADSbUxxVDADNkZ1TppSEAIWHeMvdzaTwEG+lJAQoO6fM6GgAAAAAnBDlCQCXc6i4RkMGeDW6NmSAlw4V1xiUCAAAAIAzozwB4HL69vTQrqLqRtd2FVWrb08PgxIBAAAAcGaUJwBcTnyMn5ZllSh/X5Vq66zK31elZVklio/xMzoaAAAAACfEhrEAXE7DprCZG0ptp+0kxfqzWSwAAACAFqE8AeCSoiJ8KEsAAAAAtAqW7QAAAAAAANhBeQIAAAAAAGAH5QkAAAAAAIAdlCcAAAAAAAB2UJ4AAAAAAADYQXkCAAAAAABgB+UJAAAAAACAHZQnAAAAAAAAdrgbHaCpJkyYoK5du0qSevfureTkZM2aNUsmk0kDBw7UvHnz1KkTXRAAAAAAAGhdTlGeVFdXS5LWrVtnu5acnKyZM2dq+PDhmjt3riwWi8aOHWtURAAAAAAA4KKcYqrG3r17VVlZqaSkJE2ZMkV///vftXv3bkVGRkqSRo0apS+//NLglAAAAAAAwBU5xcwTb29vTZs2TXfeeae+++47PfDAA7JarTKZTJIkHx8fnTp1yuCUAAAAAADAFTlFedK/f38FBQXJZDKpf//+8vf31+7du223V1RUyM/Pz8CEAAAAAADAVTnFsp1NmzZp8eLFkqRjx46pvLxc1113nXJzcyVJ27Zt07Bhw4yMCAAAAAAAXJRTzDyZNGmSZs+erbi4OJlMJi1cuFAXX3yxMjIy9MILLyg4OFjR0dFGxwQAAAAAAC7IKcoTT09PPf/882ddz8rKMiANAAAAAADoSJxi2Q4AAAAAAIBRnGLmCQAAAAAAkCx5FcrOKdOh4hr17emh+Bg/RUX4GB3L5VGeAAAAAADaBC/0W5clr0Jrt5xUSkKAhgzw0q6iai3LKpEkntc2RnkCwCXxixoAAMBYvNBvfdk5ZUpJCFB4iLckKTzEWykJAcrcUMpz2sYoTwC4HH5RAwAAGI8X+q3vUHGNhgzwanRtyAAvHSquMShRx8GGsWgRS16FksxHNWb6ISWZj8qSV2F0JMDmzF/U7m4m2y/q7Jwyo6MBAAB0GLzQb319e3poV1F1o2u7iqrVt6eHQYk6DsoTNFvDu/qP3HWxclb00SN3Xay1W05SoMBh8IsaAADAeLzQb33xMX5allWi/H1Vqq2zKn9flZZllSg+xs/oaC6PZTtoNqbfwdE1/KJu+G9U4hc1AABAe2t4of/zpdRJsf5GR3NaDa+3MjeU2vb2S4r153VYO6A8QbPxrj4cHb+oAQAAjBcV4aPd31Zr/is/qPzHevl26aSoiC680P+FoiJ8eA4NQHmCZuNdfTg6GnkAAADjWfIqlFtQqfkPBDZ6Q8sSXMHfZXA6lCdotvgYPy15o0RpU/7zrv6SN0r0wAR/o6MBNjTyAAAAxmK5P1wJ5Qla5HRtvZ7PPqHiklr1DHDX6dp6oyM5NUtehbJzymyzJOJj/PiFAgAAAKfGcn+4EsoTNFt2Tpkykro3WraTv6+KBrmFGk4v+vn+HJJ4PgEAAOC0WO4PV8JRxWg2GuTWdeZ0Rnc3k206Y3ZOmdHRAAAAgBbjWF24EmaeoNlokFsXZRQAAABcEZv4w5VQnqDZOAa2dVFGAQAAwFWxiT9cBeUJmo0GuXVRRrUNNuEFAAAA0FqcqjwpKSnR7bffrrVr16qqqkrJycnq16+fJCkuLk7jx483NmAHQoPceiijWh+b8AIAAABoTU5TntTU1Gju3Lny9v5paUNhYaGmTp2qpKQkg5MBcDRnbsIrybYJLydCAQAAAGgJpylPlixZosmTJ+vll1+WJBUUFOjAgQOyWCwKCgpSenq6fH19DU4JNB+zJFofm/ACAAAAaE1OcVTx5s2b1a1bN40cOdJ2LSwsTKmpqcrOzlafPn20atUqAxMCLcdRxa2vYRPeM7EJLwAAAFyBJa9CSeajGjP9kJLMR2XJqzA6UofgFOXJ22+/rS+//FKJiYnas2eP0tLSNGrUKIWGhkqSxo4dq8LCQoNTAi3DLInW17AJb/6+KtXWWZW/r0rLskoUH+NndDQAAACgxRpmrT9y18XKWdFHj9x1sdZuOUmB0g6cYtlOdna27ePExETNnz9fDz/8sDIyMhQWFqbt27dr8ODBBiYEWo6jilsfm/ACQPvihDMAaB/s7WccpyhPzmX+/Pkym83y8PBQYGCgzGaz0ZGAFuGo4rbBiVAA0D4seRVatemEuni5SZJO11i1atMJSezdBQCtjVnrxnG68mTdunW2j9evX29gEqB1MEsCAODMXn7npLw8OunJ+G62NwF+t65EL79zkt9lANDKmLVuHKcrTwBXxCwJAICzKq+s17PJ3RtNIU9NDNDTfzhucDIAcD3MWjcO5QkAAABarKraes4p5FXVVoMSAYDrYta6cZzitB0AAAA4ph4Bbuc8Hr5HgJtBiQAAaH2UJwAAAGixabH+WvJG4+Phl7xRomlMIQeAVsdRxcZh2Q4AAABa7FxTyB+YwBRyAGgLHFVsHMoTAAAA/CJsfA4A7YOjio1DeQIAAIBfxJJXoeycMtvMk/gYP8oUAJIYH1obRxUbh/IEAAAALdaw/v7nx2ZK4gUS0MExPrQ+jio2DuUJWoQGGQAASKy/B3B+2TllGh3h02hPpNERPsrOKWN8aCGOKjYO5QmajQYZAAA0OFRcox9O1irJfNT2h3zcuK6svwegg8U1qjpdr6fOeN2wNKtEx07UGR3NqbHPlDE4qhjNduY7TO5uJts7TNk5ZUZHAwAA7SzgIjet/fO/Gx+b+ed/K+AiN6OjATCYb+dOeupnrxueSgiQb2dehsL58F8tmo0dngEAQAM3Nyk1sfGLo9TEALnRnQAdXnll/TlfN5RX1huUCGg5yhM0W8MOz2dih2cAADqmf52oO+eLo38xLR/o8ILO87ohiNcNcELseYJmi4/xk3ntcXXxclNxSa16Brjrx+o6TZ/UzehoAACgnXFsJoDz4WQYuBLKE7SIp3snPRnfzTYILnmjxOhIAADAALw4AnA+nAwDV0J5gmbLzilT2pTGRxKmTeFIwl9ixZ9OyJL3o8p/rJdvl06KiuiiR+9mJg8AwPHx4giAPZwMA1dBeYJmY8PY1rXiTye0Lf9HzX8g0PaO3XOv/yBJFCgAAKfAiyMAgKtzqg1jS0pKdMMNN2j//v06ePCg4uLidM8992jevHmqr2fH5vbChrGty5L3o+ZMDWx0SsGcqYGy5P1odDQAAAAAgJyoPKmpqdHcuXPl7f3TUpFFixZp5syZevPNN2W1WmWxWAxO2HE0rG3O31el2jqr8vdVaVlWieJj/IyO5pTKfzzPEW4/UggCAAAAgCNwmvJkyZIlmjx5si655BJJ0u7duxUZGSlJGjVqlL788ksj43UoURE+Sor1V+aGUsU8+k9lbihlbfMv4Nul0zln8vh2cZr/ewIAAACAS3OKPU82b96sbt26aeTIkXr55ZclSVarVSaTSZLk4+OjU6dOGRmxw2Ftc+uJiuii517/QXOmNt7zJCqii9HRAAAAAABykvLk7bfflslk0vbt27Vnzx6lpaXpxIkTttsrKirk58eSETinhk1h57/yA6ftAAAAAIADcoryJDs72/ZxYmKi5s+fr6VLlyo3N1fDhw/Xtm3bNGLECAMTAr/Mo3d3oywBAAAAAAfltJsqpKWlKTMzU3fffbdqamoUHR1tdCQAAAAAAOCCnGLmyZnWrVtn+zgrK8vAJAAAAAAAoCNw2pknAAAAAAAA7cHpZp4AAAAAQFux5FUoO6dMh4pr1Lenh+Jj/DhlEgDlCQAAAABIPxUna7ecVEpCgIYM8NKuomotyyqRJAoUoINj2Q4AAAAASMrOKVNKQoDCQ7zl7mZSeIi3UhIClJ1TZnQ0AAajPAEAAAAASYeKazRkgFeja0MGeOlQcY1BiQA4CsoTAAAAAJDUt6eHdhVVN7q2q6hafXt6GJQIgKOgPAEAAAAASfExflqWVaL8fVWqrbMqf1+VlmWVKD7Gz+hoAAzGhrEAAAAAoP9sCpu5odR22k5SrD+bxQKgPAEAAACABlERPpQlAM7Csh0AAAAAAAA7mHkCAAAAAP/Hkleh7Jwy27Kd+Bg/ZqIAoDwBAAAAAOmn4mTtlpNKSQjQkAFe2lVUrWVZJZJEgQJ0cCzbAQAAAABJ2TllSkkIUHiIt9zdTAoP8VZKQoCyc8qMjgbAYJQnAAAAACDpUHGNhgzwanRtyAAvHSquMSgRAEdBeQIAAAAAkvr29NCuoupG13YVVatvTw+DEgFwFOx5AgAAAACS4mP8ZF57XF283FRcUqueAe76sbpO0yd1MzoaAINRngAAAADA//F076Qn47vZNoxd8kaJ0ZEAOACnKE/q6ur09NNP68CBA3Jzc9OiRYt06tQpJScnq1+/fpKkuLg4jR8/3tigAAAAAJxWdk6Z0qb8tGGsJIWHeCttSoAyN5Ry2g4cBsdpG8MpypNPP/1UkrR+/Xrl5uZq0aJFGj16tKZOnaqkpCSD0wEAAABwBWwYC0fHcdrGcYryZMyYMbrxxhslSUeOHFFgYKAKCgp04MABWSwWBQUFKT09Xb6+vsYGBQAADo937ACcT8OGsQ0zTyQ2jIVjOfM4bUm247SZHdX2nOa0HXd3d6WlpclsNis6OlphYWFKTU1Vdna2+vTpo1WrVhkdEQAAOLiGd+weueti5azoo0fuulhrt5yUJa/C6GgAHEB8jJ+WZZUof1+Vauusyt9XpWVZJYqP8TM6GiCJ2VFGcpryRJKWLFmirVu3KiMjQ9dff71CQ0MlSWPHjlVhYaHB6QAAgKM78x07dzeT7R277Jwyo6MBcABRET5KivVX5oZSxTz6T2VuKFVSrD/v6MNhcJy2cZyiPHn33Xe1Zs0aSVLnzp1lMpk0Y8YM7dy5U5K0fft2DR482MiIAADACfCOHQDAmTE7yjhOsefJuHHjNHv2bMXHx6u2tlbp6em69NJLZTab5eHhocDAQJnNZqNjAgAAB8d+BgDsYTNOOLqG/w4zN5Ta9u5idlT7cIrypEuXLlq+fPlZ19evX29AGgAA4Kwa3rH7+QujpFh/o6MBcABsxglnEBXhw3+PBnCK8gQAAKA18I4dAHtY2gfgfChPAABAh8I7dgDOh6V9AM6H8gQAAHQolrwKZeeU2WaexMf4UaYAkMTSPgDnR3kCAAA6DDaDBGAPS/sAnA/lCQAA6DDYDBLAhbC0D8C5dDI6AAAAQHthM0gAANASlCcAAKDDaNgM8kxsBgkAAC6E8gQAAHQYDZtB5u+rUm2dVfn7qrQsq0TxMX5GRwMAAA6MPU8AAECHERXho93fVmv+Kz+o/Md6+XbppKiILuxvAAAA7KI8AQAAHYYlr0K5BZWa/0Bgo9N2LMEVFCgAAOC8KE8AAECHkZ1TptERPo2OIR0d4aPsnDLKEwAAcF6UJwAAoMM4WFyjqtP1eiohwDbzZGlWiY6dqDM6GgAAcGCUJwAAODhLXoWyc8psMyXiY/yYJdFCvp076amEAIWHeEuSwkO89VRCgOa/8oPByQAAgCPjtB0AAByYJa9Cqzad0OkaqyTpdI1VqzadkCWvwuBkzqm8sl5DBng1ujZkgJfKK+sNSgQAAJwB5QkAAA7s5XdOysujk56M76acFX30ZHw3eXl00svvnDQ6mlMK6umhXUXVja7tKqpWUE8PgxIBAABnQHkCAIADK6+sV2riT8tM3N1MCg/xVmpiADMlWig+xk/LskqUv69KtXVW5e+r0rKsEsXH+BkdDQAAODD2PAEAwIFVVVvPucykqtpqUCLn1rBXzJmn7STF+rOHDAAAsMspypO6ujo9/fTTOnDggNzc3LRo0SJZrVbNmjVLJpNJAwcO1Lx589SpExNpAACupUeAm3YVVds2OJV+WmbSI8DNwFTOLSrCh7IEAAA0i1O0DZ9++qkkaf369Xr00Ue1aNEiLVq0SDNnztSbb74pq9Uqi8VicEoAAFrftFh/LXmj8TKTJW+UaFqsv9HRAAAAOgynmHkyZswY3XjjjZKkI0eOKDAwUH/5y18UGRkpSRo1apT++te/auzYsQamBACg9Z1rmckDE1hmAgAA0J6cojyRJHd3d6Wlpemjjz7SihUr9Omnn8pkMkmSfHx8dOrUKYMTAgDQNlhmAgAAYCynKU8kacmSJUpJSdFdd92l6ur/HDNYUVEhPz92yQcAuCZLXoWyc8psM0/iY/woUwAAANqRU+x58u6772rNmjWSpM6dO8tkMik0NFS5ubmSpG3btmnYsGFGRgQAoE1Y8iq0dstJPXLXxcpZ0UeP3HWx1m45KUtehdHRAAAAOgynKE/GjRunwsJCxcfHa9q0aUpPT9fcuXOVmZmpu+++WzU1NYqOjjY6JgAArS47p0wpCQEKD/GWu5tJ4SHeSkkIUHZOmdHRAAAAOgynWLbTpUsXLV++/KzrWVlZBqQBAKD9HCqu0ZABXo2uDRngpUPFNQYlAgAA6HicYuYJAAAdVd+eHtpVVN3o2q6iavXt6WFQIgAAgI6H8gQAAAcWH+OnZVklyt9Xpdo6q/L3VWlZVoniY9goHQAAoL04xbIdAAA6qoZTdTI3lNpO20mK9ee0HQAAgHZEeQIAgIOLivChLAEAADAQy3YAAAAAAADsoDwBAAAAAACwg/IEAAAAAADADsoTAAAAAAAAOyhPAAAAAAAA7KA8AQAAAAAAsIOjigEAcHCWvApl55TpUHGN+vb0UHyMH0cXAwAAtCPKEwAAHJglr0Jrt5xUSkKAhgzw0q6iai3LKpEkChQAAIB2wrIdAAAcWHZOmVISAhQe4i13N5PCQ7yVkhCg7Jwyo6MBAAB0GJQnAAA4sEPFNRoywKvRtSEDvHSouMagRAAAAB0P5QkAAA6sb08P7SqqbnRtV1G1+vb0MCgRAABAx0N5AgCAA4uP8dOyrBLl76tSbZ1V+fuqtCyrRPExfkZHAwAA6DDYMBYAAAfWsCls5oZS22k7SbH+bBYLAADQjpyiPKmpqVF6erq+//57nT59Wg899JB69uyp5ORk9evXT5IUFxen8ePHGxsUAIA2EBXhQ1kCAABgIKcoT7Zs2SJ/f38tXbpUpaWlmjhxoqZPn66pU6cqKSnJ6HgAAAAAAMCFOUV5EhMTo+joaNvnbm5uKigo0IEDB2SxWBQUFKT09HT5+voamBIAAAAAALgip9gw1sfHR76+viovL9ejjz6qmTNnKiwsTKmpqcrOzlafPn20atUqo2MCAAAAAAAX5BTliSQdPXpUU6ZM0W9+8xvddtttGjt2rEJDQyVJY8eOVWFhocEJAQAAAACAK3KKZTs//PCDkpKSNHfuXF177bWSpGnTpikjI0NhYWHavn27Bg8efMHv8/333+v2229v67gA4NK+//77Jt+PMRcAfhnGXABoP/bGXJPVarW2Y5YWefbZZ/XBBx8oODjYdm3mzJlaunSpPDw8FBgYKLPZzJ4nAAAAAACg1TlFeQIAAAAAAGAUp9nzBAAAAAAAwAiUJwAAAAAAAHZQngAAAAAAANhBeQIAAAAAAGAH5QlabMeOHUpMTDQ6hkuoqanRU089pXvuuUeTJk2SxWIxOpLTq6ur0+zZszV58mTFx8fr0KFDRkcCfhHG3NbDmNv6GHPhahhzWw9jbutjzDWGu9EB4JxeeeUVbdmyRZ07dzY6ikvYsmWL/P39tXTpUpWWlmrixImKiooyOpZT+/TTTyVJ69evV25urhYtWqSXXnrJ4FRAyzDmti7G3NbHmAtXwpjbuhhzWx9jrjGYeYIW6du3rzIzM42O4TJiYmL02GOP2T53c3MzMI1rGDNmjMxmsyTpyJEjCgwMNDgR0HKMua2LMbf1MebClTDmti7G3NbHmGsMZp6gRaKjo3X48GGjY7gMHx8fSVJ5ebkeffRRzZw509hALsLd3V1paWn66KOPtGLFCqPjAC3GmNu6GHPbBmMuXAVjbutizG0bjLntj5kngIM4evSopkyZot/85je67bbbjI7jMpYsWaKtW7cqIyNDP/74o9FxADgIxty2wZgL4FwYc9sGY277ojwBHMAPP/ygpKQkPfXUU5o0aZLRcVzCu+++qzVr1kiSOnfuLJPJxDRRAJIYc9sCYy6A82HMbX2MucagPAEcwB/+8AeVlZVp9erVSkxMVGJioqqqqoyO5dTGjRunwsJCxcfHa9q0aUpPT5eXl5fRsQA4AMbc1seYC+B8GHNbH2OuMUxWq9VqdAgAAAAAAABHxcwTAAAAAAAAOyhPAAAAAAAA7KA8AQAAAAAAsIPyBAAAAAAAwA7KEwAAAAAAADsoTwAAAAAAAOygPAEAAAAAALCD8gQAAAAAAMAOyhMAAAAAAAA7KE8AAAAAAADsoDwBAAAAAACwg/IEAAAAAADADsoTAADgMEaPHq2QkJBz/vv44481evRobdy4sc1zfPXVV/rHP/4hSdq8ebNGjRrV5o8JAAAcl7vRAQAAAM40a9Ys3XrrrWddv+iii7Rw4cJ2yXDvvffq9ddf16BBgzR+/HjdeOON7fK4AADAMVGeAAAAh+Lr66vu3bsbHcPG29tb3t7eRscAAAAGYtkOAABwSlarVatXr9bIkSN1zTXXaNq0afruu+8kScuWLVNcXFyj+7/88su64447JEn79+/X/fffr/DwcA0ZMkRxcXH65ptvJP20dEiSpk6dqszMzLOW7ezfv1/Tpk3T1Vdfreuvv16ZmZmqr6+XJGVmZurxxx/XM888o2uuuUY33XST1qxZ09ZPBQAAaGOUJwAAwCllZWXpvffe0+9+9ztt2LBBQUFBuu+++1RZWalbbrlF+fn5+te//mW7f05Ojm655RZZrVY9/PDD+tWvfqX33ntP69evV319vX73u99JkjZt2iRJ+v3vf6+kpKRGj3nixAndc889uuSSS7Rx40bNnz9f2dnZWrt2re0+H330kdzc3LRx40bdfffdeuGFF1RUVNQOzwgAAGgrlCcAAMChPPPMMwoPD2/071wbtr766qtKSUnRtddeq8suu0wZGRlyc3PT1q1bdcUVVyg4OFgfffSRJOmf//yn9uzZo5tvvlmVlZWaNGmS0tLS1LdvXw0ePFgTJ060FRzdunWT9NMeKz4+Po0e87//+7/l5eWlZ555RpdddpnGjBmjxx57TK+++qrtPl27dtWsWbMUHBys5ORk+fv7q6CgoK2eLgAA0A7Y8wQAADiUGTNmKCYmptG1Tp0av99TUVGh4uJipaSkNLqturratnRn/Pjx+uijjxQfH6+cnByFh4fr0ksvlSTdc889eu+991RQUKBvv/1WhYWF8vf3v2C2/fv368orr5SHh4ftWnh4uEpLS3XixAlJUq9eveTm5ma73cfHRzU1Nc16DgAAgGOhPAEAAA6lW7duCgoKsnufuro6SdILL7ygAQMGNLqta9eukqRbbrlFL730kkpLS5WTk6Pbb79d0k/Fy6RJk3TRRRdpzJgxuvXWW/Xtt9/q5ZdfvmA2Ly+vs6417HfS8L9nFisAAMA1UJ4AAACn4+fnp4CAAB0/flxRUVGSfipUnnjiCU2ePFnXXnut+vfvr0GDBmnjxo3au3evoqOjJUn/+7//q+LiYm3ZssVWdHzxxReyWq0XfNzLLrtMOTk5qqmpsX1tfn6+/P39bct9AACA62HPEwAA4FDKy8t1/Pjxs/6Vl5c3ut99992n5cuX6+OPP9bBgwe1YMECffnllwoODrbdp2H2SWRkpAIDAyVJ/v7+qqys1EcffaTDhw9r48aNys7O1unTp21f16VLF33zzTc6depUo8e89dZbVV9fr7lz52r//v2yWCzKzMzU5MmTz1paBAAAXAczTwAAgENZvHixFi9efNb1n598M23aNFVWVmrBggUqKyvTFVdcoddee009evSw3eeWW27RsmXLNH78eNu18PBwzZgxQ2azWdXV1Ro0aJDmzZun2bNn68iRI/rVr36l++67T88//7yOHDmikJAQ29f6+Pjo1Vdf1bPPPqsJEyaoW7dumjJlipKTk9vgmQAAAI7CZG3KHFUAAAAAAIAOivmlAAAAAAAAdlCeAAAAAAAA2EF5AgAAAAAAYAflCQAAAAAAgB2UJwAAAAAAAHZQngAAAAAAANhBeQIAAAAAAGAH5QkAAAAAAIAdlCcAAAAAAAB2/H+Qds19P1VumgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x540 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 3, figsize=(16, 7.5), sharey=True, sharex=True)\n",
    "fig.subplots_adjust(left=0.075, right=0.975, bottom=0.09, top=0.925, wspace=0.03)\n",
    "\n",
    "axes_flat = axes.ravel()\n",
    "\n",
    "for i, pop_plot in enumerate([4, 5, 6, 1, 2, 3]):\n",
    "    ax = axes_flat[i]\n",
    "    idx = data.index[data[\"pop\"] == pop_plot].tolist()\n",
    "    elev_plot = data.loc[idx, \"elev\"].values\n",
    "    wing_plot = data.loc[idx, \"wing\"].values\n",
    " \n",
    "    # Plot observed data points\n",
    "    ax.scatter(elev_plot, wing_plot, color=\"white\", ec=\"royalblue\")\n",
    "    ax.set_xlim(0, 4)\n",
    "    # Add a title\n",
    "    ax.set_title(f\"Population: {pop_plot}\", fontsize=14)\n",
    "\n",
    "ax.xaxis.set_ticks([1, 2, 3])\n",
    "fig.delaxes(axes[0][2])\n",
    "fig.text(0.5, 0.02, \"Elevation\", fontsize=14)\n",
    "fig.text(0.03, 0.5, \"Wing length\", rotation=90, fontsize=14, va=\"center\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42056ce0-1fa5-46bd-b3a3-6b85d7a2c80f",
   "metadata": {},
   "source": [
    "### 10.5 Analysis using PyMC3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a192d39-c5e8-488d-a3e3-d7f5f1721a1a",
   "metadata": {},
   "source": [
    "#### 10.5.1 Main-effects ANOVA Using PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "85123816-7737-4772-a7c5-64ec322b52d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_pop = pd.Categorical(pop).codes\n",
    "groups_pop = len(np.unique(idx_pop))\n",
    "idx_elev = pd.Categorical(elev).codes\n",
    "groups_elev = len(np.unique(idx_elev))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a5711269-f814-4c0b-9b9c-8f8f8a19868f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y, X = dmatrices('wing ~ pop + elev', data, return_type = 'dataframe')\n",
    "X = X.rename(columns = {'pop[T.2]': 'pop2', 'pop[T.3]': 'pop3', 'pop[T.4]': 'pop4', 'pop[T.5]': 'pop5', 'elev[T.2]': 'elev2', 'elev[T.3]': 'elev3'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "15e74ce5-67f9-47d6-b6d9-e84b6bd4dcb4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, beta_elev3, beta_elev2, beta_pop5, beta_pop4, beta_pop3, beta_pop2, alpha]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [4000/4000 00:09<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 25 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as anova_2w:\n",
    "    \n",
    "    alpha = pm.Normal('alpha', 0, tau = 0.001)\n",
    "    beta_pop2 = pm.Normal('beta_pop2', 0, tau = 0.001)\n",
    "    beta_pop3 = pm.Normal('beta_pop3', 0, tau = 0.001)\n",
    "    beta_pop4 = pm.Normal('beta_pop4', 0, tau = 0.001)\n",
    "    beta_pop5 = pm.Normal('beta_pop5', 0, tau = 0.001)\n",
    "    beta_elev2 = pm.Normal('beta_elev2', 0, tau = 0.001)\n",
    "    beta_elev3 = pm.Normal('beta_elev3', 0, tau = 0.001)\n",
    "    \n",
    "    sigma = pm.Uniform('sigma', 0, 100)\n",
    "    \n",
    "    mu = alpha + beta_pop2*X['pop2'] + beta_pop3*X['pop3'] + beta_pop4*X['pop4'] + beta_pop5*X['pop5'] + beta_elev2*X['elev2'] + beta_elev3*X['elev3']\n",
    "    wing_new = pm.Normal('wing_p', mu = mu, sigma = sigma, observed = wing)\n",
    "    trace_2w = pm.sample(return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "faa5f281-a8c0-4e37-a0c9-24abf863e1d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>sd</th>\n",
       "      <th>hdi_3%</th>\n",
       "      <th>hdi_97%</th>\n",
       "      <th>mcse_mean</th>\n",
       "      <th>mcse_sd</th>\n",
       "      <th>ess_bulk</th>\n",
       "      <th>ess_tail</th>\n",
       "      <th>r_hat</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>41.473</td>\n",
       "      <td>1.293</td>\n",
       "      <td>38.970</td>\n",
       "      <td>43.713</td>\n",
       "      <td>0.049</td>\n",
       "      <td>0.035</td>\n",
       "      <td>704.0</td>\n",
       "      <td>926.0</td>\n",
       "      <td>1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_pop2</th>\n",
       "      <td>-10.650</td>\n",
       "      <td>1.591</td>\n",
       "      <td>-13.556</td>\n",
       "      <td>-7.488</td>\n",
       "      <td>0.052</td>\n",
       "      <td>0.037</td>\n",
       "      <td>929.0</td>\n",
       "      <td>1366.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_pop3</th>\n",
       "      <td>-4.577</td>\n",
       "      <td>1.550</td>\n",
       "      <td>-7.569</td>\n",
       "      <td>-1.618</td>\n",
       "      <td>0.048</td>\n",
       "      <td>0.034</td>\n",
       "      <td>1052.0</td>\n",
       "      <td>1215.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_pop4</th>\n",
       "      <td>5.562</td>\n",
       "      <td>1.549</td>\n",
       "      <td>2.872</td>\n",
       "      <td>8.527</td>\n",
       "      <td>0.051</td>\n",
       "      <td>0.036</td>\n",
       "      <td>933.0</td>\n",
       "      <td>1217.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_pop5</th>\n",
       "      <td>10.284</td>\n",
       "      <td>1.542</td>\n",
       "      <td>7.516</td>\n",
       "      <td>13.206</td>\n",
       "      <td>0.052</td>\n",
       "      <td>0.037</td>\n",
       "      <td>897.0</td>\n",
       "      <td>1198.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_elev2</th>\n",
       "      <td>3.797</td>\n",
       "      <td>1.154</td>\n",
       "      <td>1.818</td>\n",
       "      <td>6.150</td>\n",
       "      <td>0.032</td>\n",
       "      <td>0.023</td>\n",
       "      <td>1281.0</td>\n",
       "      <td>1376.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta_elev3</th>\n",
       "      <td>8.842</td>\n",
       "      <td>1.161</td>\n",
       "      <td>6.763</td>\n",
       "      <td>11.035</td>\n",
       "      <td>0.031</td>\n",
       "      <td>0.022</td>\n",
       "      <td>1405.0</td>\n",
       "      <td>1277.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>3.698</td>\n",
       "      <td>0.363</td>\n",
       "      <td>3.049</td>\n",
       "      <td>4.376</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.006</td>\n",
       "      <td>1800.0</td>\n",
       "      <td>1453.0</td>\n",
       "      <td>1.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  \\\n",
       "alpha       41.473  1.293  38.970   43.713      0.049    0.035     704.0   \n",
       "beta_pop2  -10.650  1.591 -13.556   -7.488      0.052    0.037     929.0   \n",
       "beta_pop3   -4.577  1.550  -7.569   -1.618      0.048    0.034    1052.0   \n",
       "beta_pop4    5.562  1.549   2.872    8.527      0.051    0.036     933.0   \n",
       "beta_pop5   10.284  1.542   7.516   13.206      0.052    0.037     897.0   \n",
       "beta_elev2   3.797  1.154   1.818    6.150      0.032    0.023    1281.0   \n",
       "beta_elev3   8.842  1.161   6.763   11.035      0.031    0.022    1405.0   \n",
       "sigma        3.698  0.363   3.049    4.376      0.009    0.006    1800.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "alpha          926.0   1.01  \n",
       "beta_pop2     1366.0   1.00  \n",
       "beta_pop3     1215.0   1.00  \n",
       "beta_pop4     1217.0   1.00  \n",
       "beta_pop5     1198.0   1.00  \n",
       "beta_elev2    1376.0   1.00  \n",
       "beta_elev3    1277.0   1.00  \n",
       "sigma         1453.0   1.00  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_2w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47e09c2-ec3b-42de-9ae0-160dc0855e05",
   "metadata": {},
   "source": [
    "#### 10.5.2 Interaction-Effects ANOVA Using PyMC3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "fa497c6e-a545-4359-928e-1c8067500179",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, mean]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='4000' class='' max='4000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [4000/4000 00:08<00:00 Sampling 2 chains, 0 divergences]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 23 seconds.\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model_2w2anova:\n",
    "    mean = pm.Normal('mean', 0, tau = 0.001, shape = [groups_pop, groups_elev])\n",
    "    sigma = pm.Uniform('sigma', 0, 100)\n",
    "    \n",
    "    wing_new = pm.Normal('wing_p', mu = mean[idx_pop, idx_elev], sigma = sigma, observed = wing)\n",
    "    trace_2w2anova = pm.sample(return_inferencedata = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "39c85831-b1a2-49e1-a941-f6d45c27f324",
   "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>mean[0, 0]</th>\n",
       "      <td>43.989</td>\n",
       "      <td>1.554</td>\n",
       "      <td>41.077</td>\n",
       "      <td>46.757</td>\n",
       "      <td>0.024</td>\n",
       "      <td>0.017</td>\n",
       "      <td>4171.0</td>\n",
       "      <td>1274.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[0, 1]</th>\n",
       "      <td>42.941</td>\n",
       "      <td>1.602</td>\n",
       "      <td>39.985</td>\n",
       "      <td>46.100</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.018</td>\n",
       "      <td>3775.0</td>\n",
       "      <td>1554.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[0, 2]</th>\n",
       "      <td>50.154</td>\n",
       "      <td>1.560</td>\n",
       "      <td>47.263</td>\n",
       "      <td>53.085</td>\n",
       "      <td>0.028</td>\n",
       "      <td>0.020</td>\n",
       "      <td>3129.0</td>\n",
       "      <td>1381.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[1, 0]</th>\n",
       "      <td>29.387</td>\n",
       "      <td>1.556</td>\n",
       "      <td>26.588</td>\n",
       "      <td>32.353</td>\n",
       "      <td>0.027</td>\n",
       "      <td>0.019</td>\n",
       "      <td>3365.0</td>\n",
       "      <td>1511.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[1, 1]</th>\n",
       "      <td>32.701</td>\n",
       "      <td>1.514</td>\n",
       "      <td>29.758</td>\n",
       "      <td>35.420</td>\n",
       "      <td>0.024</td>\n",
       "      <td>0.017</td>\n",
       "      <td>3796.0</td>\n",
       "      <td>1399.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[1, 2]</th>\n",
       "      <td>42.776</td>\n",
       "      <td>1.452</td>\n",
       "      <td>39.895</td>\n",
       "      <td>45.399</td>\n",
       "      <td>0.024</td>\n",
       "      <td>0.017</td>\n",
       "      <td>3788.0</td>\n",
       "      <td>1393.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[2, 0]</th>\n",
       "      <td>35.968</td>\n",
       "      <td>1.580</td>\n",
       "      <td>32.998</td>\n",
       "      <td>38.941</td>\n",
       "      <td>0.024</td>\n",
       "      <td>0.017</td>\n",
       "      <td>4138.0</td>\n",
       "      <td>1709.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[2, 1]</th>\n",
       "      <td>44.582</td>\n",
       "      <td>1.480</td>\n",
       "      <td>41.994</td>\n",
       "      <td>47.359</td>\n",
       "      <td>0.025</td>\n",
       "      <td>0.017</td>\n",
       "      <td>3706.0</td>\n",
       "      <td>1203.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[2, 2]</th>\n",
       "      <td>42.507</td>\n",
       "      <td>1.510</td>\n",
       "      <td>39.843</td>\n",
       "      <td>45.332</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.018</td>\n",
       "      <td>3424.0</td>\n",
       "      <td>1646.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[3, 0]</th>\n",
       "      <td>46.535</td>\n",
       "      <td>1.553</td>\n",
       "      <td>43.564</td>\n",
       "      <td>49.384</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.018</td>\n",
       "      <td>3513.0</td>\n",
       "      <td>1611.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[3, 1]</th>\n",
       "      <td>48.738</td>\n",
       "      <td>1.509</td>\n",
       "      <td>45.892</td>\n",
       "      <td>51.555</td>\n",
       "      <td>0.024</td>\n",
       "      <td>0.017</td>\n",
       "      <td>3905.0</td>\n",
       "      <td>1416.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[3, 2]</th>\n",
       "      <td>58.198</td>\n",
       "      <td>1.481</td>\n",
       "      <td>55.332</td>\n",
       "      <td>60.879</td>\n",
       "      <td>0.023</td>\n",
       "      <td>0.016</td>\n",
       "      <td>4194.0</td>\n",
       "      <td>1668.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[4, 0]</th>\n",
       "      <td>51.891</td>\n",
       "      <td>1.573</td>\n",
       "      <td>48.915</td>\n",
       "      <td>54.843</td>\n",
       "      <td>0.023</td>\n",
       "      <td>0.016</td>\n",
       "      <td>4895.0</td>\n",
       "      <td>1531.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[4, 1]</th>\n",
       "      <td>57.530</td>\n",
       "      <td>1.526</td>\n",
       "      <td>54.691</td>\n",
       "      <td>60.496</td>\n",
       "      <td>0.025</td>\n",
       "      <td>0.018</td>\n",
       "      <td>3626.0</td>\n",
       "      <td>1356.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean[4, 2]</th>\n",
       "      <td>58.080</td>\n",
       "      <td>1.468</td>\n",
       "      <td>55.467</td>\n",
       "      <td>61.096</td>\n",
       "      <td>0.026</td>\n",
       "      <td>0.018</td>\n",
       "      <td>3194.0</td>\n",
       "      <td>1160.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sigma</th>\n",
       "      <td>3.057</td>\n",
       "      <td>0.338</td>\n",
       "      <td>2.448</td>\n",
       "      <td>3.684</td>\n",
       "      <td>0.009</td>\n",
       "      <td>0.006</td>\n",
       "      <td>1553.0</td>\n",
       "      <td>1458.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",
       "mean[0, 0]  43.989  1.554  41.077   46.757      0.024    0.017    4171.0   \n",
       "mean[0, 1]  42.941  1.602  39.985   46.100      0.026    0.018    3775.0   \n",
       "mean[0, 2]  50.154  1.560  47.263   53.085      0.028    0.020    3129.0   \n",
       "mean[1, 0]  29.387  1.556  26.588   32.353      0.027    0.019    3365.0   \n",
       "mean[1, 1]  32.701  1.514  29.758   35.420      0.024    0.017    3796.0   \n",
       "mean[1, 2]  42.776  1.452  39.895   45.399      0.024    0.017    3788.0   \n",
       "mean[2, 0]  35.968  1.580  32.998   38.941      0.024    0.017    4138.0   \n",
       "mean[2, 1]  44.582  1.480  41.994   47.359      0.025    0.017    3706.0   \n",
       "mean[2, 2]  42.507  1.510  39.843   45.332      0.026    0.018    3424.0   \n",
       "mean[3, 0]  46.535  1.553  43.564   49.384      0.026    0.018    3513.0   \n",
       "mean[3, 1]  48.738  1.509  45.892   51.555      0.024    0.017    3905.0   \n",
       "mean[3, 2]  58.198  1.481  55.332   60.879      0.023    0.016    4194.0   \n",
       "mean[4, 0]  51.891  1.573  48.915   54.843      0.023    0.016    4895.0   \n",
       "mean[4, 1]  57.530  1.526  54.691   60.496      0.025    0.018    3626.0   \n",
       "mean[4, 2]  58.080  1.468  55.467   61.096      0.026    0.018    3194.0   \n",
       "sigma        3.057  0.338   2.448    3.684      0.009    0.006    1553.0   \n",
       "\n",
       "            ess_tail  r_hat  \n",
       "mean[0, 0]    1274.0    1.0  \n",
       "mean[0, 1]    1554.0    1.0  \n",
       "mean[0, 2]    1381.0    1.0  \n",
       "mean[1, 0]    1511.0    1.0  \n",
       "mean[1, 1]    1399.0    1.0  \n",
       "mean[1, 2]    1393.0    1.0  \n",
       "mean[2, 0]    1709.0    1.0  \n",
       "mean[2, 1]    1203.0    1.0  \n",
       "mean[2, 2]    1646.0    1.0  \n",
       "mean[3, 0]    1611.0    1.0  \n",
       "mean[3, 1]    1416.0    1.0  \n",
       "mean[3, 2]    1668.0    1.0  \n",
       "mean[4, 0]    1531.0    1.0  \n",
       "mean[4, 1]    1356.0    1.0  \n",
       "mean[4, 2]    1160.0    1.0  \n",
       "sigma         1458.0    1.0  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace_2w2anova)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2e2c575-ead7-418f-99cc-432740cf6bd3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}