diff --git a/Rethinking/Chp_07.ipynb b/Rethinking/Chp_07.ipynb
new file mode 100644
index 0000000..ccb72e3
--- /dev/null
+++ b/Rethinking/Chp_07.ipynb
@@ -0,0 +1,1923 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import arviz as az\n",
+    "import bambi as bmb\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import scipy.stats as stats\n",
+    "import statsmodels.formula.api as smf\n",
+    "\n",
+    "from scipy.special import logsumexp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext watermark\n",
+    "az.style.use(\"arviz-darkgrid\")\n",
+    "np.random.seed(1211)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "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>species</th>\n",
+       "      <th>brain</th>\n",
+       "      <th>mass</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>afarensis</td>\n",
+       "      <td>438</td>\n",
+       "      <td>37.0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>africanus</td>\n",
+       "      <td>452</td>\n",
+       "      <td>35.5</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>habilis</td>\n",
+       "      <td>612</td>\n",
+       "      <td>34.5</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>boisei</td>\n",
+       "      <td>521</td>\n",
+       "      <td>41.5</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>rudolfensis</td>\n",
+       "      <td>752</td>\n",
+       "      <td>55.5</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>5</th>\n",
+       "      <td>ergaster</td>\n",
+       "      <td>871</td>\n",
+       "      <td>61.0</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>6</th>\n",
+       "      <td>sapiens</td>\n",
+       "      <td>1350</td>\n",
+       "      <td>53.5</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "       species  brain  mass\n",
+       "0    afarensis    438  37.0\n",
+       "1    africanus    452  35.5\n",
+       "2      habilis    612  34.5\n",
+       "3       boisei    521  41.5\n",
+       "4  rudolfensis    752  55.5\n",
+       "5     ergaster    871  61.0\n",
+       "6      sapiens   1350  53.5"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sppnames = [\"afarensis\", \"africanus\", \"habilis\", \"boisei\", \"rudolfensis\", \"ergaster\", \"sapiens\"]\n",
+    "brainvolcc = [438, 452, 612, 521, 752, 871, 1350]\n",
+    "masskg =  [37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5]\n",
+    "\n",
+    "d = pd.DataFrame.from_dict({\"species\": sppnames, \"brain\": brainvolcc, \"mass\": masskg})\n",
+    "d"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d[\"mass_std\"] = (d[\"mass\"] - d[\"mass\"].mean()) / d[\"mass\"].std()\n",
+    "d[\"brain_std\"] = d[\"brain\"] / d[\"brain\"].max()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.3"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "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: [brain_std_sigma, mass_std, Intercept]\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:01<00:00 Sampling 2 chains, 1 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 2 seconds.\n",
+      "There was 1 divergence after tuning. Increase `target_accept` or reparameterize.\n"
+     ]
+    }
+   ],
+   "source": [
+    "priors = {\n",
+    "    \"Intercept\": bmb.Prior(\"Normal\", mu=0.5, sd=1),\n",
+    "    \"mass_std\": bmb.Prior(\"Normal\", mu=0, sd=10),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}\n",
+    "\n",
+    "model_7_1 = bmb.Model(\"brain_std ~ mass_std\", d, priors=priors)\n",
+    "results_7_1 = model_7_1.fit()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_7_1_OLS = smf.ols(\"brain_std ~ mass_std\", data=d).fit()\n",
+    "mean = model_7_1_OLS.params.values\n",
+    "cov = np.diag(model_7_1_OLS.cov_params())\n",
+    "post = stats.multivariate_normal.rvs(mean=mean, cov=cov, size=10000)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def var2(x):\n",
+    "    return np.mean(x ** 2) - np.mean(x) ** 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/pymc3/sampling.py:1689: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample\n",
+      "  warnings.warn(\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='1000' class='' max='1000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
+       "      100.00% [1000/1000 00:01<00:00]\n",
+       "    </div>\n",
+       "    "
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "model_7_1.posterior_predictive(results_7_1, 1000)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.5062199789556512"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "r = results_7_1.posterior_predictive[\"brain_std\"].values.reshape(-1, 7).mean(axis=0) - d[\"brain_std\"].values\n",
+    "resid_var = var2(r)\n",
+    "outcome_var = var2(d[\"brain_std\"])\n",
+    "1 - resid_var / outcome_var"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.6"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def R2_is_bad(model, result):\n",
+    "    model.posterior_predictive(result, 1000)\n",
+    "    r = result.posterior_predictive[\"brain_std\"].values.reshape(-1, 7).mean(axis=0) - d[\"brain_std\"].values\n",
+    "    return 1 - var2(r) / var2(d[\"brain_std\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.7\n",
+    "\n",
+    "Here we use the `common` key in the priors dictionary. This means that all the common predictors (fixed effects) are going to receive the same prior, a normal with mean 0 and sd of 10, as our inspection below confirms. We also start with a baseline formula, and add terms to the formula as we create models."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "formula = \"brain_std ~ mass_std + I(mass_std ** 2)\"\n",
+    "priors = {\n",
+    "    \"Intercept\": bmb.Prior(\"Normal\", mu=0.5, sd=1),\n",
+    "    \"common\": bmb.Prior(\"Normal\", mu=0, sd=10),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Formula: brain_std ~ mass_std + I(mass_std ** 2)\n",
+       "Family name: Gaussian\n",
+       "Link: identity\n",
+       "Observations: 7\n",
+       "Priors:\n",
+       "  Intercept ~ Normal(mu: 0.5, sd: 1)\n",
+       "  mass_std ~ Normal(mu: 0, sd: 10)\n",
+       "  I(mass_std ** 2) ~ Normal(mu: 0, sd: 10)\n",
+       "  sigma ~ Exponential(lam: 1)"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model_7_2 = bmb.Model(formula, d, priors=priors)\n",
+    "model_7_2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "formula += \" + I(mass_std ** 3)\"\n",
+    "model_7_3 = bmb.Model(formula, d, priors=priors)\n",
+    "formula += \" + I(mass_std ** 4)\"\n",
+    "model_7_4 = bmb.Model(formula, d, priors=priors)\n",
+    "formula += \" + I(mass_std ** 5)\"\n",
+    "model_7_5 = bmb.Model(formula, d, priors=priors)\n",
+    "formula += \" + I(mass_std ** 6)\"\n",
+    "model_7_6 = bmb.Model(formula, d, priors=priors)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "I'm skipping this for now, until we have better Laplace method."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.12"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.6108643020548935"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "p = np.array([0.3, 0.7])\n",
+    "-np.sum(p * np.log(p))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABvFklEQVR4nO3dd3gU1f7H8fekkWyAkCA1CUUlAaVJkyogoEgTsSE2RGw/QMWKDbmWi1ev3usVe0MREVBRQVSa0pESiiBFQSCEKimQQur5/TEkgrSw2c22z+t5eHZ2ZjL5hpPsfnJy5hzLGGMQERERERGXC/J0ASIiIiIi/kphW0RERETETRS2RURERETcRGFbRERERMRNFLZFRERERNxEYVtERERExE0UtkVERERE3ERhW0RERETETUI8XYC3S0tLc9u1o6KiyMjIcNv1xfMCpY2zsrKIj48HIDk5mcjISA9XVD4CpX0DmdrY/6mN/Zu72zc6OvqM56hn24OCgvTf7+/Uxv5N7ev/1Mb+T23s37yhfT1fgYiIiIiIn9IwEhEps7CwMF588cWSbREREbEpbItImYWGhjJ06FBPlyEiIuJ1NIxERERERMRN1LMtImVWWFjI0qVLAWjXrh3BwcEerkhERMQ7KGyLSJkdOXKEfv36AYE19Z+IiMiZaBiJiIiIiIibKGyLiIiIiLiJwraIiIiIiJsobIuIiIiIuInCtoiIiIiImyhsi4iIiIi4iab+E5EyCw0NZcyYMSXbIiIiYlPYFpEyCwsL49577/V0GSIiIl5Hw0hERERERNxEPdsiUmaFhYWsXbsWgGbNmmm5dhERkaMUtkWkzI4cOUL37t0BLdcuIiLla/0Gw7790K2r5elSTkphW0RERER81ugxhv0H4KJmEBPjfYFbY7ZFRERExCfl59tBG8BbRzAqbIuIiIiIT8rIsB+Dg6BSJc/WcioK2yIiIiLik9LS7ceoKAgK8r4hJKCwLSIiIiI+Kj3dfoyO9mgZp6WwLSIiIiI+KS3NfqxSxaNlnJZmIxGRMgsNDeWRRx4p2RYRESkPxcNIvLlnW2FbRMosLCyMUaNGeboMEREJMGlpBoDoKp6t43Q0jEREREREfFJxz3aVKt55cySoZ1tEXKCoqIjNmzcDkJiYSFCQfo8XERH384UbJBW2RaTMcnJy6NChA6Dl2kVEpPwU3yCpYSQiIiIiIi6mGySdtG/fPr777jsWLFjAtm3b+PPPP4mKiqJFixYMHTqUZs2aleo6P//8M7fccsspj48dO5YBAwa4qmwRERERKUfpPtCz7ZVhe8KECbz77rvUqVOHDh06EBMTw44dO5gzZw5z5szh5ZdfplevXqW+Xps2bWjTps0J+xs1auTKskVERESknOTkGHKO2NuaZ/ssNW3alAkTJpwQkFeuXMngwYMZM2YM3bt3JywsrFTXa9OmDSNGjHBHqSIiIiLiAcU3R4aFgsPh0VJOyyvHbF922WUn7Ylu1aoVF198MRkZGSUzH4iIiIhI4EnPsB+rRINlaeo/lwkJCTnusTS2b9/O+PHjyc3NpUaNGrRr144aNWq4q0QRERERcTNfmIkEfCxs7969myVLllCtWjUSEhJK/XEzZsxgxowZJc9DQkK46aabeOSRRwgODnZHqSIBJTQ0lOHDh5dsi4iIuJsvzEQCPhS28/PzeeSRR8jLy+Ohhx4qVUiOiYnhwQcfpGvXrsTGxpKTk8Pq1at5+eWXGT9+PJZlnXGJ6aioKLcu0BHt7d8hUmaB0savvfaap0vwiEBp30CmNvZ/amPfdORIDpBN9eoViI6ueMrzPN2+ljHGeLSCUigqKuLhhx9mxowZXHfddTz77LNlut6BAwfo168fhw4dYsGCBVStWvWU56YV/43CDaKjo916ffE8tbF/U/v6P7Wx/1Mb+67X3ihi8hS44XoYds/JO0bd3b6lCfJeeYPksYqKinj88ceZMWMG/fr14x//+EeZr1mtWjW6detGQUEBa9eudUGVIoGtqKiInTt3snPnToqKijxdjoiIBICSObajvffmSPDyYSRFRUU89thjfPXVV/Tp04cXXnjBZUM6in8TycnJccn1RAJZTk4OzZs3B7Rcu4iIlI+SMdtVPFnFmXltz/axQbtXr168+OKLLr2ZsbhHOzY21mXXFBEREZHyUTzPdhUvH3LvlWG7eOjIV199Rc+ePXnppZdOG7RTU1PZunUrqampx+1fv379Sc//6KOP+Pnnn6lXrx5NmjRxae0iIiIi4n6a+q8MXn/9daZNm4bD4aBevXq8+eabJ5zTvXv3kuXWJ06cyLhx4xg+fPhxK0Xee++9hISE0LhxY2rUqEFOTg5r167l119/pXLlymcM8SIiIiLifYwxJcNIvHmpdvDSsJ2SkgJAdnY2b7311knPiY2NLQnbpzJw4EAWLVrEihUrSE9PJygoiNq1a3PrrbcyZMgQatas6fLaRURERMS9MjOhoMDe9vaebZ+Y+s+TNPWflEWgtHFWVhbx8fFAYN0gGSjtG8jUxv5PbeybdiYbBt1scDhg1sxTj4rW1H8iIiIiImep+OZIb+/VBi8dRiIiviUkJITbb7+9ZFtERMSdfGWpdlDYFhEXqFChAi+99JKnyxARkQBRPDLE22+OBA0jEREREREfo2EkIhJQjDEcPHgQgKpVq2JZ3r10roiI+La0NHt+D29f0AYUtkXEBbKzs0lISAACazYSERHxjL+Wavf+zh0NIxERERERn1IyjMQHerYVtkVERETEp5TcIBnl2TpKQ2FbRERERHyKL039p7AtIiIiIj6jsNCQkWFv+8JsJArbIiIiIuIzDh0CY09GQpSGkYiIiIiIuE7xeO2oyhAS4v2zkWjqPxEps5CQEG644YaSbREREXdJLx5C4gPjtUFhW0RcoEKFCrz++uueLkNERAKALy3VDhpGIiIiIiI+pHgmEl8J2+rZFpEyM8aQnZ0NgMPh0HLtIiLiNsVLtfvKMBL1bItImWVnZxMfH098fHxJ6BYREXEHX1qqHRS2RURERMSH+NJS7aCwLSIiIiI+xJeWageFbRERERHxIb60VDsobIuIiIiID0k/2rPtC0u1g8K2iIiIiPiIvDxDZpa9rZ5tEREREREXKr45MjgYKlb0aCmlpnm2RaTMgoOD6devX8m2iIiIOxQv1V4lCoKCfGPqP4VtESmz8PBwxo8f7+kyRETEz5XMROIjQ0hAw0hERERExEf8taCNJ6s4OwrbIiIiIuITinu2feXmSFDYFhEXyMrKIiYmhpiYGLKysjxdjoiI+Km0dAOoZ1tERERExOVK5tiO9o2bI0FhW0RERER8RPHUf76yVDsobIuIiIiIj/C1pdpBYVtEREREfETJ1H9VPFrGWVHYFhERERGvZ4xRz7aIiIiIiDvk5EBurr3tS7ORaAVJESmz4OBgevToUbItIiLiasW92mFhEBHh0VLOisK2iJRZeHg4kydP9nQZIiLix4pnIomOBsvS1H8iIiIiIi7ji0u1g8K2iIiIiPiAdB9cqh0UtkXEBbKysoiLiyMuLk7LtYuIiFv4as+2xmyLiEtkZ2d7ugQREfFjaWkG8K05tkE92yIiIiLiA4p7tqtU8Z2bI0FhW0RERER8wLGzkfgShW0RERER8XrFS7X72phthW0RERER8Xq+uFQ7KGyLiIiIiJfLzzclU/9VjfFsLWdLs5GISJkFBQXRoUOHkm0RERFX2rMHCosgPByqVvV0NWdHYVtEyiwiIoLp06d7ugwREfFTybvsx7hY31qqHVwQtvft28eqVavYv38/ANWrV6dFixbUrFmzzMWJiIiIiJSE7TjP1uEMp8P2vn37ePbZZ5k3bx7GmOOOWZZF165deeqppxS6RURERKRMdu2ys2Z8vIcLcYJTYXvfvn0MHDiQPXv2EBERQYcOHYiNjQVg9+7dLFq0iLlz57JhwwYmT55MjRo1XFq0iHiXrKwsmjdvDsCaNWuIjIz0bEEiIuJXinu242N9awgJOBm2//Of/7Bnzx769u3LE088QZW/rZuZkZHBP//5T77++mv++9//MnbsWFfUKiJe7ODBg54uQURE/FRJ2PbBnm2npg1YsGABcXFxvPDCCycEbYCoqCj++c9/EhcXx08//VTGEkVEREQkUB05Yjh6a6BPjtl2KmxnZ2fTrFkzgoODT3lOcHAwzZo1Iycnx+niRERERCSwpey2HytWhCpRnq3FGU6F7XPPPbdk9pHT2b9/P+eee64zn0JEREREhF3FQ0jifG/aP3AybN96662sXLmShQsXnvKcRYsWsXLlSm655RanixMRERGRwLYz2X6M98EhJODkDZKtW7dm0KBB3HPPPfTq1YtevXpRu3ZtwJ6N5LvvvmPmzJnceOONtGnTht27dx/38cXnioiIiIiczq4Ue9q/uDjf69UGJ8P2pZdeimVZGGOYPn36SVeOM8YwceJEJk6ceNx+y7L49ddfnatWRLxSUFAQF110Ucm2iIiIqyQX92z74EwkUIaebRGRYhEREcydO9fTZYiIiB8qGbMd69k6nOVU2J4wYYKr6zjOvn37+O6771iwYAHbtm3jzz//JCoqihYtWjB06FCaNWtW6msVFRUxceJEpkyZwo4dO3A4HLRv356RI0cS76u/IomIiIgEgKwsQ2qave2L0/6BkzdIutuECRMYO3YsycnJdOjQgdtuu42WLVsyd+5cBg4cyMyZM0t9rdGjR/Pcc89hjOHmm2+mU6dOzJo1i2uuuYbt27e774sQERERkTIpXswmOhoqVgygMdvu1rRpUyZMmECbNm2O279y5UoGDx7MmDFj6N69O2FhYae9zrJly5g6dSqtW7fmgw8+KDm/T58+3HnnnTz77LO8//77bvs6RAJFdnY27dq1A2Dp0qU4HA4PVyQiIv4g+Zhp/3yVV/ZsX3bZZScEbYBWrVpx8cUXk5GRwebNm894nalTpwJw3333HRfMO3fuTJs2bVi0aNEJM6WIyNkzxpCcnExycjLGGE+XIyIifqJ4vLavDiEBLw3bpxMSEnLc4+n8/PPPOBwOWrRoccKxTp06AbB8+XLXFigiIiIiLpGcbHfgxPvotH/gY2F79+7dLFmyhGrVqpGQkHDac7Ozszlw4ABxcXEnXVa+bt26AOzYscMttYqIiIhI2SSn2I++PIzEK8dsn0x+fj6PPPIIeXl5PPTQQycN0Mc6fPgwABUrVjzp8eL9xeedSlRUlFvnDY6OjnbbtcU7BEIbHztMKzo6msjISA9WU74CoX0DndrY/6mNvZMxhpRdaYDhwgsrEx3tXGz1dPv6RNguKipi1KhRrFixguuuu47+/fuX2+fOyMhw27Wjo6NJS0tz2/XF8wKljbOyskq209LSyMvL82A15SdQ2jeQqY39n9rYe6WnGw4dtoeRVKp4iLS0sx9K4u72LU2Q9/phJEVFRTz++OPMmDGDfv368Y9//KNUH1epUiUAMjMzT3q8eH/xeSIiIiLiPXYdHUJSvTqEh/vumG2v7tkuKiriscce46uvvqJPnz688MILpR7S4XA4qFatGrt27aKwsPCEYSfFY7WLx26LiPMsyyIxMbFkW0REpKxKlmn34fHaUIaw/dtvv/HBBx+wfPlyDhw4QH5+/knPsyyLX3/99ayvf2zQ7tWrFy+++OIZx2n/XZs2bfj2229JSko6YYn5hQsXAlp6XsQVHA4HS5cu9XQZIiLiR5J3Fc9E4uFCysipsL18+XLuuOMOcnNzsSyLqKgoly5iUTx05KuvvqJnz5689NJLpw3aqamppKWlER0dTUxMTMn+6667jm+//ZZXX331uEVt5s+fz/Lly+nYsSOxsbEuq1tEREREXCO5ZI5t3/6LqVNh+6WXXiI3N5d77rmH22+//ZQzfjjr9ddfZ9q0aTgcDurVq8ebb755wjndu3enUaNGAEycOJFx48YxfPhwRowYUXJO27Ztufbaa5k6dSoDBgygc+fOHDhwgJkzZ1KlShWefPJJl9YtIiIiIq6xyw9WjwQnw/bmzZtp3rw59913n6vrASAlxR4Rn52dzVtvvXXSc2JjY0vC9uk888wzJCQkMGXKFD7++GMcDgc9evRg5MiR1KlTx6V1iwSq7OxsunXrBsDcuXO1XLuIiJSJMcZvwrZlnFhb+ZJLLqF169a8/PLL7qjJq7h7uhhNN+TfAqWNs7KyiI+PByA5OTlg5tkOlPYNZGpj/6c29k5//mnof40hOAjm/GARGurcUBKfnfqvS5curFmzhsLCQmc+XERERETklHYenYmkZk2cDtrewqmwff/992NZFk888cQZV2AUERERETkbxXNsH/2jqU9zasx2TEwMU6dO5eabb+bSSy+lcePG1KhR46Tz61qWxT//+c8yFyoiIiIigSE52R7lHOfj47XBybCdmZnJfffdx9atWzHGnHZ+XYVtERERETkbf90c6dtDSMDJsP2vf/2L5cuX06BBA6677jri4+M1+4CIiIiIuESyn8xEAk6G7blz51KrVi0mT56skC0iWJZVMhuJlmsXEZGyKCw0pOy2twM2bB85coQ2bdooaIsIYC/XvnbtWk+XISIifmD/fsjPh9BQqF7d09WUnVOzkTRq1Ig///zT1bWIiIiISIArHkISGwvBwb7/11KnwvawYcNYvXo1CxYscHU9IiIiIhLASsZrx3q2DldxahhJaGgogwYN4p577qFv3760b9+eGjVqEBR08uzeunXrMhUpIt4tJyeHPn36ADBjxgwiIiI8XJGIiPiq5F32tH/+MMc2OBm2b775ZizLwhjDV199xddff33a8zdu3OhUcSLiG4qKili9enXJtoiIiLOKp/2L84Np/8DJsN2/f3/NOCAiIiIiLpd8dKl2f5iJBJwM2y+88IKr6xARERGRAJeba9i71972l2EkTt0gKSIiIiLialt+g8IiiImGqjGersY1nOrZPlZeXh6bNm1i3759ANSoUYOGDRsSFhZW5uJEREREJHBs3GQ/NmzoP4ukOR22c3NzefXVV5k8eTLZ2dnHHXM4HAwcOJB7772XChUqlLlIEREREfF/mzbZM5Fc0Mg/gjY4Gbbz8vIYPHgwa9asASAxMZHY2FgsyyIlJYVNmzbxwQcfkJSUxEcffaRebpEAULVqVU+XICIiPm7jZvuxYaJn63Alp8L2+PHjWb16NS1btmT06NEkJh7/P7JlyxaeffZZVq5cyfjx47nzzjtdUqyIeKfIyEh+++03T5chIiI+7PBhUzITiT+FbadukJwxYwYxMTG8/fbbJwRtgISEBN566y2io6OZPn16mYsUEREREf+2eYv9WKsWVKniP8NInArbO3fupE2bNlSsWPGU50RGRtKmTRuSi39FERERERE5heKbIxs19GwdruZU2A4ODubIkSNnPO/IkSMEBwc78ylExIfk5OTQt29f+vbtS05OjqfLERERH7Rps31zZMNE/+nVBifHbCckJLBs2TKSk5OJP8WM48nJySxbtowLLrigTAWKiPcrKipi8eLFJdsiIiJnSz3bx7j++us5cuQIN998M1OnTj2ul/vIkSN88cUX3HLLLeTm5jJw4ECXFSsiIiIi/ic11bB/P1gWJCZ4uhrXcqpnu3///iQlJTFlyhRGjx7N6NGjiY6OBiAtLQ0AYwzXX389/fr1c121IiIiIuJ3inu169YFh0PDSAB45pln6NChAxMmTGDt2rWkpqYCEBoaSvPmzbnpppu4/PLLXVaoiIiIiPin4vHajfxoyr9iToXtTZs2ERQUxOWXX87ll19OQUEB6enpAFSpUoWQkDKvAi8iIiIiAeKvZdr9q1cbyjCMpHXr1kyYMMG+SEgI55xzjksLExERERH/Z4xhk5/eHAlOhu2oqCiqV6/u6lpExIc5HA5PlyAiIj5oz15Iz4CQEDj/PE9X43pOhe3mzZuzZcsWV9ciIj4qMjKSXbt2eboMERHxQcVDSM47F8LC/G8YiVNT/w0bNow//viDDz74wNX1iIiIiEgA2bTp6M2RfjiEBJzs2d62bRv9+vXjpZde4ptvvqFz587Url2bChUqnPT8/v37l6VGEREREfFT/nxzJDgZtkeNGoVlWUcHtG9i06ZNWNaJ/0HGGCzLUtgW8XNHjhzh1ltvBeCjjz4iPDzcwxWJiIgvKCw0bD46Mlk928cYNmzYScO1iASmwsJCZs+eXbItIiJSGjuTIScHIsKhXl1PV+MeToXtESNGuLoOEREREQkwxVP+JSRAcLB/duSW6gbJRo0a8fjjj5c8HzduHHPnznVbUSIiIiLi/zYevTmyoZ8OIYFShm1jDMaYkufjxo1jzpw5bitKRERERPzfxs32Y6NE/+zVhlKGbYfDQWpqqrtrEREREZEAkZ9v+P13e9tfb46EUo7ZTkxMZMmSJYwbN464uDgAdu7cyVdffVWqT6LZSERERETkWFu3QX4+VK4MtWt7uhr3KVXYHj58OMOHD2fcuHEls5AkJSWRlJR02o/T1H8iIiIicjLFN0c2TMSvZ7krVdju0KEDM2fOZMmSJezZs4dx48bRsGFDunXr5u76RMQHREZGaqiZiIiclY2bj94cmejhQtys1FP/1apVi6uvvhqgJGwPHz7cbYWJiIiIiP8q7tlu5KcrRxZzap7tjz/+mHPOOcfVtYiIiIhIAMjKMvyx3d7255sjwcmw3aZNG1fXISI+7MiRI9x9990AvPXWW1quXURETitpNRQVQVwcnHOOf/dsl2rqPxGR0yksLOSbb77hm2++0XLtIiJyRitX2eO1W7fycCHlQGFbRERERMrVipX2Y6uW/t2rDQrbIiIiIlKO9u037EyGoCBo0dzT1bifwraIiIiIlJviXu1GDaFSJfVsi4iIiIi4TCCN1waFbREREREpJ0VFhpWr7O1AGK8NCtsiIiIiUk62boX0dIgIhwsv8HQ15cOpebZvueUWpz+hZVl89NFHTn+8iHgfh8NBcnJyybaIiMjJrDjaq31RcwgNDYyebafC9vLlywE7OBtjTnrOqY5ZVmD8x4oEEsuyiIyM9HQZIiLi5YrHa7dqFTh50KmwPXfuXD788EMmTZrE5ZdfTu/evYmNjQVg9+7dfPvtt/zwww8MHDiQ2267zaUFi4iIiIjvyc01rFlrb7dq6dlaypNTYXvVqlV8+umnvPHGG3Tp0uW4Yw0bNuTSSy+lX79+3HPPPTRp0oQrr7zSFbWKiJfKzc3lgQceAOCVV16hQoUKHq5IRES8zS/rIS8PzjkH6tfzdDXlx6kbJMePH0/Lli1PCNrH6ty5My1bttT4bJEAUFBQwKRJk5g0aRIFBQWeLkdERLzQipVHh5C0DKxhxU6F7W3btlG9evUznle9enW2bdvmzKcQERERET8SaFP+FXMqbEdERPDLL79QVFR0ynOKior45ZdfiIiIcLo4EREREfF9GRmGLb/Z24E0XhucDNsdO3YkOTmZp59+mszMzBOOZ2VlMWbMGJKTk+nYsWOZixQRERER37UyCYyBc+vDOVUDq2fbqRskH3zwQZYuXcrnn3/O999/T6dOnahduzZgz0aycOFCMjMzqVq1aslNU2fr66+/ZtWqVaxfv54tW7aQn5/P2LFjGTBgQKmv8fPPP592TvCzvZ6IiIiInL2VKwNrifZjORW2a9asyWeffcaYMWNYtGgRM2fOPOGcDh06MGbMGGrVquVUYa+++iopKSlER0dTvXp1UlJSnLoOQJs2bWjTps0J+xs1auT0NUVERETkzIwxrFhpbwfaeG1wMmwDxMXF8d5775GcnMyqVavYv38/YN8U2bJlS+Lj48tU2HPPPUfdunWJjY3lnXfe4eWXX3b6Wm3atGHEiBFlqkdEREREzl5KCuzdByEh0LyZp6spf06H7WLx8fFlDtYn0759e5dfU0Tcw+FwsGXLlpJtERGRYsVLtDe+ECIi1LPtl7Zv38748ePJzc2lRo0atGvXjho1ani6LBG/YVkW55xzjqfLEBERL1S8RHvrAFqi/VgBEbZnzJjBjBkzSp6HhIRw00038cgjjxAcHOzBykRERET8V0GBYVWSvR1oU/4V8+uwHRMTw4MPPkjXrl2JjY0lJyeH1atX8/LLLzN+/Hgsy2LUqFGnvUZUVBRBQU7NkFgq0dHRbru2eIdAaONAXq49ENo30KmN/Z/a2H2W/ZxPZuYhqlSxaNc2muDg8u/d9nT7+nXYbtCgAQ0aNCh57nA46N69O82aNaNfv35MmDCBO+64g6pVq57yGhkZGW6rLzo6mrS0NLddXzwvUNo4KyuLN954A4DHHnuMyMhID1dUPgKlfQOZ2tj/qY3da/oMewHETh0Mhw6ll/vnd3f7libIu6/L1otVq1aNbt26UVBQwNq1az1djoiIiIjfKSw0zF9ob3ftEpjjtaGUYXv37t2kp6e7uZTyVfybSE5OjocrEREREfE/a9dBWhpUrgwtLvJ0NZ5TqrDdrVs3XnzxxZLnjz32GJ9//rnbiioPxT3asbGxHq5ERERExP/8ON+ehaRTRwgJUc/2aRljMMaUPJ82bRqrVq1yW1FnKzU1la1bt5Kamnrc/vXr15/0/I8++oiff/6ZevXq0aRJk/IoUURERCRgFBYaFiywt7t0DtygDaW8QbJSpUrs2bPH3bUcZ+rUqSWBvnixjKlTp7J8+XIAWrZsybXXXgvAxIkTGTduHMOHDz9upch7772XkJAQGjduTI0aNcjJyWHt2rX8+uuvVK5cmZdeeklT/4mIiIi42C/r4WAqVKwIrVp4uhrPKlXYbtKkCcuWLeOxxx4rGXaxadMmxo0bd8aPtSyLYcOGnXVhq1atYtq0acftS0pKIikpqeR5cdg+lYEDB7Jo0SJWrFhBeno6QUFB1K5dm1tvvZUhQ4ZQs2bNs65LRERERE7vx5+ODiHpAKGhgd2zbZljx4ecwoYNG7j77rs5cODA2X8Cy2Ljxo1OFecN3D1djKYb8m+B0sZFRUXs2rULgLi4OLfOTe9NAqV9A5na2P+pjV2vqMhw1bWGgwfhxbEW7dt5Lmx7w9R/perZvvDCC/n+++/55Zdf2Lt3L6NGjaJly5Zcc801ZS5SRHxfUFAQderU8XQZIiLiBX5ZDwcPQmRk4K4aeaxSL2oTGRlJ27ZtARg1ahR16tThqquuclthIiIiIuJ7fjo6C0nHDhAWFthDSMDJFSTnzp2Lw+FwdS0i4qPy8vJ47rnnAHjyyScJCwvzcEUiIuIJRUWGn+bb210DfBaSYk6F7b/PTX3w4EH27dsHQI0aNU67/LmI+J/8/PySG6YfffRRhW0RkQD160Y48Cc4HNC6laer8Q5Ohe1iEydO5OOPP2bnzp3H7a9bty633HILgwYNKlNxIiIiIuI7imch6dAeKlRQzzY4GbaLioq4//77mT17NsYYKleuTO3atbEsi927d7N9+3aeffZZli1bxquvvopl6T9bRERExJ8ZY/hRQ0hO4FTYnjx5MrNmzaJ+/fo88sgjdO3a9bjjP/30Ey+++CKzZ89m8uTJDBw40CXFioiIiIh3+nUj7N8PERFwcRtPV+M9nJoM98svv6RixYpMmDDhhKAN0KVLFz766CMcDgdffPFFmYsUEREREe9WPAtJ+3YaQnIsp8L277//Ttu2bTnnnHNOeU61atVo164dv//+u9PFiYiIiIj3M+avWUgu7aKgfSynl3krzThsjdUWERER8X/rfoE9eyEiXENI/s6pMdv169dn2bJlpKamEhMTc9JzUlNTWbZsGfXr1y9TgSLi/SIiIli8eHHJtoiIBJZvZthDSLpdCuHh6mw9llM921dddRWHDx9m8ODBLF269ITjy5YtY8iQIWRmZjJgwIAyFyki3i0oKIhGjRrRqFEjgoKc/oOZiIj4oMOHDT/+ZG/37aOg/XdO9WwPGjSIhQsXsmDBAoYMGUJMTAy1a9cGYPfu3aSmpmKMoXPnzpprW0RERMSPzZoDeXlwbn24oJGnq/E+ToXt4OBg3nrrLcaPH8+ECRPYs2cPBw8eLDleu3ZtbrrpJgYPHqxeLpEAkJeXxyuvvALAAw88oBUkRUQChDGG6UeHkPTpbel+vZOwjDGmrBfZs2cP+/fvB6B69erUqlWrzIV5i7S0NLddOzo62q3XF88LlDbOysoiPj4egOTkZCIjIz1cUfkIlPYNZGpj/6c2LptNmwxD7zaEhcJXX1hUruxdYdvd7RsdHX3Gc8q0XHuxWrVq+VXAFhEREZEzm/6t3Wd7ySV4XdD2FhrjISIiIiJnLTvbMHuuvd1PN0aeksK2iIiIiJy1H+dDdjbE1obmzTxdjfdS2BYRERGRs3bsjZFBQerZPhWFbRERERE5K9v+MKzfAMFBcEVPT1fj3RS2RUREROSszDh6Y2T79nBOVfVqn45LZiMRkcAWHh7OnDlzSrZFRMR/5eUZvp9lb/ftraB9JgrbIlJmwcHBtGjRwtNliIhIOViwEA4dgmrnwMVtPF2N9ytT2P7999+ZMmUK69atIy0tjW7duvHII48AkJSUxPr16+nXrx9VqlRxRa0iIiIi4mHFc2v37gXBwerZPhOnw/aHH37Iyy+/TEFBAQCWZZ2wQs/YsWMJCwtj4MCBZatSRLxaXl4eb731FgB33323lmsXEfFTO3caViWBZUHvKxS0S8OpGyR/+ukn/vWvf1GzZk3GjRvHkiVL+Puq7y1atCAmJoa5c+e6pFAR8V75+fmMGTOGMWPGkJ+ff8bzc3Jy+Mc//kHXrl3p3bs3EydO5J577uE///mP22p89913ufnmm5k2bRr9+vWjc+fOPPHEE2RmZrrtc4qI+JvPptp5r11bqFVLYbs0nOrZ/vDDD4mIiODDDz8kPj7+lOc1bNiQP/74w+niRMQ/vfbaa6xevZoXX3yR6Oho3nzzTTZv3kxCQsIpP2bNmjWMHDnytNd99NFH6dnz1HNQ7dq1i7lz5/Lvf/+brKwsnn/+eV588UWeeeYZp78WEZFAkZpq+P57e3vQQAXt0nIqbG/YsIHmzZufNmgDREdHs2rVKqcKExH/lJ2dzfTp0xkzZgytW7cGYPTo0fTr1++0H9ewYUM+/vjj054TExNz2uN5eXmMHj2a6tWrA/Dggw/y4IMPct9991G1atWz+CpERALPl18Z8vKhUSNo1tTT1fgOp8J2fn4+kZGRZzwvNTWV4OBgZz6FiPiplJQU8vPzufDCC0v2RUVFUbdu3dN+XHh4+Bl/wT+TGjVqlARtgCZNmlBUVMSOHTsUtkVETiMnx/DlV/b2oOstLEs926XlVNiOi4tj06ZNpz0nLy+PzZs3U69ePWc+hYjIcVwxjERERJwz8zt7ur/ateGSTp6uxrc4FbYvvfRS3nvvPT788ENuu+22k57z3nvvkZqayi233FKmAkXEv8TGxhISEsKGDRuoWbMmAIcOHWLnzp1cdNFFp/w4Vwwj2bdvHwcOHKBatWoArF+/nqCgoDP2qouIBLLCQlNyY+T111qa7u8sORW2hw4dyvTp03nxxRdZu3YtPXr0AODgwYPMnj2b2bNnM336dOLi4rjxxhtdWrCI+DaHw0Hfvn157bXXiIqKIjo6mrfeeougoNNPjuSKYSRhYWE888wz3HvvvWRlZfHKK6/QrVs3DSERETmN+Qthzx6Iqgy9r/B0Nb7HqbAdFRXFhx9+yL333sv333/PDz/8AMDChQtZuHAhxhjOP/98Xn/9dSpWrOjSgkXE+4SHh/PNN9+UbJ/JiBEjyMnJ4aGHHsLhcDBo0KBymYIvLi6OLl268MADD3Do0CE6dOjAww8/7PbPKyLiq4wxTPrM7tW+qj+Eh6tX+2w5vahN/fr1+frrr5k3bx6LFy8mJSWFoqIiatasSfv27bn88st1c6RIgAgODqZjx46lPt/hcDBmzJjj9i1evNjFVZ3c1VdfzdVXX10un0tExNetXQcbN0FYGFx9lYK2M8q0XHtQUBDdu3ene/furqpHRERERLzEpMl2r/YVl0N0tMK2M8oUtkVEwJ4O9KOPPgLg1ltvJTQ01MMViYhIWW3fYVi8xF6a/frrFLSd5dRy7bNmzeKqq65i6dKlpzxnyZIlXHXVVcyZM8fp4kTEN+Tl5fHII4/wyCOPkJeX59Q13nzzzTNO7VcWd9xxBxMmTHDb9UVE/M1nR3u1O3aAOvEK285yKmx/+eWX7N69m5YtW57ynFatWpGSksIXX3zhdHEiIiIiUv7+/NPww2x7+4brFbTLwqmwvWnTJhITEwkLCzvlOWFhYTRs2PCMi9+IiIiIiHcZP8GQnw9Nm0DTJgrbZeFU2D548OBxSx6fSrVq1Th48KAzn0JEREREPGD3HsP0Gfb2HbcraJeVU2G7cuXK7Nmz54zn7d27F4fD4cynEBEREREP+HC8obAQWreCi5orbJeVU2G7SZMmrFmzhs2bN5/ynM2bN7NmzRqaNGnidHEiIiIiUn627/hrrLZ6tV3DqbA9aNAgCgsLueuuu/j+++9POP79999z1113UVRUxKBBg8pcpIiIiIi43wfjDUVF0KkDXNBIYdsVnJpn+5JLLmHw4MGMHz+ekSNH8vTTTxMXFwfArl27OHToEMYYbr75Zrp27erSgkXE+1SoUIHPPvusZFtERHzPb78Z5v1oz6t9+xAFbVdxelGbUaNG0ahRI95++222bdtGRkZGybHzzjuPO+64g/79+7uiRhHxciEhIVx22WWeLkNERMrg3Q/sebUv7Qrnn6ew7SplWkHyyiuv5Morr2T//v3s3bsXgJo1a5ZqphIRERER8Q7rNxiWLIXgILj9NgVtV3LJcu3Vq1dXwBYJYPn5+UydOhWAa6+9Vsu1i4j4mHfft3u1e/bUapGu5pKwLSKBLS8vj+HDhwP2X7wUtkVEfMeqJMOqJAgJgdtuUdB2NafD9sGDB/n0009ZsWIFBw4cIC8v76TnWZbFnDlznC5QRERERNzDGFPSq31lX6hZU2Hb1ZwK21u3buWmm24iPT0dY4yraxIRERGRcjB/AazfABUqwM03KWi7g1Nh+8UXXyQtLY3LLruMu+66i3r16hEZGenq2kRERETETXJzDePesDtNb7gezqmqsO0OToXtlStXUr9+fV599VUsSw0jIiIi4ms+/Qz27oPq1eGmQcpz7uLUCpLGGBo1aqSgLSIiIuKD9u4zfPKp3as97G6L8HBlOndxKmw3btyY3bt3u7oWERERESkHb7xlyM2F5s3sRWzEfZwK2yNGjOCXX35h3rx5rq5HRHxQhQoV+OCDD/jggw+0XLuIiJdbvcZelj0oCO4bbmmkgps5PfXfLbfcwogRI+jTpw/t27enZs2aBAWdPLu3bt3a6QJFxPuFhITQv39/T5chIiJnUFBg+O9r9vCRfn2gQQMFbXdzKmzffPPNWJaFMYavv/6ab7755rTnb9y40aniRERERMR1pn8LW7dCxYowdIiCdnlwKmz3799ff3IQkRIFBQXMmDEDgD59+hASosVpRUS8zaFDfy1gc8cQiypVlOXKg1PviC+88IKr6xARH5abm8uQIUMASE5OVtgWEfFC731gOHQIzq0PV/bzdDWBw2vfEb/++mtWrVrF+vXr2bJlC/n5+YwdO5YBAwac1XWKioqYOHEiU6ZMYceOHTgcDtq3b8/IkSOJj493U/UiIiIi3mPLb4avjo76vW+ERUiIerXLS5nDdnp6Ohs2bCAtLY3atWvTokULV9TFq6++SkpKCtHR0VSvXp2UlBSnrjN69GimTp1KgwYNuPnmm9m/fz/fffcdixcvZvLkydSrV88l9YqIiIh4o4ICw9h/GYqKoGsXaNlCQbs8OTX1H0BqaioPPvggHTt2ZOjQoTz88MNMnTq15PjUqVNp06YNK1eudOr6zz33HPPmzWPZsmUMHDjQqWssW7aMqVOn0rp1a7788ksefvhhXnrpJV5//XXS09N59tlnnbquiIiIiK+YNBl++x0qVYL7RyholzenwnZ6ejoDBw7k22+/pUGDBgwaNAhjzHHn9OjRg6ysLH744QenCmvfvj2xsbFOfWyx4vB/3333ERYWVrK/c+fOtGnThkWLFmlxHhEREfFbO3caPhxvZ7R7h1tUraqwXd6cCttvvfUWO3fuZNiwYUybNo2nnnrqhHOqVKlCYmIiK1asKHORzvr5559xOBwnHdrSqVMnAJYvX17eZYmIiIi4XVGR4YWXDHn50KY19LzM0xUFJqfC9pw5c6hXrx4jRow47Xnx8fHs27fPqcLKKjs7mwMHDhAXF0dwcPAJx+vWrQvAjh07yrs0EREREbf76htY9wtEhMMjD2qlSE9x6gbJffv20a1btzOeZ1kWmZmZznyKMjt8+DAAFStWPOnx4v3F551KVFTUKVfGdIXo6Gi3XVu8QyC0ccWKFfnwww8BqFGjBqGhoR6uqPwEQvsGOrWx//PHNt69p5C33kkHYOT9Dho1ivBsQR7k6fZ1KmxXrFiRAwcOnPG8nTt3EhMT48yn8BoZGRluu3Z0dDRpaWluu754XiC18ZVXXgngsV+wPSGQ2jdQqY39nz+2sTGGp542ZGdDk8bQ87Ic0tKOeLosj3B3+5YmyDvVZdukSRN++eUXkpOTT3nOpk2b2LRpk8umAjxblSpVAk79xl+8v/g8EREREX8wazYs+xlCQ2HUwxZBQRo+4klOhe2bbrqJvLw8hg8fztatW084vmPHDh5++GGMMdx4441lLtIZDoeDatWqsWvXLgoLC084XjxWu3jstog4r6CggFmzZjFr1iwKCgo8XY6ISMA6eNDw6jh79pHbbrWoW1dB29OcCtuXXHIJQ4cOZfPmzfTp04eePXtiWRaLFi2iX79+9OrVi99++4277rqLVq1aubrmUmvTpg3Z2dkkJSWdcGzhwoUAtG7durzLEvE7ubm5DBw4kIEDB5Kbm+vpckREAlJRkeGf/7KXZD//PBjk3DIl4mJO3/n30EMP8Z///IeEhAS2b9+OMYYDBw6wZcsW6taty7///W/uv/9+F5Z6aqmpqWzdupXU1NTj9l933XWAvRplXl5eyf758+ezfPlyOnbsWOa5vEVERES8wedfws/LISwMRj+pJdm9RZmWa7/iiiu44oorSE1NZdeuXRhjqFmzJjVq1ChzYVOnTmXVqlUAbNmypWRf8bzYLVu25NprrwVg4sSJjBs3juHDhx83HWHbtm259tprmTp1KgMGDKBz584cOHCAmTNnUqVKFZ588sky1ykiIiLiab9vNbz5tj18ZNg9FufWV9D2FmUK28ViYmJcPuvIqlWrmDZt2nH7kpKSjhsSUhy2T+eZZ54hISGBKVOm8PHHH+NwOOjRowcjR46kTp06Lq1ZREREpLzl5hrGPGvIz4f27WBAf09XJMeyzN/XWZfjuHu6GH+bbkiOFyhtnJWVRXx8PADJyclERkZ6uKLyESjtG8jUxv7PH9r4lf8W8eVXEBMNH31gER2tXu1i3jD1n1M927fcckupzgsNDaVKlSo0atSI3r17U6tWLWc+nYiIiIicxOIlhi+/srefeExB2xs5FbaLx01blsWpOsaPPfbtt9/y3//+l4ceeojBgwc7V6mIiIiIlDh40DD2RTtrXXcNXNxGQdsbORW2586dy0cffcSnn37KFVdcQa9evUp6rffu3cvMmTOZOXMmAwcOpFevXqxcuZK3336bf/3rX5x//vl07NjRpV+EiHhWWFgYL774Ysm2iIi4lzH2NH/p6XDeeXDXHQra3sqpsL127Vo++eQT3n33XTp06HDcsYYNG9KlSxeuvPJK7rzzTpo3b86dd95J06ZNGTx4MJ988onCtoifCQ0NZejQoZ4uQ0QkYEyc9Nc0f2OesqhQQWHbWzk1z/b7779Py5YtTwjax+rQoQMtWrTggw8+AOxp+Bo2bMi6deucq1REREREWLnK8M579vCR+++1qF9PQdubORW2t23bRvXq1c94XvXq1fnjjz9KntetW5dDhw458ylFxIsVFhayaNEiFi1aRGFhoafLERHxW3v3GcY8Yygqgt69oG9vT1ckZ+LUMJLw8HDWr1+PMQbLOvlvU8YY1q9fT3h4eMm+3NxcKlas6FylIuK1jhw5Qr9+/YDAmvpPRKQ85eUZnnrakJ4BCQnwwH3WKXOYeA+nerbbt2/Pzp07efbZZ8nJyTnh+JEjR3j++efZuXPncUNNduzYoen/RERERJzw6muGjZugcmV4/h8ap+0rnOrZfuCBB1iyZAmTJk3i22+/pWPHjiUhes+ePSxatIhDhw4RExPDyJEjAdi6dSt//PEHt99+u+uqFxEREQkA335n+Ho6WBaMfsKiVi0FbV/hVNiOjY3ls88+4+mnn2bZsmV8++23J5zTrl07xowZQ2xsLADx8fEsWrSISpUqla1iERERkQCyeYvh5VfsGyJvv82i7cUK2r7EqbAN9s2O48ePZ+fOnSQlJbF//37Avinyoosuom7dusedHxYWxjnnnFO2akVEREQCSEaG4cnRhrx8aN8WbrnJ0xXJ2XIqbA8fPpxq1arx9NNPU6dOHerUqePqukREREQCWn6+4YnRhj17oXZtePIJi6Ag9Wr7GqdukJw/fz7p6ekuLkVEREREwJ7V7cV/G9ashchIeOF5i8qVFLR9kVM923FxcSedhUREAlNoaChjxowp2RYRkbL55FP47gcIDoJnnrY4t76Ctq9yKmz37t2bDz74gAMHDlCtWjVX1yQiPiYsLIx7773X02WIiPiFH38yvP2ufUPkffdaXNxGQduXOTWM5K677qJVq1bcdNNNzJ49m/z8fFfXJSIiIhJwft1oePafdtC+5moY0F9B29c51bPds2dPjDHs2bOHe++9F8uyiImJoUKFCieca1kWc+bMKXOhIuK9CgsLWbt2LQDNmjUjODjYwxWJiPievfsMox435OXZM4+M+D8FbX/gVNhOSUk57rkxhj///NMlBYmI7zly5Ajdu3cHtFy7iIgzsrIMjz5mSE2D886DMaMtgoMVtv2BU2F706ZNrq5DREREJCDl5hoee9KwdRvERMO//mnhcCho+wunxmyLiIiISNkVFhqeed6QtBocDnjxBYuaNRS0/YnCtoiIiIgHGGP49yuG+QsgNNSeS7thooK2vylT2F60aBHDhg2jU6dONG7cmMcff7zk2MKFCxk7diz79u0rc5EiIiIi/uad9wzTv4WgIBjzlEWLixS0/ZFTY7YBnnvuOSZOnIgxBofDQUFBAcaYkuPVqlXjo48+olatWgwePNgVtQaETz8z7NtnuG+ElmQVERHxV5OnGiZMtLcfesCi8yV6z/dXTvVsf/XVV3zyySdceOGFTJs2jaSkpBPOadiwIbVq1WLevHllLjKQfP6l4YtpsHmLpysRERERd/h+luG11+0OyrvusOjXR0HbnznVsz1p0iQqV67MO++8Q0xMzCnPS0xMZMsWpcazUa8u7N8PmzZDo4aerkakdEJDQ3nkkUdKtkVE5OQWLDSMfcEO2tdfCzcN8nBB4nZOhe0tW7bQpk2b0wZtgIoVK2r+7bPUMBGWr4BNmw2g33TFN4SFhTFq1ChPlyEi4tUWLTGM/oehsAh6Xg7D7rGwLL3X+zunb5AszTfH/v37CQ8Pd/ZTBKTiu5A3bfZwISIiIuIyS5YanhxtKCiA7t3gsUd0b1agcCps16tXjw0bNpCfn3/KczIzM9m0aRPnn3++08UFooaJ9uP2P+DIEXP6k0W8RFFRERs3bmTjxo0UFRV5uhwREa+y9GfDE0eD9qVd4cnHtDpkIHEqbPfs2ZMDBw7w8ssvn/KcV155hcOHD9O7d2+niwtE1apB1RgoLILffvd0NSKlk5OTQ4cOHejQoQM5OTmeLkdExGv8vNzwxJOG/HzocgmMfsIiJERBO5A4NWb71ltv5dtvv+Wjjz5i9erVdOvWDYDk5GTGjx/P7NmzWbVqFRdccAHXXnutSwv2d5ZlkZhoWLIUNm6CJo09XZGIiIg4Y8VKexn2vHzo1BHGjFbQDkROhe3w8HDGjx/PqFGjWLBgAevWrQNg5cqVrFy5EoAOHTrw0ksvERYW5rpqA0SjhhZLlho26yZJERERn7RipeHRxw15edChPTzztIJ2oHJ6UZuYmBjeeecdNm3axKJFi0hJSaGoqIiaNWvSoUMHmjZt6so6A0rxuG3dJCkiIuJ75i80jHnGHjrSvi08O8YiNFRBO1A5HbaLNWzYkIYNNSG0KxWH7Z3JkJVliIzUD6iIiIgv+GGW4Z8v2NP7dbkEnn5KQTvQOXWD5CeffEJqaqqra5GjoqMtatQAY7SSpIiIiK+Y9rXh2X/aQbtXT3uMtoK2OBW2n3vuOS655BLuuusuZsyYwZEjR1xdV8DTUBIRERHfMXGS4eX/2FP2XjMARj2iMdpic2oYyeDBg5k5cybz589nwYIFRERE0L17d/r27UuHDh0ICnJ6rRw5qmGixfwFRitJik8IDQ1l+PDhJdsiIoHCGMO77xs+/sR+fstNcMftWhlS/mIZY5xaOcUYw7Jly/jmm2+YPXs2mZmZWJZFTEwMvXr1om/fvn5xk2RaWprbrh0dHX3K669YaRj5kKF2bZjyqX558VWna2PxfWpf/6c29n9laeOCAsN/XjV8Pd1+fs9dFjfeoJDtTdz9MxwdHX3Gc5wO28fKy8tj3rx5TJ8+nQULFpCfn49lWdSpU4e+ffuW9Hj5Ik+F7UOHDb362k3z7dcWUVH64fVFeqP2b2pf/6c29n/OtnFOjmHMs4bFS8Cy4MH7Lfpfqfdqb+MNYdslXaZhYWH07NmT119/ncWLF/Pss8/SqlUrduzYweuvv+6KTxFwKleyiIu1t3WTpHi7oqIidu7cyc6dO7Vcu4j4vbQ0w70P2EE7LAyef0ZBW07N5eMTDh06xMGDB9UT4AKJuklSfEROTg7NmzenefPmWq5dRPzarl2Gu4cZNm6EypXh1VcsLumkoC2nVuZ5tgFSU1P57rvvmD59OmvXri3Z37p1a/r27euKTxGQGjW0mDvPsHGTbpIUERHxtF83Gh55zJCeDrVqwcv/sqhTR+/PcnpOh+2cnBxmz57N9OnTWbp0KYWFhRhjaNCgAX379qVv377UqlXLlbUGnOLp/zarZ1tERMSjFiw0PPO84cgRSEyAF8daVK2qoC1n5lTYfuCBB/jxxx85cuQIxhhq1qxJ79696devH4nFYx+kzBIa2Ddd7D8ABw8a/VCLiIiUM2MMn3wKb79rT1pwcRt7+XWHQ+/JUjpOhe2ZM2dSqVIlrr76avr160fr1q01n6QbOBwWdesatm+3x213aO/pikRERAJHbq7hXy8ZZs2xn199FYwYpsVq5Ow4FbZfe+01OnfuTFhYmKvrkb9pmMjRsG3o0F4/3CIiIuXh4EHDY08aft0IwUEw8j7NOCLOcWo2kh49eihol5OGifYPtmYkERERKR9bfjPccbcdtCtVglf+raAtznPJbCTiPg2Pmf7PGKPhOuKVQkJCuP3220u2RUR81dx5hrEv2jdC1q0D//qnRVyc3nvFeaV6V+zWrRuWZfHhhx8SHx9Pt27dSv0JLMtizpw5ThcY6BqcD8HBkJZm3yhZo7qnKxI5UYUKFXjppZc8XYaIiNMKCgxvvm2YPNV+3qY1/GO0RaVKCtpSNqUK2ykpKQAUFBQc91zcr0IFi3PrG377HTZtUtgWERFxtYMHDaP/YVi7zn5+840wdIhFcLCCtpRdqcL2pk2bTvtc3KthInbY3mzofIl+8MX7GGM4ePAgAFWrVtVwJxHxGUmr87n/QcPBg+BwwBOPWXTWipDiQi5frl1cL/HoTZIb9TuOeKns7GwSEhJISEggOzvb0+WIiJyRMYbPvzTcNvQQBw9CvXrw3lsK2uJ6upPJBzRqaD/qJkkREZGyy8w0vPSyYe6P9vNLu8Koh7VQjbiHU2E7KSmJn3/+ma1bt3Lo0CEsyyIqKorzzz+fiy++mGbNmrm6zoB2bn0IC4XMTNi+A+rX83RFIiIivmnTJsPoZwy7d9sTEDxwv4N+fXLUkSVuc1Zhe9OmTTz++ONs3LgRsHtZj1X8jdq0aVOef/55zj//fBeVGdhCQy1atjQsXQbzfjTcfpteEERERM6GMYYpn8ObbxsKCqBmDRgz2qJTxwjS0o54ujzxY6UO2+vWrePWW28lJyeHiIgILrnkEho1akR0dDTGGNLS0ti4cSMLFy5k7dq1XH/99UyYMIELLrjAnfUHjB7dLZYuM/wwG4YM1lASERGR0srIMDz/gmHJUvt550vg0YctKmtaPykHpQrbhYWFPPzww+Tk5HDNNdcwatQoKlaseNJzMzMzGTt2LF988QUPPvggM2fOVDB0gU4dICIcdu+GDb9C4ws9XZGIiIj3W73G8Ozzhv0H7CGZw4dZXHUlyiZSbko1G8ncuXPZsWMHvXr14rnnnjtl0AaoWLEizz//PD179mT79u3MmzfPZcUGsogIi0susbd/mG1Of7KIiEiAy8szvPFWEfeOtIN2fDy8/YbFgP6WgraUq1KF7R9//JGgoCBGjhxZ6gs/+OCDAFo90oUu72G/OMybZ690JeItQkJCuOGGG7jhhhu0XLuIeNy2bYY77zF8+hkYA317w/tvWzRooJAt5a9U74rr16+nfv36xMfHl/rC8fHxnHvuuWzYsMHp4tatW8drr73G6tWrKSgoICEhgcGDB9OrV69SffyXX37JY489dsrjH3/8MRdffLHT9ZW3FhdB1Rg4mAo/L4cO7T1dkYitQoUKvP76654uQ0QCXFGRYeoX8PY7hrx8qBJlj83u1FEhWzynVGH7wIEDtGzZ8qwvXq9ePVatWnXWHwewbNkyhg4dSlhYGL179yYyMpJZs2YxcuRI9u7dy5AhQ0p9rW7dutGoUaMT9sfGxjpVm6eEhFh072aYPBVmzTZ0aK8XDxEREYB9+w3/fMGwKsl+3r4tjHrEIiZG75XiWaUK25mZmVSqVOmsL16xYkUyMzPP+uMKCgp46qmnsCyLiRMnlgTlYcOGcc011/DKK69w+eWXlzosd+/enQEDBpx1Hd7osh4Wk6caFi6GrCxDZKReRMTzjDElK0c6HA6NhxSRcmOMYca38NobhuxsCA+H4f9ncWVf3QQp3qFUY7YLCgqc+oYNCgqisLDwrD9u2bJl7Ny5kz59+hzXI12pUiXuvvtu8vPzmTZt2llf1x8kNIB6dSEvD+Yv8HQ1Irbs7Gzi4+OJj4/Xcu0iUm727jU88LDhX/+2g/aFF8AH71r076ebIMV7eOWdTMuXLwegY8eOJxwr3rdixYpSX+/XX38lPT2dgoIC4uLiaNeuHdHR0a4ptpxZlsVlPeCd9ww/zDb0ukIvJiIiEliMMXw9HV5/05CTA2FhcOdQi2uvhuBgvS+Kdyl12P7qq6/46quv3FjKX7Zv3w5A3bp1TzhWrVo1HA4HO3bsKPX1JkyYcNzz8PBwhg0bxp133lmmOj2lRzd45z1IWg0HDhiqVdMLi4iIBIbdewz/eumvsdlNGsNjj1rUidd7oXinUoftvy/NXlrO/BmneJz3qcaJV6xYkcOHD5/xOnFxcTz11FN07NiRmjVrkpGRwdKlS3nllVd4+eWXiYiI4Oabbz7tNaKioggKKtVoG6c408MeHQ0tW2SwKqmAxUsjuO3WCDdUJq7iq39FORthYWEl29HR0URGRnqwmvIVCO0b6NTG3iE/3/DxJ0d4461sjhyxx2bff6+DQQPDy9ybrTb2b55u31KF7U2bNrm7Drdo06YNbdq0KXkeHh5O//79ufDCC7n66qsZN27cGecFzsjIcFt90dHRpKWlOfWxl3a1f6v/6uts+vc74uLKxFXK0sa+JCsrq2Q7LS2NvLw8D1ZTfgKlfQOZ2tg7/LrRHpe9dav9vMVF8MiDFnFxRzh0qGzvgWpj/+bu9i1NkHdfl20ZFK9Qearea2dnRynWoEEDWrZsSXp6OluLf3J9TNcuEBoKv2+1J+8XERHxN1lZhv/+r4i7/s8O2pUrw+OPWrz6ikVcnIaNiG/wyrBdr149gJOOyz5w4ADZ2dknHc99Nop/E8nJySnTdTylciWLtkfX45k1R2FbRET8hzGGn+YbbrrV8PmX9iqQl18GEz+26HWFZhoR3+KVYbt169YALFq06IRjxfuKz3FGYWEh69evB6B27dpOX8fTipdvnzXHXjVLxFOCg4Pp168f/fr1Izg42NPliIgP27nTns7vyacNB/6E2Nrwn39bPPV4ENFVFLLF93jl1H/t2rUjPj6eGTNmcMstt5TMtX348GHeeustQkND6d+/f8n5+/fv5/Dhw1SvXv244SXr16+ncePGx127sLCQf//73+zYsYOLL76Y6tWrl8vX5A7t2kLFSNi/HxYthks6eboiCVTh4eGMHz/e02WIiA/LzjZ8NMFeJbmgAMJCYdANcPONFhUqKGSL7/LKsB0SEsJzzz3H0KFDufHGG49brj0lJYVHH32UuLi4kvNfeeUVpk2bxtixY49bKfLqq68mMTGRxMREatSoQUZGBsuXL2f79u3UrFmT559/3hNfnstUqGBx1VWGCZ/Aex8YOrTX/KIiIuJbjDH8OB/GvW7Yf8De164t3D/CIjZW72ni+7wybAO0bduWTz/9lP/973/MnDmTgoICEhISeOihh+jVq1eprjFkyBDWrFnDkiVLyMjIIDQ0lDp16nDPPfdw2223ERUV5eavwv1uuN5i2jTDtj9g7o9wWXdPVyQiIlI6v/1ueO11Q9Jq+3mtmnDfCIsO7bXUuvgPyzg7gXaAcPd0Ma64/kcTDO++b4iLhU8+sggJ0QuUtwiUKaWysrKIj48HIDk5OWDm2Q6U9g1kamP3SE01vPuBYca39s2PYaFw4yC4aVD5DxlRG/s3Tf0nLnHt1VClCuxKge++93Q1IiIiJ5eXZ5g4yTDwJsP0GXbQ7tbVnmXk9tuCNDZb/JLXDiOR0nM4LG65Ef73uuHDjw2XXwZhYXrBEhER72BP5QdvvmPYvdvel5gA9w63aNZU71fi3xS2/cSV/eCzKfbMJF9Pt3u7RUREPG3NWsMbbxl+3Wg/r1oV7r7D4vLLIChIQVv8n4aR+IkKFSxuvcV+0fr4E0NOjobii4iI52z7w/DIY0UMv88O2hHhcNutMGmCxRU9LQVtCRjq2fYjva+ATydBym74/Eu4+UZPVyQiIoFm/37D++MN330PRUUQHAR9+8Jtt1hUraqALYFHPdt+JCTEYshg+4Xs088Mhw+rd1tERMpHWprhf+OKGHij4duZdtDucglMGG/x0MggBW0JWOrZ9jPdu8GET2H7dpg81TB0iF7cxP2Cg4Pp0aNHybaIBI5Dhw2TJhs+/xxyjtj7mjeDu++0aHyh3oNEFLb9THCwxR1D4InR9pK3/a80nKPeBHGz8PBwJk+e7OkyRKQcZWcbpn4Bkz4zZGbZ+xo1hDuHWrRqqUVpRIopbPuhSzpBo0awcSO8+JLhX2P1oiciIq6RlWX4YhpMnmLIOGTvO+9cGDrEomMHvd+I/J3Cth+yLIvHHobb7zIsWQYzvoW+fTxdlYiI+LLMTMPnX9pDFA8ftvfFx8OQwRbdumoaP5FTUdj2U+eea3HnUHj9TcP/Xje0aAGxtfVCKO6RlZVFYmIiAJs3bw6Y5dpFAsHhw3+F7MxMe1+deBh8i0W3S+3hiyJyagrbfuz6a2HxElizFp4fa3jtv3pRFPfJzs72dAki4kKpqYbJUw3TvobiH+96deHWWywu7aL3E5HSUtj2Y0FBFk+MgltvN6z7xV5h8sYbPF2ViIh4s917DJ9+Zpg5E/Ly7X3169k92V06K2SLnC2FbT9Xq5bFfSNg7L8M731guLgNnH+eXihFROR427YZJk4yzJkLhUX2vgsvgJtvtGjfTmOyRZylsB0AevWEhYtg0WJ49p+Gd9+EsDC9aIqIBDpjDKuSYNJkw8/L/9rfupUdsi9qrtlFRMpKYTsAWJbFow/B+g2GrVvhg/GGu+/Ui6eISKAqKDDM+wk+m2zY8pu9LyjInjr2phssGjbUe4SIqyhsB4joaItHHoTHnzJMnATNmhnaXawXUxGRQHL4sGHGTJj6hWH/fntfeDj0vgKuu8YiNlbvCyKuprAdQC7pZNG3j2H6DBg9xvDG/6BBA72wStkFBQXRoUOHkm0R8S47dhimfmn4/gc4cnRJ9ZhouHqARf9+EBWl9wIRd1HYDjAP3Gexe7c9Ru/hxwxvvwE1qutFVsomIiKC6dOne7oMETlGUZHh5xUw9XPD8hV/7T+3vt2L3aM7VKig138Rd1PYDjChoRbP/QP+b4Thj+3wyCjDG69BZKRecEVE/MHhw4bvvodp3xiSk+19lgUd2tshWzc9ipQvhe0AVKmSxUsvwF3/Z9i6DZ582vDSCxASohdfERFf9dtvhi+/MsyaA7m59j6HA/r0gquv0nhsEU9R2A5QNWtavPgCDL/XsGIl/PsVw6MPq7dDnJOVlUXz5s0BWLNmjZZrFyknubmGn+bDtK8N6zf8tf/c+jCgv8VlPcDh0Ou6iCcpbAewxASLMaPhsSftu9Nr1YJbb/Z0VeKrDh486OkSRALGH9sN30w3fD8LDh+29wUHQ5dLYMBVFk2bqPNExFsobAe4Du0t7r8XXvmv4d33DRER9pg+ERHxLrm5hh/nw9ffGH5Z/9f+GjWgb2+Lvr2halW9fot4G4VtYUB/iwMHDBMmwv/GGbKz7R5u9YqIiHiWMYZNm+HbmfYy6plZ9v7gIGjfHvr1tWjTCoKD9Xot4q0UtgWAO4daVKgA731geO8DQ1YW/N/dCtwiIp6Qlm6YNdsO2dv++Gt/zRrQt49F7yvgnHP0+iziCxS2BbBD9eBb7DvX/zfOMGkyZGUbHrxfPSYiIuUhP9+wdBl8P8uweAkUFtr7w0Khc2fofYVFi4sgKEivySK+RGFbjnPdNRYOB7z4b8M30yEry/DU45oWUETEHYwxbNwE3/9gmDMPDh3661jDROjdy6LbpVC5kl6DRXyVwracoE8vO3A/85xh7jzIyTH8YzREROjFXk4uKCiIiy66qGRbRE4vZbc9BvuHWYadyX/trxoDl/WAnpdbnHeuXnNF/IHCtpzUpV0sIiLgiacMS5bC3cMN/3wGLYogJxUREcHcuXM9XYaIV0tNNcz7CWbPMWz49a/9FSpA505w+WUWLVvoL4ki/kZhW06p3cUW/30Znhxt2LoVbr/L8PRT9n4RETmzw4cNCxfBnHmGlaugqMjeHxQELVtAj24WXTpr4RkRf6awLafVtInF++/YS7pv+BUeGWW4/Ta45SbdpCMicjJZWYZFi2Huj4blK6Cg4K9jFzSCHt0tLu2iObFFAoXCtpxRtWoWr/3XnqXkq2/s6QE3b4YnHoOKFfVmIZCdnU27du0AWLp0KQ6Hw8MViZSvrCzD4qXw40+Gn3+GvPy/jp1bHy7tatH9UoiL02umSKBR2JZSCQuzeOgBi4YNDa/8x7BwMdxxj+HpJ6Fhot48Ap0xhuTk5JJtkUCQkWFYtATmzzesWAX5xwTsunXg0q7QtYvFufX1GikSyBS25az06WVx3rnwxGhDcjLcdY/hlpsNt9xkERqqNxQR8W8HDhwN2AsMq1dDYdFfx+rWgc6XQLdLLc6tr0XBRMSmsC1nrVFDiw/egVdeNcz7ET78CBYuMjzxGDQ4X28uIuI/jDFs+8Meg71wkT0n9rHOPw+6dLbofAnUr6fXPxE5kcK2OKVKFYtnnrbofIk9rOT3rTD0LsNtt8JNgzR1lYj4roICw7pfYPFSw5Kl6SQn/zU0yrLgwgugU0eLzp00BltEzkxhW8qkW1eLi5rBS6/Y01u994H9+NADdg+4iIgvyMgw/LwcFi+xHzOzio8UERoKrVraAbtDO80iIiJnR2FbyiwmxuKfz8LsOfCf/xk2b4E77zFc0dNw11BLb0wi4nWMMfz2Oyz7GZb9bFi/4a85sAGqREHbtnBZj4o0viBL82CLiNMUtsUlLMvish7QogW8+bbhh1kw8zv4ab5h8C1w7dXoBko/ZlkWiYmJJdsi3ujwYXvWkGU/29PzHUw9/vh550H7dtChnUWjhhAcbBEdXYG0tGzPFCwifsEymqfrtNLS0tx27ejoaLde35PWbzD893+GTZvt5/HxcO8wi7YXB1YY8+c2FrWvtysosF+Dfl5uWLESft14fO91RDi0bAlt29ivTTVrnvjapDb2f2pj/+bu9o2Ojj7jOerZFrdofKHFO2/Cdz/A2+/Y0wQ+PMrQrCncfhu0uChwAreIlA9jDLtSYNUqWLHKsGrVsWOvbfXqwsUXQ7uLLZo2sdcQEBFxJ4VtcZugIIveV0CXS+CjCYbPv4C16+DekYaLmhtuv82ieTO90YmI8/48aEhKgpVJhlVJsG/f8ccrVYLWraBNK4tWraBmDb3miEj5UtgWt4uMtPi/uy2uvdowYaJh+reweg0Mv8/QsoVhyGCLZk31BujLsrOz6datGwBz587Vcu3iNmnphjVrYPUaQ9Ia2L79+OMhIdD4QmjV0qJNa0hMsMdei4h4isK2lJtq1SweuN/ixkGGCZ8YZsyEVUmwKsnQpLHh+mstOnXUG6MvMsawefPmkm0RV0lLM6xZC2vWGpJWwx/bjz9uWdCgAbRqAS1b2ENDIiL0GiIi3kNhW8pdjeoWDz1gcdMgw0efGL77Hn5ZD7+sN9Sqac9c0ruX3SMuIoFl717DmnWwdp1h7VrYmXziOeedCxc1h4uaWzRvBlFReq0QEe+lsC0eU7OmxaMPWdx+m2HaV4avvoY9e+F/rxveHw99ehmu7GdRJ15vpCL+qKDAsO0P+OUXWLfe8Mt62L//xPPOOxeaNbVvrG7ezF7BVkTEVyhsi8edU9Xijtstbr7R8MNsmDLVsGMnTJ4Kk6camjcz9Olt0eUSCA/Xm6yIrzp0yPDrRtjwq70c+q8bISfn+HOCg+1x1s2aQrNmFk0bQ+XK+rkXEd+lsC1eIzzc4sq+0Lc3/LwCpk0zLFtOyXjN/7wKPbob+vSySEwIrPm6RXxNQYHhjz9gw692uN7w68mHhERGwoUXQNMmFk0awwWNNOZaRPyLwrZ4naAgi3ZH58E9cMAw83uYMdOwZw989TV89bWhbh3o3g26dYU6dfTGLOJJxhh274ZfN8HGjYaNm2DLb5Cbe+K5cXF2uG5yoUWTJva817opWkT8mcK2eLVq1SxuvRluvhGSVtuhe8EC2LET3v/Q8P6HkNDA0O1Si0u7QK1aetP2BMuyiI+PL9kW/2WMYc9e2LwZNm+xV2jc8hscOnTiuZGR0KihHa4bX2gvga7x1iISaLRc+xlouXbvk5VlWLgY5sy1l2AuLPzrWEICdGxv0aE9JDTwfPBTG/s3f2/fwkJD8i747Xf47TfDlt9OHaxDQ+H88+GChtCokUWjRIiPt/9S5cv8vY1FbezvtFy7iBMiIy16XgY9L7NITzfMXwjzfrTn4N2yBbZsMXwwHqpXg/btDR3bW1zUHCpU8O03fRF3ysqyZwbZuhV+32r47XfYug2OHDnx3JAQe4aQxERITLBIbADnnQehofoZExH5O4Vt8WlVqtg3VV7Z1yItzbB0GSxaYli+AvYf+GuMd1goNG1qaN3KonUrOP883+9xE3FGQYHdW73tD/jjD2OH622wZ8/Jzw8Pt39eGjSABufbNyefW1/BWkSktDSM5Aw0jMQ35ebaPd2LFhuWLIUDfx5/vEoUtGxhL4rRtKl9k5Y7wnegtHFOTg59+vQBYMaMGURERHi4ovLhze1bUGBISYHtO+0lzf/Ybs8OsmMnFBSc/GOqnWP3UJ93rh2sExpAbGxg38DozW0srqE29m8aRiLiJhUqWLRrC+3aWhhjz9u9YiWsWGlYvQbSM2DujzD3R/t3zajK0LSJoWlTi2ZNocH56rk7G0VFRaxevbpkW8pPVpZhZzLs3Ak7kw3bd8COHbAr5dShOiIC6teD+vXhvHMtzj8asLUSo4iI6ylsi9+zLIt6de3e62uvtigosOf8XbnKXlhj/QbIOAQLF8PCxXb4DguF8883XNAIGjW0aNQI4mI19EQ8IzfXkLIbUlIgeRfsSjHs3AnJyXAw9dQfFxEOderY3/v161ucW98eAlK9ur6XRUTKi8K2BJyQELv3ullTO2zk5xs2b4F1v8DadfaS0YcO2avb/boRwA7gFSva0wyef579J/bzz7dDjHrAxRUOHbbnqk7ZDbt3w+499vNdKfYS5qcb8Fc1xp75o04dqFfHom5dqFvXvklYoVpExLMUtiXghYZaNL4QGl8IgwZaJ12gY/MWyMy05/pOWg3FATwkBOrVNdSvD/XqWtSrB/Xq2ONcQ0IUcsRmjOFwJuzbB3v3wp69sHevOfpoP8/MPP01IiPtv67ExdmPdeIt6tSB+DioWFHfayIi3sqrw/a6det47bXXWL16NQUFBSQkJDB48GB69epV6mvk5eXxzjvv8M0337Bnzx6ioqLo2rUr999/P1WrVnVj9eKrLMsiNtYOzD262SGmoMCeFu333+G33+1p0X7fagek37fa/4oDONghPC7OcF79Q9SoUURsrFUSlKqdo95Gf1IcpP88YN+Ie+BP+PNP2LfPsG+/3Su9bx/knGQKvb+rGgO1ax/9Vwtq17aIj7O/F6tEeX7eeBEROXteG7aXLVvG0KFDCQsLo3fv3kRGRjJr1ixGjhzJ3r17GTJkyBmvUVRUxD333MOiRYto3rw5l112GTt27GDq1KksXbqUKVOmEBMTUw5fjfi6kBB7ZoaEBgB24DHGsG+fHbT/2A47dhq2b7dvTss5Ys8AsX17/tEr/BXEw0KhZk1DzZrY/2pY1Kxhb1evBlWramiKNygoMBw6BKlpkJYGfx6EgwchNdVw8KA9VjotLY19+81J56I+mSpRlLR7rZpQs6ZlP9aAWrUgIkLtLiLib7xy6r+CggKuuOIK9u7dy5QpU2jUqBEAhw8f5pprriElJYUffviB2NjY017niy++4PHHH6dPnz78+9//LukVmjRpEmPGjOH666/nmWeeOe01NPWfnK2iIsP+A3bYTk2LYMtv2aSk2GNvd+8+fsXLU4mOtnvAz6kK51Szezyjoy1iou1j0dEQXcUeR+4NvZ1ZWVk0b94cgDVr1hAZGenZgv6mqMiQlQWHM+3x+IcP2zPSZGRARoYp2U5Pt4N1Wrr9/GxeHStVstusWjX7sXp1ixo17F+gih/Dwz3fVnJ29Drt/9TG/s0bpv7zyrC9aNEibr/9dgYMGMDYsWOPOzZt2jRGjRrFiBEjGD58+GmvM3DgQFavXs28efOOC+bGGHr06MHBgwdZunQp4eHhp7yGwraUxd/buKDAHlqwdy/s3WeP2y0es7tvnz0E4VTTtZ1McDBUrmxPXVip0l/bFSvaY3wrVrSIdBRv21O+RUTYC5U4jj5WqOCdw1oKCw25uZCba69ieOToY3Y25OTYj9k5xduGzEzIyoLMrKOPmfa/w0cfnZmR0LIgKsr+xaZq1WP+xVhUjYF69SoRXuEw55yjIO2v9Drt/9TG/s0bwrZXDiNZvnw5AB07djzhWPG+FStWnPYaubm5rF27lvr165/QA25ZFu3bt2fy5MmsX7+eVq1auahykdMLCbGIrQ2xtYv3HB/QiooMGYf+Pv7XkJb213CG4p7XrCy7l7x438mV7nfpChUMYWH2EJewsL/+hYRCSLAd6oOD7bHowcEQHARWkB1Gi/8FWfZnMwZM0fHbhUV2rYWF9i8Tx27n5UN+HuTnH93Oh7yjz12tQgWoXMn+xaRKFTtIR0XZwzuioiyiovjrrwdHj59uQZfo6FDS0hSyRUTk1LwybG/fvh2AunXrnnCsWrVqOBwOduzYcdpr7Ny5k6KiIurVq3fS48X7t2/fftqwHRUVRVBQUKnqdkZpfiMS33a2bVy1qj0X8pkcOWLIyCgiPcOQnm5IzygiI8OQkWE4dKiIrCzD4UxDZqbh8GH7MTvbkJNz9N8x44yLe5C9VXi43XMcHm4R6bBwOCwiIy0cDuznkRaVKlpUrBhEpUoWFSsWP7eOhuggoipbhIW5Z5VQ8W9qY/+nNvZvnm5frwzbmUfnwKpUqdJJj1esWJHDhw+f9hrFxytWrHjKaxz7uU4lIyPjtMfLQn+68n/ubuOwMHsscPVqZ/NRFmBRVGRKhmjk5tq9ybl5f/Us5+baPc8Ff++RLrB7qo2x/xUVQV5uDv/97/VYwL33T6ZChYiS3m4suye8pFf82B7yELs3PTT0r1710DD7eUS4HbLDwv4+Lt1Q2h77Y2UdHV7iSvoZ9n9qY/+nNvZvGkYiIh4TFGSVjOEuq6wsw7D/WwzAlX0NkZEaWiEiIgLgvvERZVDc63yq3uvMzMxT9noXKz5+qp7r4v2n6vkWERERESkrrwzbxeOpTzYu+8CBA2RnZ590PPex4uPjCQoKKhn//XfF+081pltEREREpKy8Mmy3bt0asKcA/LvifcXnnEp4eDhNmzbljz/+ICUl5bhjxhiWLFmCw+GgcePGLqpaREREROR4Xhm227VrR3x8PDNmzGDjxo0l+w8fPsxbb71FaGgo/fv3L9m/f/9+tm7desKwk+uuuw6AV155hWOnE//ss89ITk6mb9++p51jW0RERESkLLzyBsmQkBCee+45hg4dyo033njccu0pKSk8+uijxMXFlZz/yiuvMG3aNMaOHcuAAQNK9l911VXMnDmTGTNmsGvXLlq3bs3OnTuZNWsWcXFx3H///R746kREREQkUHhlzzZA27Zt+fTTT2nRogUzZ85k0qRJVK1alf/85z8MGTKkVNcICgrizTffZMSIEaSmpjJ+/HiSkpK45pprmDx5MjExMW7+KkQCh8PhwOFweLoMERERr+KVy7V7Ey3XLmWhNvZval//pzb2f2pj/+YN82x7bc+2iIiIiIivU9gWEREREXEThW0RKbMjR45w/fXXc/3113PkyBFPlyMiIuI1vHI2EhHxLYWFhcyePbtkW0RERGzq2RYRERERcROFbRERERERN1HYFhERERFxE4VtERERERE3UdgWEREREXEThW0RERERETfRcu0iIiIiIm6inm0RERERETdR2BYRERERcROFbRERERERN1HYFhERERFxE4VtERERERE3CfF0Af5k3bp1vPbaa6xevZqCggISEhIYPHgwvXr1KvU18vLyeOedd/jmm2/Ys2cPUVFRdO3alfvvv5+qVau6sXo5k7K0rzGGBQsWMG/ePJKSkti9ezcFBQXUrVuXXr16cdttt1GhQoVy+CrkdFzxM3ysjIwM+vTpw/79++nYsSPvv/++iyuWs+WqNj548CBvv/02P/30E3v27MHhcFCvXj2uvPJKBg0a5Kbq5Uxc0b779u3j3XffZcmSJezevRuHw0HdunW5/vrr6du3L8HBwW78CuR0vv76a1atWsX69evZsmUL+fn5jB07lgEDBpzVdYqKipg4cSJTpkxhx44dOBwO2rdvz8iRI4mPj3d53Zr6z0WWLVvG0KFDCQsLo3fv3kRGRjJr1ixSUlJ49NFHGTJkyBmvUVRUxB133MGiRYto3rw5rVu3ZseOHcyePZu4uDimTJlCTExMOXw18ndlbd/c3FyaNm1KWFgYbdq0ISEhgby8PBYtWsT27dtp0qQJEyZMICIiopy+Ivk7V/wM/92DDz7IvHnzyM7OVtj2Aq5q440bNzJkyBAOHTpE586dOe+888jOzmbr1q2Ehoby7rvvuvkrkZNxRfsmJydz7bXXkp6eTseOHUlMTCQzM5O5c+dy4MABBgwYwNixY8vhq5GTufTSS0lJSSE6OhqHw0FKSopTYfvJJ59k6tSpNGjQgM6dO7N//36+++47IiMjmTx5MvXq1XNt4UbKLD8/33Tv3t00btzY/PrrryX7Dx06ZC677DJz4YUXml27dp3xOp9//rlJSEgwDzzwgCkqKirZ/+mnn5qEhATz1FNPuaV+OT1XtG9eXp554403THp6+gn777rrLpOQkGDeffddt9QvZ+aqn+Fjff/99yYhIcF88sknJiEhwQwZMsTVZctZcFUbHz582HTp0sW0bdvWbNy48aSfR8qfq9r36aefNgkJCWb8+PHH7c/IyDBdunQxCQkJZ/1aIK6zePHikv//t99+2yQkJJgvvvjirK6xdOlSk5CQYG688UaTm5tbsv+nn35y22u1xmy7wLJly9i5cyd9+vShUaNGJfsrVarE3XffTX5+PtOmTTvjdaZOnQrAAw88gGVZJfsHDhxIfHw806dP58iRI67/AuS0XNG+oaGh3HPPPURFRZ2w/6677gJgxYoVri9eSsVVP8PFUlNTGTNmDFdeeSWdO3d2R8lyllzVxp9++im7d+/mwQcfpGHDhiccDwnR6ExPcFX7JicnA5zwc1u5cmVatGgBQFpamgsrl7PRvn17YmNjy3SN4qx13333ERYWVrK/c+fOtGnThkWLFrF79+4yfY6/U9h2geXLlwPQsWPHE44V7ztTkMrNzWXt2rXUr1//hG8ky7Jo37492dnZrF+/3kVVS2m5on1Pp/jNWeMAPcfVbfz0008THBzME0884ZoCpcxc1cYzZ87Esiwuv/xytm3bxoQJE3j33XeZO3cueXl5ri1aSs1V7ZuQkADA/Pnzj9t/6NAhVq9eTbVq1Tj//PPLWq540M8//4zD4Sj55elYnTp1Av76fnIV/QruAtu3bwegbt26JxyrVq0aDoeDHTt2nPYaO3fupKio6JTjhIr3b9++nVatWpWlXDlLrmjf0/niiy8A6NChg9PXkLJxZRt//fXXzJo1i9dff52oqCgOHz7sylLFSa5o47y8PLZs2UJMTAwTJkzgtddeo6ioqOR4fHw8r7/+OomJiS6tXc7MVT/Dt99+O/PmzWPs2LEsXLjwuDHb4eHhjBs3jvDwcFeXL+UkOzubAwcOkJCQcNIOruLvn7K8p5+MerZdIDMzE7D/XHUyFStWPOMbbvHxihUrnvIax34uKT+uaN9TmT9/PpMnT+a8887j2muvdbpGKRtXtfG+fft4/vnn6dOnD927d3dpjVI2rmjjjIwMCgsLSU9P54033uDhhx9myZIlLFiwgP/7v/9j165d3HPPPeTm5rq8fjk9V/0Mn3POOUyePJlOnTqxcOFC3nvvPT777DMOHz5M//79Tzp0SHxHabOWqztJFLZFPGTdunWMHDmSSpUq8eqrrx43dkx805NPPklISIiGj/ip4l7swsJCbrjhBoYMGULVqlWpUaMG9913Hz179iQlJYXvv//ew5WKs3bs2MENN9xAamoqEydOJCkpifnz5zNs2DDeeOMNBg8eTGFhoafLFB+jsO0CZ/pNKDMz85S/bRcrPn6qnuvi/af6bUzcxxXt+3e//PILt99+O0FBQbz33ns0aNCgzHWK81zRxtOmTWPBggWMHj1aU3R6IVe+ToM9BdnfFe/TvTXlz1Wv06NGjWL37t289dZbtGrVisjISGrWrMmdd97JTTfdxOrVq/n2229dWruUn9JmrbN9Tz8ThW0XKB5PfbIxPgcOHCA7O/uk48iOFR8fT1BQUMm4s78r3u/yuR/ljFzRvsf65ZdfGDJkCEVFRbz//vs0bdrUVaWKk1zRxr/++itg3+GemJhY8q9bt24ALFq0iMTERK688krXFi+l4oo2djgc1KhRA7Bnp/i74n0aRlL+XNG+mZmZJCUlcd5551GtWrUTjl988cWAPc+6+CaHw0G1atXYtWvXSf9CUfz9czbv6aWhsO0CrVu3Buw3078r3ld8zqmEh4fTtGlT/vjjD1JSUo47ZoxhyZIlOBwOGjdu7KKqpbRc0b7FioN2YWEh7733Hs2aNXNdoeI0V7TxRRddxDXXXHPCv+KV62rWrMk111xDjx49XFy9lIarfo7btm0LwO+//37CseJ9ZZ2aTM6eK9o3Pz8fOPXUfqmpqQAa8ufj2rRpQ3Z2NklJSSccW7hwIVD69/RSc/nM3QEoPz/fdOvW7bST6ScnJ5fs37dvn/n999/NoUOHjruOFrXxTq5q319++cW0atXKNG/e3KxcubLc6pczc1Ubn0xycrIWtfECrmrjVatWmYSEBNO7d2+TkZFRsn///v2mU6dOpmHDhmbbtm3u/4LkOK5q38svv9wkJCSYKVOmHLc/IyPD9OzZ0yQkJJjFixe794uRUjnTojYHDx40v//+uzl48OBx+z2xqI2Wa3eRs1kmdtSoUUybNu2EJUZPtlz7zp07mTVrFrGxsUydOlVjQT2krO2bnp7OZZddRkZGBp06dTppj3alSpUYPHhweX1J8jeu+Bk+mV27dtGtWzct1+4FXNXGL7zwAh9++CG1atWia9euFBQUMHfuXA4ePMgDDzxQslCVlC9XtO/8+fP5v//7PwoKCmjXrh2NGjXi0KFDzJs3j9TUVC6//HL+97//eeLLE+wFaVatWgXAli1b2LBhAy1atCgZ9tGyZcuSmb1ee+01xo0bx/DhwxkxYsRx1/n7cu0HDhxg5syZREZG8tlnn1G/fn2X1q15tl2kbdu2fPrpp/zvf/9j5syZFBQUkJCQwEMPPVTyZ+QzCQoK4s033+Sdd97h66+/Zvz48VSpUoVrrrmG+++/X0Hbg8ravpmZmWRkZAD2n6mK/1R1rNjYWIVtD3LFz7B4N1e18ahRo0hISGDixIlMmzYNy7Jo1KgR//jHPzRMyINc0b6dO3dm0qRJvP/++6xatYoVK1YQFhbGeeedx7Bhw7jhhhvc/FXI6axateqElUCTkpKOGxJSmml0n3nmGRISEpgyZQoff/wxDoeDHj16MHLkSOrUqePyutWzLSIiIiLiJrpBUkRERETETRS2RURERETcRGFbRERERMRNFLZFRERERNxEYVtERERExE0UtkVERERE3ERhW0RERETETRS2RURERETcRGFbRERERMRNtFy7iIiUSExMPO75q6++Ss+ePUuep6amMnfuXNatW8e6dev47bffKCwsZOzYsQwYMOCk15wzZw7Dhg07bt/mzZtdX7yIiBdS2BYRkeM4HA4uv/xyAGJjY487lpSUxJNPPnlW16tVqxZXXXUVAD/88APZ2dmuKVRExAcobIuIyHGio6N54YUXTnqsatWqDBo0iMaNG9OkSRMmTJjAlClTTnu9Cy+8sOR6y5cvV9gWkYCisC0iIqV20UUXcdFFF5U8tyzLg9WIiHg/3SApIuLn5s6dy/XXX0+zZs24+OKLGTFiBH/88QevvfYaiYmJfPnll54uUUTEb6lnW0TEj02aNIkxY8ZgWRatWrWiWrVqrF27lmuvvZauXbt6ujwREb+nsC0i4qdSUlIYO3YsoaGhvPnmm3Tq1AmA/Px8HnvsMb755hsPVygi4v80jERExE998cUX5Obm0rt375KgDRAaGsoTTzxBRESEB6sTEQkMCtsiIn5q5cqVAPTq1euEY9HR0XTo0KG8SxIRCTgK2yIifmr//v3AiXNlFzvVfhERcR2FbRERERERN1HYFhHxU9WqVQPsGyVPZvfu3eVZjohIQFLYFhHxU61atQLg+++/P+FYeno6ixcvLu+SREQCjsK2iIifGjBgAGFhYUyfPp0lS5aU7M/Pz2fs2LFaNl1EpBxonm0RET8VHx/PqFGjeOaZZ7j99ttLFrVZs2YNhw4dom/fvkyfPv2sr3vdddeVbO/atQuAN954g88++wyACy64gDFjxrjkaxAR8XUK2yIifuzGG2+kRo0avPPOO6xbt44KFSrQqlUrHnzwQWbOnOnUNdeuXXvCvuTkZJKTkwGoUKFCmWoWEfEnCtsiIn6ue/fudO/e3WXX27x5s8uuJSLi7xS2RUTkOGlpaYwaNQqwe8abNGlSputt2LCBCRMmlFxbRCSQKGyLiMhxsrOzmTZtGgBdunQpc9jes2dPyfVERAKNwraIiJRwxxCR7t27a+iJiAQsyxhjPF2EiIiIiIg/0jzbIiIiIiJuorAtIiIiIuImCtsiIiIiIm6isC0iIiIi4iYK2yIiIiIibqKwLSIiIiLiJgrbIiIiIiJuorAtIiIiIuIm/w/g+7sZ0X+2igAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 720x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# Figure 7.5\n",
+    "p = np.array([0.3, 0.7])\n",
+    "q = np.arange(0.01, 1, 0.01)\n",
+    "DKL = np.sum(p * np.log(p / np.array([q, 1 - q]).T), 1)\n",
+    "\n",
+    "plt.plot(q, DKL)\n",
+    "plt.xlabel(\"q[1]\")\n",
+    "plt.ylabel(\"Divergence of q from p\")\n",
+    "plt.axvline(0.3, ls=\"dashed\", color=\"k\")\n",
+    "plt.text(0.315, 1.22, \"q = p\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.13 & 7.14\n",
+    "\n",
+    "We don't have a `lppd()` function, but we can extract the log-likelihood and finish the computation ourselves. T"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ 0.46031432,  0.49022383,  0.39503769,  0.47387998,  0.33126112,\n",
+       "        0.28998899, -0.67546589])"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "logprob = results_7_1.log_likelihood[\"brain_std\"].stack(sample=[\"chain\", \"draw\"]).values.reshape((7, 2000))\n",
+    "\n",
+    "n = logprob.shape[0]\n",
+    "ns = logprob.shape[1] \n",
+    "\n",
+    "lppd = np.zeros(n)\n",
+    "for i in range(n):\n",
+    "    lppd[i] = logsumexp(logprob[i]) - np.log(ns)\n",
+    "lppd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.15 to Code 7.18\n",
+    "\n",
+    "We don't have the means to do these tasks well in Bambi yet."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.19"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "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>speed</th>\n",
+       "      <th>dist</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>4</td>\n",
+       "      <td>2</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>4</td>\n",
+       "      <td>10</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>7</td>\n",
+       "      <td>4</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>7</td>\n",
+       "      <td>22</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>8</td>\n",
+       "      <td>16</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "   speed  dist\n",
+       "0      4     2\n",
+       "1      4    10\n",
+       "2      7     4\n",
+       "3      7    22\n",
+       "4      8    16"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "data = pd.read_csv(\"data/cars.csv\")\n",
+    "data.head()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "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: [dist_sigma, speed, Intercept]\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:03<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 3 seconds.\n"
+     ]
+    }
+   ],
+   "source": [
+    "priors = {\n",
+    "    \"Intercept\": bmb.Prior(\"Normal\", mu=0, sd=100),\n",
+    "    \"mass_std\": bmb.Prior(\"Normal\", mu=0, sd=10),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}\n",
+    "\n",
+    "m = bmb.Model(\"dist ~ speed\", data, priors=priors)\n",
+    "results_m = m.fit()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.20"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "post = results_m.posterior.stack(sample=[\"chain\", \"draw\"])\n",
+    "a = post[\"Intercept\"].values\n",
+    "b = post[\"speed\"].values\n",
+    "sigma = post[\"dist_sigma\"].values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_samples = 1000\n",
+    "n_cases = data.shape[0]\n",
+    "logprob = np.zeros((n_cases, n_samples))\n",
+    "\n",
+    "for s in range(n_samples):\n",
+    "    mu = a[s] + b[s] * data[\"speed\"]\n",
+    "    logprob[:, s] = stats.norm.logpdf(data[\"dist\"], loc=mu, scale=sigma[s])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.21\n",
+    "\n",
+    "We have an array with shape `(50, 1000)`. We apply the `logsumexp()` funcion on `axis=1` to return an array of shape `(50,)`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lppd = logsumexp(logprob, 1) - np.log(n_samples)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.22\n",
+    "\n",
+    "Similar idea here"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pWAIC = np.var(logprob, axis=1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.23"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "421.24274061937064"
+      ]
+     },
+     "execution_count": 23,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "-2 * (np.sum(lppd) - np.sum(pWAIC))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.24"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "16.317728486436113"
+      ]
+     },
+     "execution_count": 24,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "waic_vec = -2 * (lppd - pWAIC)\n",
+    "np.sqrt(n_cases * np.var(waic_vec))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Setup for Code 7.25+\n",
+    "\n",
+    "We need to obtain models `model_6_6`, `model_6_7`, and `model_6_8` from previous chapter."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "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",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>h0</th>\n",
+       "      <td>9.903</td>\n",
+       "      <td>2.003</td>\n",
+       "      <td>6.476</td>\n",
+       "      <td>14.276</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>h1</th>\n",
+       "      <td>14.191</td>\n",
+       "      <td>2.418</td>\n",
+       "      <td>8.693</td>\n",
+       "      <td>18.207</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>treatment</th>\n",
+       "      <td>0.500</td>\n",
+       "      <td>0.503</td>\n",
+       "      <td>0.000</td>\n",
+       "      <td>1.000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>fungus</th>\n",
+       "      <td>0.250</td>\n",
+       "      <td>0.435</td>\n",
+       "      <td>0.000</td>\n",
+       "      <td>1.000</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "             mean     sd  hdi_3%  hdi_97%\n",
+       "h0          9.903  2.003   6.476   14.276\n",
+       "h1         14.191  2.418   8.693   18.207\n",
+       "treatment   0.500  0.503   0.000    1.000\n",
+       "fungus      0.250  0.435   0.000    1.000"
+      ]
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# number of plants\n",
+    "N = 100\n",
+    "\n",
+    "# simulate initial heights\n",
+    "h0 = np.random.normal(10, 2, N)\n",
+    "\n",
+    "# assign treatments and simulate fungus and growth\n",
+    "treatment = np.repeat([0, 1], N / 2)\n",
+    "fungus = np.random.binomial(n=1, p=0.5 - treatment * 0.4, size=N)\n",
+    "h1 = h0 + np.random.normal(5 - 3 * fungus, size=N)\n",
+    "\n",
+    "# compose a clean data frame\n",
+    "d = pd.DataFrame.from_dict({\"h0\": h0, \"h1\": h1, \"treatment\": treatment, \"fungus\": fungus})\n",
+    "az.summary(d.to_dict(orient=\"list\"), kind=\"stats\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "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: [h1_sigma, h0]\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:01<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 2 seconds.\n"
+     ]
+    }
+   ],
+   "source": [
+    "priors = {\n",
+    "    \"h0\": bmb.Prior(\"Lognormal\", mu=0, sigma=0.25),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}\n",
+    "model_6_6 = bmb.Model(\"h1 ~ 0 + h0\", d, priors=priors)\n",
+    "results_6_6 = model_6_6.fit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "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: [h1_sigma, h0:fungus, h0:treatment, h0]\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:02<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 2 seconds.\n"
+     ]
+    }
+   ],
+   "source": [
+    "priors = {\n",
+    "    \"h0\": bmb.Prior(\"Lognormal\", mu=0, sigma=0.2),\n",
+    "    \"h0:treatment\": bmb.Prior(\"Normal\", mu=0, sd=0.5),\n",
+    "    \"h0:fungus\": bmb.Prior(\"Normal\", mu=0, sd=0.5),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}\n",
+    "model_6_7 = bmb.Model(\"h1 ~ 0 + h0 + h0:treatment + h0:fungus\", d, priors=priors)\n",
+    "results_6_7 = model_6_7.fit()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "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: [h1_sigma, h0:treatment, h0]\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:02<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 3 seconds.\n"
+     ]
+    }
+   ],
+   "source": [
+    "priors = {\n",
+    "    \"h0\": bmb.Prior(\"Lognormal\", mu=0, sigma=0.2),\n",
+    "    \"h0:treatment\": bmb.Prior(\"Normal\", mu=0, sd=0.5),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}\n",
+    "model_6_8 = bmb.Model(\"h1 ~ 0 + h0 + h0:treatment\", d, priors=priors)\n",
+    "results_6_8 = model_6_8.fit()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.25"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "Computed from 2000 by 100 log-likelihood matrix\n",
+       "\n",
+       "              Estimate       SE\n",
+       "deviance_waic   321.09    12.22\n",
+       "p_waic            3.73        -\n",
+       "\n",
+       "There has been a warning during the calculation. Please check the results."
+      ]
+     },
+     "execution_count": 29,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "az.waic(results_6_7, scale = \"deviance\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.26"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>rank</th>\n",
+       "      <th>waic</th>\n",
+       "      <th>p_waic</th>\n",
+       "      <th>d_waic</th>\n",
+       "      <th>weight</th>\n",
+       "      <th>se</th>\n",
+       "      <th>dse</th>\n",
+       "      <th>warning</th>\n",
+       "      <th>waic_scale</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>model_6_7</th>\n",
+       "      <td>0</td>\n",
+       "      <td>321.094660</td>\n",
+       "      <td>3.733966</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>1.000000e+00</td>\n",
+       "      <td>12.218020</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>True</td>\n",
+       "      <td>deviance</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>model_6_8</th>\n",
+       "      <td>1</td>\n",
+       "      <td>376.411366</td>\n",
+       "      <td>3.221266</td>\n",
+       "      <td>55.316706</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>14.339417</td>\n",
+       "      <td>11.273896</td>\n",
+       "      <td>True</td>\n",
+       "      <td>deviance</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>model_6_6</th>\n",
+       "      <td>2</td>\n",
+       "      <td>384.832482</td>\n",
+       "      <td>2.265256</td>\n",
+       "      <td>63.737822</td>\n",
+       "      <td>8.302248e-13</td>\n",
+       "      <td>15.524433</td>\n",
+       "      <td>13.298590</td>\n",
+       "      <td>True</td>\n",
+       "      <td>deviance</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "           rank        waic    p_waic     d_waic        weight         se  \\\n",
+       "model_6_7     0  321.094660  3.733966   0.000000  1.000000e+00  12.218020   \n",
+       "model_6_8     1  376.411366  3.221266  55.316706  0.000000e+00  14.339417   \n",
+       "model_6_6     2  384.832482  2.265256  63.737822  8.302248e-13  15.524433   \n",
+       "\n",
+       "                 dse  warning waic_scale  \n",
+       "model_6_7   0.000000     True   deviance  \n",
+       "model_6_8  11.273896     True   deviance  \n",
+       "model_6_6  13.298590     True   deviance  "
+      ]
+     },
+     "execution_count": 30,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "compare_df = az.compare(\n",
+    "    {\n",
+    "        \"model_6_6\": results_6_6,\n",
+    "        \"model_6_7\": results_6_7,\n",
+    "        \"model_6_8\": results_6_8,\n",
+    "    },\n",
+    "    ic=\"waic\",\n",
+    "    scale=\"deviance\",\n",
+    ")\n",
+    "compare_df"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.27"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 31,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array(11.27389615)"
+      ]
+     },
+     "execution_count": 31,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "waic_m_6_7 = az.waic(results_6_7, pointwise=True, scale=\"deviance\")\n",
+    "waic_m_6_8 = az.waic(results_6_8, pointwise=True, scale=\"deviance\")\n",
+    "\n",
+    "# pointwise values are stored in the waic_i attribute.\n",
+    "diff_m_6_7_m_6_8 = waic_m_6_7.waic_i - waic_m_6_8.waic_i\n",
+    "\n",
+    "n = len(diff_m_6_7_m_6_8)\n",
+    "\n",
+    "np.sqrt(n * np.var(diff_m_6_7_m_6_8)).values"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.28\n",
+    "\n",
+    "We use the values from our models, 11.27 for the standard error and 55.35 for the difference between the WAIC of the two models."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 32,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([26.018, 84.622])"
+      ]
+     },
+     "execution_count": 32,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "55.32 + np.array([-1, 1]) * 11.27 * 2.6"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 33,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAE3CAYAAADmP0YBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzKElEQVR4nO3deVSV5d7/8Q/IoMYog6Z1UCl24TxlIg6Iq9LCIU9qJzBTG6yOmGmWA5JUnsw6KJVDRh1Re6AMi8fh8VhZmTaY08HKNDU1JyAEnJDh/v3Rz31CEYEN3Gx4v9ZyLbmn/WF/3fjluq59bwfDMAwBAADAFI5mBwAAAKjPaMYAAABMRDMGAABgIpoxAAAAE9GMAQAAmIhmDAAAwEQ0YwAAACaiGQMAADCRk9kB6ovs7GyzI9g1T09P5eTklNi2fft2SVLnzp3NiASVXheYi5rUTtSldqqJunh7e1/zGEbGYBccHa/8p3rw4EEdPHjQhDS4pLS6wFzUpHaiLrVTbalL7UgBAABQT9GMAQAAmIhmDAAAwEQs4Ifdqi1z/QAA2IJmDHZr6NChZkcAAMBmDC0AAACYiGYMdmvbtm3atm2b2TEAALAJzRjs1q+//qpff/3V7BgAANiEZgwAAMBENGMAAAAmohkDAAAwEbe2gN1ydnY2OwIAADajGYPdGjRokNkRAACwGdOUAAAAJqIZg9365ptv9M0335gdAwAAm9CMwW4dPXpUR48eNTsGAAA2oRkDAAAwEc0YAACAiWjGAAAATMStLWC3XF1dzY4AAIDNaMZgt+655x6zIwAAYDOmKQEAAExEMwa7tXXrVm3dutXsGAAA2IRpStitY8eOmR0BAACbMTIGAABgIpoxAAAAE9GMAQAAmIg1Y7BbjRs3NjsCAAA2oxmD3RowYIDZEQAAsBnTlAAAACaiGYPd2rx5szZv3mx2DAAAbMI0JezWyZMnzY4AAIDNGBkDAAAwEc0YAACAiWjGAAAATMSaMZQpPz9fn3zyiU6ePKmmTZsqPDxcrq6uZsdSfn6+fvrpJ2VnZ6tRo0a1JhcAABVVK0fGZs+eraioqAqdY7FYtH79+mpKVD8lJSUpODhYkZGRevrppxUZGang4GAlJSWZnqtdu3Z69dVXlZiYqMjISLVr1870XAAAVEatbMZqwpdffqkRI0aoQ4cO6tq1q0aNGlWu8z788ENZLJZS/+zevbuaU9ecpKQkRUdHKzs7u8T27OxsRUdHm9b4XMoVFBSklJQU7dmzRykpKQoKCjI1FwAAlVUvm7GNGzdq0qRJGjRokFJTU5WcnKy//vWv5Tp34MCB1vtbXfozaNAg3XjjjWrXrl01J68Z+fn5io2NLfOY559/Xvn5+TUT6P/Lz89XXFycQkJClJqaKhcXF+3bt0/9+/dXamqqQkJC9MILL9R4LgAAbFGhNWNRUVFq3bq1GjVqpA8//FCOjo4aP3687r//fs2ZM0dpaWlyc3PTxIkTNWTIEEnS3r17NWfOHG3fvl0NGzZUv379NH36dLm7u0uSioqKNG/ePH3wwQeSpCFDhqioqKjE4xqGoaVLlyo5OVmnTp1SQECAxo0bp8GDB1f4Gy4qKtILL7ygKVOmaPjw4dbtgYGB5Tq/YcOGatiwofXr8+fP67PPPtO4cePk4OBQ4TxVJSIiosqulZWVdcWI2OV+//13hYWFycfHp8oetyxOTk46efKkMjMz5e/vr3vvvVddunSRJL3yyiuSpNOnTysjI6NGc5klLS3N7AgAgCpS4QX8aWlpeuihh5SSkqJPP/1UL730kr788kv16tVLq1atUmpqqmbMmKGQkBC5ublp7Nixat++vd5//33l5ORo5syZmjZtmhISEiRJiYmJSklJUVxcnCwWi1auXKm0tDS1adPG+pjx8fFav369YmJi1KpVK+3cuVMzZ86Up6en+vbtW6H8e/bs0fHjx+Xs7KyhQ4fq1KlTslgsmjx5soKDgyv6dGjdunU6f/68hg0bVuZxnp6ecnSsvoFIJ6eqey/G5c1wWcdV5eOW5/EkycPDQ05OTtbm91IGDw8PU3KZwdvb2+wIVrUpC/5ATWon6lI71Ya6OBiGYZT34KioKF28eFHJycmS/hix6tGjhzp27KhFixZJkgoKCtSpUyfNmzdPubm5evnll/X555/Lzc1NkvTNN99o1KhR2rBhgwICAhQaGqoHHnhA48ePlyQVFxdrwIAB8vf3V1JSks6dO6fbb79diYmJ6tq1qzXLiy++qEOHDumtt96S9McC/vnz5+uuu+4q83tYs2aNJk2apObNm2vq1Km64YYbtGLFCm3YsEHr1q2Tv79/BZ4+aeTIkfLx8dEbb7xR5nHXGmmqTdauXavIyMhrHrd8+XINHDiwBhL98WJZsWKFIiMjlZKSov79+2vVqlWSZG2EN27cqOHDh9dorvrO29vbrv5t1wfUpHaiLrVTTdSlPM1ehYcPLBaL9e8ODg7y8fEpsc3Z2VkeHh7KysrS4cOHZbFYrI2YJHXq1EmOjo7av3+/mjRpooyMDHXs2NG639HRUe3bt9eJEyckSfv371d+fv4V04AFBQVq0aJFReOruLhYkvTYY49ZG7e4uDht3bpVq1ev1iOPPFLua+3bt087duzQkiVLKpyjNgsPD7/mP9AmTZooPDy8BlP9kcvX11fx8fHq06dPiX0FBQWKj4+Xn59fjecCAMAWFW7GLp/+cXBwKHXbtQbcyru+6tJ1Fi5cqObNm5eZpTz8/PwklVwj5uTkpICAAB0/frxC10pOTtb111+vXr16VThHbebq6qrY2FhFR0df9ZhZs2bV+H29XF1dNXPmTEVHR2vo0KHq27ev/vKXv2jjxo2Kj4/Xli1bNH/+fO43BgCwK9W6sCYwMFCrVq3SmTNnrKNjO3bsUHFxsQIDA+Xu7i4/Pz/t2rVLPXr0kPRH87V7927rdGFgYKBcXFx07Ngx6zG2aNu2rVxcXHTw4EHrtGdxcbGOHDmi0NDQcl8nPz9fH3/8saKioqp1LZhZLt3nLTY2tsQImbe3t2JjYyt8H7iqzhUXF6ctW7ZYt/v6+mr+/Pmm5QIAoLKqtRmLiIjQggULNHXqVE2YMEG5ubmKiYnRHXfcoYCAAEnSqFGjtGTJErVs2VJBQUFauXKlMjIyrM2Ym5ubxowZo7lz58owDHXr1k3nzp3Tzp075ejoqBEjRlQok5ubm0aOHKmEhAQ1a9ZMLVq00IoVK5STk1Ohd2euX79eeXl511y4b8+ioqI0fPjwWncH/tqaCwCAyqjWZqxRo0Z6++239dJLL+m+++6Tq6urwsPDNX36dOsxY8aMUWZmpmbMmCFJGjx4sCIiInTgwAHrMRMnTpSvr68SExMVGxsrNzc33XrrrRo3blylcj3zzDNydnbWs88+q/Pnz6tNmzZatmxZhRbvv//++woNDb1i6rSucXV1rZWL4V1dXdWoUSO1bNlSYWFhZscBAKDSKvRuSlQe76KxTWlvKLj83ZSoebxDrPahJrUTdamdasu7KeveYicAAAA7UufujLlo0SItXry41H1dunTR0qVLyzw/Jibmqnc3j4iI0OzZs23OCAAAcEmda8ZGjhypAQMGlLrvzx9jdDXR0dEaO3Zsqfv+fL80AACAqlDnmjEvLy95eXlV+nwfH586/7mGdUWTJk3MjgAAgM3qXDOG+oN3UQIA6gIW8AMAAJiIZgx2a+PGjdq4caPZMQAAsAnTlLBbOTk5ZkcAAMBmjIwBAACYiGYMAADARDRjAAAAJmLNGOyWr6+v2REAALAZzRjsVp8+fcyOAACAzZimBAAAMBHNGOzWhg0btGHDBrNjAABgE6YpYbfy8vLMjgAAgM0YGQMAADARzRgAAICJaMYAAABMxJox2K2mTZuaHQEAAJvRjMFuhYaGmh0BAACbMU0JAABgIpox2K1169Zp3bp1ZscAAMAmTFPCbp07d87sCAAA2IyRMQAAABPRjAEAAJiIZgwAAMBErBmD3br++uvNjgAAgM1oxmC3QkJCzI4AAIDNmKYEAAAwEc0Y7NaaNWu0Zs0as2MAAGATpilhty5cuGB2BAAAbMbIGAAAgIloxgAAAExEMwYAAGAi1ozBbt1www1mRwAAwGY0Y7Bb3bt3NzsCAAA2Y5oSAADARDRjsFsff/yxPv74Y7NjAABgE6YpYbcKCgrMjgAAgM0YGQMAADARzRgAAICJaMYAAABMxJox2K2//OUvZkcAAMBmNGOwW926dTM7AgAANmOaEgAAwEQ0Y7BbqampSk1NNTsGAAA2YZoSdqu4uNjsCAAA2IyRMQAAABPRjAEAAJiIZgwAgFIcOXJER44cMTsG6gHWjMFutWrVyuwIAOqovLw8rV69WpI0evRoubu7mxsIdRojY7BbnTt3VufOnc2OAaAO2rZtmwzDkGEY2rZtm9lxUMfRjAEA8Cd5eXlKT0+3fp2enq68vDwTE6GuY5oSdmvVqlWSpGHDhpmcBKh/Lr3+6qLTp0/LMAzr14ZhKCUlRV5eXpW+ppOTkwoLC6sgXdXjZ6j5aMZqiKenpxwdGYi0hbe3d4mvnZycSt2OmsXzX/vURE0uvf7qmsLCQp09e/aK7WfPnpWXl5dN33dtfc7q+2u4Nnz/tfNfRh2Uk5NjdgS75u3trezs7BLbLv2Wefl21JzS6gJz1VRNBg8eXO2PYYbPPvtMv/32W6n7mjRporCwsEpdtza/VmprrppQE3UpT7PHUA0AALpyrdjlWDuG6kIzBgCApN27d5dYK3Y5wzC0e/fuGkyE+oJpStitm266yewIAOoQi8WiEydOXPMYoKrRjMFudejQwewIAOoQX19f3lkIUzBNCbtVVFSkoqIis2MAAGATmjHYrdWrV1s/rgQAAHtFMwYAAGAimjEAAAAT0YwBAACYiGYMAADARNzaAnaL+/0AAOoCmjHYrbZt25odAQAAmzFNCbt18eJFXbx40ewYAADYhGYMdistLU1paWlmxwAAwCZMUwIAapX8/Hx98sknOnnypJo2barw8HC5urrWuwyoP2plMzZ79mzt27dPSUlJ5T7HYrFo/vz5uuuuu6oxGQCgOiUlJSk2NlbZ2dnWbd7e3oqNjVVUVFSNZYiLi1NmZqZ1m6+vr2bOnFljGVC/1Ntpyi+//FIjRoxQhw4d1LVrV40aNarc5+7evVujR49W165d1bVrVz344IPavXt3NaYFgLovKSlJ0dHRJRoxScrOzlZ0dHSFfkG3NUNQUJBSUlK0Z88epaSkKCgoqMYyoP6pl83Yxo0bNWnSJA0aNEipqalKTk7WX//613Kde/bsWT388MPy9/dXcnKykpOT5efnp7Fjx+rMmTPVnBwA6qb8/HzFxsaWeczzzz+v/Pz8as0QFxenkJAQpaamqn///rr++uvVv39/paamKiQkRC+88EK1ZkD9VKFpyqioKLVu3VqNGjXShx9+KEdHR40fP17333+/5syZo7S0NLm5uWnixIkaMmSIJGnv3r2aM2eOtm/froYNG6pfv36aPn263N3dJUlFRUWaN2+ePvjgA0nSkCFDVFRUVOJxDcPQ0qVLlZycrFOnTikgIEDjxo3T4MGDK/wNFxUV6YUXXtCUKVM0fPhw6/bAwMBynX/gwAGdPn1af//733XjjTdKkqKjo5WWlqaDBw+qXbt2Fc6Eyrn11lvNjgDYtYiICLMjWGVlZV0xIna533//XWFhYfLx8am2DJmZmfL399e99957xf7Tp08rIyOjUhmcnJxUWFhYVVFrFd5IZbsKrxlLS0vTQw89pJSUFH366ad66aWX9OWXX6pXr15atWqVUlNTNWPGDIWEhMjNzU1jx45V+/bt9f777ysnJ0czZ87UtGnTlJCQIElKTExUSkqK4uLiZLFYtHLlSqWlpalNmzbWx4yPj9f69esVExOjVq1aaefOnZo5c6Y8PT3Vt2/fCuXfs2ePjh8/LmdnZw0dOlSnTp2SxWLR5MmTFRwcfM3zW7VqpSZNmuiDDz7QE088IUl6//331bx5c918881XPc/T01OOjvVyILLKeHt7l/i6Z8+eJiXBn11eF5ivvDVxcqo9y4Yv/yW8rOOqK/elDB4eHqU+hoeHh00ZatPzXZXs/WdAbcjvYBiGUd6Do6KidPHiRSUnJ0v6Y8SqR48e6tixoxYtWiRJKigoUKdOnTRv3jzl5ubq5Zdf1ueffy43NzdJ0jfffKNRo0Zpw4YNCggIUGhoqB544AGNHz9eklRcXKwBAwbI399fSUlJOnfunG6//XYlJiaqa9eu1iwvvviiDh06pLfeektS+Rfwr1mzRpMmTVLz5s01depU3XDDDVqxYoU2bNigdevWyd/f/5rPw/79+/XEE0/o8OHDkqQWLVrorbfeUqtWra56zrV+40PZvL29r3gOL1y4IElq2LChGZGg0usCc9lrTdauXavIyMhrHrd8+XINHDiwWjOkpKSof//+V+zfuHGjhg8fXqkM9lqXuq4m6lKeZq/CQzV//ggaBwcH+fj4lNjm7OwsDw8PZWVl6ZdffpHFYrE2YpLUqVMnOTo6av/+/crLy1NGRoY6duz430COjmrfvr316/379ys/P1/jxo1Tp06drH/ee+89azNUEcXFxZKkxx57THfddZfatm2ruLg4ubu7a/Xq1dc8/8KFC5o2bZo6dOig5ORkvffeewoODtbjjz+uc+fOVTgPKm/NmjVas2aN2TEAVIHw8PBr/qfVpEkThYeHV2sGX19fxcfHq6CgoMS+goICxcfHy8/Pr1ozoH6q8Jjp5cOsDg4OpW671oCbg4NDuR7v0nUWLlyo5s2bl5mlPPz8/CSVXCPm5OSkgIAAHT9+/Jrnp6Wl6fDhw3rvvffUoEEDSdK8efN022236d///nel1rEBQH3n6uqq2NhYRUdHX/WYWbNmVeu9vlxdXTVz5kxFR0dr6NChmjhxooKDg/XDDz8oPj5eW7Zs0fz587nfGKpctU5gBwYGatWqVTpz5ox1dGzHjh0qLi5WYGCg3N3d5efnp127dqlHjx6S/mi+du/ebZ0uDAwMlIuLi44dO2Y9xhZt27aVi4uLDh48aJ32LC4u1pEjRxQaGnrN8y9cuCAHB4cS678cHR3L1YACAK7u0j28zLzP2KXHiIuLK/EmL19fX82fP5/7jKFaVGszFhERoQULFmjq1KmaMGGCcnNzFRMTozvuuEMBAQGSpFGjRmnJkiVq2bKlgoKCtHLlSmVkZFibMTc3N40ZM0Zz586VYRjq1q2bzp07p507d8rR0VEjRoyoUCY3NzeNHDlSCQkJatasmVq0aKEVK1YoJyenXKNaISEhmjt3rmJjYzVq1CgVFxdryZIlcnR01O23317xJwkAYBUVFaXhw4ebevf72pAB9Uu1NmONGjXS22+/rZdeekn33XefXF1dFR4erunTp1uPGTNmjDIzMzVjxgxJ0uDBgxUREaEDBw5Yj5k4caJ8fX2VmJio2NhYubm56dZbb9W4ceMqleuZZ56Rs7Oznn32WZ0/f15t2rTRsmXLyrV4PzAwUIsWLdLrr7+uESNGyMHBQbfccoveeustNWvWrFJ5AAD/5erqWm2L9O0pA+qPCr2bEpXHu2hsU9o7Xvbu3Sup5JtKULN4h1jtQ01qJ+pSO9WWd1PWzZueoF6gCQMA1AV1rhlbtGiRFi9eXOq+Ll26aOnSpWWeHxMTc9W7CUdERGj27Nk2Z0TVOHv2rCTpuuuuMzkJAACVV+emKU+fPq2cnJxS9zVs2FBNmzYt8/ysrKyrfsakm5tbpT+Gg+Fp25Q2lLxq1SpJ0rBhw8yIBDH1UhtRk9qJutROTFNWEy8vL3l5eVX6fB8fn2r73DMAAIDL8WGJAAAAJqIZAwAAMBHNGAAAgInq3Jox1B9//kB5AADsFc0Y7NbNN99sdgQAAGzGNCXsVl5envLy8syOAQCATWjGYLc2bNigDRs2mB0DAACb0IwBAACYiGYMAADARDRjAAAAJqIZAwAAMBG3toDd6tSpk9kRAACwGc0Y7Fbr1q3NjgAAgM2YpoTdOn36tE6fPm12DAAAbEIzBrv1ySef6JNPPjE7BgAANqEZAwAAMBHNGAAAgIloxgAAAExEMwYAAGAibm0Bu9WlSxezIwAAYDOaMditli1bmh0BAACbMU0Ju5WVlaWsrCyzYwAAYBOaMditTZs2adOmTWbHAADAJjRjAAAAJqIZAwAAMBHNGAAAgIloxgAAAEzErS1gt2677TazIwAAYDOaMditG2+80ewIAADYjGlK2K2MjAxlZGSYHQMAAJvQjMFuffHFF/riiy/MjgEAgE1oxgAAAExEMwYAAGAimjEAAAAT0YwBAACYiFtbwG716NHD7AgAANiMZgx2q3nz5mZHAADAZkxTwm6dOHFCJ06cMDsGAAA2oRmD3frqq6/01VdfmR0DAACb0IwBAACYiGYMAADARDRjAAAAJqIZAwAAMBG3toDd6tmzp9kRAACwGc0Y7FazZs3MjgAAgM2YpoTdOnbsmI4dO2Z2DAAAbEIzBru1detWbd261ewYAADYhGYMAAA7d+TIER05csTsGKgkmjEAAOxYXl6eVq9erdWrVysvL8/sOKgEmjEAAOzYtm3bZBiGDMPQtm3bzI6DSqAZAwDATuXl5Sk9Pd36dXp6OqNjdohbW8Bu9e7d2+wIAEy0atUqsyOUm5OTkwoLC6v8uqdPn5ZhGNavDcNQSkqKvLy8qvyx7MmwYcPMjlAhNGM1xNPTU46ODETawtvbu8yvYQ7qUPvUl5o4OdnXf2FVnbewsFBnz569YvvZs2fl5eVld89PVarIa6A2vF4cjD+31Kg22dnZZkewa97e3lc8h5feOXTjjTeaEQkqvS4wFzWpnaqjLp999pn+85//lLqvXbt2CgsLq9LHq4tq4vVSnmaPoRrYrW+//Vbffvut2TEAoMZdvlbscqwdsy80YwAA2Jndu3errIktwzC0e/fuGkwEW9TfCWUAAOyUxWLRiRMnrnkM7APNGAAAdsbX19fu3jGIq2OaEgAAwESMjMFu9e3b1+wIAADYjGYMdsvHx8fsCAAA2IxpStitQ4cO6dChQ2bHAADAJoyMwW59//33kqSWLVuaGwQAABswMgYAAGAimjEAAAAT0YwBAACYiGYMAADARCzgh90KDw83OwIAADarlSNjs2fPVlRUVIXOsVgsWr9+fTUlQm3k5eUlLy8vs2MAwDXl5+dr7dq1euedd7R27Vrl5+fbxbVRM+rtyNiXX36p119/XT/99JOcnZ0VHBysZcuWlfv8jz76SImJiTpw4IAaN26sPn36aO7cudWYGJc7cOCAJKl169YmJwGAq0tKStLs2bOVlZVl3ebt7a3Y2NgKDzyUdu24uDhlZmZat/n6+mrmzJk2Xxs1p142Yxs3btRzzz2niRMnas6cOTIMQ3v27Cn3+cuWLdOSJUs0ZcoUdezYURcuXODmoybYsWOHJJoxALVXUlKSoqOjr9ienZ1t3V7ZpunStUNCQvTmm2+qTZs22rNnj+Lj422+NmqWg2EYRnkPjoqKUuvWrdWoUSN9+OGHcnR01Pjx43X//fdrzpw5SktLk5ubmyZOnKghQ4ZIkvbu3as5c+Zo+/btatiwofr166fp06fL3d1dklRUVKR58+bpgw8+kCQNGTJEFy9e1IEDB5SUlCRJMgxDS5cuVXJysk6dOqWAgACNGzdOgwcPtmazWCyaP3++7rrrrjK/h6KiIoWHh+vxxx/X8OHDK/RkSVJubq569+6tN954Qz179iz3ednZ2RV+LPyXt7f3Fc/hqlWrJEnDhg0zIxJUel1gLmpSe+Tn5ys4OLjMejRp0kR79uyRq6trha/drl07BQUFKTU1Vc7OztZ9BQUFGjp0qPbt26fdu3dX+Nr1SU28Xry9va95TIVHxtLS0vTQQw8pJSVFn376qV566SV9+eWX6tWrl1atWqXU1FTNmDFDISEhcnNz09ixY9W+fXu9//77ysnJ0cyZMzVt2jQlJCRIkhITE5WSkqK4uDhZLBatXLlSaWlpatOmjfUx4+PjtX79esXExKhVq1bauXOnZs6cKU9Pzwp/WPSePXt0/PhxOTs7a+jQoTp16pQsFosmT56s4ODga56/efNmFRUVKTMzUwMHDlReXp7at2+vZ599VjfeeGOFsgAASoqIiDA7QpXJysq65n/0v//+u8LCwir8WbtZWVnKzMyUv7+/7r333iv2nz59WhkZGZW6dn3i5OSk1NRUs2NUvBm7+eab9fe//12S9NBDD2nJkiVycnLSgw8+KEl64okntHTpUm3fvl25ubk6f/685s6dKzc3N0l/LM4fNWqUfv31VwUEBOhf//qXxo0bp4EDB0qSpk+frs2bN1sf79y5c3rnnXeUmJiorl27SpJuvPFG7d69WytWrKhwM3bkyBFJ0oIFCzR16lTdcMMNWrFihaKiorRu3Tr5+/uXef7Ro0dlGIYWLlyoadOmycvLS2+88YZGjRqltWvXqlGjRqWe5+npKUfHWvl+Cbtx+W8XTk5OpW5HzeL5r33suSaXXtd1QVFRUbmPq+j3fenaHh4epZ7r4eFR6WvXN7Xh9VLhClksFuvfHRwc5OPjU2Kbs7OzPDw8lJWVpcOHD8tisVgbMUnq1KmTHB0dtX//fjVp0kQZGRnq2LGjdb+jo6Pat2+vEydOSJL279+v/Px8jRs3Tg4ODtbjCgoK1KJFi4rGV3FxsSTpscces05pxsXFaevWrVq9erUeeeSRa55fUFCgGTNmKDQ0VJI0b9489ezZU5999pm1qbxcTk5OhbPiv0obSi4sLJTEFLCZmBKrfey9JrVhlKKqrF27VpGRkdc8btasWVf9v+Na1540aZL69+9/xf6NGzdq+PDhlbp2fWK305SXd9gODg6lbrvWUrQ/N1ZluXSdhQsXqnnz5mVmKQ8/Pz9JUmBgYInrBAQE6Pjx4+U+/6abbrJuc3d3l7+/v44dO1bhPKi8O+64w+wIAHBV4eHh1/zPvkmTJpW6Z2J4eLh8fX0VHx+vPn36XLFmLD4+Xn5+ftyP0U5U67xZYGCgfv75Z505c8a6bceOHSouLlZgYKDc3d3l5+enXbt2WfcbhqHdu3eXuIaLi4uOHTumgICAEn8qMzLWtm1bubi46ODBg9ZtxcXFOnLkyBXNXmk6d+4sSSXOP3v2rDIyMsp1PqqOu7u79Y0gAFDbuLq6KjY2tsxjZs2aVakF9q6urpo5c6a2bNmioUOHauPGjTp27Jg2btyooUOHasuWLZoxYwaL9+1EtU4kR0REWNdmTZgwQbm5uYqJidEdd9yhgIAASdKoUaO0ZMkStWzZUkFBQVq5cqUyMjKsa7fc3Nw0ZswYzZ07V4ZhqFu3bjp37px27twpR0dHjRgxokKZ3NzcNHLkSCUkJKhZs2Zq0aKFVqxYoZycnBLvzryaVq1aKTw8XC+++KKef/55eXp6asGCBfLx8VFYWFjFnyRU2r59+yT9sY4RAGqjS7eWqI77jF06Ny4ursTdAXx9fTV//nxua2FHqrUZa9Sokd5++2299NJLuu++++Tq6qrw8HBNnz7desyYMWOUmZmpGTNmSJIGDx6siIgI6w09JWnixIny9fVVYmKiYmNj5ebmpltvvVXjxo2rVK5nnnlGzs7OevbZZ3X+/Hm1adNGy5Ytu+bi/Uvmzp2rf/zjHxo/frwMw1CXLl307rvvXnXxPqrHpRFUmjEAtVlUVJQeffRRffDBBzp58qSaNm2q8PDwKhm1ioqK0vDhw/XJJ59U+bVRcyp0nzFUnj0vqK0NuM9Y7WTvi8XrImpSO1GX2qm2LODnXgsAAAAmqnM3H1m0aJEWL15c6r4uXbpo6dKlZZ4fExOjtLS0UvdFRERo9uzZNmcEAAC4pM41YyNHjtSAAQNK3dewYcNrnh8dHa2xY8eWuu/P90sDAACoCnWuGfPy8pKXl1elz/fx8eGjI+zEtT6HFAAAe1DnmjHUH9ddd53ZEQAAsBkL+GG39u7dq71795odAwAAmzAyBruVnp4uqeTnpQIAYG8YGQMAADARzRgAAICJaMYAAABMRDMGAABgIhbww27dfffdZkcAAMBmNGOwW+X5RAUAAGo7pilht3744Qf98MMPZscAAMAmNGOwWz/++KN+/PFHs2MAAGATmjEAAAAT0YwBAACYiGYMAADARDRjAAAAJuLWFrBbERERZkcAAMBmNGOwWy4uLmZHAADAZkxTwm6lp6crPT3d7BgAANiEZgx2a+/evdq7d6/ZMQAAsAnNGAAAgIloxgAAAExEMwYAAGAimjEAAAATcWsL2K0hQ4aYHQEAAJvRjMFuNWjQwOwIAADYjGlK2K1du3Zp165dZscAAMAmNGOwW/v379f+/fvNjgEAgE1oxgAAAExEMwYAAGAimjEAAAAT0YwBAACYyMEwDMPsEAAAAPUVI2MAAAAmohkDAAAwEc0YAACAiWjGAAAATEQzBgAAYCKaMZhixYoVioiIUOfOndW5c2eNGDFCmzZtkiQVFBTolVdeUUREhDp27KjQ0FA9/fTTOnbsWIlrXLx4UXFxcerevbs6duyoxx57TCdOnDDhu6k7yqrL5WJiYmSxWPT222+X2E5dqlZ5anLw4EE9+eST6tq1qzp06KChQ4fql19+se6nJlXvWnU5e/as4uLi1Lt3b7Vv31533nmn3n333RLXoC7Va/HixbJYLJo9e7Z1m2EYSkhIUGhoqNq3b6+oqCjt27evxHk5OTmaMmWKunTpoi5dumjKlCnKzc2t3rAGYIJ///vfxqZNm4xDhw4ZBw4cMF577TUjODjY+PHHH43c3Fxj9OjRxpo1a4xffvnF2LVrl3H//fcbAwYMMAoKCqzXiImJMXr27Gls3rzZSE9PNyIjI41BgwYZhYWFJn5n9q2suvzZunXrjMGDBxuhoaHG0qVLS+yjLlXrWjU5fPiw0b17d2POnDlGenq6cfjwYWPTpk3GsWPHrNegJlXvWnWZMWOG0a9fP2Pr1q3GkSNHjNTUVKNNmzZGamqq9RrUpfrs2LHDCAsLMyIiIoznn3/eun3x4sVGx44djfXr1xt79+41JkyYYPTs2dPIy8uzHjN27Fhj4MCBxvbt243t27cbAwcONB599NFqzUszhlqjW7duxnvvvVfqvn379hlBQUHGTz/9ZBiGYeTm5hpt2rQxPvroI+sxx44dMywWi/HFF1/USN764vK6HD161AgNDTX2799vhIWFlWjGqEvN+HNNJk2aZEyaNOmqx1KTmvPnutx9993G/PnzS+x/4IEHrI0Bdak+ubm5Rnh4uLF161YjMjLS+pwXFxcbPXv2NN58803rsefPnzc6duxordv+/fuNoKAgY9u2bdZjvvvuOyMoKMj45Zdfqi0z05QwXVFRkdasWaNz586pU6dOpR5z5swZSZKnp6ckKT09XQUFBQoNDbUec/311yswMFA7duyo/tD1QGl1KSws1NNPP63x48crMDDwinOoS/W6vCbFxcX69NNPddNNN2ns2LG6/fbbNWzYMK1du9Z6DjWpfqW9Vjp37qzPPvtMx48flyRt375dP/74o3r16iWJulSnmTNn6s4779Ttt99eYvvRo0eVkZGhnj17Wrc1bNhQ3bp1sz7nO3bsUOPGjdW5c2frMV26dFHjxo2rtS5O1XZl4Br27t2rkSNHKj8/X40bN9brr78ui8VyxXEXL17UP/7xD4WFhalZs2aSpMzMTDVo0EDe3t4ljvXx8VFmZmaN5K+ryqpLQkKCvLy89Le//a3Uc6lL9bhaTTIyMnTu3DktXrxY0dHRmjx5sr7++mtNnjxZjRs3Vt++falJNSrrtTJjxgzNmjVLffv2lZOTk3VbWFiYJF4r1SUlJUWHDx/WK6+8csW+jIwMSZKvr2+J7T4+Pjp16pSkP+rSpEkTOTg4WPc7ODioSZMm1VoXmjGYplWrVlq9erXy8vL0f//3f5o6daqSkpIUFBRkPaawsFBTpkxRXl6eFi5caGLa+uNqdcnOztaHH36ojz76yOyI9c7VanJppDg8PFwPPfSQJOnWW29Venq6li9frr59+5qYuu4r62fY8uXLtX37di1cuFDNmzfXtm3bNHfuXLVo0UK9e/c2O3qddODAAb322mtauXKlnJ2dzY5TIUxTwjQuLi4KCAhQ27Zt9fTTT+vWW28t8W6jwsJCTZo0SXv37tW7775b4jdIX19fFRUVKTs7u8Q1s7KyrvitBxVztbp8++23ysjIUGhoqIKDgxUcHKzffvtN8+bNs/7nQl2qx9Vq4u3tLScnpyumjFu3bm2dHqMm1edqdblw4YJee+01TZkyRf369dMtt9yiyMhIDRw4UImJiZKoS3XYuXOnsrOzdc8991h/Rn377bdauXKlgoOD5eXlJUlXjHD9+Tn39fXV77//LuNPH9ttGIZ+//33aq0LzRhqjeLiYl28eFHSH7e3eOqpp7R3714tW7ZMfn5+JY5t27atnJ2d9dVXX1m3nThxQr/88stV152hci7V5W9/+5s+/vhjrV692vrH399fo0ePtjbR1KVmXKqJi4uL2rVrp4MHD5bYf+jQITVv3lwSNalJl+pSWFiogoICNWjQoMT+Bg0aqLi4WBJ1qQ79+/dXWlpaiZ9Rbdu21d13363Vq1erVatW8vPz05YtW6zn5Ofna9u2bdbnvFOnTjp37lyJ9WE7duwoc01zVWCaEqaYN2+e+vbtq2bNmuns2bP63//9X3377bdavHixCgsLFR0drf/85z9atGiRHBwcrHP97u7uatiwodzd3TVs2DC98sor8vHxkZeXl+bMmSOLxaKQkBCTvzv7VVZdfHx85OPjU+J4Z2dn+fr6qnXr1pJEXapBWTWRpHHjxmnixInq2rWrbr/9dn3zzTdau3at3njjDUnUpLqUVRc3NzfddtttevXVV9W4cWM1b95c3333nVavXq0pU6ZIoi7VwcPDQx4eHiW2NW7cWJ6entblL6NGjdLixYvVunVrtWzZUgsXLlTjxo11zz33SJICAwPVq1cvzZo1y3p/slmzZiksLMz6c6460IzBFJmZmZoyZYoyMjLk7u4ui8Wit956S7169dLRo0f1ySefSJLuvffeEufNmTPHum369OlycnLSU089pQsXLqhHjx6aO3fuFb+NovzKqkt5UZeqda2a9O/fX7Nnz9bixYv14osvKiAgQC+//HKJ9WLUpOpdqy6vvfaaXnvtNU2ePFk5OTlq3ry5oqOjFRkZab0Gdal5Dz/8sPLz8zV79mzl5OSoQ4cOSkxMlJubm/WYV199VXFxcRo7dqwkqV+/foqJianWXA7GnydGAQAAUKNYMwYAAGAimjEAAAAT0YwBAACYiGYMAADARDRjAAAAJqIZAwAAMBHNGAAAgIloxgAAAExEMwagzkpISJDFYpHFYtEtt9yibt26adiwYfrnP/9p/YitqmaxWLR8+fJquTaAuomPQwJQp7m7u2vp0qWSpLy8PP3www967733lJycrKVLl6pt27ZV+njJycm64YYbqvSaAOo2Pg4JQJ2VkJCg5cuX65tvvimxPTc3Vw888IDy8/O1bt06PgsQgKmYpgRQ73h4eGjKlCn69ddf9dVXX0mS8vPzNXfuXPXp00dt27bVoEGD9Pnnn1vPefbZZzVs2LArrrVixQp16NBBZ86ckXTlNOWmTZv00EMPqUePHurcubOGDx+uzZs3l7hGQkKCunfvrh9++EHDhw9Xhw4dNGTIEG3btu2Kx0tJSVFERITatWunkJAQTZgwQXl5edb927ZtU2RkpDp06KDu3btrxowZ1mwAaieaMQD1Uvfu3eXk5KRdu3ZJkiZMmKDU1FQ9+uijWrRokdq1a6fx48frxx9/lCQNGDBA6enpOnLkSInrrF27Vn369JGbm1upj3P06FGFhYVp7ty5SkhIUKdOnfTwww/r+++/L3HchQsXNHXqVI0YMUILFiyQi4uLnnzySZ0/f956zJtvvqmYmBh169ZNb7zxhmJjY+Xm5qZz585Jkr7//nuNHj1avr6+WrBggZ577jl9/vnnmjZtWpU9bwCqHmvGANRLrq6u8vb2VmZmprZu3apNmzYpKSlJt912myQpNDRUhw4d0sKFC7VgwQL17NlTXl5eWrdunR555BFJ0smTJ/X9998rPj7+qo8TGRlp/XtxcbG6d++u/fv364MPPlCXLl2s+y5cuKBp06apR48ekiR/f38NGTJE3333nXr37q3c3FwtXrxYDz74oJ577jnreXfccYf176+++qo6depUIk/Tpk01evRo/fzzzwoKCrLpOQNQPRgZA1BvXVoyu2XLFvn5+alz584qLCy0/unRo4fS09MlSU5OTrrjjju0du1a6/nr1q1To0aN1Ldv36s+xokTJzR16lT16tVLwcHBatOmjTZv3qxDhw6VOM7Z2Vndu3e3fh0YGCjpj4ZPknbs2KELFy7o3nvvLfVxzp8/r507d2rAgAElvocuXbrI2dlZe/bsqfDzA6BmMDIGoF7Kz8/X6dOn5evrq5MnTyojI0Nt2rS54rg/L+4fOHCgUlJSdPDgQbVq1Urr1q1Tv3791LBhw1Ifo7i4WOPHj9fZs2c1YcIEBQQEqFGjRlqwYIGysrJKHHvdddfJ0fG/vx+7uLhYc0rS6dOnJUl+fn6lPlZubq6Kior0/PPP6/nnn79i//Hjx8t4NgCYiWYMQL309ddfq7CwUB07dtTXX3+tpk2b6o033ijznNtuu02+vr5au3athgwZop07d1qnLEvz66+/6ocfftBbb72l3r17W7dfuHChwnm9vLwkSRkZGWrSpMkV+93d3eXg4KAnn3xSffr0uWK/v79/hR8TQM2gGQNQ7+Tm5mrevHkKCAhQSEiIHBwc9M4776hx48bW6cHSNGjQQHfddZfWrVsnV1dXeXh4qFevXlc9/tKo1qVRLkn67bfftGPHjgqv3+rUqZMaNmyo1atXa+rUqVfsb9y4sTp27KiDBw/qySefrNC1AZiLZgxAnVZUVKSdO3dKks6ePas9e/bovffe0/nz57V06VI1aNBAPXv2VGhoqMaMGaOHH35YN910k86cOaOffvpJ+fn5evrpp63XGzBggJYvX653331X/fv3L9FoXa5169Zq1qyZXn75ZUVHR+vs2bNasGBBpUapPDw89Pjjj+uf//ynCgoK1Lt3b128eFGff/65nnzySTVt2lSTJ0/W6NGj5ejoqDvvvFPXXXedjh8/rk2bNumpp55Sq1atKvy4AKofzRiAOi0vL08jRoyQg4OD3Nzc9Je//EWDBg1SZGSkdf2Vg4ODXn/9dS1atEj/+te/dPz4cXl6euqWW25RVFRUiet16dJF119/vY4fP66BAweW+dguLi5KSEjQ7NmzNWHCBDVr1kyPPfaYvv32W/38888V/l4effRReXp6atmyZfqf//kfeXp6qmvXrrruuuskSV27dtWKFSu0YMECPfPMMyouLlbz5s3Vq1cv+fr6VvjxANQM7sAPAABgIm5tAQAAYCKaMQAAABPRjAEAAJiIZgwAAMBENGMAAAAmohkDAAAwEc0YAACAiWjGAAAATEQzBgAAYCKaMQAAABPRjAEAAJjo/wEmHrfFwM5pcwAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 600x300 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "az.plot_compare(compare_df);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.30"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "array(6.39613523)"
+      ]
+     },
+     "execution_count": 34,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "waic_m_6_6 = az.waic(results_6_6, pointwise=True, scale=\"deviance\")\n",
+    "diff_m6_6_m6_8 = waic_m_6_6.waic_i - waic_m_6_8.waic_i\n",
+    "n = len(diff_m6_6_m6_8)\n",
+    "np.sqrt(n * np.var(diff_m6_6_m6_8)).values"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.31\n",
+    "\n",
+    "dSE is calculated by compare above, but rethinking produces a pairwise comparison. This is not implemented in arviz, but we can hack it together:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n",
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>model_6_7</th>\n",
+       "      <th>model_6_8</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>model_6_8</th>\n",
+       "      <td>11.273896</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>model_6_6</th>\n",
+       "      <td>13.298590</td>\n",
+       "      <td>6.396135</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "           model_6_7  model_6_8\n",
+       "model_6_8  11.273896   0.000000\n",
+       "model_6_6  13.298590   6.396135"
+      ]
+     },
+     "execution_count": 35,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "dataset_dict = {\"model_6_6\": results_6_6, \"model_6_7\": results_6_7, \"model_6_8\": results_6_8}\n",
+    "\n",
+    "# compare all models\n",
+    "s0 = az.compare(dataset_dict, ic=\"waic\", scale=\"deviance\")[\"dse\"]\n",
+    "# the output compares each model to the 'best' model - i.e. two models are compared to one.\n",
+    "# to complete a pair-wise comparison we need to compare the remaining two models.\n",
+    "\n",
+    "# to do this, remove the 'best' model from the input data\n",
+    "del dataset_dict[s0.index[0]]\n",
+    "\n",
+    "# re-run compare with the remaining two models\n",
+    "s1 = az.compare(dataset_dict, ic=\"waic\", scale=\"deviance\")[\"dse\"]\n",
+    "\n",
+    "# s0 compares two models to one model, and s1 compares the remaining two models to each other\n",
+    "# now we just nee to wrangle them together!\n",
+    "\n",
+    "# convert them both to dataframes, setting the name to the 'best' model in each `compare` output.\n",
+    "# (i.e. the name is the model that others are compared to)\n",
+    "df_0 = s0.to_frame(name=s0.index[0])\n",
+    "df_1 = s1.to_frame(name=s1.index[0])\n",
+    "\n",
+    "# merge these dataframes to create a pairwise comparison\n",
+    "pd.merge(df_0, df_1, left_index=True, right_index=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.32"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def standardize(x):\n",
+    "    return (x - np.mean(x)) / np.std(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = pd.read_csv(\"data/WaffleDivorce.csv\")\n",
+    "d[\"A\"] = standardize(d[\"MedianAgeMarriage\"])\n",
+    "d[\"D\"] = standardize(d[\"Divorce\"])\n",
+    "d[\"M\"] = standardize(d[\"Marriage\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "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: [D_sigma, A, Intercept]\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:01<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 2 seconds.\n",
+      "Auto-assigning NUTS sampler...\n",
+      "Initializing NUTS using jitter+adapt_diag...\n",
+      "Multiprocess sampling (2 chains in 2 jobs)\n",
+      "NUTS: [D_sigma, M, Intercept]\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:01<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 2 seconds.\n",
+      "Auto-assigning NUTS sampler...\n",
+      "Initializing NUTS using jitter+adapt_diag...\n",
+      "Multiprocess sampling (2 chains in 2 jobs)\n",
+      "NUTS: [D_sigma, A, M, Intercept]\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:02<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 2 seconds.\n",
+      "The acceptance probability does not match the target. It is 0.8802196963185166, but should be close to 0.8. Try to increase the number of tuning steps.\n"
+     ]
+    }
+   ],
+   "source": [
+    "priors = {\n",
+    "    \"Intercept\": bmb.Prior(\"Normal\", mu=0, sd=0.2),\n",
+    "    \"A\": bmb.Prior(\"Normal\", mu=0, sd=0.5),\n",
+    "    \"M\": bmb.Prior(\"Normal\", mu=0, sd=0.5),\n",
+    "    \"sigma\": bmb.Prior(\"Exponential\", lam=1),\n",
+    "}\n",
+    "\n",
+    "model_5_1 = bmb.Model(\"D ~ A\", d, priors=priors)\n",
+    "results_5_1 = model_5_1.fit()\n",
+    "\n",
+    "model_5_2 = bmb.Model(\"D ~ M\", d, priors=priors)\n",
+    "results_5_2 = model_5_2.fit()\n",
+    "\n",
+    "model_5_3 = bmb.Model(\"D ~ M + A\", d, priors=priors)\n",
+    "results_5_3 = model_5_3.fit()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.33"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:146: UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>rank</th>\n",
+       "      <th>loo</th>\n",
+       "      <th>p_loo</th>\n",
+       "      <th>d_loo</th>\n",
+       "      <th>weight</th>\n",
+       "      <th>se</th>\n",
+       "      <th>dse</th>\n",
+       "      <th>warning</th>\n",
+       "      <th>loo_scale</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>model_5_1</th>\n",
+       "      <td>0</td>\n",
+       "      <td>126.958182</td>\n",
+       "      <td>3.721284</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.894053</td>\n",
+       "      <td>12.741410</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>False</td>\n",
+       "      <td>deviance</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>model_5_3</th>\n",
+       "      <td>1</td>\n",
+       "      <td>128.774804</td>\n",
+       "      <td>4.838442</td>\n",
+       "      <td>1.816623</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>12.855868</td>\n",
+       "      <td>0.735610</td>\n",
+       "      <td>False</td>\n",
+       "      <td>deviance</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>model_5_2</th>\n",
+       "      <td>2</td>\n",
+       "      <td>140.313353</td>\n",
+       "      <td>2.952110</td>\n",
+       "      <td>13.355171</td>\n",
+       "      <td>0.105947</td>\n",
+       "      <td>9.795754</td>\n",
+       "      <td>9.225215</td>\n",
+       "      <td>False</td>\n",
+       "      <td>deviance</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "           rank         loo     p_loo      d_loo    weight         se  \\\n",
+       "model_5_1     0  126.958182  3.721284   0.000000  0.894053  12.741410   \n",
+       "model_5_3     1  128.774804  4.838442   1.816623  0.000000  12.855868   \n",
+       "model_5_2     2  140.313353  2.952110  13.355171  0.105947   9.795754   \n",
+       "\n",
+       "                dse  warning loo_scale  \n",
+       "model_5_1  0.000000    False  deviance  \n",
+       "model_5_3  0.735610    False  deviance  \n",
+       "model_5_2  9.225215    False  deviance  "
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "az.compare(\n",
+    "    {\"model_5_1\": results_5_1, \"model_5_2\": results_5_2, \"model_5_3\": results_5_3},\n",
+    "    scale=\"deviance\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.34"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tomas/anaconda3/envs/bmb/lib/python3.8/site-packages/arviz/stats/stats.py:1406: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. \n",
+      "See http://arxiv.org/abs/1507.04544 for details\n",
+      "  warnings.warn(\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAA9hAAAPYQGoP6dpAABI6ElEQVR4nO3df3zN9f//8fvZWJofM2xoiVJDfm/yc6UQ70qIlCKhUqSf+qGU3nq/fane1VtSfRSJ9E4qraQitKXCmN8/kjS0ZMOM2ZvN2ev7x95bZj/O62yv1/m12/Vy6SLnvM7rPM5jr2P38zzP1/PlMAzDEAAAAADLBXm7AAAAACBQEbYBAAAAmxC2AQAAAJsQtgEAAACbELYBAAAAmxC2AQAAAJsQtgEAAACbELYBAAAAm1TxdgG+ICMjw9slWC4sLEyZmZneLsPn0SfX6JFr9Mg1emQOfXKNHrlGj8yxok/h4eEut2FkO0AFBfGjNYM+uUaPXKNHrtEjc+iTa/TINXpkjqf6xE8DAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsEkVbxcAAAAAlIfTaWjzFunIUaluHaltGyk42OHtsoogbAMAAMDvJCQa+vcMQ+npf90WESE9/IDU/SrfCdxMIwEAAIBfSUg0NHFS0aAtSenp0sRJhhISDe8UVgLCNgAAAPyG05k/ol2W6a8bcjp9I3ATtgEAAOA3Nm9RsRHtc6Wl5W/nCwjbAAAA8BtHjlq7nd0I2wAAAPAbdetYu53dCNsAAADwG23b5K86UpbIyPztfAFhGwAAAH4jONihhx8oe2m/h8Y5fGa9bcI2AAAA/Er3qxya8ryj2Ah3ZKQ05XmHT62zzUVtAAAA4He6X+VQXDdxBUkAAADADsHBDsW093YVZWMaCQAAAGATwjYAAABgE5+dRhIfH68NGzZo27Zt2r17t3JzczV16lQNHDiwyHa5ublauXKlVq5cqS1btujPP/+UJF166aW66aabdOuttyo4ONgbLwEAAACVnM+G7enTpys1NVXh4eGKjIxUampqidvt379fDz74oEJDQ9WlSxf16NFDJ06c0KpVqzR58mQlJibqzTfflMPhW5PlAQAAEPh8Nmz/85//VOPGjRUVFaVZs2bp5ZdfLnG7GjVqaNKkSbrpppsUGhpaeHt2drbuuOMOrVq1Sl9//bWuu+46T5UOAAAASPLhOdtdu3ZVVFSUy+3q16+voUOHFgnakhQaGqqRI0dKkpKSkmypEQAAACiLz4ZtK1Spkj9wz5xtAAAAeENAh+1PPvlEkhQXF+flSgAAAFAZ+eyc7YpauHChEhMT1blzZ3Xv3r3MbcPCwhQUFHifO8LDw71dgl+gT67RI9fokWv0yBz65Bo9co0emeOJPgVk2F61apX+8Y9/KCoqSi+99JLL7TMzMz1QlWeFh4crIyPD22X4PPrkGj1yjR65Ro/MoU+u0SPX6JE5VvTJTFgPuOHchIQEPfjgg6pbt67ee+89RUZGerskAAAAVFIBFba/++47jRs3TuHh4Zo3b54aNWrk7ZIAAABQiQVM2P7uu+/0wAMPKCwsTPPmzVPjxo29XRIAAAAquYAI2wkJCUWCdpMmTbxdEgAAAOC7J0guWrRIGzZskCTt3r278LZ169ZJkmJjYzV48GD9+uuvGjdunHJyctSxY0d9+eWXxfYVFRWlgQMHeq54AAAAQD4ctjds2KDFixcXuS05OVnJycmFfx88eLAOHz6snJwcSSoxaEtSx44dCdsAAADwOJ8N29OmTdO0adNcbtepUyf9/PPPHqgIAAAAcE9AzNkGAAAAfBFhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsAlhGwAAALAJYRsAAACwCWEbAAAAsEkVbxdQmvj4eG3YsEHbtm3T7t27lZubq6lTp2rgwIElbp+VlaUZM2Zo2bJlSk9PV2RkpPr06aNx48apevXqHq4eAAAA8OGwPX36dKWmpio8PFyRkZFKTU0tddvs7GwNGzZMO3fuVFxcnG644Qbt3LlTc+bMUVJSkhYsWKDzzjvPg9UDAAAAPjyN5J///KdWrlypNWvWaMiQIWVu+84772jnzp265557NHv2bD322GOaPXu27rnnHm3dulVz5871TNEAAADAWXw2bHft2lVRUVEutzMMQ4sWLVJoaKjGjh1b5L6xY8cqNDRUixYtsqtMAAAAoFQ+G7bNSklJUVpammJiYhQaGlrkvtDQUMXExOjAgQM6ePCglyoEAABAZeX3YXvfvn2SpCZNmpR4f8HtKSkpnikIAAAA+B+fPUHSrBMnTkiSatSoUeL9BbdnZWWVuo+wsDAFBfn9545iwsPDvV2CX6BPrtEj1+iRa/TIHPrkGj1yjR6Z44k++X3YtkJmZqa3S7BceHi4MjIyvF2Gz6NPrtEj1+iRa/TIHPrkGj1yjR6ZY0WfzIR1vx/OrVmzpqTSR64Lbi9t5BsAAACwi9+H7caNG0sqfU52we2lzekGAAAA7OL3YbtJkyaKjIxUcnKysrOzi9yXnZ2t5ORkXXjhhWrYsKGXKgQAAEBl5fdh2+FwaPDgwcrOztYbb7xR5L433nhD2dnZuuWWW7xUHQAAACoznz1BctGiRdqwYYMkaffu3YW3rVu3TpIUGxurwYMHS5LuvvturVixQm+//bZ27typyy+/XDt27NDq1avVunVr3Xnnnd55EQAAAKjUfDZsb9iwQYsXLy5yW3JyspKTkwv/XhC2Q0ND9f7772vGjBlatmyZ1q5dq4iICI0aNUr333+/qlWr5tHaAQAAAElyGIZheLsIbwvE5XFY9scc+uQaPXKNHrlGj8yhT67RI9fokTks/QcAAAD4OcI2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYJMq3i7AKoZhaPny5Zo/f75+++03nThxQg0aNFCnTp10zz33qFGjRt4uEQAAAJVMwITtF154Qe+++64iIiLUs2dP1ahRQ7t27dJHH32kJUuW6MMPP1R0dLS3ywQAAEAlEhBhOz09Xe+9956ioqIUHx+vmjVrFt43d+5cTZ06Ve+++66mTp3qxSoBAABQ2QTEnO3U1FTl5eWpffv2RYK2JF199dWSpIyMDC9UBgAAgMosIMJ248aNVbVqVW3cuFFZWVlF7vvuu+8kSZ07d/ZCZQAAAKjMHIZhGN4uwgpz587VtGnTVK9evSJztteuXavBgwdr4sSJqlKl5FkzeXl5CgoKiM8dAAAA8CEBE7YlaenSpXrmmWd08uTJwttiY2M1fvx4xcbGlvq4QJxiEh4eHpCvy2r0yTV65Bo9co0emUOfXKNHrtEjc6zoU3h4uMttAmY49/XXX9cTTzyh++67TwkJCUpOTtaCBQt0+vRpDR8+XCtWrPB2iQAAAKhkAiJs//jjj5oxY4aGDh2q0aNHq0GDBqpevbo6dOigt956S1WqVNELL7zg7TIBAABQyQRE2E5MTJQkderUqdh9ERERuuSSS7Rv374i00sAAAAAuwVE2M7NzZUkHT16tMT7jx49qqCgIFWtWtWTZQEAAKCSC4iwHRMTIyl/RZITJ04Uue8///mP/vzzT7Vr104hISHeKA8AgErH6TSUvNHQ8hX5fzqdAbMeA+CWgLiC5N/+9jf95z//UVJSkvr06aMePXqoZs2a2rFjh9asWaNq1arpqaee8naZAABUCgmJhv49w1B6+l+3RURIDz8gdb/K4b3CAC8IiJHt4OBgzZkzR+PHj1f9+vW1ZMkSzZs3T7/99pv69eunTz/9VG3atPF2mQAABLyEREMTJxUN2pKUni5NnGQoIZERblQuATGyLUkhISEaPXq0Ro8e7e1SAAColJzO/BHtskx/3VBcNyk4mBFuVA4BMbINAAC8b/MWFRvRPldaWv52QGVB2AYAAJY4UvKiYOXeDggEbk0jOXbsmLKzs1WzZk3VrFmzzG1PnDihEydOKDQ0VLVr165IjQAAwA/UrWPtdkAgMD2ynZ2drb59+6pfv35Kd/UdkaT09HT169dPAwYM0KlTpypUJAAA8H1t2+SvOlKWyMj87YDKwnTY/uyzz3T48GGNGTNGl1xyicvtL7nkEo0bN05//vmnPv/88woVCQAAfF9wsEMPP1D2iY8PjXNwciQqFdNhe8WKFTr//PM1bNgw0zu/7bbbdP7552vZsmXlKg4AAPiX7lc5NOV5R7ER7shIacrzDtbZRqVjes72rl271LZtW5133nmmd37eeeepbdu22rVrV7mKAwAA/qf7VQ7FdctfdeTI0fw52m3bsNwfKifTYTszM1P16tVz+wnq1aun9evXu/04AADgv4KDHYpp7+0qAO8zPY2kevXqysrKcvsJsrKyFBoa6vbjAAAAAH9nOmxHRUVp27Ztbj/Btm3bFBUV5fbjAAAAAH9nOmx36dJFR44ccWtlkfj4eB0+fFjdunUrV3EAAACAPzMdtocNG6aqVatq8uTJSkpKcrn9unXrNHnyZJ133nlurWACAAAABArTYbthw4aaMGGCTp48qTvvvFPjx4/XqlWrdOjQIeXm5io3N1eHDh3SqlWrNH78eI0YMULZ2dmaMGGCGjRoYOdrAAAAAHySW5drv/322xUUFKSpU6fqyy+/1NKlS0vczjAMhYSE6Omnn9aQIUMsKRQAAADwN26FbUkaMmSIrr76as2fP1+JiYnas2ePDMOQJDkcDl166aXq3r27hg0bxog2AAAAKjW3w7YkNWjQQI8//rgef/xx5ebm6vjx45KkWrVqqWrVqpYWCAAAAPircoXts1WtWlV169a1ohYAAAAgoJg+QRIAAACAe0yPbA8fPrzcT+JwOPTee++V+/EAAACAPzIdttetW1fuJ3E4HOV+LAAAAOCvTIftefPm2VkHAAAAEHBMh+2OHTvaWQcAAAAQcGw/QTIjI0MLFiyw+2kAAAAAn1Phpf9KkpOToxUrVujzzz/X999/L6fTqaFDh9rxVAAAAIDPsjRsJyUlKT4+Xt98842ysrJkGIaqVKmi7t27W/k0AAAAgF+ocNj+9ddf9fnnn+uLL77QwYMHJUmGYahjx47q27ev+vTpo7CwsAoXCgAAAPibcoXtI0eOaMmSJfr888+1Y8cOSfkB+7LLLlNGRoaOHDnC6iUAAACo9EyH7dOnT+vbb79VfHy8fvzxRzmdThmGoYiICPXt21f9+/dX8+bNdfvtt+vIkSN21gwAAAD4BdNhu2vXrsrOzpZhGAoNDVXv3r3Vr18/denShYvWAAAAACUwHbZPnjwph8OhyMhITZ48Wddcc42ddQEAAAB+z/Q621deeaWCgoKUlpamsWPH6uqrr9ZLL72kXbt22VkfAAAA4LdMj2y//fbbOnLkiL744ovCEyNnz56tOXPm6NJLL9WAAQPUt29fO2sFAAAA/IrDMAyjPA/89ddfFR8fryVLluiPP/6Qw+GQw+FQlSpVlJubq6SkJNWoUcPqem2RkZHh7RIsFx4eHpCvy2r0yTV65Bo9co0emUOfXKNHrtEjc6zoU3h4uMttyn259qZNm+rRRx/VypUrNW/ePA0cOFDVq1dXTk6ODMNQ165ddf/992vp0qU6depUeZ8GAAAA8FvlHtkuScFl2uPj47V69WqdOXNGDodD559/vpKTk616GssF4qc/PtWaQ59co0eu0SPX6JE59Mk1euQaPTLHUyPbpuds33jjjerQoYNiY2PVoUMHNWjQoNg2ISEhuu6663TdddcpIyNDX375peLj47V161b3KgcAAAACgOmw/csvv2jPnj368MMPJUkNGzZUhw4dCgN406ZNi2wfHh6uYcOGadiwYUpJSbG0aAAAAMAfmA7b77//vtavX6/169dr06ZN+uOPP/T555/riy++kCTVrl27cNQ7NjZWLVu2VFBQ/pTwJk2a2FI8AAAA4MtMh+2CUWxJMgxDu3bt0vr167VhwwZt2LBB6enp+vbbb7VixQpJ0vnnn6927doVBvBOnTrZ8woAAAAAH2U6bJ/N4XCoRYsWatGihe644w5J0oEDB7R+/XolJSUpOTlZKSkp+umnn/TTTz/J4XBox44dlhYOAAAA+Lpyhe2SNGrUSI0aNdKAAQO0a9cuLVu2TB988IEyMzOtegoAAADAr1Q4bOfk5Gjz5s2F00k2bdqkrKwsGYah4OBgXX755YqJibGiVgAAAMCvuB22MzMzlZycXDhfe8eOHcrNzZVhGKpRo4batm2rmJgYxcTEqG3btgoNDbWjbgAAAMDnmQ7bkyZNUnJysvbu3au8vDxJ0gUXXKDevXsXhutmzZrJ4XDYViwAAADgT0yH7Y8++kgOh0PNmzfXnXfeqS5duqh+/fp21gYAAAD4NdNhu0aNGsrKytLOnTv17LPPqmXLloqNjVVsbKzat29v6nKVAAAAQGViOmwnJSXp559/LjwRsuDiNnPmzJGUf+GagvAdGxuriy66yLaiAQAAAH9gOmwXTCFp3ry5hg4dKil/be2C4J2cnKyPP/5YH3/8sRwOh+rWrauYmJjC8N2qVSvbXgQAAADgiyq09N/Za2tL0tGjR5WcnFwYwFeuXKnly5dzURsAAABUSkFW7qxGjRqqVauWatWqpbCwMIWEhMgwDBmGYeXTAAAAAH6hQiPbx44dK1xzOzk5Wdu3b9eZM2ckqTBgN27cWLGxsRWvFAAAAPAzboXt33//vfAEyQ0bNui3334rMnIdFBSk5s2bKzY2Vh06dFCHDh1Ur149WwoHAAAAfJ3psN29e3elpaVJ+mvUOiQkRK1bt1aHDh0UGxurmJgY1ahRw55KAQAAAD9jOmwfOnRINWrUUPv27QvDdZs2bRQSEmJnfQAAAIDfMh22Fy9erObNm3M5dgAAAMAk02G7RYsWdtZhmeXLl+uDDz7Qjh07lJ2drYiICLVr106PP/64GjZs6O3yAAAAUIlUaDUSX2IYhp577jktXLhQF110ka6//npVr15daWlpSkpKUmpqKmEbAAAAHhUwYXvevHlauHChbr/9dj3zzDMKDg4ucn/BkoQAAACAp1h6URtvOXXqlGbOnKlGjRpp4sSJxYK2JFWpEjCfKwAAAOAnAiKBrl69WpmZmRo4cKDy8vK0bNkypaSkqGbNmuratasaN27s7RIBAABQCQVE2N6+fbuk/Ivq3HjjjUpJSSm8LygoSCNGjNCTTz7ppeoAAABQWQVE2D5y5Igkae7cubr88su1aNEiNW3aVDt37tSzzz6rOXPmqFGjRrr99ttLfHxYWJiCggJiRk0R4eHh3i7BL9An1+iRa/TINXpkDn1yjR65Ro/M8USfHEbB5SD92LPPPquPPvpI1apV07Jly1S/fv3C+3bv3q3+/fvrwgsv1PLly0t8fEZGhqdK9Zjw8PCAfF1Wo0+u0SPX6JFr9Mgc+uQaPXKNHpljRZ/MhPWAGM4tuER8q1atigRtSYqOjlajRo20f/9+HT9+3BvlAQAAoJIKiLB9ySWXSJJq1qxZ4v0Ft586dcpjNQEAAAABEbY7deokSdq7d2+x+3Jzc7V//36FhoaqTp06ni4NAAAAlVhAhO2LLrpIcXFx2rdvnxYtWlTkvlmzZun48ePq1asXa20DAADAowImfT733HMaMmSInnnmGX377be65JJLtGPHDq1Zs0ZRUVF64oknvF0iAAAAKpmAGNmW8ke3P/nkEw0cOFDbt2/X/PnztW/fPg0dOlSLFi1SRESEt0sEAABAJRMwI9uS1LBhQ02dOtXbZQAAAACSAmhkGwAAAPA1hG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCZVvF0AAOs5nYY2b5GOHJXq1pHatpGCgx3eLgsAgEqHsA0EmIREQ/+eYSg9/a/bIiKkhx+Qul9F4AYAwJOYRgIEkIREQxMnFQ3akpSeLk2cZCgh0fBOYQAAVFKEbSBAOJ35I9plmf66IaeTwA0AgKcQtoEAsXmLio1onystLX87AADgGYRtIEAcOWrtdgAAoOII20CAqFvH2u0AAEDFEbaBANG2Tf6qI2WJjMzfDgAAeAZhGwgQwcEOPfxA2Uv7PTTOwXrbAAB4EGEbCCDdr3JoyvOOYiPckZHSlOcdrLMNAICHcVEbIMB0v8qhuG7iCpIAAPgAwjYQgIKDHYpp7+0qAAAA00gAAAAAmxC2AQAAAJswjQQAUCqn0/Da/H9vPjcAWIWwDQAoUUKioX/PMJSe/tdtERHSww/I9pVtvPncAGClgJ5GMmvWLDVr1kzNmjXTpk2bvF0OAPiNhERDEycVDbuSlJ4uTZxkKCHRCMjnBgCrBWzY3r17t2bMmKHQ0FBvlwIAfsXpzB9VLsv01w05ndaHXm8+NwDYISDDdm5uriZMmKAWLVqoV69e3i4HwFmcTkPJGw0tX5H/J6HJ92zeomKjyudKS8vfzh1mfvZ2PTcAeEtAztl+66239Msvv2jx4sV65513vF0OgP9hHq5/OHLU2u0kafm3pzVlmuufvR3PDd/ECbCoLAIubG/fvl1vvfWWHnzwQV166aXeLgfA/xTMwz1XwTzcKc8TuH1F3TrWbpeQaOiZ57JknPPjL+lnb/VzwzfxwRuVSUBNI8nJydGTTz6p5s2b6+677/Z2OQD+h3m4/qVtm/zgU5bIyPztXCn42Z8btM929s/eyueGb+IEWFQ2ATWyPX36dKWkpOjTTz9VcHCw6ceFhYUpKCigPndIksLDw71dgl+gT65VtEfrknKVnn68zG3S0qS9v9VUxyuqVui5vCXQjqOJE07rkceyJKlIUHb8b9Dx6SdrqF6981zupzw/e6ue218F2rF0NqfT0Gszj0kqOVA7HNKMmQ7d2Ld2mVNKArlHVqFH5niiTwETtjdu3Kg5c+Zo3Lhxio6OduuxmZmZNlXlPeHh4crIyPB2GT6PPrlmRY9S9pkbqUrZd0KXXep/XyEH4nHUIVb652RHiV/1PzTOoQ6x2crIyHa5n/L87K16bn8UiMfS2ZI3Gjp0qPRjwjCkPw/l6buEDMW0L/nfgkDvkRXokTlW9MlMWA+IsH3mzBlNmDBBzZo10+jRo71dDoBzMA/XP3W/yqG4bqrQSWzl/dlb8dzwPZwAi8ooIMJ2dna2UlJSJEmtWrUqcZtbb71VkjRz5kyWAwQ8rGAebllLujEP1zcFBzsU0778jy/42R8+rFLnbZf2s6/oc8P38MEblVFAhO2QkBDdfPPNJd63fv16paSkqEePHqpTp46ioqI8XB2A4GCHHn5AJa5GUuChcQ5GLQNQwc/+mef42YMP3qicAiJsV6tWTVOmTCnxvgkTJiglJUX33nuv2rVr59nCABTqfpVDU55XsXm4kZH5YYvlvgJX96scevVf1TVlWhY/+0qOD96ojAIibAPwD8zDrbyu7XWe2rc7yc8efPBGpUPYBuBRzMOtvPjZowAfvFGZBHzYnjZtmqZNm+btMgAAwFn48IXKIvCu5AIAAAD4CMI2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgE8I2AAAAYJMq3i4AAHyZ02lo8xbpyFGpbh2pbRspONjh7bJQSXD8Af6PsA0ApUhINPTvGYbS0/+6LSJCevgBqftVBB7Yi+MPCAxMIwGAEiQkGpo4qWjQkaT0dGniJEMJiYZ3CkOlwPEHBA7CNgCcw+nMH1Esy/TXDTmdBB5Yj+MPCCyEbQA4x+YtKjaieK60tPztAKtx/AGBhbANAOc4ctTa7QB3cPwBgYWwDQDnqFvH3Ha//87X+LCe2ePP7HYAvIuwDQDnaNsmf9UHV2a/q8IT1ZxOQ8kbDS1fkf8n82lRXmaOv8jI/O0A+D6W/gMqIdbuLVtwsEMPP5C/6oMr0183dP75pzT1RZZogzXMHH8PjXPwngX8hMMwjEo//JKRkeHtEiwXHh4ekK/Lar7aJzvDcOlr9zpKDIa+2iNPePe9PM1+t2L7mPJ8yX2tbCrzceSOs/tU0ns1MjI/aFfmY4pjyTV6ZI4VfQoPD3e5DSPbgI+x80IWBWv3nqtg7d4pzzMSe7YLL3RIqth4xPTXDcV145sDuK/7VQ7FdRPfQgF+jjnbgA+x80IWrN3rPitOQGOJNlREcLBDMe0durZn/p8EbcD/ELYBH2F3GGbtXveZPVHSFZZoA4DKi7AN+Ai7wzBr97ov/0S1io8kskQbAFRehG3AR9gdhlm7t3y6X+XQlOcdxUa4IyOlf/w9f+TbUUYeZ4k2AKjcOEES8BF2h+GCKRFljZ4TDEtW1olqQUGGnnmOJdoAACVjZBvwEXZfyMLMlAiCYelKO1Gt+1UOvfqvGiWOfLPsHwCAkW3AR3jiQhb5UyLE2r0Wu7bXeWrf7iRLtAEAiiFsAz7EE2HYn9fu9eUrX+aPfHu7CgCAryFsAz7GE2HYH4OhnRf7AQDALoRtwAf5Yxi2ky9d+bKk0XUAAEpD2Abg08xe7McTl0QvbXR94oTT6hBr61MDAPwUq5EA8Gm+cuXLgtH1c2tJT5ceeSxLCYlc5h4AUBxhG4BP84UrX5odXXc6CdwAgKII2wB8mi9c+dLV6LpheGZ0HQDgfwjbAHya3Rf7McMXRtcBAP6JsA3Ap/nClS99YXQdAOCfCNsAbOF0GkreaGj5ivw/KzKfOf9iPw6vXRLd1ei6w2H/6DoAwD+x9B8Ay9lxARpvXvkyf3RdJa71XcDu0XUAgH9iZBuApcpaIm/iJKNCS+TlX+zHoWt75v/pyXBb1uj6q/+qwVUsAQAlYmQbgGVycvL04stlb+OpC9DYobTR9Xr1zlNGRra3ywMA+CDCNgBLJCQaeullKTOz7O0Klsjz18vR54+ue7sKAIC/IGwDqLCCqSNmsUSee5xOwytz1QEAFUfYBuC2s8NfeG1D/57h3uNZIs88O042BQB4DmEb8DG+PopZUvhzB0vkmVfaNwYFJ5tOeZ7ADQC+jrAN+BBfH8V0d7pISR4Y61sfHnyV05l/LJTFn082BYDKgqX/AB9h55J5VjAT/sx4baa8/lr8weYtcvntQcHJpgAA30XYBnyA2VHMilyFsaLMhD8zfOXDg68zexIpJ5sCgG8LiGkkhw4d0ldffaXExETt3btXhw8fVlhYmGJiYnT33Xerbdu23i4RKJPZUcw5cw3FxnhnHrfVoY4pEGUzexIpJ5sCgG8LiJHt+fPna+rUqTpw4IC6deumkSNHKjY2VitWrNCQIUO0dOlSb5cIlMlskH1vvvTgI4ZuHuL5kWGrQx1TIMrWto2KXa3yXJxsCgC+LyBGttu0aaP58+erY8eORW5fv369RowYob///e/q1auXQkJCvFQhUDZ3g6w3VqMoCH9ljcBHREgTJ0iJq6VPF7veZ1kfMnx9VRa7BQc79PADKvOE1IfGefaS9YAZlf29C5wrIMJ27969S7y9Q4cO6tSpk1avXq2ff/5ZrVu39nBlgDlmgmxJzEzFsOoXn5nw9/ADDrVvJ/2619yoe3htQ8kbi9fm66uyeEr3qxya8ryK9SIyMj9oV6ZewD/w3gWKC4iwXZYqVaoU+RPwRWaCbEnS0qSNmwwFBf11gRlJyjjmUN06UmamoddmyrJffK7CnyTdPMTcGty1aklTpknp6X+95ogI6dqehj74sPj2lXVt6e5XORTXTYwUwuexLjxQModhGAG7JMAff/yhPn36KCwsTAkJCQoODi5xu4yMDA9XZr/w8PCAfF1W87U+leeCMbVqSsdPuP9cU543NzJaUo9KGi1f/YP7HxbKIzJSWvQf35o+4WvHkS+iR+b4a5+cTsPlB22r3rv+2iNPokfmWNGn8PBwl9sEbNjOzc3VyJEjlZSUpBdeeEEDBgwoddu8vDwFBQXEuaJ+w+k0tCH5jNIP5ymiXpBiY6r4VHjypoLerFmbq/97+7+2PIfDIdWPDNKyr2pb0nen09C11x3ToUN5LretH+nQ6dPSsczy/9Pz7ju11PGKquV+PABrrUvK1ci7j7vcjvcuKqOAnFuRl5enCRMmKCkpSbfcckuZQVuSMjMzPVOYB/nyp9rS5/R5fg6qr/bpskulSy429Oln1qxtfS7DkP48lKfvEjIU077snpvpUfJGQ4cOuQ7PD9wvNb3E0MPj3Sq3mJR9J3TZpb7z4cxXjyNfQo/M8dc+pewz9+HZiveuv/bIk+iROZ4a2Q644dy8vDw9/fTTWrJkifr166fJkyd7uyScxdevkmglp9NQ8kZDy1fk/+nuBWny53HbGyitWjvb7H7q1HEo41jFXxNrS3tXRY9tBB7WhQdKF1Aj23l5eXrqqaf02WefqW/fvpo2bRrTQ3yI2askBsKFTqw6I7+0ExJr1ZKOu/7G1iWrfvF58hcta0t7F6tNoCRmVlTivYvKKmCS6NlB+/rrr9eLL75Y6gmR8A6zV0n09wudWD163/0qhz7+0KHXXnXouWfz/3z+uYrXaeUvPncuwGJm27KwtrT3VKZvprzB6TS0LinXL78xMPNNHO9dVFYBEbYLpo589tln+tvf/qaXXnqJoO2DzE41sPqy4J5kdvS+PFNKYto7dG3P/D/bt3NUKLBK1v7ic+cXrZltbx9SPJBHRppfQQXWs+vYRr6ExPzVPEbefVyT/2F47UqxFZH/TVzxf5t476KyC4hpJDNnztTixYsVGhqqJk2a6M033yy2Ta9evdSiRQsvVIcCZqcQHD1qyOn0z6kk7ozex7Qv//OUd11uyb4LorhzARYz2957D1eh8yWeOrYro0Ban5p14YHiAiJsp6amSpKys7P11ltvlbhNVFQUYdvLzF4lccZM6cOPDL+cA+rJ0fvSAmtZ7hopDR/214i21ZdVducXratt80fzy10KLFYZvpnyhkA8l4X3LlBUQITtadOmadq0ad4uAy64MxrrjyM6kufPyD87sKYfNvTaDCmzjBMnv/hSGj4s///tOtHNnV+0/FL2H6w2YQ++MQACX0DM2YbvK1gqLCc3f3S1Xj1zj3vhX4bWb8izdR6olcuYuXOioFUK5nNH1HOUGbSlv35pc6Ib3OWNY7sy4BsDIPAFxMi2v7H6q3tfV9oI6nV/k776uuzHHj8uPTxeioiwZ1rJ8m9Pa8o060Z3zYzemz0x0d3jxOwv4/TDht6aVfY2/va1dXlUtvdhRVl5bOMvfGMABD7Ctof58xq15QknZZ344ypon7u91dNKEhINPfNcloxzyqvoc7lzomBZtbl7nJj9ZXzsWGB/bW3mOPXn96E3WXFsoyjWpwYCH2Hbg/z5jPPyhBMzJ/64q2DEVTJ/tntJ4UvKDwznBu2Snqs8I3UVOSO/vMeJ2V/atWubew2e/NrailFmp9PQvPcNffSxdOLEX7efe5z68/vQF7DahLX4xgAIfIRtD/HnM87LG07mvW9+lQyz0tLy9/v5EpkK/qV9SOjX1/7R3fKc/FeR48TsL+2aNSXJ9YcgT31tbcUoc0KioRf+ZZR4Vc2zj9O4bvLb96Ev4cRWa/GNARDYOEHSQ/z16onlvZBFQqKh2e/aU9Psd4v3sqQT+8o6CdBsbZ4+Kamix4mZi0r40oluVpyoWbAPV5evn/66oY2bXH8A9MX3IQJfwZVi332nVuGVYhf9h6ANBAJGtj3EX884L8+yVO5OH7lrpIqNVJfX2dNMrJjC4umTkqw4TsysX+0LX1tb8W2PO8daWpq0cZO52nztfYjKITjYoY5XVNVllxKwgUDCyLaH+OsZ5+UJf2YCeoHIyPyLrHz8oUP/flmqVdP9Gs9WEPzdqaGs2jx9UpJVx8m5l3c/N6z6wmWVrfi2x4qfc0l87X0IAPBfjGx7iL+ecV6e8OfOqODZI6gdYh168vGS54e7w6pRSW+clOTJ48TbJ7pZMYrv7s+6fTtp6df+9z4EAPgvRrY9JP+r+7JDjC+ecV6e+b1mA/pdI4ufAFfWiOtdI83tt24d8zXcP+Z8r47unsvTx4mrEXA7WTGK784IdGSk1L6df74PAQD+i5FtD/LHM87LM7/XzOhsRET+9JGSlDbiKkmfLyn7BLezg7+ZEeJ77zlft9x8yqeWMfPH46Q8rBjFN7OPAgXHaferVCn6CwDwDQ7DKGul4cohIyPDo8/niSvXhYeHW/q6SlqeraxwUtpygQXKO3Lszn7NbDugfx2P//zN8pUrHFp9LJ3NiuPE1T5q1ZKefKz4fqzsr509ChT0yBz65Bo9co0emWNFn8LDw11uQ9iW58O2J9jxRnM3nLgb0M1yZ7+utuUfJNfs7pEVx0lJ+6hVSxo8KP8bFLs/pHAcuUaPzKFPrtEj1+iROYRtDwrEA9JX3mh2jc66s9+ytvWVPvkyT/TIqitIeuubAI4j1+iROfTJNXrkGj0yx1NhmznbsJVdV5pzZ79c7c73WfEz4ucMAPBFrEYCAAAA2ISwDQAAANiEsA0AAADYhLANAAAA2ISwDQAAANiEsA0AAADYhLANAAAA2ISwDQAAANiEsA0AAADYhLANAAAA2ISwDQAAANiEsA0AAADYhLANAAAA2MRhGIbh7SIAAACAQMTINgAAAGATwjYAAABgE8I2AAAAYBPCNgAAAGATwjYAAABgkyreLgAVl5WVpRkzZmjZsmVKT09XZGSk+vTpo3Hjxql69eqm9pGSkqKvv/5a33//vfbt26djx46pbt266tSpk+699141bdrU5ldhLyt6JEnz58/Xjh07tG3bNv36669yOp2aN2+eOnXqZGP11tqyZYtmzJihjRs36syZM4qOjtaIESN0/fXXm95HTk6OZs2apc8//1wHDx5UWFiYrrnmGj388MOqW7eujdV7RkV7tH//fsXHx2v79u3avn270tLSFBUVpZUrV9pcuedUpEeGYSgxMVErV65UcnKy/vjjD505c0aNGzfW9ddfr5EjR+q8887zwKuwX0WPpYSEBH322WfauXOnDh8+rNzcXDVs2FAxMTG65557dPHFF9v8Cuxnxb9JZ8vMzFTfvn2VlpamuLg4zZ492+KKPa+iPfr000/11FNPlXq/v/0eK41Vx9KRI0f0f//3f/ruu+908OBBhYaGqkmTJurfv79uv/12t+sibPu57OxsDRs2TDt37lRcXJxuuOEG7dy5U3PmzFFSUpIWLFhg6pfW9OnTtXTpUkVHR6tnz56qUaOGdu/erfj4eH3zzTd65513dMUVV3jgFVnPqh5J0j//+U9JUkREhOrUqaP09HQ7S7fcmjVrdPfddyskJEQ33HCDqlevrmXLlumRRx7Rn3/+qVGjRrncR15ensaMGaPVq1erXbt26t27t/bt26dFixbpp59+0kcffaQ6dep44NXYw4oerV+/Xq+//rqCg4PVtGlTHT582AOVe05Fe5STk6PRo0crJCREHTt2VFxcnHJycrR69Wq9+uqr+vbbbzV//nydf/75HnpF9rDiWEpMTNTmzZvVpk0bRUZGqkqVKtq7d68+++wzffHFF5o1a5a6dOnigVdjDyt6dK7nn39eWVlZNlTrHVb2qGfPnmrRokWx26Oioqws2Sus6tPOnTs1atQoHT9+XN27d1efPn2UnZ2tX3/9VatWrSpX2JYBvzZ9+nQjOjraeOmll4rc/tJLLxnR0dHGW2+9ZWo/n3zyibF9+/Zity9ZssSIjo42rr/+ekvq9QaremQYhrFq1SojLS3NMAzDePbZZ43o6GhjzZo1ltZrl9zcXKNXr15Gq1atjB07dhTefvz4caN3795Gy5Ytjd9//93lfj7++GMjOjraePTRR428vLzC2z/44AMjOjraePbZZ22p3xOs6tH+/fuNjRs3Gv/9738NwzCMVq1aGddcc41tdXuSFT3Kyckx3njjDePYsWPFbr/33nuN6Oho4+2337alfk+x6lg6depUibf/+OOPRnR0tDFw4EDLavY0q3p0tq+//tqIjo423n//fSM6OtoYNWqU1WV7lFU9+uSTT4zo6Gjjk08+sbNcr7GqTydOnDCuvvpqo3PnzsbOnTtLfJ7yYM62HzMMQ4sWLVJoaKjGjh1b5L6xY8cqNDRUixYtMrWvgQMH6vLLLy92+w033KAmTZpoz549Onr0qCV1e5KVPZKkq6++WhEREVaX6RFr1qzR/v371bdv3yIjGzVr1tR9992n3NxcLV682OV+Cvr16KOPyuFwFN4+ZMgQNWrUSF988YVOnTpl/QvwAKt61KhRI7Vr107VqlWzs1yvsKJHVatW1ZgxYxQWFlbs9nvvvVeSlJSUZH3xHmTVsVTat25dunRRWFiY9u/fb1nNnmZVjwocPXpUf//739W/f391797djpI9zuoeBSqr+vTBBx/ojz/+0Pjx49W8efNi91epUr4JIYRtP5aSkqK0tDTFxMQoNDS0yH2hoaGKiYnRgQMHdPDgwQo9T9WqVSWV/yDzJk/1yB+sW7dOkhQXF1fsvoLbXAWc06dPa/Pmzbr44ouLfe3ocDjUtWtXZWdna9u2bRZV7VlW9CjQ2d2jgn9ngoODy70PX2B3nzZu3KjMzExddtll5d6Ht1ndo+eee07BwcGaOHGiNQX6AKt7tGPHDs2ZM0ezZs3S0qVLlZGRYU2hXmZVn5YuXSqHw6E+ffpo7969mj9/vt5++22tWLFCOTk55a7P/9ITCu3bt0+S1KRJkxLvb9KkiVavXq2UlBQ1bNiwXM+xZcsW/fLLL2rdurVq1apV3lK9xhM98hcpKSmSpMaNGxe7LyIiQqGhoYX9Ks3+/fuVl5dXZj8LnqtDhw4VKdcrrOhRoLO7R5988okkqVu3buXehy+wuk+rV6/Wxo0blZOTo3379mnVqlUKDw8v86Q3X2dlj+Lj47Vs2TLNnDlTYWFhOnHihJWleo3Vx9H8+fOL/L1atWq6//77NXr06ArV6W1W9CknJ0e7d+9WnTp1NH/+fM2YMUN5eXmF9zdq1EgzZ85Us2bN3K6PsO3HCv4xqVGjRon3F9xe3hNFTpw4oSeffFJBQUF6/PHHy1ekl9ndI39S8Bpr1qxZ4v01atRw+Qsq0PtpRY8CnZ09SkhI0MKFC9W0aVMNHjy43DX6Aqv79MMPP2jOnDmFf2/cuLFeeeUVtWrVqmKFepFVPTp06JCmTJmivn37qlevXpbW6G1W9ejCCy/Us88+q7i4ODVo0ECZmZn66aef9Morr+jll1/W+eefrzvuuMPS2j3Jij5lZmbK6XTq2LFjeuONN/T444+rf//+OnPmjD788EO9+eabGjNmjL766iu3V0sibPuAadOmufX1xPDhw0sdWbTKqVOndP/992vv3r165JFHvL4kkC/2CIB1tmzZokceeUQ1a9bU9OnTFRIS4u2SfMqTTz6pJ598UidPntSvv/6qmTNn6rbbbtP/+3//TzfeeKO3y/OqZ555RlWqVAmo6SNW69ixozp27Fj492rVqmnAgAFq2bKlBg0apNdff1233XabX04XtUrBKLbT6dTQoUOLrF7y0EMP6bffftNXX32lr7/+Wv3793dr35W3qz5k4cKFys7ONr19nz591KRJk8JPcKWNIhbcXtooZGlOnz6tsWPHau3atbr33nt13333ufV4O/haj/xRwWss7dN9VlZWsRPWzhXo/bSiR4HOjh5t3bpVd911l4KCgvTOO+/49TzkAnYdS9WrV1ebNm00c+ZMDRo0SJMmTVK3bt38crlNK3q0ePFiJSYmavr06X7ZA1fs/jfpsssuU2xsrH788Uf9+uuv5Zoi4Qus/P0mST169Ch2f48ePfTVV19p27ZthG1/tHHjxnI9rmBuUsFcpXMV3O7OCO+pU6c0duxY/fDDD7r77rv16KOPlqs2q/lSj/xVwWvct29fsa+e09PTlZ2drTZt2pS5j0aNGikoKChg+2lFjwKd1T3aunWrRo0apby8PM2ZMydg+mv3sVSlShV16tRJu3bt0tatW/1y9Q0rerRjxw5J+SOPJVm9erWaNWum5s2bKz4+vuJFe5gn/k0KDw+XJP33v/+t0H68yYo+hYaGqn79+jp06FCJ56gV3Hb69Gm362M1Ej/WpEkTRUZGKjk5udiob3Z2tpKTk3XhhReaPvHv7KA9atQov52nfTare+TPCi5KtHr16mL3Fdzm6sJF1apVU5s2bfTbb78pNTW1yH2GYejHH39UaGio384jtaJHgc7KHhUEbafTqXfeeUdt27a1rlAv88SxlJaWJumvFaP8jRU9at++vW6++eZi/xVcMbBBgwa6+eabde2111pcvWfYfRw5nc7C1aMuuOCCcu/H26zqU+fOnSVJe/bsKXZfwW3lugBQuVbnhs9w94It2dnZxp49e4zU1NQit586dcoYOXKkER0dbUydOtX2uj3Jqh6dyx8vatOzZ88yF/0/cOBA4e2HDh0y9uzZYxw/frzIfgL9ojZW9OhcgXZRGyt6tHXrVqNDhw5Gu3btjPXr13usfk+xqk9btmwpcf+JiYlGy5YtjQ4dOhgnT56050XYzK73m2EYxoEDBwLmojZWvd/OdebMGWPatGlGdHS0cccdd9j3IjzAqj5t2LDBiI6ONm644QYjMzOz8Pa0tDTjyiuvNJo3b27s3bvX7fochmEY7kd0+Irs7Gzddttt2rVrl+Li4nT55Zdrx44dWr16tVq3bq3333+/yIU11q5dq+HDh6tjx45FlgCaMGGCFi9erIiICN16660lPtdNN92kCy+80PbXZDWreiRJs2bN0t69eyXlT21JSUlRXFxc4YVuevXq5dNnw5d2OdvU1FQ9+eSTRU4IKTgmpk6dqoEDBxbenpeXp3vuuafwcu1XXHGF9u/fr2XLlikqKkqLFi3y67mTVvTo6NGjevHFFwv/Hh8fr2rVqqlPnz6Ftz3xxBN+26eK9ujYsWPq3bu3MjMzdeWVV5Y4ol2zZk2NGDHCUy/JFlYcS82aNVN0dLSio6PVoEED/fe//9XPP/+s9evXq2rVqnr11Vf9dtRWsqZHJfn999/Vs2dPxcXFafbs2Xa/DFtZdRwV/Fe/fn1lZmZq3bp1SklJUYMGDfT++++rUaNG3nh5lrHqWJo2bZreffddNWzYUNdcc43OnDmjFStW6MiRI3r00UcLL7zlDuZs+7nQ0FC9//77mjFjhpYtW6a1a9cqIiJCo0aN0v3332/6CnYFUwLS09P1+uuvl7hNx44d/TJsW9UjSfr+++8LF88vcPbXVlFRUT4dtjt37qwPPvhAr732mpYuXaozZ84oOjpajz32WOHXrq4EBQXpzTff1KxZsxQfH6+5c+eqdu3auvnmm/Xwww/7bYAsYEWPsrOzi12t7Nzbxo0b57e9qmiPsrKylJmZKSn/PfX9998X2yYqKsrvw7YVx9Kjjz6qtWvXKikpSUePHlVQUJAaNmyoW2+9VXfeeaeaNm1q86uwlxU9CnRW9GjUqFHatGmTfvzxR2VmZqpq1aq66KKLNGbMGI0cOTIgTvy26liaMGGCoqOjtWDBAi1evFgOh0MtWrTQ5MmTy/3BlpFtAAAAwCacIAkAAADYhLANAAAA2ISwDQAAANiEsA0AAADYhLANAAAA2ISwDQAAANiEsA0AAADYhLANAAAA2ISwDQAAANiEy7UDQAmaNWtW5O8Oh0M1atRQdHS0brrpJt18881yOBxFtsnKytLcuXO1YsUKpaSkKDc3V3Xq1FHDhg0VExOjHj166Iorrijcfu3atRo+fLg6duyo+fPnF9lXXl6ePv30U8XHx+vnn3/WyZMnFRYWpnr16qlNmzbq2LGj+vXrZ+q1fPrpp3rqqaeK3FalShXVrVtXHTp00F133aWWLVu6056AM2HCBC1evFjz5s1Tp06dvF0OgABC2AaAMtx0002SJKfTqQMHDig5OVkbNmzQTz/9pFdeeaVwuz/++EPDhg1TamqqQkND1aZNG9WrV0/Hjh3Ttm3btGnTJu3evbtI2C5NTk6Oxo4dq++//15BQUFq06aNLrjgAuXk5GjXrl1atGiR4uPjTYftAhdddJFiY2MlSdnZ2dq2bZu+/PJLLVu2TNOnT1fPnj3d2p+dfv/9d/Xs2bPEDyIA4E8I2wBQhmnTphX5+w8//KDRo0fryy+/1I033qhrrrlGkvT8888rNTVVcXFxevnll1W7du3Cx+Tl5Wnt2rX6+eefTT3nggUL9P3336thw4aaPXu2mjZtWuT+X375RYsXL3b7tcTGxhZ5Pbm5uZo8ebIWLVqk5557TldeeaVCQkLc3i8AoHTM2QYAN3Tr1q1wRPnbb7+VJJ06dUrff/+9JGnSpElFgrYkBQUFqUuXLhoxYoSp5/jmm28kSffff3+xoC1Jl112mZ544olyvoK/VK1aVU8//bRCQ0OVnp6uLVu2VHifAICiCNsA4KbLL79ckvTnn39KkjIzM3XmzBlJUp06dSq8/4yMDElSeHh4hfflSmhoqC6++GJJf72enTt36sUXX9TAgQPVuXNntWrVSj179tTf//53HTp0qNg+fv/9dzVr1kx33HGHsrKyNHXqVPXo0UMtW7bUlClTCrc7duyYXn75ZV1//fVq06aNYmNjNXz4cK1atarI/mbMmFE4pWXdunVq1qxZ4X8TJkwosu2ePXs0fvx4xcXFqVWrVrryyiv1xBNPaO/evZb0JycnRw8++KCaNWumcePGKScnx5L9Aqg8CNsA4KaTJ09Kyh8ZlvJD8XnnnSdJ+uCDDyq8/wYNGkiSPv74Y+Xm5lZ4f64UvJ6CKSSzZs3Se++9Jyl/6kn37t1lGIb+85//aNCgQSUGbil/hH/YsGFavHixWrRooR49eigsLEyS9Ntvv2nAgAGaNWuWTp06VRiOt2zZovvuu0+zZ88u3E+LFi3Up08fSVK9evV00003Ff5XMOdckn766ScNGjRIS5YsUUREhHr37q26desqPj5egwYN0vr16yvcl3vvvVfffPONBg0apOnTpzPNBoDbmLMNAG4wDEPfffedpL9WLAkJCdGAAQO0cOFCvfLKK1q2bJmuueYatW7dWq1bt3Z7tPuWW27RmjVrtGrVKvXu3VvXXnut2rVrp1atWumiiy6y9PXs2bNHBw4cKPJ6br31Vk2cOFH16tUr3C4vL09vvPGGZsyYoX//+9+aOnVqsX1t2bJF7du317fffqtatWoV3u50OvXggw/q4MGDevzxxzVq1CgFBeWP9ezbt0+jRo3Syy+/rCuvvFLR0dHq1auXmjdvrm+++UaXXHJJsXnzUv4Jno899phOnTqlSZMmaejQoYX3zZ07V1OnTtX48eO1bNmywg9C7sjIyNDo0aO1ZcsW3XXXXZZM2wFQOTGyDQAmOJ1OpaSk6Omnn9bGjRsVEhKiQYMGFd7/9NNPa+DAgXI4HNq2bZtmzJih0aNHq2vXrho8eLCWLl1q+rluuOEGPfXUUwoNDdUff/yh9957T4888oiuvfZa9ejRQ7NmzdLp06cr9Hqys7P1008/ady4cXI6neratasaN24sSercuXORoC3lzzsfN26c6tevr5UrV5a634kTJxYJ2pK0atUq7d69W3369NHdd99dGLQlqXHjxpowYYKcTqc++ugj0/V/9dVXOnz4sNq3b18kaEvSiBEj1LJlS/3555+F89/dcejQIQ0bNkxbtmzR+PHjCdoAKoSRbQAow7nrbUtS9erV9cILLxQZZa5WrZqmTp2q++67T998842Sk5O1detWHT58WFu2bNEjjzyi5ORkPfPMM6aed8SIERowYICWLVumtWvXatu2bUpJSVFqaqpefvllrVixQu+9956qVatm+rUsXry4xFVMWrVqpRdffLHIbRkZGVq5cqV++eUXHT9+XHl5eZKkM2fO6NixYzp27FixE0EjIiLUunXrYvtfvXq1JOnaa68tsa6CqSFbt241/VoKpojceOONJd7fr18/bd++XevXr3dricR9+/bpqaee0sGDB/WPf/xDt9xyi+nHAkBJCNsAUIaCdbbPvqhN7969C+cin6tx48YaPXp04d+3b9+uGTNmaNWqVZo/f76uu+66IvOOy1K7dm3dcssthYEvNTVVCxYs0Ny5c7Vp0ya9++67GjNmjOnXcvY622df1KZbt25FRpuXLFmiZ599VtnZ2aXu6+TJk8XC9gUXXFDitqmpqZKkxx57TI899lip+yw4MdSMtLQ0SVJUVFSJ91944YVFtjNr8uTJOnPmjB577DGCNgBLELYBoAwlzRd2R8uWLfXGG29o8ODB2rZtmxISEkyH7XNFRUXpiSeekNPp1Ny5c5WQkOBW2D53ne2SpKamFq748fTTT+vqq69W/fr1C0fQhwwZoo0bN8owjGKPLW1udMGo+JVXXllsesrZPLH6iivXXXedlixZovfee089evQocelFAHAHYRsAbBYUFKSOHTtq27Ztbo3elqZz586aO3euJfs6V0JCgnJzczVq1Cjdeeedxe4vOJnSHQWrqwwePLhwlZGKioyMlPTXqPm5Cm4v2M6swYMHKyYmRpMnT9add96p+fPnFy6NCADlwQmSAOAB+/btkyTVr1/f5bYljRqXtC93g6QZx48fl/RXQD5bUlKSDh8+7PY+u3XrJklavny56ccULKtYsH75uTp06CBJ+vLLL0u8//PPPy+ynTtuv/12TZo0Senp6Ro+fHhhvwGgPAjbAFBBx48f180336yvv/662EVP8vLytGjRIq1cuVJBQUHq1auXy/2NGTNG8+bN07Fjx4rdt3nzZr355puSpL/97W+W1H+2Jk2aSMoPq2fP2T506JCee+65cu2zd+/euvTSS/XFF19o5syZxXpkGIY2bNigDRs2FN4WHh6uqlWr6sCBA3I6ncX2ed1116levXrasGGDFi5cWOS+efPmadu2bapfv365R9KHDh2qiRMnKi0tTcOHD9f+/fvLtR8AYBoJAFhg69ateuihhxQaGqpWrVopMjJSJ0+e1O7du5WamiqHw6FHHnlEzZs3d7mvgwcPasqUKXrhhRfUvHlzXXjhhTIMQ/v379fOnTslSddcc41uvfVWy19Hjx49dNlll2nbtm3q3bu3YmJidPr0aa1du1bNmzdX+/bttXHjRrf2WaVKFc2cOVN33XWXXnvtNS1YsEDNmjVTnTp1dOzYMe3cuVNHjhzRU089VTifPSQkRHFxcVq1apX69++vyy+/XFWrVlVMTIwGDRqk0NBQ/etf/9J9992nSZMmaeHChbr44ou1d+9e7dixQ6GhoXrllVfKtcZ2geHDhysvL09Tp07V8OHDNX/+fDVq1Kjc+wNQORG2AaCCatasqYULF2r16tVau3atUlNTtXnzZkn500b69++vIUOGKCYmxtT+XnvtNSUmJuqHH37Qb7/9psTEROXm5qp27dq6+uqrdeONN+qGG26Qw+Gw/LWEhIRowYIFevXVV5WYmKhVq1apfv36GjZsmO6///4iK624o0mTJvrss8/0/vvva/ny5dq0aZOcTqfq1atXeLXJ6667rshjCj5w/Pjjj1qyZImcTqecTmfh+uZdunTRxx9/rLfeektr1qzR7t27Vbt2bfXr109jxozRJZdcUuF+jBgxQnl5eXrhhRd05513at68eYUrnQCAGQ7D1eRAAAAAAOXCnG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAmhG0AAADAJoRtAAAAwCaEbQAAAMAm/x/jnPrU49cMJQAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 720x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "psis = az.loo(results_5_3, pointwise=True, scale=\"deviance\")\n",
+    "waic = az.waic(results_5_3, pointwise=True, scale=\"deviance\")\n",
+    "\n",
+    "plt.scatter(psis.pareto_k, waic.waic_i)\n",
+    "plt.xlabel(\"PSIS Pareto k\")\n",
+    "plt.ylabel(\"WAIC\");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Code 7.35\n",
+    "\n",
+    "Bambi does not support t distribution yet."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Last updated: Fri May 21 2021\n",
+      "\n",
+      "Python implementation: CPython\n",
+      "Python version       : 3.8.5\n",
+      "IPython version      : 7.18.1\n",
+      "\n",
+      "statsmodels: 0.12.2\n",
+      "scipy      : 1.5.4\n",
+      "matplotlib : 3.3.3\n",
+      "arviz      : 0.11.2\n",
+      "numpy      : 1.20.1\n",
+      "pandas     : 1.2.2\n",
+      "bambi      : 0.5.0\n",
+      "\n",
+      "Watermark: 2.1.0\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "%watermark -n -u -v -iv -w"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "bmb",
+   "language": "python",
+   "name": "bmb"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/Rethinking/data/cars.csv b/Rethinking/data/cars.csv
new file mode 100644
index 0000000..6908989
--- /dev/null
+++ b/Rethinking/data/cars.csv
@@ -0,0 +1,51 @@
+speed,dist
+4,2
+4,10
+7,4
+7,22
+8,16
+9,10
+10,18
+10,26
+10,34
+11,17
+11,28
+12,14
+12,20
+12,24
+12,28
+13,26
+13,34
+13,34
+13,46
+14,26
+14,36
+14,60
+14,80
+15,20
+15,26
+15,54
+16,32
+16,40
+17,32
+17,40
+17,50
+18,42
+18,56
+18,76
+18,84
+19,36
+19,46
+19,68
+20,32
+20,48
+20,52
+20,56
+20,64
+22,66
+23,54
+24,70
+24,92
+24,93
+24,120
+25,85