05-unpooled_model.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",
    "# 5. An Unpooled Model"
   ]
  },
  {
   "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",
    "\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building the model in `Stan`\n",
    "\n",
    "We construct a model with no pooling, where we treat each county as having a prevailing radon level that is independent of all other counties: \n",
    "\n",
    "$$y_i = \\alpha_{j[i]} + \\beta x_i + \\epsilon_i$$\n",
    "\n",
    "* $y_i$: measured log(radon) in household $i$\n",
    "* $\\alpha_{j[i]}$: prevailing radon level in county $j[i]$, which is the county $j$ to which household $i$ belonds\n",
    "* $\\beta$: effect on measured log(radon) in moving from basement to ground floor measurement\n",
    "* $\\epsilon_i$: error in the model prediction for household %i%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specifying the unpooled model in Stan\n",
    "\n",
    "For the unpooled model, we have the parameter $\\alpha_{j[i]}$, representing a list of (independent) mean values, one for each county. Otherwise the model is the same as for the pooled example, with shared parameters for the effect of which floor is being measured, and the standard deviation of the error.\n",
    "\n",
    "We construct the `data`, `parameters` and `model` blocks in a similar way to before. We define the number of samples (`N`, `int`), and two vectors of log-radon measurements (`y`, length `N`) and floor measurement covariates (`x`, length `N`). The main difference to before is that we define a list of counties (these are the indices 1..85 defined above, rather than county names), one for each sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "unpooled_data = \"\"\"\n",
    "data {\n",
    "  int<lower=0> N;\n",
    "  int<lower=1, upper=85> county[N];\n",
    "  vector[N] x;\n",
    "  vector[N] y;\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define three parameters: $\\alpha_{j[i]}$ - one radon level per county (`a` - as a vector of length 85, one value per county); change in radon level by floor, $\\beta$ (`beta`, a real value), and the Normal distribution scale parameter `sigma`, as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "unpooled_parameters = \"\"\"\n",
    "parameters {\n",
    "  vector[85] a;\n",
    "  real beta;\n",
    "  real<lower=0, upper=100> sigma;\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also define a *transformed parameter*, for convenience. This defines a new variable $\\hat{y}$ (`y_hat`, a `vector` with one value per sample) which is our estimate/prediction of log(radon) value per household. This could equally well be done in the model block - we don't *need* to generate a transformed parameter, but for more complex models this is a useful technique to improve readability and maintainability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "unpooled_transformed_parameters = \"\"\"\n",
    "transformed parameters {\n",
    "  vector[N] y_hat;\n",
    "  \n",
    "  for (i in 1:N)\n",
    "    y_hat[i] <- beta * x[i] + a[county[i]];\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this transformed parameter, the model form is now $y \\sim N(\\hat{y}, \\sigma_e)$, making explicit that we are fitting parameters that result in the model *predicting* a household radon measurement, and we are estimating the error of this prediction against the observed values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "unpooled_model = \"\"\"\n",
    "model {\n",
    "  y ~ normal(y_hat, sigma);\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the unpooled model in `Stan`\n",
    "\n",
    "We again map Python variables to those used in the `Stan` model, then pass the data, parameters (transformed and untransformed) and the model to `Stan`. We again specify 1000 iterations of 2 chains.\n",
    "\n",
    "** Note that we have to offset our Python indices for counties by 1, as Python counts from zero, but `Stan` counts from 1. **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map data\n",
    "unpooled_data_dict = {'N': len(clean_data.log_radon),\n",
    "                      'county': clean_data.county + 1,  # Stan counts start from 1\n",
    "                      'x': clean_data.floor_measure,\n",
    "                      'y': clean_data.log_radon}\n",
    "\n",
    "# Fit model\n",
    "unpooled_fit = pystan.stan(model_code=unpooled_data + unpooled_parameters +\n",
    "                                       unpooled_transformed_parameters + unpooled_model,\n",
    "                            data=unpooled_data_dict,\n",
    "                            iter=1000,\n",
    "                            chains=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting the fit\n",
    "\n",
    "We can extract the sample from the fit for visualisation and summarisation. This time we do not use the `permuted=True` option. This returns a `StanFit4Model` object, from which we can extract the fitted estimates for a parameter using indexing, like a dictionary, e.g. `unpooled_fit['beta']`, and this will return a `numpy` `ndarray` of values. For $\\alpha$ (`a`) we get a 1000x85 array, for $\\beta$ (`beta`) we get a 1000x1 array. Mean and standard deviation (and other summary statistics) can be calculated from these.\n",
    "\n",
    "When extracting vectors of $\\alpha_{j[i]}$ (radon levels per county) and the associated standard errors, we use a `pd.Series` object, for compatibility with `pandas`. This allows us to specify an index, which is the list of county names in `clean.data.mn_counties`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AITKIN      0.897083\n",
       "ANOKA       0.931966\n",
       "BECKER      1.541253\n",
       "BELTRAMI    1.603539\n",
       "BENTON      1.459515\n",
       "dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Extract fit of radon by county\n",
    "unpooled_estimates = pd.Series(unpooled_fit['a'].mean(0), index=clean_data.mn_counties)\n",
    "unpooled_se = pd.Series(unpooled_fit['a'].std(0), index=clean_data.mn_counties)\n",
    "\n",
    "# Inspect estimates\n",
    "unpooled_estimates.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To inspect the variation in predicted radon levels at county resolution, we can plot the mean of each estimate with its associated standard error. To structure this visually, we'll reorder the counties such that we plot counties from lowest to highest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAFtCAYAAADBM4kgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYlNXiB/DvDAMqYKIIZaaYu+bV61LXfTe3zD00oZuZ\npak/86aoSTqKqKDlggsuj92rlrnhdu1eNZfqkomVu2LdMr0qKSq4gAIzc35/EJPIO8zAbO+Z+X6e\npyfnnZn3Pecd4DvnvOc9RyOEECAiIiIpaN1dACIiIrIdg5uIiEgiDG4iIiKJMLiJiIgkwuAmIiKS\nCIObiIhIIgxuksJrr72GVatWFdm+du1avPPOOyXa15IlS7Bz585iX3P//n389a9/NT/u378/7t+/\nX6LjWBIdHY0jR444ZF+O9sEHH+DcuXMuOVafPn1w7NgxlxxLyZdffoklS5YAAA4ePIjY2FiH7duV\n55G8D4ObpDBs2DAkJSUV2b5lyxZERkaWaF//93//h759+xb7mszMTJw+fdr8ePv27QgMDCzRcSyZ\nPXs2WrVq5ZB9OVpycjK8ZWqH06dP4+7duwCAzp07Y9q0aQ7btzedR3I9nbsLQGSLrl27Ys6cOfj+\n++/RvHlzAEBKSgoAoFWrVhBCYM6cOTh16hSysrIghMDs2bPRtGlTTJ06FZmZmbhy5Qo6duyImzdv\nom7duhg+fDi2bt2KzZs3w2AwIDMzE2+99RaGDBmC999/Hw8fPkT//v2xbds2NGzYEN9++y2CgoKw\nbNkyfP7559DpdKhRowamT5+O4OBgREZGomnTpvjhhx9w7do1tGjRAvHx8UXqEhkZicjISDz33HN4\n/fXX0aFDB5w8eRJ3797Fu+++i549exZ5z6FDh7B48WIIIVCuXDno9XrUr18fX3zxBZYtWwaTyYTA\nwEBMnjwZjRs3xtKlS5GZmYno6GgAKPT48XI+//zzmDdvHhYtWoQbN25g4sSJmDVrFt5++2189dVX\n5i8s3bt3x5IlS1CvXj1zubZv346tW7fiwYMHKF++PBITE6HX63Hp0iVkZmYiICAAH374IWrUqIGf\nf/7ZfF5r1KiB7Oxs836Kq8fVq1dx48YNXLt2DcHBwVi4cCFCQkKKnKPExETs27cPQghUrVoVM2bM\nQEhICPbt24fExERotVr4+Phg0qRJ8PPzw2effWY+XlhYGPbu3YvExERERkaiUaNG+Pbbb3H79m1E\nRkbi1q1bSElJwcOHD7Fo0SLUqVMHJ06cwIIFC5CXl4f09HS0bt0as2fPxsKFC83nMS4uDjVr1kRs\nbCx+/PFHGAwGtGrVClFRUdBqtViyZAkOHDgAX19fBAUFYd68eahcubIdvynkFQSRJBISEsSUKVPM\nj9977z2xbt06IYQQx48fF+PHjzc/t3LlSjFq1CghhBBTpkwRw4cPNz83ZcoUsXbtWpGVlSXCw8NF\nZmamEEKIEydOiKZNmwohhLhy5Yr530IIUb9+fZGRkSG2bt0qhgwZIh4+fGgu04gRI4QQQkRERIh3\n331XCCHE/fv3Rbt27cTRo0eL1CMiIkLs3btXXLlyRdSrV08cPnxYCCHE3r17RadOnYq8/ubNm6JF\nixYiNTVVCCHEvn37xMiRI8XPP/8s2rRpI65cuSKEEOLIkSOiTZs24v79+yIhIUHExMQUOncFj4sr\nZ6dOncTZs2eFEEKMGTNGfPrpp0IIIb755hsRHh5epGxJSUnihRdeEFlZWUIIIf7973+L2bNnm5+f\nPn26+bj9+vUT27ZtM5/rhg0bipSUFKv16Natm3n/o0aNEgkJCUXKsX37djFhwgRhNBqFEEJs2rRJ\njBw5UgghRNeuXcXJkyeFEEIkJyeLZcuWFTknSUlJ4u233zafn3HjxgkhhDh58mShz2jOnDnigw8+\nEEII8be//U2kpKQIIYTIysoSLVu2NJ+7R8/j1KlTxYYNG4QQQhiNRjFp0iSxZs0akZaWJpo3by5y\nc3OFEEJ8/PHH4osvvihSN6LHscVN0ggPD8dLL72E7Oxs5ObmIjk5GXq9HgDw5z//GePHj8fGjRtx\n+fJlpKSkFOrabtasWZH9+fv7IzExEYcOHcKlS5dw/vx5PHjwoNgyfP311xgwYADKlCkDIP/ae2Ji\nIgwGAwCgU6dOAICAgACEhYXhzp07xe7P19cXHTp0AAA0bNhQ8fU//PAD6tata27pduvWDd26dcOn\nn36KVq1aoWrVqgCAli1bonLlyjh79myxx7RWTvF7F++rr76KBQsWYOjQodi8eTOGDh2quK969erB\n398fQH6rvFq1atiwYQMuXbqElJQUNG3aFJmZmbhw4YL5EkWTJk3M9fn222+LrccLL7xg3n/Dhg2R\nmZlZpAyHDx/G6dOnMWDAAACAyWRCTk4OAKB3795455130LFjR7Ru3Rpvvvmm1fPz4osvAgCqVasG\njUaDtm3bAgCqV69u7umZN28evvzyS6xcuRK//PILHj58WKgXoeA8FpRty5YtAICcnBxotVo89dRT\naNCgAfr374927dqhffv2qr2EQurC4CZphISEoHXr1tizZw+ys7PRvXt3czgfPnwYc+bMwRtvvIGu\nXbuiZs2a2L17t/m9AQEBRfZ3/fp1hIeHIzw8HC1atED37t3x5ZdfFlsGk8lU6LHRaITRaDT/kS5b\ntmyh54WV65y+vr7mf2s0GsXX+/j4QKPRFNr2008/Kb7WaDSav0Q8+nxeXl6h19lSztatW+PBgwc4\ncuQIvvvuO8TFxSnWoSBUAeDTTz/Fli1bEBERgT59+qBChQq4evUqNBpNkfr5+PhYPPaj9Xi0rI+f\nhwImkwkjR47EkCFDzPXNyMgAALz77rsYNGgQkpOTsX37dqxevRrbt29X3E8BPz+/Qo8LyvqoV199\nFQ0aNED79u3Rs2dPnDx50mJdFi9ejJo1awIA7t27Z67H+vXrcebMGXzzzTeYO3cu/vKXvzj0Wjt5\nJg5OI6kMHToUu3btws6dOzFs2DDz9m+++QadO3fGkCFD8Kc//QkHDhwoErKPO336NCpVqoTRo0ej\nTZs2OHToEID8INHpdIXeX/AHuV27dkhKSjK3zNevX4/nn3++UACXxON/6JX+8Ddp0gQ///wzfv75\nZwD514MnTpyIli1bIjk5GVeuXAEAHDlyBNevX0fjxo1RqVIlc4s1Ozsb//nPf2wqj06nMwcmkH++\no6Oj0adPnyJhpiQ5ORkDBgzAwIEDUaNGDRw6dAgmkwkVKlTAc889Z251nj9/HufPnweAYuthq7Zt\n22LLli3mkf8JCQmYOnUqjEYjOnfujOzsbISHh2PGjBn45ZdfkJeXBx8fnyJfaJQofSZ3797FuXPn\nMGnSJHTt2hW//fYbLl++DKPRCKDweWzXrh3+/ve/QwiB3NxcjBkzBp988glSU1Px0ksvoVatWnjr\nrbfw+uuv48KFCzbXmbwXW9wklRdeeAGZmZmoWLEi6tSpY94+ZMgQTJw4Ef3798cTTzyBLl26YO3a\ntcXuq127dti2bRu6d++OypUro0uXLqhcuTIuXbqE6tWro0GDBujVqxc+/fRTcwtp0KBB+O233zB4\n8GAIIVC9enXMnz8fQNHWoKXW4aPbbXlPcHAwFixYgKioKPNgqoULF6JmzZqYMWMGxo4dC6PRiHLl\nyiExMRGBgYF4+eWX8fXXX6N79+4IDQ1F06ZNbTpmly5dMGHCBMyePRutW7dG//79ER8fb27JWvPG\nG29g+vTp2LFjB4KCgtC1a1dzL8aHH36IqVOnYuPGjQgLC0OtWrUAALVq1bJYD1sNHjwYN27cQHh4\nOLRaLapUqYK5c+fCx8cH06ZNw3vvvQdfX19otVrMnTsXvr6+aNWqFcaNGwdfX18899xzNp2fAk88\n8QTeeust9OvXD08++SRq166N9u3b4/Lly2jZsmWh8xgdHY3Zs2fj5ZdfhsFgQJs2bfDmm2/Cx8cH\nPXv2xIABA+Dv749y5cqZBxMSFUcjrPXlEZHX2rNnD3bu3Kl4Dz0RuYfbuspv3bqFjh074uLFi+4q\nAhEVIzIyEsuXL8fkyZPdXRQieoRbusoNBgNmzJhRZIAMEanH+vXr3V0EIlLglhZ3XFwchg4ditDQ\nUHccnoiISFouD+6kpCQEBwejTZs2nBKQiIiohFw+OC0iIsI8SjM1NRXPPvssVqxYgeDgYIvvMRiM\n0OmK3kdJRETkbdw6qjwyMhKzZs3Cs88+W+zr0tPvFdkWElJecbuMWBf18qT6sC7q5Un1YV0cd2xL\n3DoBi6X7XImIiEiZWydgWbdunTsPT0REJB1OeUpERCQRBjcREZFEGNxEREQSYXATERFJhMFNREQk\nEQY3ERGRRBjcREREEmFwExERSYTBTUREJBEGNxERkUQY3ERERBJhcBMREUmEwU1ERCQRBjcREZFE\nGNxEREQSYXATERFJhMFNREQkEQY3ERGRRBjcREREEmFwExERSYTBTUREJBEGNxERkUQY3ERERBJh\ncBMREUmEwU1ERCQRBjcREZFEGNxEREQSYXATERFJhMFNREQkEQY3ERGRRBjcRERELqDXl4FeX8bu\n/TC4iYiIXGD3bh1279bZvR8GNxERkUQY3ERERBJhcBMREUmEwU1ERCQRBjcREZFEGNxEREQSYXAT\nERFJhMFNREQkEQY3ERGRRBjcREREEmFwExERSYTBTUREJBEGNxERkUQY3ERERBJhcBMREdnBUets\n24rBTUREZAdHrbNtKwY3ERGRgzmzFc7gJiIicjBntsIZ3ERERBJhcBMREUmEwU1ERCQR1w2De4TJ\nZEJ0dDQuXrwIrVaLmTNnonbt2u4oChERkVTc0uI+ePAgNBoNNm7ciPHjx+Ojjz5yRzGIiIik45YW\nd9euXdG5c2cAwNWrV1GhQgV3FIOIiEg6bgluANBqtZg6dSr279+PJUuWuKsYREREUtEIIYQ7C3Dr\n1i0MHjwYn3/+OcqWLav4GoPBCJ3Ox8UlIyIibzVpUv7/58+3/toaNfL//+uvJd9WGm5pce/YsQPX\nr1/H22+/jTJlykCr1UKrtXy5PSMju8i2kJDySE+/58xiugzrol6eVB/WRb08qT6eUpdNmwKg1WoR\nFWW9LiZTAAAgPT2rxNssCQkpb/E5twxO69GjB86fP4+IiAi8+eabmDZtGvz8/NxRFCIiIqe7fTsT\nN29m47ffHmLkyCRkZGSWel9uaXGXLVsWixYtcsehiYiIXG7y5EN48CASgAY7d74GYD1Wr+5fqn1x\nAhYiIiInu3TpCQCa3x9pfn9cOgxuIiIiJwsLuwOgYCy4QFjY3VLvi8FNRETkZPHxnVGuXBZ8fbPR\nt+96xMd3KvW+3HYfNxERkbeoWDEIlSv7A0Cpr20XYIubiIhIIgxuIiIiiTC4iYiIJMLgJiIikgiD\nm4iISCIMbiIiIokwuImIiCTC4CYiIpIIg5uIiMiBHLkSmBIGNxERkQPlrwQWgLw8f+zc+Rqiog45\ndP8MbiIiIhvo9WWg15ex+jpHrgSmhMFNRERkg927ddi92/oSH45cCUwJg5uIiMiBHLkSmBKuDkZE\nRORAjlwJTAlb3ERERBJhcBMREUmEwU1ERCQRBjcREZFEGNxEREQSYXATEZFXs3ViFbVgcBMRkVez\ndWIVtWBwExERSYTBTUREVErOXglMCYObiIiolJy9EpgSBjcREXkkVww6c/ZKYEoY3ERE5JFcMejM\n2SuBKZFnGB0REZHKxMd3xr59WTAYtOjVa2uxK4H16WNwyDEZ3ERERKVUkpXA9PochxyTXeVEREQS\nYXATEVGx1DazmNrK42oMbiIiKpbaZhZTKo83hTmDm4iInMZVgersLxfumGjFEgY3ERE5jdpa66X9\nIuGOiVYsYXATEZFLubNbu7RfJNwx0YolDG4iInIptbXCbeGOiVYsYXATEbmRXl8Gkya5uxRkTXx8\nZ5QrlwVf32z07bu+2IlWnI3BTUTkRrt367Bli7tLQY8qGIiWlpZtHohWMNHKU0+VxerV/VGxYpDb\nyidXXwURkZcquCbsqNm3yLL8gWiRADTYufM1AOutzor2OEdNb6qEwU1EJIGCa8KOCm5nfBHQ68vA\n3x+IinLYLt3CEQPRnPkFi13lREReyNIAMXtGfDuj298dI9DVNBBNCYObiIjM3DXi21JAu2NilYKB\naH5+7h+IpoTBTUREbueuLwxKE6sUDESrUsXf7QPRlDC4iYioxDxlbnA1TaxiKwY3EZECTwkmZ5Fx\nEhUlar+erYTBTUSkwFOCiYqnpolVbMXgJiJyEVe04tlTYJnSQDQ1TaxiKwY3kQrxj6862fu5uKIV\nz54Cy9S0wpc9rAZ3bm4uVqxYgaioKNy7dw9Lly5Fbm6uK8pG5LX4x9d57L1P2RWfC7+42U+pdS3j\nQDQlVoN71qxZePDgAc6dOwcfHx9cvnwZ06ZNc0XZiIgczhXhq7aWuTd+EVBqXcs4EE2J1eA+e/Ys\n/va3v0Gn08Hf3x9xcXE4f/68K8pGRB7K04NEbT0maiuPo9naupZxIJoSq8Gt0WiQm5sLjSb/BGRk\nZJj/XRoGgwFRUVEYNmwYXnnlFRw8eLDU+yIiOXl6kJBjKAWy0jZbW9cyDkRTYvU357XXXsPw4cOR\nnp6O2NhYfPHFF3jnnXdKfcBdu3ahYsWKiI+Px507d9CvXz907ty51PsjIvJWnrJiWH4Ya2AwaDFy\nZBLi4zujYsUgxVW6ABTZptS63rSpOfbty4LBoEWvXlulbV0rsRrc/fr1Q6NGjXD06FEYjUasWLEC\n9evXL/UBe/bsiR49egAATCYTdDp+6yYiebkzPB29Ypg9CsLXaPwjfIVAkUBW2mZpGU3Lg8kKbwsL\nu4MTJ8Tv2wu3rgGUeElOtbOamuPGjUNCQgJq165t3vbXv/4V//jHP0p1wHLlygEA7t+/j/Hjx2PC\nhAml2g8RuYantOqcRU3h6SpKLWRbW8dK2ywFtFIgA6LItvj4zh7bulZiMbjHjBmD1NRU3LhxA126\ndDFvNxqNeOqpp+w6aFpaGsaOHYuIiAj06tXL6usrVvSHTudTZHtISHm7yqEmrIt6uaM+Wq1zjl2a\n/e3Zk///Zcv8HFYOR9SvtO+19dhKr7NnmzP26Yoy3rqViVu3smEwaDF27G6sWNEL06f/Bw8eDEFB\n+JYp8xmuXauIR8M3/zFs2la37r1CYVy37gOEhJTH2rV9Ua1afiD367cTK1a8DABFtlWqFITQ0Pw9\n7tjxmtvPmbNZDO64uDhkZmYiNjYW0dHRf7xBp0NwcHCpD3jz5k2MGDEC06dPR8uWLW16T0ZGdpFt\nISHlkZ5+r9TlUBPWRb3cVR+TKQAAkJ6e5bB9lrYuziiLvfu053Ox9dhKr7NnW3Gv1Wq1herj6GPb\nU8aRI3chOzu/hbx581Dk5BRtIf/4YzmEhd1G/mCw/PB9+umM3x9b3xYT0wm7d//RYo6J6fT7+fBB\ncHB+d/fSpS/BaMw/4uPb0tPvOeX8PP659O5d5vfXFe5dcdbvqyUWgzswMBCBgYFYsWIFzp07h+zs\nbAghYDQa8fXXX2PQoEGlKszKlStx9+5dLF++HMuWLYNGo8GaNWvg5+e4b/NEROQYSt3YSl3YBd3V\nRqMWPXv+0V2t1IX9+DZZrker5XKI1WvckydPxvHjx3Hnzh3UrFkTqamp6NSpU6mDe9q0aZzAhciJ\n7LkmzevZ9LjiQlopfLVabaHwVQpkGUJazazex33s2DHs2bMH3bt3R0xMDDZv3gwhhLW3EZGb2HOP\ntLfeX+3pE8LYQ2nSEk+5H1pWVoM7NDQUvr6+qFWrFi5cuIA6derg2rVrrigbEUlG1gD01i8stmBI\nq4/V4H7yySexcuVKNG3aFJ999hn27NmDO3fuuKJsRCQZewJQKfTt/SIg6xcJR7J19jFLryX1sRrc\nsbGxeOaZZ9C4cWO8+OKL+Oc//wm9Xu+CohGRN1EKfaVten0ZTJpU9P1KIe1tLWlbpwO1tLylpyx7\n6ems/kQHBgaiQ4cOuHbtGrp06VLonm4i8l7uGsi2e7cOWi0QFVV0uzvKoyZKk6DYOvsYoDyC3NP0\n6WNwdxHsZjW44+LisHnzZgQF5V/XEEJAo9HgwIEDTi8cETmGXl8G/v5Fw84etgYlR6rbT2mmMqWp\nQ229dUtp9jHA0kxlnsUTfg6tBveBAwfw1VdfISAgwBXlIXILTw8XS61UVx0b8Nxz62j2TCcaFmb7\ndKBK29w5daillrAntJAdzWpw16tXD7m5uQxu8mhK4eKMVip5r+JazY8uzGFPd7fSiliWJjdR2ubO\niVAsfbFzxRc+2b4cWA3uvn374sUXX0TdunXh4/PHfOHr1q1zasGI3M2drVSSg61d2PYuUWlrd7cs\nM5CpjWy9QVaDe86cOZg2bRqefvppV5SHiEiV7OnCtneJypJ0d5PrubrFbjW4y5cvj379+rmiLOSh\nPP36sTvx3DqHrSFdkhHb9ixRWZLubleRrXvZmVz9+2c1uJs3b45x48ahffv28PX1NW9nmJOtODjJ\neXhu7WdPSJdkxHZxreZHF+aQpbtb6WeuTx8D/P25YJSzWQ3uBw8eIDAwED/88EOh7QxuopJh61id\n7AnpknRhF9dqfnxhDlnp9TkICfFDerq7S+LZrAb33LlzXVEOIo/H1rFrOfreZ1m6sL2VN3XdWwzu\nt99+GytXrkTnzp2h0WiKPM8JWIhILVxx77MsXdjeypu+EFsM7piYGADA+vXrXVYYosfZ2r3MbuiS\nk/GcFQT0o/c9W7rVyhn3PpNjuKp17KmtcIuLjISGhgIA5s2bh6pVqxb67/3333dZAcl72LNIhLct\nJuEIajpntq5gVbAIRm5u4UUwLHV3A+L3bfktaaVtXLbS9fT6HLu+MPbpY7AplO09jlpZ/K0dM2YM\nUlNTcePGjUILixgMBlSpUsUlhSPv4uhrwJZalDK2ND2JPd3allrRvPfZu3j7767F4I6Li0NmZiZi\nY2MRHR39xxt0OgQHB7ukcCQXtU0RaumLAAeJOYetg8Hs6da2tAgGB44Vz51dxp7aXe1OFoM7MDAQ\ngYGBWLx4MX755RfUr18fu3fvxrlz5zB8+HBzVzpRAVdNEcoWszq5YkrPgoB+9L5nwL1zbLuTraHo\nzt8V/p46ntULXJMmTULNmjWRk5ODhIQE9O3bF1OmTMHatWtdUT6iIthiVidXTOlZENCect+zvfg7\n4J0sDk4rcOXKFYwfPx579+7FoEGDMGbMGNy5c8cVZSMilVIaOGbrYLD4+M4oVy4Lvr7Z6Nt3faFA\nfnSAGAeNESmz2uI2Go24ffs2Dhw4gISEBKSnp+Phw4euKBt5MHZ3y0HpurWl269K2moGvKtb21V4\nTdnzWW1xjxgxAq+88go6dOiAunXrIiIiAmPGjHFF2ciDqelWpJJQumVNRiW9/Sovz/rtV2w1q4On\n3gJFf7D6l7NPnz7o06ePuXv8X//6F7Raq3lPbsTWrPN4yvV1Z9x+RUSuYTWBU1NT0aNHD/Tt2xdp\naWno3r07zp4964qyUSnJ2pol11EKZFsnMQGgeJ1aZrZO6OHo98rCG+ooE6t/3WNiYrBs2TK89957\nqFKlCmbOnIkZM2Zg69atriifx1Pbvc9UcjL2cNizNjTgebdf2fPZyfS520IpoD2tjrZQ8xKlNi3r\nWatWLfPj1q1bY968eU4tlDdx1b3P5DxK3efOmGO9tF8QlOb35kAydXNF69bSMbwxpJWoeYlSq8Ed\nFBSE1NRU8wphu3btQoUKFZxeMCKZ2XotvCTXzG15ra3Tia5e3d/jZhVTCiJZu3ddEZ4MaHlZDW69\nXo/Jkyfjp59+QosWLRAWFoYFCxa4omxEVAxbQ9ryxCieRSmI3BVOsn5hIDlYDe7q1atj48aNyM7O\nhslkQmBgoCvKRSoi4zVcb2BrSHvaCHA1X3sswN8Vciab7+vy9/dnaLuI2u4V5ih191O6x9rWUeAF\nI8D9/DxjBLhen4P5891ditLh6GxyBN6QrUIlCUpbQ17pdWr7guDJbJ3wxNJrlSZCKel0olWq+HMS\nFDfj5CjkCGxGSc6eQVD2jIYmy+xZb3r16v42d4Fv2tTco0aByzC4TG3lIe9kNbjPnTuHxMRE3Llz\nB0II8/Z169Y5tWDkHp4yM5ir2D9ArOg2W69T2xvSrgohe5aeVNvPoTvLwy8NVMBqcE+ePBnh4eGo\nU6eO+ZYwIk+ndO+zECh1SNs64QmgPDmKpfuu7eGqEFJb+MqK55EKWA3usmXLIiIiwhVlIXI6pRay\nrYEM2DaXtz3rTQNQfK3MXeBE5FhWg7tt27ZYv3492rZtizJl/hjI9PTTTzu1YKRuar8WXpLlKAHb\nF9cobUhbCl6lbZZeK0NXqQxlJJKd1eDeuXMnAODjjz82b9NoNDhw4IDzSkWq565r4fa0mFev7m9X\nINs6l7ezWsdq/ZL0KBnKaCt+CSG1shrcBw8edEU5yAq1t3Dt5YoubKDki2sYjVr07Mm5vL2Rp/6u\nkfysBvft27cxa9YsHDlyBEajES1btoRer0flypVdUT7VcVeAetJob3tul7Kvxax8/RiwHMhardap\nc3mzVed6POckO6vBPX36dDRt2hSzZ8+GyWTCpk2bMG3aNKxcudIV5VMdTwpQd7HndilnLUfprkCW\n+edI1gBU2zmXYQpXUherwf2///0PS5cuNT8eOXIkdu3a5dRCkedQal3bc7uULMtRqi0cnMEb6ugK\nal4+ktTJ6pSnGo0GaWlp5sfXrl2DTscJ17xZcdN3pqVlF5q+056pOoubvvOpp8qap+9U2kZE5Kms\nJvD48eMRHh6OJk2aQAiBkydPIiYmxhVlIxWwd/pOe6fq5MAvIqLCrAZ3p06d0KRJE5w6dQomkwkz\nZ85EcHCwK8pWIjKMulZ7GZ0xfaczpupUG1mv9boTzxlR6VkM7kevaz/q3LlzAICxY8c6p0SlJMOg\nMbWX0RnTdzpjqk61UevnqWY8Z0SlZ/Eat7+/P/z9/XHu3DkcPnwYgYGBCAoKwtGjR3Hx4kVXlpGc\nwBHrOyu1Hj0VAAAUwUlEQVRdj358zWdefyYiciyLLe433ngDALB371588skn5ulOX3nlFQwbNsw1\npXMjtXdrl4StXeBhYfbNDKZ03zMRETmW1WvcmZmZMJlM5se5ubm4e/euUwvlakohrfZu7ZLw1vWd\niYg8kdXgDg8Px4ABA9CxY0eYTCYcOnQIw4cPd0XZXEbWkLZ1mlBXre/sDTx9UJWn14/IE1gN7jfe\neAMvvPACUlJSoNFosGTJEtSvX98VZXMKT+oCt/W2LFu7wMk6tf3c2Bq0ts7OpVQ/hjmRulgN7tzc\nXKSlpaFSpUoAgLNnz2Lv3r0YP3680wtnL0/pAre0RKWtt2WxC9xz2fpzbM/sXDL9rhB5A6szp40d\nOxbr1q3DwoUL8fXXX2Px4sXIzMy0+8AnT55EZGSk3fspzu7dOnNQy0JptLfS7GMAFEd8K21zxsju\nPn0MbIkREbmB1VS7ePEi9u3bh9jYWAwcOBBRUVGYMWOGXQdds2YNdu7ciYCAALv2o2al7ZIvyYQn\nts7b7QxshRERuYfVFndwcDA0Gg2effZZXLhwAU8++SRu3rxp10HDwsKwbNkyu/ahdqVt7dt6LzWg\nfI8075smIvJsVoO7Tp06iImJwV/+8hf8/e9/x6pVq5CTY19rq1u3bvDx8bFrH57K1glPSD68vEBE\njmC1SajX63H8+HHUrl0b48aNw5EjR/DRRx+5omxmFSv6Q6crGvQhIeXN/9Zqnb/N0vZJk/L/P39+\n8a97fNutW5m4dSsbBoMWY8fuxooVvbB2bV9Uq5bf1d2v306sWPEyKlUKQmho/nt37HjNanlcdS4e\nV5Jz5s5yO7retvqjk8n66G57jl3cex1VFzXwpLoAnlUf1sW5ig3uX375BQEBAWjRogUAoEuXLmjS\npAkWL17skBXChBDWXwQgIyO7yLaQkPJIT79nfmwy5V8vT0/Pcto2S9s3bcrfFhVl+f23b2fixo2K\nvwfyFsTHd0ZU1CFkZ+dfz968eShycvJX1AoOzh/tvXTpSzAagfT0eyUqj6vOxeNMpgBotdpCn4s7\ny2ip3L17l/l9e06xr7VUH1ew9ZyX5L2P/87IzJPqAnhWfVgXxx3bEotd5QkJCRg4cCB69OiBb775\nBkajEatWrUL37t1x7do1hxRMo9FYf5GEbB0Zbvl2LnImvT6Hg+uISFoWW9w7duzA3r17cePGDSxZ\nsgSrV6/GzZs3sWjRIrRr187uA1etWhWfffaZ3ftxN3uWwlReZcu9eA2WiEjdLAZ3QEAAQkNDERoa\nilOnTqFfv35Ys2aNVw8qc3RIF9zOZTRq0bOnOmYvc9XMWfbssySzhXkjb603kbewGNxa7R+96BUr\nVsSUKVNcUiA1c0RIK81epvYVtZzRrWzPPksyW5g38tZ6E3kLi9e4H73+XLZsWZcURi2UrlEDtt9j\nrXT7Fu+vdi7eakVE3sJii/unn35Cly5dAADXr183/1sIAY1GgwMHDrimhG6g1LJevbp/iVvSAOcB\ndxW2MonIW1gM7r1797qyHKpSkilGZQ5ptbVQ1VYeIiI1shjcVatWdWU53EZpwJml0d6WQlrWwFFb\nK1VtS0rauhQmEZErybV0lhModYuXdK1qtQWgEn65KN2xS7sUJhGRs3hVcCu1rpW6xR3R/e3ooLR3\nfzJ8uSAiIuu8KriVWtdhYcIpk6A4OigZvM4la48EEXkfrwpupdb1pk3NXbJ+NakbvxgRkSy8KriV\nBp3JPCocUFdLkYO5iIicz6uCu6SDzmSgppYiB3MRETmflMGt15eBvz8QFWX5NUoD0dzdulZT65iI\niORkccpTNdu9W4ctW/54bOsymu7G5SSJiMheUra4H2fr4h9ERESyk7LF/ThbF/8gIiKSnUcEt60r\ndJF7cOUuIiLH8Yiuci7+4dj3Oxqv6xMROY5HBLfMIa3E3qBjUMpPbV++iEg9PCK4iTwNv3wRkSUe\ncY2biIjIW7DFbQG7KomISI0Y3Bawq5KIiNSIXeVEREQSYXATERFJRLrgLpiXPC0t2zwvuStxMhEi\nInIn6a5xK81Lbu+92yUJYl77JiIid5IuuO1dPEQppBnGREQkC+mCOyzsDk6cEMgP75IvHsKQJiIi\nmUkX3AXzkhuNWvTsuVX6xUP69DHA39/P3cUgIiJJSBfcBfOSa7Vaj5mXPCTED+np7i4JERHJQLpR\n5URERN6MwU1ERCQRBjcREZFEGNxEREQSYXATERFJRLpR5SXBqUmJiMjTeHRwc7IVIiLyNOwqJyIi\nkgiDm4iISCIMbiIiIolIeY2b83sTEZG3kjK4Ob83ERF5K3aVExERSUTKFrcS3rNNRETewGOCm/ds\nExGRN2BXORERkUQY3ERERBJhcBMREUmEwU1ERCQRBjcREZFEGNxEREQSYXATERFJxOX3cQshoNfr\nceHCBfj5+SE2NhbVqlVzdTGIiIik5PIW9xdffIHc3Fx89tlneO+99zB37lxXF4GIiEhaLg/u77//\nHu3atQMANGnSBGfOnHF1EYiIiKTl8uC+f/8+ypcvb36s0+lgMplcXQwiIiIpufwad2BgILKyssyP\nTSYTtNrivz9UrOgPnc6nyPaQkPIKr5YT66JenlQf1kW9PKk+rItzuTy4mzVrhkOHDqFHjx44ceIE\n6tata/U9GRnZRbaFhJRHevo9ZxTR5VgX9fKk+rAu6uVJ9WFdHHdsS1we3N26dUNycjKGDBkCAByc\nRkREVAIuD26NRoOZM2e6+rBEREQegROwEBERSYTBTUREJBEGNxERkUQY3ERERBJhcBMREUmEwU1E\nRCQRBjcREZFEGNxEREQSYXATERFJhMFNREQkEQY3ERGRRBjcREREEmFwExERSYTBTUREJBEGNxER\nkUQY3ERERBJhcBMREUmEwU1ERCQRBjcREZFEGNxEREQSYXATERFJhMFNREQkEQY3ERGRRBjcRERE\nEmFwExERSYTBTUREJBEGNxERkUQY3ERERBJhcBMREUmEwU1ERCQRBjcREZFEGNxEREQSYXATERFJ\nhMFNREQkEQY3ERGRRBjcREREEmFwExERSYTBTUREJBEGNxERkUQY3ERERBJhcBMREUmEwU1ERCQR\nBjcREZFEGNxEREQSYXATERFJhMFNREQkEQY3ERGRRBjcREREEmFwExERSYTBTUREJBEGNxERkUQY\n3ERERBJhcBMREUmEwU1ERCQRtwX3/v378d5777nr8ERERFLSueOgsbGxSE5ORoMGDdxxeCIiImm5\npcXdrFkz6PV6dxyaiIhIak5tcW/duhX/+Mc/Cm2bO3cuevbsiZSUFGcemoiIyCNphBDCHQdOSUnB\npk2b8OGHH7rj8ERERFLiqHIiIiKJMLiJiIgk4rauciIiIio5triJiIgkwuAmIiKSCIObiIhIIm6Z\nOa20hBDQ6/W4cOEC/Pz8EBsbi2rVqrm7WKVy8uRJLFiwAOvXr8fly5cxZcoUaLVa1KlTBzNmzHB3\n8WxiMBjw/vvv4+rVq8jLy8OoUaNQu3ZtKesCACaTCdHR0bh48SK0Wi1mzpwJPz8/aesDALdu3cLA\ngQPx8ccfw8fHR9q6DBgwAIGBgQCAZ555BqNGjZK2LqtWrcLBgwdhMBgQERGBZs2aSVuX7du3Iykp\nCRqNBjk5OUhNTcUnn3yCOXPmSFcfIQSmTZuGixcvwsfHBzExMer9nRES2bdvn5gyZYoQQogTJ06I\n0aNHu7lEpbN69Wrx0ksvifDwcCGEEKNGjRLHjh0TQggxffp0sX//fncWz2bbtm0Tc+bMEUIIcefO\nHdGxY0dp6yKEEPv37xfvv/++EEKIo0ePitGjR0tdn7y8PDFmzBjRvXt38csvv0hbl5ycHNG/f/9C\n22Sty9GjR8WoUaOEEEJkZWWJxYsXS1uXx82cOVNs3rxZ2vp89dVX4t133xVCCJGcnCzGjRun2rpI\n1VX+/fffo127dgCAJk2a4MyZM24uUemEhYVh2bJl5sdnz55FixYtAADt27fHkSNH3FW0EunZsyfG\njx8PADAajfDx8cG5c+ekrAsAdO3aFTExMQCAa9euoUKFClLXJy4uDkOHDkVoaCiEENLWJTU1FdnZ\n2RgxYgRef/11nDx5Utq6/Oc//0HdunXxzjvvYPTo0ejcubO0dXnU6dOn8d///heDBw+W9u9ZmTJl\ncO/ePQghcO/ePeh0OtV+NlJ1ld+/fx/ly5c3P9bpdDCZTNBqpfr+gW7duuHq1avmx+KRO/ICAgJw\n7949dxSrxMqVKwcg/3MZP348JkyYgLi4OPPzMtWlgFarxdSpU7F//34sXrwYycnJ5udkqk9SUhKC\ng4PRpk0bJCYmAsi/FFBAprqULVsWI0aMwODBg/Hrr79i5MiR0v7OZGRk4Nq1a1i5ciX+97//YfTo\n0dJ+Lo9atWoVxo0bV2S7TPVp3rw5cnJy0KNHD2RmZiIxMRHfffed+Xk11UWq4A4MDERWVpb5sYyh\nreTROmRlZeGJJ55wY2lKJi0tDWPHjkVERAR69+6N+fPnm5+TrS4F5s6di4kTJ2LQoEHIyckxb5ep\nPgXXHZOTk3HhwgVMnjwZGRkZ5udlqkuNGjUQFhZm/ndQUBDOnTtnfl6mugQFBaFWrVrQ6XR49tln\nUaZMGVy/ft38vEx1KXDv3j38+uuveP755wHI+/dszZo1aNasGSZMmIDr168jMjISeXl55ufVVBep\nUq9Zs2b48ssvAQAnTpxA3bp13Vwix2jYsCGOHTsGAPjqq6/QvHlzN5fINjdv3sSIESMwadIk9O/f\nHwDQoEEDKesCADt27MDKlSsB5HebabVaNGrUyLwgjkz12bBhA9avX4/169ejfv36iI+PR7t27aT8\nbJKSkjBv3jwAwPXr13H//n20adNGys+lefPm+PrrrwHk1+XBgwdo2bKllHUpcOzYMbRs2dL8WNa/\nAdnZ2eYBkOXLl4fBYEDDhg1V+dlI1eLu1q0bkpOTMWTIEAD5LSNPMHnyZHzwwQfIy8tDrVq10KNH\nD3cXySYrV67E3bt3sXz5cixbtgwajQbTpk3D7NmzpasLAPTo0QNTpkxBREQEDAYDoqOjUbNmTURH\nR0tZn8fJ+nM2aNAgvP/++xg2bBg0Gg3mzZuHoKAgKT+Xjh074rvvvsOgQYPMd8lUrVpVyroUuHjx\nYqG7e2T9ORsxYgSmTp2KV199FUajERMnTsRzzz2nys+GU54SERFJRKquciIiIm/H4CYiIpIIg5uI\niEgiDG4iIiKJMLiJiIgkwuAmIiKSCIObSCWys7Mxa9YsvPjii+jXrx8iIiJsnht57969mDp1qsPL\ndObMGURGRjp8v9YkJCTg+++/d/lxiWTA4CZSiVGjRsHPzw+ff/45duzYgWnTpiEqKso8C5W7aDQa\nlx8zJSWl0BzeRPQHqWZOI/JUKSkpSEtLw7p168zbGjRogNGjR2P58uX4+OOPERkZiaCgIPz3v//F\nwoUL8eOPP2LFihUICAhA9erVUaZMGQDAqVOnMG/ePDx8+BAVK1bErFmzULVq1SLvT09Px5IlS2A0\nGvHMM88gJiYGFSpUQHJyMubNmwc/Pz/Url1bsbypqamYPn06Hj58iAoVKmDBggV48sknkZiYiN27\nd8PHxwdt2rRBVFQUrl27hsjISBw8eBAAsHTpUgDA2LFj0bZtW/To0QPff/89dDodFi1ahGPHjuHM\nmTOIjo5GQkIC3n77bRw6dAhA/vSaq1atwurVq535cRCpGlvcRCpw+vRpNGrUqMj2559/HqdPnzY/\nrlevHv71r3+hUqVKiI+Px4YNG7BlyxY8fPgQAJCXl4cPPvgAH330EZKSkjB8+HBER0cXeX9oaCg+\n/PBDrF27FklJSWjTpg3mz5+P3NxcTJ48GYsWLcK2bdvMczc/buLEiRgzZgx27dqF3r17Y926dfjy\nyy9x+PBhbN++HTt27MClS5ewceNGAJZb7Tdv3kTr1q2xfft2tGjRAhs2bEC/fv3QqFEjxMbGom7d\nunjmmWdw9OhRAMD27dsxYMCA0p1kIg/BFjeRCmg0GhiNxiLbH12dCMhfhx4Ajh8/jmbNmiE4OBgA\nMGDAABw4cAC//vorLl++jNGjR5uXvszOzi7y/lOnTiEtLQ2vvfYahBAwmUwICgrCjz/+iNDQUNSq\nVQtA/jzhj68JkJGRgfT0dHTo0AEAzGsHxMXFoXfv3vDz8wMADBw4EDt37jS/zpK2bdsCAOrUqVNo\nGcWC8hfsp0mTJvj2228xc+bMYvdH5OkY3EQq0LhxY2zYsAFGoxE+Pj7m7cePH8ef/vQn8+OyZcsC\nyA/6R68BF7zHaDSievXq2L59O4D88EtPTy/yfqPRiObNm2P58uUAgNzcXNy/fx9paWmF1rp+tCwF\nfH19Cz3Ozc3F9evX8fiyB0IIGAwGaDSaQs/l5eUV2kdB0D/+ugI9evTAwoUL8e9//xsdOnQocnwi\nb8OuciIVaNGiBWrXro05c+bAYDAAyB/RnZiYiDFjxhR5ffPmzXHixAlzYP7zn/8EANSsWRN37twx\nt1y3bt2KiRMnFnl/kyZNcOLECfz6668AgGXLlmH+/PmoV68ebt++jfPnzwOAeb+PCgwMRJUqVcwj\n3nfs2IGEhAS0atUKe/bsQU5ODgwGA5KSktCyZUs88cQTuHv3LjIyMpCbm2te1rI4Op3OfB7Kli2L\n9u3bY9GiReblY4m8GVvcRCqxdOlSfPTRR3jppZeg0+nMg75atGgBoPB14uDgYOj1egwfPhzlypVD\nvXr1IISAn58fFi1ahNjYWOTm5iIwMBBxcXFF3l+5cmXMmTMH7777LkwmE5566inMnz8fOp0OCxcu\nxJQpU+Dj41Ootf+o+fPnY8aMGYiPj0fFihURHx+PypUr4/z58xg4cCCMRiPatm2LiIgIaLVajBgx\nAgMHDsTTTz9t7q5/vEyPateuHfR6PeLi4vDnP/8ZvXr1wvHjx9G4cWO7zzOR7LisJxGpmslkwsKF\nCxEcHIzXX3/d3cUhcju2uIlI1QYOHIhKlSphxYoV7i4KkSqwxU1ERCQRDk4jIiKSCIObiIhIIgxu\nIiIiiTC4iYiIJMLgJiIikgiDm4iISCL/D3j3aUg34JEiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1097b2c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get row order of estimates as an index: low to high radon\n",
    "order = unpooled_estimates.sort_values().index\n",
    "\n",
    "# Plot mean radon estimates with stderr, following low to high radon order\n",
    "plt.scatter(range(len(unpooled_estimates)), unpooled_estimates[order])\n",
    "for i, m, se in zip(range(len(unpooled_estimates)),\n",
    "                    unpooled_estimates[order],\n",
    "                    unpooled_se[order]):\n",
    "    plt.plot([i,i], [m - se, m + se], 'b-')\n",
    "plt.xlim(-1, 86)\n",
    "plt.ylim(-1, 4)\n",
    "plt.title('Variation in county radon estimates')\n",
    "plt.xlabel('Ordered county')\n",
    "plt.ylabel('Radon estimate');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this visual inspection, we can see that there is one county with a relatively low predicted radon level, and about five with relatively high levels. A single state-wide estimate of prevailing radon level could not represent this variation well, and this reinforces our suspicion that a pooled estimate is likely to exhibit significant *bias*."
   ]
  }
 ],
 "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
}