BSM/Chapter_03_04_Gibbs_sampling_for_concussions_data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "az.style.use(\"arviz-white\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gibbs sampling for the concussions data\n",
    "\n",
    "Let $Y_i$ be the number of concussions (aggregated over all teams and games) in season $i$ (1=2012,…,4=2015). We model these counts as $Y_i\\sim\\mbox{Poisson}(N\\lambda_i) \\mbox{  where  } \\lambda_i \\mid \\gamma\\sim\\mbox{Gamma}(1,\\gamma),$ $\\lambda_i$ is the concussion rate in year $i$ and $N$ is the number of games in each year. The prior for $\\gamma$ is $\\gamma\\sim\\mbox{Gamma}(a,b)$. The objective is to deterimine if the concussion rate has changed over time by comparing the posteriors of the $\\lambda_i$.\n",
    "\n",
    "Gibbs sampling cycles through the parameters and updates each using a draw from its full conditional distributions. The full conditional distributions are:\n",
    "\n",
    "$$\\lambda_i \\mid \\mbox{rest}\\sim\\mbox{Gamma}(Y_i+1,N+\\gamma)$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\gamma \\mid \\mbox{rest}\\sim\\mbox{Gamma}(a+4,b+\\sum_{i=1}^4\\lambda_i).$$\n",
    "\n",
    "This produces draws from the joint posterior of $(\\lambda_1, \\dots,\\lambda_4,\\gamma)$\n",
    "\n",
    "To evaluate whether the rate has changed between years $i$ and $j$, we approximate the posterior probabilities that $\\lambda_i > \\lambda_j$ using the proportion of the MCMC samples for which this is the case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = np.array([171, 152, 123, 199])\n",
    "n = 4\n",
    "N = 256"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty matrix for the S MCMC samples\n",
    "S = 2500\n",
    "samples = np.zeros((S, 5))\n",
    "# colnames(samples) = c(\"lam1\", \"lam2\", \"lam3\", \"lam4\", \"gamma\")\n",
    "\n",
    "# Initial values\n",
    "λ = np.log(Y / N)\n",
    "γ = 1 / λ.mean()\n",
    "\n",
    "# priors: lambda|gamma ~ Gamma(1,gamma), gamma ~ InvG(a,b)\n",
    "a = 0.1\n",
    "b = 0.1\n",
    "\n",
    "# Gibbs sampling\n",
    "for s in range(S):\n",
    "    λ = stats.gamma(Y + 1, scale=1 / (N + γ)).rvs()\n",
    "    gamma = stats.gamma(a + 4, scale=1 / (b + λ.sum())).rvs(1)[0]\n",
    "    samples[s] = (*λ, gamma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summarize the posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAGbCAYAAAC1c/byAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X9c1fX9//87HPyBHgISLWYqKGA4teav4Q+mViunJIi0dxdyNfulae/ZQiv9uL3bZupa2tJlWW7rh9pWiFgoRmUtMNQpM6MwEUUTQdLAOAIGh/P9wy9nnvgh6OGcF5zb9XLp4jmv55Pn63EM5e7z9Xy+Xl42m80mAAAAeARvdxcAAAAA1yH8AQAAeBDCHwAAgAch/AEAAHgQwh8AAIAHIfwBAAB4EMIfAACAByH8AQAAeBDCHwAAgAch/AEAAHgQH3cX0JQtW7Zo3759ys3N1aFDh1RTU6Nly5YpPj6+VeN89913eumll/TOO+/o5MmT8vf31/jx4/XII4+oZ8+ebVQ9AACAMRk2/D333HMqKipSYGCgevXqpaKiolaPUVdXp4ceekhZWVm64YYb9NOf/lTHjx9XSkqKPvnkE7355psEQAAA4FEMe9l3yZIl2rFjh3bt2qU777zzssbYvHmzsrKyNGXKFP3zn//U/PnztWrVKi1ZskQnT57UM8884+SqAQAAjM2w4W/MmDHq3bv3FY3x1ltvSZKSkpLk5eVlPz59+nQNGDBA27Ztk8ViafRrKysrlZOTo8rKyiuqAQAAwEgMG/6u1Pnz5/Xpp58qNDS00RA5duxYfffdd/r0008b/fqDBw9q+PDhOnjwYFuXCgAA4DIdNvwdP35cdXV1CgkJabS9/nhhYaHLagIAAHC3Dhv+KioqJElms7nR9u7du0tSk5d9AQAAOqIOG/4AAADQUIcNf35+fpKantk7d+6cpKZnBgEAADqiDhv++vbtK29v7ybX9NUfb2pNIAAAQEfUYcNfly5dNHToUB09erTRG0Tv3LlTnTt31g033OCG6gAAANyjQ4S/0tJSFRQU2Dd51Pv5z38uSVqxYoVsNpv9+KZNm1RQUKDJkydz2RcAAHgUwz7e7a233tK+ffskSYcOHbIf27NnjyTplltu0S233CJJWrlypTZv3tzg2b9xcXHatm2btm7dqhMnTmjUqFH66quvlJGRoeDgYM2fP9/FnwoAAMC9DBv+9u3bp82bNzscy8nJUU5OjiSpd+/e9vDXFJPJpBdeeEEvvfSS3n77bb3yyivy9/fXtGnT9Mgjj/BcXwAA4HG8bBdfD4VdTk6Ohg8frn379mnYsGHuLgcAAMApDDvzBwAA8H1Wq1WZmZkqLi5WcHCwoqOjZTKZ3F1Wu9IhNnwAAICOLyUlRWFhYZo4caISExM1ceJEhYWFKSUlxd2ltSuEPwAAYHgpKSlKSEjQkCFDlJ2drYqKCmVnZ2vIkCFKSEggALYC4Q8AABia1WpVUlKSYmJitGnTJlVXV+udd95RdXW1Nm3apJiYGM2fP19Wq9XdpbYLhD8AAGBomZmZKiws1JgxYxQeHu5w2Tc8PFyjR4/W0aNHlZmZ6e5S2wU2fAAAAEMrLi6WJC1cuFC+vr4ObaWlpVq0aJFDPzSP8AcAAAytV69e9tc33XSTJk+eLF9fX1VVVdkf5vD9fmga4Q8AABha/Vo+s9msAwcO2MOeJPXp00dms1kWi4U1fy3Emj8AAGBo9Wv5LBaLioqKHNqKiopksVgc+qF5hD8AAGBodXV19tddunRxaOvatWuj/dA0LvsCAABDCwgIkCRdddVVKi0tVXZ2tv0JH6NHj1bPnj1VUVFh74fmMfMHAAAMrby8XJL07bffavr06fr8889VVVWlzz//XNOnT1dFRYVDPzSPmT8AAGBo3t7/navaunWrw4aPpvqhafwuAQAAQ5swYYJT+3k6wh8AADC06OhoeXl5NdvHy8tL0dHRLqqofSP8AQAAQ8vMzJTNZmu2j81m41YvLUT4AwAAhvb+++9Lkrp3766+ffs6tPXr10/du3d36IfmEf4AAICh7du3T5I0a9Ys5efn69lnn9XDDz+sZ599VocOHdIDDzzg0A/NY7cvAAAwtG7dukmStmzZopSUFBUWFtrbnnvuOZlMJod+aB4zfwAAwNDqN3IUFBTo+PHjDm3Hjx9XQUGBQz80j/AHAAAMbc6cOfbX33+E28XvL+6HphH+AACAoX3yySdO7efpCH8AAMDQPvroI6f283Rs+AAAAIZmtVrtrzt16qSf/OQnCg4OVnFxsT7++GPV1NQ06IemEf4AAIChnTp1yv7622+/VdeuXe3vq6ur5evr26AfmsZlXwAAYGh79uyxv77jjjuUnZ2tiooKZWdn64477mi0H5rGzB8AAHC5I0eOqLy8vEV9Kyoq7K/fe+89paWl2d9fPAtYUVGhnJycS44XEBCg/v37t6LajoXwBwAAXOr06dMKDw9vcNuWljh//rzD++rqavvrY8eOafjw4Zccw2QyqaSkREFBQa0+f0dA+AMAAC4VFBSk/Pz8Fs/8VVVVady4cZIuzNqNHj1a6enp+tnPfqbs7Gz7OFlZWfb1f80JCAjw2OAnEf4AAIAbtPaya2xsrLZs2aLy8nKlp6dLkv3X+vaxY8c6tcaOig0fAADA8FJTUxUbG9toW2xsrFJTU11cUftF+AMAAO1CamqqKisr7Tt877jjDlVWVhL8WonwBwAA2g1fX1898cQTkqQnnniiRWv84IjwBwAA4EEIfwAAAB6E8AcAAOBBCH8AAAAehPAHAADgQQx9k+cDBw5o9erV2r9/v2pqahQWFqZ77rlHt99+e4vHKCgo0Jo1a7Rr1y6dPXtWPXv21M0336yHH35YAQEBbVg9AACA8Rg2/O3evVv33XefOnXqpClTpsjPz08ZGRmaP3++ioqKNHv27EuOsX//fs2cOVPV1dW6+eab1adPHx08eFCvv/66MjMz9Y9//EOBgYEu+DQAAADGYMjwV1tbq8WLF8vLy0sbNmzQoEGDJElz587VnXfeqdWrV2vSpEkKCQlpdpzFixersrJSa9as0c0332w/vm7dOv3pT3/Ss88+q9///vdt+VEAAAAMxZBr/nbt2qXjx48rJibGHvwkyWw2a86cOaqtrVVKSkqzYxw7dkz5+fkaMmSIQ/CTpHvvvVcBAQF6++23VVlZ2SafAQAAwIgMGf727NkjSRo3blyDtvqHNtf3acrp06clSdddd12DNm9vb/3gBz9QVVWVPv300ystFwAAoN0wZPgrLCyUJPXr169Bm7+/vwIDA3Xs2LFmx7j66qslSSdOnGjQVldXp5MnT0qSjh49eoXVAgAAtB+GDH8Wi0WS5Ofn12i72WxWRUVFs2OEhoaqb9+++uyzz/TRRx85tL3yyisqLy+XpEuOAwAA0JEYcsOHs/z2t7/VQw89pIceeki33HKL+vTpoy+//FJZWVmKiIjQoUOH5O1tyPwLAADQJgyZfMxms6SmZ+UsFkuTs4IXi46O1oYNG/STn/xEu3bt0uuvv66ysjI9//zzGjVqlKT/Xh4GAADwBIac+au/hcuxY8c0ePBgh7azZ8+qrKxMP/rRj1o01g033KC1a9c2OP7qq69KUoPxAQAAOjJDzvyNHDlSkpSVldWgbefOnZJkn7m7HEVFRdq3b5/CwsI0cODAyx4HAACgvTFk+Bs9erT69OmjtLQ05eXl2Y9bLBatWbNGPj4+mjZtmv14aWmpCgoKGlwmPnfunGw2m8OxiooKPfbYY7JarXr00Ufb9oMAAAAYjCEv+/r4+GjJkiW6//77lZiYqJiYGJnNZmVkZOjEiRN65JFHFBoaau+/cuVKbd68WcuWLVN8fLz9+Pvvv69nn31WUVFR6tWrl86cOaMdO3bom2++0bx58xrc/BkAAKCjM2T4k6SoqCht3LhRq1atUnp6umpqahQWFqZ58+Zp6tSpLRpj4MCBuv7665WVlaXy8nKZzWbdeOON+uUvf6moqKg2/gQAAADG42X7/nVRSJJycnI0fPhw7du3T8OGDXN3OQAA4P/Hz+grY8g1fwAAAGgbhD8AAAAPQvgDAADwIIQ/AAAAD0L4AwAA8CCEPwAAAA9C+AMAAPAghD8AAAAPQvgDAADwIIQ/AAAAD0L4AwAA8CCEPwAAAA9C+AMAAPAgPu4uAO2H1WpVZmamiouLFRwcrOjoaJlMJneXBQAAWoGZP7RISkqKBgwYoIkTJyoxMVETJ07UgAEDlJKS4u7SAABAKxD+cEkpKSmaPn26SktLHY6XlpZq+vTpBEAAANoRwh+aZbVaNXv27Gb7PPTQQ7JarS6qCAAAXAnCH5r10Ucf6euvv5Yk3XzzzcrOzlZFRYWys7N18803S7owA/jRRx+5sUoAANBShD80a8eOHZKkqKgovfXWW9q1a5cWLlyoXbt26a233lJUVJRDPwAAYGzs9kWzvvrqK0lS9+7d5efnp9raWnvbggULFB0d7dAPAAAYG+EPzerTp48k6YMPPlCvXr00aNAg1dXVydvbW1988YU+/PBDh34AAMDYCH9o1rhx4+yvS0tLG+z4bawfAAAwLtb8oVlffPGFU/sBAAD3IvyhWYcPH7a/9vX1dWi7+P3F/QAAgHER/tCszz77TJL0ox/9SD169HBoCwoK0o033ujQDwAAGBtr/tCs+tm93Nxc1dTUOLR99dVXKikpcegHAACMjZk/NGvgwIGSZA9+kyZN0s6dOzVp0iSH4/X9AACAsRH+0Kzf//739tfe3t7avn27xo4dq+3bt8vb27vRfgAAwLgIf2jWb3/7W/trLy8vDRgwQBERERowYIC8vLwa7QcAAIyLNX9oVn5+viSpb9++On78uAoKChza64/X9wMAAMZG+EOzwsPDlZGRoePHjzfaXn88PDzclWUBANwsPz9fFRUVbjl3Xl6ew6+u5ufn165/7nnZbDabu4swopycHA0fPlz79u3TsGHD3F2O25w9e1YBAQGSLlz2vfjb5eL35eXl8vf3d0uNAADXys/PV0REhLvLcKtDhw612wDIzB+a9dJLL9lf22w2zZgxQ48++qhWrlyp9evXO/RbsGCBO0oEALhY/Yzf+vXrFRkZ6ZYaysrKFBgY6PLz5uXlacaMGW6b9XQGwh+atWXLFkmS2WyWxWLR+vXrHUJf/fEtW7YQ/gDAw0RGRnr01bH2it2+aNbZs2clSTNnzlRlZaXmzp2rW2+9VXPnzlVlZaXuvvtuh34AAMDYmPnzQEeOHFF5eXmL+l533XXKzc3VX//6V82YMUP33nuvfar9s88+0yuvvGLvl5OT06IxAwIC1L9//8stHwAAXAHCn4c5ffq0wsPDVVdX16qvq6ys1I9//OMm27dv367t27e3aCyTyaSSkhIFBQW1qgYAAHDlDB3+Dhw4oNWrV2v//v2qqalRWFiY7rnnHt1+++0tHuPUqVN6+eWX9cknn+jkyZPq1q2b+vXrp//5n//R7bffLpPJ1IafwHiCgoKUn5/f4pk/q9WqiRMn6ty5c032MZvN2rFjR4t/LwMCAgh+AAC4iWHD3+7du3XfffepU6dOmjJlivz8/JSRkaH58+erqKhIs2fPvuQYX331le644w6Vl5dr3LhxmjhxoiwWiz744AM9/vjj2r17t5YtW+aCT2Msrb3k+tprr2n69OlNtr/66qsaOXLklZYFAABcwJAbPmpra7V48WJ5eXlpw4YNWrJkiR5//HFt2bJF4eHhWr16tQoLCy85zl//+leVlZVp0aJFWrdunRYsWKDf/e532rZtm3r37q2UlBQVFRW1/Qdq5+Lj47Vp0yb17dvX4Xi/fv20adMmxcfHu6kyAADQWoYMf7t27dLx48cVExOjQYMG2Y+bzWbNmTNHtbW1SklJueQ4X331lSRp/PjxDsevuuoq/ehHP5IkffPNN06svOOKj4/XkSNHtHbtWknS2rVrVVBQQPADAKCdMWT427NnjyRp3LhxDdrGjh3r0Kc59XfezszMdDhusVj0n//8R0FBQQoLC7vScj2GyWTSiBEjJEkjRozwuPWSAAB0BIZc81d/Sbdfv34N2vz9/RUYGKhjx45dcpz77rtPO3bs0FNPPaWPP/5YERER9jV/nTp10qpVq+Tr6+vs8gEAAAzLkOHPYrFIuvDg5MaYzWaVlJRccpyePXvqH//4hxYsWKB//etf+te//iVJ6tKli+6//36HS8oAnMtqtSozM1PFxcUKDg5WdHQ0s8UAYACGDH/Ocvz4cc2aNUvdunXThg0bFBkZqYqKCr3zzjv685//rMzMTG3cuFGdOnVyd6lAh5KSkqKkpCSHjVkhISFasWIF60QBwM0MuebPbDZLUpMPTbZYLE3OCl5s4cKFOnnypF588UWNGDFC3bt317XXXqsHHnhAM2fO1IEDB/T22287tXbA06WkpCghIUFDhgxRdna2KioqlJ2drSFDhighIaFFm7UAAG3HkOEvJCREkhpd13f27FmVlZU1uh7wYhaLRXv37tWAAQPUs2fPBu1RUVGSpM8///zKCwYg6cKl3qSkJMXExCg1NVVRUVEym82KiopSamqqYmJiNH/+fFmtVneXCgAey5Dhr/6GwVlZWQ3adu7cKUkaNWpUs2PU1NRIksrKyhptr7/FS+fOnS+7TgCOMjMzVVhYqEWLFsnb2/GvF29vby1cuFBHjx5tsAMfAOA6hgx/o0ePVp8+fZSWlqa8vDz7cYvFojVr1sjHx0fTpk2zHy8tLVVBQYHDZeLAwECFhobq5MmTeuuttxzGt1gsWrdunSQ1+7xaAK1TXFwsSRo8eHCj7fXH6/sBAFzPkBs+fHx8tGTJEt1///1KTExUTEyMzGazMjIydOLECT3yyCMKDQ2191+5cqU2b96sZcuWOSwmX7RokR566CEtXrxYaWlpGjRokCoqKvThhx/q9OnTmjhxoiZMmOCGTwh0TMHBwZKk3Nxc+9KKi+Xm5jr0A9B+hQaY1LX8sGwnDRkl2kzX8sMKDWjfdy4w7P+xqKgobdy4UatWrVJ6erpqamoUFhamefPmaerUqS0a4yc/+Yn++c9/at26ddq7d6/27t2rzp07a8CAAXrwwQd11113ycvLq40/CeA5oqOjFRISoqVLlyo1NdXh0m9dXZ2WLVum0NBQRUdHu7FKAFfKdP6svpjTQ6YPf6XaD91djWuFS/p8Tg8dOn/W3aVcNsOGP0kaOnSo/fJsc5YvX67ly5c32jZ48GD9+c9/dnZpABphMpm0YsUKJSQkKC4uTgsXLtTgwYOVm5urZcuWKS0tTcnJydzvD2jnrF38NWjNGaVt2qjI6693dzkulXfwoGKmJyp5ur+7S7lshg5/ANqf+Ph4JScnKykpSWPGjLEfDw0NVXJyMvf5AzqIo+VWVQeEyesHQ91diktVl9TqaHn7vmMB4Q+A08XHxys2NpYnfACAARH+ALQJk8nEhioAMCDCH4A2wbN9AcCYDHmfPwDtW0pKisLCwjRx4kQlJiZq4sSJCgsL49FuAGAAzPwBcKr6Z/tOmTJFCxYskK+vr6qqqpSenq6EhAQ2fQCAmxH+ADhN/bN9hw8frs8++0xpaWn2tn79+mn48OGaP3++YmNjuQQMAG5C+APgNPXP9i0sLJSvr69DW2lpqY4dO2bvx2YQAHAPwh8ApykqKrK/vummmzR58mT7Zd9t27Zp69atDfoBAFyL8AfAaU6dOiXpwiXe3Nxce9irP9avXz8dO3bM3g8A4Hrs9gXgNGfOnJEkHTt2TEOGDFF2drYqKiqUnZ2tIUOG2C/71vcDALge4Q9Am7HZbPb/AADGwGVfAE5z9dVXS7pwifezzz5zeLZvSEiI/bJvfT8AgOsR/gA4zbXXXivpwmXfyZMnKzY2VtXV1eratasOHz6sbdu2OfQDALge4Q+A0/Tu3dv+uj7oXaofAMC1WPMHwGmio6PVs2fPZvv06tVL0dHRLqoIAPB9hD8ATnXu3Llm2y0Wi4sqAQA0hsu+AJzmgw8+UGVlpSSpa9euqq6utrfVv6+srNQHH3ygW2+91V1lAnCSnJwct527rKxMgYGBLj9vXl6ey8/pbIQ/AE7z2muvSZL69u0rSTp+/Li9rVevXqqrq9OJEyf02muvEf6Adqy2tlaS9MADD7i5Evfx8/NzdwmXjfAHwGkKCwslXQh9Xbp0cWg7deqUzp8/79APQPs0atQo7d69Wz4+7okReXl5mjFjhtavX6/IyEiXn9/Pz0/h4eEuP6+zEP4AOE3fvn21c+dOSdJ3333n0Hbx+/qZQQDt16hRo9xdgiIjIzVs2DB3l9HusOEDgNPceOON9tfff6rHxe8v7gcAcC3CHwCnaekze3m2LwC4D+EPgNPs3bvXqf0AAM5H+APgNKdOnZIkBQUFqV+/fg5tISEhCgoKcugHAHA9NnwAcJqAgABJ0unTp/XjH/9YCxYskK+vr6qqqpSenq6tW7c69AMAuB7hD0Czjhw5ovLy8hb1HT58uH23b0ZGhj3sSVLnzp0d+rXk5rABAQHq379/KysGADSH8AegSadPn1Z4eLjq6upa/bU1NTUO7y++1cuqVau0atWqS45hMplUUlJiv1wMALhyhD8ATQoKClJ+fn6LZ/4k6bnnnrM/6aMxd999t+bNm9eisQICAgh+AOBkhD8AzWrtZddXX31V11xzjVasWOEwY+jt7a2kpCQ9/fTTzi4RANAK7PYF4HRPP/20qqqq9Oijj0qSHn30UVVVVRH8AMAACH8A2kTnzp111113SZLuuusuhw0fAAD3IfwBAAB4EMIfAACAByH8AQAAeBDCHwAAgAfhVi/tVH5+vioqKlx+3ry8PIdf3cHPz0/h4eFuOz8AAO0Z4a8dys/PV0REhFtrmDFjhlvPf+jQIQIgAACXwdDh78CBA1q9erX279+vmpoahYWF6Z577tHtt9/eoq+/6aabVFRU1GyfDRs2aMSIEc4o12XqZ/zWr1+vyMhIl5+/rKxMgYGBLj+vdGHGccaMGW6Z9QQAoCMwbPjbvXu37rvvPnXq1ElTpkyRn5+fMjIyNH/+fBUVFWn27NmXHOPuu+9uNCSUlZVpw4YN8vf315AhQ9qifJeIjIzUsGHD3F0GAABoRwwZ/mpra7V48WJ5eXlpw4YNGjRokCRp7ty5uvPOO7V69WpNmjRJISEhzY7zy1/+stHjf/vb3yRJU6dOVZcuXZxZOgAAgKEZcrfvrl27dPz4ccXExNiDnySZzWbNmTNHtbW1SklJuezxk5OTJUkJCQlXXCsAAEB7Ysjwt2fPHknSuHHjGrSNHTvWoU9r5eTkqKCgQIMHD9b1119/+UUCAAC0Q4YMf4WFhZKkfv36NWjz9/dXYGCgjh07dllj18/63XHHHZddHwAAQHtlyPBnsVgkXbifW2PMZvNl7fY8d+6c0tPT5evrq5iYmCuqEQAAoD0yZPhrK9u2bVNlZaUmTZoks9ns7nIAAABczpDhrz6YNTW7Z7FYmpwVbA4bPQAAgKczZPirv4VLY+v6zp49q7KyskbXAzbn8OHD2r9/v/r379/ubuoMAADgLIYMfyNHjpQkZWVlNWjbuXOnJGnUqFGtGpNZPwAAAIOGv9GjR6tPnz5KS0tTXl6e/bjFYtGaNWvk4+OjadOm2Y+XlpaqoKCgycvENTU12rJlizp16qS4uLg2rx8AAMCoDBn+fHx8tGTJEtlsNiUmJuo3v/mN/vjHPyo2Nlb5+fl6+OGHFRoaau+/cuVKTZ48We+9916j4+3YsUPffPONJk6cqB49erjqYwAAABiOIR/vJklRUVHauHGjVq1apfT0dNXU1CgsLEzz5s3T1KlTWzUWl3wBAAAuMGz4k6ShQ4dq3bp1l+y3fPlyLV++vMn2l19+2ZllAQAAtFuGvOwLAACAtkH4AwAA8CCGvuyLpoUGmNS1/LBsJz3rf2HX8sMKDTC5uwwAANotz0oOHYTp/Fl9MaeHTB/+SrUfursa1wqX9PmcHjp0/qy7SwEAoF0i/LVD1i7+GrTmjNI2bVTk9de7uxyXyjt4UDHTE5U83d/dpQAA0C4R/tqpo+VWVQeEyesHQ91diktVl9TqaLnV3WUAANBuseEDAADAgxD+AAAAPEibh7/c3FwtWbKkrU8DAACAFmiTNX+nTp3S22+/rdTUVB05ckSStHjx4rY4FQAAAFrBaeHv/PnzysjIUGpqqnbt2qW6ujrZbDbdcMMNmjZtmrNOAwDoYKxWqzIzM1VcXKzg4GBFR0fLZOJ+nkBbueLwt2fPHqWmpurdd99VZWWlbDabrrnmGk2dOlXTpk1T//79nVEnAKADSklJUVJSkgoLC+3HQkJCtGLFCsXHx7uvMKADu6zwd/z4caWmpmrLli06efKkJKlz586aPHmypk2bprFjx8rLy8uphQIAOpaUlBQlJCQoJiZGb7zxhgYPHqzc3FwtXbpUCQkJSk5OJgACbaBV4W/r1q1av3699u/fL0my2WySpODgYK1Zs0aRkZHOrxAA0OFYrVYlJSUpJiZGqamp8va+sP8wKipKqampiouL0/z58xUbG8slYMDJWrzb9+WXX9b8+fP1n//8R9ddd53mzp2ru+66S76+viouLtYdd9yhpKQk5eXltWW9AIAOIDMzU4WFhVq0aJE9+NXz9vbWwoULdfToUWVmZrqpQqDjavHM3xtvvKGuXbvqT3/6k2655Rb78V/96lfasGGDNm7cqK1bt2rbtm0aM2aM7r//fo0ePbpNigYAtG/FxcWSpMGDBzfaXn+8vh8A52nxzF9JSYkmTpzoEPwkyd/fX3PmzNGOHTv0hz/8Qf3799fOnTt17733Kj4+Xtu2bbNfHgYAQLqwXEi6cC/YxtQfr+8HwHlaHP5iY2M1ZsyYJts7d+6sO+64Q1u3btWLL76oESNG6IsvvlBSUpJuvfVWpxQLAOgYoqOjFRISoqVLl6qurs6hra6uTsuWLVNoaKiio6PdVCHQcbU4/C1btkwJCQkt6jthwgS9/vrrSklJ0c9+9jOm7QEADkwmk1asWKG0tDTFxcUpOztbFRUVys7OVlxcnNLS0vTMM8+w2QMNWK3L0wsmAAAgAElEQVRW7d27V5K0d+9eWa1WN1fU/rTJEz7qDRo0SCtXriT8tZGcnBy3nLesrEyBgYFuOTcbioCOIz4+XsnJyUpKSnK4shQaGsptXtCo798XctasWVq2bBn3hWylNg1/9Viz4Vy1tbWSpAceeMDNlbiPn5+fu0sA4ATx8fGKiYnRmjVrVFBQoAEDBmjOnDnq3Lmzu0tDGzty5IjKy8tb3H/Hjh167LHHFB0drV/+8pd68skn9eSTT2rHjh1KSEjQ008/rZtuuqlFYwUEBHj0Qyi8bOzGaFROTo6GDx+uffv2adiwYe4up4E9e/bIx8cl2d1BXl6eZsyYofXr17vtvo5+fn4KDw93y7nROkb/cwT34wkfnun06dO65pprGqz3dBWTyaSSkhIFBQW55fzu5vr0AKcYNWqUW88fGRnJD3MAV4QnfHiuoKAg5efnt3jmb+/evZo1a5ZeeeUVDRkyRJLjEqQDBw5o5syZWrt2rUaMGHHJ8QICAjw2+EmEPwCAG1z8hI9NmzZp586deueddxQcHKxNmzZp+vTpPOGjg2vNZdcvv/xSkjR9+nSZzeYG7eHh4Zo5c6b8/PyYmGgBwh/QweXn56uiosIt567foOOujTosETCu+id8zJo1SxEREQ0u+z744IN65513lJmZqQkTJritThjDxfeFjIqKatDOfSFbh/AHdGD5+fmKiIhwdxmaMWOG28596NAhAqAB1d8FYtGiRY1e9v1//+//OfSDZ7v4vpAXPwta4r6Ql4PwB3Rg9TN+7tyg465bA9VvTnLXrCea16tXL0nS2LFjHX6YR0VFKTU1VePHj1dWVpa9Hzxb/X0hExISFBcXp4ULF9r/sbBs2TKlpaUpOTmZJQItRPgDPAAbdGBkVVVVevzxx5Wfn6/w8HD98Y9/5LGgaID7QjoP4Q8A4HKlpaWSpKysLIcF/BkZGXr++ecb9APqff8fBu66XUx71uLHuwEA4CwtXZjPAn7Uq7810NChQx0eBzh06FAlJCQoJSXF3SW2G4Q/AIDLXXwvNl9fX4e2i9+35J5t6PguvjVQamqqoqKiZDab7WtEY2JiNH/+fJ7z20KEPwCAyy1YsMD++vz58w5tF7+/uB88V/2tgRYtWuSw01eSvL29tXDhQh09elSZmZluqrB9Yc0fAMDl/v3vf9tff3/N1sXvL+4Hz1V/y5/BgwfLarUqMzNTxcXFCg4OVnR0tAYPHuzQD80j/AEAXO7iRftTpkxRWFiYqqqq5Ovrq8OHD2vr1q0N+sFz1a/9/Mtf/qK1a9c2elPwi/uheYQ/AIDLhYeHKycnR9KF57LWhz1J6tOnj0M/IDo6Wr169dLChQsb3BT8qaee0qJFi9SrVy9u8txCrPkDALicj89/5x6++uorh7aL31/cD57t4llgm81m/w+tR/gDALjcxbN7zuiHji0zM1Nff/21li1bptzcXI0ZM0ZXXXWVxowZo88//1xLly5VaWkpGz5ayND/pDpw4IBWr16t/fv3q6amRmFhYbrnnnt0++23t2qcM2fO6KWXXtKHH36o4uJidevWTSEhIYqNjVViYmIbVQ8AaMrVV1/t1H7o2Oo3cjz88MNasGBBgw0flZWVWrRoERs+Wsiw4W/37t2677771KlTJ02ZMkV+fn7KyMjQ/PnzVVRUpNmzZ7donLy8PN1777369ttvNX78eN12222qrKxUQUGBPvzwQ8IfALhBS5/cwRM+IP13I0dubq6ioqI0YcIEh/bc3FyHfmieIcNfbW2tFi9eLC8vL23YsEGDBg2SJM2dO1d33nmnVq9erUmTJikkJKTZcSwWi+bMmSNJ2rRpk66//voG5wEAuN6WLVvsr729vR1u73Lx+y1btuhPf/qTy+uDsURHRyskJERLly5Vamqqw73+6urqtGzZMoWGhrLho4UMueZv165dOn78uGJiYuzBT5LMZrPmzJmj2traFj3GZePGjTp58qSSkpIaBD+JhcQA4C5nz56VJHXp0qXBur6+ffuqS5cuDv3g2Uwmk1asWKG0tDTFxcU5PN4tLi5OaWlpeuaZZ2QymdxdartgyPSzZ88eSdK4ceMatI0dO9ahT3PS09Pl5eWl2267TUeOHNHOnTtVXV2t/v37Kzo6Wp07d3Zu4QCAFrnqqqtUWlqq8+fPa/DgwVqwYIF8fX1VVVWl9PR0+33crrrqKvcWCsOIj49XcnKykpKSNGbMGPvx0NBQJScnKz4+3o3VtS+GDH/1f+j79evXoM3f31+BgYE6duxYs2N89913+vLLL3X11Vdr/fr1WrVqlcNlhT59+uj555/XwIEDnVo7YDShASZ1LT8s20lD/nFvM13LDys0gFkAVzpy5IjKy8tb1HfMmDE6fPiwJOn99993uM9f/axffb/6+wFeSkBAgPr379+KitHexMfHKzY2tsGGD2b8WseQPw0sFoskyc/Pr9F2s9mskpKSZsc4e/asrFarysvL9fzzz2vBggWKjY1VbW2t/vGPf+iFF17QQw89pPT0dIe/aICOxHT+rL6Y00OmD3+l2g/dXY1rhUv6fE4PHTrPZUNXOH36tMLDwxs8qq0lmnu272uvvabXXnutReOYTCaVlJQoKCio1TWg/TCZTA02fKB1DBn+nKH+LyCr1aq77rpL9957r71t3rx5Kiws1LZt27R9+3bFxsa6q0ygTVm7+GvQmjNK27RRkY2se+3I8g4eVMz0RCVP93d3KR4hKChI+fn5LZ75s1qtuu2221RWVtZkn6uvvlrbt29v8axOQEAAwQ9oAUOGP7PZLEmqqKhotN1isTQ5K1jv4vabbrqpQfvEiRO1bds25ebmEv7QoR0tt6o6IExePxjq7lJcqrqkVkfLre4uw6O09pLrunXrlJCQoC5duqi6utp+3NfXV9XV1Xr55Zc1cuRIZ5cJeDxD7vatv4VLY+v6zp49q7KyskbXA16sW7duuuaaayQ1vmC4/tj3LzcAAFyjfgH/tdde63D82muvZQE/0IYMGf7q/6WXlZXVoG3nzp2SpFGjRl1ynKioKEmyLyq+WP2x3r17X3adAIArEx8fr8OHD2vt2rWSpLVr1yo/P5/gB7QhQ4a/0aNHq0+fPkpLS1NeXp79uMVi0Zo1a+Tj46Np06bZj5eWlqqgoKDBZeI777xTkvTyyy/r22+/tR//+uuv9dprr8nb21u33nprG38aAEBzTCaTRowYIUkaMWIEOzeBNmbINX8+Pj5asmSJ7r//fiUmJiomJkZms1kZGRk6ceKEHnnkEYWGhtr7r1y5Ups3b9ayZcsc/rU4bNgwzZw5U3//+981depUTZw4UbW1tfrggw905swZPfroow7jAAAAdHSGDH/ShUu2Gzdu1KpVq5Senq6amhqFhYVp3rx5mjp1aovHeeKJJxQREaENGzZo8+bN8vLyUmRkpH73u9/ppz/9aRt+AgAAAOMxbPiTpKFDh2rdunWX7Ld8+XItX768yfb4+HjWjwAAAMiga/4AAADQNgh/AAAAHoTwBwAA4EEIfwAAAB6E8AcAAOBBCH8AAAAehPAHAADgQQx9nz8AzpGTk+O2c5eVlSkwMNDl57340ZAAgP8i/AEdWG1trSTpgQcecHMl7uPn5+fuEgDAUAh/QAc2atQo7d69Wz4+7vmjnpeXpxkzZmj9+vWKjIx0+fn9/PwUHh7u8vMCgJER/oAObtSoUe4uQZGRkRo2bJi7ywAAiA0fAAAAHoXwBwAA4EEIfwAAAB6E8AcAAOBBCH8AAAAehPAHAADgQQh/AAAAHoTwBwAA4EEIfwAAAB6EJ3x4oCNHjqi8vLzVX2e1WpWamipJ2rRpk6xWq0wmU6vHCQgIUP/+/Vv9dQAA4MoR/jzM6dOnFR4errq6uisaZ+nSpVq6dOllfa3JZFJJSYmCgoKuqAYAANB6hD8PExQUpPz8/FbN/O3YsUOPPfaYoqOjde+996pHjx46c+aM/va3vykzM1NPP/20brrpphaPFxAQQPADAMBNCH8eqDWXXK1Wq6ZPn66YmBi9+eabevHFF7V3714NGDBA7777rn7+859rzZo1+vWvf31Zl4ABAIBrEf7QrMzMTBUWFmrs2LHy8/NTbW2tvW3BggX6+c9/rqNHjyozM1MTJkxwX6EAAKBFCH9oVnFxsSRpw4YNuuaaa/SLX/xC/fv315EjR/T6669r48aNDv0AAICxEf7QrPq1ed27d1eXLl30zDPP2Nv69u2r7t2769y5c6zhAwCgneA+f2jWZ599Jkk6d+6cSktLHdpKS0t17tw5h34AAMDYCH9oVkFBgf21v7+/XnrpJZ08eVIvvfSS/P39G+0HAACMi8u+aJbNZpMkXXvtteratasefPBBe1toaKiuueYanTp1yt4PAAAYG+EPzaqf3Tt//rwKCwuVnZ2t4uJiBQcHa/To0QoODnboBwAAjI3LvmiWj8+Ffx+UlZWpX79+OnTokMaPH69Dhw6pX79+Kisrc+gHAACMjfCHZtXfu6937946ffq0Zs2apd69e2vWrFk6c+aMevfu7dAPAAAYG9M1aNaECRPUq1cvFRUVafLkyQoLC1N1dbW6du2qw4cPa9u2berVqxfhDwCAdoLwh2aZTCa98MILSkhI0Icffqht27bZ27p16yYvLy+98MILPNoN6CDy8/NVUVHh8vPm5eU5/OoOfn5+Cg8Pd9v5AVch/OGS4uPjlZycrKSkJBUWFtqPX3PNNXrmmWcUHx/vvuIAOE1+fr4iIiLcWsOMGTPcev5Dhw4RANHhEf7QIvHx8YqNjVVmZqZ9t290dDQzfkAHUj/jt379ekVGRrr8/GVlZQoMDHT5eaULM44zZsxwy6wn4GqEP7SYyWRibR/gASIjIzVs2DB3lwGgjRg6/B04cECrV6/W/v37VVNTo7CwMN1zzz26/fbbW/T1KSkpWrhwYZPtWVlZ6tmzp7PKBQAAMDzDhr/du3frvvvuU6dOnTRlyhT5+fkpIyND8+fPV1FRkWbPnt3isW6++eZGL2F069bNmSUDAAAYniHDX21trRYvXiwvLy9t2LBBgwYNkiTNnTtXd955p1avXq1JkyYpJCSkRePdcsstbEoAAACQQW/yvGvXLh0/flwxMTH24CdJZrNZc+bMUW1trVJSUtxYIQAAQPtkyJm/PXv2SJLGjRvXoG3s2LEOfVoiLy9Pf//731VbW6u+fftq7NixMpvNzikWAACgHTFk+Ku/l1y/fv0atPn7+yswMFDHjh1r8Xivvfaaw3s/Pz/95je/UWxs7BXVCQAA0N4YMvxZLBZJF0JaY8xms0pKSi45Tp8+ffR///d/GjdunHr16qUzZ87oo48+0nPPPafHH39cAQEBGj9+vFNrBwAAMDJDrvlzlpEjRyoxMVF9+/ZV165d1bt3b911111auXKlbDabVq1a5e4SAQAAXMqQ4a9+PV5Td1q3WCxNzgq2xLhx4xQcHKzPP/9c33333WWPAwAA0N4YMvzV38KlsXV9Z8+eVVlZWaPrAVsjMDBQNptN1dXVVzQOAABAe2LI8Ddy5EhJF57A8X07d+6UJI0aNeqyx7dYLDpy5IiuuuqqK5pBBAAAaG8MueFj9OjR6tOnj9LS0nT33Xfbn85hsVi0Zs0a+fj4aNq0afb+paWlqqioUK9evRzC3L59+zR8+HCHsaurq7V48WJVV1crLi5OXl5ervlQANAOhAaY1LX8sGwnDfnjoc10LT+s0ACTu8sAXMKQf7p9fHy0ZMkS3X///UpMTFRMTIzMZrMyMjJ04sQJPfLIIwoNDbX3X7lypTZv3qxly5Y5PMkjMTFRYWFh+uEPf2jf7Zudna3i4mJdf/31evTRR93x8QDAkEznz+qLOT1k+vBXqv3Q3dW4Vrikz+f00KHzZ91dCtDmDBn+JCkqKkobN27UqlWrlJ6erpqaGoWFhWnevHmaOnVqi8aYOXOm9u/fr6ysLH377bfq1KmTBgwYoLvuuku/+MUv1LVr1zb+FADQfli7+GvQmjNK27RRkddf7+5yXCrv4EHFTE9U8nR/d5cCtDnDhj9JGjp0qNatW3fJfsuXL9fy5csbHH/iiSfaoiwA6LCOlltVHRAmrx8MdXcpLlVdUquj5VZ3lwG4hCE3fAAAAKBtEP4AAAA8COEPAADAgxD+ALQJq9WqvXv3SpL27t0rq5X1VABgBIQ/AE6XkpKisLAwzZo1S5I0a9YshYWFKSUlxc2VAQAMvdsXgPsdOXJE5eXlLe6/Y8cOPfbYYxo3bpzGjx+vV199Vffcc4+OHDmihIQEPf3007rppptaNFZAQID69+9/uaUDABpB+APQpNOnTys8PFx1dXWt/trMzExlZmZKkl599VX78QULFrR4DJPJpJKSEgUFBbX6/ACAxhH+ADQpKChI+fn5LZ7527t3r2bNmiUvLy9FR0fr3nvvVY8ePXTmzBn97W9/U2Zmpmw2m9auXasRI0ZccryAgACCHwA4GeEPQLNac9k1Ly9PkjRp0iSlpaXJ2/u/y4p/8YtfKCYmRunp6erevbuGDRvm9FoBAJfGhg8ATvP1119LkuLj4x2CnyR5e3srLi7OoR8AwPUIfwCcpmfPnpIu7Pb9/jrBuro6paamOvQDALgel30BOE3v3r0lSenp6YqLi9PChQs1ePBg5ebmatmyZUpPT3foB2PKyclxy3nLysoUGBjolnPXL1kAPAHhD4DTREdHKyQkREFBQTpw4IDGjBljbwsJCdGIESN05swZRUdHu7FKNKW2tlaS9MADD7i5Evfx8/NzdwlAmyP8AXAak8mkFStWKCEhQVOmTNGCBQvk6+urqqoqbd++XVu3blVycrJMJpO7S0UjRo0apd27d8vHx/U/GvLy8jRjxgytX79ekZGRLj+/dCH4hYeHu+XcgCsR/gA4VXx8vJKTk5WUlKS0tDT78dDQUCUnJys+Pt6N1eFSRo0a5dbzR0ZGshMcaGOEPwBOFx8fr9jYWGVmZqq4uFjBwcGKjo5mxg8ADIDwB6BNmEwmTZgwwd1lAAC+h1u9AAAAeBDCHwAAgAch/AEAAHgQwh8AAIAHIfwBAAB4EMIfAACAByH8AQAAeBDCHwAAgAch/AEAAHgQwh8AAIAHIfwBAAB4EMIfAACAByH8AQAAeBDCHwAAgAch/AEAAHgQwh8AAIAHIfwBAAB4EMIfAACAByH8AQAAeBDCHwAAgAcxdPg7cOCAHnjgAY0cOVI33nijEhIS9M4771z2eN9++62io6M1cOBA3XfffU6sFAAAoH3wcXcBTdm9e7fuu+8+derUSVOmTJGfn58yMjI0f/58FRUVafbs2a0ec8mSJbJYLG1QLQAAQPtgyJm/2tpaLV68WF5eXtqwYYOWLFmixx9/XFu2bFF4eLhWr16twsLCVo25Y8cObdmyRY8++mjbFA0AANAOGDL87dq1S8ePH1dMTIwGDRpkP242mzVnzhzV1tYqJSWlxeOVl5frt7/9rWJiYjRx4sS2KBkAAKBdMGT427NnjyRp3LhxDdrGjh3r0Kcl/vCHP8hqtWrx4sXOKRAAAKCdMuSav/pLuv369WvQ5u/vr8DAQB07dqxFY7333ntKS0vTs88+q8DAQJ07d86ZpQIAALQrhpz5q9+U4efn12i72WxWRUXFJcf55ptv9H//93+65ZZbNHnyZKfWCAAA0B4ZMvw5y+9+9zvV1tbqySefdHcpAAAAhmDIy75ms1mSmpzds1gsTc4K1nv//fe1fft2LV++XD179nR6jQAAAO2RIWf+QkJCJKnRdX1nz55VWVlZo+sBL5aXlydJeuKJJzRw4ED7fzfffLMkKSsrSwMHDlRsbKxziwcAADAwQ878jRw5UmvXrlVWVpamTJni0LZz505J0qhRo5odY9CgQUpISGhwvLKyUtu2bdO1116rcePGKTg42HmFAwAAGJwhw9/o0aPVp08fpaWl6e6771ZkZKSkC5d716xZIx8fH02bNs3ev7S0VBUVFerVq5f9cvDNN99sn+W72IkTJ7Rt2zaFhYXpqaeecs0HAgAAMAhDXvb18fHRkiVLZLPZlJiYqN/85jf64x//qNjYWOXn5+vhhx9WaGiovf/KlSs1efJkvffee26sGgAAwPgMOfMnSVFRUdq4caNWrVql9PR01dTUKCwsTPPmzdPUqVPdXR4AAEC75GWz2WzuLsKIcnJyNHz4cO3bt0/Dhg1zdzkA0KHxdy7gOoa87AsA8BxVVVVavny5JGn58uWqqqpyc0VAx0b4AwC4TVxcnLp166a33npLkvTWW2+pW7duiouLc3NlQMdF+AMAuEVcXJy2bNnSaNuWLVsIgEAbMeyGDwBA+3LkyBGVl5e3qG9VVZU9+HXu3Fnfffedva3+/ZYtW7Rz5075+vq2aMyAgAD179+/9YUDHobwBwC4YqdPn1Z4eLjq6upa/bUXB7/vvx83blyLxzGZTCopKVFQUFCrawA8CeEPAHDFgoKClJ+f3+KZv7vuuksHDx5U7969ZbVaVVJSYm+79tpr5e3trZMnT+r666/Xhg0bWjRmQEAAwQ9oAcIfAMApWnPJtXv37pKkoqKiBm0XB8Hu3btz6xfAydjwAQBwuZtuusn+ulOnTkpMTNTKlSuVmJioTp06NdoPgHMw8wcAcLmrr77a/rqmpkYbN27Uxo0bm+0HwDmY+QMAuNwHH3zg1H4AWo7wBwBwubKyMvtrLy8vh7aL31/cD4BzEP4AAC7Xu3dv++uuXbs6tF38/uJ+AJyDNX8AAJeLi4vT22+/LUmaMGGCpkyZIl9fX1VVVWnr1q1KT0+39wPgXIQ/AIDLnT171v46PT3dHvaa6wfAObjsCwBwuZ49ezr82tp2AJePmT8AgMvVr+X7+uuv9bOf/UwWi0VnzpxRjx49ZDab7TOBrPkDnI/wBwBwuejoaIWEhMhkMikjI0NWq9XeZjKZNGDAANXV1Sk6OtqNVQIdE5d9AQAuZzKZdMcdd6igoEBBQUFKSkrS888/r6SkJAUFBamgoEAJCQkymUzuLhXocLxsNpvN3UUYUU5OjoYPH659+/bxXEkAcDKr1aqwsDAFBQXp9OnTKiwstLeFhoaqR48eOnPmjPLz8wmAgJNx2RcA4HKZmZkqLCzUG2+8oZEjRyozM1PFxcUKDg5WdHS09uzZozFjxigzM1MTJkxwd7lAh0L4AwC4XHFxsSRp8ODBMplMDQLe4MGDHfoBcB7W/AEAXC44OFiSlJub22h7/fH6fgCch/AHAHC5+t2+S5cuVV1dnUNbXV2dli1bptDQUHb7Am2A8AcAcDmTyaQVK1YoLS1NcXFxys7OVkVFhbKzsxUXF6e0tDQ988wzbPYA2gBr/gAAbhEfH6/k5GQlJSVpzJgx9uOhoaFKTk5WfHy8G6sDOi7CHwDAbeLj4xUbG9tgty8zfkDbIfwBANyqsd2+ANoOa/4AAAA8COEPAADAgxD+AAAAPAjhDwAAwIMQ/gAAADwI4Q8AAMCDEP4AAAA8COEPAADAg3CT5yZUVVVJkvLy8txcCQAAQMtcf/316tatW7N9CH9NKCwslCTNmDHDvYUAAAC00L59+zRs2LBm+3jZbDabi+ppV06fPq13331XISEh8vX1dXc5AAAAl9SSmT/CHwAAgAdhwwcAAIAHIfwBAAB4EMIfAACAByH8AQAAeBBu9dJBnTp1Sunp6fr444915MgRnT59Wv7+/ho2bJjuv/9+3XDDDQ2+xmKxaPXq1crIyNDXX3+tnj176tZbb9X//u//ymw2O/Q9c+aMkpOT9fnnnys3N1dFRUWSpC+//LLRegoLC7V9+3ZlZmbq2LFjKi8vV48ePfTjH/9Ys2bN0oABA5z/m4AWMdr3SlFRkdauXavc3FydPHlSFotFPXv2VHh4uGbOnKnRo0c7/zcBLWa075fGPPnkk3rjjTckSVlZWerZs+cVfGJcLiN+rwwcOLDJtqSkJD344IOX+WnbF3b7dlDPPPOMXn75ZfXt21cjR45Ujx49dOzYMb3//vuy2WxasWKFJk+ebO9fWVmpxMRE5eXlaezYsRo0aJAOHjyozMxMRUZGauPGjQ5bx3fv3q27775bXl5e6tevn06dOqWqqqom/9D9+te/1rZt2xQREaFhw4bJbDbr0KFD+vjjj9W1a1f99a9/1YgRI9r89wUNGe175ZNPPtG8efN044036rrrrpPZbNapU6f0wQcfyGKx6Ne//rVmz57d5r8vaJzRvl++Lzs7WzNnzpSvr68qKysJf25kxO+VgQMHqnfv3po2bVqDttGjR3vOzyEbOqR3333X9u9//7vB8X//+9+2H/7wh7ZRo0bZzp8/bz/+3HPP2SIiImxPP/20Q//6488995zD8a+//tq2Z88eW0VFhc1ms9luu+02W0RERJP1bNq0yZaXl9fgeFpami0iIsI2efLkVn0+OI/RvlfOnz9vs1qtDY6XlJTYxowZY/vhD39oO3v2bKs+I5zHaN8vF6uoqLBNnDjR9vDDD9tmzJhhi4iIsJWWlrb2I8JJjPi9EhERYZsxY8blfqQOg/Dnge69915bRESE7cCBAzabzWarq6uzjRs3znbjjTfazp0759C3urraNnLkSFt0dLStrq6uyTFb8xf099166622iIgI25kzZy7r69F2jPa9MnfuXFtERESj/5CA+7n7+2Xx4sW2UaNG2b7++mvCn8G563uF8HcBGz48kI+Pj8OvhYWFKi0t1bBhwxrcFbxLly4aMWKETp06pWPHjrVJPZ06dXKoB8ZhpO+VsrIyffrpp/L19dV1113n9PFx5dz5/ZKVlaU333xTCxcuVFBQ0BWPh7blzu+Vb7/9Vm+99ZZefPFFvfnmm/bHuXoSftp6mJMnT+qTTz5Rz549FRERIUn2P0whISGNfk2/fv3s/Zrqc7kOHDig/Px8DRkyRFdddZVTx8aVcbYHgxwAAAVWSURBVPf3yokTJ7R582ZZrVaVlpZqx44dqqio0JNPPtlg4Tfcz53fLxaLRYsXL9b48eMVFxd32ePANdz9d8vBgwe1ePFi+3svLy/dfvvt+v3vf+8xj3Ml/HmQmpoaPfbYY/ruu+80f/58mUwmSVJFRYUkNfkDtf54fT9nqaio0OOPPy5vb28tWLDAqWPjyhjhe6WoqEh/+ctf7O+7deumpUuXKjY29orHhnO5+/vlqaeeUkVFxf/X3r2EQvfGcQD/elO85DJym5WRDSk7DSEb5Zoi47IgJqVQSjGuCzIRxQILVqIsZ0STJqSRhVymXJpIWSiLcSsRjXF5Fxp//v/hff8zOOd1vp+yMOf0PGem7zn9zjzPmQednZ1utUOfT+isqNVqZGZmIiIiAh4eHrBYLBgYGMD09DTu7+/R39/vVvt/Cw77SsTDwwNaWlqwtraGwsJCwe+ObTYbamtrcXBwgLq6OiiVSkGPh/4hlqwolUrs7e1hZ2cHRqMRxcXF0Gg06OrqEuR4yDmh82IymaDT6dDQ0IDw8PAv7Zv+H6GzAgAajQZxcXEICAiAv78/EhISMDY2hoiICBgMBuzv73/5MQmBxZ8EPD4+oq2tDdPT08jNzUVHR8er7X5+fgCehk6ccbzu2M9dt7e3qKmpwcrKCqqqqvizHSIitqwAT3NCFQoFNBoNSkpKMDExAZPJ9GHtk+uEzsvNzQ3a29uhVCpRVFTkUhv0NYTOynt+/vyJ7OxsAIDZbP7w9sWIw77f3MPDA1pbW6HT6ZCTk4Oenh78+PG65nfMpXhr0qtjLoZjP3fYbDZUV1djeXkZlZWVqK+vd7tN+hhiy4ozSUlJmJycxOrqKlJTUz+lD/ozYsjL2dkZrFYrrFYroqOjne6TnJwMAJiamkJMTIxL/ZB7xJCV35HJZACebiikgMXfN/byhMvKykJvb+/z/IqXFAoFQkNDYTabcX19/epJK5vNhvX1dYSGhrp90r0s/NRqNef5iYjYsvKW4+NjAHwyXGhiyYuvry8KCgqcbjOZTDg5OUFOTg68vb0RGBjoUh/kHrFk5Xc2NzcBQDK/JMBh32/q5QmXkZGBvr4+pycc8PSkk0qlwvX1NYaHh19tGxkZwcXFBVQqFTw8PFw+npeFX0VFBTQajctt0ccSW1a2trZgs9n+8/rR0RFGR0cBACkpKS63T+4RU15kMhm0Wq3Tv8jISABAU1MTtFot5HK5S32Q68SUFQCwWCxOv9mbnZ2FwWCATCaTzPKRXN7tmxocHMTQ0BB8fHxQVlbm9JuStLS052GQfy+rExsbi93dXSwtLTldVgd4uqg6zM3N4erq6tWSOY2NjQgKCnreV6/XIyQk5M25OXl5eZK56xITsWWluroaGxsbiI+Ph1wuh6enJw4PD2EymWC321FeXo7m5ubP+CjoD4gtL28pLS3F6uoql3cTkNiy0tTUhPn5eSQmJkIul+Px8REWiwXr6+vw8vLC4OCgZKaTsPj7phzF1nu6u7uRn5///P/l5SWGhoZgNBpxenqK4OBgpKeno7a21ukk2/cWyAaAhYWF52LOcSF+z/j4OJ/6FYDYsrK4uIiZmRlsb2/j9PQUdrsdQUFBiIuLg0qlkszFWazElpe3sPgTntiyMjc3B71ej93dXZyfn+Pu7g5hYWFISEiAWq1GVFSUC+/y78Tij4iIiEhCOOePiIiISEJY/BERERFJCIs/IiIiIglh8UdEREQkISz+iIiIiCSExR8RERGRhLD4IyIiIpIQFn9EREREEsLij4iIiEhCWPwRERERSQiLPyIiIiIJYfFHREREJCG/AKQf0Mz4MDpmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boxplot(samples[:, :4], labels=[y for y in range(2012, 2016)])\n",
    "plt.ylabel(\"λ\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "trace = {k: v for k, v in zip((\"λ1\", \"λ2\", \"λ3\", \"λ4\", \"γ\"), samples.T)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLsAAADTCAYAAABp7hHfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XeYG+XVNvB71LbvuuHecMGADcYVjOktQHhDMCWBFEIJIXQSEiCB0EJCINRAgIQPyJuEN5AEQgkYbIzt9br3bu+6rtf22t71rrZK074/xvKOpJE0I400kvb+XRfXZbZIz4xmZuc5c855BFVVVRAREREREREREeUBl9MDICIiIiIiIiIisguDXURERERERERElDcY7CIiIiIiIiIiorzBYBcREREREREREeUNBruIiIiIiIiIiChvMNhFRERERERERER5g8EuIiIiIiIiIiLKGwx2ERERERERERFR3mCwi4iIiIiIiIiI8gaDXURERERERERElDcY7CIiIiIiIiIiorzBYBcREREREREREeUNBruIiIiIiIiIiChvMNhFRERERERERER5g8EuIsq4Bx54ABMmTMAdd9yBpqammD/37W9/G2PGjMFHH32UwdERERER5S7eZxERAR6nB0BE3c/w4cNRWFiIWbNmQVVVvPLKK1E/s3PnTqxatQqlpaW46KKLHBglERERUe7hfRYRETO7iMgBt956K2bOnIl+/frhyy+/RH19fdTPvP/++wCAr3/96ygsLMz0EImIiIhyEu+ziIgY7CIih1RUVOCaa66BqqqYOXNm2PcURcGHH34IAJgxY4YTwyMiIiLKWbzPIqLujsEuInLMxRdfDABRN2FVVVXYv38/RowYgVNOOcWJoRERERHlNN5nEVF3xmAXETlm1KhRGDlyJFatWhWWYv/BBx8AAK644gqnhkZERESU03ifRUTdGYNdROSoiy66KCzF3u/3Y/bs2XC73bj88ssdHh0RERFR7uJ9FhF1Vwx2EZGjIlPsP/nkEwQCAZxxxhno16+fk0MjIiIiymm8zyKi7orBLiJy1PHHH4/hw4cfTbEPpdazYSoRERFRanifRUTdFYNdROS4UIr9K6+8grVr16JHjx4477zznB4WERERUc7jfRYRdUcMdhGR4772ta8BAN59910AwGWXXQafz+fkkIiIiIjyAu+ziKg7YrCLiBw3btw4DBo06Oj/M7WeiIiIyB68zyKi7ojBLiLKCqGnjmPGjMHYsWMdHg0RERFR/uB9FhF1Nwx2EVFWGD16NADgxBNPdHgkRERERPmF91lE1N0w2EVERERERERERHmDwS4iIiIiIiIiIsobDHYREREREREREVHeYLCLiIiIiIiIiIjyhqCqqur0IIiIiIiIiIiIiOzAzC4iIiIiIiIiIsobpoJd7e3tWLlyJdrb29M9HiIiIqJuhfdZRERERPYyFezavHkzJk2ahM2bN6d7PERERETdCu+ziIiIiOzFMkYiIiIiIiIiIsobDHYREREREREREVHeYLCLiIiIiIiIiIjyBoNdRERERERERESUNzxOD4BSI8sqamsVbKmWsG2bDH+LivZ2FYoMVFQI6NHDhaFDXTh+jAdDBrvgcglOD5mIiIiIiHKMKKpYvFRE/74ujB7NaSQRZTdepXKQLKtYvETErNlBfDkniMbDqqnfKy0VcNqpXpx1hhdnn+1Dn95M7CMiIiIiosRq9yjw+1X4/TKDXUSU9XiVyiGiqOKTTwN44/91YOcu5ejXCwuBUSPdOG60B716CSguFiAIQHOzisZGBdt2yKiultHaqmL2l0HM/jIIt7sNZ57hxYxvFuLss7zweJjxRURERERExhTZ3AN2IqJswGBXjqhaGMQTT7ahdo8W5CovE3Dx13y48AIfpkz2wuuNH6ySJBWbNsuoXBDE3PlBbNggY+48EXPniRg00IXrv1+IGd8sRFERg15ERERERERElLsY7MpyTc0Knv59Oz78KAAA6N1LwA+uL8K3ri5ESYn5wJTHI+CkcR6cNM6D224txrbtEj74MID/fBhA3V4Fv3mqHa+93oFbbynC1VcXwpcgeEZERERERERElI3YtCmLbdok4epvN+PDjwIQBOC71xXis//2xI0/KLIU6DIycoQH991bgtkze+LhX5ZgyGAXGg+r+M3v2vE/lzfh81kBqCpTlYmIiIiIiIgotzDYlaU++TSA7/6gGXv3Khg61IW//285Hry/BCXF9mZcFRYK+PY1hfj4Pz3wyMMl6NNHwJ46BT+5rxU/uq0Fu3bLtr4fEREREREREVE6MdiVhd54swP3P9iKzk7grDO9ePfvFRh/sjet7+n1CrjmqkJ89nFP/PhHRfB6gaqFIr55ZRPe+ksHZDakJCIiIiIiIqIcwGBXFlFVFa+82o7nX2wHANx0QyFefrEM5eWZ+5iKiwXccVsxPvx3D5x+mhfBIPD759pxw81+1O5hlhcRERERERERZTcGu7LISy934I+vdQAA7rmrGD+5pwRutzON4ocNc+NPr5XhsUdKUFwMrFgp4cprmvHZzIAj4yEiIiIiIiIiMoPBrizxv3/rwJ/e0AJdD/y8GD+8qcjhEQGCIOCqGYX4z796YOIED9raVNx3fysef7IVgQDLGomIiIiIiIgo+zDYlQU+mxnA757RShfvvbsY3/uO84EuvUGD3HjrjXLccrM2rnffC+AHN/lx4IDi8MiIiIiIiIiIiMIx2OWwFStFPPhQKwDgO9cW4qYbCh0ekTGPR8Dddxbj9T+WobxcwNp1Eq65rglr14pOD42IiIiIiIiI6CgGuxx04ICCe+9rgSgCF17gw/0/K4YgONOjy6wzpvvw7jsVGDXSjYMHVVx/kx9fzGYfLyIiIiIiIiLKDgx2OSQoqrj3Zy1oaFBx3Gg3fvNEqWPN6K0aOsSNd/5agbPP0lZr/Ml9rfjLXzucHhYREREREREREYNdTnnm9+1YvVpCWZmAF54rQ3FxbgS6QkpKBPzhhTJ8+1sFUFXg6d+34/fPt0FV2bieiIiIiIiIiJzDYJcD5nwVxDv/6AQAPPVkKYYNdTs8ouS43QIeerAEP7mnGADw1tudeOrpdga8iIiIiIiIiMgxDHZl2MFDCn71mNaQ/obrC3HO2T6HR5QaQRBw0w1FeOShEgDA397pxBNPtkFRGPAiIiIiIiIiosxjsCuDVFXFrx5txeHDKsaMceOuO4qdHpJtrrm6EL9+rASCALz7zwAeeawNssyAFxERERERERFlFoNdGfTePwOYXynC5wN+92QpfL7c6tOVyBXfLMRvf10Klwt4/z8BPPSrVga8iIiIiIiIiCijGOzKkH37Zfz++TYAwD13FWP0aI/DI0qP/7msAE//thRuN/DRJ0H84qFWljQSERERERERUcYw2JUBqqriiV+3ob0dOGW8B9+9rtDpIaXVJRcX4NmnS+HxAJ98GsRvf8em9URERERERESUGQx2ZcBnM4OYVynC6wUee6QEbnd+lS8aufCCAjz5RCkA4J1/dOLV1zscHhERERERERERdQcMdqVZU7OC3/5OK1+85eYijBqZn+WLRi67tAC/uF9rwv/Kqx34v3c7HR4REREREREREeU7BrvS7KWX29F4WMXIEW7cfFOR08PJuO9cV4TbbtW2+8nftuHTzwIOj4iIiIiIiIiI8hmDXWm0cZOE9/6pBXce/mUJfN78L180ctutRbj2WwVQVeDBh1qxoCro9JCIiIiIiIiIKE8x2JUmiqLiyd+2QVWBSy/xYcpkr9NDcowgCPjFAyW49BIfJAm456ct2LRZcnpYRERERERERJSHGOxKk48+CWD1GglFRcB99xY7PRzHuVwCnnyiFNNO86KjA7j9zhYcOKA4PSwiIiIiIqKsdqhBwcZNEmSZK9wTmcVgVxq0tat4/sV2AMCPf1SMfv3cDo8oO/i8Ap57phQjRrhRf0DBHXf70d7OCzYREREREVEsS5aK2LFTxo4dstNDIcoZDHalwZtvdeDQIRVDhrjwve8WOj2crFJe7sIf/1CGnj0FbNgo44FftEBRGPAiIiIiIiKKh4kCROYx2GWz/fUy3v7fDgDAT+4p7rZN6eMZMtiNl54vg9cLfPmVeDQLjoiIiIiIiIgoVQx22eyllzvQ2QlMnODBhef7nB5O1po4wYtfP14KAHjz7U786/1Oh0dERERERESUxZhHQWQag1022rRJwkcfBwAAP/tpMQSBV6N4Lru0ALfdWgQAeOLJNixfITo8IiIiIiIiIiLKdR6nB5BPnn+pHaoKXHKxDyef5HV6ODnhtluLsGOnjM9mBnHvT1vw7v9VYOAANvQnIqLcdN5556Gurs7we9/61rfw+OOPZ3hERERERN0Pg102WbJURNVCER4PcNcdxU4PJ2cIgoAnHi3Fjh3N2LxFxt33tuB/36pAURGz4oiIKDeVlZXh+uuvj/r6uHHjHBgNERERUffDYJcNVFXFCy9pTdavvrIAQ4cwM8mKoiIBf3ihDNdc14yNm2Q8+ngrnvpNKctAiYgoJ5WXl+POO+90ehhERERE3RZ7dtlgzlci1q6TUFQI/OgWZnUlY+BAN557pgxuN/DJp0H85X/ZsJ6IiIjST1VVp4dARERENmOwK0WKouLFP2hZXd/7XhGO6cNdmqypU7y4/z4tWPjsC+2oWhh0eERERETWBYNBfPDBB3jttdfwzjvvYPPmzU4PiWLo6FAx+8sgNm+RnB4KERER2YhljCmaNTuIbdtllJcJuOH7hU4PJ+ddd20hNm2W8cGHAdx3fyvefaeCZaFERJRTDh48iAceeCDsa2eeeSaefvpp9OrVy6FRkZGabTKCIrBtu4zjx/C2mIiIKF8wDSkFiqLitT91AAC+991ClJdzd6ZKEAQ8/MsSnDTOA79fxV33tKCtneUFRESUG2bMmIG//vWvWLRoEVasWIH33nsPZ511FiorK3HbbbexZI6IiIgoAxidScFXc0VsrZZRUiLgO9cxq8suBQUCXnq+DH36CKiukfGLh1qhKJwcEBFR9rvjjjswdepU9OrVC6WlpRg/fjxef/11TJo0CatWrcK8efOcHiIRERFR3mOwK0mqquLV17VeXd+5rhAVzOqyVd++Lrz4XBk8HmD2l0G8/ucOp4dERESUFJfLhRkzZgAAVq5c6fBoiChZoqiiulpCa6vi9FCIiCgBRmiSNL9SxKbNMoqKgO9/h1ld6XDKeC8e/mUJAODlP3ZgzldsWE9ERLmpZ8+eAICODj68IcpVmzZL2FojY16l6PRQqJsSBMHpIRDlDAa7kqBldWk3q9d+qxA9e3I3pstVMwpx7bcKAAD3/6IFNTVcLYmIiHLP2rVrAQCDBg1yeCRElKymJrbVSIUsq+js5D4kosxglCYJCxeJWLdeQmEh8IPvFzk9nLx3/89KMGWyB+3twB33tKCpmanjRESUfWpqauD3+6O+vnz5crz11lvw+Xy46KKLHBgZERnx+xXU1cmmf57rS6Rm7jwRX34VRFsbdyQRpR/XWLZIVVX88UhW1zVXFaJ3b8YL083rFfD878twzXXNqK1V8LOft+LVV8rg8TCNl4iIssdnn32GN954A9OmTcOgQYPg8/mwdetWVFVVweVy4bHHHsPAgQOdHibpsCKoe6us0soRvV4Bfftau6c/fFhhdYdFnQEtyHXwoIKSErfDoyEntbWpUFUVpaU8hyh9eHRZtGSphNWrJfh8wI0/YFZXpvTs6cIfXihDUSGwcLGI515od3pIREREYU499VSce+652L59Oz744AP89a9/RU1NDS699FL84x//wNVXX+30EClCNmXqNDYqzHhxiN9vvWogm44dolyiKCrmzg9iXqUIWeaJFKmtXUXVwiD27jOfdUrGmNll0Wt/0oIsV80oxDHHMFaYSceP8eDJX5fiJ/e14i9/7cSYMR5c/j8FTg+LiIgIADB16lRMnTrV6WFQDmptVbBoiZZl9PVLeG9D2YkBvu5LUVS4XPakwsq6GI4oAm4m+YVZv15CU7OKVaslDBzAnZMKRmssWLFSxLLlEjwe4MYbuAKjE752YQFuvUXLqHv08VasXcfVcIiIiCg52VLG6G9JTxSBWRO5qaFBQXMS2WZ2qK+XsX1H7mWU7K+X0dbO4z0dmpoVfPZ5EJu32L9QGAOo0USRO8UuDHZZ8NqRXl0zvlmAAf0ZZXXK7T8uwrnneBEMAnff24KDB9mwnoiIqLvavEXCwkVBBnYi7NgpY+YXQeyvz73ARS5I1yS9s1PF4qUiFlR1PdBtbVWwbr2Ejo70H+PLV0rYtFnC4cO5c39dXy9jxUoJc+cFnR5KXtq8WbuGbNvOawnlFga7TFqzVsTCxSI8HuDmG9mry0kul4CnnizFyBFuHDio4q57WxAI8AaXiIioO9q2XcbhJhV79+XO5DxKGm5jNm7SsjDWrLE/GyNdDhxQ0NnZve/pOg3uaasWidhdK2PZchE1NVJGViY3Gke2amzMnbGmKluyUYlyAYNdJr32Jy2r6xuXFWDQIGZ1Oa201IWXXypDebmAteskPPp4K1TmwRIREXVbvA3IbXv3yli2QsRXc7t5do7BcSwdiVe2tKrYUi2jaiHbeFDmcI4VjfskNzDYZcL6DRLmV4pwuYCbb2JWV7YYOsSN554phdsNfPRJEK8eKTMlIiIi++3aLaOujmUslB4HD2nZSkqOzCGZYUPdRbbFdVpaFLQ72J+tpUXB518EsbU6d7JmuysGu0x4/c/aCoyXXerDsKHM6som007z4eFflgAAXnm1Ax9/EnB4RERERPmns1PF+g0SVq+V+ESbKJ0YRLOMV6T0yqZLfiCgYv4CEV852J9t8xYZsgJU1/DhT7ZjsCuBzVskzPlKhCAAt/yw2OnhkIGrryzEjT/QVsd8+NFWrFjJ1G4iIqIQUVShpJguI0lZNNshaxg8SYtsCgAQpVM2Heutrc4PJt37g1mj9mGwK4HX/6yVxl38NR+OHc6srmx1793FuPACH0QRuPOeFuzaxUg7ERFRUFTxxewg5s7jg6BUtLWpaG3N4Qb4FrS1q5hfGUTtnuTupTo7VVRXS1w8iFLW2qpga7UEUczeYykYVLFwURC7a/N37pFNez8bypyzKfhH8THYFUd1tYQvZmkpkrf+kL26sllohcaTT/KguVnFrXf40dTUPW5KiYiIYmls1P4Wdti4wl13u9FXFBVz5wcxr1KELOfgxlsc8saNElpaVaxdl1w/mmXLRWytkbF8BQOslJp5lSKqa2Rs2pS9vZG2Vmurwa5bn71jzDYpXUUduARHlu53t7+BuYzBrjhef0PL6rrwAh9GjfI4PBpKpLBQwB9eLMOggS7s3q3g9rtb0NHBqxERERElT9IlbIjdIH4jp5ig4m/R7r2amrP7HiyZbKGMTHJ175HP/fG2VktYsVI0tY2Nh7N3P3SLEu8s2sRUS/Kt2rRZS35ps6EhflubipYWJmNkEoNdMWzfIWPm51pW14+Y1ZUz+vR24Y8vl6G8TMDq1RJ++vOWrE59JiIiygWqavxvirZpk4TlK8xN4kPSuk/Z/8VQsplrIeyrk5rqGhn76xUcakh88Cc6P6ycP52dKnbtlnMzS9Mh+v3b2Kgg6ODcKtNljNt3yJBkoKam63qRTBBaVbUM4fkLRM5NM4jBrhj+/EY7VBU49xwvTjieWV25ZNRID175QxkKCoB580U88nhbXj8ZIyIiStbevTIqq+x5ak2a7Ttl1B9QcDiLs1EI2F9vPcMiXtBXlrUgCqsKrFFMBJ3s3KMLF4lYv0HC5i3522PLbvpjfdESEQsWOJfiqupOW6fmd8m8raIbd7AbZAhnCwa7DOzaLeO/nx3p1XULV2DMRRMnePHcM2Vwu4EPPwrg2RfanR4SERFRxiVKPlm1RoLfr2I9+83YzigDob5exoqVIoLB7AmIyLKKpiYFqqqt2tnQyDKbZGytlrF+g4TKqmBqL8SMsWg2ni6h/oUHDthznHfH5+kp94BM4ded2t9Ws5sPNSixHyJ1w2PGKQx2GXjj/3VAloEzz/Bi3FhmdeWqc8724bFHSgAAb73dibf/0uHwiIiIiFInyyqamu0NSFgp6el2kzvd9qZaurZ8pYT99Qo2b86e4OKy5RKqFonYtUvBps3Mdom0u1bGfz8LJFxt79Ah7Zy0s69btzvXYkhHn6Zc3bVOlM9m03GomAg6tbend/XcRPujqUnBkqUi5s5LMfBNKWOwK0JdnYyPPgkAAG69hb26ct0VlxfiJ/do2XnPPNeODz8OODwiIiKi1CxbLqFqoYhdu9MbmNi3Xw7rU0KJmc3Y6syi25FQJteu3TJ27mKwK1Jolb1Eq+05ERBQVRVbqyUcasjvbDw7e3aRddm0f2XdaRhrXF/N01bPTVcGbaJXbWpK8BMJApa50g9w23YJC6qCjvZwS4TBrghvvNkBSQKmnebFKeO9Tg+HbHDjDwpx/fcKAQAPP9KKz2dl0R0mERGRRaHgxO40B7tWrpKwpVrGoQaFDepN2rCRwcF8lpFD38Kb7KlTUF0jY8nS/G4ClM2XnO5wPUznJlrdf50B87/QbrIXZVu7ivUbpLi9K8PGmcQOSfX37aCqKg41KLYtzrB5i4xmv4qdO7P3IQmDXTr79st4/z9aIOTHP2JWV74QBAH3/aQYV1xeAFkGfv5AK+Z8xbRSIiIiM9pau8FsLg79JCXRE/fWHN5XTk/a7UhmUBQ1vQ3is+zjbWvLzICcXugpLW+f4U1atDiIxUtyNChp876K93lKktZDMBZ9tpZdGX9LlmqZ0suXx/589u5TsHataOl1s82WrVpgfNVqex/KKFmcWMpgl86bb3VCkoApkz2YNJFZXfnE5RLw2CMluOxSHyQJuPe+FsyvZMCLiIgoES6T3iU8wy16v+TqJCgXiaIaVT6zoErEnLlBNPszM/vq6FBty5Iwkg3H0549MmZ9GUSjjQsX6D83U5uYBfshFZ2dKhoPa4s/ZHPJVyzpPA4jX3rxEhFVi8SEPfIA+4JdoQB5a4LgcW2dgqZmJeHrGj0UyYZzedeRMvX6JBZnaGtXsXev7Hjg2yoGu444eFDBv97vBMAVGPOV2y3gySdK8bWLtIDX3T9pwcJFDHgRERHFE9lwO7dudZMjyyoWLgpi8xbJcJJSXS3hi1lB25og59j8wXGqquKL2UHMmh0MCza1HMms27cv/cGu1lYVc+YGMb8ytWydhgYFS5eJKWdpfTkniLo6+8uJ1qyTIIrAoiUitu+w5/Wra6y9TjZnjlg9d3OkHZNjmv3aDrXjWE7HdVVVcvd6ncq4584LYtUaCXv3ZvHJaIDBriPefLsDwSBwyikenDqVKzDmK49HwO9+U4rzz/UiGATuvKcFy+KkrBIREWWrXL3hTqeODhV+G7J69u9XcLhJxbbtxhOurTUyJFnrWZIv4vWrSadAQLW82p4++OFUs/99+7VBtB/JCkn2fFy8VMTBQwpWrj5yP5pkNKQzoGL12vT2jJMcOtzTErTI4OGe69dqu7N57Hq5TZuluJmVyb5NW7uK3bWy4XVJEJIbf44fAkc1HlbR2akiYKF3mpMY7IL2ROW9f2lZXT++pQhCriyBQEnxegX8/ukynHWmF52dwI9v9+d9Y08iIqJkRd3SZvE97py5QVRWiejsTG2Qsi6Ywub86dPSomD2nCCqFubGfVg6P//OUK+xHDveJElNKRgSNFFkkejVk3n7zoC2kmUira0KamokR8u5nQwsNDQoCNhUCNPWrkaXGMfYNP18PFZAq3aPgu0xHkgAyZ8/8+YFsW69hG3bDF67m4cJZFnFl18FMXtOblRHMdgFLaursxMYN9aN6aezV1d34PMJeOHZMkw/3YuOTuDHd/hRWZUbJy0REVEycqnXxqrVIiqrgkmPucWhRvE5tIuzQig7yt+Swo5zcGWz7i4YVPH5rCAqFyQfrFy3XrKc2WeX6hrj7B0A8PsVLKgKYl6liC3VckpZnKkeKlaDwcGg9V5yoqgaBvSWLrMvED13XhALqkQEgubHdqhBwcwvgtiy1TgwGdlnK+y8THK/h36toTFGZleiF8jSgJgdl6zOztRfI5O6fbDr4EEF//eu9qndcVsxs7q6kYICAS+/WIZzzvYiEADuvLuFqzQSEVG3FO8mWFUzH0/Yu0+B36+i0WCyYYadt3PJThCs/F6m9m+zX8GePflTekn2shLAq9km4VCDFqw0Ci7LsoqGBsXUawYcKkWNZ/kK6Wj/KABxVwhMJNXzu8NCpmogoGLWl9YybxRF64H3xezoBwzpiENa6U+3caMW5KoxyrJKwMzhLEnWNlAAkvtAY/xOOhe4ALTs2Q0bJa1cvBvG57t9sOuNNzsQCADjT/bgjOnM6upufD4Bzz9bhosu9EEUtVUaZ36ehX9xiYgo71nNFGlpjV6NLv7rWx1Rel8nnex8dKnfXCe2fdMmyVS5lRkLqkSsWSfh4KH0NBm2WuqVC8eSWU70lUrrKnkJXnvL1vjBh5WrJCxeKmJrdW4GV4MWso+sSPchf/hIUE6ycMnQl5KK6W37FiXR372ob1vYgZKsYvUaEXv3GR+Dhw8r+HyWtUQHqw9S4mUt7t0nY+YXQezcmb5zZP4CETt3yVi3PsMfbJbo1sGu/fUy3v2nltV15+3s1dVd+bwCnnmqFP/zdW2Vxp890IoPP2bAi4iIMiuZMpkVK+zvdbR3n4xt27tujFU1fZPqvftkzP4ygMOHs2eFJ1VVsX+/vmlX+PetBBgBoKZGQq2FbKqOjq7X7+xUsX2nHLfcyoxAQA1b3aw1TWWetXsy8zk6FSTTT8xzMVDX2qpgxUoxum9SGhw4qL3Hrt2Jj/1U92UqJaWmfzWVaWKWHSsdHVq/slAvsLApsMmxhvb57loZS5eJcTOkAgEVNduSC7ZY/Wj1P797t4K6vQpWrTZ+b6urggLWG9Tv3KVEjSskNK4Nm9IfiEolM1HPKFySzSGUbh3s+tMbHRBFYPIkD047lVld3ZnHI+DJJ0px5YwCKArwy4db8be/dzg9LCIi6kZCJUFWNB62N7NLVVWsWi1lbJXBVaslBILAMl3QLnLi2tZ5sHpUAAAgAElEQVSmBWqsZACkcvNdV6fEzHxqalYwa3bsTIDIIba0KNhSLWPtOvOTmcVLuvaFvlF+KoGkqoVi2Ep9Uoz+PJEkyXwpGhA/i8FO4YsGZFkk4QhJUrG/Xk57mZIVy5ZL2F+voKrKOEjuWBDRmbc11NamZQPJJk+3TO8zOwILS5aKqK6RsXJV4oclid5u3XotUzRedtKy5WJ4JqCFfRZ5SbGyuxNd4zLx2W3aLFla6XZ3rYzDTfYPLNEiA03NimMr8qZTtw127d0r49/va9k77NVFAOB2C3j04RJ897pCqCrw26fb8fyLbVl7E0VERPln334Z23dok4KGBgWbtyRu3jy/Mnj0d1IVs/lsmleIU2OsfigIwNz5QaxeK2Hv3gSzT5vGFRl01N8HWC03EWPMJUMveahBwY6I12zvMN6Q9RuSf/of2fNna41s2J8n0vIVWinajp1mZ/7JjjA3WDn2V6+RsGKllNLnZvfuDB1bZl7XrqmRXTOsFSvFpB4IJBIZ2F+6TERdomuNRXY0TTdSUyNh7rxgzJLLvXuNr1ehoEboYYn+s448xmMNN/LnYl3rAIT1Pkv0OtE/YO3n7Xroob1W9Jut35D4wUukBQuCpgP0TpQbtrerqFooYu68/Otd3W2DXa/9qQOSBJx2qhdTJjOrizQul4AHfl6Mu+8sAgC88WYnHn6kzXLzQiIiomSsXCVh02YJzX4Fi5eK2LZdjgqGRGppVbFpc+IbZDP350ZPwjPyzCfOZCsk0dNuuyY5Vn43amJocV8tWSpi4yYJhyIyyWL1mLFbvAkqADQ0auPabaIUzSnOlTTG/379AW3f7akzFzhJ13P3AwcU0yVkqezLPXtkbEq2HMvE++6vV7Bkqf1l2zt3yWFl1LGCzfqPZ+9eGQuqgmh3IBNmx86ujKwt1TLa2lVs3yFDFNWoxu+r1kiWy65TzZpsa1PR1Bz/mFdMnBKh8yGVbFFXGiIdDY0K5ASXw8hzWcreyycAa2Xt6Qg4p5PH6QE4YXetjP98pGV13Xl7kcOjoWwjCAJuubkYvXu58OgTbfjgwwAaDyt49ukyFBUxA5CIiNIvqGsdaWXlKr1mvxJWAmKqjDHW1yNeR1VVrN8go3dvAQMHuJMaXyxmJkJG7Jp2Rv6lT0cwJXIyFDlpXrVaQnGxAK83vfcddgcPMhV4cuoRZFC3u6wGOrdslVBW1nW+NDYqaE3y3O5608Q/EioR7tHD2sw/ZpZnDGuOlOr2OcaFY/qkFmXw+5WMrtBopZk7oAWRAGDDRnO/mMxKgkFRhdcDw+qjffuVsICSqgBzvgpCkoExx4Vfj2UJQIK8jljHrumMTnQdinPna9lBF5znQ0GB8fXLzMOZWGOzco1xuQSk42qReo85e8Zhl3wucOuWmV2vvt4BWQbOPMOLU8Yzq4uMXTmjEC8+X4aCAmDefBE3/8hvW3M/IiKiePRPT5NptbBzp4wFVSJWrkq9JEJV1aib8z11CnbXyjEb/6b0fgm+f/CgglWrxaiMBf0YJRmorAqiOsEqhobNs7Pkxr/dIBASK8tBFGOvzBnv8PG3ZNmsK08dPKigZlv4+bJoiZjRkqVAp7XPeuMmKXF5tMFLihEldUET8dTI47qySsTS5fZncdnNbOWHvizSTKDD79d6Ay5bHvv4UCOmJKHsoUOHrH3Obe1q2AOG2tquz7wuRhkkkDgIFS9bKF7QOOp94n877rj0mV129s5LJli1LAeO53zU7YJd23fI+OS/2qOC239c7PBoKNudd44Pb7xejvJyAavXSLj2e83YtSvLc1GJiCjn6SeZVmNd6zdIhqs7mbk/N3sTH1rFywpZVnGoQTEM2JhZDSw0tqXLRezdp2BLRK8dfcnN0mUi/H4VW5NabSt8h1uZ2MT72XiliWbfw6h8RlVVfDE7iFmzgxlrEN/WZtzg3ujd45VRJZ3hkGMxulilcQlleDs7IsZpJQMnntbW+A+L5y8QE5a+xZPyKqCCFjBevcZaQMLMSrXJnJO7dmvbE2uhDADYtl13MYjzN0KStAU+YjVrnzsviKqFXb2akrlmhui3NdkM3UiJgmpmryFG40n2+pPM7+kfLGRbZlc+63bBrpf/2A5FAc49x4uTxnXLKk6yaOIEL/76djkGDnRh924F136vGcvTsNQ7ERGRHQyzlRKId/Ndu0dBfb21zAQja9ZKWLJUTLjSoxKjWX2kzk41LMAVmemQrMjgon4IVrZ9zx4Zm7d0BQsig3PJMHp/fQlWME39hfWrdPn9CubOD+IrE82Mq6slzJodDGuUHQioWLNWDOuTlIqabbLpcjI7xTsUIrdNf2y2tauG95GBILC/PjKAa99qk2aC5vMqrQd7zHzNTCn2egcac+tt3yHHbUxvZv81+xX4/eGvEa9fk9+vRAUC6+tl7K5NfK0I9YRLZPUaCavXSlizNvb+TbRSnxGjoJP+2i2nGHhv9iuYNTtgucTU6UBSTY2UVMamLKu2neuWZUk2czp0q2DX2rUiPv8iCEEA7rydWV1k3qiRHvzf3ypw0jgPmptV3HSLHx9/ksFmAkRE1G2ZzeyqrpbQ0hJnAmQqtcv4y9sjmuQn009n335tbIka7usnelZu/ePNE4Kiin37ZVOlLKn0L9GPYe8+JWZTfbMrnkUyzOwy+bt2OXBkkm3U3D5yu0JZIut1wagNGyXsqVOwcHHswIqVptp79ynYuUuOf+xngH7EclSJWdd3V60SYwYqVqyUwr731dwg5lsMQO3bL6PZb/D6Dk5ozXya6WgmboUdgeIFVSIqq8SwDCcp4uMLnSOiqKKySsS8SjEsyLF8pb1Bv1BGkdngGKBlaK1YqWXHWqGPbykmY/ux3kEUjUtgk87GSu7XLL/WlurEG270+zO/iF+2mkl798lRQdt4srnnV7cJdqmqimeeawcAfPMbBRhzHLO6yJo+vV14641yXHiBD5IEPPDLVvzxtXbnovBEREQ6W2tkzF+QWuaxmT9pjY1acCFdqnVlNPHGEy8DS6+zU8Ws2UGsXCUlzCozfF2b/syn0nsmJNdWwgrRb4uZFexmze6KPHR0RGc8GO2bWCukqaqKVavFxP2n0kif8RJZKhjp8OGu78sKLDexX7lKwoKq6OtAOiakZs8NMz/ndjs3YxYEwJNgnQ2j/Rdru/TZXJEB9kWLRRw8pIQF11K9xujHFlpBNRX79inYXx//dYwyu4yya+2cJ23ZKiEYNP96Zsrj48n0FC9e2WoyklncpqFBwarVEioNriG5qNsEu76aqzVpLSwE7uAKjJSkoiIBzz1Tiht/UAgAeOXVDjz4y1ZLF14iIiIr7JqkqqrWg2ft2iPL1m9N7ilyrMwsVVVRUyOhMYnJVsxN1P15NZo0Ra4SaeSruV2zyp275ITlc/F2d8JmyvF+IGpyaP3eYe06KWryHDTon6b/DLwem4MIBi8XFFUoSvRCBkbcFhbvPNSgYM7cIJYuS3ysxnrvAwcU7N2nmOo/1dKimP9c4vzYvn2xjzEzDdvzjZksn6CoJtUL0K4+de4k8yAMD5c416X2DhVLl4loaU29NDxd4pVexqKqaljAObRNdm5bzTYZjYeN/g5oAfH99XJYgCfR3wezCwzYLVOft9mm/Pp7jGQWLdm5y1zWtBO6RbBLFFU890IbAOD73y1C/372LpFN3YvLJeCn95bg0V+VwO0GPv5vED+8lSs1EhFR9tu+Q0ZtndYnpmZb9IzGzE2412scPNlTp2BLtYxFS+ybzVuZFMQqQ4qcCy9cLCIYjNMfJTKzy2QPsUTiNla28LqRjZZXrekK4oQmLfrPoKDA/GsnIxjUMufmzjP3uVtZXTS0KFBkRpvZXlEATPf82blTy4xMeoVR3fub6bmUCfpjPC2ZXSZ/zkzwxO9XMXuO9VpCOwIHNTWyrWWUYUGWGD+zQleumKngx86dsqmAYrL7wui1Ez8gSD37q3KBiE9nBrFipYSqRcbXoci3qKuT0WyxTNNJdXUy1m+Q0l5RlMzLiyKww8HM2Xi6RbDrH+91YsdOBT17CkczcohSdfWVhXjtlTKUlgpYvkLCdd9rdjRFnoiI8pNdc1T9TWwqK2V5YmRAxFtqPpFY2S6JbuzVJAMMlXHKPa1O9Natl7DmSLaclYlCsp/BpoiVNvV9dYze32oZnFWNRzLlOjrjZ7WFVoPTZ3aFjdfgQLcSoIm17yP7ZwHAwYPRXwytbhfqLWfl/erq5Ljb70SBnpbtkt4xGO3zhkYlamVFVTEf0LA6mbdyHsVa8bGhUUl4rFkpY7R6bquq+SwcI2bPkw2bJCxdnjgo7TLxepGjPXBACVsowmxm18KFIubOE5PO0AsGgZYYf3viZXbFatififY0yQSGVq+VsGu3bPr6lPD11ohHg5Nh1Z667beSyBHKCDt8WMGKldmTvpr3wa6GBgWvvNoBALjrjmKUleX9JlMGnT7Nh78dWalx124F1363GfMr07QUEhERUQrM3MPbdZ+fqC+RWWrM/9EkGzDqDJgruQMSNL4PqthdK2NPnZJwmyNbHpjJeDBSW2dto40ayackYlxu3cw4Xg+jL2YH0damhk3MEz0ktBLs0pdQ6csIjTI6jCb8qWT2rI6z0p1TFiwUHSmPq92joGph+P61sjJfMoEicz+nRo0rXcyuKBvS1Kxi5hfBqEB2OphpOp9MFmAgiLBMqdC/El7noJV2pvKwJBb9NTlyHLG28eDB9J80e+OUOSci2jTNrNurGK4aWVvbNTYrSRw+n7ZDFy4WE/Z7y6S8j/y8+HI7WlpUnHiCG1dekeY8buqWRo/24B9/q8DECR60tqq47c4WvPFmBxvXExGRPWxKyTCT4WPmL1fMP2+6r69cZc+k8kCCFcR27U5/RnW8eXpTc0QAy8LrmLpNsHgr4cSthz5IFKtJfEjtHmufV6wJqdF2Vtd0TdzaO1TU12vvFdk0O9ZEzOoE38q+3m9hJTwr94/xftLvV8MWA3BFpOu0mVgoIJFtBqXQRsyuzAdYD3ZF/rwkaZ99ZKZUot1qZvEKIHqBBaPXnTs/2LUipondvPlIP7nIVW+dItgQITia2ZX6SyVNf4ybHcf+/dnxGcTSGVCxY6d8NFM2FaFjWX/tS/a6kK0rMuZ1sGv9Bgnvf6Ctjf3g/SWOrvJB+a13bxf+35/Lcc1VBVBV4PkX2/HzB1tte7JNRETdl103kUuXdQWgUnnNWJNG/ZebmlW0tiqor5cxd14QLS3JPemt3RP+e5HZUY2Nyf+dTSazK96vpBIsMZrgxw2cKarh0vBbq9OXGRKrBExflpio4XOiFe8imSmnCokMPC5fKaGpSYnaj7GCXZGZXZ2dqm3Nz81k01i1fYecMDAUb3XWufNSTxFpN3mfKyvmzzerwa6WiIban88KYvlKCes3dJ0Lfr+C3bvtyTZZsrRrn8bbpA0bzJ+Ldq7GaAdTr5dgzIcOKZAkFQdMZvmkku1kisl9rM+ctfK5BAJqwsVP7FCzTcbGTeHHd7JaWlVUV0sxH+hY2f5sDXYlue5E9pNlFU882QpVBf7n6z5MnOB1ekiU53xeAY88XIrjx3jwm9+14dPPgti5sxkvvVCGAf25KAIRESWv2a9AAFBWZs8dZazSNjszg+ZVdr1JvEm3FbO/7JqgC4K9DeP318vGixhZeI9ke3YZ9ZXSXs/4BVetlrC/XsHYE8Jv5ffUKRg4MD0TroULRVx6SXiVhLZ6Ydf/J2oG7/YIiLVD24wyDy1kdhmp3aOguNjcz+qzWZqatVK8HhWxzzenE/jNrC4ZJsam7K6VUVsrY/Kk9M2VrASwrPysKKoxe1DtqVPQt6+MAf3dqKxK/foTWljBbIAvJJPl43ax46/Mvv0KDh4Kml4gYtt2GX36WM/DsSsgnaov5wQzmsVWbxBEDARU7NwlY+gQN4qKzH2KW2tkAKlnswnQVlTNNnmb2fW3dzqxfoOMsjIBP7mnxOnhUDfyrWsK8cbr5ejZU8DGTTK+dW2zbeUcRETU/cgysKBKRGWViAU29Zyxc8XEkHRP2FREh0nsfM+164xnZWGTqTjvZ+XJdkODEnZvoFhsTh3KTjIqe6qpSU9219ER6rZz/oLwvlCJPo94u6jxcPQvW1m5MZbIMsZY9FlkoV42kdliYa+bffO6pKxbL6GpWQ3LWLJbovJWvYOHzEe7OuMtigBg5SrJtiqLhkbFsAzXyiIaqfxMJqVUzq5jNtAVEpmlZ4bZfafvP3XggBK3PN2q+no54+WaoaDwps0S1q7Tzt15lUHUbJMxZ27m+0dv3ylj1uzs61udl8GuPXtk/OHldgDAT+8tRt++ebmZlMWmTPbi3b9XYMwYNxoaVdxwsx/v/auTfbyIiMgy/XLu6SiHAnSTmyQmZk49zd1fr6ChMfkspqiGxTF+ztKkyGxpJMIDKbEzuyy89xFGQSNbRbx8p+74TJSVY2lXqmrMpvFmXycYNL8Qgf69zJxnVjN8HJdguLFWtLPlrRXzn4NR0+xUJAqIWWEUEE9lZduQbJoeDB3izqrxJGJ2rKHVeoNBFctW2BfYnTNXK5u1KtU5Yei3t++QUbtHQVubGpaxvXBREHV7Uz84GxoUNKbwd9ZpeRcFUlUVjz7Rho5OYMpkD66awab05IxBg9z429sV+NpFPkgS8NgTbXjoV23s40VERFnLyqppIQsWiGhqVnJqggQgavWvmM3QTd7nq2ryk1arGRAh2dAmpbOz69/JrjBpZH6lmHIfmP31ypEyHc3Bg7E/TDuyyLKZZJA9aHXBgGTFCuZmwsLFaayuMHHOZyKzy84jN9Xy8FRsTGI1SitjVVXV9tVpk53XiTbEdPUBs8hyzsNNqqXVFI2IIrB4qZiWTPBMybtg1z//FcCixSIKCoDHflWa93+4KLsVFwt49ulS3HNXMVwu4D8fBfCd65szsnoUERGRFYqiYoWJJ9SRN9UdnSqqFoqW/rYZNVdPxO4J2Jq14dsa655R/77+JMpskpVK8Cxd2tpUdAbCv2alBFOfpZhIa5saVv7W3q5i8RLx6CqLyYjV3wnI3gbLdtlisNpgrNJdu9mR/RRipcwxE2Kdo4ebVNPHe7ad58lk+DrFyjV5a7XctUqmCencRjtK/qyUkCdDTLDgSC7Iq2BXdbWEp55pAwDcdUcxhg1jU3ByniAI+OFNRXjj9XL07iVgyxYZ11zbjC/nZF9dMxERZZ9MTSraTS45bkeZnJmgWrpFTsBjBTvMZqU0Hs58dlums8VXrooOaoa1NEuwA7bvkC2tWKYvLVy3XkJDo6KVDNm82S2tKg41ZEcQJV1lqG0mz+90sNKzK5Gly0TIRwKs2RBwiTeGVaslk/2vsmBDjsjGILtdarbJWLXa+b89dpk3P70ZV2azmrNZ3gS7OjpU3Hd/KwIB4IzpXnz/u4VOD4kozKlTvfjnPypwyiketLaquOveFjz7fFvCZbqJiKh7y8eJR3uHiro6Z7Oco/ZrimWMTkyiMn1oGGVR6IOGZo7VoIX5mb5pvL730tZqe/f13n15MKvLYopi73Vsz57s+LxUNX7WmtmegvGqx80sOGFnVuLuWq0HFGW/dPcNTGcfv0zJm2DXU8+0oWabjN69BTz5RClcrjzPRaac1K+fG2+/UX40GPvm2524+Uf+rEvJJiKi7FF/IDN/IzLd+mF1RBnhosXBuOWN6Q76dXSoWLgoGNXTZdt280E5O8eoqpkPZiVDHwy0s1wNAFy6Io3Wtq69weBUboksfU6Vkz3A9MxkQZm5JsTrI7WlOvH1x2xWrllmgnT5+BAmUnOaFoSxQpbNPRjS906kLnkR7Pr7Ox34178DEATgqSdL0ad3XmwW5SmvV8D9PyvBc8+UorgYWLZcwtXfbsKyOH0kiIiI8k2BL/z/Gw+rqKzK3N/CToN+Ooeb1KgyvQ4bV3PLR2FljDa/Np9d5wdZtj84srVaClvR1Al2Z6wlq7YuS6J/ZLtFi8WoB0NG7FxhMp/kfFSockEQTz3TDgC49+5inD7Nl+A3iLLD1y4qwLvv9MDIEW4cPKjihpv9eOnldpY1EhFRt+D1WotkHE5TLyM72T3CbJhIJ9LWlr5oVw5sPpmgqsCOnfaVLVdXS6iukbFuvbP9l1SYOEfz9CBWVaCpmUG2dMuG7LJcltPBrk2bJfz0561QFGDGNwtw4w/Yp4tyy4hj3fi/v1fgissLoKrA63/uwPU3+h3vY0JERN1Pppskuy2uI8QMq+ykb8Vg9yGUC8E+SkxVgeoa++5tpWy5TVYTl2jm6zHc1q6twkuUzXI22LVxk4SbbvGjrU3FlMkePPxQScZ7TRDZoaRYwK8fL8UzT5WitFTA6jUSZlzTjE8/CyT+ZSIiIjs4sAKX1WBXt5RjE2W7A5J29wAjZ+j7reUTM331ZJv7lWWL5iaenJT9cjLYtW69hJt+6Edzs4qTT/LgpefL4LOYCk+UbS69pAD/frcCp4zXVmv82QOtePCh1rjNeomIiOygqvFXBEsHt5v3bnHl4BxZsrmqTBRzcCdQt6GqasLVWleucrbUMl3MNM4nclrOBbs++jiAG25uhr9FxSmnePDn18pQXp5zm0FkaPBgN/7yZjluvaUILpd2vF9xVTOqFgadHhoREeUxVUXCSZvdPJ7Mvl+uaWtX87YEyqzaPXzgR9lLUZl9SJTNciZK1Nqq4MGHtEyXjg5g2mle/OnVcpSW5swmEJni8Qi48/Zi/OXNcgwd6sL+egW3/LgFjz/ZijablxYmIiICtDKjqkWZ7b/i8XC1vXh27JQNV4wkouzgREYsEZmX9ZGi1lYFr/2pHRde0oSPPg7A5QLuuK0Ir/+xDCXFvEOi/DVxghf/frcHrvu2tvDCu+8FcOXVTVjOpWWJiCgPHDigcKKYQL6WQBHlA7vLdonIXlmTQK4o2ooOhxq0XNCDBxUsWyZh1WoRHZ3azxw73IVHf1WKyZO8Do6UKHOKiwX88sESnH+eDw890oraPQquv9GPK2cU4Kf3FqOCJbxERJSjAqzQJyIiojTJmmDX0mUSbr29xfB7I0a4cestRbj4Ih+bmVK3dNqpXvznXxV45rl2/OvfAfz7/QDmzg3i/p+X4NKLfVyJlIiIiIiIiOiIrAl2nXiCG1dcXnA0s6ukWMDECR5MneLFyJFuuNjUgbq50lIXHvtVKb5xWQEefaIN27fL+PkDrfjwIy8e/mUJhgzmGu5ERJQfCnz5k/l11hlezF/AFgRERESZlDXBrvJyF379eKnTwyDKepMmevHv9yrw5lsdeP3PHahaKOLyGU24/vtFuPnGIvayIyKinHfmGT7MnpMf0S4+sCUiIso8NvwhykE+r4BbbynGB//sgdNO9SIQAP705w58/RuH8cGHnVDY8ZeIiHKYz+f0COzDTgNERN2PzwuMHMHKGycx2EWUw4YPd+ON18vw0vNlGDLEhYMHVTz0qzZ8+zvNWLGSJRNERJSbBEHAqJH5MUlw8W6bqFvo05snO3UZPNiNslI+7XASz0iiHCcIAs4/z4eP3u+B+35SjNJSARs2yvj+DX7cersfGzZyXWQiIso97jy5S2Wwi6h74Lmeuny57gNaVm+QuQeOyqPDiah78/kE3HB9ET79qAeuuaoAbjdQuUDENdc24657/NiylUEvIiLKIfnyQDxftiMNfF6nR0BkHwa7Uher7Nub5LXizOnJX2TGHJdadrEgALLE1jJO4ilJlGd693bhkYdL8fF/euAbl/ngcgFffiVixtXN+OnPWrB5C4NeRESU/fKl11U+ZSoMGWzPxpSXCxg9yo1TpzoX7Ro6JDfLZE84Pr3ri5WX58mJ54BsvGadNM6D6dPMnWdlpQJOPMHZ9euEGJeYZHtflZe7kj7XUy1BFARAVlJ6CUeNGJ6b10i9PPrzS0R6w4a68dsny/Cff1fg4ou0Tr8zvwjiymua8cNb/Vi8RISq8mkDERFlJ1k2/np5WRbOKGOYPNGTlRNgvXFjzU9um/323De4BOC40R4UFTm3c0Yc68Y5Z+XeSggjjnVj6uT0BAldAqDEOO8oMU8WxgaGDnGjR4/EU/7hw9w4Y7o3ZiLq1y8psHdgMaTjijBoYHIhDyHFi7eA9ATVk90eqwZk6H3SKfe3gIjiGjnCg2efKcP7/6zAJRdrmV4LF4m46RY/rv52Mz79LABRZNCLiIgSKyrUnvy7MhCjCAaNvz55UvREf+AA529pB/R3YfxJXYGjC87zoV8/N1wJdpbH2UQKa6VXMW4XrG6DmA19bASguNjpQSTnmGPSc7wLgnGQOZcCzE5yuXNzP/XqKWDsiZ6Y16pMZnvpx6AvXTx0yPpcJRR87NUr+fNl8kQPJoz3oDiJwLzgQloC+iyXNY+7iqibGHOcB7//XRk++7gHvnNtIQoLgU2bZfzsgVZcdMlhvPzHduzbz8d5REQUW79+LhybodIGRTae3BgFVtKRPdW/n7XbZEEIH0dBQeJB9eopOFrKB8BS4DLWfj5zug9TDIKQsYQesjmd9ZZq5kamJSorKi1JMRPFBcgG511plqwoN+EU5yLDZjIg3RnO7DppXPwxmS1fTHQeDBuqXQutZIEm249PH8g5blTX+wWD1oNd005L/drar58bAwe6k3ooka7LS2RQsiBLElSTLTVNJwa7iLqZwYPd+MUDJZg9sydu/3ERevUUcOCgildf78BFlzThznv8qKwKQlGY7UVE1F2df64Px41yo6gw/KY6VP1u5018v77Gt6Ox/gx5vUJUgMZoPKkO0eu1XrpjtTvAtNN86FHhwtQpXvS1mK1jV9Ax1mfp82rNnfv07hpXrAyPoiKgb4zP0YiSxj42fXq7MGmiB8OGxt8/Zo6P0aMS7+OiQiFjfdlOiJNhM2yoO2E/r94JMlwEARh1JMAweJDr6Lk5fJjzk9hpp3oxcIAbo0a6cUwfV9T50iiX7dYAACAASURBVKe3C+NODN9+u3rMAUh4PAHmgl0DB7gw7VRvws8ikX59XVHXZ72ePQRT5YsAwtqaGF5Lj3zNzD4AgHPO8uHMM7oiML16mr8a69/fpXu7ZIJN5eXG25/M+ZpMNlU6guk+b/RYysoyG9IZM9od9nchJFPllVZk34iIKCN69nThtluL8eWsnvj970oxZbIHigLM+UrErbe14MKLm/Ds823YWs2G9kRETkt11bpzzjb/6HdAfxcKCwWMHu1BWVmMYJeFO8h4E84hg10xv28UOLroAm07zARWUn6qbzFwFW9ek2jKc0wfF45LceWvZMULPJWXu1Bc3DX6igrjLbE6qYvVj80Op071on8/d8IsFDNDPm504hn2yJFunHZq7GPtdBPH4VlnpJ6BMnKEO+q8KIzILuyZIOggQAuinnWGFyef5MGkiR5ceL4PPXvaN2W8+CItkG5VqBRtzHEeTJ3ixYkRga0ePYSo5vonn5TZrEl9ECJWgNDt1rZlwikejDlOW6ThgvN8loPzxw53xzwfAUCycPuuv9YmWvUwdA2Op7g4/HWs9K3S70P9OXr8GPuy+goLBS0gl2ClRp9uU91JlKjaHevyeLQyfl/ERxDrgVE69Owh4NhjjRcXMZPNnGkMdhF1cz6vgEsuLsDb/68CH31Qge9+pxDlZQL21yt48+1OXHFVM664uglvvt2BvftY5khE5ITBg91hT02tZgWUFAsYa7Lviv4GPfJmPZQBMKC/ufefdqoXxx/vwYD+LvQwmJgJgoB+/dw45WQPhgwKf02jIIzXGyPYYvC1khIB557twzF9UrvdNduvKN7ExkzczGofNLvWmIk1KTZ6ef1E1OgYOC1OSeb007u+F3rtXFgnJ9FE0uVCzGhmr55CVKDIKLCVjsyMAf1dOP88X1j2z7HHulFeJsTsIxTKFiwrc0EQBAiCAJ8v8YFpZbLtdgthAVQjE8Z7wl7T6Bz0RlzO3C5rQfh08OgCIkOHGA8mVIJWUCBg1EgP+vR2oaBAsBQYGTrEjd69XfD5hJjBVCt98cKDXV0DGTnCjTHHucOC2V6vkHCVQkEQwq4VkSsSxnsAov8M9RlY+mNm0sTYf8tCx3uiLMeSEiFhea4+My6ZEtXIa/rgQS7TQc3IfTRhvAcXXaAFniMDqYIATJwQf3uPG+VOOajes4eA06f5DAN/U6d44fMJmDLJ63hpvh6DXUR01MgRHjz48xLM/bInXni2FBec54PHA2zdKuPZ59tx4cVNuObaJrz+53Zs286MLyKiTHG5wp8sFySxMNbw4cZPYyPFLTE7MilKNJEI6dXLBZ9XwMQJXkw/PTojIDSHGjTIjVILE36zQZLiYgFTpyR34933SM+uM6ZrN++J+iWl+kw7USP7SJGLy4QCAlbLG30WjiX9pNfoM+jVK/Y29KiI/nwjJ5ChbJDhw9xRn5vd/ZoSHUOhY/yU8R5MjjO5FhAe6Dz5pNg/e9K49K1AGRrDhPEe9O7lwtgTo8fh8wo48wwfzjnL+JwYOTK57MLJk7wYYyEzMVFgp29fV9hCFIa9gCJeY9Agt+Hrjo/zeVh1boIMWX3GTczAfIxtT7RP9Fl6RYVdX+/Z04Xzz/Vh3FgtEy/USL1Pn/AXPPvM2NdB/bnQq5eAAp8WqD1+jAejRkbvv9OnecOC10b01wr9Cp9jT/BgZMRrhj280Y1FiJHl5fXE3lmnjNcy/44dnvjvif6aGxmwjSybjLxWmfq7cuTlz5juxbixHow/2YtLLjY+hiIfypx8khdfv6QAl17swzln+zBwYFfQ0eMRwq4zBYXAgP7xzz+vVzAVVI/X9y9edl1o/H37ugxLHJ2SPSMhoqxRUCDgwgsK8OLzZZg3pyceebgEk44sn75ho4yXXu7AN65oxmWXH8azz7dh8RIxqcaRREQUW88eXTedLsF6z5BQ6aP+1rVPb1dYqc9xo9zo30/rIRMSCMS+noe+4/VqmWKJnvBHigwE6LcpMtYTL+gWGagQbF4e8rQjpXCANmnr09uF0aMTTOZTHYLF34/c5FOnejFpogcnHO8Oe/Lfr68rbklf32NchoGK0D5WY0SFkum7NepIMCVUxuZyCWHvPf10H8af5MHxY9yoiChJs7svVmirYvVKC2VWeDxa9mFISURWkiAgLHtqyODofTn2BC1TadBAl6kMngnjEwdopp3qxTF9XCgvE9Cvr1Z6DAADB2pllaGSIqNeR7ECq6n0GBo6xB1Vbp1shoeZa51+pNNO9aKoKLyXX+gzGWzweSQrXkaa26VN9Hv1FDBkUNfnESnWpSpWVtq4sVr219m6AGXkqo+FhQKGDXWjfz83pp3mxYkneKL668UKvkVth1vAeef64pbmejyCYfD6wvN9uPii6GBOZP/FkmIhrBxSX/Io6p6l6/e32UPT49ECL2aP5XPP9uH007yYPEkrJ40lMtgV2Ycs8noFdAU/K8pdR/udxRpX7CCoEHXNidTzSAaalSzm0hIBp031RmV8xzr3RhzrTmlVS6c4vNgxEWW7HhUuXHNVIa65qhCHGhR8NTeIL+cEsXiJiB07Fex4uxNvvt2JoiLtKccZ032YPs2LoUPN/6EhIqJogUDXvwXBet+uqVO92LhRwpjjwm/3hg11Y916Cb17uTDaoCdRvPIXfdxj+HA3vF5g9VptdjLmODe2bO16hG9UejRksBuqCqxbr/2O/sZ60CAXdtcKR2/YrQS7Ro9y4+BBBf36ubBzV3Il95MmerBipTau3gZPphNNwLWyr6TeWnt9i787dJgbtXVdO8nnE44G6PS91kIZMus3dM0ix4x2Y0u1tp9cLmDUSE/YZ2eG0eeT6O/+mOM8GDzIjRJd9kCvni4A2nsXFwElR4ITcsQMOdaCBQU+7ZgOBLXm36Ydeb2KCgEHDkZ/O9bn3bOngCmTvZg7PwjgyLnpEzB9mvfo75SVCmhpVTFwgLYtw4e7MfxIxp3RaoeRinX7Z8J4D/bUKVHN8nv1cmGqicnnKeM9WL1GigrW+rxA0EKpW0hhgYBOg4C4zyfggvN9+HSmtl/69XUZBgAAJAzsRpdPx/8Ztyf6a1MmJzfNHTPajbIyActXmq9gKC4ScNaZXrjdAqadFj/7qzhG5ozbBRi947Ch7qjG8IVxsjELCwXD7M7IoGdJsYC2dvXI98LHZDXLNCRWyassRX+A+uCbqgLl5QL8fhU9KgQcf5wbsqIFicae4IHHE34+qtAWUZEVYO68YNjrmukbpT+eiou7ymrj/a7+e5HZnAAw9kQP6g8oGDDABb9fRVOTGnNF34pyAc3+8H1i9W+HrDtYQsHCsWM9WLlKhM8r4PBhJax8NPITGHKkFLZ3b2DDpq4XKyoS4G+J/rzijc/mZ022YrCLiEzr09uFq68sxNVXFqK1VcH8ShGVVSKqFgbR0KBi3nwR8+Zrd059jxEwebIXU478N3wYg19ERFaMGOE+GqAQXMCwIW4calDQv78bTU2J02oqyl2GE6+hQ7TGxrHKFfTlcZGX7cgAh37SMHKEGzt2yJBlLRshVhP5IYNdWLde+7d+AuP1CjjrzK7xWll9q7BQy0YIBNSkg12RE75IRhNAr7crOGj2T1xRoYCOzujJRLyeMAU+IBDUyvtC2+dxa9l/h5uiX2vIYBdk2RNWVjhooAt1e7UPsKKiK8AUGvf4kzzYtFnqCoAkiMkMGODCoQYFnohxf+1CHz6fFTT+JSAs0KV/f+3f4dmMevpjbfrpXhw8oMDfomLMcW6Ulhofa8OGurFrd/zjIVbfnnifp34bQhk5+v4+p53mhd+vordBWWes150+zYtly0UMG+ZGjwoXRo5wo7NTxcCBbgwcmHx2Unm5K+y8Cpk4wYvFS61Hu04+yYOly41/T//5xWuGr98HF57vg88noK5OPho4N3O/aLhqoO6gSTZgE1qNUh960vdPmjzRg02bZZwy3oMly0RIEtCvnytuA/MB/V0YMsSNhgYlZi+vSRO9WLFSxIkneFCzTUZLqxq12uL4kzxoPKwFVKxyuwWcMd2LBVXaZ1dY2BXsSmalQSsSLfCuKMBJYz3Yt1/BiGPdYYGlUJBYv0q8qsIwc27YUHfamqSPHOHWAkgyMGGCN+r4KyoSjpb5VZQDQwbHfq2iovjBLjM9MYuKuv4dOtZLigWcOV071w8cULBsRfR5OmWSF/UHFAwbavwe48Z6UH+g6/odCkIaPUiYNNGDDRu0cyFbZe/IiCirlZa6cOklBbj0kgIoiootW2UsOBL4Wr1GwoGDKj79LIhPP9MumL17a/Xl2n9ejD3RnfGlcomIsl1JsQBRVFFRoZU9hIJdLpcAr7cra2B3LbB3XxJ1ZEdUxFiSHYif2RWZYaC/XRcEAecfKQOJN9GMF9DQO+F4D9o7RAwepGVtDdFNElNpbH72mVogYvUaSVeWaa4R/dTJXkiyipWrtM9lki5gUFgY7ze7FBUBHZ3RX/f5tPInFcAeXcbWOWf5UFIioLNTRUEBjga74k2uBcE4uyPWzwJaudfgwW789zMtpTC0j3v1dKF2T9d4Joz3oKVV7WqSHZFxmChoGKlHD61MNFHT8rDfqXAZllFFOvEEN/r3d2GJQVAnVGI0cIAboqhleFUt7Po5s8FLo2PY5xXQp7e5Fwj1UurRw4ULL+hK2bFz9TkjZbGyrmI492wfWltVHBNR9hmrxCpyJUg9fSN3Mw3wE2V2Hf1azP/p+tKUKV40NythmYwejzYmo75wQwa7wlZ17NfPfbSs9czpPtQfUOI2XAe0YNIxfVxxy8x69nThgvO1z79HDxd27JSPBnpCQudosvTX/dJSAQ2N2r+Tab4eafCgONvWoyu4bkRVtW3WB4wjhX3eBsfD1CneuPv31Klew+uAWfq/v0YiV0iMx+jY1ZdyTpyQOI27Xz83xp4AFBYZf1+JOGlCwee+fV1RD6JOP82LbdtlHD/GHRVEnD7Ni2DQOLjYv5/7aDZxtmKwi4hS5nIJOOF4D0443oMf3lSEzk4Va9ZKWLZcxPIVItasldDQoOKruSK+misC6IAgaKsCnTROC4CNG+vB6FHpeyJDRJTtjh2u9Vsxug5GZjnps6OMRAalrNA/hddnsJx/ri9mH5oQq9kUkb1n9Ip1T6kjAzf6+/hE/UwA7an8tu3aZKu01IXSUqB/fxfa21VsrZYxaqRbK8M6zxc3y6Froq8FuwoKtEyP+gMKjh3uxv79iQOQWnDJOFp38snaJGdPXVcNa+gzCO37cWM9kCQVRUUChg5x43CTZHl1TisGDXJBEDxHJ6L6LKNE+97n1SaJPXrE/jlBEAx7O0VOwPv11fpTRa5uGM//b+/eo6Is8ziAf9+5ocOMhCIqCHJzUERA5KKCaYaXBTXvWXmyU66lmUdbM7fa02l32zbbykBLsa3cLQtRPFpZkp1NE8VbGaUieeN+Vy7DIDDzvvvHOOMMM4MzMDPvMPw+53CAd17e+c08zPs+z+99LgKBNulkOETV30+AkWEio2FUIwKFJnOT3at3UUiQEI1NnMUejJZ0PqylXmmOJhFr/99/KVCjto61OORKx3C4l6H4ThN1R40V4dZNFn5+AotJaR8fBsOGCqxe6dRahp9dwyMPGypAZRWL4GChPul09aoGag0wdIgA42NN//8mJIhRXq7BqC4W45BKrU8q20IqZcwuMGAPSZPEqKjQDovV9Xr06eGKteNjRWaTHtOmStDcfO/PiDU3L3o6OsRw0nRxF7kkXQ9Ys4siWDDEV2DTtc/w/1QiBuLjtHPOqVTqeyZODXVOhhrqPBdhVzcHvL0FiBtvefVQa2/kuCJKdhFC7K5fP23FVVd5bWvjcOGiGgW/qvHbb9rv5RUsrl3T4No1DQ4c1FbqhUJtg2ZUuBCjRokQrhBiVLioV06ISAghthrub5rwHxUuRH09B38/4/OgpYp/+Egh/P2F3VrxzW+YABWVrNFzhYUKwbLaBqG5RFdPelgBppP8WsuwsdDVRMo6/n4CfbLr7jG0S8/Hjrt7MGtvuIyPFaHt9t3Ema6nhzXvx5gIIU7k29aYMmSYyBw+XIgBAxiLQ/EssSZO3T4Mw8Dfv3uxCoUMpt5jBTtLDP/HxWLtsSYnd+9Y2oa4NtnVrx9jMpRS93zS/gxUrV2/ObqwRo/uXjOKYRhMTBTj5KkO/UT9fPHwYBA7ToTaOtamya11+vcznTw7YLhQP0k/w2iTKyoVh6oqVp8YYhjGpPeKuf/JEYHaoX/mhu0ZngJ1f2uU7DJ4PCZahKARnNHQyuQkCUrLNAgJNl8G2vmM7FP/7Ol50p4Me0VOnSJBwy1tYrInLJ19+vdnrLoWWVoEw/L+Nu2uFxerHSIaHWX5sxsdJUJIMIcBXfR+1hkbKUJHB4fQENvOBYbn6+Bgof5GgjXXMmsJBAxCgoS4dqN7w/rdBSW7CCEOp61MiY0qNvX1LH79TY1f7yS/Ll5So6GBw5WrGly5qsFXh+6OF/cdzGDkSBFCQ4QIDRVqv4cIrboQEUKIMxUUFCAjIwPnz59HR0cHwsLCsHz5csyZM6dbxwsNESE05N77JU0Uw9OTsXrFLXPGRoowbBhrdAdcJGJMVvWyh2lTJVCpuC6HrXRlTIQIp053IDTEeNiF4R17w95wcrkAiQniLid2tkVPhm7I5QLMSJHYbR5La6+F1j7dYB8BautYBPagd6COLUN7zJmYKEbhZTUiHNTLpbOpU8Sov8mZ7fkRMVqE0lINFGYWdbDVwIECpM6y3/9AT4hEDIYN7WZZWxG+NrkC/WT9tuhqFVGz752FeAQCxmj+OkDbY9LRQ0V1nyXD1VFdiaeUgae057FZM3RZ5slA2WK+l5etyStzybHOcweaYzgE1RKGYYxWLO5KYED33ruQYKHNC4J0R1iYEK23uW7N8QZYv4KnK6NkFyGEF4MGCTB1ikR/x5fjONTUsCi8rMHlIjUKL2tQWKhGSSmLmloONbUdRnNpANokWGioSJ/8Cg0VIiRYiPvuY1yiAkkI6VtOnTqFp556CmKxGGlpaZDL5cjNzcWGDRtQXl6OZ555xiHP278f0+2kkSGRiHHa/BvW3vG3xNNTOyF9ZwIBo18+vvOwEh879dLoiqVLj26CeV2DjI9rVGiIEGXlrNkhe4Zix4lw8yYHH5/uxxg/XowrV9WIjupZT4WBAwWYNLGHGTMbMIzlubaCg4R2HbLmDvUUe74CW+cQM0dgprcXn+LjRGhvt77HaG8TrhCipYWzqgfc5GQx1Grzc7R1tfKuOYZFG64QQtXC2TS8mW+G1yZHDkEXi017UFojdpwIHR3WTRHg6ijZRQhxCQzD6O+4TLn/bsW2RcWhqEiNq9e0Pb6uXtUOfayqvpsEO5lvnASTyxkEBggQGCBEQKAQgcMFCAwUIjBACB8fSoQRQuxPrVbjlVdeAcMw+OyzzxAREQEAePbZZ7F06VJkZGRg1qxZCAoKsttz6lbmGxXOX68BV2hQdsZnw9JSgysxQYzLRRooRtpWVvZ8JTKZAH+YKYFAwOhX3JSbGf4oEjHw9e3ZM2snQXZeksoWXl5UB7AXe76XXgMESIgT2zQ/kPd9DG7fvru4hNHQRrtF1n0Mw8DDTr1JXVFYqPWpBIGAsdjT09aykhrcKLElBlfy4AMStLa6ZpKu2z09XVDv/O8ghPQZnlIG42LEGBdjfGdCqdTO+aVPgt35XlnJormZw4WLGly4aNpFuH8/3EmACTF0mADDhggwdOjdr8E+XS8fTQgh5uTn56OkpAQLFizQJ7oAQCaTYfXq1Vi/fj1ycnLw/PPP2+05x0SIEBbK78IePR2m5m6kUu0cVZ1XKJTLLU8AbM6EBDEuXFR3OYyrO3Q9CsRiBjOnS+yyCltvMWWyGA0NnFs15Phyf7IYZeUsQkPt+152XunxXiZO0H7QdDcxO0/8TVzbQG8GN29xCLByTsBJE8RobbVuPi1X168fc88FX0jPUbKLENIryWQCREUJ9CtX6dy+zaGsTIOSMhYlJRqUlGpQWsKipFSDikoWrbeBoiINiiyMlRcKtZWtoUMEGDiQwcCBAgwcKMAgb4OfB2qHDMnl2jt21FOMEHL69GkAQHJyssljSUlJRvvYE9/DY4YO0fai9e5itb2+xh5DPwYNEuD+yY7NJFozz4470S0mQHpOLhdg9Cj+Ew6d618MwyA5SQxW4x7zDbm7+Dgxmpu5LldrNeTtLYC3t4ODIm6Fkl2EELfSrx+DsDARwsJMH2vv4FBRzqK0TJsEq6piUVnFoqqK1Q6LrGGh0UD7uxXLxwPayZAHyBnI5Qzkcm0CbICcgadMu0KRpycDmScDqefdnz09tY9JDX6npBkhvduNGzcAACNGjDB5zMvLC97e3iguLnZyVI7HMAzGRlJ1khDiGrzcoNdPXyESMUYrZBJib1Q7IYT0GRIxg6AgIYIsTDCr0XCoq9MmvqprWNy6yeHmLRY3b3K4eZNF/U3tz/U3WTQ1cWBZoKMDqL/Jof4mB8DGGTYNCIXaCZd1iTCZjIFUapwgk3ZKlnVOqOn2kfanxBkhzqZUKgEAcrnc7OMymQxVVVXODIkQQgghpM+iZBchhNwhFDJWLUsMaFePVKmApmbtHGFNTRyamjk0N3NobmKhbOHQ0sKhRcWhRXnne4vBl4qDUsmhtVV7PI0G2mM09XxKVYFAO2+MzFPX44yBXGbw850eaJ23D5ALILvTM83cajmEEEIIIYQQ0htQsosQQrqBYRh4egKenkIMG9r942g0HFpbubvJMaOEmHYifpUKaGnhoFSyaFEBqhbOYjKN47RLOCuV2mRaVXX34pJItAmz/v11X0D/foa/393m4cFALNbOj2HyXQSIxIBYZPyYSAiA0U6WLBBoV1ASCLTLhjMC7XaGMfj9zuPMnX0NU3GGvdiMOrQxMNlu/Hf3/tncMaw/jukf2Pr8FmPpwXFs7fXHcdr/K40G+v8vltX9zIHlAO7Otrs/a7ezGu02lr2z/c7f3P3Z8JgcRowQ9tohKLI7kwE1NzebfVypVFrs9UUIIYQQQuzLqmRX652uB5cuXXJoMIQQQgAPifZroI2TcHIch9u3tZP0t7ZqE2C6RFqr6m5yTKUCVCoWLS0cVCptUk2luvNzi7ZnWasKaGxwwIsjLslSYozjtF/O4jWAwY4PBjhl4uxRo0ZBKpXa7XhBQUEAgOLiYkRGRho91tjYiFu3bmHcuHFm/5bqWYQQQghxJ/auZ3WHVcku3aSry5Ytc2QshBBCCOnjEhOd8zznzp1DbGys3Y4XHx+PHTt24Pjx40hLSzN6LC8vDwCQkJBg9m+pnkUIIYQQd2LvelZ3MBx373u2dXV1OHz4MIKCgtC/f39nxEUIIYQQ4jD2vuOoVqsxa9YsVFdXY8+ePRg9ejQA7fDFpUuX4vr16/jqq68QHBxs8rdUzyKEEEKIO3GFnl1WJbsIIYQQQkjX8vPzsWLFCojFYsyePRsymQy5ubkoKyvDunXrsGrVKr5DJIQQQgjpEyjZRQghhBBiJwUFBUhPT8f58+fR0dGBsLAwLF++HHPnzuU7NEIIIYSQPoOSXYQQQgghhBBCCCHEbfTO9b0JIYQQQgghhBBCCDGDkl2EEEIIIYQQQgghxG2InPEkBQUFyMjIMJm/Ys6cOc54epd34MABnDt3Dr/99huKiorQ0dGBN954AwsWLOA7NJdSXV2Nb775BseOHcO1a9dQV1cHLy8vxMbGYsWKFYiOjuY7RJfQ1NSE9PR0/PrrrygrK0NjYyO8vb0RHByMxx57DDNmzADDMHyH6XJ27tyJf/3rXwCArKwsxMTE8ByR65g2bRrKy8vNPvbwww/jr3/9q5Mjcm3fffcddu/ejYsXL6K1tRU+Pj6IiYnBCy+8gGHDhvEdHq9ycnLw5z//uct9JkyYgF27djkpIn5R/ci12Hquq62txZYtW3D06FE0NjbCz88Pc+bMwcqVKyGRSEyO0d7ejszMTHz55ZeoqKiAl5cXpkyZgnXr1mHw4MEOeU3uyta6szPK6ssvv8SuXbtw5coViMVixMTEYO3atRg7dqxdX7u7sKUMN23ahP3795s9TnBwML799luT7VSGjtWddplSqURGRgZyc3NRW1uLwYMHY8aMGXjuuecgk8lM9mdZFrt370ZWVhaKi4shlUqRmJiI9evXIygoyGxcP/74IzIzM3HhwgVwHIfIyEisXLkSkydPtvdb0OvZWoYZGRnYunWr2WNJJBL8+uuvJtv5LkOHJ7tOnTqFp556CmKxGGlpaZDL5cjNzcWGDRtQXl6OZ555xtEhuLz33nsP5eXl8Pb2hq+vr8WKVl/33//+Fzt37kRgYCAmTZqEQYMGobi4GEeOHMGRI0fw9ttvIzU1le8weXfr1i3s27cP0dHRePDBB3Hfffehvr4e//vf/7B27VosWbIEf/vb3/gO06VcvXoV6enpkEqlUKlUfIfjkuRyOZYvX26yPTIykodoXBPHcXj11VeRlZWFwMBApKamwtPTEzU1NThz5gzKy8v7fLJr9OjRWLNmjdnHDh8+jN9//x3JyclOjoofVD9yTdae62pra7FkyRJUVlYiJSUFQUFBOHfunD55mZmZCYHg7gAKlmWxatUqHD9+HNHR0Zg+fTpKSkqQk5ODEydOYM+ePZTwsoEtdWdnlNX27dvx7rvvws/PD0uXLoVKpcLXX3+NRx55BP/+97+RmJjosPeit+pO++fxxx/HgAEDjLZ5e3ub7Edl6Hi2tstUKhWWLVuGS5cuISkpCWlpaSgsLMQnn3yCU6dOYffu3ZBKpUbP8eqrr2LPnj0ICwvDsmXLUF9fj0OHDiEvLw9ffPEFwsLCjPY/ePAgXnjhBXh7e2P+/PlgGAbffPMNVqxYgbfeeosWiumku23r+fPnw9/f32ibUCg0+xy8lyHnQB0dHVxKWstXfgAAENZJREFUSgoXGRnJXbhwQb+9ubmZS0tL4yIiIrjr1687MoReIS8vjysrK+M4juN27NjBKRQKbt++fTxH5XoOHz7MnTlzxmT7mTNnuDFjxnAJCQlcW1sbD5G5FrVazXV0dJhsb25u5lJTUzmFQsEVFRXxEJlrUqvV3MKFC7lFixZxGzZs4BQKBffzzz/zHZZLeeCBB7gHHniA7zBc3q5duziFQsG99tprnFqtNnnc3OeSaLW1tXEJCQlcREQEV1tby3c4Dkf1I9dky7lu48aNnEKh4D777DP9NpZluRdffJFTKBTc3r17jfbfu3cvp1AouPXr13Msy5ps37hxo31eRB9hS93Z0WV1/fp1LiIigpsxYwbX1NSk315UVMRFR0dzKSkpdP43w5Yy1JVVaWmpVcemMnQ8W9tl7733HqdQKLjNmzcb7a/b/t577xltP3nyJKdQKLhHH33U6DgnTpzgwsPDuccee8xo/4aGBi4uLo5LTEzkKioq9Nurq6u5pKQkLi4ujmtoaOjRa3Y3tpZheno6p1AouPz8fKuO7wpl6NA5u/Lz81FSUoLZs2cjIiJCv10mk2H16tVQq9XIyclxZAi9wqRJk0yyo8TUjBkzEBcXZ7I9Li4OiYmJaGhowOXLl3mIzLUIhUKIRKadNmUymb7HRHFxsbPDclk7d+5EYWEh/vGPf1i8K0HIvdy+fRvbtm1DQEAAXnrpJbP/S+Y+l0TryJEjaGhowNSpU+Hj48N3OA5H9aPeTalU4tChQwgICMAjjzyi384wDJ5//nkIBAJkZ2cb/Y3u9z/96U9GUwksXLgQoaGhOHToEJRKpXNegBuwtu7sjLLKycmBWq3GqlWrIJfL9dtHjhyJhx56CCUlJcjPz+/2a3VXjmz/UBk6ni3tMo7jkJ2dDalUimeffdZo/6effhpeXl7Yu3cvOI7Tb9eV4bp164yGGk+cOBHJyck4c+YMrl+/rt/+7bffoqmpCcuWLTPqRe/r64vHH38cTU1NZoe79mWOblu7Qhk6NNl1+vRpADA7JCEpKcloH0J6QteIpMakZW1tbcjPzwfDMCZdRvuqoqIibN26FatWrcLIkSP5Dseltbe3Y//+/di+fTt2796NwsJCvkNyKXl5eWhoaEBKSgpYlkVubi4yMzPx+eefU3LZCnv37gUALF68mOdInIPqR67LmnPd+fPn0d7ejkmTJpnMgenr6wuFQoFffvkFbW1tALTX319++QXBwcFmG/dJSUlob2/HL7/84pgX1Yc5o6x0n1XdZ9eQbo6ZM2fO2O019WVHjx5FZmYmPvnkE5w8eRIajcZkHypD/nVul924cQM1NTWIjY01Garo4eGBuLg4VFdXG9WXTp06BalUitjYWJPj666dhmXS1XVVV4Z0XbVeV23rs2fPYufOnfjoo4/www8/oL293ewxXKEMHZoZuHHjBgBgxIgRJo95eXnB29ubGgGkxyoqKnDixAkMHjwYCoWC73BcRlNTE3bt2gWWZVFfX49jx46hsrISa9assTghYF+iVquxadMmhIaGYuXKlXyH4/Jqa2uxadMmo22TJ0/G5s2bMXDgQJ6ich2//fYbAG3Pyrlz5xrdqRIIBHjiiSfw4osv8hWeSysvL8fJkycxZMiQPjOBLNWPXJc15zpd2Vi6lo4YMQKFhYUoLS1FWFgYSkpKwLKsxf1122/cuGG2sU26zxlldePGDUilUrNzruk+47rPPOmZzotEBAUF4Z133sGYMWP026gM+WWuXWbN51C3X1BQEFQqFWpra6FQKMz2lDcsQ52urquGxyf3dq+2dXp6utHvgwcPxptvvml0/XKVMnRoskvXPdSwO6ghmUyGqqoqR4ZA3FxHRwc2btyI9vZ2bNiwgYahGWhqajJaMUMsFmPjxo148skneYzKdWzfvh2XL1/Gnj17IBaL+Q7HpS1YsAAJCQkICwuDRCLB1atXsXXrVhw7dgyrV6/G559/3udX+KyvrwcAfPzxx4iIiEB2djZCQ0Nx6dIl/OUvf8FHH32EgIAAPProozxH6npycnLAsiwWLFjQZ87hVD9yTdae65qbmwF0XX4A9PvpvptbbQwAPD09AYCGMTqAM8pKqVRavOmjOw6Vbc/Ex8dj2rRpiIqKgre3N8rKypCVlYVPP/0UTz75JA4ePIghQ4YAoDLkk6V22b3KxNbPobky6eq6KpVKIRQK9ccllnXVth49ejTefPNNxMfHw8fHB1VVVfj666+xY8cOrFq1Cnv27MGoUaMAuE4ZOnQYIyGOxLIsXnrpJZw5cwZLlizBvHnz+A7JpQwfPhyXL1/GxYsX8f3332Pt2rV499138dxzz0GtVvMdHq8KCwuxfft2PPnkk0Z3A4l5a9asQUJCAgYOHAiZTIbo6Gjs2LED48ePx88//4yjR4/yHSLvdPNMiMVibNu2DVFRUfD09ERcXBzS09MhEAjw8ccf8xyl62FZFjk5OWAYBgsXLuQ7HNLH0bmOENe0cOFCzJgxA0OHDoWHhwdCQ0Px0ksv4Y9//CMaGhrwySef8B1in0ftst7vXmWYkpKCefPmwd/fHx4eHhgxYgRWr16Nl19+GW1tbXj//fd5itwyhya7OmdpO1MqlRbvtBDSFY7j8Morr+DgwYOYO3cuXnvtNb5DcllCoRDDhw/HypUrsW7dOnz33XfYs2cP32Hx6sUXX0RAQACee+45vkPptQQCARYsWAAA+Omnn3iOhn+6611kZKT+7rLOyJEjERAQgJKSEjQ1NfERnsvKy8tDRUUFJkyYgICAAL7DcRqqH/Ue5s51urLpqvwM99N9t9QzpKWlBYDlO+Ck+5xRVjKZ7J7Hp7J1jEWLFgEAfv75Z/02KkPnu1e77F5lYuvn0FyZdHVdValU0Gg0dF3tQk/a1vPmzYNIJDJqD7hKGTo02aUbi2lubGVjYyNu3bpldkwmIV3RZZ337duH2bNn45///CcEAuqkaA3dhH99fYLGwsJCXLt2DWPHjkV4eLj+a//+/QCAhx9+GOHh4Thy5AjPkbo2b29vAEBrayvPkfAvJCQEgOWhMrrtt2/fdlpMvUFfm5heh+pHvUvnc9295vApLi6GQCDQJ3ADAwMhEAgs7q/bTvNp2p8zyspwjiFzx++8P7Efc/UQKkPnsqZdZs3n0HA/3fxpZWVlZhchsFSGhsfq6vjEWE/b1hKJBJ6enkZ1XFcpQ4dmCOLj4wEAx48fN3ksLy8PAJCQkODIEIibYVkWL7/8MnJycpCamorNmzf3mTle7KG6uhoA+vx7tmjRIrNfupPstGnTsGjRIoctie0uCgoKAIDeJwCJiYkAgGvXrpk81tHRgZKSEkilUprM38CtW7fw/fff47777sP06dP5DsepqH7Uu3Q+18XExEAikeDEiRP6Icw6NTU1KCoqQnR0NDw8PABoVxuLiorC9evXUV5ebnL8vLw8SCQSREdHO/iV9D3OKCvd51n32TX0448/Gu1D7Eu3oqJhPYTK0HmsbZcFBQXB19cXP/30E1QqldFjbW1tOHv2LHx9fY0SGQkJCVCpVGZHD+iunYZl0tV1VVeGdF01ZY+29Y0bN9DY2GjSHnCFMnRosmvixIkICAjAV199hUuXLum3K5VKvP/++xCJRJg/f74jQyBuxPDDOGvWLLz11lt9PmljzqVLl8x2/2xoaMC7774LALj//vudHZZLef31181+jRs3DgDw9NNP4/XXX8fo0aN5jpR/V65cMTv07uzZs/j4448hkUgwY8YMHiJzLYGBgUhOTkZxcTGys7ONHsvMzERTUxNSUlLMLuHcVx04cAAdHR2YM2cOJBIJ3+E4FdWPXI8t5zqZTIbU1FSUlpbi888/1+/LcRzeeecdsCxr0ltxyZIlAIC3337bKOmyb98+XL16FampqTRMygGcUVYLFiyASCTCBx98YFT/+v3333HgwAEEBgZiwoQJjnqJbq+2thYlJSUm26urq/H3v/8dADB79myjx6gMHc+WdhnDMFi8eDFUKhW2bdtm9NiOHTvQ2NiIxYsXGy12pCvDLVu2oL29Xb/95MmTOH78OOLj4xEcHKzf/oc//AFyuRyffvopKisr9dtramrwn//8BwMGDMCsWbPs8trdhS1lqFQqUVhYaLK9sbERL7/8MgAgLS3N6DFXKEOG63ybw87y8/OxYsUKiMVizJ49GzKZDLm5uSgrK8O6deuwatUqRz59r5CdnY1z584BAIqKinDhwgXExsbqs9spKSlISUnhM0SXkJGRga1bt0IqleLxxx8322hMSUnp8wmK119/HXv37kViYiL8/PzQv39/VFRU4IcffoBKpcLMmTOxZcsWGvppxqZNm7B//35kZWUhJiaG73BcQkZGBj788ENMnDgR/v7+kEgkKCoqQl5eHgQCAV577bU+NwTNkpKSEixduhT19fWYOnUqQkJCcPHiReTn58Pf3x9ZWVlmlzXvq+bMmYOioiIcPHgQ4eHhfIfjdFQ/ci22nutqamqwZMkSVFVVYfr06QgKCsLZs2fx008/ITk5GTt37jS6zmo0GqxcuRLHjx9HdHQ0EhISUFpaitzcXAwZMgTZ2dl0frCBLXVnZ5TVBx98gC1btsDPzw8zZ86ESqXC119/jba2Nnz44YeUKDHD2jI8deoUli9fjvHjxyMkJAReXl4oLy/X12vnz5+PN954wyhRQmXoeLa2y1QqFR599FFcunQJSUlJGDNmDAoLC3Hs2DGMHj0au3fvhlQqNfr7V155BdnZ2QgLC8OUKVNQX1+PQ4cOwcPDA1988QXCwsKM9j9w4AA2btwIb29vpKWlgWEYfPPNN6irq8PmzZvx0EMPOe4N6YVsKcOysjI8+OCDiIyMhEKhwKBBg1BdXY1jx46hoaEBSUlJ2L59u8nNS77L0OHJLkDb/Ts9PR3nz59HR0cHwsLCsHz5csydO9fRT90r6BrYlqxZs4Ym0sa93ycAeOONN/QTyfZVZ8+exb59+3D+/HnU1NTg9u3b8PLyQkREBObNm6c/cRBTlOwydfr0aezevRsXL15EXV0d2tvbMWjQIIwfPx5PPPEEoqKi+A7RpVRWViI9PR0//vgjGhoa4OPjg2nTpuHZZ5/FoEGD+A7PZRQUFGDx4sWIiooy6QnXl1D9yHV051xXU1ODLVu24OjRo2hsbISfnx/mzp2LlStXmu2t2N7ejszMTBw8eBAVFRXw8vLClClTsG7dOvj6+jrjZboNW+vOziirgwcPYteuXbhy5QrEYjFiYmKwdu1auk5aYG0ZVlZWYtu2bSgoKEBVVRVaWlogk8kwZswYLFq0CKmpqWb/nsrQsbrTLmtubsbWrVtx+PBh1NXVwcfHBzNnzsSaNWvMznnKsiw+++wzZGVlobi4GFKpFImJiVi/fr1RjyBDx44dQ2ZmJi5cuAAAGDNmDJ5++mlMnjy5B6/WPdlShkqlEu+88w7Onz+PiooKNDc3o3///lAoFJg7dy4WL15stlcY32XolGQXIYQQQgghhBBCCCHOQOOYCCGEEEIIIYQQQojboGQXIYQQQgghhBBCCHEblOwihBBCCCGEEEIIIW6Dkl2EEEIIIYQQQgghxG1QsosQQgghhBBCCCGEuA1KdhFCCCGEEEIIIYQQt0HJLkIIIYQQQgghhBDiNijZRQghhBBCCCGEEELcBiW7CCGEEEIIIYQQQojboGQXIYQQQgghhBBCCHEblOwihBBCCCGEEEIIIW6Dkl2EEEIIIYQQQgghxG38H8Ur0F+PaFzdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1200x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "az.plot_trace(trace, var_names=[\"γ\"]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>mc error</th>\n",
       "      <th>hpd 3%</th>\n",
       "      <th>hpd 97%</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>λ1</th>\n",
       "      <td>0.68</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.58</td>\n",
       "      <td>0.77</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>λ2</th>\n",
       "      <td>0.60</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.51</td>\n",
       "      <td>0.69</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>λ3</th>\n",
       "      <td>0.49</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.41</td>\n",
       "      <td>0.57</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>λ4</th>\n",
       "      <td>0.79</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.69</td>\n",
       "      <td>0.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>γ</th>\n",
       "      <td>1.58</td>\n",
       "      <td>0.80</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.34</td>\n",
       "      <td>3.09</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mean    sd  mc error  hpd 3%  hpd 97%\n",
       "λ1  0.68  0.05      0.00    0.58     0.77\n",
       "λ2  0.60  0.05      0.00    0.51     0.69\n",
       "λ3  0.49  0.04      0.00    0.41     0.57\n",
       "λ4  0.79  0.05      0.00    0.69     0.89\n",
       "γ   1.58  0.80      0.01    0.34     3.09"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "az.summary(trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximate $\\text{Prob}(\\lambda_i > \\lambda_j \\mid Y)$ for all pairs of $i$ and $j$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9244"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Is the rate higher in 2015 than 2012?\n",
    "np.mean(samples[:, 3] > samples[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0028"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Is the rate higher in 2014 than 2012?\n",
    "np.mean(samples[:, 2] > samples[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1376"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Is the rate higher in 2013 than 2012?\n",
    "np.mean(samples[:, 1] > samples[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.93"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Is the rate higher in 2015 than 2013?\n",
    "np.mean(samples[:, 4] > samples[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0444"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Is the rate higher in 2014 than 2013?\n",
    "np.mean(samples[:, 2] > samples[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Is the rate higher in 2015 than 2014?\n",
    "np.mean(samples[:, 3] > samples[:, 2])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}