03-bias_variance_tradeoff.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",
    "# 3. The Bias-Variance Tradeoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "import clean_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modelling distribution of radon in MN\n",
    "\n",
    "Visual inspection of the variation in (log) observed radon levels for the households shows a broad range of values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAFkCAYAAADSRRn0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHXi//H3AEJyU0Ss1fUWkWul7QqZjywz08LtsWlb\nJiqYZbvdNBNJJFTUVlltTbeN1OJRW1gp3by0Vg/J1M1spTb1W5r7a5HwtlyEEPDC7fz+8OFsbuXA\nDDMDn3k9/2Iu55z3BxjenDNnzsdmWZYlAABgBD9vBwAAAC2HYgcAwCAUOwAABqHYAQAwCMUOAIBB\nKHYAAAzi9mLfs2ePkpKSzrtv48aNSkhIsN/Ozc3VnXfeqYSEBG3dutXdkQAAMFaAO1eenZ2t9evX\nKyQkxH7fvn379NZbb9lvl5WVKScnR++8845Onz6tcePGafDgwWrXrp07owEAYCS37rH37NlTWVlZ\n9tsVFRVavny50tPT7fft3btXsbGxCggIUGhoqHr16qUDBw64MxYAAMZya7GPGDFC/v7+kqTGxkbN\nnj1bs2bNUvv27e3Pqa6uVlhYmP12cHCwqqqq3BkLAABjufVQ/Pd99dVXKioq0rx583TmzBn9+9//\nVmZmpq699lpVV1fbn1dTU6Pw8HCH67MsSzabzZ2RAQBoczxS7JZlqV+/ftq4caMk6ciRI5oxY4bS\n0tJUVlam5cuXq7a2VmfOnFFBQYFiYmIcrtNms6m01Hf37KOiwhi/j47fl8cuMX7G77vjj4oKc/wk\neajYL7Rn3blzZyUlJWn8+PGyLEvJyckKDAz0RCwAAIxja8uzu/nqf22Sb//XKvn2+H157BLjZ/y+\nO/6m7rFzgRoAAAxCsQMAYBCKHQAAg1DsAAAYhGIHAMAgFDsAAAah2AEAMAjFDgCAQSh2AAAMQrED\nAGAQih0AAINQ7AAAGIRiBwDAIBQ7AAAGodgBADBIgLcDAEBr09DQoMLCApfW0avXpfL392+hREDT\nUewA8D8KCws07akNCu7QxanlT1aW6M+P367o6JgWTgY4RrEDwI8I7tBFoRHdvB0DaDbeYwcAwCAU\nOwAABqHYAQAwCMUOAIBBKHYAAAxCsQMAYBCKHQAAg1DsAAAYhGIHAMAgFDsAAAah2AEAMAjFDgCA\nQSh2AAAMQrEDAGAQih0AAINQ7AAAGIRiBwDAIBQ7AAAGcXux79mzR0lJSZKk/fv3a8KECZo4caLu\nv/9+lZeXS5Jyc3N15513KiEhQVu3bnV3JAAAjBXgzpVnZ2dr/fr1CgkJkSQtWrRIc+fOVZ8+fbR2\n7Vq98MILmjx5snJycvTOO+/o9OnTGjdunAYPHqx27dq5MxoAAEZy6x57z549lZWVZb+9bNky9enT\nR5JUX1+vwMBA7d27V7GxsQoICFBoaKh69eqlAwcOuDMWAADGcmuxjxgxQv7+/vbbnTt3liT985//\n1GuvvaZJkyapurpaYWFh9ucEBwerqqrKnbEAADCWWw/F/5hNmzZp1apVev755xUREaHQ0FBVV1fb\nH6+pqVF4eHiT1hUVFeb4SQZj/L47fl8eu+T+8VdUhLq8jk6dQt2Wk5+/b4/fEY8W+/r165Wbm6uc\nnBx7effv31/Lly9XbW2tzpw5o4KCAsXExDRpfaWlvrtnHxUVxvh9dPy+PHbJM+MvL692/KQmrMMd\nOfn5++74m/oPjceKvbGxUYsWLVLXrl31yCOPyGazaeDAgZoyZYqSkpI0fvx4WZal5ORkBQYGeioW\nAABGcXuxd+vWTWvWrJEk/eMf//jR54wZM0ZjxoxxdxQAAIzHBWoAADAIxQ4AgEEodgAADEKxAwBg\nEIodAACDUOwAABiEYgcAwCAUOwAABqHYAQAwCMUOAIBBKHYAAAxCsQMAYBCKHQAAg1DsAAAYhGIH\nAMAgFDsAAAah2AEAMAjFDgCAQSh2AAAMQrEDAGAQih0AAINQ7AAAGIRiBwDAIAHeDgAALa2hoUGF\nhQVOL19U9G0LpgE8i2IHYJzCwgJNe2qDgjt0cWr544f3K/LnfVs4FeAZFDsAIwV36KLQiG5OLXuy\nsriF0wCew3vsAAAYhGIHAMAgHIoH0Opc6OS3iopQlZdXX3B5Tn6DL6PYAbQ6nPwGOI9iB9AqcfIb\n4BzeYwcAwCAUOwAABqHYAQAwCMUOAIBBKHYAAAzi9mLfs2ePkpKSJElFRUUaP368EhMTNX/+fPtz\ncnNzdeeddyohIUFbt251dyQAAIzl1mLPzs7W7NmzVVdXJ0nKzMxUcnKyVq9ercbGRuXl5amsrEw5\nOTlau3atsrOztXTpUvvzAQBA87i12Hv27KmsrCz77a+++kpxcXGSpCFDhuiTTz7R3r17FRsbq4CA\nAIWGhqpXr146cOCAO2MBAGAstxb7iBEj5O/vb79tWZb965CQEFVXV6umpkZhYWH2+4ODg1VVVeXO\nWAAAGMujV57z8/vv/xE1NTUKDw9XaGioqqurf3B/U0RFhTl+ksEYv++O3/SxV1SEejuCyzp1CnXb\nz8n0n78jvj5+Rzxa7FdccYXy8/N1zTXXaPv27Ro0aJD69eunZcuWqba2VmfOnFFBQYFiYmKatL7S\nUt/ds4+KCmP8Pjp+Xxi7o0le2oLy8mq3/Jx84ed/Ib48/qb+Q+PRYk9NTdWcOXNUV1en6OhoxcfH\ny2azKSkpSePHj5dlWUpOTlZgYKAnYwEAYAy3F3u3bt20Zs0aSVKvXr2Uk5Pzg+eMGTNGY8aMcXcU\nAACMxwVqAAAwCMUOAIBBKHYAAAxCsQMAYBCKHQAAg1DsAAAYhGIHAMAgFDsAAAah2AEAMAjFDgCA\nQSh2AAAMQrEDAGAQih0AAINQ7AAAGIRiBwDAIBQ7AAAGodgBADAIxQ4AgEEodgAADEKxAwBgEIod\nAACDUOwAABiEYgcAwCAUOwAABqHYAQAwCMUOAIBBKHYAAAwS4O0AAFqnhoYGFRYWuLSOXr0ulb+/\nfwslAtAUFDuAH1VYWKBpT21QcIcuTi1/srJEf378dkVHx7RwMgAXQrED+EnBHbooNKKbt2MAaAbe\nYwcAwCAUOwAABqHYAQAwCMUOAIBBKHYAAAxCsQMAYBCPf9zNsiylp6fr4MGD8vf315NPPil/f3/N\nmjVLfn5+iomJUUZGhqdjAUCrcaGLA1VUhKq8vPqCy3NhIN/m8WL/+OOPderUKb3++uv65JNPtGzZ\nMtXV1Sk5OVlxcXHKyMhQXl6ehg8f7uloANAquHJxIC4MBI8Xe1BQkKqqqmRZlqqqqhQQEKA9e/Yo\nLi5OkjRkyBB98sknFDsAn8bFgeAsh8W+d+9e9e/fv8U2GBsbqzNnzig+Pl7fffedVq5cqc8++8z+\neEhIiKqqqpq0rqiosBbL1RYxft8dvyfGXlER6vI6OnUKdSprS2zb25wdu+T6+F3Zdltg8thagsNi\n/9Of/qSKigqNGjVKo0aNUlRUlEsbzM7O1oABAzR9+nQVFxcrKSlJdXV19sdramoUHh7epHWVljbt\nHwATRUWFMX4fHb+nxu7ofdymrsOZrC2xbW9zduznlvXWtls7X3/tN4XDs+JfeeUVrVy5UrW1tZo8\nebIeeOABvf/+++eVcXOcPHlSoaFn/xsNCwtTfX29rrjiCu3atUuStH37dsXGxjq1bgAAfF2T3mPv\n1q2bRo8erYCAAK1Zs0Y5OTlatmyZUlJSNGLEiGZtcPLkyUpLS9P48ePV0NCglJQUXXnllZo9e7bq\n6uoUHR2t+Ph4pwYDAICvc1jsubm52rBhg0pLSzV69Gi99tpruuSSS1RcXKw77rij2cUeHh6urKys\nH9yfk5PTrPUAAIAfcljsn332mR599FENHDjwvPsvvvhiPm8OAEAr4/A99hkzZmjbtm2SpEOHDmnm\nzJkqKyuTJN16663uTQcAAJrFYbGnpKSoe/fuks7upcfFxWnmzJluDwYAAJrPYbF/9913SkhIkCQF\nBgbq7rvvVkVFhduDAQCA5nNY7O3bt7cfipeknTt3qn379m4NBQAAnOPw5Ln58+fr8ccftx9+/9nP\nfqYlS5a4PRgAAGg+h8Xet29fvfvuu6qoqFC7du3sF5cBAACtj8Ni37dvn1auXKnKykpZlmW//5VX\nXnFrMAAA0HwOiz01NVVjx45VTEyMbDabJzIBQJtmNTaqqOhbp5d3ZVnAYbFfdNFFSkxM9EQWADDC\nqapSLV1bpuAOx5xa/vjh/Yr8ed8WTgVf4bDYr7/+euXk5Oj6669XUFCQ/f6uXbu6NRiAts2VvVYT\n9lhdmU/9ZGVxC6eBL3FY7OvXr5ckvfTSS/b7bDabPvzwQ/elAtDmubLXyh4r4DyHxb5lyxZP5ABg\nIGf3WtljBZzn8AI1lZWVmj17tiZOnKjy8nKlpaXpxIkTnsgGAACayWGxz5kzR/369dN3332n0NBQ\ndenSRSkpKZ7IBgAAmslhsR8+fFhjx46Vn5+fAgMDNX36dP3nP//xRDYAANBMDovd399fVVVV9s+w\nFxYWys/P4WIAAMALHJ48N3XqVCUlJenYsWN6+OGHtXv3bi1atMgT2QAAQDM5LPYhQ4boqquu0t69\ne9XQ0KAFCxaoc+fOnsgGAACayWGxP/vss+fd3r9/vyRpypQp7kkEAACc1qw3y+vq6rRlyxYdP37c\nXXkAAIALHO6x/++e+SOPPKL77rvPbYEAAIDzmn16e01NjY4ePeqOLAAAwEUO99iHDRtm/6ibZVk6\nceIEe+wAALRSDos9JyfH/rXNZlN4eLhCQ0PdGgoAADjHYbHn5+df8PHRo0e3WBgAAOAah8X+97//\nXfn5+brlllvUrl07ffTRR4qMjNRll10mm81GsQMA0Io4LPaSkhKtW7dOnTp1knT2rPjf/e53mj9/\nvtvDAXBeQ0ODCgsLnF6+qOjbFkwDwFOaVOwdOnSw3w4MDFRVVZVbQwFwXWFhgaY9tUHBHbo4tfzx\nw/sV+fO+LZwKgLs5LPahQ4dq0qRJuvXWW2VZlv72t79p1KhRnsgGwEXBHbooNKKbU8uerCxu4TQA\nPMFhsaelpem9995Tfn6+goKCNHXqVA0ePNgT2QAAQDM16QI1Xbp0UUxMjB577DEFBga6OxMAAHCS\nwz32l19+WXl5eSopKdGtt96quXPn6q677tLkyZM9kQ8A0AxWY6PLJz726nWp/P39WygRPM1hsb/z\nzjvKzc3V3XffrU6dOunNN9/UmDFjKHYAaIVOVZVq6doyBXc45tTyJytL9OfHb1d0dEwLJ4OnOCx2\nPz+/8w6/BwUF8Z8cALRirpw0ibbPYbEPHDhQixcv1qlTp5SXl6e1a9dq0KBBnsgGAACayWGxz5w5\nU7m5uerTp4/WrVunG2+8UQkJCS5t9Pnnn9eWLVtUX1+vxMREDRgwQLNmzZKfn59iYmKUkZHh0voB\nAPBVDov9/vvv14svvuhymZ+za9cuffHFF1qzZo1Onjyp7OxsffDBB0pOTlZcXJwyMjKUl5en4cOH\nt8j2AADwJQ4/7nb69GkdO+bcSRg/5uOPP9bll1+uhx9+WA899JCGDRumffv2KS4uTpI0ZMgQ7dy5\ns8W2BwCAL/nJPfZNmzbp17/+tUpKSnTTTTepc+fOCgoKkmVZstls+vDDD53aYEVFhY4ePapVq1bp\n0KFDeuihh9TY2Gh/PCQkpMmXrI2KCnMqgykYv++Ovyljr6hgemU4p1On0Fb9+mrN2VqDnyz2Z555\nRrfccosqKyu1ZcsWe6G7qmPHjoqOjlZAQIB69+6toKAgFRf/99KVNTU1Cg8Pb9K6Skt995r1UVFh\njN9Hx9/UsZeXV3sgDUxUXl7dal9fvv7ab4qfLPZf/epX6tevnyzL0s0332y//1zB79+/36lgsbGx\nysnJ0aRJk1RcXKxTp05p0KBB2rVrlwYOHKjt27dz1j0AAE76yWLPzMxUZmamHnroIa1YsaLFNjh0\n6FB99tlnuuuuu2RZlubNm6du3bpp9uzZqqurU3R0tOLj41tsewAA+BKHZ8W3ZKmfk5KS8oP7cnJy\nWnw7AAD4miZNAgMAANoGih0AAINQ7AAAGIRiBwDAIBQ7AAAGodgBADAIxQ4AgEEodgAADEKxAwBg\nEIodAACDUOwAABiEYgcAwCAUOwAABqHYAQAwCMUOAIBBKHYAAAxCsQMAYBCKHQAAg1DsAAAYhGIH\nAMAgFDsAAAah2AEAMAjFDgCAQSh2AAAMQrEDAGAQih0AAINQ7AAAGIRiBwDAIBQ7AAAGodgBADAI\nxQ4AgEEodgAADEKxAwBgEIodAACDUOwAABjEa8V+/PhxDR06VAcPHlRRUZHGjx+vxMREzZ8/31uR\nAABo87xS7PX19crIyNBFF10kScrMzFRycrJWr16txsZG5eXleSMWAABtnleKffHixRo3bpy6dOki\ny7K0b98+xcXFSZKGDBminTt3eiMWAABtnseL/e2331ZkZKQGDx4sy7IkSY2NjfbHQ0JCVFVV5elY\nAAAYIcDTG3z77bdls9m0Y8cOHThwQKmpqaqoqLA/XlNTo/Dw8CatKyoqzF0x2wTG77vjb8rYKypC\nPZAEJurUKbRVv75ac7bWwOPFvnr1avvXEydO1Pz587VkyRLl5+frmmuu0fbt2zVo0KAmrau01Hf3\n7KOiwhi/j46/qWMvL6/2QBqYqLy8utW+vnz9td8UHi/2H5Oamqo5c+aorq5O0dHRio+P93YkAADa\nJK8W+yuvvGL/Oicnx4tJAAAwAxeoAQDAIBQ7AAAGodgBADAIxQ4AgEFaxVnxAIDWwWpsVFHRty6t\no1evS+Xv799CidBcFDsAwO5UVamWri1TcIdjTi1/srJEf378dkVHx7RwMjQVxQ4AOE9why4Kjejm\n7RhwEu+xAwBgEIodAACDUOwAABiEYgcAwCAUOwAABqHYAQAwCMUOAIBBKHYAAAxCsQMAYBCKHQAA\ng1DsAAAYhGIHAMAgFDsAAAah2AEAMAjFDgCAQSh2AAAMQrEDAGAQih0AAINQ7AAAGIRiBwDAIBQ7\nAAAGodgBADAIxQ4AgEECvB0AwI9raGhQYWHBD+6vqAhVeXm1w+WLir51RywArRzFDrRShYUFmvbU\nBgV36OLU8scP71fkz/u2cCoArR3FDrRiwR26KDSim1PLnqwsbuE0ANoC3mMHAMAgFDsAAAah2AEA\nMIjH32Ovr6/XE088oSNHjqiurk4PPvigLrvsMs2aNUt+fn6KiYlRRkaGp2MBAGAEjxf7hg0bFBER\noSVLlujEiRMaNWqUfvGLXyg5OVlxcXHKyMhQXl6ehg8f7uloAAC0eR4/FD9y5EhNmzZN0tnP6fr7\n+2vfvn2Ki4uTJA0ZMkQ7d+70dCwAAIzg8WJv3769goODVV1drWnTpmn69OmyLMv+eEhIiKqqqjwd\nCwAAI3jlc+zHjh3TlClTlJiYqNtuu01PPfWU/bGamhqFh4c3aT1RUWHuitgmMH6zx19REertCIBT\nOnUKdevr0/TXvqs8XuxlZWWaPHmy5s6dq0GDBkmS+vbtq/z8fF1zzTXavn27/X5HSkt9d88+KiqM\n8beB8f/UZWGbgkvCoq0qL6922+uzrbz23aGp/9B4vNhXrVqlEydO6LnnnlNWVpZsNpvS09P1hz/8\nQXV1dYqOjlZ8fLynYwFu4cplYbkkLABneLzY09PTlZ6e/oP7c3JyPB0F8AhnLwvLJWEBOIML1AAA\nYBCKHQAAg1DsAAAYhGIHAMAgFDsAAAah2AEAMAjFDgCAQSh2AAAMQrEDAGAQih0AAIN4ZXY3oK1w\nZRIXiYlc4HusxkaXf+979bpU/v7+LZTI91DswAW4MomLxEQu8D2nqkq1dG2Zgjscc2r5k5Ul+vPj\ntys6OqaFk/kOih1wwNlJXCQmcoFvcuU1A9fxHjsAAAah2AEAMAjFDgCAQSh2AAAMQrEDAGAQih0A\nAINQ7AAAGIRiBwDAIBQ7AAAGodgBADAIl5QFALQajiaRqagIVXl59U8+zgQyFDsAoBVxZRIZJpA5\ni2IHALQqTCLjGoodRmM+dcB3MBf8WRQ7jMZ86oDvYC74syh2GI/51AHfwWF8Pu4GAIBRKHYAAAxC\nsQMAYBCKHQAAg3DynIdkLl+lb46dcXr5a/t21u8njW/BRAAAE1HsnuLfXlaE8x+hqGssacEwnuXq\nZ8lN+FwpgNbP1c/BNzQ0SLLJ39/5g+Et8feOYofbufJZclM+Vwqg9XP1c/DHD+9X+7BIp6+b0VJ/\n71pNsVuWpXnz5unAgQMKDAzUwoUL1b17d2/HQgvhs6UA2gJXr3vRGv7WtZpiz8vLU21trdasWaM9\ne/YoMzNTzz33nLdjGcHVQ+GS9w6H/9ShMUczPJ3DJWEB+JpWU+yff/65brjhBknS1VdfrS+//NLL\niczh6mVVvXk4vCUOjXFJWAC+pNUUe3V1tcLCwuy3AwIC1NjYKD8/Mz6R11hXo8bj/+f08hV+Dfr3\nv//ff283cY9Vapm9VlfWUVT0rU5WOnfy36mqcrUPi3R625Kc3va57Us2ryzvzW17e/m2nN3by7fl\n7K4u35azS679rfo+m2VZVousyUV//OMf9ctf/lLx8fGSpKFDh2rr1q3eDQUAQBvTanaHBwwYoG3b\ntkmSdu/ercsvv9zLiQAAaHtazR7798+Kl6TMzEz17t3by6kAAGhbWk2xAwAA17WaQ/EAAMB1FDsA\nAAah2AEAMAjFDgCAQdpssZ86dUoPP/ywEhMTdd9996mkpO3OfuaM6upqPfjgg0pKSlJCQoJ2797t\n7Uget3nzZs2YMcPbMTzGsixlZGQoISFBEydO1KFDh7wdyeP27NmjpKQkb8fwuPr6es2cOVMTJkzQ\n3XffrS1btng7kkc1NjbqiSee0Lhx4zRhwgR988033o7kFcePH9fQoUN18ODBCz6vzRZ7bm6urrrq\nKq1evVq/+c1v9MILL3g7kke99NJLuu6665STk6PMzEwtWLDA25E8auHChVq2bJm3Y3jU9+dTmDFj\nhjIzM70dyaOys7M1e/Zs1dXVeTuKx23YsEERERF69dVX9cILL+jJJ5/0diSP2rJli2w2m15//XVN\nmzZNTz/9tLcjeVx9fb0yMjJ00UUXOXxuq7mkbHPdc889OvdJvaNHj6pDhw5eTuRZ9957rwIDAyWd\n/YEHBQV5OZFnDRgwQCNGjNDatWu9HcVjfH0+hZ49eyorK0szZ870dhSPGzlypP2qnI2NjQoIaLN/\nup0yfPhwDRs2TJJ05MgRn/t7L0mLFy/WuHHjtGrVKofPbRO/HW+++aZefvnl8+7LzMzUVVddpUmT\nJulf//qXXnzxRS+lc78Ljb+0tFQzZ85Uenq6l9K510+NfeTIkdq1a5eXUnmH6fMpODJixAgdOXLE\n2zG8on379pLO/g5MmzZN06dP93Iiz/Pz81NaWpo2b96sZ555xttxPOrtt99WZGSkBg8erJUrVzp8\nvhEXqCkoKNADDzygzZs3ezuKRx04cEApKSlKTU3V9ddf7+04Hrdr1y6tXbtWS5cu9XYUj2A+hbN7\nazNmzNCaNWu8HcXjjh07pilTpigxMVF33HGHt+N4zfHjxzVmzBht2rSpSYelTZCYmCib7ezkMl9/\n/bV69+6tFStWKDLyxyfIahN77D9m1apVuvjiizV69GgFBwd7Za5wb/rmm2/02GOPafny5erTp4+3\n48ADBgwYoI8++kjx8fE+PZ+CAfsizVZWVqbJkydr7ty5GjRokLfjeNy6detUXFysBx54QEFBQfLz\n8/OZI1WStHr1avvXSUlJWrBgwU+WutSGi/2uu+5Samqq3nrrLVmW5XMnEj399NOqra3VwoULZVmW\nwsPDlZWV5e1YcKMRI0Zox44dSkhIkCSf+50/59yeiy9ZtWqVTpw4oeeee05ZWVmy2WzKzs62n2dj\nuvj4eM2aNUuJiYmqr69Xenq6z4z9fzXl99+IQ/EAAOAs3zmWAQCAD6DYAQAwCMUOAIBBKHYAAAxC\nsQMAYBCKHQAAg1DsgGF27drl9AxoxcXFSktLc3rbH3zwgdPL//Wvf/W5K+kB7kCxAwZy9iIuCxcu\n1O9///sWTtM0EyZM0IoVK3xy9jagJbXZK88BuLDCwkLNmTNHlZWVCg4OVnp6uvr166fi4mKlpKTo\nxIkTiomJUX5+vrZt26aioiKVlpaqd+/ekqRhw4bp6quv1tdff61XX31VL7/8sj799FNVVlYqIiJC\nzz77rCIjI7VhwwatWLFCISEh6tGjh32mwd27d2vRokWqra1VRESEFixYoO7duyspKUn9+/fX559/\nroqKCs2ePVs33HCD2rVrp9jYWL377rs+fS10wFXssQMGsixLjz/+uO655x5t2LBBaWlpmjZtmurq\n6rRw4ULddtttWr9+veLj41VSUiJJ+uijjxQbG3veem688Ua99957qq6u1sGDB7V27Vq9//776tGj\nhzZu3KiSkhItWbJEq1ev1htvvKHTp09Lkurq6pScnKyMjAytW7dOY8eOPW9Gsvr6eq1Zs0azZs3S\nsmXL7PfHxcVpy5YtHvgOAeai2AEDnTx5UkVFRRo+fLiks/O3d+zYUQUFBdqxY4duv/12SWfnuQ4P\nD5ckffvtt7rkkkvOW0///v0lST169FBqaqpyc3O1ePFi7d69WzU1Nfriiy80YMAARUZGymaz6be/\n/a2ks0cLOnbsqCuvvFLS2Wt9Hzp0SNXV1ZJkn1c+JiZGJ06csG+va9euKiwsdNN3BfANFDtgoB+b\nAqKxsVENDQ32edz/l5+f3w9mSTw3LeZXX32l++67T5ZlKT4+XsOHD5dlWbLZbOet69zyjY2NP8hg\nWZb9uecO19tstvOe165dO5+atQtwB15BgIFCQkLUvXt3bd68WdLZ97vLysp0+eWX67rrrtPGjRsl\nSdu2bbPvMXfv3l1Hjx790fXl5+fr2muv1dixY3XppZdqx44damxsVGxsrHbv3q3i4mJZlqV3331X\nktS7d29VVlbqyy+/lCRt2rRJXbt2tR8d+L7vF/vhw4fVs2fPlvtGAD6Ik+cAA9lsNj311FOaO3eu\nnnnmGQUpz7uKAAABGElEQVQFBSkrK0sBAQFKS0tTamqq3njjDfXp08detjfddJNSUlLOW8c5I0eO\n1NSpU3XHHXcoIiJCN954ow4fPqzIyEjNmzdP9957r9q3b68+ffrIsiwFBgbq6aef1oIFC3Tq1Cl1\n7NhRy5cv/8F6//f2p59+qmHDhrnzWwMYj2lbAR+Tk5Oj6667TtHR0dq3b5/mzJmjt956S5L06KOP\naurUqYqJifF4rtraWo0fP16vv/662rVr5/HtA6ag2AEfs337di1dulR+fn4KCgpSRkaG+vbtK0k6\nevSo/vKXvygzM9PjuV566SX16NFDN998s8e3DZiEYgcAwCCcPAcAgEEodgAADEKxAwBgEIodAACD\nUOwAABjk/wO/+DpXTJQZmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bf98358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = clean_data.srrs_mn.activity.apply(lambda x: np.log(x + 0.1)).hist(bins=25)\n",
    "fig.set_xlabel('log(radon)')\n",
    "fig.set_ylabel('frequency');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We aim to determine which contributions of the prevailing radon level, and the floor at which radon level is measured,  produce this distribution of observed values, by modelling the relationship between those quantities and the measured radon level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conventional approaches\n",
    "\n",
    "Two conventional alternatives to modelling, *pooling* and *not pooling* represent two extremes of a [tradeoff between *variance* and *bias*](http://scott.fortmann-roe.com/docs/BiasVariance.html).\n",
    "\n",
    "* *variance*: the difference between observations and the mean behaviour of the model\n",
    "* *bias*: the difference between the mean behaviour of the model, and the mean behaviour of the \"true\" relationship\n",
    "\n",
    "### The bias-variance tradeoff\n",
    "\n",
    "Where the variable we are trying to predict is $Y$, as a function of covariates $X$, we assume a relationship $Y = f(X) + \\epsilon$ where the error term $\\epsilon$ is distributed normally with mean zero: $\\epsilon \\sim N(0, \\sigma_{\\epsilon})$.\n",
    "\n",
    "We estimate a *model* $\\hat{f}(X)$ of $f(X)$ using some technique. This gives us squared prediction error: $\\textrm{Err}(x) = E[(Y − \\hat{f}(x))^2]$. That squared error can be decomposed into:\n",
    "\n",
    "$$\\textrm{Err}(x)=(E[\\hat{f} (x)] − f(x))^2 + E[(\\hat{f}(x) − E[\\hat{f}(x)])^2] + \\sigma^2_e$$\n",
    "\n",
    "where\n",
    "\n",
    "* $E[\\hat{f} (x)] − f(x))^2$ is the square of the difference between the model $\\hat{f}(x)$ and the 'true' relationship $f(x)$, i.e. the square of the *bias*\n",
    "* $E[(\\hat{f}(x) − E[\\hat{f}(x)])^2]$ is the square of the difference between the mean behaviour of the model and the observed behaviour of this model, i.e. the square of the *variance*\n",
    "* $\\sigma^2_e$ is the noise of the 'true' relationship that cannot be captured in any model, i.e. the *irreducible error*\n",
    "\n",
    "With a known true model, and an infinite amount of data, it is in principle possible to reduce both *bias* and *variance* to zero. In reality, both sources of error exist, and we choose to minimise *bias* and/or *variance*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The trade-off in the radon model\n",
    "\n",
    "Taking $y = \\log(\\textrm{radon})$, floor measurements (basement or ground) as $x$, where $i$ indicates the house, and $j[i]$ is the county to which a house 'belongs'. Then $\\alpha$ is the radon level across all counties, and $\\alpha_{j[i]}$ is the radon level in a single county; $\\beta$ is the influence of the choice of floor at which measurement is made; and $\\epsilon$ is some other error (measurement error, temporal variation in a house, or variation among houses).\n",
    "\n",
    "We can take two contrasting, and extreme approaches:\n",
    "\n",
    "* **Complete pooling** - treat all counties the same, and estimate a single radon level: $y_i = \\alpha + \\beta x_i + \\epsilon_i$\n",
    "* **No pooling** - treat each county independently: $y_i = \\alpha_{j[i]} + \\beta x_i + \\epsilon_i$\n",
    "\n",
    "When we do not pool, we will likely obtain quite different parameter estimates $\\alpha_{j[i]}$ for each county - especially when there are few observations in a county. As new data is gathered, these estimates are likely to change radically. This is therefore a model with *high variance*.\n",
    "\n",
    "Alternatively, by pooling all counties, we will obtain a single estimate for $\\alpha$, but this value may deviate quite far from the true situation in some or all counties. This is therefore a model with *high bias*.\n",
    "\n",
    "So, if we treat all counties as the same, we have a biased estimate, but if we treat them as individuals, we have high variance - the *bias-variance tradeoff*. It may be the case that neither extreme produces a good model for the real behaviour: models that minimise bias to produce a high variance error are *overfit*; those that minimise variance to produce a strong bias error are *underfit*."
   ]
  }
 ],
 "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
}