07-partial_pooling_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Multilevel Modelling using PyStan\n",
    "\n",
    "This is a tutorial, following through Chris Fonnesbeck's [primer on using PyStan with Bayesian Multilevel Modelling](http://mc-stan.org/documentation/case-studies/radon.html).\n",
    "\n",
    "# 7. Partial Pooling - Introduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/site-packages/Cython/Distutils/old_build_ext.py:30: UserWarning: Cython.Distutils.old_build_ext does not properly handle dependencies and is deprecated.\n",
      "  \"Cython.Distutils.old_build_ext does not properly handle dependencies \"\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pystan\n",
    "import seaborn as sns\n",
    "\n",
    "import clean_data\n",
    "import unpooled_model\n",
    "\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pooling and Multilevel/Hierarchical Models\n",
    "\n",
    "### pooled model\n",
    "\n",
    "When we *pool* data, we imply that they are sampled from the same model. This ignores all variation (other than sampling variation) among the units being sampled. That is to say, observations $y_1, y_2, \\ldots, y_k$ share common parameter(s) $\\theta$:\n",
    "\n",
    "![pooled model](images/pooled_model.png)\n",
    "\n",
    "### unpooled model\n",
    "\n",
    "If we analyse our data with an *unpooled* model, we separate our data out into groups (which may be as extreme as one group per sample), which implies that the groups are sampled independently from separate models because the differences between sampling units are too great for them to be reasonably combined. That is to say, observations (or grouped observations) $y_1, y_2, \\ldots, y_k$ have independent parameters $\\theta_1, \\theta_2, \\ldots, \\theta_k$.\n",
    "\n",
    "![unpooled model](images/unpooled_model.png)\n",
    "\n",
    "### partial pooling/hierarchical modelling\n",
    "\n",
    "In a *hierarchical*, or *partial pooling* model, model parameters are instead viewed as a sample from a population distribution of parameters, so the unpooled model parameters $\\theta_1, \\theta_2, \\ldots, \\theta_k$ can be sampled from a single distribution $N(\\mu, \\sigma^2)$.\n",
    "\n",
    "![partial pooled model](images/partial_pooled_model.png)\n",
    "\n",
    "One of the great advantages of Bayesian modelling (as opposed to linear regression modelling) is the relative ease with which one can specify multilevel models and fit them using Hamiltonian Monte Carlo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partial Pooling\n",
    "\n",
    "### A simple model\n",
    "\n",
    "The simplest possible partial pooling model for the radon dataset is one that estimates radon levels, with no other predictors (i.e. ignoring the effect of floor). This is a compromise between pooled (mean of all counties) and unpooled (county-level means), and approximates a *weighted average* (by sample size) of unpooled county means, and the pooled mean:\n",
    "\n",
    "$$\\hat{\\alpha} \\approx \\frac{(n_j/\\sigma_y^2)\\bar{y}_j + (1/\\sigma_{\\alpha}^2)\\bar{y}}{(n_j/\\sigma_y^2) + (1/\\sigma_{\\alpha}^2)}$$\n",
    "\n",
    "* $\\hat{\\alpha}$ - partially-pooled estimate of radon level\n",
    "* $n_j$ - number of samples in county $j$\n",
    "* $\\bar{y}_j$ - estimated mean for county $j$\n",
    "* $\\sigma_y^2$ - s.e. of $\\bar{y}_j$, variability of the county mean\n",
    "* $\\bar{y}$ - pooled mean estimate for $\\alpha$\n",
    "* $\\sigma_{\\alpha}^2$ - s.e. of $\\bar{y}$\n",
    "\n",
    "### Specifying the model\n",
    "\n",
    "We can define this in `stan`, specifying `data`, `parameters`, `transformed parameters` and `model` blocks. The model is built up as follows.\n",
    "\n",
    "Our observed log(radon) measurements ($y$ approximate an intermediate transformed parameter $\\hat{y}$, which is normally distributed with variance $\\sigma_y^2$:\n",
    "\n",
    "$$y \\sim N(\\hat{y}, \\sigma_y^2)$$\n",
    "\n",
    "The transformed variable $\\hat{y}$ is the value of $\\alpha$ associated with the county $i$ ($i=1,\\ldots,N$) in which each household is found.\n",
    "\n",
    "$$\\hat{y} = {\\alpha_1, \\ldots, \\alpha_N}$$\n",
    "\n",
    "The value of $\\alpha$ for each county $i$, is Normally distributed with mean $10\\mu_{\\alpha}$ and variance $\\sigma_{\\alpha}^2$. That is, there is a common mean and variance underlying each of the prevailing radon levels in each county.\n",
    "\n",
    "$$\\alpha_i \\sim N(10\\mu_{\\alpha}, \\sigma_{\\alpha}^2), i = 1,\\ldots,N$$\n",
    "\n",
    "The value $\\mu_{\\alpha}$ is Normally distributed around 0, with unit variance:\n",
    "\n",
    "$$\\mu_{\\alpha} \\sim N(0, 1)$$\n",
    "\n",
    "In `data`:\n",
    "* `N` will be the number of samples (`int`)\n",
    "* `county` will be a list of `N` values from 1-85, specifying the county index each measurement\n",
    "* `y` will be a `vector` of log(radon) measurements, one per household/sample.\n",
    "\n",
    "We define parameters:\n",
    "\n",
    "* `a` (`vector`, one value per county), representing $\\alpha$, the vector of prevailing radon levels for each county.\n",
    "* `mu_a`, a `real` corresponding to $\\mu_{alpha}$, the mean radon level underlying the distribution from which the county levels are drawn.\n",
    "* `sigma_a` is $\\sigma_{\\alpha}$, the standard deviation of the radon level distribution underlying the county levels: variability of county means about the average.\n",
    "* `sigma_y` is $\\sigma_y$, the standard deviation of the measurement/sampling error: residual error of the observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "partial_pooling = \"\"\"\n",
    "data {\n",
    "  int<lower=0> N;\n",
    "  int<lower=1,upper=85> county[N];\n",
    "  vector[N] y;\n",
    "}\n",
    "parameters {\n",
    "  vector[85] a;\n",
    "  real mu_a;\n",
    "  real<lower=0,upper=100> sigma_a;\n",
    "  real<lower=0,upper=100> sigma_y;\n",
    "}\n",
    "transformed parameters {\n",
    "  vector[N] y_hat;\n",
    "  for(i in 1:N)\n",
    "    y_hat[i] <- a[county[i]];\n",
    "}\n",
    "model {\n",
    "  mu_a ~ normal(0, 1);\n",
    "  a ~ normal(10 * mu_a, sigma_a);\n",
    "  \n",
    "  y ~ normal(y_hat, sigma_y);\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We map Python variables onto the model data (remembering to offset counts/indices by 1, as Stan counts from 1, not from 0):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "partial_pool_data = {'N': len(clean_data.log_radon),\n",
    "                     'county': clean_data.county + 1,\n",
    "                     'y': clean_data.log_radon}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we fit the model, to estimate $\\mu_{\\alpha}$, and $\\alpha_i, i=1,\\ldots,N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "partial_pool_fit = pystan.stan(model_code=partial_pooling,\n",
    "                               data=partial_pool_data,\n",
    "                               iter=1000, chains=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're interested primarily in the county-level estimates of prevailing radon levels, so we obtain the sample estimates for `a`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sample_trace = partial_pool_fit['a']  \n",
    "means = sample_trace.mean(axis=0) # county-level estimates\n",
    "sd = sample_trace.std(axis=0)\n",
    "samples, counties = sample_trace.shape\n",
    "n_county = clean_data.srrs_mn.groupby('county')['idnum'].count()  # number of samples from each county"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to compare the results from our partially-pooled model to the unpooled model from `05-unpooled_model.ipynb`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz8AAAF+CAYAAACszqSjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xt81NWd//H39zsTEpIBk2gSbzRaq3jZbqvY/ryBFEHk\nkhUVBBWwxct662qrooDaKFpttK1bQUvtrnTpRYSqiJeKLspWvLFYXIstIipWEAiEAJOEkJnv+f0x\nzGQScpn79fV8PHiEZGa+3zOTyZzz+X7O+RzLGGMEAAAAADnOTncDAAAAACAVCH4AAAAA5AWCHwAA\nAAB5geAHAAAAQF4g+AEAAACQFwh+AAAAAOQFgh/E5Pjjj1djY2OHn7388suaMmVKRrSlN9dcc42e\nffbZJLWo3c6dO3X88cf3er8ZM2boiSeeSMg5Fy1apD/84Q+SpCeffFKPP/54Qo4rSVdccUXUrzUA\nRGvTpk068cQTdcEFF+iCCy7QuHHjNG7cOP3xj3+M+lgrVqzQL37xC0nS8uXLdd999/V67pNPPlmS\nNGfOHN17773RP4EEiaVvSGVfHElf+u6776qmpiYh5/viiy/0b//2b5Kkbdu26ZJLLknIcaWOfSdy\nmzvdDUB2siwrqp8nUzrOGSljTMrb99577+m4446TJE2aNCmhx165cmVCjwcA3SkqKtIzzzwT+n7r\n1q2qqanR17/+9dBnXCQ++OAD7d69W5I0bNgwDRs2rNfHZHK/Eolsb393Nm3apE8//VSSVFlZmdBg\nJbzvRG4j+EFMetsbd86cOdq0aZO2bdumzZs3q7y8XA8//LAqKio0bNgwDR8+XKtXr5bX69V3v/vd\n0NWbhQsX6re//a1cLpcOPvhg3XnnnTrqqKPk9Xp199136+9//7ssy9LgwYN18803y7btDm1ZvHix\nfv/730uSSktLdccdd+irX/2qtm3bpttvv1319fWqqqrSzp07u2z3jBkzVFhYqA8++EA7duzQeeed\np/Lyci1fvlw7duzQvffeq//3//5fj+1ZtmyZHn74YfXt21cnnXRSh+N3bt+dd96po48+utvXsa2t\nTQ899JBWrVolx3F0wgkn6I477lBJSYl+//vfa+HCherTp48KCwt1991369NPP9Xy5cv15ptvqrCw\nUA0NDWpsbNQdd9yhYcOGqaamRq+//rp27dqlG264Qe+9957Wrl2rgoICPfbYY6qoqNBrr72mefPm\nyefzqaGhQePGjdO//du/acaMGZKkqVOnhrJJs2fP1pdffimfz6cxY8bo6quvlt/v1+zZs/Xee++p\noKBAAwYM0P3336++ffv2+J4BgJ5UVVWpurpan332mQYMGKDa2lpt3LhRjY2NKikp0U9/+lMdddRR\nmjJlikpLS/Xpp5/qvPPO05NPPinHceTxeFRdXa2XX35Zv/zlL7VmzRo99NBDamtrU319vc4444xu\nszzvvfeefvjDH+r111+XJO3du1fDhg3T888/r/Ly8tD95syZo48++kgNDQ2qr6/XCSecoPvuu08l\nJSVav369Zs+ercbGRtm2re9+97saN26cpAP7vrvuukvV1dUd2rBhwwb9+Mc/VmNjoxzH0eTJk3XR\nRRdJkv793/9dzz//vMrKyjRgwIAun8O7776rn/3sZ6qsrNT69evVt29fff/739eCBQv02WefacSI\nEaHP+e764p760s7tmzJlii688MIef6evvfaaHnvsMfl8PhUVFWn69On65je/qU8++USzZs3Svn37\nZIzRhAkTNHHiRN15553atm2brrzySt19990aO3as/vKXv2jOnDn6/PPP9fnnn6u+vl7//M//rDPP\nPFPPPvusNm3apFtvvVWjR4/Wjh07dNddd2nHjh3avn27Dj/8cD388MN67733OvSdl156qX75y19q\n2bJlMsboiCOO0I9+9CNVVFRo2bJl+uUvfynbtuVyuXTrrbfq1FNP7fF5IsMYIAYDBw40O3fu7PCz\nP/3pT2bKlCnGGGMeeeQRM2LECNPU1GSMMeaaa64xjzzyiDHGmO985ztm1qxZxhhjtmzZYk477TTz\n0Ucfmbfeesuce+65oeM+/fTTZvTo0cYYY6ZPn27uu+8+Y4wx+/btM9OmTTO/+tWvOrTl3XffNZdd\ndpnZu3evMcaYN954I/T46667zvz7v/+7McaYzz//3Jx88snmmWeeOeB53X777WbixInG7/eb+vp6\nM3DgQPPb3/7WGGPMb37zGzNt2rQe27N9+3Zz6qmnmg0bNhhjjPn1r39tjj/+eGOMMe+880637bv9\n9tvNf/7nfx7Qnjlz5pi6urrQ9z/72c/M3Xffbfx+v/mnf/onU19fb4wxZsmSJeapp5464FiPPPKI\nmT17duh1f+CBB4wxxrzwwgvmhBNOMOvWrTPGGHP99debefPmGWOMmTp1qtm4caMxxpitW7eaE088\nMfQ7GThwoGlsbAzd77XXXjPGGNPa2mqmTp1qXnrpJbNq1SozatSoUJsfeugh85e//OWA5wYA3fni\niy/MySef3OFn7733nvn2t79ttmzZYv70pz+Ze++9N3TbXXfdFfqsmzx5cqiPMabj5+DTTz9t/vVf\n/9UYY8wPf/hD8+677xpjjGlqajKnnXaaWbt2bYdzhz923LhxZsWKFcYYYxYvXmx++MMfHtDuRx55\nxAwePNjs2LEjdI6f/OQnxufzmeHDh5tXXnnFGBP4bB0yZIhZs2ZNj31f8PPc5/OZMWPGmA8//NAY\nY8yePXvM6NGjzfvvv29effVVM3bsWNPc3Gz8fr+59tprQ31xuHfeececdNJJ5m9/+5sxxpgrr7zS\nTJo0yfh8PtPQ0GBOOukks23bNvPmm292257u+tKe2vfOO++YsWPHHtCezz77zIwdOzbUp6xfv96c\neeaZpqWlxcycOTPUx9fX14de6/Bjdf49nXPOOcbr9Zq9e/eab3/726H+7tVXXzUjR440xgT68ccf\nfzzUhquuuso88cQTHV5rY4x55plnzA9+8APj9/uNMcYsXLjQXHXVVcYYY4YPH27ef/99Y4wxK1eu\nNHPnzj3guSGzkflBTLpKqTuOI5fLFfr+29/+toqLiyVJJ554Yoe1IpdddpmkwJW8IUOG6I033tD2\n7ds1atQolZaWSpIuuOAC/fjHP9YXX3yhP//5z3ryySclSQUFBbrkkkv0m9/8RldddVXomK+//ro+\n//xzTZo0KZQN2r17t3bt2qW33npLt99+uyRpwIABOvPMM7t9bt/5zndk27YOOeQQ9e3bV4MHD5Yk\nfeUrX9GuXbskqdv2VFdXa+DAgfrqV78qSZo4caIeeughSYF55121Lzgdoyuvv/669uzZE5pu5vP5\ndPDBB8u2bY0aNUoTJ07U0KFDdeaZZ0Y0p/rcc88NPZeKiopQin/AgAGh389jjz2m119/Xc8995w+\n+eQTSVJLS0vo92KMUUtLi1atWqXdu3fr4YcfDt3nb3/7m6688kq5XC5NmDBBZ511lkaMGKF//ud/\n7rVtABBu7969uuCCC2SMkd/vV1lZmX7605+qqqpKI0eO1IABA/Tb3/5WGzdu1LvvvhtapyMpoivx\nDzzwgFasWKF58+bpk08+0d69e9Xc3KyDDjqoy/tfeumlWrRokYYMGaKFCxdq+vTpXd5v5MiRoWzQ\n+PHjdf/99+uiiy7Svn37NHz4cEmBKVsjR47U//zP/2jv3r1d9n2bNm0KHfOzzz7T559/rpkzZ4b6\nj9bWVn344Yf6+OOPNWLEiFB2ffz48Zo/f36XbTviiCNC61C/8pWvqF+/fnK5XCorK1O/fv20a9cu\nvfHGG932xd31pT21L9gfdrZy5Upt375d3/3ud0OPcbvd2rhxo0aMGKHbbrtN//d//6fTTz9ds2bN\n6vIY4c444wyVlJSEXt8hQ4aEnmew7546dar+93//V/Pnz9dnn32mjz/+WN/4xjcOONbrr7+uDz74\nIJS5chxHra2tkqQxY8bouuuu09ChQ3XGGWfoyiuv7LVtyCwEP4hJeXm5GhsbQx+OkrRjx44O3xcV\nFYX+3zlYCg+S/H6/XC6XHMc54DzBTs90mmbnOI58Pt8BPzv//PN18803h3725Zdf6qCDDpJtd6zt\nEX7+zvr06dPhe7f7wD+Tzm0Ntse27Q63hZ+nu/b179+/27b4/X7NmjUrFIA1NzeHPoDr6ur08ccf\n680339Tjjz+uxYsX69FHH+32WJ2fW1fPq6WlRePGjdO5556rU089VePHj9err77a4fW3LEt+v1+S\nQtPuJKmhoUF9+/ZV3759tWTJEr333nt6++239YMf/EBTp07V5Zdf3mPbACBc5zU/4X7/+99r0aJF\nmjx5smpqanTQQQd1CBaCF956cumll+qEE07QkCFDNGrUKL3//vs9TumuqanRz3/+c73zzjtqaWnp\nNsAK/2w1xoT6t+76sa7O2bmP8/v96t+/f4fXo76+Xv3799eDDz7Y4RiJ7t+Cz8Pv93fbl3bXvoMO\nOkhr1qzpsi2O4+j000/Xz372s9DPNm/erMMOO0wDBw7UsmXLtHLlSr399tuaO3du6IJjPM/twQcf\n1F//+ldddNFFOu2003p8/a+66qrQutm2trbQFL+bbrpJ48eP18qVK/XMM8/o8ccf7/Z9isxEtTfE\nZMiQIVqwYEHoQ2PXrl165plndPbZZ0f0+GB1mM2bN+vNN9/UkCFDNHjwYL300ktqaGiQJP3xj39U\naWmpqqurdeaZZ+p3v/udJGnfvn1auHDhAdmbM888Uy+88ILq6+slBQbm06ZNkyQNHjxYCxculCRt\n2bJFb731VlzP/6yzzuqyPaeccoo2bNigdevWSZKefvrpiNrXncGDB+t3v/ud9u3bJ8dxdNddd+nn\nP/+5du7cqaFDh6q0tFRTp07VTTfdFDqny+VSW1tbTM9r48aNam5u1k033aShQ4fqnXfeUVtbWyjY\ncbvdamtrk8fj0Te+8Q39x3/8hyRpz549mjx5sv77v/9br7/+ui6//HKdfPLJuuGGGzRu3Dj9/e9/\nj6k9APJXT4HIypUrdeGFF+qiiy7SUUcdpddee63LQbvU9Wfi7t279eGHH+rWW2/V8OHDtWXLFn3+\n+eehz7quFBUVqaamRjNnzuyxmMzy5cvl9XrlOI6eeuopDRs2TEcffbT69OmjV199VVKgeMPLL7+s\nM888U2edddYBfV9ZWVmHNT9HH320CgsL9dxzz0kK9GPnn3++PvzwQw0ePFh/+tOftGfPHjmOE7pP\ntIKvd099cXd9aXftW7t2bbfnO+2007Ry5crQDIM///nPuuCCC9Ta2qqbb75ZL7zwgkaPHq277rpL\nHo9HW7ZskcvlOuDCZzRWrlypyy+/XP/yL/+isrIyvfnmm6H3Tfj75KyzztKiRYvk9XolSY888ohm\nzJghv9+vYcOGqbm5WRMnTtSPfvQjffLJJzH3uUgPMj+IycyZM/XAAw9o7NixcrvdMsaEypFG4ssv\nv9SFF16o1tZWzZo1S0cddZSOOuooXX755aEMQVlZmebNmydJuuOOOzR79mzV1NSora1NgwcP1jXX\nXCOpPat01lln6corr9S0adNk27Y8Ho/mzJkjSbrzzjs1c+ZMjRkzRoceemhE5afDj91Z5/YMGTJE\n11xzjdxutx566CHdcsst6tOnj771rW+FHtNT+7pz3XXXqa6uLjT144QTTtBtt92mkpISXXfddbr8\n8stVWFiogoKCUPnWIUOGaPbs2RE/l3DHH3+8zj77bI0aNUqVlZU65ZRTdNJJJ+nzzz/XgAEDNHz4\ncF166aV69NFH9dOf/lT33HOPampq5PP5NHbsWI0dO1aO4+jPf/6zxo4dq+LiYpWWlnbZHgDoSU+f\nWdOmTdNdd92lZ599VqWlpRo+fLhWrFjR5eNOP/10ff/731dBQUGoCE3//v119dVXa9y4caqqqtLX\nvvY1DRkyJPRZ150LL7xQTz31lM4///xu73PIIYfo6quvVkNDg0499VT967/+q9xut+bOnat7771X\nv/jFL+Q4jr7//e/r29/+tiR12/cFFRQU6NFHH9W9996rX//61/L7/brppptCU/3Wr1+viy66SAcd\ndFDE/VtnwdftjDPO6LY93fWlPbXv3Xff7fJ8X/va13TPPffohz/8oaRA8PHYY4+pqKhI1113ne64\n4w499dRTsm1b5557rr71rW9p9+7dsm1bF198cYeMUaSuv/56/eQnP9G8efNUXl6ukSNHauPGjZI6\n9p1XX321tm7dqokTJ8q2bR122GG6//775XK5NGvWLN18880qKCiQbdu6//77VVBQEHVbkD6W6enS\nCpAEw4YN08MPP8w6EABAVvnVr36lL7/8Uj/60Y+6vH3OnDnavn27amtrU9swABGLOfPjOI7uuOMO\nffrpp7JtW3fffbe+9rWvhW6fP3++Fi9eHFr0d8899+ioo46Ku8HIfrm6/wAAIHedc845Ovjgg/XY\nY4+luykA4hBz5ufVV1/Va6+9pvvuu0/vvvuu5s+f32Gx9a233qrvfe97OvHEExPWWAAAAACIVcyZ\nn+HDh4d2Sd60adMBpSHXrl2refPmqb6+XkOHDtXVV18dX0sBAAAAIA5xFTywbVszZszQK6+8ol/8\n4hcdbhszZowuu+wyeTweXX/99VqxYkWPlcCMMUyHAgBkLPopAMh+CSl4sGPHDk2YMEEvvvhiaG8X\nr9crj8cjKVCPf9euXbr22mt7PE59/Z54mwKkVUVFP97HyGoVFf3S3YSMxt83sh39FLJdvP1UzPv8\nPPvss6HSh4WFhbJtO7T5ldfrVU1NjVpaWmSM0dtvvx0qLwkAAAAA6RBz5mfv3r26/fbbtX37dvl8\nPl199dVqbm5WS0uLJkyYoBdeeEFPPPGECgsLdfrpp+uGG27o9ZhciUC244oash2Zn57x941sRz+F\nbBdvP5VR+/zwx4hsR6eCbEfw0zP+vpHt6KeQ7dI27Q0AAAAAsgnBDwAAAIC8QPADAAAAIC8Q/AAA\nAADICwQ/AAAAAPICwQ8AAACAvEDwAwAAACAvEPwAAAAAyAsEPwAAAADyAsEPAAAAgLxA8AMAAAAg\nLxD8AAAAAMgLBD9JUFtbqNrawnQ3AwAAAEAYgp8kWLrUraVL3eluBgAAAIAwBD8AAAAA8gLBDwAA\nAIC8QPADAAAAIC8Q/AAAAADICwQ/AAAAAPICwQ8AAACAvEDwAwAAACAvEPwAAAAAyAsEPwAAAADy\nAsEPAAAAgLxA8AMAAAAgLxD8AAAAAMgLBD8AAAAA8gLBDwAAAIC8QPADAAAAIC8Q/AAAAADICwQ/\nAAAAAPJC1gU/tbWFqq0tTHczAAAAAGSZrAt+li51a+lSd7qbAQAAACDLZF3wAwAAAACxiDn4cRxH\nM2fO1CWXXKLLLrtMH3/8cYfbly9frvHjx2vSpElatGhR3A0FAAAAgHjEHPwsX75clmXpD3/4g268\n8Ub97Gc/C93m8/n0wAMPaP78+VqwYIEWLlyohoaGhDQYAAAAAGIRc/AzfPhwzZ49W5K0adMmHXTQ\nQaHbNmzYoOrqank8HhUUFGjQoEFatWpV/K0FAAAAgBjFVTnAtm3NmDFDr7zyin7xi1+Efu71etWv\nX7/Q9yUlJdqzZ0+vx6uo6NfrfWw78vumSza0EcnD7x3IXfx9IxfwPkY+i7ts2v33369bbrlFEyZM\n0IsvvqiioiJ5PB55vd7QfZqamtS/f/9ej1Vf33uA5Dgl++/bFHujkywb2ojkqKjoF9H7GMhUDIp6\nxt83sh39FLJdvP1UzNPenn32Wc2bN0+SVFhYKNu2Ze9PeRxzzDHauHGjdu/erX379mnVqlX65je/\nGVdDAQAAACAeMWd+zjvvPN1+++2aPHmyfD6fZs6cqWXLlqmlpUUTJkzQjBkzNG3aNBljNGHCBFVW\nViay3QAAAAAQlZiDn6KiIj388MPd3j506FANHTo01sNnNceR2tokr1fyeFJ77trawv1fW1N7YgAA\nACDDsclpgnm90tatlrZtszVyZLHClj6lxNKlbi1dGvdSLgAAACDnEPwk2Lp1tnw+S5K0fr1L69bx\nEgMAAACZgJF5gg0c6MjtNpKkY4/1a+BAJ80tAgAAACAloNQ1OvJ4pKoqo7Y2o5dfbk75mh8AAAAA\nXSPzkwS2LRUWpr7YAQAAAIDuEfwAAAAAyAsEPwAAAADyAsEPAAAAgLxA8AMAAAAgLxD8AAAAAMgL\nBD8AAAAA8gLBDwAAAIC8QPADAAAAIC8Q/AAAAADICwQ/AAAAAPICwQ8AAACAvEDwAwAAACAvEPwA\nAAAAyAsEPwAAAADyAsFPjnEcqbVV8nrT3RIAAAAgsxD85BCvV9q61dK2bbZGjiwmAAIAAADCEPzk\nkHXrbPl8liRp/XqX1q3j1wsAAAAEMTrOIQMHOnK7jSTp2GP9GjjQSXOLAAAAgMzhTncDkDgej1RV\nZdTWZvTyy83yeNLdIgAAACBzkPnJMbYtFRaKwAcAAADohOAHAAAAQF4g+AEAAACQFwh+AAAAAOQF\ngh8AAAAAeYHgBwAAAEBeIPgBAAAAkBfY5ycJamp86W4CAAAAgE4IfpKgtrY13U0AAAAA0ElMwY/P\n59PMmTO1adMmtbW16ZprrtGwYcNCt8+fP1+LFy9WeXm5JOmee+7RUUcdlZAGAwAAAEAsYgp+nnvu\nOZWVlamurk67du3SuHHjOgQ/a9euVV1dnU488cSENRQAAAAA4hFT8DNq1Cidd955kiTHceR2dzzM\n2rVrNW/ePNXX12vo0KG6+uqr428pAAAAAMQhpuCnb9++kiSv16sbb7xRP/jBDzrcPmbMGF122WXy\neDy6/vrrtWLFCp199tnxtxZpVVtbuP8ra5oAAACQfSxjjInlgV9++aVuuOEGTZ48WRdccEGH27xe\nrzwejyTp97//vXbt2qVrr702/tZKCi4d+uyzhBwu5yTz9eG1BwAAQDaLKfOzfft2XXHFFbrrrrt0\n2mmndbjN6/WqpqZGL774ooqKivT2229r/PjxER23vn5Pr/dxnJL9922KvuF5IJmvD6997yoq+kX0\nPgYyVUVFv3Q3IaPx941sRz+FbBdvPxVT8DNv3jzt3r1bjz76qObOnSvLsnTxxRerpaVFEyZM0C23\n3KIpU6aosLBQp59+uoYMGRJXIwEAAAAgXjFPe0uGSK5EDBoUyD6sXk32oSvJfH147XvHFTVkOzI/\nPePvG9mOfgrZLi2ZH2SumhpfupsAAAAAZCSCnxxDJTYAAACga3a6GwAAAAAAqUDwAwAAACAvEPwA\nAAAAyAsEPwAAAADyAsEPAAAAgLxA8AMAAAAgLxD8ICK1tYVqbLTS3QwAAAAgZgQ/iMjSpW41N0ut\nrZLXm+7WAAAAANEj+EFEHEfy+6Vt22yNHFlMAAQAAICsQ/CDiLS1SVJg2tv69S6tW8dbBwAAANmF\nESwiUlAgSUaSdOyxfg0c6KS1PQAAAEC0si74cRzWnaSDbUsul1RZ6ejll5vl8aS7RQAAAEB0sir4\n8XqlrVst1p2kUWGhCHwAAACQlbIq+Fm3zpbPx7oTAAAAANHLquhh4EBHbjfrTtKhsdGSw8sNAACA\nLOZOdwOi4fFIVVVGbW2GdScp1twsGZPuVgAAAACxy6rMjxRYeM+6EwAAAADRyrrgBwAAAABiQfCT\ng2prC1VbW5i1xwcAAACSgeAnTpkYCCxd6tbSpclZzuU48R8/E1+z3mRjmwEAANARwU+ckhloZAqv\nV/L7JcnS1q3xV33LxtcsG9sMAACAjjIq+OHqemYK7KcU2F/J57PU1hYofc3vCgAAANkko4Ifrq5n\npsB+SoE61263UUFBoPQ1vysAAABkk4wKfpCZPB7J5ZIko6oqI5t3DQAAALIQw1hEzLKk88/3Jex4\nTHMEAABAKhH8IGK2LdXWtibseNk0zdFxpNbWQPEHAAAAZCeCH6AXXq+0daulbdtsjRxZTAAEAACQ\npQh+gF6sW2fL5wtUu1u/3rW/+h0AAACyDaM4RMyY/Jz2NXCgI7c7UO3u2GP9+6vfpQ9rpQAAAGJD\n8JMH4h0sz5pVKL9fchxLI0cWx73JabbxeKSqKqPKSkcvv9wsjye97cmmtVIAAACZhOAnD8Q7WF6y\nxK3gJqfr17vU1paghmWQ3gJE25YKC5X2wAcAAACxI/hBrwoKpOAmp8ce69//feyCldMyKYOUTdmU\nxkZLjY1WupsBAACQdWIKfnw+n6ZPn67LLrtMF198sZYvX97h9uXLl2v8+PGaNGmSFi1alJCGIn3a\nNzU1evnl5rg2OQ2vnLZ1q5VRAVC2aG6Wmpryc/0VAABAPGK61P3cc8+prKxMdXV12rVrl8aNG6dh\nw4ZJCgRGDzzwgJ5++mkVFhbqkksu0TnnnKPy8vKENhzp4fEEMjbGxJa5Ca+c5vNZamszCW5hbvN6\nJb9fkgLrrzJhDRIAAEC2iOka/qhRo3TjjTdKkhzHkdvdHkNt2LBB1dXV8ng8Kigo0KBBg7Rq1arE\ntBZpF8zcOI6lLVssbd0a3ePDK6e53SbuKXSpUlPjU02NL93N0K23Fip8/RVltwEAACIXU+anb9++\nkiSv16sbb7xRP/jBD0K3eb1e9evXL/R9SUmJ9uzZE2czkQksq2Pmxu+3NHp0sVasiDz7EKyc1tZm\n5HZLbW2BgCrTsxe1ta3pboIk6a233Aqsv7JUUGB05JHMGwQAAIhUzCu8v/zyS91www2aPHmyRo8e\nHfq5x+ORN2wxQlNTk/r37x/RMe39i0kqKvr1cB/1ep9UyrT2SAe2Kd42Bh9v29JZZ5XI5QpOvZL+\n8Q+Xtm3rp6OPjvx4brfkckmbN0s+nzR6dD+tWpXeACiRv8dkvhfCpxq2tVlqauqnioqknQ5AJ5n0\nWQ/Eivcx8llMwc/27dt1xRVX6K677tJpp53W4bZjjjlGGzdu1O7du1VUVKRVq1bpiiuuiOi4zv6R\nXX19Uw/3Ken1PqmUae2RDmxTQ4Nn//exrZAPHC+Q7Wlp8aqqqkSbN1uSLB17rF+Vlc2qr4/ueK2t\nks8XiDgBPaOZAAAgAElEQVT+/nfpjTeaNGhQ+rIYPl+J2tqkTz9tiisIq6jop/r65GU6Xa7g7yK2\n1x7oDYOiniXz7xtIhWT3U0CyxdtPxRT8zJs3T7t379ajjz6quXPnyrIsXXzxxWppadGECRM0Y8YM\nTZs2TcYYTZgwQZWVlXE1EtFxnI7TyZqbE3t8lyvw7+CDY9/0s6AgsObH5wsM4gcOTF/gE1zH5PNl\nRxEBy5KKikzGtxMAACDTxBT8zJo1S7Nmzer29qFDh2ro0KGxtglx6GognyzxbPpp2+1rf5I5iA9u\nXNrTmp3wdUzBIgLpzEJ1J/i7NYYqeQAAALGgVFSO6TyQD1QHS6zGxsTsz2Pb8QVQkYhk89LwCnTp\nzkL1pHOZcCq9AYhHbW1h6AJRPpwXACSCn5zTuZT022/HXNOiW15vYJ+fXBGsQFdZGfs0vt4korPP\nliANQHaI5OJQNp2XoApAJFL/qYekCi8lXVDQXsUskeIJfGprC9XYaKm0NPqDRDKFLVbJzkIFO/p4\n2h7+u2W9DwB0tHSpW42Ngex4pmxPkK+S2V8D8SLzk4OCA/lkBD7xWrrUHXMBhnRdpcwkqZgqCADZ\nqrlZedtPZFLmi/4amYx3JnoVXN9jTGDKG9KnpsaX7iYAADJQImYYAPmA4Ac9ClYYkwJFDkaMKI77\nmMXFgUE8V4W6192UATo1AOhaoorxSEzbAnJZBk6MQiYJrzAmSRs2uOI+ZmmpSVmH4jhSa2vvGaua\nGl9SsyqRtiOIKQMAEJ3m5sQV4+EzOHbBtb1ApiL4QY/CK4xJ0jHH+NPYmugEs1bbttkaObK4x8Cj\ntrY1aQFZNO0AAKArmbSmpyfxrO0FUoHgJ8fFOw0gWGHMto0qKhy98kr2fKJ1tXlptrQj2kwRACC3\nZVM2ynFE9gcZKzv+ijKY40htbYFBaiZW4ErENADblixLKipK7HNM9uL9YNbK57Pi2hcn3rnf0bYj\nmCny+SyNHFlMWWsAQNIlcp2TMSL7g4xF5icO2TSdKRMrtSVzqpmUuM1L473aFm07MiVjBQD5Llum\nmiVCNmWWgHgwqopDNg1SHcdKWIBm7c9kp7rscixTwTJlX5xo2hG+ziqejBUAID4EBEDuydzRehbo\nbpCaSVeKAut9Ehug2bbkcqW2BGg2ZdnilaiMFQBksmT2lZnUDwPILAQ/cehukJoJV4oaGy01Nlr7\n1/tkThYh1oX82ZRl60q0pbRtW2ppsfTQQ3TeQLZh4B2ZZPaVmdAPR4r3C5Ba2TWCzECZMq2qs+bm\njosNbdskNIsQyxqieLI32T4VLJb1Tc3NyprOG0C7bBp4Z4p8DgB4vwCpRfCTJywrMQGa1yv5/YE1\nRCNGRBfAxJO9SedUsHSUna6p8am4OHXnA4B0SnQAEJz9AACdEfzkiURVe1uzxlZwDdGGDa7930cm\n3uxNOrJs6VprVFvbqtLSBG1VDgB5pvPsB6SO3+/s/2p01VVPa+fOxjS3COiI4CdPJLLaW6zSkb2J\ndq1NZ+lca8RGpwCASNTWFmZMpmvnzr37/2dryZKpmj79tbS2B+iMSaZ5oePgfdCg2NfLfPObjgIF\nFCwdc4x///eRS3X2JtJ1Nt1t7paojVKjFb7R6SmnlOi995pCr1kiN6IDAHQUDCSyKfu+dKlbzc1S\naWm6WyL5fO0zRCRLGzf2T2dzgAOQ+ckLiSsU4PEEylzbttErr+ROGebu5puna61ReMapsdHukHFi\ncSwA9M5xFFMRhWAgkU38fkeOY7Rly960TjVraGiU4+xTYNzhSNqp6urdaWkL0B2Cnxzm90vBqy+J\nrvaWqAIK2WD3bkstLVZKn28gSA0ErW63ybrqdgAyR75UUnM6fUwakz8VM3fu3CtjLLW1Ffc41ay2\ntlC33pqYc9bWFurcc4t17rnFoffXbbe9Jr//IAXGHpYOP/wx1dV9JzEnBBKE4CeJ0tXh1NT4VFgo\nhaedpdiDlZoan77+dSeutTPxSHfHnY6Fs8EMm2RUVWU6/O5YCwQgGvmQLW5oaNy/r520fXuWpW0S\nINKpZkuXurVoUWLOuXSpWx98YOuDD+zQ+ytw3vZ2VFaeqLKyDJiLB4Qh+EmidHU4tbWtKiszCmYO\nJCMrjnWQtbWtWrasOW1rTGJ5HdMRMCXyilq4trb2QCdd1ecAIJPddlsw02GppaVEfn/2rNdJBLe7\nfbaAZNI21ay6eleHdmzb9iHV3pBxCH6SrLHRSkvWwrbbMwfJEOtc6nhEU7ktkYFncbEi2nMnkVfU\nwoUHOumsPgcgdzQ0NGr79uaUrxEJP+/27c1yOs9Vi1HnjEO6pGumQnl5kSzLqKCgWeefvyBtU83q\n6obJ5QoGQEabN19LtTdkHEZOSdbcnJtzjtMxl7q2tjVp2aeeppKVlpqUV/3xeoPz1zsGOvHulQQA\nUiBT0tJS0usakWSet6WlRFu3JuYzrHPGIR7hxQOiDdCSNeOjt2DVtm3ZtqVDDy3S449fkLapZmVl\npbLtPgr0XbakMqq9IeMQ/CDv9TaVLB1rbNats2VM+9XLAQMcDRzopK36HIDc0jlTkqoBaufz+nyJ\nCRTq6obt/19wynfs2Z/w4gEtLSVqaNjb+4OSLNnBaiIzgZkyBQ/oDsEPUireTUeToaepZOlaYxOe\n4XG5jF58sX2fn1TvlQSgXboLsCRK50xJqgaoB2ZorNC+OvFs0mk6JXssK/bsT+fiAYHv0yvZwWoi\ng6vy8qL9/3PSOgUP6E76/6KRApb8/syoDpbMqWux6mkqWTRrbBKZIQpmeGzb6NBDjaqq4j8mgPjl\nSuW0urph6tu3KeVrRMLPGwxQgvvqxFNVM7zgQbxrfjpnLgLfp1eyg9V4givHUYfA1bYD/aTLZaV1\nCh7QHYKfHBa4Etb+Ycbi+K71NJUs0jU2wQzRli1KWIbItgP7Kdmdfm2ZmD0DkF3Kykp1yCHFKirq\nqyOOmJSyAWrwvIceWiTbTlxhgs6D986ZoGgEiwdIjizLhGUy0ifZwWo8wZUxqd8OAogHo+EcFihv\n3f5hlk+L46PNwnQ3lSzSNTbJqsJWXKwDAp1MzJ4ByE7NzdKCBQVZP5Wvq+l0sQoWD3C5AsfYtm1f\nSividSU8aExGNiVdmUAgHQh+8oKRy5W9a0SiDWQSvU4nkjU2yajCVlPj05QpbQQ6AJIqF6qShhc8\niGe9TzjHMaHCB6msiJcOyQ6ugExC8IOoFRcrrk1TpcgXDccSyKRjL5xghujQQ5WwKmxkeACkUjYX\ncwgfrCdqOl3nqePJqIgXTZW1eItCAAgg+EmRVHcqjhNe/Sa+gged215aag5YhxKtSBcNxxLIpGsv\nHKqwAchmqSzmENxLR5K2b0/egpEtW2LfaLzz1PFkVMSLpspavEUhUsmyItscHEgHgp8USWWnEsyW\nOE541ZvYCx50bntNjS9lH2qxBDLshZN5svmKMpCLguWl0yW4l45kqaWlRH5/cs7T1iYtWRJb32fb\nlizLJHUdTCJKWAezR36/0fbtzSlbmxQewDqO6XBe21bKNwcHIhVX8PP+++9rypQpB/x8/vz5Gjt2\nrKZOnaqpU6fqs88+i+c0iFJ4tiTIshJX8KC2tjVlH2qxBjLRZmF6qqBGdbX45Up5YCBXBMtLp0vn\nvXQSKTggDx43lk1KHSfwz7atpK6DSUQJ62D2SLLV0lKSkrVJDQ2N2rrVG5phYoz0ne/8V1qLQsSL\ni3T5I+bRyK9//WstWbJEJSUlB9y2du1a1dXV6cQTT4yrcYjNwIGB8pyBq2qBT6bDDzdZmwFJxXSy\nntbWRLruJpAR65OoJuUUxwlcgfV6mRYIILCXTltbsCpbYi+mdV6rE7gYGN054imVHY26umFatqxJ\nPp+t0aMXx5RdijR7FMgQWfL7bV111dOqqxsWc0B3222vye+fovbX1dbmzTM1ffoCSZNjOma6BS/Q\nsdY298Wc+amurtbcuXO7vG3t2rWaN2+eLr30Uv3qV7+KuXGIjccT3Bum/dM73jU66F1tbasefDDd\nrcg8ia6+ByC7NTZacrmCc6fb+ym/3ySwpHR7NsVx2hJwvORIRJW1SLNHwQzRvn3xV6/rHHAFv27c\n2D+lU+OBWMQ8JB4xYoRcLleXt40ZM0Z33323/uu//kurV6/WihUrYm4ggOyWjup7QLpEU70rXzU3\nSy0tXU11sxNcUjoQENh2aqbcpmsdVXCPHttu09e+9mG32aNErC8KOnBfpcDX6urdKZ0aD8QiKZ8I\nl19+uTz757acffbZ+vDDD3X22Wf3+jh7f3qioqJft/eZOFG93ieVghmV8PZ0zrKEl91MR7stS5o4\n0Yr53F09R8uS9u2T+vbtF9M0pq6OmYj7xvOYRMmU92amOOssye2WfD7p+OOls84qYeobslZv/dQN\nNzyvlpZJkiwtWTJVhYVPauHCS1LYwgNF00/t3CnV1fVLcRY7fFNSS5s3l0X9ObpjR1dBZmC6m+P4\n5HJJ5eW9Z1V27GjUjh3Nktqn9Nu2HVF7XnghENiVl0fX30bbX3W+X0VFP1VWBv6/fv1J3T7uuOOa\ntWZN+1TD445r6fI9EUk7/vM/z9eAAU1qbg6meBxdfPFCPfbYv6i8vJ+qqqJ7TpkgneMGpFbcwY/p\nNDHW6/WqpqZGL774ooqKivT2229r/PjxER3LcQIL8uvrm7q9z/Tp2n+f2NqbaI4T+IAMb3PwZwFW\n2OLLnp9bYrWPLm1bmj7dG/Nr5jglamy0dP31gQ03vV5p82aPfD5Lp5zij6miWlevWyLuG89jEqGi\nop/q6/ek9JzZ4HvfK1RDg/Tgg61qaZFaWtLdInSHjr9nvfVTH33UV+ED+Y8+6pv2z4Tu+ykr7PtA\nP9XUZGnhQqPp05P52dlVh9E+KD/88J1Rv2ZnnLFJ0kHqWDwh8Jz8/n6aNu13evzxC3o9zlVXPafm\n5ilqD5wC/yLvq6yI7y8FskUNDQUqLY3kMYHXrb7+wLnDDQ3d3xY0e/ZgLV3aJL/f1qhRizV79nc6\nvM7R9ZsuHXxwsVpbA1tpuFyW5swZK79fqq/fo61b09MHxyNd4wZEL95+Ku75J9b+3S6ff/55LVq0\nSB6PR7fccoumTJmiyZMn67jjjtOQIUPiPQ3SLHwH8HinMaViagAV2jLLffe16rHHWsn4IOclonpX\n/nJiLim9cePR6moNSuCfK+IpXp2nhkVT+MDvd+T3S5s2+SOe8pioqnuR7AEUXF902GHFCa9el9j1\nWkByxZX5OeKII/Tkk09KksaOHRv6+ZgxYzRmzJj4Woa4WUmKL4J77/h8ltzu6EtoBz/sSxNfNTSE\nai0A0iER1bsSKVjhy+eLvMJXY6O1v+xvaj9HXS4rouxMV7quHtf+faRBaHX1rg5Tw6LpR3fu3Cup\nRI5ToCVLpkpaEPPzyRaBjKGl4HqtfHjOyH6sPM5htp2cKm/BvXds26iqKvYS2q2tiqjyF1kcANki\nEdW7EilY4autLfIKX+GZ/mxRXl60/3+dUzVGffs2RRyEBosHtBdLiDz66bx3UTwFBbJF57Li+fCc\nkf2y69MtA3U1KA/uaVJQkIYGpUBtbaF277ZkWdEHV7W1hdq5MzBHOFj6uLc1Q2RxACA2iazwlcns\nHjqjQw4pVllZZPPXgsHrP/4RfRs6Z5+ybcpjLBcZLSu49js7nzPyE8FPnDoPzIN7mgSnhKW+Pcnf\nnTieOcrtj+24ZmjQoOimzgEAetd5Gld+DE4Tv3FqJMrLi7RpkxRYu/TbhE55rK0tlG1LRUW93zf2\nc0R/odG2Axczk/GckyE4RuKiauRy8TXLqGlvjhP5VKhMFV4MwOeLbrFkIixd6s746QqBOdSBF+bY\nY/1RrxkCAEQmOI2roKA55mICydBeBTXRi1M7Fz2Ij99vtH17c0QL+YPZp+DapUROeVy61C3HUcbu\nn5OM55wM3Y2Rsnl/rsD6vORd+M6GcWW0Mib4yZVd4IPFACTJ7Y5usWS+cJxAAFRZ6cRUJhvIFbNm\nFeraawuz9vMOmS/T1iAFGRPsHNsH84mpGNZ548142WppKUngxqvJ4zipmf2Ri2JZG5cpcjE4SbaM\nCX7WrlVO7AIfLAZQWemoqiozr9BEy3ECixqdBCVojAn8KywUgU8Pkn01B71LZnDi9UpPPFGgP/6x\nT1Zf8AESx07QwDNRswmSv1aqsdGKqG/1+x05jukxM2FM9hWqyBT5sjYOARkTYZx0kkIZk2yfCmXb\ngYF9MiqtRau4OPAvVsGMnONY2rrVYoCWQtl2NSfXgrVkByfx7peF1NuxI3unxmSS4uLwrRgO3Jsn\n2oFn573jXK6OUy6cmK/cJX+/puZmRTQ9fufOvTLGyrjMRKIuiqYb+3Pll4zpbcMzJkyFSpzSUhPX\nHOHOa5gYoKE7S5a4tXixO2cC5GQHJwMHOiotDYwcsv2CT7647rqXsnZqTFeKiowsqy20riX2ICE6\npaUm7OJg52lq0Q88Oxfh6byeqKFhb1TH67gJtxNVqexkydQy2qle15wsmbo2DsmRUZeVgxkTj0ca\nNOifurzP6tV/7fLnmXb/zZs3SbLl9x+qzZu3KJiGr62t2v+1Y9WMRLVH+jSq+/d2fMcplvR/kgok\ntXU7QIv2+FLgNRo06Kyo2hPp8dN1/88/35iw42/e/KaMcevkk4fKtptjak+q7h9e5XDkyGI1N/9T\nqM3paE8i7u84xbKsN2VMmdzuj3XllTX6y1/eTWh7SkqK1afPcWpu/khnn92cEa9Pd+9hSJ9+6lH4\nAPTFF2/RoEF3dnXPLh+frvd5d+1pbv5CxtiSDlVLS6M2bXIkHS6Xy5WW9kg+9e17mP73f4s6vK69\nHX/z5jfk9x8mKdDuznvPBAIHf8TtDw+kLMsrt/sPGj78x722J9DvHym/v+O5emq/379BkqtDf9jV\n/QNltH2Stkg6S+vWNWnQoOJej39gG9/QkUcOiPj+PR2/sfFxSd+V1KfDc86Ez7FY7r9v3xuS1O0m\nrZna/s2b39j/fffvn1iP7/c7+vLLlTJGqq4+ReXlffWXv3yY0PbHcv94+yku4yeNLb+/UpJLfn+V\ngh/KyZ7K1NhodbpqFR/bbpbLtU2SXy7XtgRn5Hj7dcdxiuX3V8pxqrR169L9QWjm6pwlaWs7Ls0t\nip9tN+vww89QZeX5qqqq6TKYS8Q5CgvXJOXYSLyjj96j8AxFOrYzSKRMvGp/yCHFPe7ZE4n2iqKB\nf21tfl1++ZMxZrYa1b//vXG1p3dGxcUv9niP9k1cA5mo9u/Tp7T0x/vHBx1lc+U0dBSYbmkrUPjD\nE3UWNVNZxmTOx99XvhL4YFq9uinNLYnPoEElam0NbOIZ5HIZHXpoe6o/Wc/xsMM6RieHH27iPt+g\nQSXavNnS4YcbrV7dFPo+ePxojn3yySXatMmSFNgH6aOPvDk1xbGiop/q6/fEfZzVq22NGlUS+v6l\nl5oyei8kr1c67jiPfD5Lxx7rZ+pqFquo6JfuJmSshoZGHXmkWz6frdGjF6uu7jtdVk8bNCjwt5sp\nfVlX7Rk0qETbtzerpaVEgSyJX5JXUj9ZlnT44Y7+8peWpLVn8+bg/jDtXC7pyy+jnzfb+XiWZcKq\nyUnBfX/OP39Bt1f1w1VVlex/fHSPGzSoRP/4hxXV8wj22b3dv3M/3Nkxx3jk9UpHHtlzn3zMMR5Z\nlqWPP46/nwq264svLHk80oYNgedw1VXPaMmSKYr29UuFnv42M+3vNlLJbPe55/631qwZF/r+m998\nVsuWnZPw80Qr3n6KS+9JUlAQCHiC/H5LbW2pb0dNjS+mXZuTJfAasIaoN+El07NhPQhr9pAPyssz\ns2x0rMrLi+Ry7VJgWrZXUn9Jtoyxsu4Kb/g6nwMDn8DXSNfJdJ42lynra3rScR1V95qbpaYuxsjx\nFKyx7Y77D1E5LXUFgJKdZcvVQhCMPJPEtrW/1HX7nj8FBak7f7DKW21ta8J25U3EBrSB16D9Ncn0\nQX26ZGMwcf75Po0f78uKtgIIbMp56KEeWZZkWS4duE4mtRzHxDx46xywhN0S+hrpwC18I+5IHxcc\nhErxPY9YBM8dyYasxcVSScmBP0/klPxcHTBHI1XVWpO9P1GuFoIg+EkilyvwL7jnTypLX/fvb9S3\nr0lI5S3Hkfx+hTagdZzAB2jfvvEEREZVVYaBcg/CC4Bkg0QG2gCSJ/xqcUPDXtm2paIiqeN6plRf\nmApkbG666eXYHt0pYOl42x6NGvUfEQ/cbLt9A9ZIK70FB6HB55HKSoDt5+59Q9bSUqOysuiOH20W\nI1MHzLm4FinZWbZM3SQ5XgQ/KZCOPX+2brVCwUq8AVD4VLXAYvZAcLV3r6I6R7AiWMercpkv1/av\nAZDfwq8Wt7SUyHGMysuL1LdvkyRHlmWStqA+PEvRrn3w9tZbsXWW4QFLZ0ceWaLf/GZiTAO3Qw4p\njuhxsQxCGxoa5ThGfr+JazCe7AFwb1mM4mJ1mF6fqQPmZGdJ0oEsW2wIfnKQMUro/iThU9WOPdYv\ntzswbzjac4RXBJMsnXpq5qxF6km6NhvNtPVaAHJD58GyMYEpcIccUiyXy5JtW3FXXOtOeIaknQn7\ntyMp5022WAaht9322v71SXZcg/F0D4BLS01WZP1zcS1SpmbZMh3BTw6yLCV0sbxtt0/fe/rpZm3b\nZmnnTlvhAVEk5whfxC9Jf/2rK2c2xEwGppEBSIbOg2Urhcn49gFox41Ig205/fQsmefbSWAQ6lXw\nee3b19prJidRg/HgADhTNmTNVL0Fidl4wTFTs2yZjuAnRyVjsXxhofTFFx2zN2VlkZ/D45HKytqD\nnw0b4s9KAQCiE361uG/fprApY8nXcQDasUCBZUkPPzw26mP6/U5Ytbf0CAw625/PSy9d2WsmJ1EZ\nm+AA2OWyIp6ml496y5JwwbFr2RgU9oaRZ4o4TqA4QEx7rMUgWYvlw7M3brdRcXF05+jTR4o2Y5TL\nWE8EINXCrxYfckhqN1Buz1JInQsUeDxWTAP3wEaMXQdwsWS1ioraA6neqqeFC1TIizyTU1c3TJZl\nJDk9TlmKpJpbMADMpcX8iUaWJDa5GBQS/KSA47QXINi61UpZAJRojY2WHnqoMJRViqWCXfD+tm2S\nXsI50sAinQFIutYTZatZswp17bWFTJcEkiA4gI6kZHKswrMUAe2BRvheMdFoDzoOnEoXy9Ilx2kO\nHae36mnhAhXyIs/klJWVyrYtuVxWj4PxSKq5BQPAXFrMDyQLwU8KtLW1Fwfw+ZK/2akxyckwNTcH\nBuvBrFI862EtK/qsVLRBSqSBRW/3C2btGHCnl9crPfFEgf74xz4JqWIIoKP2DErvJZMzScego6u9\nfqITbQYnqL1Cnkno4vNI1gbF2uZ45GLpaOQHgp8UKChQh6liydrs1OsN7MfjOFZWZ5i6k44sSbA8\nd6LKhiN24dUCE1HFEEBH7RubOpJatGLFlqwY0JaXF+2fPiZ1t9dPNKLN4ASFV8iLZFpVpKWuI1kb\nFGub45GLpaMTIdVBYS6uyUk2Rg8pYNuKa6pYpAKDwdRlmPIBA+7MEb7ejPViQOK1D6AtScVqbLw1\nKwa0tm2Hpo8FxFf8IDyDk8zqaZGWuo6kmlswAExlyeNcLB2dCKkOCnNxTU6ysdggRYJTxZIpMBgM\ndFzJzDDlk+CA2+ezkjLgdpzAtEivN/HFKXKNxyN973ttamiQHnywldcLGSvTrsJG2p7y8iJt2iTl\n+4A2PIMTqJ6WnEpykQYPwXVSmzdb3bYnGAAGF/OnQnX1Lq1ZEwyW2WAziKAw8xH85BCPJ7AfjzEm\noRmmmhqfFiwoUHNz7/eNhGUFdoTOBh5PIGvX1pb4Ag3BKXU+n6WRI4uTXgAiF9x3H1e3kPky7Sps\npO2xbVuWFehDkjmgDRRWCF+jkjjB6m4mAfGKZSU3kI02eCguzqzAuq5umJYta5LPZ2v06MXsL7Qf\nQWHmYw5PDrKs+IoRdFZb2xpzFZ6u2HbsVX2SobeCBskqG86UOgCZxratpE+fai+skPjgx7bj6/+C\n6zUCjG68cWtC2tWVSEtdB5WWmrgC6+Bz+/LL5oSsRempdHQ+r0PpbT8hpB+jLeS1dBY0YA0LgEwU\nPn0qGXuhdK5MJiU/yxKp9rLSloyxkrpeI9JS14kSfG779nVci5KMBfr5vA6F/YQyH8EPclYkJaoT\nmX2prS3UrbdGfv/glLrKSocpbwDyRufKZFIgWxPLYDl8A1DHiX9GQSrXa9TWFqa0Kmt3z62nBfq5\nVs46nzNSaEfwg5SqqfGlZL1PpBmdSLIvkX5YLl3q1qJF0bUzWVPqACBTtZemduR2t4XW6cQifANQ\nY6wDAiDHMXKiiDAiKSudKEuXuhOyNilS3T23ngK+XCtnnc8ZKbQj+MlBydrkNBFqa1s1ZUpbTFde\notlsNNKMTiTZFz4sASBxwktTH3aYO641Op0H7oH+z4R9b6mhYW/Ex2svK21kWSZr1mtEcmEx+Nz6\n9Om4FqWngI/KZchFVHvLIcFNTiVLW7cGBvXJFGvqOJZAItrKaNGUqE5n9oX0OwDErnNlrUC1Oil8\nwB5YY+SP6HjB9RpffCF5PFbWrNeorW3tdRPw4HOzbbtDOeyeqrZRuQy5iOAnhxy4yWn8wU9tbWBz\nomC56+D/A7elLhvSVSZn0KDuA5pklqhOJDJKQHbJtQsW4Z/t2Sg4cG9pKd5f6TQw9S28XHdgjVF0\nMq0qaTIFgyJJB+wRRDlr5CKCnxySjE1Og1eSVq9uCv0/HQP2WDYbTXVGhw1LgdyXaxcswp9Pb5mD\nTBS+AWhw2ldzcyAAcrtbVF39qYYMGaBklNXOBz0FRkC2Ys1PCjU2WmpsTN4HcHCTU9tO7CanQdGs\nuUm0TK+MFpyWt2WLUl4yGwDikSvrGktLjfr18+/P/ASqyr3wQn898EBqAp8RI9o/+HOhMhpyS21t\nYcNHynQAABKcSURBVGg2T76La3j8/vvva8qUKQf8fPny5Ro/frwmTZqkRdGWv8phzc2Bf8mW6E1O\npfTuhxOUrExOIkpfdldggQ8bANkiVdU4k6l9A1VbLS0lKa1OtnnzU/v/Z2nJkqm66aaXe32MZSnr\nX3McKBNLai9d6s7K7G4yxPwq/PrXv9aSJUtUUlLS4ec+n08PPPCAnn76aRUWFuqSSy7ROeeco/Ly\n8rgbmy2Cb/h0vskS/UcX7ZqbbJKIK57dTctL51RBAIhGcNF8c3Nqp4glcvDfeQPVVFYne/vtjud+\n663er0Lm09qifEKfn9liHp1XV1dr7ty5mj59eoefb9iwQdXV1fLsvzw/aNAgrVq1SiNHjoyvpVkk\n+KZPZ/CT6D+8WNbc5JPgtDy/38rIaXkAkKkSOfh3u539xX7SUZ1su4LrbgNfd6Tw3OHFiAr3f2UA\nDnQl5tH5iBEjtGnTpgN+7vV61a9fv9D3JSUl2rNnT6ynQYbIlupp4VKdcrZtye2m2AEARCORn9Xl\n5UXavDlQ7a1v3+aUVic7/XSPXnpJCgQ+f9DppyeuM4jkNQoGO4MGlXT4PpkybWoXEImEpyY8Ho+8\nYQtCmpqa1L9/ZGlne/9ClYqKfr3cMzsE193YttXp51bSn2Oijh98DhUV/eR2Bwb3Rx+dnt9PeFsi\nMXdu8H99ktKezrpqX7RtBpDZcuVveeLEwNfOzye830rWc504UfrVr6SmpsB55s4NfkbH9lkd3ubA\nBqqB7ysrPTruuOgDkFhfg9/+9lIdckig2MLFF1t67LFLVF7e/eOjOU80/VlP/U53t8XymGjbhfRi\nPNIu7uAnUEu/3THHHKONGzdq9+7dKioq0qpVq3TFFVdEdCzHCUylqq9virdZGWHMmEItWFDQYbdp\nSWppMfr006YkZQgCB62vT0xFAscp2X+8pg7/T4dozp+OtL/jlMi2bdXXt2c6fb4StbUpib9vILHo\nGHsW/vedzYIz1uvrO/488DkbKBWdrM/66dOlhQtLZIylMWPaVF8f3+d0eJsDrP0/j+05OE6Jiotj\naZsrdLFzzpyx8vt7fr8k67Xuqa/sqp+K5DHd3YbskUu/x3j7qbhrgllW4A/9+eef16JFi+R2uzVj\nxgxNmzZNl1xyiSZMmKDKysp4T5OVamtbD5jL7PcrrRXT8kU6qprU1Pg0YUL795lQIQ8AMlVpqUnY\nBari4sROwUpk24B0a2ho1PbtzdqyZS9l2BVn5ueII47Qk08+KUkaO3Zs6OdDhw7V0KFD42pYrgiv\n/BZIkuVmxTQEgt2Kij6hK6m5XCEPHc2aVaiGBunBB1vJ8AFpEB6sLFhQkJJtJRIh0UEb0JXbbntN\nLS1TFCzDLi3I601rKfidZOGV3wJJskAlmGRVTMvl/QKyrYOgQl5+8HqlJ54okM9n6f/+z5U1BUGA\nXJSuct2xIsOEVAiUfE9PCfhMRPCTYi6XdPDBTtIGSLm6X0A2lu7Mxgp5iB4ZPgDZIrCRLcUJ8k11\n9S6tWZOuEvCZh+AnDQoLKYccrWzdLNS2+X3nOjJ8ALJF5+nZQT3NrMi2WRc4UF3dMC1b1iSfz9bo\n0YtTWgI+ExH8oEeOI7W1KbRYv7HRUm1tYVqCkHSeG+iOxyN973ttrPlBTqip8enNN1064wx/upuC\nFOqpX6XPzX5lZaU65JDAuoh8XusTRPCDbgWrlfl8lkaOLJbjSM3NgSxMOj4Mozl3eNCWzsEoV8zy\nw333MThAbmCgmxqB6luWfD5bV131tOrqhqmsrDTdzQLyAsEPutV5LUNlZXZM5+kctKVzvQ0DCQA4\nUCZfGEpF26i+BaQPwQ+61XktQ7ZUz2EBOgBktky+MJSKtlF9C0ifuDc5Re4KViurrAxUpzv/fF9W\nlNIOBm2SWIAOAMg41dW7FNj6QqL6FpBaBD8pEigvme5WRC+8WlltbWtWlNLuHLSxAB0AcldNje+A\nqWqZvnloXd0w9e3bpIKCZp1//oK8r74FpBLT3lIkVRuvZfKHfSpRYhoA8kNX09QyffNQqm8h1Rgf\ntiP4yTGZ/GEfD8eRjAl8BQAAQORydXwYC6a9IeMFq7c5jqWtW63QnkMAAKBdoIR2s7Zs2aurrnpa\nO3c2prtJQMYh+EHGC6/e5vNZWreOty0AAJ0FSmiXqK2tWEuWTNX06a+lu0lAxmHaGzJeeMltl8vo\nyCOZ+wYA6Fo61zZ8/evp7Z8ooQ30juAHGc/jkSorjTZvlvx+SxdemN6NSwEAmSudaxuWLWtO27ml\nQAntNWuMAgEQJbSBrjB/CFnB55OCV7OCG5f2pKvSpwAA5DJKaAO9I/ODrFBQIAU2hLMi2riUqiYA\ngHxDCW2gdwQ/yAq2Lblc0sEHs3EpAAAAYsO0N0TFcaTWVqWt3DQblwIAACBWBD8plO0bdQb329m2\nzdbIkcXstwMAAICswrS3FOm4UWfg+2zIYIQXDQjfbydYdGDQoCyN5AAAAJB3CH5SpKuNOrMhcAgv\nHBC+304kRQcAAACATMK0txQJBg6S5HabrAwcPB6pqsqosjL1RQdqanwqLk7d+QAAAJB7CH5SJBg4\n2LZRVZXJiilvXbHt9BQdqK1tVWmpSe1JAQAAkFMIflLItiXLCnwFAAAAkFqs+UHWCC++AABANqNP\nA9KD4AdZI7z4AgAA2Yw+DUgPJmABAAAAyAsZlfkhBQwAAAAgWTIq+CEFDAAAACBZMir4yQfFxWS4\nAABAcjDGAHpG8JNipaWGDFcWq60t3P+V3yEAIPPQPwE9o+ABEIWlS91aupRrBgAAANmI4AcAAABA\nXojpErYxRrW1tVq3bp369Omj++67TwMGDAjdPn/+fC1evFjl5eWSpHvuuUdHHXVUQhoMAAAAALGI\nKfh59dVXtW/fPj355JN6//33df/99+vRRx8N3b527VrV1dXpxBNPTFhDAQAAACAeMQU/q1ev1uDB\ngyVJ3/jGN/TXv/61w+1r167VvHnzVF9fr6FDh+rqq6+Ov6XICPleRcZxpLY2yeuVPJ50twYAAADR\niCn48Xq96tevX/tB3G45jiPbDiwhGjNmjC677DJ5PB5df/31WrFihc4+++zEtBhplc9VZLxeaetW\nSz6fpZEji/Xyy80EQAAAAFkkpuDH4/Goqakp9H144CNJl19+uTz7R4Vnn322Pvzww4iCn4qKfr3e\nJ5sFX6Jcf5656pNPJN/+xNf69S5t29ZPRx994P34/QK5i79v5ALex8hnMQU/p5xyil577TWdd955\nWrNmjY477rjQbV6vVzU1NXrxxRdVVFSkt99+W+PHj4/ouPX1e2JpTtZwnBJJUn19Uy/3RCaqrJTc\nbo98PkvHHutXZWWz6us73qeiol/Ov4+R2xgU9Yy/b2Q7+ilku3j7qZiCnxEjRmjlypWaNGmSJOn+\n++/X888/r5aWFk2YMEG33HKLpkyZosLCQp1++ukaMmRIXI0EMoHHI1VVGbW1Gaa8AQAAZCHLGGPS\n3YigXL8SMWhQIPOzejWZn2zV2++QK2rIdmR+esbfN7Id/RSyXbz9FJucAgAAAMgLBD8AAAAA8gLB\nDwAAAIC8EFPBAyBf5fsmrwAAANmM4AeIQj5v8goAAJDtmPYGAAAAIC8Q/AAAAADICwQ/AAAAAPIC\na35SiMXyAAAAQPoQ/KQQi+UBAACA9GHaGwAAAIC8QPADAAAAIC8Q/AAAAADICwQ/AAAAAPICwQ8A\nAACAvEDwAwAAACAvEPwAAAAAyAsEPwAAAADyAsEPAAAAgLxA8AMAAAAgLxD8AAAAAMgLBD8AAAAA\n8gLBDwAAAIC8QPADAAAAIC8Q/AAAAADICwQ/AAAAAPICwQ8AAACAvEDwAwAAACAvEPwAAAAAyAsE\nPwAAAADyAsEPAAAAgLxA8AMAAAAgLxD8AAAAAMgLBD8AAAAA8kJMwY8xRj/60Y80adIkTZ06Vf/4\nxz863L58+XKNHz9ekyZN0qJFixLSUAAAAACIR0zBz6uvvqp9+/bpySef1M0336z7778/dJvP59MD\nDzyg+fPna8GCBVq4cKEaGhoS1mAAAAAAiEVMwc/q1as1ePBgSdI3vvEN/fWvfw3dtmHDBlVXV8vj\n8aigoECDBg3SqlWrEtNaAAAAAP+/vbsJjepswwB8R1OtZFKqNi51UUgXUirBhZtKXWRbhKpEIYLi\nsuKqlCJU8Ido0WWzakFxUZW6ENRF6SaLoCDxByq6E1tcxXahmaqpZL5FMd8XPjRxojOZea9rd+ac\nmTwMz5kn95yXOdSps54nTUxMpLu7+78v0tmZqampLFq06P/2dXV15fHjx3N63Z6e7tkPggVOH0P7\ncn7TDvQxJasr/FQqlVSr1entF8Hnxb6JiYnpfdVqNe+9996cXnd8fG4hCRaqnp5ufUxL80/Rqzm/\naXXmFK1uvnOqrmVvfX19GRkZSZLcvHkzvb290/s+/PDD3L9/P48ePcrk5GSuXbuWdevWzatIAACA\n+arryk9/f39GR0czMDCQJBkaGsrFixfz5MmTbN26Nd988012796dWq2WrVu3ZtWqVW+0aAAAgNfV\nUavVas0u4gWXYWl1lhPQ6ix7ezXnN63OnKLVNWXZGwAAQKsRfgAAgCIIPwAAQBGEHwAAoAjCDwAA\nUAThBwAAKILwAwAAFEH4AQAAiiD8AAAARRB+AACAIgg/AABAEYQfAACgCMIPAABQBOEHAAAogvAD\nAAAUQfgBAACKIPwAAABFEH4AAIAiCD8AAEARhB8AAKAIwg8AAFAE4QcAACiC8AMAABRB+AEAAIog\n/AAAAEUQfgAAgCIIPwAAQBGEHwAAoAjCDwAAUAThBwAAKILwAwAAFEH4AQAAiiD8AAAARRB+AACA\nInTW86Rnz57lq6++yp9//plKpZKjR49m+fLlM445cuRIrl+/nq6uriTJ8PBwKpXK/CsGAACoQ13h\n56effkpvb2++/PLLXL58OcPDw9m/f/+MY27fvp0ff/wx77///hspFAAAYD7qWvY2NjaWjRs3Jkk2\nbtyYK1euzNhfq9Vy//79fPvtt9m+fXvOnz8//0oBAADmYdYrPz///HNOnTo147EPPvhgeglbV1dX\nJiYmZuz/+++/Mzg4mF27duX58+fZuXNnPv744/T29r7yb/X0dL9u/bDg6GNoX85v2oE+pmSzhp8t\nW7Zky5YtMx7bu3dvqtVqkqRaraa7e+ZJtGzZsgwODmbp0qVZunRpNmzYkLt3784afgAAAN6Wupa9\n9fX1ZWRkJEkyMjKS9evXz9h/79697NixI7VaLf/880/Gxsaydu3a+VcLAABQp45arVZ73Sc9ffo0\nX3/9dcbHx7NkyZKcOHEiK1euzMmTJ7NmzZps2rQpJ0+ezKVLl/LOO+9k8+bN2bZt29uoHwAAYE7q\nCj8AAACtxk1OAQCAIgg/AABAEYQfAACgCMIPAABQBOEHAAAowqw3OW2GGzdu5OzZs+no6Mj+/ftT\nqVSaXRLU5erVq7l48WIOHz7c7FKgLleuXMnly5fz9OnT7NmzJx999FGzS1owzCrahVlFK3vdObUg\nr/ycO3cuBw8ezBdffJFLly41uxyoy++//547d+5kcnKy2aVA3Z49e5ZDhw5l9+7dGR0dbXY5C4pZ\nRTswq2h1rzunGh5+bt26lcHBwSRJrVbLgQMHMjAwkJ07d+aPP/5IkkxNTWXJkiXp6enJ+Ph4o0uE\nWc2lj1evXp1du3Y1s0x4pbn08WeffZYnT57k9OnT2bx5czPLbSizinZgVtHq3sacauiytx9++CEX\nLlxIV1dXkuTXX3/N5ORkzpw5k1u3bmVoaCjDw8N59913Mzk5mfHx8fT09DSyRJjVXPv4BfcRZiGa\nax//9ddfOX78ePbt25cVK1Y0uerGMKtoB2YVre5tzamGXvlZs2ZNvv/+++ntsbGxfPrpp0mSTz75\nJLdv306SbNu2LQcOHMjZs2fz+eefN7JEmNVsffzbb7/NOL6jo6Oh9cFczPXz+NixY3n48GFOnDiR\nX375pSm1NppZRTswq2h1b2tONfTKT39/fx48eDC9PTExke7u7untxYsXZ2pqKmvXrs3Q0FAjS4M5\nm62POzs7MzU1lUWL/v1u4bvvvmt4jTCbuX4eHzt2rBnlNZVZRTswq2h1b2tONfUHDyqVSqrV6vT2\n/56E0Cr0Me1AH7+c94Z2oI9pdW+qh5va9X19fRkZGUmS3Lx5M729vc0sB+qij2kH+vjlvDe0A31M\nq3tTPdzU+/z09/dndHQ0AwMDSWL5AC1JH9MO9PHLeW9oB/qYVvemerij5uc9AACAAljsCQAAFEH4\nAQAAiiD8AAAARRB+AACAIgg/AABAEYQfAACgCMIPAABQBOEHAAAown8A46i6f9IdMw8AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a1c5cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Obtain unpooled estimates\n",
    "unpooled = pd.DataFrame({'n': n_county,\n",
    "                         'm': unpooled_model.unpooled_estimates,\n",
    "                         'sd': unpooled_model.unpooled_se})\n",
    "unpooled['se'] = unpooled.sd/np.sqrt(unpooled.n)\n",
    "\n",
    "# Construct axes for results\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14,6),\n",
    "                         sharex=True, sharey=True)\n",
    "jitter = np.random.normal(scale=0.1, size=counties)  # avoid overplotting counties\n",
    "\n",
    "# Plot unpooled estimates\n",
    "axes[0].plot(unpooled.n + jitter, unpooled.m, 'b.')  # means\n",
    "for j, row in zip(jitter, unpooled.iterrows()):\n",
    "    name, dat = row\n",
    "    axes[0].plot([dat.n + j, dat.n + j], [dat.m - dat.se, dat.m + dat.se], 'b-')\n",
    "\n",
    "# Plot partially-pooled estimates\n",
    "axes[1].scatter(n_county.values + jitter, means)\n",
    "for j, n, m, s in zip(jitter, n_county.values, means, sd):\n",
    "    axes[1].plot([n + j, n + j], [m - s, m + s], 'b-')\n",
    "\n",
    "# Add line for underlying mean\n",
    "for ax in axes:\n",
    "    ax.hlines(sample_trace.mean(), 0.9, 100, linestyles='--')  # underlying mean from partial model\n",
    "\n",
    "# Set axis limits/scale (shared x/y - need only to set one axis)\n",
    "axes[0].set_xscale('log')\n",
    "axes[0].set_xlim(1, 100)\n",
    "axes[0].set_ylim(-0.5, 3.5)\n",
    "\n",
    "# Set axis titles\n",
    "axes[0].set_title(\"Unpooled model estimates\")\n",
    "axes[1].set_title(\"Partially pooled model estimates\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By inspection, there is quite a difference between unpooled and partially-pooled estimates of prevailing county-level radon level, especially as smaller sample sizes. The unpooled estimates at smaller sample sizes are both more extreme, and more imprecise."
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}