04-pooled_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",
    "# 4. A Pooled 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 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 complete pooling, where we treat all counties the same, and estimate a single radon level: \n",
    "\n",
    "$$y_i = \\alpha + \\beta x_i + \\epsilon_i$$\n",
    "\n",
    "* $y_i$: measured log(radon) in household $i$\n",
    "* $\\alpha$: prevailing radon level across the state\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 pooled model in `Stan`\n",
    "\n",
    "To build a model in `Stan`, we need to define `data`, `parameters`, and the `model` itself. This is done by creating strings in the `Stan` language, rather than having an API that provides a constructor for the model.\n",
    "\n",
    "We construct the `data` block to comprise the number of samples (`N`, `int`), with vectors of log-radon measurements (`y`, a `vector` of length `N`) and the floor measurement covariates (`x`, `vector`, length `N`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Construct the data block.\n",
    "pooled_data = \"\"\"\n",
    "data {\n",
    "  int<lower=0> N;\n",
    "  vector[N] x;\n",
    "  vector[N] y;\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we initialise parameters, which here are linear model coefficients (`beta`, a `vector` of length 2) that represent both $\\alpha$ and $\\beta$ in the pooled model definition, as `beta[1]` and `beta[2]` are assumed to lie on a Normal distribution, and the Normal distribution scale parameter `sigma` defining errors in the model's prediction of the output (`y`, defined later), which is constrained to be positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initialise parameters\n",
    "pooled_parameters = \"\"\"\n",
    "parameters {\n",
    "    vector[2] beta;\n",
    "    real<lower=0> sigma;\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we specify the model, with log(radon) measurements as a normal sample, having a mean that is a function of the choice of floor at which the measurement was made:\n",
    "\n",
    "$$y \\sim N(\\beta[1] + \\beta[2]x, \\sigma_e)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pooled_model = \"\"\"\n",
    "model {\n",
    "  y ~ normal(beta[1] + beta[2] * x, sigma);\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the pooled model with `pystan`\n",
    "\n",
    "We need to map Python variables (from `clean_data`) to those in the `Stan` model, and pass the data, parameters and model strings above to `Stan`. We also need to specify how many iterations of sampling we want, and how many parallel chains to sample (here, 1000 iterations of 2 chains).\n",
    "\n",
    "This is where explicitly-named local variables are convenient for definition of Stan models.\n",
    "\n",
    "Calling `pystan.stan` doesn't just define the model, ready to fit - it runs the fitting immediately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pooled_data_dict = {'N': len(clean_data.log_radon),\n",
    "                    'x': clean_data.floor_measure,\n",
    "                    'y': clean_data.log_radon}\n",
    "\n",
    "pooled_fit = pystan.stan(model_code=pooled_data + pooled_parameters + pooled_model,\n",
    "                         data=pooled_data_dict,\n",
    "                         iter=1000,\n",
    "                         chains=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting the fit\n",
    "\n",
    "Once the fit has been run, the sample can be extracted for visualisation and summarisation. Specifying `permuted=True` means that all fitting chains are merged and warmup samples are discarded, and that a dictionary is returned, with samples for each parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Collect the sample\n",
    "pooled_sample = pooled_fit.extract(permuted=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is an `OrderedDict` with two keys of interest to us: `beta` and `sigma`. `sigma` describes the estimated error term, and `beta` describes the estimated values of $\\alpha$ and $\\beta$ for each iteration. These can be obtained through key values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.38870138, -0.60100033],\n",
       "       [ 1.45212072, -0.58113051],\n",
       "       [ 1.33121831, -0.69810621],\n",
       "       ..., \n",
       "       [ 1.39174434, -0.65172936],\n",
       "       [ 1.34036587, -0.51205498],\n",
       "       [ 1.34826004, -0.55500121]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Inspect the sample\n",
    "pooled_sample['beta']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While it can be very interesting to see the results for individual iterations and how they vary (e.g. to obtain credibility intervals for parameter estimates), for now we are interested primarily in the mean values of these estimates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: 1.3639661176627098, beta: -0.5912920402438805\n"
     ]
    }
   ],
   "source": [
    "# Get mean values for parameters, from the sample\n",
    "# b0 = common radon value across counties (alpha)\n",
    "# m0 = variation in radon level with change in floor (beta)\n",
    "b0, m0 = pooled_sample['beta'].T.mean(1)\n",
    "\n",
    "# What are the fitted parameters\n",
    "print(\"alpha: {0}, beta: {1}\".format(b0, m0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a scatterplot with the mean fit overlaid, to visualise how well this pooled model fits the observed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFtCAYAAAAeffM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0VHX+//HnTCoQSEIIYBBiLAiCgrIqq7IoC7ZVJIKK\nIoq6dlFX/BqKAooxErsrsoprAV1YV0FYLKwKunv8WRANgkqRDkEgCSGkl7m/PyZlxhRCyv0wd16P\nczx6J3dm3p8T83nfT3dZlmUhIiIiAc9tOgARERFpGUrqIiIiDqGkLiIi4hBK6iIiIg6hpC4iIuIQ\nSuoiIiIOoaQucgTr1asXw4cPZ8SIEdX/PPTQQwAkJyeTn59Pfn4+119/ffV7Zs2axfLlyw/7u049\n9VQyMzNbLPb6LFu2jLFjxx7yviFDhvDjjz+2ejwiThJqOgARqZ/L5WLevHlER0fX+tmiRYsA2Llz\nJ2vWrKl+/auvvuKEE05o0nfZxc7vEgkmSuoiRzDLsqhvf6hevXrx1VdfMXnyZIqLi0lOTmbUqFGs\nXbuW9PR03G43gwcP5sknn2TlypV4PB569+7Ngw8+SLt27fj222959NFHcbvd9OnTp97vGTJkCJde\neimfffYZBw4c4K677uK7777jxx9/JCwsjNmzZxMfH8/GjRuZMWMGubm5uN1uxo0bx4gRIwB47rnn\nWLp0KbGxsXTv3r36s8vKyuqNT0QOn7rfRY5w1113HcnJyYwYMYLk5GRycnKAmtZuWloakZGRLFq0\niDFjxtC3b19SUlIYOnQoL7/8MqGhoSxcuJD33nuPzp0789RTT1FWVsa9997LpEmTWLhwIWeffTbF\nxcX1xlBaWsrixYt54IEHmDp1KuPGjWPx4sV07dqVRYsWUVFRwR133MF1113HkiVLePnll3nmmWdY\nvXo1n376KZ988glLlixhwYIFFBUVVX9uffGJSNOopS5yhKuv+72hHZ6rfvbZZ59x8OBBvvjiCwDK\ny8uJi4tjw4YNhIWFceaZZwJw4YUXMm3atHo/7/zzzwegR48exMfH07NnTwC6d+9Obm4uW7dupbS0\nlKFDhwLQuXNnLrjgAv773/9y4MABhg0bRps2bQAYNWoUr7/+eoPxiUjTKKmLHOGaczxDRUUFU6ZM\nYdCgQQAUFhZSUlLC7t27a31uSEhIvZ8THh5e/d+hobWrDY/HU+vzPB4P5eXluFwuv5/5fk998YlI\n06j7XSTAhYaG4vF4/K7Ly8sBGDRoEG+99RalpaV4PB6mTp3KM888Q8+ePbEsi//+978AfP7559Xd\n+k2RlJREeHg4n3zyCQB79uxh2bJlnH322Zxzzjl89NFHHDx4EI/Hw5IlS6rfV198ItI0aqmLHMEa\nmiVe9bP4+Hh69+7NxRdfzPz58znvvPOYOXMmpaWl3HnnnTz++OMkJydjWRa9e/cmJSWF0NBQZs2a\nVZ1Ee/fuTadOnQ47hiqhoaG88MILpKam8vzzz+PxeBg/fjxnnHEGABs3bmTkyJFER0fTq1ev6vfd\ncccdpKen14qvsd8rIv5cOnpVRETEGYy11C+//HKioqIAOProo3nsscdMhSIiIuIIRpJ6aWkpAHPn\nzjXx9SIiIo5kZKLcunXrKCws5KabbmLcuHGsXr3aRBgiIiKOYmRMfcOGDaxevZorrriCrVu3cvPN\nN7Ns2TLcbk3GFxERaSoj3e/HHHMMiYmJ1f8dExPDvn376NKlS533W5almbAiIiKHYCSpL1y4kPXr\n1zNt2jT27NlDQUEB8fHx9d7vcrnYt++gjRHaKz6+vWPL5+SygcoX6FS+wOXksoG3fE1hJKmPGjWK\nyZMnM2bMGFwuF4899pi63kVERJrJSFIPDQ0lPT3dxFeLiIg4lprHIiIiDqGkLiIi4hBK6iIiIg6h\npC4iIuIQSuoiIiIOoaQuIiLiEErqIiIiDqGkLiIi4hBK6iIiIg6hpC4iIuIQSuoiIiIOoaQuIiLi\nEErqIiIiDqGkLiIi4hBK6iIiIg6hpC4iIuIQSuoiIiIOoaQuIiLiEErqIiIiDqGkLiIi4hBK6iIi\nIg6hpC4iIuIQSuoiIiIOoaQuIiLiEErqIiIiDqGkLiIi4hBK6iIiIg6hpC4iIuIQSuoiIiIOoaQu\nIiLiEMaSenZ2Nueeey5btmwxFYKIiIijGEnq5eXlTJs2jcjISBNfLyIi4khGkvrMmTO5+uqr6dy5\ns4mvFxERcSTbk/rChQuJi4vj7LPPxrIsu79eRETEsVyWzZn12muvxeVyAbBu3TqSkpKYPXs2cXFx\ndoYhIiLiOLYndV9jx47lkUceISkp6ZD37tt30IaIzIiPb+/Y8jm5bKDyBTqVL3A5uWzgLV9TGF3S\nVtViFxERkeYLNfnlc+fONfn1xuXk5JKSsoLMzFgSEnJITx9CbGyM6bBERCRAGU3qwS4lZQWLF48F\nXIAFzGPOnGTDUYmISKDSjnIGbdvWAW9CB3BVXouIiDSNkrpBiYkH8LbQASwSE/NMhiMiIgFO3e8G\nTZo0gJUr08jNPZqYmJ1MnjzcdEgiIhLA1FI3KC3tOzIzJ1FYeB2ZmZN47LHvTIckIiIBTEndII2p\ni4hIS1JSN0hj6iIi0pI0pm5QevoQYF7lOvX9pKefZzokEREJYErqBsXGxjBnTrLjtzsUERF7qPtd\nRETEIZTURUREHEJJXURExCGU1EVERBxCSV1ERMQhlNQN2rRpG/37/5V27ebSv//zbNmyzXRIIiIS\nwLSkzaCRI5eQmTkJcFFYaJGcnEZGxnjTYYmISIBSS92g/fuPxnebWO+1iIhI0yipGxQbuwPfbWJj\nY3eaDEdERAKcut8NWrToMpKTa45eXbRIR6+KiEjTKakblJSUSEbGeG0TKyIiLULd7yIiIg6hpC4i\nIuIQ6n43KCcnl5SUFZVHr+aQnj6E2NgY02GJiEiAUlI3KCVlBYsXj8W7rM0C5jFnTrLhqEQk0KnB\nELyU1A3atCkEWABEAQcrr0VEmkcNhuClpG5QTs42wLujHFjk5KQZjkhEnGDbtg74bmzlvZZgoIly\nBkVGdsf3D69Nm+4mwxERh0hMPIDvxlaJiXkmwxEbqaVu0NatP+H9w/O21Lds+RkYYTYoEQl46elD\ngHmVY+r7SU8/z3RIYhMldYM8nm7UjKnn4/EkGI5IRJwgNjaGOXOStbFVEFJSNygy8leKi++kqqUe\nGTnVdEgiIhLAjCR1j8fDgw8+yJYtW3C73Tz88MMcf/zxJkIxat68QVx11RQ8nuNxu3/hzTcHmw5J\nREQCmJGJcsuXL8flcjF//nzuuecenn76aRNhGPfmm3vweFKBG/F4Upk3b4/pkETEAXJycrn55kWc\ncca/ufnmhezfn2s6JLGJkZb60KFDGTJkCAC7du0iOjraRBjGaZ26iLQGrVMPXsbG1N1uN5MmTeLj\njz/m+eefNxWGUVqnLiKtQevUg5fRdeppaWksW7aMBx98kOLiYpOhGBEVdQy+f3jeaxGR5jnqqN3A\nP4B/A29VXkswMNJSf++999izZw+33norERERuN1u3O6Gny/i49vbFJ19tm/3X6e+fftPxMdfYziq\nlufE350vlS+wObF8ERFtgaupqlsiIuY6spxOLFNzGUnqF154IRMnTuTaa6+lvLycKVOmEB4e3uB7\nnLjWsrg4Dt916sXFcY4rp9PXyap8gc2p5du6NRrfXsCtW6MdV06n/u6qNPWBxUhSj4yM5NlnnzXx\n1UeUsLACyspGU/U0HRb2mOmQRMQBEhMPkJFR0wuobWKDhzafMWjgwDb873/zgfbAQX7/+zamQxIR\nB9A2scFLSd2gqKi2ftft2impi0jzaZvY4KWkbtCOHfHAqOrrnTvfMReMiIgEPCV1g/bs+Q4oxdv9\nnsevv/4IXGA2KBERCVhK6gbl5HjwXXaSkzPJcEQiIhLIjG4+E+wqKo7Bd9mJ91pERKRp1FI3aje+\nm8/Ar2bDERFHyMnJJSVlReXs9xzS04cQGxtjOiyxgZK6UTnA40A3YBeQbTYcEXGEe+9dxkcf3URV\ng6G09O+88cZVpsMSGyipGxUPTKSmpT7DbDgi4ghffeXGd2jvyy810hoslNSNaovv0augdeoi0hKy\n8B/aUy9gsFBSN+pXYAI1f3gPmA1HRBzh97+P4sMPfXerjDIdkthESd2o4/DtIoNjDcYiIk7x7LOX\nEh6+gsxMSEgoJz39EtMhiU000GLUXrwtdCr/vc9gLCLiFJbld2UqDDFALXWj8oCaLjI4YDYcEXGE\nlJQVLF48lpqhvXnMmZNsOCqxg1rqRhUe4lpE5PBt29YB36E977UEAyV1o7oB24H9wDbgKLPhiIgj\nJCYewHdoT+epBw91vxuVB8ykpossxWw4IuIIkyYNYOXKNHJzjyYmZieTJw83HZLYREndqGPR7HcR\naWnTp39FZuYkwEVhocX06X/njTcSTYclNlD3u1Ga/S4iLU87ygUvtdSNysZ/9nuW2XBExCG0o1yw\nUlI3KhbvRLkEvAe66BQlEWm+AQNC+PTTmsOiBgwINx2S2ERJ3ahiYBo1T9OTzYYjIg7RDriPqrrF\n5XrJcDxiFw20GJWE/0Q5TWQRkeZbuTIc37rlm2/UUg8WaqkbtQ74B94x9TxgvdlwRMQRiooy8R1T\nLyrabTgisYuSulFtgaup6X7/yWw4IuIIbdpUUFZWMwm3TZsK0yGJTdT9blQ3/LvfuxmMRUSc4vTT\no/DuUundrfL009sZjkjsopa6Udlo2YmItLTw8GjgJqrqlvDwvxuOSOyipG5UFv7r1HPMhiMijrB7\ndyd8ewG91xIM1P1ulOc31+VGohARZ9GBLsFLLXWjTgKu8bnWeeoi0nzp6UOAeWRmxpKQsJ/09PNM\nhyQ2UVI3ah/+Y+raJlZEmi82NoY5c5KJj2/Pvn0HTYcjNrI9qZeXlzN58mR27dpFWVkZt912G0OG\nDLE7jCPEHmAq3k1oNuOdqSoi0jw5ObmkpKyobKnnkJ4+hNhYbUMdDGxP6kuWLCE2Npb09HQOHDjA\niBEjgjiptwMeQeepi0hLSklZweLFY6mpW+YxZ06y4ajEDrYn9YsuuogLL7wQAI/HQ2hoMI8AJAAL\ngCi8s9+PMhuOiDjCpk0h+NYt3msJBrZn1DZt2gCQn5/PPffcw1/+8he7QziC/ArcQ83T9CSz4YiI\nI+TkbMNbn3jrlpycNMMRiV2MNJN3797NXXfdxbXXXsvFF1/cqPfEx7dv5ahMOA7/HeWOc2Q5nVgm\nXypfYHNi+aKijsG3bmnf/hhHltOJZWou25N6VlYWN910E1OnTmXgwIGNfp8zZ3D+jP/s9/WOK6fT\nZ9+qfIHNqeXbuHENvnXLhg1rHVdOp/7uqjT1gcX2pP7SSy+Rl5fHiy++yKxZs3C5XLzyyiuEhwfj\n0YCF+O8ol282HBFxBMvyna+Tj2Vpvk6wsD2pT5kyhSlTptj9tUeoE6jpIgM43lQgIuIgLtceLOsu\nqlrqLtdDpkMSmwTz1PMjwH7gPmq63/WHJyLNd8457fjf/x7He/LjLs45J8p0SGIT7f1uVGf8J8p1\nMRiLiDjFlCmDiYwsxOUqJzKygIce+oPpkMQmaqkbtRP/iXI7zYYjIo5w442fUVzs3diquNhi3Lg0\nMjJONh2W2EBJ3ah4fCezeK9FRJonKysG381nvNcSDJTUjSoGRlPTUn/UbDgi4gilpduB26mqW0pL\ntbFVsFBSNyoX/yVtOnpVRFrCMfjP1znGWCRiL02UMyoLWIv3CNY16OhVEWkZVfN1QPN1gota6kZ1\nAvribalHAKVmwxERR+jYsYicnJolbR07FpsOSWyilrpRbX5zHWkkChFxlg8/vIGEBIu2bSEhweLD\nD8eZDklsopa6UeHA1dRMlJthNhwRcYSkpEQyMsY7fn90qU0tdaPi8J/MEmcwFhFxik2bttG//19p\n124u/fs/z5Yt20yHJDZRS92orfhvPqM/PBFpvssuW8TevQ8CLgoLLS699FHWrr3XdFhig0Mm9bKy\nMpYuXcry5cvZunUrbrebxMREhgwZwp/+9CfCwsLsiNOhCoCaySw6pU1EWsLevf6bz+zdG204IrFL\ng0n9s88+Y/bs2QwYMIDk5GQSEhIICwtj586dfPXVV7z55pvcfvvt/PGPf7QrXofpBkykpqX+mNlw\nRMQhMoGaU9pAJ2MGiwaT+tatW3nzzTdrtcaPO+44Bg8eTGlpKW+99VarBuhsMfiPqWsrRxFpCS5g\nGt5NZzYDHqPRiH0anCg3bty4BrvXw8PDueGGG1o8qODxC/4bRGwyGIuIOEch0AvveRK98W5JLcGg\nURPl/ve///HMM8+Ql5eHZVlYloXL5eLTTz9t7fgcLg6YCSTgHVPvaDYcEXGITvgvl91iNhyxTaOS\n+qOPPsrEiRM54YQTcLlch35DAAv7bDltXn0ZT1wnrLhOeOI64YmLo6L3SZSf3K+Fvy0Pb1Kv+sNL\naeHPF5HgpOWywapRST02NpbzzjuvtWM5IoRuWEfERx/Uer3o+pvIf+KZWq+Hf/g+bd74e2Xy74Sn\nk/dhoPyUfo14COiO/x9e9+aGLyJCWNhmyspqlsuGhamlHiwaldQHDBhAWloagwYNIiIiovr1008/\nvdUCM6Xoz7dRfMVo3NnZuLKycOdk487OouL4E+q8P2TjesKXf1Lr9cLbx9eZ1CPefZvIf8zDExfH\nX9nEPj5jH535im58T3aLl0dEgk/Xrm3ZsWMS0BPYQNeu2oI6WDQqqf/www8A/PTTT9WvuVwu5s6d\n2zpRmeR2Y8V2pCK2I9STyH0Vjf8LRTfegjs7y/tPjvdhoKLniXXeH7J5E+H/+xzwLjiB1QA8zEN8\nT+1JiZFzXyPyXwtqegLi4rDi4ij7/dmUn9K/qaUUEQfbtasMSKOqpb5rl5a0BYtGJfV58+YBkJ+f\nj8fjoUOHDq0aVEBxuSAqCk9UFJ7EYw55e+H/TaLw7vtw789hyCnzieds4tnHT5wErKp1v3vnDkK/\n+QqXZfm9nv/QI3Um9TazXyBi8bveOQEd46ofBsrOPa8V5gSIyJHI4+mG7+YzHk+C4YjELo1K6jt2\n7OAvf/kLO3bswLIsEhISePbZZznmmGNaOTyHiojA0/UofqAQ+CM1E+Xm17q1cPJUClOm4MrNre4N\ncGVnU3Firzo/2r13D6Fr1+Aq9T/G9WDbJ+pM6m2fSCPig6WV8wHiqh8GSs6/iIqTT2l+WUXEgF+A\nquFRLZcNJo1K6lOnTuXPf/4zF154IQAffPABDz30UHULXpqqCKgZ9/KuLa1DSAhWXBwVcXFUUHe3\nfpWCaTMomPoIroJ875yAygeB8hN713m/K+8A7m1bCf1xjd/rFQnd6kzq7aY/SPgny36zOqAjJZcm\nU9Gn76GLLCI26ID/krYNZsMR2zQqqe/fv786oQNcfPHFzJ49u9WCCh5h+I57wf0t87EuF1ZUe6yo\n9niOSWrw1oIZj1Mw43EoLq6eFOjKyqq3J8BVXIQ7ax8hGzf4DQlU9Dqp7qR+6610/M/H1SsDPB3j\nsOI6UXzFaCp6n9SsYopIfbriv7Kmq8FYxE6NSurh4eH8+OOP9OnTB4C1a9fSpk2bVg0sOByL/x/e\nseZCiYzEk9ANT0K3Bm/Lf/wp8h9/CsrLce3fX9MT0LPuhwBcLigqIvSHDFxlZdUvlw38fZ1Jvf2t\nNxD2zdeVyT+u+mGgeOwN9U4+FJHf2oi3wVB1WNQ+s+GIbRqV1CdPnsz48eOJiYnBsiwOHDjA008/\n3dqxBYFM4B9Ae7wb0WSaDedwhIZixcdTER9PRUP3/e1v5Mx4AiwL18G86iGB+pYIWm3bgdtN6KaN\nuNasrn699IKL60zqHa4ZRdjqDDxxNZMCrbg4im65nYrj6vgOjwfcDe6OLOIA4fgfFjXJbDhim0Yl\n9f79+7Ns2TK2bt2Kx+MhKSmJ8PDw1o4tCJTiP+71gNlwWpPLhdUhGqtDNJ5jj6v3tvxnXqi5KCys\nXiZYUc97rE7xeNq3x717N6Hrfq5+vfiqa+q8P+ayiwhZ/3P1MEDVMsHCeybUPVRRVgY6XlgCjNt9\nAh5PTS+g233o5bniDA0m9UmTGn66S0tLa9Fggs/x+He/H28wliNQ27Z42vbA071HvbccfN5nbkdZ\nmXefgOxsKuqZS1DRIxHXgVzv5kLbtuKq8PYzFN1yR533x55/LiFbNuGpXB5YNSRQMHkqnm5H135D\nfj60a+cddhAxJCxsMyUl2lEuGDWY1M8991wAli1bRnFxMZdffjkhISEsXbrUb2c5aaqf8O9+/6nh\n26VhYWF4unSFLvVPCjo46+WaC4+nOsFX1PPgUN77JKyQENzZWYSu/xlXsfe0q4IHJtd5f+zQQYRs\n31bTE9DJ2xOQP2MmVpcute537c/Bat8BQhvVaSbSKHFx7cnMrDp6dQudOrUzHJHYpcGa5IILLgBg\nzpw5/Otf/6o+zOXcc89l1KhRzfri1atX8+STTwb5srgw/LvfvzcbTrDx3T2wHgdfnFNzYVlQUIA7\nO6vuVjpQ3v80rNiO3mGDnTsI/flHAPJTn6jz/tjBv8e951esmBi/ZYIHn34eq2PtQzjce37FEx0D\nkdr2U+qXmZkNpFOzo5yDh/bET6OaBwUFBeTk5BAX561k9u7dS1FRUZO/9JVXXmHx4sW0axfsT4+9\n8O9+r2cGuRwZfHYPrM/Bv/3d/4XSUtw52VhxdZ+SVXbW2bh3765ZTrh5Ey6Ph7wXXqrz/tizT8ed\ndwBPu6jKngDvcEDeS69BHXG5t2/Dio3FimqvIYGgcgStrBFbNSqp33777Vx22WWcdtppeDweMjIy\nmD59epO/NDExkVmzZvHAA8H+9LgZbwu9qqW+1Wg00grCw/F0PareHx/826v+L1RU4MrN9Y7L/5Zl\nUTr0/OpdBd3ZWYSuXeOdzFfXElPLouPvT8NVVoYVHl6zOqBjHAf+8S+oY7JryC8b8cR2xIqN1SqB\ngLYP/7oly2w4YptGJfXhw4czcOBAvv/+e1wuFw8//HB1q70phg0bxq5du5r8fufIAx6nZi1prtlw\nxLzK3QPr5HLV7gmwLFwF+RASUvv+0lJKLr8CV2UvgDsrm5Atm3Ft2lhnQqe0lI5nDfB+rNuN1bEj\nxMcT3bETBxYurd3StyxCflxbPXmwzs8UQ3bhX7cE0HJZaZZGJfXs7Gw++OADCgoKsCyLdevWsXPn\nTtLT01s7vmrx8e1t+y77hACJeCfKhQJ7HFlOJ5bJl/ny1XfAUntY8Fbtl0tKiK9romt+Ptx8M2Rl\n4dq3D1dWFuzdS3hWFvGd6/iO3FwYcrZPGB0gPh4SE+HTT2vfX14O333nvSc+/ohZJWD+99fyuneP\nZMeOqrolhO7df3VkOZ1YpuZqVFK/66676NGjBxkZGQwdOpQvvviCwYMHN/vLrd+cPNaQffsONvv7\njjzHAb7rqbMcV874+PaOK5OvwC1fad0vpz7ldxkf3559ew5AHWV05ebR7qZbcGVn4c7OqRwWyILN\nW8ip4373r7uJO/PM6msrMhJPXCcqjjuBA+8srh1LSQmh33+HVbmCwIqOafEhgcD9/TVsx444fCfh\n7tixyXHldOrvrkpTH1gavff7/PnzmTlzJueffz633XYb99xzT5O+0JfrCHhKN+sH4CEgCe/4+q9m\nwxGpSz2J1IqJJT/tyTp+UPfDuhUaRuFtd3mTf9WQQHY2rgN1DzuF7NpB7PALat4fEoLVMY7yPn05\n8PZ7td9QUEDYd9/W7CfQMS6INw6Kx3+iXLzBWMROjUrq0dHRACQlJbFu3Tr69evH/v37m/XF3bp1\nY8GCBc36jMDXCXgEbeUojlLPw7rVqRMFjzzW6I+x2raj8O77KnsCsmt6Airq3pg4ZPMmYkZe6vea\nJzqGst+dTt78d2uHmXcAVn1BaEib6kmEtG3b6PiObD/jvwfGOrPhiG0aldQHDhzI3XffTUpKCjfe\neCM//vgjYUH7BNySeuL/NN3TYCwiRxZP16MoeHB6o++3OnakYEJK9a6CVYcNEVH3mv6Qn3+GSy8i\n1vcz2ral9Jw/kPfm27Xud2VnE7ay8rChyqWEVofoI2JewG9FRMRRUlLT/R4R8YvpkMQmjUrq119/\nPfn5+XTr1o2nn36alStXcuedd7Z2bEFgNf5P06uBK4xGJBKoPN2OpjBlSuPvT0iA1FQKt+/CXXnQ\nkCsnB6uezYhCf8gg+rrRfq9ZoaGUnn8Rea/XnpDo3vMroV9/6XPGQCfvUkEbdg90uzvj22DwXksw\naNT/XWPGjOHDDz8EoE+fPtVHsEpzefDfUe4rs+GIBBFP9x4weTIFjZxsVZF0LPnTHq3sCciqXCaY\nRUVCQp33h676lug/X+/3muVyUTI8mYNzXq91v3vHdsK+/tKb/DtVPgh0jGvS7oFFRVvwXadeVLT1\nsD9DAlOjknqvXr147733OOWUU4j0+R8soZ7/maWxTsK/+732+eIicmTwHJNE0Z13N/r+8t4ncfCx\n9MpegJo5AXWeBgiEff0lHe64udbrxSOv5ODsV2q9HrJpI6HffF29T4B3WKCTd/dAjgYWAFFAPt71\n6hIMGpXUV69ezerVq/1ec7lcfFrXWlQ5DFn47/qUbTYcEWkxnqRjKf7zbY2+v+y033HwiWerk3/V\nhkHlvXrXeX/Y//uC9hNqP2QUXX0t0AMYT03d8hAhP/9E2Ddf+Z02WD0koN0DHaNRSX358uWtHUeQ\nygLm4x1TP4i2chQJXp5jj6P42OMafX/Z2eeQ99yL3p6AnJqegPKTT4H5m/HfUS6H8P+uIOqh2its\nisbdRH76M7VeD834jtBvv8Hq6PMA0KlySEC7Bx6xDnme+i233EJSUt3dRRs3buTVV1/VuepNVoh3\nrXpPYAPebjIRkUOrOPZ4Ko49vs6fhU17jbKyeLzd72GEhbWjdNgg8uI6+Zwd4H0QKO9/Wp2fEf7Z\ncto99kit1wvvuJuC6Y/Wej30qy8JW7XSu1GQz5wAT+cudZ9NIK2iwaR+7733kpqayr59+xgwYABd\nu3YlJCRiDP4tAAAV3ElEQVSEzMxMvv76a7p27crEiRPtitWB4oA0tE5dRFpSfPxBMjNrut/j49Ma\nfAioS8nwEZQfd3z1yoDqnoD+p9Z5f/jyj2n3bO3NiAompNS5KiHss+WEZXxXk/wrewIquh5V54mD\n0jgNJvUuXbrw/PPPs337dlasWMHmzZtxu910796dJ598kh49etgVp0PFATOBBLxdZLEN3y4i0ggz\nZvTippsmUdULmJp6+mF/xuE+BBSPHkP5ab+rmRNQ+TBQb0/AJ8to+/LsWq/nP/QIRePvrX3/Rx8Q\nuvp7PJ06YcV1gj+dD2Ha+/23GjWm3qNHD66//vpD3yiHKRu11EWkpd15Zwa+dcvtt09l+/bzW/U7\nPcceR+lhzAko+vNtlJ33x5qhgMqlguX9+td5f/h/PqTNm2/UvPDee3DWkOaG7TiNSuqDBw9m7969\ndOjgPakpLy+PDh06cPTRR/Poo4/Su3fdszPlUI7Hf0lb45+KRUTqU1zcA9+6xXt9ZPEck0RpPcv7\n6lJ47/2UJI+q7glof2rdwwDBrlFJ/fTTT+fCCy9k6NChAHz++ed89NFHjB07locfflh7uDfZRvyX\ntGkrRxFpCb/gX7dsMhtOC/D0SMTTI7H6un18+zpPDwx2jVqcuHHjxuqEDt6W+/r16znppJMoKSlp\nteCcLwbvmPrcyn9Hmw1HRByiG97NZ/5d+e+jzIYjtmlUS71Dhw4sWLCA4cOH4/F4+Pe//010dDSb\nNm3C4/G0dowOlot3LanG1EWkJe3Cf/MZ1S3BolFJ/cknnyQ1NZUnnniC0NBQzjrrLGbOnMmyZcuY\nMGFCa8foYJ3x38pRZx6LSEuIpWZlTSbeXkEJBo1K6l26dOGpp55i8+bNVFRU0LNnT0JDQxk7dmxr\nx+dw66jpFvMA6w3GIiLO8StwBt7dKkOAb82GI7ZpVFJfs2YN99xzDzExMXg8HrKyspg1axb9+vVr\n7fgcLhb/U9pWN3y7iEgjxMREkptbU7fExKwxHZLYpFFJPTU1lWeeeaY6iWdkZDBjxgzeeeedVg3O\n+Xrhv6Stl8FYRMQpCgv9l8sWFWm5bLBo1Oz3wsJCv1Z5//79Neu9RazH20Kn8t8bDMYiIk5RWlq1\npA3AoqREy2WDRaNa6tHR0XzyySfVy9o+/vhjYmI08aJlTMG76YzvH6GISHP89jz1o82GI7ZpVFKf\nMWMG//d//8eUKd5N+bt3784TTzzRqoEFhzy85x6H4p3Qst1sOCLiCGFh2ykru5OqMfWwMB28FSwa\nTOpjx47F5fKOy0RGRnL00UdjWRZt2rRh2rRpzJ0715YgnSsKmEjNRLkHzIYjIo5QVlYKzMfbWDhI\nWZmGS4NFg0l9/PjxdsURpHrjP1FOe+iLSEs4BbjG57rYVCBiswaT+hlnnGFXHEFqHf77M68zG46I\nOMRneOfpHI/3jIndwBUmAxKbNGpMXVpL1Tax3fBu65hrNhwRcYgEIBVtExt8lNSNGgjc6HP9qqlA\nRMRReuI/tNfTYCxiJyV1o74E9lDTUt+CushEpPl+RkN7wUlJ3ahY/Ge/p5gNR0QcIh//oT2dOx4s\nlNSN0jaxItIaTkdDe8FJSd2o1cA/8K4lzQPWoO53EWm+n1D3e3BSUjfKwv+Utq/NhiMiDnEQ381n\n4IDZcMQ2tid1y7KYPn0669evJzw8nNTUVLp37253GEeI/vh3v/c3GIuIOIXLdSaWdY3PtTafCRaN\nOqWtJX3yySeUlpayYMECJkyYQFpamt0hHEE2oFPaRKSlRURsxrduiYjYYjIcsZHtSX3VqlUMGjQI\ngH79+rF27Vq7QzhijB5djHdTiFeBSVxzTYHhiETECS67rICaumUil1+uja2Che1JPT8/n/bt21df\nh4aG4vF47A7jiLBgQSSQhneWahr/+Ec7wxGJiBP8858R1NQtj6tuCSK2j6lHRUVRUFDTIvV4PLjd\nh362iI9vf8h7Ak/tXZ+cWE4nlsmXyhfYnFk+1S3Byvakftppp7FixQouvPBCMjIy6NmzcdsX7tvn\nxM0T1uO/7GSD48oZH9/ecWXypfIFNueWT3VLoGvqA4vtSX3YsGF88cUXjB49GiCoJ8olJ+exaNEk\nvE/V6xk5UuNeItJ8f/rTXt5/v6ZuueSSvaZDEpvYntRdLhcPP/yw3V97RFq0KBrvuJf3afrddycx\ne7bhoEQk4L3/fhd865alS3VKW7CwfaKc+NJJSiLSGlS3BCsldaOqxr1A69RFpOVoD4xgpaRu0G23\nheG7Tv2OO0IMRyQiTvCXv3TAd536hAlRhiMSuyipG/Tyyx5816n/7W/WId4hInJozz2Xj+869Wee\nKTQckdhFSd0gj+d4fMe9vNciIs2juiV4Kakb5HZvxHfcy+3+xWQ4IuIQqluCl5K6QW+/fS5u9xTg\nVdzuKbz99mDTIYmIA6huCV46T92gbt260bVrFLm5ocTEtKN7926mQxIRB+jb9yQuvXQPmZmxJCSE\ncPLJJ5kOSWyipG7QyJFLyMycBLgoLLRITk4jI2O86bBEJMClpKxg8eKx1GwTO485c5INRyV2UPe7\nQfv3H43vZBbvtYhI82zb1gHfusV7LcFASd2g2Ngd+E5miY3daTIcEXGIxMQD+NYtiYl5JsMRG6n7\n3aBFiy4jOTmN3NyjiYnZyaJFw02HJCIOkJ4+BJhXOaa+n/T080yHJDZRUjcoKSmRjIzxjj9CUETs\nFRsbw5w5yapbgpC63w1avPg/dO6cisv1Kp07p/L++/8xHZKIOIDqluCllrpBN9+8Et/jEW+4YRJ7\n955vOCoRCXSqW4KXWupG6XhEEWkNqluClZK6UToeUURag+qWYKWkbtBrr52O79Gr3msRkeZR3RK8\nXJZlBcR5n06ewenkGapOLhuofIFO5QtcTi4beMvXFGqpi4iIOISSukGrVq2hR4903O5X6dFjJhkZ\na0yHJCIOoLoleGlJm0HJyR9SXPwI4KK42GL48Kls336y6bBEJMCpbgleaqkbVFKShO+yE++1iEjz\nqG4JXkrqBkVEbMZ32UlExBaT4YiIQ6huCV7qfjdoyZKLGT58KiUlSUREbGHJkotMhyQiDqC6JXgp\nqRvUv//JbN9+suOXZoiIvVS3BC91v4uIiDiEkrqIiIhDKKmLiIg4hMbURUQcJicnl5SUFWRmxpKQ\nkEN6+hBiY2NMhyU2MJbUP/74Yz766COeeuopUyGIiDhSSsoKFi8eS9V56jCPOXOSDUcldjCS1FNT\nU/niiy/o3bu3ia8XEXG0bds64Lv5jPdagoGRMfXTTjuN6dOnm/hqERHHS0w8gO/mM4mJeSbDERu1\nakv9nXfe4Y033vB7LS0tjYsuuohvvvmmNb9aRCRopacPAeZVjqnvJz39PNMhiU2Mnaf+zTff8M9/\n/lNj6iIiIi0kYGa/O3lXJCfv+uTksoHKF+hUvsDl5LKBt3xNoXXqIiIiDmGspX7GGWdwxhlnmPp6\nERERx1FLXURExCGU1EVERBwiYCbKOZG2chSR1qC6JXgpqRukrRxFpDWobgle6n43SFs5ikhrUN0S\nvJTUDdJWjiLSGlS3BC91vxukrRxFpDWobglexraJPVxO3znIqeVzctlA5Qt0Kl/gcnLZQDvKiYiI\nBD0ldREREYdQUhcREXEIJXURERGHUFIXERFxCCV1ERERh1BSFxERcQgldREREYdQUhcREXEIbRNr\nkI5HFBGRlqSkbpCORxQRkZak7neDdDyiiIi0JCV1g3Q8ooiItCR1vxuk4xFFRKQlKakbFBsbw5w5\nyY4/QlBEROyh7ncRERGHUFIXERFxCCV1ERERh1BSFxERcQgldREREYdQUhcREXEIJXURERGHsH2d\nen5+Pvfffz8FBQWUlZUxceJE+vfvb3cYIiIijmN7Un/ttdc466yzuO6669iyZQsTJkxg4cKFdoch\nIiLiOLYn9RtuuIHw8HAAysvLiYiIsDsEERERR2rVpP7OO+/wxhtv+L2WlpZG37592bdvHw888ABT\npkxpzRBERESChsuyLOvQt7Ws9evXc//995OSksI555xj99eLiIg4ku1J/ZdffmH8+PE8++yznHji\niY1+n5MPPHHygS5OLhuofIFO5QtcTi4beMvXFLaPqT/99NOUlpaSmpqKZVl06NCBWbNm2R2GiIiI\n49ie1F988UW7v1JERCQoaPMZERERh1BSFxERcQgldREREYdQUhcREXEIJXURERGHUFIXERFxCCV1\nERGHWbVqDT16pON2v0qPHjPJyFhjOiSxie3r1EVEpHUlJ39IcfEjgIviYovhw6eyffvJpsMSG6il\nLiLiMCUlSYCr8spVeS3BQEldRMRhIiI2A1XHelhERGwxGY7YSN3vIiIOs2TJxQwfPpWSkiQiIraw\nZMlFpkMSmyipi4g4TP/+J7N9+8mOP8lMalP3u4iIiEMoqYuIiDiEkrqIiIhDKKmLiIg4hJK6iIiI\nQyipi4iIOISSuoiIiEMoqYuIiDiEkrqIiIhDKKmLiIg4hJK6iIiIQyipi4iIOISSuoiIiEMoqYuI\niDiEkrqIiIhDKKmLiIg4hJK6iIiIQyipi4iIOISSuoiIiEOE2v2FRUVFTJgwgby8PMLDw3n88cfp\n3Lmz3WGIiIg4ju0t9bfffpu+ffvy5ptvcumllzJnzhy7QxAREXEk21vq119/PZZlAZCZmUl0dLTd\nIYiIiDhSqyb1d955hzfeeMPvtbS0NPr27cu4cePYsGEDr776amuGICIiEjRcVlWz2YDNmzdz6623\n8vHHH5sKQURExDFsH1N/6aWXeO+99wBo27YtISEhdocgIiLiSLa31LOzs0lJSaGkpATLspgwYQKn\nnnqqnSGIiIg4ktHudxEREWk52nxGRETEIZTURUREHEJJXURExCGOuKReUlLC3XffzZgxY7j11lvZ\nv39/rXtef/11rrzySq666ipmzZplIMrDZ1kW06ZNY/To0Vx33XXs2LHD7+fLly9n1KhRjB49mn/9\n61+Gomy6Q5Vv6dKlXHnllVxzzTVMnz7dTJBNdKiyVZk6dSpPP/20zdE136HK98MPPzBmzBjGjBnD\nfffdR1lZmaFIm+ZQ5fv4448ZOXIkV1xxBfPnzzcUZfOtXr2asWPH1no90OsWqL9sgVyv+KqvfFUO\nq26xjjCvvfaa9de//tWyLMt6//33rUcffdTv59u3b7dGjhxZfT169Ghr/fr1tsbYFP/5z3+siRMn\nWpZlWRkZGdbtt99e/bOysjJr2LBh1sGDB63S0lJr5MiRVnZ2tqlQm6Sh8hUXF1vDhg2zSkpKLMuy\nrPvuu89avny5kTiboqGyVZk/f7511VVXWU899ZTd4TXbocp32WWXWdu3b7csy7Lefvtta/PmzbbH\n2ByHKt95551n5eXlWaWlpdawYcOsvLw8E2E2y5w5c6xLLrnEuuqqq/xed0LdUl/ZAr1eqVJf+aoc\nbt1yxLXUV61axR/+8AcA/vCHP/Dll1/6/TwhIYFXXnml+rq8vJyIiAhbY2yKVatWMWjQIAD69evH\n2rVrq3+2adMmEhMTiYqKIiwsjAEDBrBy5UpToTZJQ+ULDw9nwYIFhIeHA4HzO6vSUNkAvv/+e9as\nWcPo0aNNhNdsDZVvy5YtxMTE8NprrzF27Fjy8vJISkoyFWqTHOr3FxYWxoEDBygpKQHA5XLZHmNz\nJSYm1tlr6YS6pb6yBXq9UqW+8kHT6hbb9373Vdc2sp06dSIqKgqAdu3akZ+f7/fzkJAQYmJiAJg5\ncyYnnXQSiYmJ9gTcDPn5+bRv3776OjQ0FI/Hg9vtrvWzdu3acfDgQRNhNllD5XO5XHTs2BGAefPm\nUVRUxFlnnWUq1MPWUNn27dvHCy+8wIsvvsgHH3xgMMqma6h8+/fvJyMjg2nTptG9e3duvfVW+vbt\ny5lnnmkw4sPTUPkAbrzxRkaOHEnbtm0ZNmxYdf0TSIYNG8auXbtqve6EuqW+sgV6vVKlvvI1tW4x\nmtRHjRrFqFGj/F4bP348BQUFABQUFPj9D1mltLSUSZMm0b59+4AZR4mKiqouF+BXqURFRfk9vBQU\nFNChQwfbY2yOhsoH3nHN9PR0tm3bxgsvvGAixCZrqGwfffQRubm53Hzzzezbt4+SkhKOPfZYRowY\nYSrcw9ZQ+WJiYujRo0d163zQoEGsXbs2oJJ6Q+XbvXs3b775JsuXL6dt27bcf//9LFu2jAsuuMBU\nuC3KCXVLQwK5XjmUptYtR1z3+2mnncbnn38OwOeff87vfve7Wvfcfvvt9O7dm+nTpwdMV5lvuTIy\nMujZs2f1z4477ji2bdtGXl4epaWlrFy5kv79+5sKtUkaKh/AQw89RFlZGS+++GJ1d1mgaKhsY8eO\n5d1332Xu3LnccsstXHLJJQGV0KHh8nXv3p3CwsLqyWWrVq3i+OOPNxJnUzVUvpKSEkJCQggPD69u\n+eXl5ZkKtdms3+wl5oS6pcpvywaBXa/81m/L19S6xWhLvS5XX301KSkpXHPNNYSHh/PUU08B3hnv\niYmJVFRU8O2331JWVsbnn3+Oy+ViwoQJ9OvXz3DkDRs2bBhffPFF9dhIWloaS5cupaioiCuuuIJJ\nkyZx4403YlkWV1xxBZ07dzYc8eFpqHx9+vRh4cKFDBgwgLFjx+JyubjuuusYOnSo4agb51C/u0B3\nqPKlpqZy3333AXDqqacyePBgk+EetkOVb8SIEYwePZrIyEh69OhBcnKy4YibrqqR46S6pcpvyxbo\n9cpv1fW7a9LnWHU9/oiIiEjAOeK630VERKRplNRFREQcQkldRETEIZTURUREHEJJXURExCGU1EVE\nRBziiFunLiKtZ9euXVxwwQWccMIJ1ZtduFwuoqOjGTFiRMBtnCMi/pTURYJMly5dWLRokd9rkyZN\nMhSNiLQkJXUR8fPuu+/y+uuv43K56NOnD1OnTqVNmzasWLGC5557Dsuy6N69O4888ggdO3ZkyJAh\n9OvXj3Xr1vHWW29VH7IhIvbTmLpIkNmzZw/JycmMGDGC5ORkXn311eqfbdiwgZdeeom33nqLJUuW\n0KZNG1544QVycnKYNm0as2fPZvHixZx66qk88sgj1e8bPHgwH374oRK6iGFqqYsEmYa631euXMmQ\nIUOqT/K68sormTx5Mqeffjr9+vXjqKOOAuCqq67i5Zdfrn7/KaecYlP0ItIQtdRFpJrH46l1WlRF\nRQWWZfm97vF4qKioqL6OjIy0LUYRqZ+SukiQaegMpzPOOIMVK1ZUHz/69ttvM3DgQE455RRWr15N\nZmYmAP/85z8ZOHCgLfGKSOOp+10kyFQd8ViXE088kVtuuYUxY8ZQUVFBnz59ePjhh2nbti0zZszg\nzjvvpLy8nISEBFJTUw/5eSJiLx29KiIi4hDqfhcREXEIJXURERGHUFIXERFxCCV1ERERh1BSFxER\ncQgldREREYdQUhcREXEIJXURERGH+P+uCNevgWYiqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bed1b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the fitted model (red line) against observed values (blue points)\n",
    "plt.scatter(clean_data.srrs_mn.floor, np.log(clean_data.srrs_mn.activity + 0.1))\n",
    "xvals = np.linspace(-0.1, 1.2)\n",
    "plt.plot(xvals, m0 * xvals + b0, 'r--')\n",
    "plt.title(\"Fitted model\")\n",
    "plt.xlabel(\"Floor\")\n",
    "plt.ylabel(\"log(radon)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The answer is: not terribly badly (the fitted line runs convincingly through the centre of the data, and plausibly describes the trend), but not terribly well, either. The observed points vary widely about the fitted model, implying that the prevailing radon level varies quite widely, and we might expect different gradients if we chose different subsets of the data.\n",
    "\n",
    "The main error in this model fit is due to *bias*, because the pooling approach is an an inaccurate representation of the underlying radon level, taken across all measurements."
   ]
  }
 ],
 "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
}