"
],
"text/plain": [
" rank waic p_waic d_waic weight se \\\n",
"m6.7 0 337.244430 3.308165 0.000000 1.000000e+00 11.844569 \n",
"m6.8 1 399.758470 3.089429 62.514039 0.000000e+00 14.941817 \n",
"m6.6 2 409.200795 1.712087 71.956364 7.406298e-12 12.402004 \n",
"\n",
" dse warning waic_scale \n",
"m6.7 0.000000 True deviance \n",
"m6.8 13.865603 True deviance \n",
"m6.6 12.789497 False deviance "
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post = m6_6.sample_posterior(random.PRNGKey(77), p6_6, sample_shape=(1000,))\n",
"logprob = log_likelihood(m6_6.model, post, h0=d.h0.values, h1=d.h1.values)\n",
"az6_6 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"post = m6_7.sample_posterior(random.PRNGKey(77), p6_7, sample_shape=(1000,))\n",
"logprob = log_likelihood(\n",
" m6_7.model,\n",
" post,\n",
" treatment=d.treatment.values,\n",
" fungus=d.fungus.values,\n",
" h0=d.h0.values,\n",
" h1=d.h1.values,\n",
")\n",
"az6_7 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"post = m6_8.sample_posterior(random.PRNGKey(77), p6_8, sample_shape=(1000,))\n",
"logprob = log_likelihood(\n",
" m6_8.model, post, treatment=d.treatment.values, h0=d.h0.values, h1=d.h1.values\n",
")\n",
"az6_8 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"az.compare({\"m6.6\": az6_6, \"m6.7\": az6_7, \"m6.8\": az6_8}, ic=\"waic\", scale=\"deviance\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.27"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"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",
"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"
]
},
{
"data": {
"text/plain": [
"DeviceArray(13.789175, dtype=float32)"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post = m6_7.sample_posterior(random.PRNGKey(91), p6_7, sample_shape=(1000,))\n",
"logprob = log_likelihood(\n",
" m6_7.model,\n",
" post,\n",
" treatment=d.treatment.values,\n",
" fungus=d.fungus.values,\n",
" h0=d.h0.values,\n",
" h1=d.h1.values,\n",
")\n",
"az6_7 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"waic_m6_7 = az.waic(az6_7, pointwise=True, scale=\"deviance\")\n",
"post = m6_8.sample_posterior(random.PRNGKey(91), p6_8, sample_shape=(1000,))\n",
"logprob = log_likelihood(\n",
" m6_8.model, post, treatment=d.treatment.values, h0=d.h0.values, h1=d.h1.values\n",
")\n",
"az6_8 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"waic_m6_8 = az.waic(az6_8, pointwise=True, scale=\"deviance\")\n",
"n = waic_m6_7.n_data_points\n",
"diff_m6_7_m6_8 = waic_m6_7.waic_i.values - waic_m6_8.waic_i.values\n",
"jnp.sqrt(n * jnp.var(diff_m6_7_m6_8))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.28"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"DeviceArray([12.960003, 67.03999 ], dtype=float32, weak_type=True)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"40.0 + jnp.array([-1, 1]) * 10.4 * 2.6"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.29"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
"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",
"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"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAE3CAYAAADmP0YBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqgElEQVR4nO3de1TVdb7/8RfKxQuIKOC18lJSKl6ZDCVvWJmpR/SklTInp/JyMi31eEuT0llaqaGWqWPTJDIlZ0TL8bJcWFqG2liaeR09YaPhBRGUSLnu3x/+3CMqXyDZfDZfn4+1XEu++3t5s99sePH9fPYHD4fD4RAAAACMqGK6AAAAgDsZYQwAAMAgwhgAAIBBhDEAAACDCGMAAAAGEcYAAAAMIowBAAAYRBgDAAAwyNN0AeUpIyPDdAm2891330mSOnTo8JuO9/f318WLF8uzJJQj+uO+6I17oz/uy916ExAQUOI+3BmDpZSUFKWkpPzm46tU4UvMndEf90Vv3Bv9cV+VsTeVr2IAAAAbIYwBAAAYRBgDAAAwyFYT+FH+KuPYOwAAlQlhDJaioqJMlwAAgK1x2wMAAMAgwhgs7dmzR3v27DFdBgAAtkUYg6WffvpJP/30k+kyAACwLcIYAACAQYQxAAAAgwhjAAAABrG0BSx5eXmZLgEAAFsjjMFS//79TZcAAICtMUwJAABgEGEMlnbv3q3du3ebLgMAANsijMHSqVOndOrUKdNlAABgW4QxAAAAgwhjAAAABhHGAAAADGJpC1jy8fExXQIAALZGGIOlvn37mi4BAABbY5gSAADAIMIYLO3cuVM7d+40XQYAALbFMCUspaammi4BAABb484YAACAQYQxAAAAgwhjAAAABjFnDJZq1KhhugQAAGyNMAZLjz/+uOkSAACwNYYpAQAADCKMwdKOHTu0Y8cO02UAAGBbDFPC0tmzZ02XAACArXFnDAAAwCDCGAAAgEGEMQAAAIOYM+aGcnJytHXrVp09e1b16tVTZGSkfHx8jNRx5MgRZWRkqHr16sbqAADAzghjbiYuLk4xMTHKyMhwbgsICFBMTIyio6MrtI5Zs2bp/PnzkqQ///nPCgwM1IwZMyq0DgAA7K5CwthXX32ld999V0eOHJGXl5dCQ0P14YcfFrt/z5499fPPP9+0vVu3blq+fLkrSzUqLi5O48aNu2l7RkaGc3tFBKFrdXTu3FlLlixRq1atdPDgQcXGxlZoHQAA3Ak8HA6Hw5UXSEpK0tSpU/Xyyy8rPDxcDodDhw8fVt++fYs95sKFCyooKHB+nJaWpoEDB2rOnDmKiooq9rjr7yZVNjk5OWrZsqXl51CnTh0dPHjQpUOFOTk5Cg0NVYsWLbR27Vrt3LlTktS1a1fl5eUpKipKx44d0/79+0tVR0BAQKXui93RH/dFb9wb/XFf7tabgICAEvcp052x6OhoNWvWTNWrV1diYqI8PDw0atQoDR06VLNnz9aGDRvk5+en8ePHq3///iooKNDs2bM1ceJEDRkyxHme5s2bW16nTp06RT7+29/+Jl9fX/Xu3bss5bpEv379XHLe9PT0Er94Lly4oB49eqhu3bouqeFaHefPn1dwcLAGDhyojh07SpLefvttSVJmZqbS0tJKXYenp6fy8/NdVq8rrV+/3nQJAIA7QJmHKdevX6/hw4crISFBSUlJmjt3rnbs2KGIiAitWbNGiYmJevXVV/XQQw/pzJkzOn36tKpVq6aBAwfq7NmzCgkJ0cSJE9WyZctSXc/hcOhvf/ub+vfvr+rVq1vu6+/vrypVXPsGUU9P14zsXn8nsKT9XFXD9XXUqlVLnp6e8vDwkPTvz7tWrVplrsOV9bpSaX6bsYM75fOsjOiNe6M/7quy9aZMw5TR0dHKzc3V6tWrJV0NSuHh4WrXrp2WLl0qScrNzVX79u01f/58FRQUaPz48WrYsKEmT56sxo0bKz4+Xlu2bNGmTZsUHBxc4jV37Nih5557TuvWrdMDDzxgua873ZYsq40bN2rYsGEl7rdq1Sr16dPH5XUkJCSoV69eWrNmjSRp0KBBkq4OOw8ePLjUdbjb7WIURX/cF71xb/THfblbb0oTDMt8GykkJMT5fw8PD9WtW7fINm9vb/n7+ys9PV2FhYWSpFGjRql3795q3bq1Zs2aJV9fX3366aelul5CQoJCQ0NLDGKVXWRkZIkNq1OnjiIjI11eR2BgoGJjY5WXl1fksby8PMXGxiooKMjldQAAcKcocxi7ccjJw8PjltscDoeCgoIkFZ0j5unpqSZNmig1NbXEa6Wnp+vzzz/X4MGDy1pmpePj46OYmBjLfWbOnOnydb58fHw0Y8YMJScnKyoqSikpKSooKFBSUpKioqKUnJys6dOns94YAADlxKWTeVq3bi1vb2+lpKQoLCxMklRYWKiTJ08qIiKixOMTExPl5eXl0mE5d3JtuQjT64xdu86sWbOUnJzs3B4YGKiFCxeyrAUAAOXIpWHM19dXTz31lBYvXqz69eurUaNGio+P18WLF9W/f3/nfr1799awYcOKzJm6NnH/iSeekK+vryvLdCvR0dEaPHiw8RX43aUOAADszuVvc5s0aZK8vLw0ZcoUXb58Wa1atdLKlStVr1495z4pKSk3TbbbvXu3Tpw44VxS4U7i4+PjFncDfXx8VL16dTVp0kQ9evQwXQ4AALbk8kVfK5I7vXvCLm58N2VZudu7WlAU/XFf9Ma90R/35W69ccm7KQEAAFB+CGMAAAAGEcYAAAAMqpx/pwYV5sa/EwoAAMoXYQyWeBclAACuxTAlAACAQYQxWEpKSlJSUpLpMgAAsC2GKWHp4sWLpksAAMDWuDMGAABgEGEMAADAIMIYAACAQcwZg6XAwEDTJQAAYGuEMVjq1q2b6RIAALA1hikBAAAMIozB0pYtW7RlyxbTZQAAYFsMU8JSVlaW6RIAALA17owBAAAYRBgDAAAwiDAGAABgEHPGYKlevXqmSwAAwNYIY7AUERFhugQAAGyNYUoAAACDCGOwtGnTJm3atMl0GQAA2BbDlLD066+/mi4BAABb484YAACAQYQxAAAAgwhjAAAABjFnDJYaNGhgugQAAGyNMAZLnTt3Nl0CAAC2xjAlAACAQYQxWNqwYYM2bNhgugwAAGyLYUpYunLliukSAACwNe6MAQAAGEQYAwAAMIgwBgAAYBBzxmCpcePGpksAAMDWCGOw1KlTJ9MlAABgawxTAgAAGEQYg6XPPvtMn332mekyAACwLYYpYSkvL890CQAA2Bp3xgAAAAwijAEAABhEGAMAADCIOWOwdPfdd5suAQAAWyOMwdLvfvc70yUAAGBrDFMCAAAYRBiDpbVr12rt2rWmywAAwLYYpoSlwsJC0yUAAGBr3BkDAAAwiDAGAABgEMOUAACU0smTJyVJAQEBhiuBnRDGYKlp06amSwAAt5CVlaV169ZJku655x6zxcBWGKaEpQ4dOqhDhw6mywAA4/bs2SOHwyGHw6GvvvrKdDmwEcIYAAAlyMrK0oEDB5wf7927V1lZWQYrgp0wTAlLa9askSQNGjTIcCUAKsq11z3+LTMzUw6Hw/lxYWGhEhISVLt2bXNFuSF+Vvw2tgpj/v7+qlKFm33lydPz6pfI7UxWZaKre6M/7stUb6697nFVfn6+srOzb9qenZ2t2rVr83xdx12+n7hLHaVlq6+gixcvmi7BdvLz8yVJGRkZv+n4gICA33wsXI/+uC+TvfmP//gPI9d1V1988YV+/vnnWz5Wp04d9ejRo4Ircl/u8P3E3b6vlSYYchsJAIBi3DhX7EYHDhxg7hhuG2EMAIBi7N+/v8hcsRs5HA7t37+/AiuCHdlqmBLl79577zVdAgAYExISojNnzty03dPT0zmNIyQkpKLLgs0QxmCpbdu2pksAAGMCAwNv+Q5Bd5uXhMqNYUpYKigoUEFBgekyAACwLcIYLK1bt8755z8AAED5I4wBAAAYRBgDAAAwiDAGAABgEGEMAADAIJa2gCXWzwEAwLUIY7DUunVr0yUAAGBrDFPCUm5urnJzc02XAQCAbRHGYGn9+vVav3696TIAALAthikBAG4vJydHW7du1dmzZ1WvXj1FRkbKx8fHWA3NmzdXp06dKrwG2BNhDADg1uLi4hQTE1Pkb0EGBAQoJiZG0dHRFVbDrFmzdP78eee2wMBAzZgxo8JqgH1VyDDlV199pSFDhqht27YKCwvT8OHDLfcvKChQbGysevbsqdDQUPXs2VPvvPOO8vPzK6JcAICbiIuL07hx4276o9wZGRkaN26c4uLiKqyGFi1aKCEhQQcPHtSmTZvUokWLCqsB9ubhcDgcrrxAUlKSpk6dqpdfflnh4eFyOBw6fPiw+vbtW+wxS5cu1QcffKA5c+bo/vvv15EjRzRlyhQNHz5cL774YrHH3fhixe1bs2aNJGnQoEG/6fiAgAD64sboj/uiN1eHBVu2bGn5PNSpU0cHDx502XBhTk6OQkND1aJFC61du1ZeXl6Srvbn3LlzioqK0rFjx7R//36GLN2Eu712AgICStynTMOU0dHRatasmapXr67ExER5eHho1KhRGjp0qGbPnq0NGzbIz89P48ePV//+/VVQUKDZs2dr4sSJGjJkiPM8zZs3t7zO3r171bNnT/Xq1UuS1LhxY/Xs2VP79+8vS7koBw888IDpEgC4UL9+/UyXUKz09PQSf6heuHBBPXr0UN26dV1Ww/nz5xUcHKyBAwc6t3t6eio/P1+ZmZlKS0tzaQ12wBvBrJV5ztj69es1fPhwJSQkKCkpSXPnztWOHTsUERGhNWvWKDExUa+++qoeeughnTlzRqdPn1a1atU0cOBAnT17ViEhIZo4caJatmxZ7DU6duyov/71r/q///s/NW/eXMePH9fu3bs1YsQIy9r8/f1VpQpvEC1PXbp0ue1zlOa3AphDf9xXRfTG09N9pw4XFBSUej9XfR7XaqhVq9ZN1/D09FStWrVcXoMdVPT3mcr2fa1Mw5TR0dHKzc3V6tWrJUkOh0Ph4eFq166dli5dKunqulTt27fX/PnzVVBQoPHjx6thw4aaPHmyGjdurPj4eG3ZskWbNm1ScHDwLa/jcDi0cOFCLV26VFWrVlV+fr5Gjhyp8ePHW9bnTrcl7eLKlSuSpGrVqv2m493tdjGKoj/ui95IGzdu1LBhw0rcb9WqVerTp49La0hISHCO1kj/7k9SUpIGDx7s0hpQNu722ilNMCzzbaTr/zyOh4eH6tatW2Sbt7e3/P39lZ6ersLCQknSqFGj1Lt3b7Vu3VqzZs2Sr6+vPv3002KvsXnzZq1du1bz5s1TYmKi5s+fr//93/9lkqQBGzZs0IYNG0yXAeAOFBkZWeIPsjp16igyMtKlNQQGBio2NlZ5eXlFHsvLy1NsbKyCgoJcWgPsr8xh7MbbsB4eHrfc5nA4FBQUJKnoHDFPT081adJEqampxV5j7ty5+sMf/qC+ffsqJCREffv21ciRI/X++++XtVwAQCXl4+OjmJgYy31mzpzp0onzPj4+mjFjhpKTkxUVFaWkpCSlpqZq8+bNioqKUnJysqZPn87kfdwWlw5wt27dWt7e3kpJSVFYWJgkqbCwUCdPnlRERESxx125cuWmuV9Vq1Z1ZakAADd0bQ0vk+uMXbvGrFmzNHjwYOf2wMBALVy4kHXGcNtcGsZ8fX311FNPafHixapfv74aNWqk+Ph4Xbx4Uf3793fu17t3bw0bNsw5N6Bnz55atmyZGjVqpBYtWujo0aNatmxZkfF6AMCdITo6WoMHDza6Av+NNbACP8qTy9/6MWnSJHl5eWnKlCm6fPmyWrVqpZUrV6pevXrOfVJSUor8xjN9+nQtWrRIs2fPVlpamoKCgtSvXz+99NJLri4XAOCGfHx8jE+Qv74Gd5skjsrN5Yu+ViReGOXv6NGjkoq+caMs+Ibl3uiP+6I37o3+uC936025L/qKO89vDWEAAKB0WCEVlrKzs5WdnW26DAAAbIswBkubN2/W5s2bTZcBAIBtEcYAAAAMIowBAAAYRBgDAAAwiDAGAABgEEtbwFKbNm1MlwAAgK0RxmDpvvvuM10CAAC2xjAlLGVlZSkrK8t0GQAA2BZhDJa2bNmiLVu2mC4DAADbIowBAAAYRBgDAAAwiDAGAABgEGEMAADAIJa2gKX27dubLgEAAFsjjMFSs2bNTJcAAICtMUwJS5mZmcrMzDRdBgAAtkUYg6WtW7dq69atpssAAMC2CGMAAAAGEcYAAAAMIowBAAAYRBgDAAAwiKUtYKljx46mSwAAwNYIY7DUpEkT0yUAAGBrDFPCUnp6utLT002XAQCAbRHGYGnbtm3atm2b6TIAALAtwhgAAIBBhDEAAACDCGMAAAAGEcYAAAAMYmkLWHrwwQdNlwAAgK0RxmDprrvuMl0CAAC2xjAlLKWlpSktLc10GQAA2BZhDJa+/PJLffnll6bLAADAtghjAAAABhHGAAAADCKMAQAAGEQYAwAAMIilLWApPDzcdAkAANgaYQyWGjZsaLoEAABsjWFKWDpz5ozOnDljugwAAGyLMAZLX3/9tb7++mvTZQAAYFuEMQAAAIMIYwAAAAYRxgAAAAwijAEAABjE0haw1KVLF9MlAABga4QxWKpfv77pEgAAsDWGKWEpNTVVqamppssAAMC2CGOwtHPnTu3cudN0GQAA2BZhDACAO8TJkyd18uRJ02XgBoQxAADuAFlZWVq3bp3WrVunrKws0+XgOoQxAADuAHv27JHD4ZDD4dCePXtMl4PrEMYAALC5rKwsHThwwPnxgQMHuDvmRljaApa6du1qugQAKBdr1qwpt3N5enoqPz+/3M7napmZmXI4HM6PHQ6HEhISVLt2bXNFuUhZejNo0CAXV1M6tgpj/v7+qlKFm33lKSAgwC3OAdehP+6L3pQvT8/y/ZFX3udzlfz8fGVnZ9+0PTs7W7Vr1640n0dZlPZzcpfXmIfj+qhcyWVkZJguwXauvevmrrvu+k3HBwQE0Bc3Rn/cF71xb5WpP1988YV++OGHWz4WGhqqHj16VHBFruVuvSlN4OM2Eix98803+uabb0yXAQD4DW6cK3Yj5o65B8IYAAA2tX//flkNgDkcDu3fv78CK8Kt2G+gGAAASJJCQkJ05syZEveBWYQxAABsKjAw0G3eMYjiMUwJAABgEHfGYKl79+6mSwAAwNYIY7BUt25d0yUAAGBrDFPC0okTJ3TixAnTZQAAYFvcGYOlb7/9VpLUpEkTs4UAAGBT3BkDAAAwiDAGAABgEGEMAADAIMIYAACAQUzgh6XIyEjTJQAAYGuEMViqXbu26RIA4I6Xk5OjrVu36uzZs6pXr54iIyPl4+NT4eeAa1RIGPvqq6/07rvv6siRI/Ly8lJoaKg+/PBDy2POnTun+fPna/v27crOztbdd9+t119/XWFhYRVRMv6/H3/8UZLUrFkzw5UAwJ0pLi5OMTExysjIcG4LCAhQTEyMoqOjS32OWbNm6fz5885tgYGBmjFjRqnPAddxeRhLSkrS1KlT9fLLL2vOnDlyOBw6fPiw5TGXLl3S008/rY4dO2r58uUKCAjQqVOnFBAQ4OpycYO9e/dKIowBgAlxcXEaN27cTdszMjKc20sKU9fO0blzZy1ZskStWrXSwYMHFRsbW+pzwLU8HA6Ho7Q7R0dHq1mzZqpevboSExPl4eGhUaNGaejQoZo9e7Y2bNggPz8/jR8/Xv3791dBQYEiIyM1evRoDRkypNRFLViwQN98840++eSTMn0y1//WgPKxZs0aSdKgQYN+0/EBAQH0xY3RH/dFb9xbRfQnJydHLVu2tLxOnTp1dPDgwWKHG3NychQaGqoWLVpo7dq18vLycj6Wl5enqKgoHTt2TPv377fNkKW7vXZKcyOpzHfG1q9fr+HDhyshIUFJSUmaO3euduzYoYiICK1Zs0aJiYl69dVX9dBDD+nMmTM6ffq0qlWrpoEDB+rs2bMKCQnRxIkT1bJly2KvkZSUpK5du2rChAlKTk5WcHCwnnzySQ0dOlQeHh5lLRkAAKd+/frd9jk8PT2Vn59fDtUULz09vcRQceHCBfXo0aPYvyOcnp6u8+fPKzg4WAMHDrzp8czMTKWlpVmeo7IpTW/Wr19fQdWUTpnD2H333aeXXnpJkvTcc89pxYoV8vLy0vDhwyVJY8aM0QcffKDvvvtOBQUFkqTY2FhNnjxZjRs3Vnx8vKKjo7Vp0yYFBwff8honT55UfHy8nn32WX3wwQc6fPiwZs+eLQ8PDw0dOrTY2vz9/VWlCqt1lCdPz6tfIrczRMzwsnujP+6L3rjGte9r7nKe4lz7GVqa/Yqr5do5atWqdct9atWqVeI5KqOSPhd3e22V+ZkPCQlx/t/Dw0N169Ytss3b21v+/v5KT093NnnUqFHq3bu3JGnWrFlKTk7Wp59+qhdeeOGW13A4HGrVqpUmTJggSWrZsqVOnDih+Ph4yzB28eLFsn46KMG13y5+6y1fd7tdjKLoj/uiN66zdu3a2z5HRfRn48aNGjZsWIn7zZw5U3369LE8x/jx49WrV6+bHk9KStLgwYMtz1HZlKY3FfnaKk3wK/NtpBvTpoeHxy23ORwOBQUFSZKaN29e5PgmTZooNTW12GsEBQXp3nvvLbKtefPmlsfANR599FE9+uijpssAgDtOZGRkiT/I69SpY7keZGRkpAIDAxUbG6u8vLwij+Xl5Sk2NlZBQUGsKWmYS8f0WrduLW9vb6WkpDi3FRYW6uTJk2rYsGGxx3Xo0KHIMZJ04sQJy2PgGn5+fvLz8zNdBgDccXx8fBQTE2O5z8yZMy0n3vv4+GjGjBlKTk5WVFSUkpKSlJqaqqSkJEVFRSk5OVnTp0+3zeT9ysqlA8S+vr566qmntHjxYtWvX1+NGjVSfHy8Ll68qP79+zv36927t4YNG+a8Hftf//Vfevrpp/X++++rT58+OnTokOLi4jR+/HhXlotbOHbsmKSrcwUBABXr2pITt7PO2LV9Zs2apcGDBzu3BwYGauHChSxr4QZcPltv0qRJ8vLy0pQpU3T58mW1atVKK1euVL169Zz7pKSkFPkia9Omjd577z0tWLBAS5YsUcOGDTVu3Dg988wzri4XN9i/f78kwhgAmBIdHa3Bgwff1ur55XEOuE6Z1hlzd0x2LX+sM2Zv9Md90Rv3Rn/cl7v1xiUT+AEAAFB+CGMAAAAGEcYAAAAMss9yu3CJa4v1AgAA1yCMwVLNmjVNlwAAgK0xTAlLR48e1dGjR02XAQCAbXFnDJYOHDggqejfJAUAAOWHO2MAAAAGEcYAAAAMIowBAAAYRBgDAAAwiAn8sPTEE0+YLgEAAFsjjMFStWrVTJcAAICtMUwJS4cOHdKhQ4dMlwEAgG0RxmDp8OHDOnz4sOkyAACwLcIYAACAQYQxAAAAgwhjAAAABhHGAAAADGJpC1jq16+f6RIAALA1whgseXt7my4BAABbY5gSlg4cOKADBw6YLgMAANsijMHS0aNHdfToUdNlAABgW4QxAAAAgwhjAAAABhHGAAAADCKMAQAAGMTSFrA0YMAA0yUAAGBrhDFYqlq1qukSAACwNYYpYen777/X999/b7oMAABsizAGS8ePH9fx48dNlwEAgG0RxgAAAAwijAEAABhEGAMAADCIMAYAAGCQh8PhcJguAgAA4E7FnTEAAACDCGMAAAAGEcYAAAAMIowBAAAYRBgDAAAwiDCGMklMTNSAAQMUFhamdu3aacCAAVq3bt0t950xY4ZCQkL0wQcfFNmem5urN954Q506dVL79u314osv6ty5cxVQvb2VpjcpKSkaM2aMwsLC1LZtW0VFRSklJcX5OL1xnZL6k52drVmzZqlr165q06aNHnvsMf3lL38pcg76UzGWLl2qkJAQvfHGG85tDodDixYtUkREhNq0aaPf//73+vHHH4scR38qxo39ycvL09tvv61+/fqpXbt2ioiI0IQJE5SamlrkOHfuD2EMZVK3bl299NJLWr16tT777DMNGDBA06ZN0/bt24vst3nzZv3www8KDg6+6Rx//OMflZSUpNjYWMXFxen8+fN68cUXxSort6ek3pw8eVJPP/20GjdurI8++kh///vfNW7cONWoUcN5DnrjOiX1Z+7cufriiy/05ptvasOGDRoxYoTmzZtXJLDRH9fbt2+fEhISFBISUmT7n/70J3300UeaOXOmEhISVLNmTT333HO6cuWKcx/643q36s+VK1d06NAhjR49WomJiVqyZIlOnz6t559/Xvn5+c793Lo/DuA2DRgwwDFv3jznx6dOnXJEREQ4jh8/7ujRo4djxYoVzscuXbrkaNWqlWP9+vXObcePH3e0aNHC8Y9//KNC674TXN+b8ePHO8aPH1/svvSm4l3fnyeeeMKxcOHCIo8/88wzjtdff93hcNCfinDp0iVHZGSkIzk52TFs2DDnc19YWOjo0qWLY+nSpUX2bd26tWPdunXOj+mPaxXXn1s5duyYo0WLFo4jR444j3Xn/nBnDL9ZYWGhvv76a6WkpCgsLEySlJ+frwkTJmj06NFq3rz5TcccOHBAeXl56ty5s3Nb8+bN1aBBA+3bt6+iSre9G3tTWFiozz//XPfdd59GjBih8PBwDRo0SBs3bnQeQ28qzq1eOx06dNAXX3yh06dPS5K+++47HTlyRA8//LAk+lMRZsyYoccee0zh4eFFtp86dUppaWlFnns/Pz+1adPG+dzTH9crrj+38ssvv0iSateuLcn9++NpugBUPqmpqXriiSeUm5urKlWq6LXXXlO3bt0kSYsXL1bt2rX1zDPP3PLY8+fPq2rVqgoICCiyPTAwUOfPn3d57XZXXG/S0tL066+/aunSpRo3bpxeeeUV7dq1SxMnTlSNGjXUvXt3elMBrF4706dPV0xMjLp37y5Pz6vfmqdNm6YePXpI4rXjagkJCfrXv/6lt95666bH0tLSJF0dar7e9c89/XEtq/7cKDc3V3PnzlWPHj1Ur149Se7fH8IYyiw4OFjr1q1Tdna2vv76a82ZM0eNGzdW1apVlZiYqE8//bTM53Q4HPLw8HBBtXeW4nrTrFkzSVJkZKSGDx8uSXrggQf0ww8/KD4+Xt27dy/2nPSm/BTXn/DwcP31r3/Vt99+qyVLlqhx48bat2+f3nzzTQUHB+uRRx4p9pz05/b9+OOPWrBggeLj4+Xt7V3sfjc+z6V57unP7SttfySpoKBAkyZN0sWLF7VkyZISz+0u/SGMocw8PT11zz33SJJatmypEydOaNmyZerYsaPS0tIUERHh3LegoEDz5s3TRx99pC+//FKBgYEqKChQRkaG6tSp49zvwoULN/3WibIrrjfLly+Xp6fnTUPHzZs3dw5V0hvXK64/7du317x58xQbG6vIyEhJUkhIiI4fP67ly5frkUceoT8utG/fPmVkZKhfv37ObQUFBfrHP/6hTz75RH//+98lXb270qBBA+c+Fy5c0H333SeJ148rldSfffv2ydvb2zlN5siRI4qLiyvSB3fvD2EMt62wsFA5OTl65pln9NhjjxV57LnnnlPfvn315JNPSpJat24tLy8v7dq1S3369JF0dbmF1NRUtWvXrqJLt71rvfH29lZoaGiRZSwk6cSJE2rYsKEkemPCtf7k5+crPz9fVaoUncZbtWpV5//pj+v06tVLrVu3LrJt6tSpatKkiUaOHKmmTZsqKChIO3fuVGhoqKSrc5K+//57vrdVgJL64+Xlpby8PL3yyis6duyYVq5cedM7+d29P4QxlMnChQvVqVMnNWrUSDk5Odq+fbs+++wzvfbaa6pbt+5Nv2F4eXkpMDDQOUzm5+enQYMG6a233lJAQID8/Pw0e/ZstWnTRh07djTxKdmGVW8k6fnnn9fLL7+ssLAwPfTQQ9q9e7c2btyo9957TxK9cTWr/vj6+urBBx/UvHnzVL16dTVu3FjfffedVq9erVGjRkmiP65Uq1Yt1apVq8i2GjVqyN/fXy1atJAk/f73v9fy5cvVrFkz3XXXXVq4cKECAwOdv4DSH9cpqT/5+fkaO3asDhw4oKVLl6pKlSrOeX5+fn6qVq2a2/eHMIYyuXTpkl599VWdO3dONWrUUNOmTfX22287f9MojWnTpqlKlSoaN26ccnNz1blzZ8XExLjFuH1lVlJvevXqpTfeeEPLli3TH//4R91zzz168803i8wXozeuU1J/FixYoAULFmjy5MnKzMxUgwYNNHLkSD3//PPOc9Afc1544QVdvnxZM2fO1KVLl9S+fXutWLFC1apVc+5Df8w4c+aMPv/8c0nSwIEDizw2Z84c5zZ37o+Hw+EOq50BAADcmVhnDAAAwCDCGAAAgEGEMQAAAIMIYwAAAAYRxgAAAAwijAEAABhEGAMAADCIMAYAAGAQYQyAW1m8eLFCQkIUEhKi+++/X7/73e80aNAgvfPOO84/cVLeQkJCtGrVKpecGwBKwp9DAuB2/Pz8tGLFCklSVlaWDh06pI8//lirV6/WihUrbvqjwbdr9erVaty4cbmeEwBKiz+HBMCtLF68WKtWrdLu3buLbL906ZKGDh2qK1euaPPmzapataqhCgGgfDFMCaBSqFWrliZOnKh//etf2rFjhyQpJydHb731lrp166bWrVurf//+2r59u/OYyZMn6z//8z9vOtfHH3+sNm3aKDs7W9LNw5Tbtm3T8OHDFR4erg4dOujJJ5/Ul19+WeQcixcvVqdOnXTw4EENHjxYbdu21YABA7Rnz56brpeQkKB+/fopNDRUnTt31tixY5WVleV8fM+ePRo2bJjatm2rTp06afr06frll19u7wkDUGkQxgBUGuHh4fL09NT+/fslSWPHjtXatWs1cuRILV26VK1atdLo0aN1+PBhSVKfPn30ww8/6OTJk0XOs3HjRnXv3l01a9a85XVOnTqlnj176q233tLixYvVoUMHjRgx4qagdeXKFU2ZMkVDhgzRokWL5O3trRdffFGXL1927rNkyRK99tprevDBB/Xee+8pJiZGvr6++vXXXyVJ3377rZ599lkFBgZq0aJFmjp1qrZt26Zp06aV2/MGwL0xZwxApeHt7a2AgACdP39eO3fu1LZt2xQXF6cHH3xQkhQREaETJ07o/fff16JFi9SlSxfVrl1bmzZt0ogRIyRJ586d0549e/TOO+8Ue51hw4Y5/19YWKhOnTrp+PHjWrNmjcLCwpyPXblyRdOmTVN4eLgkKSgoSFFRUfrmm2/UrVs3Xbp0ScuWLdPw4cM1efJk53GPPvqo8//z589X+/btFRsb69wWHBys4cOH65///KdatGhxe08aALfHnTEAlcq1aa7JyckKCgpShw4dlJ+f7/zXpUsXHThwQJLk6empRx99VBs3bnQev3nzZlWrVk3du3cv9hqpqan6n//5Hz388MNq2bKlWrVqpR07dujEiRNF9vPy8lKnTp2cH997772SpLNnz0qS9u7dqytXrmjAgAG3vM7ly5e1b98+Pf7440U+h7CwMHl5eengwYNlfXoAVELcGQNQaeTk5CgzM1OBgYE6e/as0tLS1KpVq5v2u35yf58+fZSQkKCUlBQ1bdpUGzduVM+ePVWtWrVbXqOwsFCjR4/W5cuXNXbsWN1zzz2qXr26Fi1apPT09CL71qxZU1Wq/Pt3Wm9vb0lSbm6uJCkzM1PS1Ttmt3Lp0iUVFBTo9ddf1+uvv37T46dPn7Z4NgDYBWEMQKWxa9cu5efnq127dtq1a5fq1aun9957z/KYTp06KSgoSBs3blRUVJT27dunF154odj9f/rpJx05ckR/+tOf1LVrV+f2K1eulLne2rVrS5LS0tJUp06dmx738/OTh4eHxowZo27dut30eHBwcJmvCaDyIYwBqBQuXbqkefPm6e6771bnzp3l4eGhDz/8UDVq1FDz5s2LPa5KlSp67LHHtGnTJvn4+MjX11cPP/xwsfvn5ORIujoEec3PP/+svXv3lnn+Vvv27VWtWjWtW7euyJyxa2rUqKF27dopJSVFY8aMKdO5AdgHYQyA2ykoKNC+ffskSdnZ2Tp48KA+/vhjXb58WStWrFDVqlXVpUsXRURE6A9/+INeeOEF3Xvvvfrll1905MgR5eTkaMKECc7zPf7441q1apX+8pe/6JFHHnEOJ95Ks2bNVL9+fc2dO1evvPKKsrOztWjRot90l6pWrVr67//+b73zzjvKy8tT165dlZubq+3bt2vMmDGqV6+eJk6cqGeffdYZGmvWrKnTp09r27ZteuWVV9S0adMyXxdA5UIYA+B2srKyNGTIEHl4eMjX11d33323+vfvr2HDhjnnX3l4eOjdd9/V0qVL9dFHH+n06dOqXbu27r///iLvhpSkjh07qkGDBjp9+rT69OljeW1vb28tXrxYb7zxhsaOHasGDRpo9OjR2rVrl/75z3+W+XMZOXKk/P39tXLlSn3yySfy9/dXWFiYc1mNsLAwxcfHa9GiRZo0aZIKCwvVsGFDPfzwwwoMDCzz9QBUPqzADwAAYBBLWwAAABhEGAMAADCIMAYAAGAQYQwAAMAgwhgAAIBBhDEAAACDCGMAAAAGEcYAAAAMIowBAAAYRBgDAAAwiDAGAABg0P8DAIqnwX2NyMgAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"compare = az.compare(\n",
" {\"m6.6\": az6_6, \"m6.7\": az6_7, \"m6.8\": az6_8}, ic=\"waic\", scale=\"deviance\"\n",
")\n",
"az.plot_compare(compare)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.30"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"DeviceArray(7.524193, dtype=float32)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post = m6_6.sample_posterior(random.PRNGKey(92), p6_6, sample_shape=(1000,))\n",
"logprob = log_likelihood(m6_6.model, post, h0=d.h0.values, h1=d.h1.values)\n",
"az6_6 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"waic_m6_6 = az.waic(az6_6, pointwise=True, scale=\"deviance\")\n",
"diff_m6_6_m6_8 = waic_m6_6.waic_i.values - waic_m6_8.waic_i.values\n",
"jnp.sqrt(n * jnp.var(diff_m6_6_m6_8))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.31"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"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",
"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"
]
},
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
m6.6
\n",
"
m6.7
\n",
"
m6.8
\n",
"
\n",
" \n",
" \n",
"
\n",
"
m6.6
\n",
"
0.0
\n",
"
12.7080345
\n",
"
7.5581884
\n",
"
\n",
"
\n",
"
m6.7
\n",
"
12.7080345
\n",
"
0.0
\n",
"
13.690906
\n",
"
\n",
"
\n",
"
m6.8
\n",
"
7.5581884
\n",
"
13.690906
\n",
"
0.0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" m6.6 m6.7 m6.8\n",
"m6.6 0.0 12.7080345 7.5581884\n",
"m6.7 12.7080345 0.0 13.690906\n",
"m6.8 7.5581884 13.690906 0.0"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post = m6_6.sample_posterior(random.PRNGKey(93), p6_6, sample_shape=(1000,))\n",
"logprob = log_likelihood(m6_6.model, post, h0=d.h0.values, h1=d.h1.values)\n",
"az6_6 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"waic_m6_6 = az.waic(az6_6, pointwise=True, scale=\"deviance\")\n",
"post = m6_7.sample_posterior(random.PRNGKey(93), p6_7, sample_shape=(1000,))\n",
"logprob = log_likelihood(\n",
" m6_7.model,\n",
" post,\n",
" treatment=d.treatment.values,\n",
" fungus=d.fungus.values,\n",
" h0=d.h0.values,\n",
" h1=d.h1.values,\n",
")\n",
"az6_7 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"waic_m6_7 = az.waic(az6_7, pointwise=True, scale=\"deviance\")\n",
"post = m6_8.sample_posterior(random.PRNGKey(93), p6_8, sample_shape=(1000,))\n",
"logprob = log_likelihood(\n",
" m6_8.model, post, treatment=d.treatment.values, h0=d.h0.values, h1=d.h1.values\n",
")\n",
"az6_8 = az.from_dict({}, log_likelihood={\"h1\": logprob[\"h1\"][None, ...]})\n",
"waic_m6_8 = az.waic(az6_8, pointwise=True, scale=\"deviance\")\n",
"dSE = lambda waic1, waic2: jnp.sqrt(\n",
" n * jnp.var(waic1.waic_i.values - waic2.waic_i.values)\n",
")\n",
"data = {\"m6.6\": waic_m6_6, \"m6.7\": waic_m6_7, \"m6.8\": waic_m6_8}\n",
"pd.DataFrame(\n",
" {\n",
" row: {col: dSE(row_val, col_val) for col, col_val in data.items()}\n",
" for row, row_val in data.items()\n",
" }\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.32"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 1000/1000 [00:00<00:00, 1334.71it/s, init loss: 2138.6682, avg. loss [951-1000]: 60.6515]\n",
"100%|██████████| 1000/1000 [00:00<00:00, 1365.18it/s, init loss: 962.7464, avg. loss [951-1000]: 67.4809]\n",
"100%|██████████| 1000/1000 [00:00<00:00, 1106.15it/s, init loss: 3201.7393, avg. loss [951-1000]: 60.7879]\n"
]
}
],
"source": [
"WaffleDivorce = pd.read_csv(\"../data/WaffleDivorce.csv\", sep=\";\")\n",
"d = WaffleDivorce\n",
"d[\"A\"] = d.MedianAgeMarriage.pipe(lambda x: (x - x.mean()) / x.std())\n",
"d[\"D\"] = d.Divorce.pipe(lambda x: (x - x.mean()) / x.std())\n",
"d[\"M\"] = d.Marriage.pipe(lambda x: (x - x.mean()) / x.std())\n",
"\n",
"\n",
"def model(A, D=None):\n",
" a = numpyro.sample(\"a\", dist.Normal(0, 0.2))\n",
" bA = numpyro.sample(\"bA\", dist.Normal(0, 0.5))\n",
" sigma = numpyro.sample(\"sigma\", dist.Exponential(1))\n",
" mu = numpyro.deterministic(\"mu\", a + bA * A)\n",
" numpyro.sample(\"D\", dist.Normal(mu, sigma), obs=D)\n",
"\n",
"\n",
"m5_1 = AutoLaplaceApproximation(model)\n",
"svi = SVI(model, m5_1, optim.Adam(1), Trace_ELBO(), A=d.A.values, D=d.D.values)\n",
"svi_result = svi.run(random.PRNGKey(0), 1000)\n",
"p5_1 = svi_result.params\n",
"\n",
"\n",
"def model(M, D=None):\n",
" a = numpyro.sample(\"a\", dist.Normal(0, 0.2))\n",
" bM = numpyro.sample(\"bM\", dist.Normal(0, 0.5))\n",
" sigma = numpyro.sample(\"sigma\", dist.Exponential(1))\n",
" mu = a + bM * M\n",
" numpyro.sample(\"D\", dist.Normal(mu, sigma), obs=D)\n",
"\n",
"\n",
"m5_2 = AutoLaplaceApproximation(model)\n",
"svi = SVI(model, m5_2, optim.Adam(1), Trace_ELBO(), M=d.M.values, D=d.D.values)\n",
"svi_result = svi.run(random.PRNGKey(0), 1000)\n",
"p5_2 = svi_result.params\n",
"\n",
"\n",
"def model(M, A, D=None):\n",
" a = numpyro.sample(\"a\", dist.Normal(0, 0.2))\n",
" bM = numpyro.sample(\"bM\", dist.Normal(0, 0.5))\n",
" bA = numpyro.sample(\"bA\", dist.Normal(0, 0.5))\n",
" sigma = numpyro.sample(\"sigma\", dist.Exponential(1))\n",
" mu = numpyro.deterministic(\"mu\", a + bM * M + bA * A)\n",
" numpyro.sample(\"D\", dist.Normal(mu, sigma), obs=D)\n",
"\n",
"\n",
"m5_3 = AutoLaplaceApproximation(model)\n",
"svi = SVI(\n",
" model, m5_3, optim.Adam(1), Trace_ELBO(), M=d.M.values, A=d.A.values, D=d.D.values\n",
")\n",
"svi_result = svi.run(random.PRNGKey(0), 1000)\n",
"p5_3 = svi_result.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.33"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"UserWarning: The default method used to estimate the weights for each model,has changed from BB-pseudo-BMA to stacking\n",
"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",
"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",
"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"
]
},
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
rank
\n",
"
waic
\n",
"
p_waic
\n",
"
d_waic
\n",
"
weight
\n",
"
se
\n",
"
dse
\n",
"
warning
\n",
"
waic_scale
\n",
"
\n",
" \n",
" \n",
"
\n",
"
m5.1
\n",
"
0
\n",
"
126.515847
\n",
"
4.108174
\n",
"
0.000000
\n",
"
8.879604e-01
\n",
"
13.580529
\n",
"
0.000000
\n",
"
True
\n",
"
deviance
\n",
"
\n",
"
\n",
"
m5.3
\n",
"
1
\n",
"
128.609601
\n",
"
5.490680
\n",
"
2.093754
\n",
"
5.354823e-15
\n",
"
14.173850
\n",
"
1.047467
\n",
"
True
\n",
"
deviance
\n",
"
\n",
"
\n",
"
m5.2
\n",
"
2
\n",
"
139.775924
\n",
"
3.282653
\n",
"
13.260077
\n",
"
1.120396e-01
\n",
"
10.458673
\n",
"
9.845431
\n",
"
True
\n",
"
deviance
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" rank waic p_waic d_waic weight se \\\n",
"m5.1 0 126.515847 4.108174 0.000000 8.879604e-01 13.580529 \n",
"m5.3 1 128.609601 5.490680 2.093754 5.354823e-15 14.173850 \n",
"m5.2 2 139.775924 3.282653 13.260077 1.120396e-01 10.458673 \n",
"\n",
" dse warning waic_scale \n",
"m5.1 0.000000 True deviance \n",
"m5.3 1.047467 True deviance \n",
"m5.2 9.845431 True deviance "
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"post = m5_1.sample_posterior(random.PRNGKey(24071847), p5_1, sample_shape=(1000,))\n",
"logprob = log_likelihood(m5_1.model, post, A=d.A.values, D=d.D.values)[\"D\"]\n",
"az5_1 = az.from_dict(\n",
" posterior={k: v[None, ...] for k, v in post.items()},\n",
" log_likelihood={\"D\": logprob[None, ...]},\n",
")\n",
"post = m5_2.sample_posterior(random.PRNGKey(24071847), p5_2, sample_shape=(1000,))\n",
"logprob = log_likelihood(m5_2.model, post, M=d.M.values, D=d.D.values)[\"D\"]\n",
"az5_2 = az.from_dict(\n",
" posterior={k: v[None, ...] for k, v in post.items()},\n",
" log_likelihood={\"D\": logprob[None, ...]},\n",
")\n",
"post = m5_3.sample_posterior(random.PRNGKey(24071847), p5_3, sample_shape=(1000,))\n",
"logprob = log_likelihood(m5_3.model, post, A=d.A.values, M=d.M.values, D=d.D.values)[\n",
" \"D\"\n",
"]\n",
"az5_3 = az.from_dict(\n",
" posterior={k: v[None, ...] for k, v in post.items()},\n",
" log_likelihood={\"D\": logprob[None, ...]},\n",
")\n",
"az.compare({\"m5.1\": az5_1, \"m5.2\": az5_2, \"m5.3\": az5_3}, ic=\"waic\", scale=\"deviance\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Code 7.34"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.\n",
"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"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAHrCAYAAAAe4lGYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABXUklEQVR4nO3deVxU9f7H8fcBRCQMsVxpUUM0MxdEXHNJuy0uqalZbi1WZtrNm5V27brkzfYszF+2aZLmkppepTSz3BWFzCX3XdQ0wQVJ2c7vj2lGEZAB5gwMvJ6PRw/0nDNnvuczY7znO9/z/RqmaZoCAAAA4HJehd0AAAAAoLgibAMAAAAWIWwDAAAAFiFsAwAAABYhbAMAAAAWIWwDAAAAFiFsAwAAABYhbAMAAAAW8SnsBhRliYmJhd2EPAsMDNTZs2cLuxklCjV3P2ruXtTb/ai5+1Fz9ysONQ8KCsr1GHq2ixkvL15Sd6Pm7kfN3Yt6ux81dz9q7n4lpeYl4yoBAACAQkDYBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAsQtgGAAAALOJT2A0AAAAA8iM93dRvW6TTCdIN5aX69SRvb6Owm5UJYRsAAAAeZ8VKUxMnmTp+4vK2KpWlwYOk1q2KTuBmGAkAAAA8yoqVpkaOMlWjhvTJx4aWRhv65GNDNWpII0eZWrHSLOwmOhC2AQAA4DHS02092s2bSePHGap7hyF/f9vP8eMMNW8mTfw/U+npRSNwE7YBAADgMX7bIh0/IfXtbcjLK/NwES8vQ317Gzp+3HZcUUDYBgAAgMc4nWD7WaN69vvt2+3HFTbCNgAAADzGDeVtP/cfyH6/fbv9uMJG2AYAAIDHqF/PNutI1HRTGRmZx2VnZJiKmm6qShXbcUUBYRsAAAAew9vb0OBBhtauk0aMNLVtu6nkZNvPESNNrV0nDX7WKDLzbTPPNgAAADxK61aGxo2RJk4yNfC5y73bVapI48YYRWqebcI2AAAAPE7rVoZathArSAIAAABW8PY2FNawsFtxbYzZBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAs4lPYDcjOH3/8oe+//14rV67U/v379eeffyowMFBhYWEaMGCA6tev7/S5MjIyNGPGDM2aNUuHDh2Sv7+/mjRpoqFDh6patWrWXQQAAABKvCLZsx0VFaXx48fryJEjat68uR5//HE1atRIP/30k3r16qXo6GinzzVq1Ci9/vrrysjIUJ8+fdS6dWstX75c3bt31969ey28CgAAAJR0RbJnu169epo+fbrCw8Mzbd+0aZMee+wxjRkzRu3bt5evr+81z7N+/XrNnj1b4eHhmjJliuP4Ll266PHHH9fo0aP19ddfW3YdAAAAKNmKZM/2P/7xjyxBW5LCw8PVpEkTnTlzRrt27cr1PHPmzJEkvfDCC5mCebNmzdSyZUtt3LhRBw4ccF3DAQAAgCsUybB9LT4+Ppl+XsuGDRvk7++vsLCwLPtatmwpSdq4caNrGwgAAAD8zaPC9rFjx7R27VpVqFBBoaGh1zw2OTlZp06d0k033SRvb+8s++03Rx48eNCClgIAAABFdMx2dlJTU/Xyyy8rJSVFw4YNyzZAX+n8+fOSpICAgGz327cnJSXleI7AwEB5eXnU5xFJUlBQUGE3ocSh5u5Hzd2LersfNXc/au5+JaHmHhG2MzIy9Oqrr2rjxo3q2bOnunTp4pbnPXv2rFuex5WCgoKUmJhY2M0oUai5+1Fz96Le7kfN3Y+au19xqLkzHxaKfLetaZoaOXKkFi5cqM6dO2vMmDFOPa5s2bKScu65tm/PqecbAAAAKKgiHbbtPdpz585Vx44d9eabbzo9rMPf318VKlTQ0aNHlZ6enmW/faw2C9sAAADAKkU2bGdkZOjf//635s2bpwceeEBvv/12ruO0rxYREaHk5GTFxcVl2bd69WpJUuPGjV3SXgAAAOBqRTJsXxm077vvPr3zzjvXDNoJCQnat2+fEhISMm3v2bOnJGnChAlKSUlxbF+3bp1Wr16txo0bq3r16tZcBAAAAEq8InmD5Mcff6x58+bJ399f1apV0//93/9lOaZ9+/a6/fbbJUnTp0/XxIkTNXjwYA0ZMsRxTNOmTdWjRw/NmTNHXbt2VevWrXX69GlFR0crICBAo0ePdtclAQAAoAQqkmE7Pj5ekm2u7E8++STbY4KDgx1h+1rGjh2rWrVqadasWYqKipK/v7/atm2roUOH0qsNAAAASxmmaZqF3YiiyhOnoykO0+h4GmruftTcvai3+1Fz96Pm7lccal4spv4DAAAAPBVhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwiE9hNyAnCxYsUGxsrLZt26bdu3crNTVV48ePV7du3Zw+x4YNG9SvX78c98+aNUsNGjRwQWsBAACArIps2P7www8VHx+voKAgVaxYUfHx8fk+V0REhCIiIrJsr1y5ckGaCAAAAFxTkQ3b48aN06233qrg4GB9+umneu+99/J9roiICA0ZMsSFrQMAAAByV2TDdvPmzQu7CQAAAECBFNmw7UoHDx7UtGnTdPHiRVWtWlXNmzdX+fLlC7tZAAAAKOZKRNhetGiRFi1a5Pi7n5+fhgwZogEDBhRiqwAAAFDcFeuwXb58eb388stq06aNqlatqnPnzmnDhg1699139c477yggIEC9evXK8fGBgYHy8vK82RGDgoIKuwklDjV3P2ruXtTb/ai5+1Fz9ysJNS/WYbtmzZqqWbOm4+9lypRR586dVbt2bXXr1k2RkZHq2bNnjoH67Nmz7mqqywQFBSkxMbGwm1GiUHP3o+buRb3dj5q7HzV3v+JQc2c+LHhet60LhIaGqn79+vrzzz916NChwm4OAAAAiqkSGbaly59ELl68WMgtAQAAQHFVIsN2Wlqafv/9dxmGoSpVqhR2cwAAAFBMFYuwnZCQoH379ikhISHT9l9//VWmaWbalpaWprffflvx8fFq2bKlypUr58aWAgAAoCQpsjdIzpkzR7GxsZKk3bt3O7bFxMRIktq3b6/27dtLkqZPn66JEydq8ODBmVaKfPHFFyVJDRs2VKVKlXT+/Hlt3LhRBw4cUNWqVTVmzBh3XhIAAABKmCIbtmNjYzV//vxM2+Li4hQXFydJCg4OdoTtnPTq1UurVq1STEyMEhMT5ePjo1tuuUUDBw7UE088ocDAQMvaDwAAABjm1eMs4OCJ09EUh2l0PA01dz9q7l7U2/2ouftRc/crDjVn6j8AAACgEOUrbNuHcgAAAADIWb7C9qOPPqqOHTvqq6++0pkzZ1zcJAAAAKB4yFfYbty4sfbt26c333xTrVq10osvvqgNGza4um0AAACAR8tX2I6KitKSJUv05JNP6vrrr9fixYv12GOP6d5779Vnn32mP//809XtBAAAADxOgWcjSU9P1/LlyzV79mytWbNGGRkZ8vHx0d13363u3bvrrrvukmEYrmqvW3niHbLF4c5eT0PN3Y+auxf1dj9q7n7U3P2KQ82dmY2kwPNse3t765577tE999yjEydO6Ntvv9W8efP0448/6scff1TlypXVvXt39ezZUxUqVCjo0wEAAAAew6VT/1WuXFlhYWGqW7euTNOUaZo6fvy4IiMj1a5dO7311ltKTU115VMCAAAARZZLVpA8efKk5s6dq7lz5yo+Pl6maap27dp6+OGH1axZM/3www/6+uuvNXXqVKWnp+vVV191xdMCAAAARVq+x2xnZGTol19+0Zw5c7Rq1SqlpaXJz89P999/v3r16qX69etnOj4pKUldunTRX3/9pTVr1rik8VbzxHFExWH8k6eh5u5Hzd2LersfNXc/au5+xaHmlo3ZnjBhgubNm6dTp07JNE2FhISoZ8+e6tq1q8qWLZvtYwICAhQeHq4FCxbk5ykBAAAAj5OvsP3JJ5/I19dXHTp0UK9evRQeHu7U4yIiIjx2ZhIAAAAgr/IVtl966SV169bNqa7zK3Xr1k3dunXLz1MCAAAAHidfs5HUq1fPqWXaDx48qI0bN+bnKQAAAACPl6+w3a9fP3322We5Hvf555+rX79++XkKAAAAwOPlK2zb59B25jgAAACgpHLpojZXO3nypPz9/a18CgAAAKDIcvoGye+++y7T3w8fPpxlm11aWpoOHDigdevWZZlvGwAAACgpnA7bw4cPd0zbZxiG4uLiFBcXl+PxpmmqdOnSeu655wreSgAAAMADOR22n3vuORmGIdM09fHHH+v2229Xu3btsj22VKlSqlixolq0aKGKFSu6rLEAAACAJ3E6bA8ZMsTx5/nz56tZs2YaPHiwJY0CAAAAioN8LWqzfPlyV7cDAAAAKHYsnY0EAAAAKMmc6tkeMWJEvp/AMAy98cYb+X48AAAA4KmcCtvz58/P9xMQtgEAAFBSORW2p02bZnU7AAAAgGLHqbAdERFhdTsAAACAYocbJAEAAACLELYBAAAAi+Rrnm3Jthz7woUL9dNPP+nQoUO6cOGCTNPMcpxhGFq2bFmBGgkAAAB4onyF7ZSUFD3zzDNav359tgFbkmNpdwAAAKCkytcwkilTpmjdunVq06aNli5dqgcffFCGYWjr1q2Kjo7W4MGDVaZMGT355JPauXOnq9sMAAAAeIR89WxHR0crMDBQ7733nvz9/eXlZcvspUqVUo0aNTR48GA1bdpU/fr1U/Xq1dW9e3eXNhoAAADwBPnq2T58+LDq1asnf39/SbYhI5KUnp7uOCY8PFxhYWGaMWOGC5oJAAAAeJ58hW0vLy8FBAQ4/m4P3QkJCZmOq1Spkg4cOFCA5gEAAACeK19hu1KlSjp+/Ljj77fccosk6bfffst03K5du3TdddcVoHkAAACA58pX2K5fv7727NmjixcvSpJat24tSfrvf/+rFStWaNeuXXr99de1b98+1atXz3WtBQAAADxIvsL2vffeqzJlymjNmjWSpFtvvVX9+/fX8ePHNXDgQHXp0kXTp0+Xn5+fhg0b5tIGAwAAAJ4iX7ORtGnTRqtXr860bfjw4brzzju1bNkynTt3TtWqVVPfvn1VrVo1V7QTAAAA8Dj5XkEyOx06dFCHDh1ceUoAAADAY+VrGAkAAACA3BW4Z/vYsWM6deqUUlJScjymcePGBX0aAAAAwOPkO2x/++23mjRpUqYpAHOyY8eO/D4NAAAA4LHyFbbnzp2rkSNHSpJCQ0NVrVo15tMGAAAArpKvsD116lT5+Pjoo48+0t133+3qNgEAAADFQr5ukDx48KDCw8MJ2gAAAMA15CtsBwYGyt/f39VtAQAAAIqVfIXtdu3aacuWLY7l2gEAAABkla+w/a9//UsBAQEaPny4zp075+o2AQAAAMVCvm6QfPPNNxUSEqIlS5ZozZo1qlu3ripXrpztsYZh6I033ihQIwEAAABPlK+wPX/+fMefz58/r3Xr1uV4LGEbAAAAJVW+wva0adNc3Q4AAACg2MlX2I6IiHB1OwAAAIBiJ183SAIAAADIXb56tu0SEhK0cOFCbd26VWfOnFHTpk311FNPSZJ2796tI0eOqHnz5ipTpoxLGgsAAAB4knyH7cWLF+u1117TX3/9JdM0ZRiGKlas6Nh/6NAhPf/88xo/fry6dOniirYCAAAAHiVfw0g2bdqkl156Sb6+vhoxYoS+/fZbmaaZ6Zg2bdqobNmy+vHHH13SUAAAAMDT5Ktne/LkyfLx8dHUqVNVu3btbI8pVaqUatSoob179xaogQAAAICnylfP9m+//ab69evnGLTtKleurJMnT+arYQAAAICny1fYvnjxooKCgnI9LikpSYZh5OcpAAAAAI+Xr7BdtWpV7dq165rHpKWladeuXbr11lvz1TAAAADA0+UrbLdt21aHDx/W9OnTczxmypQp+vPPP9W+fft8Nw4AAADwZPm6QfKpp57S4sWLNW7cOG3evFnt2rWTZJt3++eff9ayZcs0f/58ValSRf369XNpgwEAAABPYZhXz9nnpH379umFF17Qnj17ZBiGY65tSTJNUzVq1NDEiRNVo0YNlzbYnRITEwu7CXkWFBTkke32ZNTc/ai5e1Fv96Pm7kfN3a841NyZexjzvajNbbfdpgULFmj58uVau3at4uPjlZ6ersqVK6t58+a699575e3tnd/TAwAAAB6vQMu1e3l5qX379ozLBgAAALKRrxskAQAAAOSuQD3bO3fu1DfffKO4uDjH4jUVK1ZUWFiYHn74YdWpU8cljQQAAAA8Ub7D9scff6xJkyYpPT090/azZ89qz549+vbbbzVw4EANGTKkwI0EAAAAPFG+wvZ3332nyMhI+fv7q3fv3urQoYNuuukmSVJ8fLwWLVqkGTNmaNKkSbr55pvVpUsXV7YZAAAA8Aj5CtvTpk2Tj4+Ppk2bprp162baV6tWLdWqVUv33nuvevXqpWnTphG2AQAAUCLl6wbJffv2qUmTJlmC9pXq1q2rpk2bat++ffluHAAAAODJ8hW2AwICFBgYmOtxZcuWVUBAQH6eAgAAAPB4+Qrbd911l2JiYnTx4sUcj7l48aI2btyoli1b5rtxAAAAgCfLV9h+8cUX5evrq8GDB+vQoUNZ9h86dEhDhgyRr6+vhg0bVuBGAgAAAJ4oXzdIvv/++6pdu7aWL1+u+++/X7fffruCg4Ml2WYj2blzpzIyMtSmTRu9//77mR5rGIbeeOONgrccAAAAKOIM0zTNvD6odu3a+X9Cw9COHTtyPW7BggWKjY3Vtm3btHv3bqWmpmr8+PHq1q1bnp4vIyNDM2bM0KxZs3To0CH5+/urSZMmGjp0qKpVq3bNxyYmJubpuYqCoKAgj2y3J6Pm7kfN3Yt6ux81dz9q7n7FoeZBQUG5HpPvqf+s9uGHHyo+Pl5BQUGqWLGi4uPj83WeUaNGafbs2QoJCVGfPn10+vRpRUdHa82aNZo5c6ZCQkJc3HIAAADAJl9hOyIiwtXtyGLcuHG69dZbFRwcrE8//VTvvfdens+xfv16zZ49W+Hh4ZoyZYp8fX0lSV26dNHjjz+u0aNH6+uvv3Z10wEAAABJ+bxB0h2aN2/uGAeeX3PmzJEkvfDCC46gLUnNmjVTy5YttXHjRh04cKBAzwEAAADkpMiGbVfYsGGD/P39FRYWlmWffUrCjRs3urtZAAAAKCHyNYzEEyQnJ+vUqVMKDQ2Vt7d3lv32myMPHjyY4zkCAwPl5eV5n0ecGawP16Lm7kfN3Yt6ux81dz9q7n4loebFNmyfP39eknJcwdK+PSkpKcdznD171vUNs1hxuLPX01Bz96Pm7kW93Y+aux81d7/iUHNnPix4XrctAAAA4CGKbdguW7aspJx7ru3bc+r5BgAAAAqq2IZtf39/VahQQUePHlV6enqW/fax2rktbAMAAADkl9Nh+8yZMzp27JhjLPS1nD9/XseOHdOZM2cK0rYCi4iIUHJysuLi4rLsW716tSSpcePG7m4WAAAASginwnZycrI6duyozp0769SpU7kef+rUKXXu3FldunTRxYsXC9zI3CQkJGjfvn1KSEjItL1nz56SpAkTJiglJcWxfd26dVq9erUaN26s6tWrW94+AAAAlExOzUby3Xff6c8//9RLL72kGjVq5Hp8jRo1NHjwYL355ptauHChI/TmxZw5cxQbGytJ2r17t2NbTEyMJKl9+/Zq3769JGn69OmaOHGiBg8erCFDhjjO0bRpU/Xo0UNz5sxR165d1bp1a8dy7QEBARo9enSe2wUAAAA4y6mw/dNPP6lMmTLq06eP0yd+5JFH9OGHH2rp0qX5CtuxsbGaP39+pm1xcXGOISHBwcGOsH0tY8eOVa1atTRr1ixFRUXJ399fbdu21dChQ+nVBgAAgKUM0zTN3A5q0aKFatasqalTp+bp5I899pj27t3rGB/taTxx7sfiMGelp6Hm7kfN3Yt6ux81dz9q7n7FoeYum2f77NmzuvHGG/PcgBtvvLHQb5IEAAAACotTYfu666675kqLOUlKSpK/v3+eHwcAAAAUB06F7eDgYG3bti3PJ9+2bZuCg4Pz/DgAAACgOHAqbDdr1kynT5/WwoULnT7xggUL9Oeff6pFixb5bhwAAADgyZwK23379pWvr6/GjBmjjRs35np8TEyMxowZo9KlS+dpBhMAAACgOHEqbFeuXFkjRozQhQsX1L9/f7344ov6+eef9ccffyg1NVWpqan6448/9PPPP+vFF1/UY489puTkZA0fPlyVK1e2+hoAAACAIsmpebYlqVevXvL29ta4ceO0ePFiRUdHZ3ucaZry9fXVq6++ql69ermsoQAAAICncTpsS1KPHj3UqlUrRUVFacWKFdq7d6/s03QbhqGQkBC1bt1affr0oUcbAAAAJV6ewrYkVapUScOGDdOwYcOUmpqqc+fOyTRNBQYGqlSpUla0EQAAAPBIeQ7bVypVqpRuuOEGV7UFAAAAKFacukESAAAAQN451bPdr1+/fD+BYRj66quv8v14AAAAwFM5FbZjYmLy/QSGYeT7sQAAAIAncypsT5s2zep2AAAAAMWOU2E7IiLC6nYAAAAAxY6lN0gmJSVpzpw5Vj4FAAAAUGQVaOq/7KSlpWnlypVasGCBfvnlF6WkpKhHjx6ufhoAAACgyHNZ2N68ebMWLlyo6OhonT17VqZpytvbW02bNnXVUwAAAAAepUBh+/Dhw1q4cKH+97//6fDhw5Ik0zTVsGFDdejQQffffz+L3gAAAKDEynPYPnPmjKKjo7Vw4UL99ttvkmwBu0aNGjp79qwSEhL0zTffuLyhAAAAgKdxKmynpKRo+fLlWrhwoVatWqW0tDSZpqkbb7xR999/vx588EHVrVtXjz76qBISEqxuMwAAAOARnArbLVq0UFJSkkzTVJkyZXTvvfeqc+fOatmypby8WPEdAAAAyI5TYfv8+fMyDEMVK1bU66+/rtatW1vdLgDIl/R0U79tkU4nSDeUl+rXk7y9WckWAFA4nArbLVu21Lp163Ty5EkNHDhQVapUUadOndSpUyeFhIRY3UYAcMqKlaYmTjJ1/MTlbVUqS4MHSa1bEbgBAO7n1BiQzz//XCtXrtTw4cNVu3ZtHTt2TJMnT1anTp3UtWtXTZ06VadOnbK6rQCQoxUrTY0cZapGDemTjw0tjTb0yceGatSQRo4ytWKlWdhNBACUQIZpmnn+DbRv3z7Nnz9fixcv1vHjx2UYhry8vOTj46OUlBRt3LhRAQEBVrTXrRITEwu7CXkWFBTkke32ZNTc/a6ueXq6qV69bUF7/DhDXl6Xe7EzMkyNGGlq/wFp5tcGQ0rygfe4+1Fz96Pm7lccah4UFJTrMfm6u/G2227TsGHDtHz5ck2dOlVdunSRn5+fLl26JNM01bx5cw0ePFjff/+9Ll68mJ+nAACn/bZFOn5C6ts7c9CWJC8vQ317Gzp+3HYcAADuVKBFbQzDUNOmTdW0aVONHj1ay5Yt04IFC7R27VotW7ZMP/30k8qUKaO4uDhXtRcAsjj994yjNapnv9++/TQzkwIA3Mxly7WXLl1aHTp0UIcOHZSQkKCFCxdq4cKF+v333131FACQrRvK237uPyDVvSPr/v0HMh8HAIC7ODWMpHv37nrrrbe0bNkypxatKV++vB577DHNmzdPixcvLnAjAeBa6tezzToSNd1URkbm21AyMkxFTTdVpYrtOAAA3Mmpnu1t27Zp+/btmjp1qiSpevXqCg8PV6NGjRQeHq7g4OAcH3vbbbe5pKEAkBNvb0ODB9lmHRkx0lTf3rahI/sP2AL42nXSuDHcHAkAcD+nZiNZsWKFYmNjFRsbq23btunSpUu2Bxu2X1yVKlVyBO/w8HDVrFnT2la7iSfeIVsc7uz1NNTc/XKqebbzbFeRBj9rMM92AfAedz9q7n7U3P2KQ82dmY0kz1P/paamauvWrY7w/euvv+rs2bO2k/0dvq+//nqFhYU5er8bNGiQ99YXAZ74BigOb1xPQ83d71o1ZwVJ1+M97n7U3P2oufsVh5pbErazs3fvXsXGxmrTpk2Ki4tTfHy8I3gbhuGxN0l64hugOLxxPQ01dz9q7l7U2/2ouftRc/crDjV3Jmy7ZDaSkJAQhYSE6OGHH9aePXsUHR2t6dOn69y5c644PQAAAOCRChS209LStG3btky92ufOnZNpmjIMQyEhIQoLC3NVWwEAAACPkqewfeHCBW3evFmbNm3Spk2btHXrVseqkX5+frrzzjsVFhamsLAwNWzYUNdff71V7QYAAACKPKfC9htvvKFNmzZp165dSk9PlyTdeOONatWqlSNc16lTRz4+LlsjBwAAAPB4TqXjadOmyTAM1a5dW71791aTJk108803W902AAAAwKM5Fbb9/Px08eJF7dixQ//973915513qlGjRo7hIgEBAVa3EwAAAPA4ToXt2NhYbd++PdPc2jExMTIMQ4ZhqGbNmgoLC3ME8KpVq1rdbgAAAKDIy/c82wcOHHCE77i4OB06dMh2QsNwrChpD+C1a9d2aaPdxRPnfiwOc1Z6GmruftTcvai3+1Fz96Pm7lccam7pPNvVq1dX9erV1b17d0lSQkKCNm3a5AjgS5YsUXR0tEcvagMAAAAUhJerThQQEKCgoCAFBQUpMDBQvr6+Mk1TLligEgAAAPBI+e7ZPnfunOLi4hy92du2bVNaWpokOQL2LbfcokaNGrmmpQAAAICHcTpsnzhxwhGsN23apH379mXqufby8lKtWrUUHh6uRo0aKTw8XBUqVLCs4QAAAEBR51TYvvvuu3X8+HFJl3utfX19VbduXUe4btSoEVMAAgAAAFdwKmwfO3ZMAQEBatCggcLDwxUeHq569erJ19fX6vYBAAAAHsupsD1v3jzVrl1bXl4uu58SAAAAKPacCtt16tSxuh0AAABAsUNXNQAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARwjYAAABgEcI2AAAAYBHCNgAAAGARn8JuwLVs2bJFkZGR2rx5s1JTUxUSEqL+/furU6dOTj1+w4YN6tevX477Z82apQYNGriotQAAAEBmRTZsb9iwQU8++aRKlSqlDh06qGzZslq6dKmGDRum+Ph4DRw40OlzRUREKCIiIsv2ypUru7LJAAAAQCZFMmynpaVp5MiRMgxD06dPV506dSRJzz33nHr16qXIyEjdd999qlatmlPni4iI0JAhQyxsMQAAAJBVkRyzvX79eh0+fFgdO3Z0BG1JCggI0KBBg5SWlqZ58+YVYgsBAACA3BXJnu2YmBhJUsuWLbPsa9GiRaZjnHHw4EFNmzZNFy9eVNWqVdW8eXOVL1/eNY0FAAAAclAkw/bBgwclSbfeemuWfYGBgQoKCtKhQ4ecPt+iRYu0aNEix9/9/Pw0ZMgQDRgwoMBtBQAAAHJSJMN2UlKSJKls2bLZ7g8ICNCJEydyPU/58uX18ssvq02bNqpatarOnTunDRs26N1339U777yjgIAA9erVK8fHBwYGysurSI60uaagoKDCbkKJQ83dj5q7F/V2P2ruftTc/UpCzYtk2HaVmjVrqmbNmo6/lylTRp07d1bt2rXVrVs3RUZGqmfPnjkG6rNnz7qrqS4TFBSkxMTEwm5GiULN3Y+auxf1dj9q7n7U3P2KQ82d+bBQJLttAwICJEnnz5/Pdn9SUlKOvd7OCA0NVf369fXnn3/maTgKAAAAkBdFMmzbp/TLLgifPXtWiYmJ2Y7nzgv7J5GLFy8W6DwAAABATopk2G7cuLEkafXq1Vn2rVmzRpKyXaTGWWlpafr9999lGIaqVKmS7/MAAAAA11Ikw3azZs108803a9GiRdqxY4dje1JSkiZNmiQfHx917drVsT0hIUH79u1TQkJCpvP8+uuvMk0z07a0tDS9/fbbio+PV8uWLVWuXDlLrwUAAAAlV5G8QdLHx0fjxo3TgAED9Oijj6pjx44KCAjQ0qVLdfToUb3wwguqXr264/jp06dr4sSJGjx4cKaVIl988UVJUsOGDVWpUiWdP39eGzdu1IEDB1S1alWNGTPG7dcGAACAkqNIhm1Jatq0qWbMmKGPPvpI33//vVJTUxUSEqJ//vOf6ty5s1Pn6NWrl1atWqWYmBglJibKx8dHt9xyiwYOHKgnnnhCgYGBFl8FAAAASjLDvHqcBRw8cTqa4jCNjqeh5u5Hzd2LersfNXc/au5+xaHmHjv1HwAAAFAcELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACL+BR2AwAAKIj0dFO/bZFOJ0g3lJfq15O8vY3CbhYASCJsAwA82IqVpiZOMnX8xOVtVSpLgwdJrVsRuAEUPoaRAAA80oqVpkaOMlWjhvTJx4aWRhv65GNDNWpII0eZWrHSLOwmAgBhGwDgedLTbT3azZtJ48cZqnuHIX9/28/x4ww1byZN/D9T6ekEbgCFi7ANAPA4v22Rjp+Q+vY25OWVebiIl5ehvr0NHT9uOw4AChNhGwDgcU4n2H7WqJ79fvt2+3EAUFgI2wAAj3NDedvP/Qey32/fbj8OAAoLYRsA4HHq17PNOhI13VRGRuZx2RkZpqKmm6pSxXYcABQmwjYAwON4exsaPMjQ2nXSiJGmtm03lZxs+zlipKm166TBzxrMtw2g0DHPNgDAI7VuZWjcGGniJFMDn7vcu12lijRujME82wCKBMI2gBKLlQc9X+tWhlq2EK8jgCKLsA2gRCqslQcJ+K7n7W0orGFhtwIAskfYBlDi2FcebN5MGvWaoRrVbbNXRE23bR83xprAXRKXFufDBYCSjrANoES5euVB+4Iode+Qxo+z3Ww38f9MtWzh2lBYWAG/MHnahws+GACwArORAChRCmPlwZK4tLj9w0WNGtInHxtaGm3ok48N1aghjRxlasXKgl1rerqpuF9N/fiT7WdBa7dipalevU09P9TUmNdtP3v1Lng7AYCwDaBEKYyVB0va0uJWf7hwdTC2+oMBgJKNsA2gRCmMlQdL2tLiVn64cHUwLonfOgBwL8I2gBKlMFYeLGlLi1v14cKKYFzSvnUA4H6EbQAlSmGsPFjSlha36sOFFcG4pH3rAMD9CNsAShzbyoOG9u+XBj5n6h8P2FYg3H/AmpUHS9rS4lZ9uLAiGJe0bx0AuB9T/wEokdy98mBJWlrc9uHCNoZ6xEhTfXsr01SHa9fZrjmvtb4yGNe9I+v+/ATjKz8YjB+nTD3mxfFbBwDuR9gGUGK5e+XBkrS0uBUfLqwIxlZ9MAAAO8M0TW6xzkFiYmJhNyHPgoKCPLLdnoyaux81d6+C1NvVC8VcuThQ395GtsE4P0E+2wV4qtiG9xTGtw68x92Pmrtfcah5UFBQrsfQsw0AsIyrvz2wajhOSfrWAYB7EbYBAB7FqmDs7mFFAEoGwjYAwOMQjAF4Cqb+AwAAACxC2AYAAAAsQtgGAAAALELYBgAAACxC2AYAAAAsQtgGAAAALELYBgAAACzCPNsAMnH18toAAJRkhG0ADitWmpo4ydTxE5e3VaksDR6kfC+DDQBAScYwEgCSbEF75ChTNWpIn3xsaGm0oU8+NlSjhjRylKkVK83CbiIskJ5uKu5XUz/+ZPuZns7rDACuRM82AKWn23q0mzeTxo8z5OVl68Wue4c0fpw0YqSpif9nqmWLkjWkpLgPqeGbDACwHj3bAPTbFun4Calv78tB287Ly1Df3oaOH7cdV1KsWGmqV29Tzw81NeZ1289evYtPDz/fZACAexC2Aeh0gu1njerZ77dvtx/nafI6VKK4B9Grv8moe4chf3/bz/HjDDVvJk38P4aUAIArMIwEgG4ob/u5/4Bt6MjV9h/IfJynSE83Ne1rU3PnSWfOXt5+raESJWFIjf2bjFGv5fRNhjTwOdsQmrCGhdRIACgm6NkGoPr1pMqVpImTTC39MSNT729Ghqmo6aaqVLEdV9Tk1Gu9YqWpLg+Z+mLK5aBdqZL0zFO6Zg91SRhSU9y/yQCAooSebQBavUa6dEnatt32n2SqUiWpS2dT27ZLa9dJ48YYhdaTe/WNim1aXw7U2d3g17aNqW9mSb6+Ur07pbfflA4eNBQ13dSnn0tjR9mOza6HuiQE0eL6TQYAFEWEbaCEs49Pbt7MFry+WyD9cVL64w9p8mdSuXK2oF1Ys1NkF6iDq55R61YZ+maW1LyZbThEjeq2kDjta1MzZko33ywdOSINGmgo4Doj0zCQSZOl116VBg1RlqESJSGI1q9n+1ASNd3U+HHK1INf1L/JAABPwzASoAS7enxy395emv2NoY8+MPSfkbawWaaM1LJF4czHnNONiiEhXpoxU6oVmvUGv57dbY89k2j7eWUP9ZXDQC5csG27uof6yiCakZH5GotLEPX2NjR4kKG162wfPrZtN5WcbPs5YqSpteukwc8W3jcZAFCc0LMNlAA5zRed3Y1y3t7G3z29hqpWMTXwOdtNhtHfy63zMV/rRsX+fctoxcrzOvWnZF6V+RPPGJJMnU+y/f3qHmp7+N65y/bz6h5qWxC1jekeMdJU395y9JpHTTcLfUiNq7RuZWjcGNs4/YHPXS5ilSqF+00GABQ3hG14jOK+wIhVrrVwSUqq7e+5jU/+YorUovnl4Rp795ma9In07/+YevJxU/36uD58XmvGjD9P28Lh6dO24+rXu/zeSEi4HBzLBWYdKmEfBrIhRjn2UJeUINq6laGWLcS/KwCwEGEbHoGV7vLnyvHYV45rjppu2/7EY7bjchqfvGevLWjaxjvbQq/ttbjcy/3FFGnx96aGuPi1uPpGxSs/bF38K8Nx3Oo1psa/lbnX3ctLysiQunWRpkyTo4e62q22DwmlS9tuBP3v2Jw/JJSUIHr5mwwAgBUM07z6S1jYJSYmFnYT8iwoKKhQ2m1lr/OVgbFvbyPbr/QLM3AXVs1zk55uW/GwRo3MwzAk29jjESNN7dtv+/ttORwzaIhtNpJJkVK9O72yvBYXL5p64UVbGN/+u2tfi7hfbas2fvKxodOnleXDlj1QG4YyvTf27Tf10itS0gXbmO42raTvFtpu+rQrV0566V/Fp4faakX1PV6cUXP3o+buVxxqHhQUlOsx9GwXU9cKv64Oxlb2OluxwEhhDkdx53M7u3DJk49LX06VBg0xFRFhqnaodN110oyZ9mkApZDbjGxfi+RkSTLVrasUGJj9VHr5Zb9RccJHpnbtvjzryK23mlq92l/vfpCsS5ckf3+p9yPKNLwl6YJ04w3Snr3Srt2Xz1munPRQV1ky7AUAgOwQtosQVwWxa4VfKWsPYZXK0qCBpgIDjUzPLeX+FXpuwxTGjXE+cGd3/dkFxvR0U7FxppYslZLOS8ePS59+YapJY9tj0tNNzV8gxR+TgqtKnTtKO3YaOvWnqU2bpJWrpAvJl5+3XKD0ULe8jTtOTzf162ZTv26W/PySdXvtDDVsYFyzZtm9LuXLS0OfN9W2jesnBnJ2vujz521B+fIc2zZBQdKTj9uGiew/IF28aGv7A/dLm3+zjZO2j3++8QbXrzro7W1o0EBTr422taX3I9KxY6b+/R/p9OnLL+CFC7Yp/KTMX9L9edr2MyJCuvceqcKNRrEcBgIAKNqK9DCSLVu2KDIyUps3b1ZqaqpCQkLUv39/derUyelzZGRkaMaMGZo1a5YOHTokf39/NWnSREOHDlW1atWu+Vh3frWRc0DO21fdm2L9NXRYUrZDLtaszfqV+/4Dtp5D+8wMdkHlbNHlzJmc2+PMMIX9B6SZX+ceYnO6/rtaSrO/lZZG26Z2W7HS1BtvmY5p267m7y9dvGgbXpBXQeWkYU4MLVix0tQ775uZaiPZeoNLlcq+ZpJtdotaodKpU1mnm3u0lzRooGsD95XDMOrekfWatm233fxnf0/0fsQWXHfust08uG279PpoadIntjB+/Lh09tzlx1euZOspPnvO9hpfuiT94wFTo14zdE871w4luaF85ppVqGDorham5n2X+fi6d0jNm0oL/sewEVcqDl/1ehpq7n7U3P2KQ82dGUZSZOfZ3rBhgx599FFt2rRJ9957rx555BElJiZq2LBh+uSTT5w+z6hRo/T6668rIyNDffr0UevWrbV8+XJ1795de/futfAKnJfTXMLXWlI6O+nppt55L9nxNf+Vcw+PG2O7KczXVxo3Ro59p0/bvmYPCpJuuEH6fpFtOevEM7bQ+MxTyrE9rlrW+lrXP/tb2zH7D9iO+/d/bEG7jJ9t+623Zj5XcrItaHfpLA1/ybbN2/vyfn9/6aZg25/LlLm8/dZbpOS/cq+3vQ1nzthWJvzwfenLz8rq1ltsQTW7mv37P6befd8WtHfukmrXvnydkyJttZ8xU/r5l3x8QriG3OaLnva1KW9vqVlT2/ul3p1eatbUS4/399KkSEMtmkuTJkutW9nanf538956Q3r5RelSim1729a23mIrFnuxB+xpU2zfAoTWlCa8J/20JEhtWl9+z/n42N7fPR6SPvtSCgmxtVGy9eD/9ZftdXD23xIAAK5SJMN2WlqaRo4cKcMwNH36dI0bN06vvPKKFixYoJo1ayoyMlIHDx7M9Tzr16/X7NmzFR4ervnz5+vll1/WW2+9pU8//VRJSUkaPXq05deSm6vHwV4ZkMePM9S8mW0crDMLiPy2RYo/lpFt+N26zdbzeOmS7c9XP/d/x9qmUdu1S1r4P1tPZ/Nm0sJFthCTXXtcsay1M9fv7S19FWUq8mNbDcqVk64PtLUvuKpUubLUJOLyOZtESOtjpKlRtmNuvOHyvjdel9LSbQH3+uttvaClS0tJSbba3FEn53qnp9vaULq07bwTPzTUKMxL4Y1K6dIlW7tKl85as7p32D68nDxpmz7vyuusd6eX/jvWdv4PPpJLF4pxZuGS9HTb+OVrfVha+qNUu5at516SXnlVevs9yc/Ptv3nFVJqaoYli73Yg/uy5VJCgvSvFwyFN/KSt7dtSIh9f1qa7fX7cKIc7+c162xT9b041Lav7h3O/1sCAMBVimTYXr9+vQ4fPqyOHTuqTp06ju0BAQEaNGiQ0tLSNG/evFzPM2fOHEnSCy+8IF9fX8f2Zs2aqWXLltq4caMOHDjg+gvIA1f1DkvXHqN7ZeC1//nK5w65zfbcv262bevXx1C/Ppmf++r2XLmsdXac6enM7fr79TGUni6tWy+d+MO2vV1b21LiSUm27UMGGbr55suPu/lm6cQJ239NIjIPJ/hpuW37U0/aztGkiS2I2WsSEaEc6/3bFlsbLl3KHFBj49J04g/p6Sdt+66umf2DQEJi9tdpr31CgnOvc17Y5os2tH+/bTz1Px6wDR3Zf0COlRZz+7CUkCi98Lyh2TMMPfm4bVvdO6SRI6SBT9uud8gLsmTVQXvv/OLorG01DKlChczHJyRILZpJ//7P5fbUDLG1p8k1XlsAAKxSJMN2TEyMJKlly5ZZ9rVo0SLTMdeyYcMG+fv7KywsLMs++7k3btxYkKYWmCt6h+2uFX6vDLz2P1/53Fc/pkb17J/7ym2uWNba2etvekXP9dz5tp8nT12eau7Spcv7U674s1/pzOfb+/dUdy2a2X6Wvmp/7dDM7cqurVe399SftvEVzZtlf2ztWtk/zu7K2jvzOudV61aGZk63LcE+6jXbz5lfG2rZIvMiL9dqV43qtp7yx/t76b9jbcOPnnteeuHvoRpH462ZgtHeO7/n7xFfy36y9c5v/i1VI0baZilp1zbzY95+z9Z2e3vs12F/HayoMQAAOSmSYds+ROTWqwfkSgoMDFRQUJAOHTp0zXMkJyfr1KlTuummm+R95aDdv9lvjnRmOIqVXNE7bFe/nhRc1Svb8HtnXdvwh9KlbX++8pz79l8Oxg0bXH7e7J77ym3ODFPIrafT2etv3PjytqZNbT//8+/LM51cGZp9r/jzxSuCt2Tr+ZRsQwwkZQrp5ctfHiqRXb2zq4MkVbjR9s9o7brsj7WfU7JNTXcl+4cS+1AXV453vpJt4RLbjYthDQ3HMIzcPixl9/pcGd7tPd2jX7NucaHWrQyNHWUbTvT2e7abMHv3O+cI1CNftb3+9n/mQ/9pu2GzdSsj04e+a722AABYpUhO/ZeUlCRJKlu2bLb7AwICdOLEiWz32Z0/f95xbE7nuPK5shMYGCgvL2s/j7RpbSq46hnNnO2tyAlls8zoMXP2ed0UnK42rcs59fX8Sy9e0tBhSfrPGB899WQZ1Qzx0Z69afrsi7906VKqDEMaM66UnnqyjBqHe6vCjWc0cpR0+rSpCe8F6O62vgquekbfzPKSISPTc2fXni4PStddd0nvvJesgc9dvsHvpmAvffCuv+5pX/oarXX++p94LFBz55/VsWMZ2r1LqlrF0OxvfXRXS9t75I8/zklKkySdOOGl4KqmTBmK+9VLVSqn6/gJ8+/jbL3dkz81VLmyqZ9XeMvLK10ZGdKrr/hrztzUHOvdprWpqlUSlZBoauYsH0V+aGtvozBTVasY+nyKKT8/6YYbvDLVbPa351W+fJoSEky9NtrQhPeuU63QUo7XZe26VNW53Vt+fqbTr7OrvPJSzu+XtetS9d7bAXrvg+RsX5+2bUzN+872+tzd1tp2d+sqXXfdRf3rpQu6/XZvdetSWh0fKK39B9I1Ztzl97ZfaemXFd56uMf12rc/3XEd778ToNnfXsrTvyVk5swd93Atau5+1Nz9SkLNi2TYLirOnj3rlucZNNDUyFGpGvhcQo4rJJ47d8apc93TPkjjxlzQxEmp6t0v1bG9ShXb0tSSsuyTbF+x+/ld0IkTF9TxAVOTP7cF52eekk6cSLxme8IbSTOiTP225cp5uk15eycrMTFZuXHm+i9cOKvnBtrmWE5IlMr4mfplRap6PpKgc+el+PjL51u9JkMPdpJq1zL11rsZmWYjqVzJNvzk4iVTZ89JJ06k/103af6C5Fzr/dyztvHAv6xMVe9+CXriMSkgoKxKlTJ17LjtmI4PZGRbs+2/SzNmmur72HnH+W68wbbK4e870vP0OrtKeCNb27J7v4wbYyiicbJL358FEdHY9h6eOCld/x2frP+OT3a01f7efvd9U79uTlfTlrappCpXlp4eIH07L8mtbS1uisP0XJ6GmrsfNXe/4lBzZz4sFMl5tp9//nktWbJEc+fOVd26dbPsb9q0qQzD0Lp167J5tE1ycrIaNmyo0NBQ/e9//8uy/5dfftEzzzyjJ598Ui+//HK25yj0ebar2IZh5OXrefsbNy8rSJ49a1t178rnDgqSTPOqOaPz0R5nOXv9Vs6z7ez15Wme7avO+fMvGfrgQ9sHhrw+r5VyW1DJVe9PV7V1/4GyOnjofLbv7Wlfm5o7TzpzxWflolBjT1YcfiF6GmruftTc/YpDzT12uXb7eOpDhw5lCdtnz55VYmKiGja89hJ1/v7+qlChgo4ePar09PQs47btY7VzW9jGXVq3MtSyRe4rNjrLNkbX2X2GWt2VNWxJrmtPbpy9fvtx9hUkLyRLN95om7KvUkUj1xUkExKkc+ckLy/b+b28pMQzRp6uz96GyytIltHttf/KdQVJSWrbxivbWhf2sIZrvV8k178/C8Lb21BE41KOWUau3vd4f0P9+hS9GgMASqYiGbYbN26syZMna/Xq1erQoUOmfWvWrJEkRUREZPfQTCIiIrR48WLFxcWp8ZV32ElavXq147mKitwCT2E8tzvb4+z128KWoYgcXjpvb0MP98i8zXZe14Utb29D4Y0MhTeSgoL8lZh4+U7L3K6hMF/ngvCkdntSWwEAxVuRnI2kWbNmuvnmm7Vo0SLt2LHDsT0pKUmTJk2Sj4+Punbt6tiekJCgffv2KcE+1cTfevbsKUmaMGGCUlJSHNvXrVun1atXq3HjxqpePYc55wAAAIACKpI92z4+Pho3bpwGDBigRx99VB07dlRAQICWLl2qo0eP6oUXXsgUkqdPn66JEydq8ODBGjJkiGN706ZN1aNHD82ZM0ddu3ZV69atdfr0aUVHRysgIKBIrCAJAACA4qtIhm3JFpRnzJihjz76SN9//71SU1MVEhKif/7zn+rcubPT5xk7dqxq1aqlWbNmKSoqSv7+/mrbtq2GDh1KrzYAAAAsVSRnIykqPPEO2eJwZ6+noebuR83di3q7HzV3P2rufsWh5s7MRlIkx2wDAAAAxQFhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCCtIAgAAABahZxsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsIhPYTcA+XPq1ClNmDBBK1as0NmzZ1W1alV16tRJTz/9tHx9fZ06x8GDB/XDDz9o1apVOnTokM6cOaMbbrhBTZo00TPPPKPbbrvN4qsomrZs2aLIyEht3rxZqampCgkJUf/+/dWpUyenz5GRkaEZM2Zo1qxZOnTokPz9/dWkSRMNHTpU1apVs67xHqig9d60aZOWLVummJgYxcfHKzk5WcHBwWrXrp2eeeYZXX/99RZfgedxxXv8Sqmpqerevbt27typ6tWr64cffnBxiz2fq2qelJSkL7/8UkuXLtWRI0dUqlQp3XzzzWrXrp0GDx5sUes9kytqfu7cOU2ZMkXLli3T0aNH5evrq5tuukldu3ZVjx49VLp0aQuvwHMsWLBAsbGx2rZtm3bv3q3U1FSNHz9e3bp1y9N5iuvvTqb+80CnTp1Sz549dfz4cbVv317VqlVTbGys4uLidNddd+nTTz+Vl1fuX1oMHTpU0dHRCg0NVVhYmAICArR7926tXLlSfn5++uKLLxQeHu6GKyo6NmzYoCeffFKlSpVShw4dVLZsWS1dulRHjx7V0KFDNXDgQKfO89prr2n27NkKCQlR69atdfr0aUVHR6t06dKaOXOmQkJCLL4Sz+CKerdo0UKJiYlq1KiRbr/9dhmGoZiYGP3++++65ZZbNHPmTN1www1uuBrP4Kr3+JU+/PBDTZ06VcnJyYTtbLiq5seOHVP//v115MgRNW/eXLfffrtSUlJ0+PBhHTt2TP/73/8svhLP4Yqanzt3Tt26ddORI0fUqFEj1a9fXykpKVq5cqUOHz6spk2basqUKU79vi3u7r77bsXHxysoKEj+/v6Kj4/PV9gutr87TXicl19+2QwNDTWnT5/u2JaRkWG+8sorZmhoqPntt986dZ65c+eaO3bsyLJ90aJFZmhoqPnAAw+4rM2eIDU11Wzfvr1Zt25dc/v27Y7t58+fNzt06GDWqVPHPHDgQK7nWbdunRkaGmo++uij5qVLlxzb165da9aqVcvs3bu3Fc33OK6q9+TJk80//vgj07aMjAxz1KhRZmhoqDl69GhXN91juarmV9q2bZtZp04dc9q0aWZoaKh57733urjVns1VNU9LSzMfeughs169eua6deuyfR7YuKrmn376qRkaGmq+8cYbmbZfunTJfOihh8zQ0FAzJibG1c33SGvWrDGPHj1qmqbt/8mhoaHm3Llz83SO4vy7k49jHiYpKUnR0dG6+eab9cgjjzi2G4ahf/3rX/Ly8tKcOXOcOle3bt1Uu3btLNs7dOigatWqae/evUpISHBZ24u69evX6/Dhw+rYsaPq1Knj2B4QEKBBgwYpLS1N8+bNy/U89vq/8MILmYb0NGvWTC1bttTGjRt14MAB11+Ah3FVvZ9++mlVrFgx0zbDMDRo0CBJ0saNG13bcA/mqprbpaSkaPjw4apfv7769OljRZM9nqtqvmTJEm3dulVPPPGEmjZtmmW/jw+jQu1cVfMjR45Iklq3bp1pu6+vr1q0aCFJOn36tAtb7rmaN2+u4ODgAp2jOP/uJGx7mM2bNyslJUXNmzeXYRiZ9lWsWFGhoaH67bffdOnSpQI9T6lSpSSVrP+Bx8TESJJatmyZZZ/9f6z2Y65lw4YN8vf3V1hYWJZ99nMTAF1X75zY37ve3t75Pkdx4+qaT5w4UYcOHdJ///vfLP8/go2rah4dHS1Juu+++3T8+HF98803+vTTT/X999/rwoULLmyx53NVzWvWrClJWrVqVabtqampWrt2rfz8/NSwYcOCNhd/K86/O0tOkiomDh06JEk53ihw6623aufOnTpy5Ei+xzZt2bJFe/bs0Z133lmibi47ePCgJFsNrxYYGKigoCBH/XOSnJysU6dOKTQ0NNuQZ3/d7M9Vkrmi3tcyd+5cSZd/ucK1Nd+yZYs+//xzDR06VNWrV3dlM4sVV9V827ZtkqTY2FiNHz9eKSkpjn3ly5fXhAkT1KRJE9c02sO5quY9evTQggUL9OWXX2rbtm2qW7euUlNTtWrVKp09e1bvvfeeKlWq5Orml0jF/XcnPdse5vz585KksmXLZrs/ICAg03H5Of8rr7wiLy8vvfTSS/lrpIdKSkqSdO3a5lZX+37765DdOa58rpLMFfXOyY4dO/Txxx/rhhtu0IABA/LdxuLGVTVPSUnRiBEjdPvtt+uJJ55waRuLG1fV3D5cYdy4cerfv79WrFihdevWaeTIkTp//ryee+45nTx50nUN92Cuqrmfn5+ioqLUuXNnxcTE6Msvv1RUVJRjiEp2PbDIn+L+u5Oe7ULSpEkTnTlzxunjp02bZnmvxaVLlzR48GDt379fQ4cOpZcEHunIkSN65plnlJ6ervfff1/ly5cv7CYVOxMmTNChQ4c0d+5chum4ifn3xGFt2rTRsGHDHNv79u2rP/74Q5999pm+/fZbx70KKLiEhAQNGjRICQkJ+vTTTxUWFqZLly5p+fLlevPNN/XLL79o7ty5CgwMLOymoogjbBeSjh075mmc3Y033ijp8if1nD6V5/aJPicpKSl67rnntH79ej3zzDP5mv7L0+X2rUBSUlKudbXvz+nTt317Tp/eSxJX1Ptq8fHx6t+/vxISEhQZGZntjWQlmStqvn37dk2dOlWDBg1SrVq1XN7G4sZV7/OAgAAlJibq7rvvzrKvbdu2+uyzzxxDTUo6V9X8zTff1K+//qoFCxY4JhMoW7asevbsqfT0dI0ePVpfffWVnn/+edc1voQq7r87CduF5LXXXsvX4+xj0HIat3To0CF5eXnp5ptvdvqcly5d0qBBg7R69WoNGDBA//rXv/LVNk9nHxN26NAh1a1bN9O+s2fPKjExMdebYfz9/VWhQgUdPXpU6enpWXr97K+bJ0/O7yquqPeVjh49qn79+unkyZOaMGGC2rZt68rmFguuqPmuXbuUnp6uyMhIRUZGZtl/4MAB1apVS2XLltWmTZtc1nZP5ar3efXq1ZWYmJjtfTT2bQW9Mb64cFXNV6xYoXLlymU7a5f9g/z27dsL3mAU+9+djNn2MA0aNJCvr6/Wrl3r+FrR7uTJk9q9e7fq16/v9KpWVwbtJ554osSN075S48aNJUmrV6/Osm/NmjWSpIiIiFzPExERoeTkZMXFxWXZZz+3/blKMlfVW8octD/44AO1b9/edQ0tRlxR82rVqql79+7Z/ifZeqi6d++uLl26uLbxHspV73N7uNu7d2+WffZtBZ16rbhwVc1TUlKUlJSU6WZUO/u0uM6u2IzcFevfnYU90TfyLq+L2iQnJ5t79+414+PjM22/ePGi+cQTT5ihoaHm+PHj3dL2oiw1NdVs166dWbduXfP33393bL9yIYT9+/c7tp8+fdrcu3evefr06UznKc4T87uSq+p95MgRs23btmadOnXMJUuWuK39nshVNc8Ji9pk5aqaHz582Kxbt67ZrFkz88SJE5nO8+CDD5qhoaHm2rVrrb8gD+Cqmtt/P37wwQeZtl+6dMmxLyoqytJr8US5LWpTEn93sly7Bzp58qR69uypEydO6J577lG1atW0adMmxcXFqWXLlvrss88yLR+7YcMG9evXTxEREYqKinJsHz58uObPn68KFSro4Ycfzva5unbtqptuusnyayoq1q9frwEDBqhUqVLq2LGjAgICHEv8vvDCC3r22Wcdx0ZGRmrixIkaPHiwhgwZkuk8I0eO1Jw5c4rfkrMu5op625cJbtCgQbbz6krK8vqUZK56j2enVq1aLNeeDVfVPCoqSuPGjVO5cuV0zz33yNfXV7/88ovi4+P18MMPa+zYse6+tCLLFTXfsWOHevfurQsXLqhevXqOGyRXr16tI0eO6I477tA333zj9DfJxdmcOXMUGxsrSdq9e7e2b9+usLAwx9DX9u3bO75xLIm/Oxmz7YEqVqyo2bNna8KECVqxYoV+/vlnVa1aVUOGDNHTTz+dKWhfS3x8vCTp1KlTmjhxYrbHRERElKiw3bRpU82YMUMfffSRvv/+e6WmpiokJET//Oc/1blzZ6fPM3bsWNWqVUuzZs1SVFSU/P391bZtW+Ykvoor6m1/H2/evFmbN2/O9hjC9mWueo/Dea6qed++fRUcHKwvvvhCixcvVnp6ukJCQjRw4ED17NnTwivwPK6o+e2336558+Zp8uTJWr9+vaZPny5vb2/dcsstGjJkiJ588kmC9t9iY2M1f/78TNvi4uIcQ0KCg4OdGt5XXH930rMNAAAAWIQbJAEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIsQtgEAAACLELYBAAAAixC2AQAAAIuwXDsA5KJWrVqZ/m4YhgICAhQaGqquXbuqe/fuMgwj0zHnzp3Tl19+qZ9//lmHDh1SRkaGypcvr+DgYIWHh+u+++7T7bff7jh+3rx5GjFihLp27ao333wz07kuXbqkqKgo/fDDD9q/f79SUlIUFBSkKlWqKCwsTO3atVPjxo2dupbhw4dnWVa5dOnSqlq1qlq1aqWnnnpKFSpUyEt5ip1atWopODhYy5cvL+ymACgGCNsA4KSuXbtKktLT03XkyBHFxcUpNjZW69at0/vvv+84Lj4+Xn369NGxY8d03XXXqWHDhipfvrzOnDmjbdu2adOmTUpMTNTYsWNzfc6zZ8+qf//+2rFjh3x9fdWgQQNVrFhRSUlJ2r59u3777Tft2bPH6bBtFxYWpltvvVWSdPr0af3222/66quvFB0drZkzZ+qmm27K0/msFBkZqYkTJ2r8+PHq1q1bYTcHAPKEsA0ATrq6x3nNmjV6+umntXjxYnXq1Elt27aVJL3++us6duyY2rZtq3feeUdly5Z1PCYtLU1r1qzR6dOnnXrOjz76SDt27FDdunU1efJk3XjjjY59pmkqNjZW27dvz/O19OjRI1NwTUhI0FNPPaVt27bp7bff1kcffZTncwIAsmLMNgDkU4sWLdS5c2dJ0rJlyyRJFy9e1KpVqyRJr776aqagLUk+Pj5q3bq10z20S5culSQNGzYsU9CWbMNZwsPD1b9//wJdhySVL19ew4cPlyT98ssvSk1NLfA5AQCEbQAokDp16kiSTpw4Ick27CMtLU2SLcAWVGJioiQpKCiowOfKjX0M+aVLlxzP+8svv2jEiBG6//77FRYWpgYNGqhz58765JNPlJKSkuUc8+bNU61atRQZGakDBw5o6NChat68uWrXru34QCJJu3bt0osvvqi77rpLdevWVcuWLTVixAgdPXo00/nuvvtuTZw4UZI0YsQI1apVy/Hfhg0bMh373Xff6ZFHHlFYWJjq16+vTp06afLkybp06ZJL6rN79261bNlSdevWVXR0tEvOCaD4I2wDQAFcuHBBklSqVClJtlBcunRpSdI333xT4PNXqlRJkjRr1iyZplng812L/VokydfXV5L073//Wz/88IPKli2ru+66S40aNdKJEyf0wQcf6KmnnlJ6enq25zpw4IC6d++uLVu2qEmTJmrRooV8fGwjF5csWaKHHnpIixYtUoUKFXT33XerQoUKmjdvnh566CHt2bPHcZ57771XtWvXlmQbZ961a1fHf1f29P/nP//RK6+8ou3btys8PFytW7fWqVOn9P7776t///66ePFigWrz66+/qk+fPkpKStKkSZP0wAMPFOh8AEoOxmwDQD6ZpqlffvlF0uUZS3x9ffXggw9q9uzZevfdd/XDDz+odevWatCggRo2bJhlWEluevTooQ8++EAzZszQ+vXr1a5dOzVo0EBhYWEu6Tm/0s8//yzJFvDLlSsnSRozZoyaN28uf39/x3FJSUkaNmyYfv75Z/3vf/9Tly5dspxr8eLF6tOnj1599VV5e3s7th85ckSvvPKK/Pz8NGXKlEw3dn733Xd65ZVXNGLECH377beSpFdeeUWRkZHauXNnlnHmdkuWLNGsWbNUqVIlRUVFOW78TEpK0tNPP63Y2Fh99NFHevnll/NVl1WrVun555+Xj4+PvvjiCzVq1Chf5wFQMtGzDQB5lJ6eroMHD+rVV1/Vr7/+Kl9fXz300EOO/a+++qoefPBBGYahbdu26eOPP9ZTTz2lJk2aqG/fvlq9erXTz/XUU0/p8ccfl4+Pj/bv36/PPvtMzz33nJo3b+7oHS6ohIQEzZ07V++8844k6ZFHHnHsa9++faagLUkBAQEaMWKEJOmnn37K9pzly5fXsGHDMgVtSZo2bZr++usvvfTSS1lmUOnSpYvat2+vrVu35ummz6ioKEnS888/7wja9naOGjVKhmFo5syZ2Q57yU10dLSeffZZXXfddYqKiiJoA8gzerYBwElXz7ctSdddd53eeust3XLLLY5tZcqU0dtvv62BAwdqyZIlio2N1datW3XmzBnFxMQoJiZGw4cP1+OPP57rc3p7e2v48OHq37+/fvjhB23cuFFbtmzRqVOntG3bNr344ovavHmzRo4cmadrGTFihCMwX6lr1656+umnM207ePCgVqxYocOHDys5OVmmaTqGtBw8eDDb8zdv3lxlypTJsn3t2rWSpHbt2mX7uEaNGmnZsmXaunWr7rjjjlyvIzU1VZs3b5ZhGOrUqVOW/fbx3Tt37tTOnTtVr169XM9p980332js2LGqWrWqpkyZkuk1BgBnEbYBwEn2ebavXNTmH//4hwIDA7M9vkaNGnr22WclSRkZGfr111/1/vvva9OmTXr33Xf1j3/8Q8HBwU49d5UqVfT44487AvrOnTsVGRmpZcuWKSoqSvfff3+eel2vnGfb19dXwcHBatWqVaaFdkzT1FtvvaWpU6fmOF78ynHeV7c3O/Hx8ZJsM7lci/0GzdycOXNGqampqlChgmOs/NWCg4O1c+dOnTx50qlzSrYbXkePHq3SpUtr2rRpTr9OAHA1wjYAOOnqebbzwsvLS40aNdLnn3+u+++/X8ePH9fq1av18MMP5+t8tWvXVmRkpHr27KmtW7dqxYoVeQrbOY1/vlJ0dLSmTJmiypUr69VXX1WDBg1Uvnx5lSpVSikpKbrzzjtzfGxOwTc9PV2GYWQ7zvtKNWvWzPUa8urqVT6vpXz58goJCdG6dev01ltv6f3333fc4AkAecH/OQDAjcqUKaN69erp+PHjTvfe5sTLy0uNGzfW1q1bC3yu7Pz444+SpNGjRzsW7LE7cuRIvs5ZuXJlHT58WCNHjlRAQECB21iuXDmVKlVKf/75py5evCg/P78sxxw7dkyS8rQMva+vrz755BM9/fTTWrJkiV566SW9++67WcagA0BuuEESANzs8OHDkqSKFSvmemxu0/3l5Vx5de7cOUnZDwn5/vvv83XOZs2aSVKmObdzY59WMbtpBkuVKqUGDRrINM1sbxbdvXu3du7cqeuuu84xhaCz/Pz8NHnyZEVERCg6Olovv/yyMjIy8nQOACBsA4ALnTt3Tj169NCPP/6YZRXG1NRUTZo0STt27JCfn5/uuuuuXM/Xq1cvzZ8/X3/99Vem7aZpau7cufrpp59kGIbuuecel16HJFWrVk1S1jm+N23apC+++CJf53ziiSfk5+en8ePHa/ny5Vn2nzlzRtOnT880L7b9g8T+/fuzPWefPn0kSZGRkZl63JOSkvT666/LNE09/PDDjrnD86JMmTKaPHmywsPDtWjRIg0fPpzADSBPGEYCAC62ZcsWDR48WAEBAbrjjjtUoUIFnTt3Tjt27NCpU6fk7e2t0aNHOzWsYd++fRo+fLhGjx6tO+64Q5UrV1ZycrL27NnjWG3xn//8Z557bZ3Rt29fzZ8/XzNmzFBMTIxq1aqlP/74Q7GxsXr88cf15Zdf5vmc1apV0zvvvKOXXnpJzz77rKpXr67bbrtNpmnq2LFj2rt3r1JTU9WpUyfHkJAWLVqodOnS+uqrr7Rnzx5VrFhRhmHoySefVI0aNXTffffp4Ycf1qxZs9SxY0c1bdpUfn5+iomJUUJCgho0aKDnn38+33Xw9/fXp59+qgEDBmjBggXy9vbWG2+8kacx4ABKLsI2ALhQ2bJlNXPmTK1atUoxMTE6evSo4uLi5O3trSpVqqhNmzbq06eP0+H466+/1sqVK7Vu3TodOXJE27dvl2maqlChgjp27KhHHnlE4eHhllxL9erV9e233+qdd97Rli1btHz5clWvXl1jx45Vz5498xW2Jekf//iHQkND9eWXX2rt2rVauXKlSpcurYoVK6pTp0669957My3+U6lSJU2aNEkff/yxYmNjlZycLEnq3LmzatSoIUkaO3aswsLCNHPmTMXExCg9PV233HKL+vfvr8ceeyzbsdx5cd111+mzzz7TgAEDNG/ePHl7e+v1118ncAPIlWFavf4vAAAAUEIxZhsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALAIYRsAAACwCGEbAAAAsAhhGwAAALDI/wOWX+AhDS4n0QAAAABJRU5ErkJggg==\n",
"text/plain": [
"