Simple Linear Regression Coding Exercises in R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I. Optimize Objective Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that data comes from a linear model $Y_i = \\beta_0+\\beta_1 X_i$, simple linear regression using least square method aims to find the line that minimizes the following objective function:\n",
    "$$Q(b_0,b_1) = \\sum_i (y_i-\\hat y_i)^2 = \\sum_i (y_i-b_0-b_1 \\hat y_i)^2$$\n",
    "We could solve this using the R builtin linear regression function, but we can also solve it by optimizing the objective function, which in this case is the least squared error function. We could imagine defining any other objective function (e.g. sum of absolute error or sum of Euclidean distance) we desire and solve for the line of best fit using this approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulate some data\n",
    "n <- 100\n",
    "x <- rexp(n,rate=2)\n",
    "b0 <- 10\n",
    "b1 <- 5\n",
    "error <- rnorm(n, sd=3)\n",
    "y <- b0+b1*x+error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>10.9985664289718</li>\n",
       "\t<li>3.56114149067457</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 10.9985664289718\n",
       "\\item 3.56114149067457\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 10.9985664289718\n",
       "2. 3.56114149067457\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 10.998566  3.561141"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<dl class=dl-horizontal>\n",
       "\t<dt>(Intercept)</dt>\n",
       "\t\t<dd>10.998271484167</dd>\n",
       "\t<dt>x</dt>\n",
       "\t\t<dd>3.5614951455151</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[(Intercept)] 10.998271484167\n",
       "\\item[x] 3.5614951455151\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "(Intercept)\n",
       ":   10.998271484167x\n",
       ":   3.5614951455151\n",
       "\n"
      ],
      "text/plain": [
       "(Intercept)           x \n",
       "  10.998271    3.561495 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# objective function\n",
    "least_square_objective <- function(params, x, y){\n",
    "    b0_hat <- params[1]\n",
    "    b1_hat <- params[2]\n",
    "    SSE <- sum((y-b0_hat-b1_hat*x)^2)\n",
    "}\n",
    "# optimize\n",
    "opt <- optim(c(10,5),least_square_objective,x=x,y=y)\n",
    "opt$par\n",
    "# compare with R's built-in linear regression function\n",
    "slr <- lm(y~x)\n",
    "slr$coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could change the objective function to anything according to our needs. For example, suppose we are interested in optimizing the least absolute value:\n",
    "$\\min_{b_0,b_1} \\sum_i |y_i-\\hat y_i|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>10.6692065697048</li>\n",
       "\t<li>3.88911602401723</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 10.6692065697048\n",
       "\\item 3.88911602401723\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 10.6692065697048\n",
       "2. 3.88911602401723\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 10.669207  3.889116"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# objective function\n",
    "least_absolute_obj <- function(params, x, y){\n",
    "    b0_hat <- params[1]\n",
    "    b1_hat <- params[2]\n",
    "    sum(abs(y-b0_hat-b1_hat*x))\n",
    "}\n",
    "#optimize\n",
    "opt2 <- optim(c(10,5),least_absolute_obj,x=x,y=y)\n",
    "opt2$par"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we could find $\\hat \\beta_0$, $\\hat \\beta_1$ using MLE, in which case the objective function will be the log likelihood function $\\sum_i \\log f(y_i | \\sigma)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>-0.130775175507615</li>\n",
       "\t<li>0.979382413011618</li>\n",
       "\t<li>1.2904521961941</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item -0.130775175507615\n",
       "\\item 0.979382413011618\n",
       "\\item 1.2904521961941\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. -0.130775175507615\n",
       "2. 0.979382413011618\n",
       "3. 1.2904521961941\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] -0.1307752  0.9793824  1.2904522"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<dl class=dl-horizontal>\n",
       "\t<dt>(Intercept)</dt>\n",
       "\t\t<dd>-0.130963338886415</dd>\n",
       "\t<dt>x</dt>\n",
       "\t\t<dd>0.979730512306197</dd>\n",
       "</dl>\n"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[(Intercept)] -0.130963338886415\n",
       "\\item[x] 0.979730512306197\n",
       "\\end{description*}\n"
      ],
      "text/markdown": [
       "(Intercept)\n",
       ":   -0.130963338886415x\n",
       ":   0.979730512306197\n",
       "\n"
      ],
      "text/plain": [
       "(Intercept)           x \n",
       " -0.1309633   0.9797305 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n <- 100\n",
    "x <- runif(n, min=0,max=1)\n",
    "beta0 <- 0\n",
    "beta1 <- 1\n",
    "error <- rnorm(n)\n",
    "y <- beta0+beta1*x+error\n",
    "\n",
    "mle_obj <- function(params,x,y){\n",
    "    b0 <- params[1]\n",
    "    b1 <- params[2]\n",
    "    sigma2 <- params[3]\n",
    "    log_likelihood <- sum(dnorm(y,mean=b0+b1*x, sd=sqrt(sigma2),log=TRUE))\n",
    "    return(-log_likelihood)\n",
    "}\n",
    "\n",
    "optim(c(0,0,1),mle_obj, x=x,y=y)$par\n",
    "\n",
    "lm(y~x)$coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. Simulate data\n",
    "Two basic ways of simulating data:\n",
    "1. create a data generation model\n",
    "2. draw from an existing sample or dataset (bootstrap)\n",
    "  \n",
    "The example below simulates an experiment in which a dice is thrown 100 times and the mean of the number is recorded. This process is then repeated 1000 times, allowing us to plot the sampling distribution of the sample mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAAAAAP9NTU1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD////xw1/KAAAACXBIWXMA\nABJ0AAASdAHeZh94AAAgAElEQVR4nO2di7aiOBAAszxERWX8/59dAfUGRMCQQDepOmfuoAEa\nmi55RTR3AFiM2XoBAPYAIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyA\nSAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8gEgAHkAkAA8gEoAHEAnAA4gE4AFEAvAA\nIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQCcAD\niATgAUQC8MCuRCoPiTFpcXWc3Bhj/Td3koasuHXm0ePgFOgwMsd5HB/5sGP/MEc/ufzSZEyb\nrlsz7BhCFvtYi5biVdZHt+kXiPSgtOfR4ZKYj6mmA72mWlBpx3q5/kT6ZY6ecvmlyZhzM3hG\nJIGc/or64jSDZSK1QYcm/nhvVqBX64JKS42x9yg/zNFXLr80PXbhzWCGSAJ5lE1RPQ4XHlsn\nd5qBwzZ9TlKVj+DJr/OdJ9ICerP4QaSQuTR/9hhEEshrk1SvgXNu6sP826vxlJr08fl6Skx2\nuVtvmsPNnoH1X/kopMPzQ/12eHyOnnoV8vcqbQ/uXmYd6w/b/Hx/77TaoVtqik6ER/ikG779\nvzNVG6E81Adq5V9ge+lefIxkL6o1x0v2Xu17VSQmeZ3k+c3ll3nnbbLKZmhgxH60oTUVxo5E\nemzF3D4OyeyDk/poonlxKzrHYXnzKnlvss5/zxOFZgte2uHsm0hlezbSvnFLnqGzXgGnr7fu\n9jJ1jwq/iPRanfw1kr10/XXO7/YsXlP8zbGwVvu1sHbufOVyaN6Pl6c2WYdm6HPEfrShNZXG\njkRqjuuTw/n2fplVzVmzVVWPEaxyfL9p0ufL7n9PmtP15O+1FfPvVdXOpH3j0JxNV4+COHUL\n2DQNHxESe15mUKT8PXI+sHRPOiP1lra/HB9rZh+a+srl0LwfL095805i8tPAiF+imd61T1ns\nSKR3GaXNkU36vMb6Lp7Hxqk3W3q9v7ZeveEeH3pl0jkssxrLxon65bl9dU7MF5H+4rR/69gd\nuV7L8BGhDn++90Tqvyzr8qseh4zmvazW0j0ZGMnOjzXHerULu2yrw+tKms9cDs770XQ91Xub\nx07+dB0Y8SPa55rKQ/TC/Up9yt9QWG++N8el89+zqdnAncMy67+6Ip5nCc+j+uaC7cfMu3Hu\n7efr60SlU8Cda+Sv1+X7o7fXar2sj4Kal8V7Wa2lezIwkp2e3nL8rdnT7s5lBS+5HJx3LdKt\nvq7+8P12HR/x25rKQ/TC/c7tfGgOEk7PV0VzuH0fkMT6rx5Ivza+1bjfO9P0X9lxjs9jke69\nJfOsloHw5t4P2Hv5nvY2OPJ7PoMj9Rb2Y85PepcdPeRycN61SI98po89T3K/Do/4Ga2/MtIQ\nvXBu3PJ2W55fn6lzNv7XcXpV/lWk7mHc635mc+b9UQpuIt1HRu4vjZNIA7WwMJeD825Eeuw8\nm33Xh0jmSzREWovk9XH8zHl9EJYeTlczuvHtXcSYSJN7pO7h4b06t5eeso85TYWf3iMlQ/O8\nj4zUW9gvaxYgl0PzbkU6N1cwzy+ROiMORkOktTi8L+u0BzbPGztjkpjncUv3JGWwkCfPkbLP\n8/vy8KUU/iK8Tiue4etSvHyLP3j601mcbyP1FnZ4zULkcmjerUhtL7vb/Tow4mA0RFqL+pJV\nczewbHtqPhM/8Slal3I9wenbOM//X1ftukcp75Eueecadvo+d37tGCpr5IHwdY+2pDmxvyRW\n/L+pvl2Q6yzOjKt21ce05/Zy2/nVbcdnLofmfW8EqlNUHzTeBkYcjIZIq/F3E6XZUllTl2Uy\nsfFbksHGu/X/1/tIb6xP0UchZrfmmkPxnLS4D4v0DF/X96EboTfV+yZlu68YLq8vIz3pL0d/\nzay7pr5yOTTvtqnoLow94mA0RFqPdxnV+45XV4T6ruHl+8a3b/GPilS2c/vs2dDx6DX262JD\n80F8MN1zJStCbk3bHu00SzQ01Wvtul+F6JXX8EhP+nPsrlnnOrevXA7N27x3seXwiIPREGlF\nmu/QmPzYnvVeD/XN+fqWRf5949dXiJJi6Kp0fwPWs8vK3vZ8lVv2nMW7uTk/ytozlvrzvXPK\nYkU4WX3emgjn92j9qZq163ej65fX4Egv8uHzq6qoewT1Tmb85HJo3m1TfV+o+jLiUDREEozD\ntqlGenlHjfA6D07Ua//DxjftFbar8/cK9g4iRcwPG//vNH7ggi4gUtRr/8PGf38xontODi8Q\nKWJ+2fjVsb7ClhzYHw2DSACwGEQC8AAiAXgAkQA8gEgAHkAkAA8gEoAHEAnAA4gE4AFEAvAA\nIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQCcAD\niATgAUQC8AAiAXgAkQA8gEgAHnAX6XJsfwA1Ly7TIwPsG1eRqvTvl7fs334HiBJXkQqTnK/N\n0K1M+PEtiB1XkRJzfQ9f+X1iiB1XkTq/zxb5j7UBsEdSzL/7fwumNvAdh3Q6bobHOVJ5a4Y4\nR9qKhSJ5W479saJI98zyN61c5wILQKRQrCnS/VI095GS/Mh9pG1ApFCsKhJsDSKFApGiApFC\nsYlIXPveCkQKBSJFBSKFYkWRFl93h8UgUihWFOmSINLWIFIo1jy0q3KTNXdkBy3CshVApFCs\ne450NuZ8n3GOxBYLBCKFYuWLDbfM5BUibQYihWL1q3ZHk5SItBWIFIr1L39f0+lzILZYIBAp\nFFvcRzog0lYgUihkdhFiiwUCkUKBSFGBSKFApKhApFAgUlR4F2n5F653gsy+dvFuj8D4F+lf\nh8ENVy4IqYYVRToh0uZsIVIaxdZc89Dumsx9vmoUqd+CLUSK43hv1XOk69xnB0WR+i34Z/5b\n8EgpRPrOuhcbTtaj7QKFgDHG90hBRHoaakyVmvxpVfv3lJrkZI9Y9x871s9te37g/o1Q5sYk\nRTvWLW/GEgZX7aJiS5Hy2g9LpPaXSTJrxGP9Rtk8ta3ojHBsd5hFM1bztTZxJiFSVGwg0tud\nrLJePHYy9RtVZsq/8R5vnJ5/k84Izbdvzs20z/Z0YlFXB5GiYkuRLvaLen9Ti1U9Dvfe412a\nv7fvI1hjTSzq6iBSVGwpUu/FxxUOe7yPEW7lMXuKZM1NEIgUFVpFyt5DiCQsRJyIEmlkvN4I\nB5OeyhsiSQwRJzJEurSnQOX38XojNG8iksgQceJfpOlbun+XD+51f6FTfSXO1Bfhkmt9Z7F7\nLcH6a41QX2G4co4kMkSceBdpBqmpL2Y/S7/pcJlb5z3J7T33nkjWCMVT0wsiiQsRJ1uIdEkt\nke7HxBz+ejaYw9ujT5GsEQ7GZJey3Tn9tUsCkaJiC5HiAJGiApFCgUhRgUihQKSoQKRQIFJU\nIFIoECkqECkUiBQViBQKRIoKRAoFIkUFIoUCkaICkUKBSFGxRafV5VQHM/f5U31W60uESFHh\nX6T/OoTZcLn7004QCUKgU6T2ixiOk/pckLFAq0wiMEScaBVpi0l/DLTKJAJDxMkWIrk++PEx\nVnq6v07DXu/Wj73Lyv603QiP18Vf032VR1EiUlRsI5LTgx+z11gdkZ6/xHDqTduN8Hyd3e/r\nPYoSkaJiG5FcHvzYfNH8mtQj2MdnSf3Q63M92ue07wjmb9rVHkWJSFGxjUguD35sH31S/u1W\nXuN2H5nSmfYVwTynzb8G9P8oSkSKio3OkTp/5z340X7qkFXMj7Oa/Pr8JYbBaa0Heb0Gf3qC\nnuujKBEpKqSJ9P3Bj19Eqp/58HxiyvC0S0VyfRQlIkWFDJHeTSMPfvwm0uOQrUjrE5cv0w6K\n5GWJxkGkqJAg0rwHP+Z/5znDJnyZtm26NNMePgIuWaJxECkqJIg078GPX67ape1VtfTrtK1I\n7bTlR8AlS/RzbgJMIjBEnGzzpNXe33kPfnzfR+q8e35P8GXaVqRm4vwz4KIl+jU3ASYRGCJO\ntvgaxUfZznzw4ylpezZ03216Nly+T/s86suf0/YDLlqi0bWcN9rCSQSGiJM4vo80t/i9xlxl\nEoEh4gSRgsVcZRKBIeIEkYLFXGUSgSHiBJGCxVxlEoEh4iQOkbYAkaICkUKBSFGBSKFApKhA\npFAgUlQgUigQKSoQKRSIFBWIFApEigpECgUiRQUihQKRogKRQoFIUYFIoUCkqECkUCBSVCBS\nKBApKhaKBN/5fWMgkl6WiQRecU/35dg+rTwvLqFCwCiIJAjXdFeptSPMRkdliwYCkQThmu7C\nJOf2qc23Mhn/TVC2aCAQSRCu6W5+kePJtf41Df8hYAJEEoRrujsXNsavcrBFA4FIgmCPpBdE\nEsSCc6Syfdol50hbgUiCcE53Zl21S6sgIWAcRBLEgvtIRXMfKcmP3EfaBkQSBD0b9IJIgkAk\nvSCSIJzTXdU/mfH8bTQuf28CIgnCuYtQ0na0a2eCSFuASIJwv/x9eth0Svo/wuYxBEyASIJw\nvyHb/HdL0hsibQQiCWJpF6Eqy4ZEWvgtKZgDIgnCNd2ped2ETTP2SNuASIJwTffJHJ5DN5Mh\n0iYgkiCc01287Sknjt7YooFAJEG4p/uav4ZuB0TaAkQSBD0b9IJIgkAkvSCSIBBJL4gkCETS\nCyIJApH0gkiCcO/ZMLvzAls0EIgkCPcbsoi0NYgkCOd0X5Px56t6CAHjIJIgFtyQHX92kI8Q\nMAoiCWJBuk/Wo+0ChYAxEEkQXLXTCyIJApH0gkiCQCS9IJIgEEkviCQIRNILIgkCkfSCSIJA\nJL0gkiAQSS+IJAhE0gsiCQKR9IJIgkAkvSCSIBBJL4gkCETSCyIJApH0gkiCQCS9IJIgEEkv\niCQIRNILIgkCkfSCSIJAJL0gkiAQSS+IJAhE0gsiCQKR9IJIgkAkvSCSIBBJL4gkCETSCyIJ\nApH0gkiCQCS9IJIgEEkviCQIRNILIgkCkfSCSIJAJL0gkiAQSS+IJAhE0gsiCQKR9IJIgkAk\nvSCSIBBJL4gkCETSCyIJApH0gkiCQCS9IJIgEEkvEyKNstpCxgIi6WXBHolt4htE0gsiCQKR\n9IJIgkAkvSCSIBBJL4gkCETSCyIJApH0gkiCQCS9IJIgEEkviCQIRNILIgkCkfSCSIJAJL0g\nkiAQSS+IJAhE0gsiCcI9o5dj3nTIz4tLqBAwCiIJwjWjVWp9uSULEgImQCRBuGa0MMn52gzd\nysQUIULAxLfzEEkQrhlNzPU9fDVJiBAwkTtEEoRrRjtfVh7/5jIbzR1E0gJ7JNEgkhYWnCOV\nt2aIc6SAIJIWnDOaWWe9aRUkBCCSGhbcRyqa+0hJfuQ+UjAQSQv0bBANImkBkUSDSFqgi5Bo\nEEkLdBESDSJpgS5CokEkLXBDVjSIpIVAXYT45QM/IJIW2COJBpG0QBch0SCSFugiJBpE0gJd\nhESDSFqgZ4NoEEkLiCQaRNICIokGkbSASKJBJC0gkmgQSQvuPRtmd15go7mDSFpwzegJkdYA\nkbTgnNFrMv7lCQ8hAJHU4J7R63jHIB8hAJG0sCCjJ6vfaqAQ0YNIWuCqnWgQSQuIJBpE0gIi\niQaRtIBIokEkLSCSaBBJC4gkGkTSAiKJBpG0gEiiQSQtIJJoEEkLiCQaRNICIokGkbSASKJB\nJC0gkmgQSQuIJBpE0gIibY0ZZWxKRBIEIm2Ne8EjkiAQaWsQaRcg0tYg0i5ApK1BpF2ASFuD\nSLsAkbYGkXYBIm0NIu0CRNoaRNoFiLQ1iLQLEGlrEGkXINLWINIuQKStQaRdgEhbg0i7AJG2\nBpF2ASJtDSLtAkTaGkTaBYi0NYi0CxBpaxBpFyDS1iDSLkCkrUGkXYBIW4NIuwCRtgaRdgEi\nbQ0i7QJE2hpE2gWItDWItAsQaWsQaRcg0tYg0i5ApK1BpF2ASFuDSLsAkbYGkXYBIm0NIu0C\nRNoaRNoFiLQ1iLQLEGlrEGkXINLWINIuQKStQaRdgEhbg0i7AJG2BpF2ASJtDSLtAkTaGkTa\nBYi0NYi0CxBpaxBpFyDS1iDSLnDP6OWYm5q8uIQKEQWItAtcM1ql5o8sSIhIQKRd4JrRwiTn\nazN0KxNThAgRCYi0C1wzmpjre/hqkhAhIgGRdoFrRo359sJbiEhApF3AHmlrEGkXLDhHKm/N\nEOdIk5hRRqcca0QkQThnNLNKIa2ChNgN7jWNSFpYcB+paO4jJfmR+0gTINL+oWfDCiDS/gmU\n0dknAVGASPvHzmh6vM2f8HYwyfF+P6UmGb3UwEZDpBjo3g4ys12qknpfczrSRWgGiLR/7IxW\n58Nsl4r6kneRmEN1rwouf4+DSPunn9HLMZ3lUtJMaExz4ZsbsuMg0v4ZyOi1Pmo7TU1n/v7S\nRWgCRNo/nxktsxnnPe89Uv23Yo80DiLtn15Gq+Njd5SW1cOmfHS61zlSUT2HZ4eIEETaP52M\nXuqLDUXbG3Xi9g9X7X4AkfZP5z7SY2d0enWbGz9c4z7SLyDS/uncR8rL0CHiBJH2T+c+UvgQ\ncYJI+6eT0aqoj+eSwq9RbDRE2j92Rm9Jc4XBmOSHPne/hYgTRNo/dkazusNPvV+auPS9IESc\nINL+6XZa7Q94DxEniLR/7Iwmbc+5e4VIfkGk/WNntDBZ/a3xSzbeU2FJiN0S6PkmiKSFTkZf\nDzSZ6Ge3JMReCeQKImmhm9Fz/TyTbKrn96IQOwWRImeFjEax0RApchDJD4gUOYjkB0SKnE5G\nj+mcy0yLQuwVRIocO6PHMI+ii2KjIVLkdG/Ier5e9xlityBS5Ax2EQoXYrcgUuTYGc1NkG8k\nRbHREClyul+jyCZ+WGJxiN2CSJHTPbTjYoMriBQ5iOQHbSI597GFYbgh6wdtIjk3wjCI5AdE\nipxu0sq83q/nXh/ZEMd2QaTI6SQtaw+QefjJ7yBS5NhJO5ms+Zb5yRxChdgtiBQ53S5C1f35\nQK5QIXYLIkVOv4sQIrmBSJFjJy197pGuJg0VYrcgUuQMnCOVnnuBR7FdEClyOknLeYqQK4gU\nOZ/3kUx+DhlipyBS5NCzwQ+IFDmI5AdEihxE8gMiRU73PhJfo3AFkSIHkfywK5H4stLvDOTl\nknn9nTFEUieSc2PEDOWlotPqzyBS5AzmhUO7n0GkyBnKy8kkoUPsDkSKnOGLDcdQIXYLIkXO\nkEip3ycXR5F6RIocbsj6AZEiB5H8gEiR8+WGrM87b1GkHpEiB5H8gEiR08nLMSkffy8JX+z7\nGUSKHDsvR3Nt/r8ar32Eokg9IkVO99CuP+A9xG5BpMix85K890g8RehXECly7LwUpjlH4ilC\nDiBS5HTykj2v1xXhQuwVRIqcbl7OzVOEypAhdgoiRQ49G/yASJGDSH5ApMjp5oUfGnMFkSLn\n82LDnR8acwCRIsfOCz805g4iRU73hiw/NOYKIkVOv4sQIrmBSJFj54UfGnMHkSJn4ByJLkIO\nIFLkdPLCD405g0iR83kfae4PjV2OrXd5cfklxE5BpMhxzUuVWl9KH9+DRZF6RIocOy/5D72+\nC5Oc228v3R7nVKMTRpF6RIqc/uXvuby+BFhzHX/EcRSpR6TI6V/+nj3dfAOjSD0iRY6dlyrP\nJq4b/MEeqQsiRU53xzL/mXb119Lbrq2cI9UgUuS4ivT+Wnrz1P3RQ8IoUo9IkeOel0vR3EdK\n8iP3kRApelbISxSpR6TIeeXF869VR/c72IgUOV2Rfij6qqgv1R1TY7KJLkVRpB6RIsdVpFvy\nGLdK6CL0BJEix1Wkg8mrx5/D7eHUgcvfiBQ7riKZuheEabtCVNyQRaTYcRfpXndvsCeeCrFr\nECly3A/tru8fVKrGT5KiSD0iRc6fSL9dsb6apLje8+RhUpma0aeFR5F6RIocV5HqJzu8Oc4K\nsWsQKXIW5OV8aL4lmx8nnssaReoRKXLoIuQHRIocRPIDIkUOIvkBkSIHkWZjRhmdMkwjIgkC\nkWazRWUikhYQaTaINNkYMYg0G0SabIwYRJoNIk02RgwizQaRJhsjBpFmg0iTjRGDSLNBpMnG\niEGk2SDSZGPEINJsEGmyMWIQaTaINNkYMYg0G0SabIwYRJoNIk02RgwizQaRJhsjBpFmg0iT\njRGDSLNBpMnGiEGk2SDSZGPEINJsEGmyMWIQaTaINNkYMYg0G0SabIwYRJoNIk02RgwizQaR\nJhsjBpFmg0iTjRGDSLNBpMnGiEGk2SDSZGPEINJsEKlpdH5M5r5BpNkg0rLGfYNIsxFXmYgk\nCESajbjKRCRBINJsxFUmIgkCkWYjrjIRSRCINBtxlYlIgkCk2YirTEQSBCLNRlxlIpIgEGk2\n4ioTkQSBSLMRV5mIJAhEmo24ykQkQSDSbMRVJiIJApFmI64yEUkQiDQbcZWJSIJApNmIq0xE\nEgQizUZcZSKSIBBpNuIqE5EEgUizEVeZiCQIRJqNuMpEJEEg0mzEVSYiCQKRZiOuMhFJEIg0\nG3GViUiCQKTZiKtMRBIEIs1GXGUikiAQaTbiKhORBIFIsxFXmYgkCESajbjKRCRBINJsxFUm\nIgkCkWYjrjIRSRCINBtxlYlIgkCk2YirTEQSBCLNRlxlIpIgEGk24ioTkQSBSLMRV5mIJAhE\nmo24ykQkQbiv+uWYNz8bmheXUCFkIa4yEUkQrqtepdZP8GZBQkhDXGUikiBcV70wyfnaDN3K\nxBQhQkhDXGUikiBcVz0x1/fw1SQhQkhDXGUikiBcV92Yby+8hZCGuMpEJEGwR5qNuMpEJEEs\nOEcqb80Q50hbNSKSIJxXPbOu2qVVkBDCEFeZiCSIBfeRiuY+UpIfuY+ESDMa9w09G2YjrjIR\nSRCBVt3YhAmxOuIqE5EEsXzVJ0XZS3bFVSYiCQKRZiOuMhFJEO43ZGcfve0lu+IqE5EE4brq\nlwSRtm5EJEE4r3qVm6y5I8uh3VaNiCSIBat+NuZ8R6TtGhFJEEtW/ZaZvEKkzRoRSRDLVv1o\nkhKRtmoUKNIoY1OqZ+HaXdPpBO0lgYgUsFE9i9fugEhbNSKSIFZYu70kUFzxIZIgEGk24ooP\nkQSBSLMRV3yIJAhEmo244kMkQSDSbMQVHyIJApFmI674EEkQiGTjfkMRkZY1qgeRbMTVFyJp\nAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUX\nImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE\n1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGR\nbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJp\nAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUX\nImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE\n1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGR\nbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJpAZFsxNUXImkBkWzE1RciaQGRbMTVFyJp\nwX3tLsfc1OTFJVSI1RFXX4ikBde1q1LzRxYkxAaIqy9E0oLr2hUmOV+boVuZmCJEiA0QV1+7\nEmmUsSk14LoCibm+h68mCRFiA+QV31ijNpGcGzXgugKdj5DxzxNFORJXQoikBfZINuJKCJG0\nsOAcqbw1Q5wjbdWISIJwXoHMOlFMqyAh1kdcCSGSFhbcRyqa+0hJfuQ+EiIFbdTACiugKEfi\nSgiRtBBoBZTeIhBXQoikBecVqA7GZOVzJnvJkbgSQiQtuK5AlbQd7dqZ7CVH4koIkbTgfvn7\n9LDplDTd7BAJkUI2asB1BZJ2wluS3hBpo0ZEEoTrCrzcqbJMnEjuvSN1NSKSIFxXIDWvm7Bp\nJk6kSBoRSRCuK3Ayh+fQzWSIhEghGzXgvALF255yyfFSCMQVAiIta9SA+wpc89fQ7SArR+IK\nAZGWNWpghRVAJERa1qgBRNLbiEiCQCS9jYgkiPhE2s9NJkQSRHwi7acRkQSBSHobEUkQiKS3\nEZEEgUh6GxFJEIiktxGRBIFIehsRSRCIpLcRkQSBSHobEUkQiKS3EZEEgUh6GxFJEIiktxGR\nBIFIehsRSRCIpLcRkQSBSHobEUkQiKS3EZEEgUh6GxFJEIiktxGRBIFIehsRSRCIpLcRkQSB\nSHobEUkQiKS3EZEEgUh6GxFJEIiktxGRBIFIehsRSRCIpLcRkQSBSHobEUkQiKS3EZEEgUh6\nGxFJEIiktxGRBIFIehsRSRCIpLcRkQSBSHobEUkQiKS3EZEEgUh6GxFJEIiktxGRBIFIehsR\nSRCIpLcRkQSBSHobEUkQiKS3EZEEgUh6GxFJEIiktxGRBIFIeht3JdIoY1MKAZH0Nu5KJOdG\nISCS3kZEEgQi6W1EJEEgkt7GeERScAKFSHob4xHJuXE9EElvIyJNNq4HIultRKTJxvVAJL2N\niDTZuB6IpLcRkSYb1wOR9DYi0mTjeiCS3kZEmmxcD0TS24hIk43rgUh6GxFpsnE9EElvIyJN\nNq4HIultRKTJxvVAJL2NiDTZuB5KRXLvxrijRkSabFwPrSLRiEgzGtcDkfQ2ItJk43ogkt5G\nRJpsXA+5InEaNNWISJON6yFYJBoRaWnjeiCS3kZEmmxcD0TS24hIk43r4b4Yl2PenK3kxcUx\nxPhZkLiNIq4RkZpGd8Zm+yuuM6tSa4EytxDjseVtMWmNiBSw8VdcZ1aY5Hxthm5lYgqnEIi0\nrBGRApk4RH4AAAZ9SURBVDb+iuvMEnN9D19N4hQCkZY1IlLAxl9xnVnnAPPzaHPWoeiCo1t4\n8M/8t/Ui6Max9oeL2XG6H/ZIEIjxPRKsyoJzpPLWDE2eI0EgEEkQzru3zNpFppXPRYKZIJIg\nFtxHKpr7SEl+nLiPBIFAJEEIuS8MDiCSIBBJL4gkCETSCyIJApH0gkiCQCS9IJIgEEkviCQI\nRNILIgkCkfSCSIJAJL0gkiAQSS+IJAhE0gsiCQKR9IJIgkAkvSCSIBBJL4gkCETSCyIJApH0\ngkiCQCS9IJIgEEkviCQIuSJt9IgmRfA4rmnWK9fVIv3KFksWS8xoVhSR9p73bWNGs6KItPe8\nbxszmhVFpL3nfduY0awoIu0979vGjGZFEWnved82ZjQrikh7z/u2MaNZUUTae963jRnNiiLS\n3vO+bcxoVhSR9p73bWNGs6KItPe8bxszmhVFpL3nfduY0awoIgGoApEAPIBIAB5AJAAPIBKA\nBxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQCcADiATgAUQC8IAskaqDMYer/U6R\nmKSo1o15SteP+eASeFt8Br3W79xWjVmF36A13VyuEVOWSEnzCwJW4rPmjXTVmEXzRhIy8R8x\nH1RJ4G3xEbRcf0Vv7RtJUHv7uVyhiISJVJhD/Sd/v3ExyfV+TcxlxZhXc3iU1ql+e7WYNXng\n3yD5DJo8klvlplgx5qGJVoRMbk0nlysU0V2YSImpPx2tJBSmfPw9m+OKMfN2MGRZf8S81ysZ\nWKSPoOemqCuTrBjThE/uvZ/LFYroLkykFmvL5qY+Brj2PrvDxny9EzwznZg3k63yq1hW0IPp\nn6WFj/k85Aop70cu1ykieSIV5vQeXucDrBuzpTLZqjEzc1tDJDtoau7HpDmOXS/m8XloF3Tv\n0MvlOkUkTaTHXtk6Zl/tSODjPOHUHA+sFvNozivsA/vJzZsT/zVj3k/11Yak/7HllX4u4xTp\nlCfWx9U6OejGbLglgY8muzGb447wIvWTW19sOAQ+d+gl99hcQQsZ8iOXcYp0r4/d3x9Yax3a\n2TFrqiT0gV03Zlpfg17lHKmT3Poc6Rb6unAn5qnePVWHjyNpj3zkMlqRrOtIyVoi9a5dZcGL\nqxPz0BxHriKStaKrfUpZMdPmMl4VUN7PXK5TRAJFsla5veByC3/VrpPmW5oFvmHYi7nmr9nb\nye29sULM8PJ+5nKdIhIlUnvbwTrYODYfL2XIe4YfMR/hQh/X9WOuItKX5N5Cru1HzHbvEPLe\n1WcuVyiiuzCRmjveVf53BL1WzwY7ZtDK+hKzYY2eDd0VTav6fOW8YszC1H3eisBF3c1lhD0b\nnl2zmkJuU5H+vbFWzMMqe4feenaHVgp6XD+5z35vwT+qOoeQKxSRNJHqfrpp+/HV5qDtLLxq\nzFXOV/rr2R1aK2iZrZ3cZ0/ssDHvPZHWKCJpIgHoBJEAPIBIAB5AJAAPIBKABxAJwAOIBOAB\nRALwACIBeACRADyASAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8gEgAHkAkAA8gEoAH\nEAnAA4gE4AFEAvAAIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOABRALwACIBeACRADyASAAe\nQCQADyASgAcQSTbG3I8mOd7vhTHtD3OfUpO0vxNe5ub5a93G3PJmrOaXyk1WbrW48YJIsjHm\naB7UerQm5fWAyR5DTUP75kOoevBh0ql987TxYscHIsnm4UxV29H8TR47nHqoykxZN53v97Mx\nf2Ol93tirvWb6dbLHR2IJBtjLs3fW/O33iFVj6HK5H8jWGPV/zis2wREkk3jifXXvKjfuJXH\n7CnSq/1xKpVfr5stbrwgkmzGRMreQ9ZYx/psKblttsCxgkiy+RTp3XQw6am89UV6nEYVKedI\nq4NIsumLlP+dAzVvDojUGYKVIOOy6Yt0Nsm1vsidt1cYrv1zpLS9lMceaW0QSTZ9kZ4nRvVJ\nUPE8W7p0RXu9B6uCSLL5EKnu2WAOzcWEgzHZpWx3Tu/2pmcDHq0OIgF4AJEAPIBIAB5AJAAP\nIBKABxAJwAOIBOABRALwACIBeACRADyASAAeQCQADyASgAcQCcADiATgAUQC8AAiAXgAkQA8\ngEgAHkAkAA8gEoAHEAnAA4gE4AFEAvAAIgF4AJEAPIBIAB5AJAAPIBKABxAJwAOIBOCB/wEC\nGw82a8HvZgAAAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title \"Sampling Distribution of the Sample Mean\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dice <- 1:6\n",
    "num_sim <- 1000\n",
    "sample_size <- 100\n",
    "means <- matrix(data = NA, nrow = num_sim, ncol = 1)\n",
    "for(i in 1:num_sim){\n",
    "    sample <- sample(dice, sample_size, replace=TRUE)\n",
    "    means[i] <- mean(sample)\n",
    "}\n",
    "# plot the sampling distribution of the sample mean\n",
    "hist(means, breaks=20, main='Sampling Distribution of the Sample Mean')\n",
    "abline(v=mean(dice),col='red')\n",
    "abline(v=mean(means),col='blue')\n",
    "legend('topright'，legend=c('true mean','mean of sample mean'),fill=c('red','blue'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create three linear data generation models below, each with a different $\\epsilon$ distributions: one unif\\~(-1,1), one N\\~(0,2), one Exp\\~(0.5). Then, we repeatedly sample from them 1000 times and plot the sampeling distribution of $\\hat \\beta_0$, $\\hat \\beta_1$ (so fun!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# three linear data generation models with different error distributions\n",
    "num_sim <- 1000\n",
    "n <- 100\n",
    "x <- runif(n, min=-5,max=5)\n",
    "b0 <- -3\n",
    "b1 <- 2\n",
    "type_of_error <- c('uniform','normal','exponential')\n",
    "output <- lapply(type_of_error, function(x) matrix(data=NA, nrow=num_sim, ncol=2))\n",
    "# simulate sampling from each model\n",
    "for(i in 1:num_sim){\n",
    "    errors <- list(unif_error <- runif(n, min=-1,max=1),\n",
    "                   norm_error <- rnorm(n, mean=0, sd=2),\n",
    "                   exp_error <- rexp(n, rate=0.5))\n",
    "    ys <- lapply(errors, function(error) b0+b1*x+error)\n",
    "    slrs <- lapply(ys, function(y) lm(y~x))\n",
    "    for(j in 1:length(slrs)){\n",
    "        output[[j]][i, ] <- slrs[[j]]$coefficients\n",
    "    }\n",
    "}       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1]]\n",
       "$breaks\n",
       " [1] -3.20 -3.18 -3.16 -3.14 -3.12 -3.10 -3.08 -3.06 -3.04 -3.02 -3.00 -2.98\n",
       "[13] -2.96 -2.94 -2.92 -2.90 -2.88 -2.86 -2.84 -2.82 -2.80 -2.78\n",
       "\n",
       "$counts\n",
       " [1]   1   0   6  14  28  40  54 100 108 145 124 123  92  81  37  26  12   6   2\n",
       "[20]   0   1\n",
       "\n",
       "$density\n",
       " [1] 0.05 0.00 0.30 0.70 1.40 2.00 2.70 5.00 5.40 7.25 6.20 6.15 4.60 4.05 1.85\n",
       "[16] 1.30 0.60 0.30 0.10 0.00 0.05\n",
       "\n",
       "$mids\n",
       " [1] -3.19 -3.17 -3.15 -3.13 -3.11 -3.09 -3.07 -3.05 -3.03 -3.01 -2.99 -2.97\n",
       "[13] -2.95 -2.93 -2.91 -2.89 -2.87 -2.85 -2.83 -2.81 -2.79\n",
       "\n",
       "$xname\n",
       "[1] \"output[[x]][, 1]\"\n",
       "\n",
       "$equidist\n",
       "[1] TRUE\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"histogram\"\n",
       "\n",
       "[[2]]\n",
       "$breaks\n",
       " [1] -3.7 -3.6 -3.5 -3.4 -3.3 -3.2 -3.1 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3\n",
       "\n",
       "$counts\n",
       " [1]   2   7  20  41  81 153 169 199 158  90  54  20   4   2\n",
       "\n",
       "$density\n",
       " [1] 0.02 0.07 0.20 0.41 0.81 1.53 1.69 1.99 1.58 0.90 0.54 0.20 0.04 0.02\n",
       "\n",
       "$mids\n",
       " [1] -3.65 -3.55 -3.45 -3.35 -3.25 -3.15 -3.05 -2.95 -2.85 -2.75 -2.65 -2.55\n",
       "[13] -2.45 -2.35\n",
       "\n",
       "$xname\n",
       "[1] \"output[[x]][, 1]\"\n",
       "\n",
       "$equidist\n",
       "[1] TRUE\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"histogram\"\n",
       "\n",
       "[[3]]\n",
       "$breaks\n",
       " [1] -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2\n",
       "[16] -0.1\n",
       "\n",
       "$counts\n",
       " [1]   1  15  46  93 175 191 193 141  68  39  21  11   4   1   1\n",
       "\n",
       "$density\n",
       " [1] 0.01 0.15 0.46 0.93 1.75 1.91 1.93 1.41 0.68 0.39 0.21 0.11 0.04 0.01 0.01\n",
       "\n",
       "$mids\n",
       " [1] -1.55 -1.45 -1.35 -1.25 -1.15 -1.05 -0.95 -0.85 -0.75 -0.65 -0.55 -0.45\n",
       "[13] -0.35 -0.25 -0.15\n",
       "\n",
       "$xname\n",
       "[1] \"output[[x]][, 1]\"\n",
       "\n",
       "$equidist\n",
       "[1] TRUE\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"histogram\"\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[[1]]\n",
       "$breaks\n",
       " [1] 1.93 1.94 1.95 1.96 1.97 1.98 1.99 2.00 2.01 2.02 2.03 2.04 2.05 2.06 2.07\n",
       "[16] 2.08\n",
       "\n",
       "$counts\n",
       " [1]   4   8  36  50  85 160 170 178 136  89  43  22  15   3   1\n",
       "\n",
       "$density\n",
       " [1]  0.4  0.8  3.6  5.0  8.5 16.0 17.0 17.8 13.6  8.9  4.3  2.2  1.5  0.3  0.1\n",
       "\n",
       "$mids\n",
       " [1] 1.935 1.945 1.955 1.965 1.975 1.985 1.995 2.005 2.015 2.025 2.035 2.045\n",
       "[13] 2.055 2.065 2.075\n",
       "\n",
       "$xname\n",
       "[1] \"output[[x]][, 2]\"\n",
       "\n",
       "$equidist\n",
       "[1] TRUE\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"histogram\"\n",
       "\n",
       "[[2]]\n",
       "$breaks\n",
       " [1] 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30\n",
       "\n",
       "$counts\n",
       " [1]   5  28  68 148 256 224 171  72  22   5   1\n",
       "\n",
       "$density\n",
       " [1] 0.10 0.56 1.36 2.96 5.12 4.48 3.42 1.44 0.44 0.10 0.02\n",
       "\n",
       "$mids\n",
       " [1] 1.775 1.825 1.875 1.925 1.975 2.025 2.075 2.125 2.175 2.225 2.275\n",
       "\n",
       "$xname\n",
       "[1] \"output[[x]][, 2]\"\n",
       "\n",
       "$equidist\n",
       "[1] TRUE\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"histogram\"\n",
       "\n",
       "[[3]]\n",
       "$breaks\n",
       " [1] 1.70 1.75 1.80 1.85 1.90 1.95 2.00 2.05 2.10 2.15 2.20 2.25 2.30\n",
       "\n",
       "$counts\n",
       " [1]   2   6  18  69 150 241 240 174  64  26   8   2\n",
       "\n",
       "$density\n",
       " [1] 0.04 0.12 0.36 1.38 3.00 4.82 4.80 3.48 1.28 0.52 0.16 0.04\n",
       "\n",
       "$mids\n",
       " [1] 1.725 1.775 1.825 1.875 1.925 1.975 2.025 2.075 2.125 2.175 2.225 2.275\n",
       "\n",
       "$xname\n",
       "[1] \"output[[x]][, 2]\"\n",
       "\n",
       "$equidist\n",
       "[1] TRUE\n",
       "\n",
       "attr(,\"class\")\n",
       "[1] \"histogram\"\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAMFBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD////QFLu4AAAACXBIWXMAABJ0AAAS\ndAHeZh94AAAgAElEQVR4nO2djXrbqBZFSXrTtJM0ev+3vf63bAMCtIEjtNY308gWHKQdLSPJ\nrusmAFiN670BACOASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQA\nAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACDIrkDtv0\n9cu5//XekI3hDP4uLfB9PJAewqmRlMHwj7v55pz76L0hGwOR/Jxy2aNIRzgo8iEzP6+5DCjS\neZeOfzr3/eHc7/Mjd2Sa/h2e+fh3euq/t/fDn1/vhxO+4x//+m52C+6BLOTw5+28dkyRDnv5\n/nWIwL1NxxOVf+e9Pz71kMstq+vK21PnY+kUTs2kDIl04mMm0s/xBM+9/ZxX/u/S5Nf5wfDc\nA4nn8MddFoYU6fu0d9/T9Ok+D///ueXy/ZzLOatb+9tTd5GqJmVIpI+f465eHx2e/n1M5uP8\nqvI5/UzHF5jv4x9fQx4zT9wDiefw5r4Or9duUJE+DlJ8ny6X391f92uaTnt/jmSeyzWrW/vH\n42mqnpQhkX6mu0Sn59+OT/0c5/T7yvPikMfME/d9Xsrh6/PXsCKdJp3Tad1xrjmfxp3SCORy\na/94PNVPypBI86XZzj+u9NyDGZRQNM85fL1driiHDOVyhnZcfDtNSAu53No/rayelAWRfgIi\nzV9xbo13KFI8h3f3+f0zsEjXpb8HB/6bbjPSmz+XW3uPSFWT6i7Sn+nnIyDS/Bz4unKPIsVz\nOJ7w/B1WpI/TpeD7dNTm83SKd7pd8HKNNE3nP2/tPSJVTapz+B+Bqfj06N/srsx15R5Fiufw\nv9NK9zNmKOe7cF9Ha34fjpbP28nbP38ut/bzCf10765yUr3D/+3cr++ASPP3Ca4r9yjSQg4f\n7u33v8NhNmYoX+/u7e9RkPv7SN/vh2PmuM6Xy6X97Kn/3k5vGUx1kxoyfBgYo68XNrcKIAQi\nAQhAJIBxQSQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCQ\nJpKbfZkLwBI7PF6S9tW9LACE2ePxgkggZ4/HCyKBnD0eL4i0hh1eC6Swx+OFmw0r2OMBk8QO\nj5c97ascRIIr6ad2O3uFSQGRAjAjhRu55Nb7AZH87DGXcpF2+KrzAhl48Yo0eFarZ6RRg4Fy\nYjPSqMdLmkiH1xFE8jD4q2wpiBRpd/o3IlYUGJE9XgskEXmBGTWq1fslDsZt6UW+07WAm1Nv\nmDoIN9hUDOZEqlO2Dp1OYVxg2Q7XQ7tuLs672AlEWgEi+blJpM8lMBv3jyH1ZkNwBt2zSJ2u\nBbYh0uS/qF4rkn+5fwyJn7VbWyCZbYkUAZEmh0gZzXYtUptrgVhpkzndfomI1LCAv1z/ZBKo\neC2wOGz9YVYQ+S0iUq0C/nL9k0mg4rXA8rDVh6kDItUq4C/XP5kEKl4LJAxbe5g6IFKtAv5y\n/ZNJoOK1QMqwlYepAyLVKuAv1z+ZFOpdC6SMWnmYOiBSrQL+cv2TWQki+UGkWgX85fonsxJE\n8lOywYEP1SFSQrn+yawEkfwUieTvjkgJ5fonsxJE8oNItQr4y/VPZiWI5AeRahXwl+ufzEoQ\nyQ8i1SrgL9c/mZUgkh9EqlXAX87OX30sBJH8IFKtAv5ymxdKncvyLWBDiP+eFiKVlwv93AyV\ncoktm8G9LHjWldRDpOxyiBQuh0jRp/uASJVApAmRmhbwl0OkcDlEij7dB0SqxJ5F4mZDjwL+\ncogULmdfpAiIVKuAvxwihcs9LJv6utELpw1hRmpcwF8OkcLlzM9ObrrK5F9XUs/THZESyiFS\nuNw2RVoxbSJSeTlECp+2bVOk2bqSep7uGxSp3VfzIpK/xLZEEv97WsOI5F4WMguks32RdFf+\nmxVpmrT/nhYiZWzH0/G3XZEu7FykIIgUaK0SKfHnZkAkP8OLFDsX2bNIpedoo4vUMpcEkbq/\nm/a4McFW1W82mBUpmkusl2DgFcsNCOcSO6irz0j9RZqm0LvR6QVWb4dBkaZeuRgXaQrnknJ2\nkzOMv7tVkVzw3/phRsrfgvFFiuUS3pThRXq5X+ZtFFmn2Q5zIsVyifYTDL1iuTotc9mUSEmN\nHnZDcm2Xetu776ldAaOLVAoiBVqvFSnz52ZAJD/ji+TCG7FrkSK5LHUTjFy6XJ+GuWxKJBfb\nino3G8yLFM1lsZ9g6LLl6rTMZRyRUgqsGx+R4iUQ6bU7Inm6I1K8BCK9drcp0lR4C250kTrl\n8lzCnEgtc9mWSLFG4o98eLov/rT2vQRLDC9SIbsWKdaq8Yy0iYPlCCL5GV+k6At+eOuGF6lw\nIhxepIa5bEqk0vFHF6lTLs8lzInUMhdESu+OSPESiPTaHZE83REpXmJbIonfwN+WSH2uBeyL\nxDWSl1gusVfn4UXqVMC8SFHq/j0t0yLFQKT2BTYtkvqAiZRApNfuRkVy12/1Ky2wavzkn83f\nmI3ksmuRGuayKZHc9b/SAuvGNzszxXLZs0jR42XXNxsQyQsi+Wl5vCBSevdNirTnmw3MSI88\nXCNt4NPfzUXi098BIrns+hqpUwH7IsXY8YwUwyvSirtEiJTefZMi7fkaKcauZ6SNfbKh3TVS\n7jv4uvvzpkXikw2PuMijggLF3c2KFBtt9zNSaMRd32wo3IpdiBQYbvciNTpeECm9+yZF4mZD\nxePFzfB3tylS6Yn96CJ1yuW5hDmRqueyrIlNkToVMC9SKaOLFMM3ldzWJRdZWkQkX3dEipfY\nlEixLRlepNirSFKBdeObFSn66qp45Y0Mnbvsv6yoQ3yk8PDDi3Ra4hrplVguilfetBL2Zqfa\nx8tWRXKFWzG6SPFcBK+8S0OvX65C9eMFkYq6b1SkpX6CodcvVwGRHpmd2nH728s2b383OLXj\n9vec1cMOL1IhhRsYumFgT6RCEKlOAUQKd0Ok6KJVkVzky0/qfRTGvkixXGLdCkersFyH2rls\nVSR3/W+hUXDduvHNihTLZaFf4XDy5SpUzwWR0sZ9muQQ6bUbIkUXEWneDZHC3RApumhUpJZf\nZjHlC9NNpMa3vzcjEre/H0kbVn+zYTsiFTK8SIXUF6nhBw59G5M96LqN3YxIpQONLlL1XAQz\nkjmRTut2OiMhkh9EeuQuUmR+cdeGO7xGKp13hxepdi5bFWmp0V5FKmV0kUpBpEBrRPKDSH7G\nFmlhxNunQfYm0opRhhZpqbriLu9mRYqPekplVTD+bvZFio1U720B5bKehVzcy4JnXcoY0cVt\nirRQoLjbpkWSHDDhbogUXdy3SKrP2CFS+rIeRHqhuUjin3XZsEhV3+BHpBcuIpVnPrZI8Vws\ni7R+UyIsHS+7vdnQroBMpA4fp3plGzcb+mZ0JfcXhkhJzbc1M+Wg+QxireUmnAZkRqpaYCiR\nJAdMuNumRXKBkeMb47+yQ6Rw82FEKjtgoiVrLjdhhUgZi4hU42cXEMkPItUvgEjLJWsuN6H8\nI2WIlNt8DJE0n0EMveezWZGmqfQjZYiU23wIkabiAybYdBSRgjQVqeab0k+0Ekn90SArIgVB\nJD+dZqRxRKr108Qbsx4QyQ8irStQTaS8zWgHIvlBpHUFECmtKSIlLyJSzZ9mQCQ/iLSuACKl\nNUWk5EVEqvnTDIjkB5HWFUCktKaIlLyISDV/mgGR/CDSugKIlNYUkZIXEanmTzMgkp9eItX+\nuBAiVWJhg1p+ULX6QZSDgRlpkyLV/ozdZkUysNwFRCor0EogM6+4VxDJDyKVFWgmUtrmtAOR\n/LxuQODME5G8qxEpvHrvIvmXEcm7GpHCqxHJt4xI3tWIFF6NSL7liiLVuImJSJVAJD8WRIpt\nTSm1RGp92/v208rdO88GNH3vKGG5C4iUV6DXTGRmZvKJZGy5C4iUV6C3QAYPGBPydP+UgzGR\nZGmkdY+MZVWkJgfLci7WTueCUkmDyjpejMxILURyLwsv67pdE6X+rHHApOTia25wWTqF5xwv\nwdeX5iIFeN0FL+UiPQ4U2gx7JAZDLuRyInefs4MBcgmwx1wQaQXk4mePuaTtauY0txvIxc8O\nc9nTvgJUQyVS60vAbET7uYbeEbzQO5AzY+y3TKRB2tUktA3Wnm9M8WaUdkSkBu1qYk0YC5lM\niFRWx3q7mlgTxkImEyKV1bHeribWhLGQyYRIZXWst6uJNWEsZDIhUlkd6+1qYk0YC5lMiFRW\nx3q7mlgTxkImEyKV1bHeribWhLGQyYRIAHADkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgE\nIACRAASsF2n2d3ejf433YWVk2IJ6knY1yd2Ghu2nri+mt6Fz/wr4Q8f0Xo85KXd8dS13rzJb\njLY77oWynqRdTXK3oWH7KfrLqM196MxtePxG36wPXLqiAdNqK6tEC952YmHY1INC3a4ms23I\naJ/eKXcfDYg0Ow7ytuHecd8iLQ6LSO1EWnxVq4l7+lnWsUAk+T4LCj6ea0YvVpKGvbXLEST+\nsnot2Osaab5PSduQmulro+R9tCVS7q9mtUjqY0E9I8XriWek2a8hpZ2JGcmlbkNqpo+Nkvdx\ndiRaEGn+oFrH5xlJueOraj3OHeFyhe2CO3tp97AqoV0XkVL3KdB+oXX+Pj5n10mk1F/hUsf0\nF47ri23uKXMi2hlpqZp0RnKRR+FBe89IyduQmuljo8T69/bVvno0iXyRntqWzkjZA6Zu0LoC\nzy+iC+0WD/uUg+7x1CfFt24i5Z5SpGb62CG9fvJrUHU4tZtVcNfz/vir271dfNhZu6Vyzt3a\nRcvN2vU5aBL3ydM+acbI3cd5+6m/SCW/mrKOub+H3A0CgHUgEoAARAIQgEgAAhAJQAAiAQhA\nJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAA\nIgEIMCNSyoY8fn1U1282bAa5+LGWi4nIE74Ecv6w0lf82YNc/FjMxUTiFoOxALn4sZhLr8Sv\n34Xu3PWP43/nR48r7k3vwdy+JnS4A4Zc/JjPpVPily+Offpx/f95xWXx9u80nEPsuPnVIBc/\n9nPpL9Lki2F6eeZxxQ4OmIlc7tjPpZdIbj4Ne4J5/Fb4/Rww5OLFfi4dZ6T7ondi9r6szBbH\nPGAeFsnliv1cOor0sv+PwThvVPfOwx4w5PKK/Vx6JT6/C/N0xXiZoK/hzJq6W1/38MRAkIsf\n87mMFzlABxAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAE\nIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQg\nAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQA\nAYgEIACRAAQgEoAAgyK5wzZ9/XLuf703ZGM4g79LC3wfD6SHcGokZTD8426+Oec+em/IxkAk\nP6dc9ijSEQ6KfMjMz2suA4p03qXjn859fzj3+/zIHZmmf4dnPv6dnvrv7f3w59f74YTv+Me/\nvpvdgnsgCzn8eTuvHVOkw16+fx0icG/T8UTl33nvj0895HLL6rry9tT5WDqFUzMpQyKd+JiJ\n9HM8wXNvP+eV/7s0+XV+MDz3QOI5/HGXhSFF+j7t3fc0fbrPw/9/brl8P+dyzurW/vbUXaSq\nSRkS6ePnuKvXR4enfx+T+Ti/qnxOP9PxBeb7+MfXkMfME/dA4jm8ua/D67UbVKSPgxTfp8vl\nd/fX/Zqm096fI5nncs3q1v7xeJqqJ2VIpJ/pLtHp+bfjUz/HOf2+8rw45DHzxH2fl3L4+vw1\nrEinSed0Wneca86ncac0Arnc2j8eT/WTMiTSfGm2848rPfdgBiUUzXMOX2+XK8ohQ7mcoR0X\n304T0kIut/ZPK6snZUGkn4BI81ecW+MdihTP4d19fv8MLNJ16e/Bgf+m24z05s/l1t4jUtWk\nuov0Z/r5CIg0Pwe+rtyjSPEcjic8f4cV6eN0Kfg+HbX5PJ3inW4XvFwjTdP5z1t7j0hVk+oc\n/kdgKj49+je7K3NduUeR4jn877TS/YwZyvku3NfRmt+Ho+XzdvL2z5/Lrf18Qj/du6ucVO/w\nfzv36zsg0vx9guvKPYq0kMOHe/v973CYjRnK17t7+3sU5P4+0vf74Zg5rvPlcmk/e+q/t9Nb\nBlPdpIYMHwbG6OuFza0CCIFIAAIQCWBcEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQ\nCUAAIgEIQCQAAYgEIACRAAQgEuhxsy//2Ql72ldohHtZGJ8d7Sq0ApEABCASgABEAlDAzQYA\nKAGR1rDDV94ULt/NuKtc9rSvcvZ4LZCCm64ydd6QhuxoV/Ugkh+vSIPP3qPuVxMQyU9sRho1\nqlH3qwmI5Of0j7aeFjzrWm9MI0bdrzYMfrqygvM/1eNb0XxT2mB0v9yM3tvSD/dM7w0SIN0H\nQ/EY/d0476I1at/mfa5rOIpHrpnUPrV7LNZVK6O/m+2IVPM271ZFumXSWKR6Ay1j9HezYZGE\nL4ibFmnyXyQhUlM2LNJsnaK+umYL3PUHIvVmIyJVvs27cZG8/yYYItXHf6VoaQtfqXqbd6si\n3U0Kr5KOs/yoOpZ+N357LG1hBrsWKQIi1Wd7IlW+zYtIWcUQ6cLmRKp9mxeRsooh0oVtilTx\nNi8iZRXbgEht3izeqEj1bvMiUlYx+yK5l4UqbFWkard5ESmr2K5FCnw6ahsi1b7Ni0hZxfYt\nkn95IyJFQCQ/iJTcunh8RFqosdEoHti1SBVvNiBSRo0R/nrSvkVqMz4iLdRw8dWbYNcind94\nZEbKAZH87F4kl9y6eHxEWqiBSNFiGxVJc6KOSBk1EClabKMi5RRIGx+RFmogUrTYBkSq9xfY\nECmjxqZFqnGrcWsiTVO1v8CGSBk1Ni3Shb2LVKsAImXUQKRoMUR6XUYkXw1EihZDpNdlRPLV\nQKRoMfsixT6Ugkh+EMnPrkWKtUKkJ9p9QeTWkjmxb5EizRDJDyL52blI1QokiLTNjzsjkh9E\nqlMga0ba0pGDSH4QqU4BRMqogUjRYoj0uoxIvhqIFC1mQ6TSC5DRReqUi7eGJZF65hIoZkOk\naSq7mB9dpE65eGtYEqlnLoFiT4/a/rX818FXFVg1vlGRph65eGvYEmnqlkugWNf56XVGyhxy\nfJG65OKtYUukfrkEitkQ6fbSgkgPdMrFW8OSSNFcIudTw4vUqYB5kUoZXaQY7mXBs046UOaj\nCiBSJRBp2qlIrmy44UXqk4u3himRIrnsWiRXON7oInXKxVvjWaSeX7waywWREOmFaC6VL6rz\nZiQ7Iu36ZgMi+en5yrtZkRb7STci/1EFZtdIZWcHo4sUy2XPIi3mstcZKd6q2lRtX6QIuxYp\nwu0+xB6vkZIa6YNBpJT6/sfbEkn0FdeR+ytGRIrs5q5Fiv36d3yzYel4qTQjlarT+mZDfqPR\nRSodc3SRosdLm6+4RqRA982JVPuieqsiTVOTr7jenkh7vtmwJFLNi+oNiyTv5+1uUqTCS8Hh\nRcq9FhB+zsC0SP2PF5sixVvtd0aKIZ+RFj70Y0mkUnYt0p6vkWLIL6oX1NiKSLGP/40v0u2o\niDTSn8LYFymSy6S+qN6SSLFcIpsyvEju+t9Co+C6leNPRkWK5bLQr2iwWA1LIsVzibzyrB5W\n8KgCiLQAIvlpmkuo+9ZE2vPNBkTyg0iPPFwj8elvDwmfEJJdVG9IJAPHi1GR+hSwL1IM9UX1\nlkQqBJHqFNi2SOqLakTK625RpOg7+G3eFzApUtN38DckEp9seCRtOPUpjL/7skjhN/3r0upa\nYEMilQ65D5EKTOo0I7W+sm7Q47WTeZG4axcaoE0wgZkFkRApr7tFkUpPmIpE8nc3KVLLXLYk\nUtNcQt0titSywJZEKmV0kUpBJGEBRErrhEgL3S2KVHo/bHSRWuayJZGa5hLqblGk0xLXSK80\nzGVLIjXNJdTdokiucLzRRWqZy5ZEappLqDsiebsjEiLldbcoUstP825JpKafct6QSNm56P9G\ntVGRGhbYlEiFDC9SIYgkLIBIaZ0QaaG7SZHiX/KRUKCoj32RcnNZcwqzJZEaHi/B7hZFctf/\nSguU9TEvUstctiRS01xC3RHJ2x2RECmvOyJ5uyMSIuV1tygSt78DcPvbD5/+1g4wvEiFDC9S\nIcOLVDrS6CK1zGVLIjXNJdQdkbzdEQmR8rqbFCn27kfkvZHhReJbhLzwLUKPJA0Qm7ZGF6mU\n0UUqBZECrRHJDyL5GVukhWF2K9KKUYYWqXUuoe4mRYoMtWuRCkcaXqSWuYS6b02k3d5sQCQ/\nLXNxj0SKbUAkb8/M+zb+NBAp0ilTpGdKtiCNpiJVeFSB9FO7tTPSshCIJJ2REGnNsLlcRIq/\nit1yW3ONtEmRyl/d03rEJxG7ItXOJdih+FHtyTqp5G5FWkGiSDkP7Yi0hPKausaMhEj5w1Zl\n9QEzpkjuZcGzLrtYtUcS0kS6/q1iQyK1uaqOs/6AQaT0YtUeSUgteTpcV03V9WYkREpp3RBE\nqlkAkWIVUh4iEiJNiBSvkPJwMyJxs6FmgQFF4mZDBqUXtIh0Wrv0cYZti+Qj64AZU6TTyMxI\nygKlQmxDpPUHzLgiucAmIFJRgfFFWnfAIFJSsbqPJCDSChDJDyLJCyDScoXkh9sRSfEGvr8D\nIg0p0uoDZkyRpknwBr6/AyKNKNK0+oAZVqQgiFRUYHiRgiCSH0QqKoBIWa0QaakDItUUycAH\nwV9AJD+IVFSg9Yxk5+BBJD+IVFQAkZ6ezfq75Yj00kH0SP8XzxGpEgGR4o8RaamD2fkJkSqB\nSH4QqahAc5HM3HcwIFKFE5j1IFJRgY4zEiKlbFBzEjYj9gKASIjkexaRFpsg0tNaRPI8i0iL\nTRDpaS0ieZ5FpMUm+xUpcIKLSJ5nEWmxyY5F8i8jkudZRFpsgkhPy4jkeRaRFps0eKR4m0Al\nUmBLdi/Swps5LUWy8baSd9y+N7wl81Nat0j47uHHc8meIjU4aLJy8TzsOCNVFWk5l9iTlh4l\nk9TNvSy8rDMo0myxklR5ufgadhSp4gyVkovZOagwl3KRHgd6+b2YJTEYciGXE7n7nB0MkEuA\nPeaCSCsgFz97zGX1zYZdQy5+dpjLnvYVoBqVRWp9ZWj9daF5HukQhJ/UHaoaV3H51v1asbB9\nq1bXrC0nNpyldcm5IFJTEClhOEvrEMkmiJQwnKV1iGQTREoYztI6RLIJIiUMZ2kdItkEkRKG\ns7QOkWyCSAnDWVqHSDZBpIThLK2zIhLAPkAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEI\nQCQAAVVEuhW9/FXdtL+yO2s065cwmreftVeI50xyVr/uYWCt7/cZ77y0Ws98tJzxFgIM9YqH\nk9Qx7ejNqJ2+DfPiLnEQd/vjujh7JqufxWn2JZOM1Z499K+d1UnsvLRaT+ir4ZP7ZW3mQjhJ\nHTMOQy1uet7njBeQeb+MyKyL9JpJzur7c6HV7rmOd3W8dnC1ksctzDkwikSajxIKZ7Fj0oBV\nT+3mQuT0XCmSPY+ml0wyVyeKFIrLjkjlg60WKX/nTIl0mE8zToY1Ihm8RJr9NqPXSJfIQg3q\nitT0GimwLYv9ggktj1gqUq9rpIeDJksIN/vRol9DEmekS2Sh9Uu9/budPiE1nZHyxlpKSD5g\nXi51RZo9SBnId0pY0C95vJakihRf7e+97sxtYWg1q0XK7Fp83po3oCWRHqfgdJGe+6WO15RV\nInn30Ns8dtqX4FFLkVbMDxl9i6fbzFyqiuTO/ycK8fRSlSxSab+2JIrkfGtf9jCwNlrZu3Kp\ndgXWiuRNaHm03BFzc6kn0mmHz1dpye/HOnfrN/tRp19j5pmEpo1ZZA/rXvcwsNYn4X217y3Q\neO0a3Pe0tF/m+7GRcNI6dntDFmB3IBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAA\nRAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCDAjUsqGuKemZja+\nIuTix1ouJjIPboT/qxBvXxVoYuMrQi5+LOZiIvOyYIz++y1CyMWPxVx6ZX7/DtbrH6c9ddfX\njfuKe9NZItnfnLkVyMWP+Vw6ZX75ytqnH9f/n1dcFq9T8+07aMc7YMjFj/1c+os0+WKYXp55\nXLGDA2Yilzv2c+klkptPw55gHr/CfD8HDLl4sZ9LxxnpvuidmEMvK2MfMA+L5HLFfi4dRXrZ\n/8dgnDeqWedRDxhyecV+Lr0yn9+FebpivEzQ13Be/5Gl29J4Bwy5BDCfy4CZA7QHkQAEIBKA\nAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEA\nBCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgE\nIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAk\nAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAi\nAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQ\nCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCA\nSAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEgPA\n2tkAAALDSURBVIAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIB\nCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJ\nQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBI\nAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABE\nAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQg\nEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAA\nkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAAB\niAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEI\nQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlA\nACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgA\nAhAJQAAiAQhAJAABiAQgAJEABCASgABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQC\nEIBIAAIQCUAAIgEIQCQAAYgEIACRAAQgEoAARAIQgEgAAhAJQAAiAQhAJAABiAQgAJEABCAS\ngABEAhCASAACEAlAACIBCEAkAAGIBCAAkQAEIBKAAEQCEIBIAAIQCUAAIgEIQCQAAYgEIACR\nAAQgEoAARAIQ8H+9QyePqReE5AAAAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title \"exponential\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "par(mfrow=c(3,3))\n",
    "lapply(1:3,function(x) hist(output[[x]][ ,1],\n",
    "                            main=type_of_error[x],\n",
    "                            breaks=15,\n",
    "                           xlab='estimated b0'))\n",
    "\n",
    "lapply(1:3,function(x) hist(output[[x]][ ,2],\n",
    "                            main=type_of_error[x],\n",
    "                            breaks=15,\n",
    "                           xlab='estimated b1'))  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Inference on SLR Parameters & Estimating Uncertainties\n",
    "The analytical solutions to $Var(\\hat \\beta_0)$, $Var(\\hat \\beta_1)$ are:\n",
    "$$\n",
    "Var(\\hat \\beta_0) = \\sigma^2[\\frac{1}{n}+\\frac{\\bar x^2}{\\sum_i (x_i-\\bar x)^2}]\n",
    "$$\n",
    "$$\n",
    "Var(\\hat \\beta_1) = \\frac{\\sigma^2}{\\sum_i(x_i-\\bar x)^2}\n",
    "$$\n",
    "Using a bootstrap procedure allows us to use $Var(\\hat \\beta_i^{B}$ (the variance of the estimated betas we get from all bootstrap samples) as an estimate for the true variance of betas. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol>\n",
       "\t<li><ol class=list-inline>\n",
       "\t<li>0.00345161807048146</li>\n",
       "\t<li>0.000512537235367081</li>\n",
       "</ol>\n",
       "</li>\n",
       "\t<li><ol class=list-inline>\n",
       "\t<li>0.0429277151157414</li>\n",
       "\t<li>0.00616828119169228</li>\n",
       "</ol>\n",
       "</li>\n",
       "\t<li><ol class=list-inline>\n",
       "\t<li>0.0415703726147198</li>\n",
       "\t<li>0.00631599062297564</li>\n",
       "</ol>\n",
       "</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate}\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 0.00345161807048146\n",
       "\\item 0.000512537235367081\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 0.0429277151157414\n",
       "\\item 0.00616828119169228\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 0.0415703726147198\n",
       "\\item 0.00631599062297564\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\end{enumerate}\n"
      ],
      "text/markdown": [
       "1. 1. 0.00345161807048146\n",
       "2. 0.000512537235367081\n",
       "\n",
       "\n",
       "\n",
       "2. 1. 0.0429277151157414\n",
       "2. 0.00616828119169228\n",
       "\n",
       "\n",
       "\n",
       "3. 1. 0.0415703726147198\n",
       "2. 0.00631599062297564\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[[1]]\n",
       "[1] 0.0034516181 0.0005125372\n",
       "\n",
       "[[2]]\n",
       "[1] 0.042927715 0.006168281\n",
       "\n",
       "[[3]]\n",
       "[1] 0.041570373 0.006315991\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>true_vars_b0</th><th scope=col>true_vars_b1</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>0.003361732 </td><td>0.0005199702</td></tr>\n",
       "\t<tr><td>0.04034078 </td><td>0.006239643</td></tr>\n",
       "\t<tr><td>0.04034078 </td><td>0.006239643</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{ll}\n",
       " true\\_vars\\_b0 & true\\_vars\\_b1\\\\\n",
       "\\hline\n",
       "\t 0.003361732  & 0.0005199702\\\\\n",
       "\t 0.04034078  & 0.006239643\\\\\n",
       "\t 0.04034078  & 0.006239643\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "| true_vars_b0 | true_vars_b1 |\n",
       "|---|---|\n",
       "| 0.003361732  | 0.0005199702 |\n",
       "| 0.04034078  | 0.006239643 |\n",
       "| 0.04034078  | 0.006239643 |\n",
       "\n"
      ],
      "text/plain": [
       "     true_vars_b0 true_vars_b1\n",
       "[1,] 0.003361732  0.0005199702\n",
       "[2,] 0.04034078   0.006239643 \n",
       "[3,] 0.04034078   0.006239643 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# variance of observed beta1 hat and beta2 hat\n",
    "observed_vars <- lapply(1:3, function(x) apply(output[[x]], 2, var))\n",
    "# true variance of beta1 hat and beta2 hat using the analytical formula\n",
    "sigma2 <- list(1/3, 4,4)\n",
    "var_b1_hat <- function(sigma2){\n",
    "    sigma2/sum((x-mean(x))^2)\n",
    "}\n",
    "var_b0_hat <- function(sigma2){\n",
    "    sigma2*(1/2+mean(x)^2/sum((x-mean(x))^2))\n",
    "}\n",
    "true_vars_b1 <- lapply(sigma2, function(s2) s2/sum((x-mean(x))^2))\n",
    "true_vars_b0 <- lapply(sigma2, function(s2) s2*(1/n+mean(x)^2/sum((x-mean(x))^2)))\n",
    "true_vars <- cbind(true_vars_b0,true_vars_b1)\n",
    "observed_vars\n",
    "true_vars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IV. Prediction & Confidence Interval\n",
    "The graph below shows the confidence band for the line of best fit created from 1000 bootstrap samples. Notice that the further away $x_i$ is from $\\bar x$, the wider the band is vertically, implying a larger uncertainty.\n",
    "Given a list of xs for which we would like to estimate their corresponding y values, we compute the variance of all the $\\hat y$s predicted by the bootstrap procedure. The second graph supports our observation from the change in the vertical width of the confidence band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD////ojgWfAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAgAElEQVR4nO2djWKkKrZGrfx1mD5JeP+nnU4pioqKuoGNrnXvJJaVLqDi\nOh9usdJYADhNU7oDAFcAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJ\nQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAE\nQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABE\nAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQA\nARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAA\nkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJ\nQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAAkQAEQCQAARAJQABEAhAgg0gN\nQBizus/420/aHc+NtBw4yuXFKdAE1EnggDXelrf9PLrbA7zdQiSAjlWPrCnnESJBRYQO12Ag\nZfcIkaAe9nlkcnqESFAN0R5Z7+DO5BEiQS3s9sjYfB4hEtTCvoJdZo8QCSrhqEeIBDCwr/Cd\n3SNEgirYOEGK8SitUvuHhEiQnX0e2W2PxK3aPyZEgtwcuRCb1yNEAv3s8Mj2R3VejxAJ9BNb\nsCvoESKBeuI9Gg7qsS6stQOILXz3utgCHiESKOdo4TuvR4gEujntESsbAATyKBP7h4ZIkI2j\nF2Kze4RIoJldBbuSHiESKKYejxAJ9LL3ApJ1HjG1y9YE6Cf2BMkdyqZcHiESqOW8Rzmd2j8+\nRIIc7PTIlPUIkUAnBzyyBT1CJNDJeqFBnUeIBCo5VLAr6BEigUaiPXL17umHnGSv3e0fIyJB\naqJPkNR4hEigj10eGRUeIRKoY2/BzirwCJFAHfsKdtZo8AiRQBtVeoRIoIw6PUIk0EWlHiES\nqCKy0NAJM1+oWsojRAJN7PNovlDVjr7lZP9YEQlSUa9HiASKiD1B6hc0qPEIkUAPOzyyJrBQ\ntTuk81tkCoj0+dI0b3+TNgF1ss+jttCgxaOMIjXPf/jaPPlI0gTUzJ55nTqPcov00Xz8WPv9\n0XymaAIqZrXQsODR6Ej2vhZg/4BPifRofn63f5qXFE1AvWwU7IaNzhdlHuUWqWm8B+JNQLXE\nFr7HHnl3mpvx4+zsH/Ipkd6dSI8UTUC1xJ0g6fUoq0hvfz7/Nv/7t/nzsV5tQKS7EeuRfyHW\nP4rzizNl/5iPi9Ty3Hz8pGgCKqV+j3JeR/r6+vx8e3uWHD5WPUKkmxHpkSnjUdzL7x81KxtA\nlrjCtxkuxB450A8T+fL7h53oKG980jQBKtlR+A55pIT94z58lP+8N81rtziI8jd0XMOjjCL9\nPJ5h89a+CCJBS3zhu/NoapIOs/YP/OhR/lwW9PP5eH2+CCLBk9MeKWH/yI8e5Y/2H34/Xr4R\nCTr2epRTo11t7R/6udXf/0Lp9RWRoGXXgoZ54VsR+8d+9Ch/adzFo5dXRIJfIgvfFXiUUaTP\n5r3b+m5eEQniPVq4EKuK/aM/fJR/9Pb83bhUhEi3IK7w7W7gU1tmaNk//ONH+deb2/p+RyRY\nP0EavtfgkZ6VDZmbgOKc8kidVfvHj0ggQVTBbuJRr486jxAJyhDnkanFI0SCIkQV7CryCJGg\nBHs8Gr64Y7akMEvsfw8QCU4TWWioxyNEggJEe2Rq8QiRID87PBpOk0z/UCX73wVEgnPEe9QX\nGtR7hEiQmzMXYtV6hEiQmbiC3cyj8h8Buc7+NwKR4ARX82i4xWMviAQniJnYVeRRz/53ApHg\nOFf1CJEgJ1EeeReQuoPU+6oJv0f73wtEgqNcy6MR+98MRIKD7PLIjg1S7hEiQTbWC3bdt1o8\nmnZo/9uBSHCIHR6Z6jxCJMjF2sSuOxCr8WjO/vcDkeAIESdIFXuESJCHC3kU7M/+dwSRYD9x\nHvXiVOfRgUMWkWA30Qsa/AtIOj0KEqykbIBIsJeogt3UI3eIqiI8rTNM7SADuzyy2wduORY9\nQiRIz2rhu/tmgh7VQDcD3Q0iwT52FuzUshxHBpEgObfwCJEgMffwCJEgLXV7tP7pRdyPBLmI\nKNgp9qhn0yNEgqRsF+yu4RErGyAl8R7p08gGtkLPmufapv3vDSJBLNt3TgweabuEtHB6FIpP\nVjZAUqI86i/E6vKoJ+L06PllN4gEcWwX7Cr1aDqta7/tBZEginiP1C0N2jmtMxQbIBnbhe9B\nH1UWeWz3y8URiQSJ2CzY1elRcFpnmNpBIuq/gGS3PfK+7X+HEAm2OeaRJqU2L2/1VYbu214Q\nCTbZcwFpeuzqYNoXO903+c/A/vcIkWCLfRdil47dgkRePBr2UbUDeS7iUcS0bni4/11CJFgn\novAd4VFBr3ZP6wxTO5Bnu/Ct2qPYad34qu3+twmRYJVoj9YWWCvzaPQDs2md4RwJxKnco2n7\nMdO6Z7buBpFghc1Cg7/Azn0v7c06owr9PI4s50ggzrZHs4VBdv5f/VJsT+vGP+dZtf+tQiRY\npG6PWoIrvkNPjYoO+98rRIIl1grfY488n1R7NHpmMq3zH1FsAEFWAqkWj5ZZiqPn0BAJ5Ngq\n2Kn2yPZfxruGB9NZ3nAhCZFAkP0eBarLpQh4NHs6PK17jg+RQIqtQoPnz9oZfEHWPArEkV9l\noNgAQhz1SANbcRSqMvTPPcdHIoEMMRdiQweiBlY8Ckw950Vvig0gxmbBLmJhUBnsVh6FnvM8\n6uqRu0EkmBPpkVXoUf9l+enAQ+ufHRkuyIII0R7Z4CFaUKnd0zobiCPOkUCEyMK3Po+8aV3E\n0tmVOEIkOE+9Hvntz7sRjKNu+zm2YVqHSHCayMK39adRejxantatPNEP0m0jEpxla6Wqt6Ah\ndJwW82hjWheKI5ddbmhDlYFiA5xlo9AQOkqVeTTtxWKVwf1LOynWUWyA08R75J+rBzbzcmRa\nN1QZZnFEIsE5Ij2yujxaiaPAnvm0zo7iiESCs1zOo+g4sqM4QiQ4Q7RHwfUMxTxabnzhAmwg\njnqpKH/DWTYK3xV5tNYXa+xaHJmFd2ILRIKOjcK3To/2Tuv60r0b1hBHfZWBcyQ4w/rELuRR\neVbiaHtaN4kjM8QRIsFhtjxy/y2vxKPlf7EVR+H3YgtEgifbhYb2oEvhw0FWbj2azDl3xRFT\nOzhOfOFbDSGpF68lud1xcYRIcIwaPYrc1z/lJqd2I47Cb8cWiATrHrn/Wz9OsxM1rRsezore\nz51upGakEFM7OMZK4XvqkQ0frblZmdYt/Yshd7bjCJHgCPGF78Ie7Wrev9w1mda57XEcDbmE\nSHCAajxyrFcZAleJd8WRRSQ4wv4LsWo9Cu3YHUcWkeAAsR6tHsgZiEvDgFGDOrM4sgvFOkSC\nvVzAo82LR8aPo+m0bhZHFpFgN5EeBc47yhh1bFrXDaddoNGP0qt3T94JRIJ9rBe+e4FKezSJ\no9Xsme6yUXE0ficQCfaxHEibHmXk6LRuqDJsxtH4jUAk2EX1Hi3+7HxaFx1HFpFgH6snSO0S\nO5Pfo51lucXdXpXBLsSRCcWRRSTYxV6PtMXRyrTuTBxZRII97M6jEuyZ1s3/nRvTrjiyiAQ7\nqMqjdju6M6txZDbiyCIS7GCj0KDCo47d07r2u12II7MeRxaRIJ5dHpXx6dS07jmSaRy1/7Pr\ncWQRCaKpzKPAtG7h0tFyHJnIOLKIBLFU5NHOHxhdgz0URxaRIJKNQsP4GC16nrSzca/KYI/G\nkUUkiGOzYKeG+Muv3n47jSO7K44sIkEU+j3aN63rp6CTOJoUvWPjyCISRLF6gqTBo3U2qwwu\njiZFb6fRVhxZRIIYavHoyLTuOQqzGkcRHiESbFOLR3GMK4s2Io4iPIpxbQoi3Ywtj3SbFFdl\n8OPI7o4jywdEwiZxHum2qcdd7RKNo9+fQSRYp2qPpt0apnVeNe50HP1+QSRYZbXwXZtHo/12\nHkdubdHeOLKIBBusBJLzyDtLUipUS7jKcC6OetUQCdZY88jo9mh5Wmfk4qjfQiRYIcIjo8Oj\nWds28HBWZQjEUff0rjiyiARrRF1A0unRwvNuVepyHNnYOBr/CCLBIvV4NGUhjrpH0zjqB7Mj\njqamIRIssc+jgkRN68xQ9JaIo+kORIIF6vVo7XmbIo4sIsEiiyIp92ixyjCKo2mVYUccBX8C\nkSDMRsFudJAW5FCVYZDhSBwF9yISBIkpfCv0aBpHo+TcjqNtj5YCC5EgRNQFpPIcmdZ5cWSm\ncbTl0fIPIBIEqMOjrcr7pMoQjCOzK44OPbUEIl2eyjwKP3YncxniyHI/EoSYHzEVejTaWoyj\n55fwmCfvwJpGTO0gwHLhW49HG9O6ocqwEUc2No6Wn3k+hUgwZcUjLRrFTutyxFH7DZFgwhU8\n6jcGj07E0cas7gkiwZj1lUGJBYlj37RuEkfdVaM9cbT8zPAUIsGImjwKP7YjzwJx1E3rTsfR\n6BlEAp/qPRrHkQnFkdkRR5EaIRKMWS582/7MQjPScRT9DCKBx3rh25jiHu2qMsziyCSKI4tI\n4KPeo3XGcWScQwfjaI9GiAQeq4VvM9iklN6j4TxpXvSOjaN9GiESDFTrUSeH97BbzmDHceRM\nWhrt+nux/gwiQcdqwU6rQz2RcdRrtOrR3jiyiASOKI/U6uR7dDaODmiESODYLnyrZGFaN4mj\ndl9r2GYcLT6x9s8QCZ4snyDp8GitD9Np3SiOrNPoZBxt1PgQCX7R7lEYV4IbHp6Mo5VZ3cYb\niEhgVz2yk+NWEeNpnasunomjxSe2PEIksOo9Cjc8rF1wG6niaFsjRAKr3qMwMXHkFRuXBjqM\n+IxGiASxHmlkuEC0HEcmMo52PzEGkWCp8K0ihOKmdStxZDPEkUUkWAwkFR6FmU/rzJk4Oq8R\nIsGqRypV8pKm2zgZR0v74zVCpNuj2aPVaZ175MXR9IaJU3G0SyNEujurhQaNeRQsendxNPEo\nJo6WZ3U730hEujXVFezG07qoOFr1aGn/Xo8Q6d4sFuyKe7RrWmfX4mjZCaFZ3RNEujN6PVph\nLY7sjjiS1AiRbs1ioUGnR6JxJHZy1IFI90XtBaS15icejePIxMfR0v6DHiHSfanRo9m0zorG\n0WGNEOm+rHikcmL3S0Qc2e04SqARIt2WTY8U2jSOI3M0jhZ2n9EIkW7LQsGudB6tTetMwCPX\n9/NxdO7tRKSbotSjFbp+dcFzMI6SzOqeINItWZjYKfbIrMZR/9TS6BaH7V7gLIh0R3R6tD2t\ns9bf6Dp+No4kNEKkW6LTo0X6To2mda7jw971OEp1ctSBSPdjySOtKxp+8ad13hh2xNHCbiGP\nEOl+qPQoZlpnzOK0zhyLIzGNEOmGhAt2SvPI86Td8kfgFb2XhuZ+bsfuY2QV6b8/b80vbx//\npWoCtqjKoyfRcbTrLEhUo6wi/bw0A69JmoBNFBYaFhvul1ikiKOT7+OUjCJ9NI//fT23vv8+\nmo8UTcAWCj1aZEihfmN3HOWY1T3JKNKj+eq3v5pHiiZgg2iPtGg19ch1erQ3PDCbUaOsIjXN\n0gOxJmAdfR6tT+u2rsGajTgK702gEYl0K1YK33m0iSc4rXN9PhdHp9/GIHnPkf5+P7c4RyrD\ncsGupDMBxlUGV/3uuhwXRxlndU9ylr9fvardy0+SJmCFRY+0IRFH4b2pNMp9HenjeR3p8faH\n60j5CU/shg9UzS1UxOmREY2jhBqxsuE+aPPIw84f2JEx4ziyG3FUQCNNIjU+aZq4M5o8Wmtr\nNY687rY7Ike6uFeOnCL9vDfN69/uRSh/Z6Uyj2ynja0jjmzeJUKPdqFd+yKIlBNNHgXMcVvu\nQR9HY4+246iURpnL35//bPp8PJfZIVJWggW7Ih5FT+vceVLXXeMiajWOgjszaJT5guzz2/fj\n5RuR8rLgUYHS96ZHfRxZuxJHQWX2nDKJU2CJ0M/rKyJlJTixK+LRmEC1zosjsxJHoREVjCOb\nVaSXxl2EfXlFpIyo8ShqWmcW4qg7ZVoYUMlZ3ZOMIn02793Wd/OKSNmowaNpHNl5HJndcZRR\no7zl74/enr8bl4oQSY7Fgl1Zj8LTuvar7Srgz75GxFF5jTJfkP16c1vf74iUiaWCXeHzowCH\n4yhu2GnRs7IhcxN3YdGjvCJFT+ts17u2q5XEkUWkq6PSo41pnbckyO1xO+aj0aERIl2chUJD\nbo/GhNo+FkdqNEKka7PkUd4TpJhpXft1lD7tHrMWR5FjzgEiXRiNHgWb7uPI1BlHFpGujEaP\nln/Cjj0yG3GkSyNEujKhuVB2j8LY6fazd6M4MutxFBpvOY0Q6cJo8Gh5IjfZdT6OimqESNcl\n3qN0ZkVO656dm1yD9XbMxhG91i4jiHRRwv/NzuvRQqV7tsOOpnVmpFHAo9BoS3uESBcl3qPS\ntF1bjKPpQBTO6p4g0iVR4FFkU13XXB/Nehxp1QiRrsmCRzlXqobbmu5zXeu+t7vMchzFjbUE\niHRFglcvM6/4jmmr65nroqk1jiwiXRINHkXgetZ9N24NoB094Q0hNCwtINL1UOpReFrX97Dd\ntRxHgXHq0QiRLkjwkMvr0d5pXftv9saR5Jt2GkS6Ggo8imDUVbdrKY6Uz+qeINLF0OlRaFrn\nVxnW4qgGjRDpapT3KKKlrlt998xqHIVGJPeGSYFI1yJUaDBW37TuRBxp1AiRLoZGj9amdW7P\nQhzprzH0INKV0OjRhK5XQ+/W4igwQqUeIdKVWDjwSrvj4/rUfm/3xMeRWo0Q6UoseVTQpI1p\n3SiOxv2vSyNEuhAbHuXXaeHGIze7cx55fV0fjOzbJQwiXYbwf8LLeTTDdan97u7qCPW+tjiy\niHQdtHvU9igmjirUCJEugzKP7PRh26Peo7U4Co5EO4h0DRaOvvBhXQB7PI5q0AiRLsK6R9mt\nmT5suxgRR6FZnfR7lQZEugKqPeq70nWsW/nnd3RjHDWASBdAl0edPP33toeH4qgajRDpEoQP\nwLIC9Q/b7mzHUdUaIdIV2Pao1E0UXU/sUPxejqPtYakGkapHk0fj9qZx5MrxXjftfHtpj24Q\nqXYCh6Ap5FE4jrxrSAtxVL9GiFQ96jyaxpHzaFjt7R7b2ebiHv0gUt2EPSp/AbbzaJhnLsfR\nfARp3qq0IFLVaPVoOq0bx5FdjaM071RqEKlqpkedKeXRyrTOzOLI2NnW8p5KQKSaUePRiCGO\nep/a/dNOX0cjRKoaXR5txNGsyjCzpmKNEKlmAgeiKe6RdVUGux5Hlzk56kCkalHlkXfDhvX+\ntxxHs9HU7REiVYtSj6wfR65D0y5fa1b3BJEqRZVHTqOhwjC6MNz3cLJhl3ZUCCJVyuTYK+nR\neFpn/WndoJFdnNVdQSNEqhV9Hh2Oo8TvVCYQqUr0eDSLo766MPJo6Ois59cAkWpk6pEt9geQ\n5nFk1uJoOo6raIRIVaLHI+N71PVl7pHr5zxHrwMi1Ycaj/xpnRk8CsfRVU+OOhCpPszkUeE8\n8uLIPe6++5298KzuCSJVhxKPFqoMcXF0NY0QqT6CHuU3yV/LYIY8ar/ZcRxdXyNEqg4lHjmN\nIqZ1d9DolEjvoh0JNgFT5h5Zk9+jQJVh0GiSQ9euMfScEKl5+RLtSqAJmKDDo8GnIXXMnePI\nnhLprWn+iPZl3gSMUeVRRBzdRqNz50ifTfP6LdmZeRMwwowfFPZoO44mvb+uRieLDd+vTfMp\n2JlAE+ChyKO22O06EhdHGd6hcpys2v35F0o/Yp0JNgEOTR711jitfI/a/ZPOX9uj0+Xvn4+m\nQ6pHsyagRZVHoTjyyws300jgOtIfRMrD3KNiK77tzji6vkZM7apBkUetPdZ6C71HcXRDjSg2\n1IImj6ZxZCdxtNLz60L5uw7MaLucR+M4Ghe97xtHlguylTD2yJpyHq3H0V01YolQHWjzyPQe\neWVwO53G3UgjFq1WgRKP+mldf+ORXwafx1Gut0cF3Eahn4BHJViJo+GmpKGj9/IIkfTjHZGm\noEetPeE4MnfXCJH0Y8abrUcF7kDqL7WaUByNu3w7jRBJP2a0VcQjl0LjOLILcXRHjRBJPQo8\namWy0zjyF9sFO3wrEEk35T2axZGdxNG4vzf1CJF0U9Ijr51nB7qt7qmuU8zqOhBJM2XzyDkz\nxJGZxBEa9SCSYkIeZaWPnO7hNI4WentHEEkvxT3q7HE3EBJHKyCSXoy/kd+jIY78G48Gj0Zd\nvblGiKSY4h6ZUBy5md6op7fXCJH0UtijTqOIOEKjXxBJKWU9ctO6mDjK+75oBZF0MvLIGndE\nZ/PIuEuvZh5Ho37iUQsiqaSsR06jQBxZZnVhEEklxv+e2aN+WrcVR2jkgUgaKevRJI68m8g5\nOVoEkRRS0KPx2dEg07MnxNEyiKSPwh7ZcByh0SqIpI6ZR3npNRrFkfHVQaM5iKSNch5Z98mT\ng0c2GEd535A6QCRllPSIODoOIinDeN+yfzZDG0furqeuG2gUAyLpopRH4zjytBp5hEaLIJIq\nJh5lM8kG48igUTSIpIlSHpmFODIjj0q8I9WASIoo5FE3rfPmkraLI0scxYJIeijn0WYcodEW\niKSGcvM6Fz+DRu16ir5naLQJIqnBWxeazyPrZnXDnm6WN7mVA9ZBJC0U82gaR7bf2/YIj2JA\nJCX49ynkkcip5K0I6q4beXGERpEgkg4mHmVxaSmOBn3QKBpEUkExj9biCI12gEgaKOBRJ9Bo\nKUO/TsjrDsSBSArI7FHXhnWXj7w4ssTRQRBJAab/km1e17oziSPb+4NGe0Gk8mT1aCmOrBdH\naLQfRCpOdo/ajZU4KvRGVA0ilabIvM5fXDfrRJn3oXIQqTA5PfJX0y3EERodBJHKMvcofSZN\n4sjPJDQ6CiKVZbhuk8sjuxhHaHQCRCqKv4wgl0eTOOqLdWh0BkQqSXaP7DiOzDCtI47OgUgF\nyehRvyIoEEcGjc6DSOXI4ZH/kk4j94mqfRyh0XkQqRi553X+nbBunbf17z2CEyBSKWYeJVDH\n+9ZP4dwHfBNHoiBSKUwuj0zn0SiOXIEBjYRApEKk98izycWRW6bqVqiikRiIVIbcHg0auTll\n+wQIgUhFSO6R96KLcVT0HbgaiFSCTHk0iyP3kDgSB5EKkNUj43vUr6pDI2EQKT+5zo/MNI5s\nH0clh39NECk/JodHwThiVV0yECk7OT2aaGSIo1QgUm6yeNTpM3j0bNESR8lApMzk8GgljoqO\n/cogUl5y5VF36XUURwaP0oFIWcnnkdeIJY7Sg0g5yXZ+NC/WGTxKCiLlZPAooUKzOLIGjZKD\nSBnJ4VErT2tUrxQeJQeR8pHOo9FrWm/qaNt1FGiUHETKRns0J19dRxwVAZFykcyjcBwNy+rQ\nKAeIlAnjvibzaBJHZlAK0oNIefCviqbxyPhnR/1yIDTKBCJlIZVHk2ldt6ub1nENNiOIlIPk\nHg23j/cqEUdZQaQM5Mgj/8O2WBKUH0RKT7LzI18j630mgyGOsoNIyUnk0XRa5zTqp3mQE0RK\nTZY8GjZcEbz0sO8GIiXGDBW0BAaNPLJ2mN1BZhApLcY79ZcWaPKgX8dAHBUAkZIy9kg2lcZx\nNJweoVEJECklnkfTHJHyqNfIDoVvyA8iJSS1R/7fgyWOyoJI6Ujo0fjF+nUMxFExECkZiT3y\n4sg4j0oP+cYgUircx2El88g9cMvriKOSIFIaTDqPvA82McSRFhApCW7xWwqPTP+y7o5y4qg8\niJSCZB4NJ0ZD+Zs40gAiJSBhHo1fzxJHWkAkedJ51E/r3F8kJ460kF+kz5emefubtImyuMWj\n0h5NX84yrVNERpGa5z98bZ58JGlCA+1dE8NSbCmPBpu8OGJap4XcIn00Hz/Wfn80nymaUMDY\no2Q381mmdbrILdKj+fnd/mleUjRRnkQeeRM649ogjjSRW6Sm8R6IN1GcpHnUu2SJI3XkFund\nifRI0URp0ng0XaFq+1VBoIasIr39+fzb/O/f5s/HerWhUpFyeWS4YUIdWUVqeW4+flI0UZYM\n87quHeJIHTmvI319fX6+vT1LDh+rHtUpUsCj8zItTetKDhQCsLJBikQeea/RNsPZkUr0iNT4\npGkiJUk8MsPLWOJINXpEytyEMJ1BVtYjf6l32wxxpBREEqFfQJqmXkccqQeRJEjjkSGO6qFA\n+TviNKgykcb3TYieGpkhhIgjzWQU6fOqIiXxyPOpbQSPVJP1OtLjNXUTJUjsUdcI0zrdZD1H\n+tq4DUmgifwk9cgwrauEvMWGz+YrdRO5SetR1wZxpB6qdufI5VHRQcI2iHQK747vdB4RRxWA\nSGdw149MCo+6JoijKkCkE7h1QfIe2cEj4qgKEOk4xt1iJ+5RLw9xVAuIdBjjPjlB2iPiqEIQ\n6ShTjwQ/MahvAY2qAZEOMvZIdMW3a4E4qghEOkZyj9CoLhDpEM4jOxz/ohrhUW0g0hF8j047\n5H+KKnFUK4h0gIBHx3UaZod4VDGItB9Jj8ZnWe3L41GFINJuRD0aCdW+PBrVCCLtxfPIStXq\nhjjCo0pBpJ0Yt4xU0KPh2isa1Qoi7cP9ORUj6ZHBo+pBpF1IeuTOr/objtCoYhBpB0bQI2+x\nK3F0ARApHmMFPeq/9x8bjkY1g0jRpPDI4NFFQKRYRh6d08i/fOT9GQuoGESKxPfoLNNVQQaP\nqgeR4pD0aBCKad1lQKQo5D2yblqHRpcAkWLoL/actsfbxqMrgUgRTG+HFcAOa4xKjw4kQKRt\n3EVToRthuxAiji4FIm0yLGfAI1gCkbaQ9cgOHhk8uhCItIGwR/0kkTi6Foi0jqhH42WveHQl\nEAmG/ycAAAv+SURBVGmVdB7xKarXApHWmHl0RqfBI0MeXQ5EWkHWIzN4RBxdDkRaxs+P03Qe\nGTy6Joi0iPNI8s9MOJPgaiDSEuk8Kj0ySAAiLZDEo/7OCbgaiBTGmN4jicK3exWmdVcFkYJ0\ny3iE8sh5RBxdF0QKIemRHb7g0XVBpAByHnUWuQ8eKjssSAgizRGd11mqDLcAkWbgEewHkaa4\nD8gSwfZnR3h0bRBpgvs4BRGNLFWGu4BIY9wtd4IeWTy6AYg0YrgMe9Il63nEtO4GIJKPvEeW\nOLoHiOQh6JHtl3nj0S1ApAGv7C3gkTvfKjIUyA0i9Qh7ZKjW3QlEcnQe9dd+zpiER7cDkTr6\ny7AiHrUvhEf3AZFavMuwJz2yg5BwHxDpibcsSMIj7ii/HYj0y+DR+WmdHbbgPiCSlfXIfegW\n3AtEwiMQAJGkPOrWReDRPbm9SO60RuDWCaoMN+buIkl7xEXYm3JzkZ5HvZG4/2j4bG+4I/cW\nCY9AiFuLJHhXef/J3nBP7iyS7KectN/gptxYJPeXXIU0snh0Z+4rkqhHFo9uzm1FEvLIv4Ep\nfadBLXcVafgzK2c96oVM3mdQzE1F6j06PbNjeR38ck+R3KEvddcEHt2eO4rUH/pCHvFZQXBH\nkfAI5LmfSINHZ7F4BI7bidQf+gIeOZnS9RZq4W4iOY9OW4RH4HMzkYy7jfU0eAQ+9xLJfar9\n6Tjqz4zwCJ7cSSSDR5CKG4mER5CO+4gk61H3kgn6CVVyF5GeCkiUvW0fR3gEAzcRScgji0cQ\n5hYimd6j02k0fN4WHoHHHUQaguScR52O3YviEXjcQCRjBD3qp3V4BD6XF6k9/M9/eN3oviM+\nlxgmXF0k425gPR9FeATLXFukzgEhjwwewRKXFsm4/8MjSMyVRcIjyMZ1RWoNwCPIwmVFMu6v\n50nhvSzAlIuKZAQ98j9GFY8gzDVF6v8KpYBEeAQRXFIktyRuyJNzHlk8gg0uKJLxPZLIIzyC\nTa4nUn/n3WmPDB5BLFcTqTvuxTyyeAQxXEyk7jYH59Hp+4+eX7tXxiNY5loiGakP3PJk6l9Z\naGxwSa4kkjudEfkzfN2Ge2k8glUuJJK7e1XGI7/MgEewxWVE6u8WklsWhEcQzVVE8j5rTsok\nPIJ4LiKS75EQ3qec4BFscQmRhmldAo8sHsE2VxBJOo7GyxnwCCKoX6TRR3Ebkb8h9vzaviQe\nQRTVi+R/hPCgwWmZhheXHx1ckNpF8j76VPxuWDyCaOoWyeAR6KBqkUZ/GQKPoCAVizT6RHsp\nh/AIDlGvSN0yb9trJPPZ3ngER6hWJG9aJ3LvER7BGSoVKUmVAY/gMHWK5E/r8AgUUKVIs2md\nhEbe3+DDI9hLhSJN1wSJfDQDHsEp6hPJm9YJLQfCIzhNbSINx7zkNVg8gpNUJtIsjmQ+vg6P\n4CR1iTS6CCt2VzkewWlqEsmf1kmWvfEITlORSOM4wiPQRDUimd6j59EvpZHvkcEjOEotIj01\nSuyRxSM4SiUieXEk/MneeAQSVCFSP60bDn48AlXUIFI/rROOIzwCMfSLNI8jPAJ1qBfJdF+M\n+1t8eAQK0S6Sm9H1S0vxCDSiW6T+7+/1K7TlPDIjjxAJTqFapO4I7/4osuxfmsAjkESxSMad\nEglfgx1/ujcegQR6RXLLGJ6ZJKkRHoE8akVy6xi6OBKd1+ERSKNUJHf91cURHoFudIrUx5Hg\nbbB4BAnRKVIz0khyYodHkASlIpmhZCd7FRaPIAk6RfLjCI+gApSKZLsLsaIa4REkQ6dIjXUr\nD/AIqkCpSO6IT+WRwSMQRadI7oBP5pH3FUAApSK1X/EIakGnSN2iVTyCWtArEh5BRWgVqT/4\n8QhqQKlIcg7hEeRAp0h4BJVxA5HwCNJzfZHwCDKgVKQ0903gEaRCp0hyxQY8gizoFAmPoDIu\nLhIeQR6yivTfn7fml7eP/9Z/EI+gMjKK9PPSDLyu/qicRHgEWcgo0kfz+N/Xc+v776P5WPtR\nIYfwCHKRUaRH89VvfzWPtR/FI6iMjCI1zdKDGRIeWTyCfFw0kfAI8pL3HOnv93Mr+TkSHkFm\ncpa/X72q3cvP2k8KemTwCDKQ9zrSx/M60uPtT/LrSBOB8AjScrmVDf3H7rcvZP1vAKnQI5I3\n7zuxZrV38PmaeASZ0COSDx5BZVxLJDyCQlxKJDyCUmRd2TBi7UfxCCojo0ifiUXCIyhHzqnd\n12P95omBY4mER1CMrOdIX+sLgwZOi4RHkJe8xYZPb93qGsclwiMowqWqdngEpbiUSHgEpbiS\nSHgExbiQSHgE5biOSHgEBaleJDt8f/7L7gXS9xHAo3aR8AhUULlIeAQ6qFukwaOnPHgEpaha\npLFHBo+gGDWLNPHI/dv03QOYUrFI1tvCIyhLvSLhESiiXpF6GtP8/n/78UMNQBH2H+UZRCrb\ndJZWaERfI5kPbUSikWs2gkgVtkIj+hpBpApboRF9jSBSha3QiL5GEKnCVmhEXyOIVGErNKKv\nEUSqsBUa0dcIIlXYCo3oawSRKmyFRvQ1gkgVtkIj+hpBpApboRF9jdxHJIDrgEgAAiASgACI\nBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCAAIgEIgEgAAiASgACIBCBAbpE+Hs3j\n42dtR4pWDn80+jqf41dMMpRJI0lG8vmS4ZcybSTJSH7em+b9y9+T5pcyJ7NIr88372VlR4pW\nvtKI9DV+xSRDmTSSZCQfz5d8DEdbipFMG0nzO3k8X9MzKc0vJUBekf5rHl/269H8t7gjSStf\nzZvk67tXfYwOgyRDmTaSYiRfzfvPb/C9ux0pRjJrJMnv5OP39T+8V07zSwmRV6SP5u+/r/9r\n/izuSNLKp+zruxd9HR3jSYYybSTFSN7aBoZ2Uoxk1kiS38mj+Rk1kuaXEiSvSG/Ntx39x2i2\nI0krn82n5Ou3NB92dIwnGcq0kSQj6Zrq20kykmkjKUfy6DcTDmXaaPIWRq1N/7s025Gklbfm\n7/u/U07JNv79dibdTjKUaSNJRvLkp3l1m0lGMm0k3Ug+PEXTDWXKTUR68rr4L841tNBqikZS\njeQ3If5OWkxw9HmNpBrJ/5rGsxORhFv537//Gn6ITybyi5RqJPb7Mcx+kh1940bSjOTz7eGd\nESGSaCstP+JV0PwitciP5OfhRUOqkYwa6XalqEy/D3ZeVaTHdFyzHUla6Uh7jCcZysILijfy\n6h/PqUbyGpAmxSH+M1QbUg1lTomq3fe0avedomo3f9HUpy8JhjJtZGXfCb5fXr+9h2lGMmmk\nI8khPi1Ayv9SAm0mb8Hnz/Nk8+9wOjjbkaSV9vJC4mM8yVCmjSQZyd/JKX+SkUwbSTIS96J9\n9KX6pcy5xcqGj9838udjKBkJkWFlw6SRFCP5npbOUoxk1kiS38lzZcPP23COdNWVDfZlqHm2\nB8hLiiLotJWfdgmW+H+W3DGecCiTRlKM5L0Zlr0lG8mskTS/k0eW4ytEZpF+notx25abyY7E\nrbzIX0kfi5RkKKFGhEfSzERKMJJwI/K/k+FFk/5S5mQWCeCaIBKAAIgEIAAiAQiASAACIBKA\nAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiA\nSAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIFKNvHd/\nzPH192+mggYQqUoez783/Nk8SncEOhCpSv5rmu/fv2ec/q91QxyIVCe/k7s3JnZ6QKRKeTR/\nmNgpApEq5d/kjomdIhCpVt6Z2GkCkWrlwcxOE4hUKe8NtQZNIFKd/PcvjzhJUgQi1cmj+R/X\nYzWBSFXyb2JnWSGkCUSqkf+a5ufft28md2pApBppl9qx2E4RiFQhbvE3kzs9IBKAAIgEIAAi\nAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKA\nAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQiASAACIBKAAIgEIAAiAQjw\nf72YmqLsykBBAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Using the data generating model from previous code chunk\n",
    "n <- 100\n",
    "x <- rexp(n,rate=2)\n",
    "b0 <- 10\n",
    "b1 <- 5\n",
    "plot(x,y=b0+b1*x,type='l',col='blue',ylab='Y',xlab='X')\n",
    "poss_x <- rnorm(n, sd=1.5)\n",
    "output <- matrix(NA, nrow=num_sim, ncol=n)\n",
    "\n",
    "num_sim <- 1000\n",
    "for(i in 1:num_sim){\n",
    "    error <- rnorm(n, sd=3)\n",
    "    y <- b0+b1*x+error\n",
    "    slr <- lm(y~x)\n",
    "    abline(a=slr$coefficients[1],b=slr$coefficients[2],col='grey')\n",
    "    output[i, ] <- slr$coefficients[1]+slr$coefficients[2]*poss_x\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///89ODILAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAcXUlEQVR4nO3di0LbuKKGUQUoUApM3v9ph9wDpfwQS7bsrHX2oem0VLl9\n2JYdu6yBwcrUdwCWQEhQgZCgAiFBBUKCCoQEFQgJKhASVCAkqEBIUIGQoAIhQQVCggqEBBUI\nCSoQElQgJKhASFCBkKACIUEFQoIKhAQVCAkqEBJUICSoQEhQgZCgAiFBBUKCCoQEFQgJKhAS\nVCAkqEBIUIGQoAIhQQVCggqEBBUICSoQElQgJKhASFCBkKACIUEFQoIKhAQVCAkqEBJUICSo\nQEhQgZCgAiFBBUKCCoQEFQgJKhASVCAkqEBIUIGQoAIhQQVCggqEBBUICSoQElQgJKhASFCB\nkKACIUEFQoIKhAQVCAkqEBJUICSoQEhQgZCgAiFBBUKCCoQEFQgJKhASVCAkqEBIUIGQoAIh\nQQVCggqEBBUICSoQElQgJKhASFCBkKACIUEFQoIKhAQVCAkqEBJUICSoYISQCj35b+o7MAcX\nvMvrhzPBEHzff1PfgRkQEpGQMiERCSkTEpGQMiERCSkTEpGQMiERCSkTEpGQsklCinuvhNQV\nIWVCIhJSNmJIPzikQkhdEVI2Ykh/VkKaJyFlY67avd6V25ftv/DZPzHwCEDaEVI27jbS71J+\nr20jzY2QspEnG15uy92rkGZGSNnos3YPZfUkpHkRUjb+9PfzTd4GElJXhJRNsR/pl5DmRUjZ\nvA4RMp83CSFlcwppW5GUxiekbFYhjTU87wkpm1FI5as/pCEhZUIiElImJCIhZTMKyTbSVISU\nzSoks3bTEFI2p5DsR5qIkLJ5hcQkhJTNLyRLpdEJKZtbSLaTJiCkbHYhjXUXOBFSNrOQ7Eua\ngpAyIREJKRMSkZCymYVkG2kKQspmF5JZu/EJKZtbSPYjTUBI2fxCYnRCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoRE\nJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyoREJKRMSERCyiYI6XFVbh7bDkFV\nQsrGDOn5rqwe1w9l47bNELQgpGzEkJ63Bd2XX6/rl7vy5TJJSF0RUjZiSL/K/Xp9X1ab26/l\npsUQNCGkbMSQyvYby93Zb97/8ZkLh6AJIWWjh/R7t063WzDVHoImhJSNumr3tnW087pdzas/\nBE0IKRsxpNfVcZWtfL1AElJfhJSNuh/p/pDP6svlkZA6I6TMkQ1EQsqERCSkTEhEQsqERCSk\nTEhEQsqERCSkTEhEQsqERCSkTEhEQsqERCSk7EpC8smMIYSUXUVI24qkdDEhZdcRUot/9IoI\nKbuGkMqHX/khIWVCIhJSJiQiIWXXEJJtpIGElF1HSGbtBhFSdhUh2Y80jJCyKwmJIYSUCYlI\nSJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlI\nSJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlI\nSJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlI\nSJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlI\nSJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISJmQiISUCYlI\nSJmQiISUCYlISJmQiISUCYlISJmQiISUCYlISNnQkB5v1uuXm3Lzp9Yd+nsIpiakbGBIT+Xt\nd6vypmpJQuqKkLKBId2W3+vncrP+XW6r3aW1kDojpGxgSJsF0nO5392oR0hdEVJWIaS78iSk\nRRNSNnjV7vmprNZW7RZNSNnwyYZSHjYLpKdqd2ktpM4IKRs8/b3abCGtb35Xuj+fDMHUhJTZ\nIUskpGxgSHf31e7Jv4ZgckLKKszaNSCkrggpGxjSTXmtdlf+McRYSqMfCgsgpGxgSK93t3WP\nsvt7iHFsK5LS54SUDV61O6p2l9aThDTRuLMgpExI7wZU0meElJn+fjegkD4jpExI7wYU0meE\nlNUK6c/d0HsSh2jKNtIXhJQNDel+GdtIZu2+IqRsYEinjmZ/0Kr9SP8kpGxgSKvye31bXl5u\nfdR8wYSUVThE6OFtafTs80gLJqSsQkhP5dEnZBdNSNnQo7/fVu1eys36j5AWTEhZjdNx3W4m\nG35Vu0trIXVGSNnQ6e+Hze9+lVL3c0lC6oqQMkc2EAkpG/oximp35J9DMDkhZUNn7W6r7oj9\nbAgmJ6Rs8Cdk3zaP6n+0T0hdEVI2dBvp5eGtpZuHyqt4QuqKkLIKkw0v96tSeRVPSF0RUlZn\n1u5x9kd/8wUhZTWWSNu1u6qnWhVSV4SUVdlGWt2/1Lo/nwzB1ISUVZi1+2XWbuGElA3ej1T3\n7PmfDMHkhJTVPLKh3oSDkLoipKzmsXZCWighZUIiElImJCIhZUIiElImJCIhZUIiElImJCIh\nZUIiElImJCIhZUIiEtLeF58WchYhIiFtfXnFkgEhlfcuvXuV7hXtCGnry2toCYlISBvlw6+f\n/+kF/+DW3WpzsoY/q6pnLBZSX4S00TSk+/K8/fW57jmLhdQVIW00DamUjze+8PrrdLahr/++\nkLoipK1W20gbq+MSaRW/73W13ZbaXbVZSDMipK1Ws3Yb92W1OWPD06o8xO+731yQ7PVxdXu8\nU98agskJaa/hfqTb/ZzdXf6+1e4bX1Y3L0KaFSFlg3fI/r7bZPSd86we2nm9vf0spGZz6Qwl\npGzEIxtuyuFUKTe3lkhzIqRsxJAej5fHfCm3QpoRIWWDQ3q62zRx951Trd4f63kKa29C6oqQ\nsiqTDW//bfWdkp6PUxIvv4Q0H0LKBob0WG5fNyE9uqr5ggkpG7xD9vXr/VSXEVJXhJRVOERI\nSEsnpGxgSDf7JdJzual2l9ZC6oyQsjrbSE+rzeE/9QipK0LKBn8eaX8owm2tO/T3EExNSFmV\n/UjlrvJVkoTUFSFlTn5CdOUhfevYz4Eh3VX9YOynQzC5qw7pm5PStT4hW5eQunLdIZ19jX/t\ngn9553REd1VC6so1h/T1mRr+/nsX/NMbr3e39a9pLqTOCKl5SM5rdw2EJCQquOaQRtpGakRI\nXbnukEaZtWtESF256pDG2Y+0P/nJr++c/OTiIZjYlYf0LSOejuvSIZiakLLhJ4jcLIwc/b1o\nQsqqnbLY55GWS0hZrUOETH8vmJCywat2hyVS1Y0kIXVFSNnQyYaH/YXGfLBvwYSU1TuyoebR\nDULqypWG9KM3tJCIrjKkH54cy5ENRNcZ0tnXb//1C0b47A8skZbpGkP67lHfH//+BUP8/QdC\nWiYhff/vXzDE338gpGUS0vf//gVD/P0HQlqmawzJNhLVXWdI083aCWmhrjKkcfcjfRj45//Y\nD4dgClca0o8IiUhImZCIhJQJiUhImZCIhJQNPWXxw3euZj5oCCYnpGz40d8tWhJSV4SUDT33\n9+9fLVoSUleElFXYRvrzcFO7JSF1RUhZncmG59XbcqniCbmE1BUhZVVCerqtfEFmIXVFSNnw\nkF4f3hZHN0+vbzVVO5OQkLoipGxoSH82kw33u3Ny2Y+0UELKhu5HelsYPR6ufllWNe7RxyGY\nnJCyofuR7ipfh+LvIZickLKh+5Gq3ZF/DsHkhJQ5HReRkDIhEQkpExKRkDIhEQkpExKRkDIh\nEQkpExKRkDIhEQkpExKRkDIhEQkpExKRkDIhDVLxcp8dE1ImpAF+eMGC2RJSJqQBfngJndkS\nUiaky/30om6zJaRMSJfb3styBZtJQsqEdLlyLVtJQsqENEDZ3tEyl7t7MSFlQhpgvzxa/Mqd\nkDIhDXGIaOElCSkT0hDlOOMw9T1pSkiZkAbZFVRmdI8vIaRMSIMcV+xmc48vIaRMSMMcZhrm\nc48vIKRMSANdw2FCQsqENJBZOzaENNjidyMJ6RuERCSkTEhEQsqERCSkTEhEQsqERCSkTEhE\nQsqERCSkTEhEQsqERCSkTEhEQsqERCSkTEhEQsqERCSkTEhEQsqERCSkTEhEQsqEVNciP+Un\npExINS30c+dCyoRU00LPhCKkTEgVlQ+/LoWQMiFVJKTrJaSKhHS9hFSTbaSrJaSazNpdLSHV\nZT/SlRISkZAyIREJKRMSkZAyIREJKRMSkZAyIREJKRMSkZAyIREJKRMSkZAyIREJKRMSkZAy\nIREJKRMSkZAyIREJKRNSSwv5UIWQMiG1s5iP+QkpE1I7i/nguZAyITWznFOhCCkTUjP7BdIC\n1u2ElAmpmcMW0vxTElImpHbK2f9mTUiZkNopu+XRevaPR0iZkFo6rNXN/PEIKRNSSwuZuBNS\nJqSmljFzJ6RMSE0dZ+1mnZKQMiE1tsmolHkvk4SUCam1cjD1HbmckDIhtbZvSEjLJqTW9ptI\nc14kCSkTUmvbbaR5z9wJKRNSa+W4SJr6nlxMSJmQmtut1s25JCFlQmpu/htJQsqE1F4pM58C\nF1ImpBEcK5ppSULKhDSC/TbSeq4PTEiZkEZx/DzFLB/ZTEMadVVaSKM4vKRCGs3Ix5MIaRyH\nD/jN8pHNM6Szr6MN1/pbOhxiZPufj/OcuJtjSGN/plJII5nzHLiQvj9e22/pcIjxzfcwcCF9\nf7y237L7vvdaDNG32c7czTGkBW8jPX4d0rcrm63ZngplniEtd9bueXXbeoiuCWlcy92P9Fzu\nWw/RM6t2CzbuZMNjeW49RM9MNiyXWbsRzXUbUEiZkEY1zw/5CSkT0thmuDdJSJmQRlbG3sFR\ngZAyIY1sjseBCykT0sh2Ac1rK0lImZDGtvvE+axKElImpLEdT6o/n5KElAlpdMcZ8Nk8SiFl\nQhrf6ZxCM3mcQsqENL5tRGshLYqQJlAOMw4z2U4SUiakCZxOYiykpRDSFI6zdvMoSUiZkKax\nn26Yx0MVUiakaZyt3U19VzIhZUKayn6v7BzmG4SUCWkq+4rmsEgSUiakycznYi9CyoQ0meOs\nnZAWQEjTOazWdX9qfSFlQprQoaPeH6+QMiFN6Hj2hs4fr5AyIU1qHhfEFFImpEnN47JJQsqE\nNLFDRT2nJKRMSNPr/vxcQsqENLny4df+CCkT0uSEtARCmpyQlkBI07ONdKmOZjuFNL3DrF1H\nb4v3Og2pq7MDCqkHx3dEL2+L9/oM6bAPbur7sSOkXnS8gtdpSPuvfTxlQupEz1MOXYbU2WU9\nhNQJIf2QkPoYojdC+qHOrtgmpF509bZ4r8uQOruGqJB60dXb4r1eQ+rpKlNC6kc/74oP+gyp\nrydMSES9htQTIREJKRNSh3paZdnoLKTenp4tIXWnv1mHrkLq7+nZElJ3yrqn2aiNvkI6+9oR\nIfWmw8NXewqp1x3XQurNbo99V4skIWVC6s3hbOAdpSSkTEjdKfv/CelTtpH6GqJfhytU9PJJ\nm95C6m0LckdI/TlcNambd0tXIXW10nsipP6czq3fyTums5C6JKQedXbNcyFlQurR6UJ+XTwT\nfYTUy/L5c0Lq0+Fd08Uz0UNInc4xHAmpT13tLukipLOvPRJSp/ZvnC7WZzoIqasfLJ8RUqdO\nqzLTpySkTEjdKt1sJwkpE1LHenn3dBCSbaSL9Pt8jUpIJ52s5f6TkDp2DGniN1APIU3/LHxN\nSD3bPg/Tn7+tj5D6JqSenV25ZMqShJQJqW/bKyftbk13J4SUCal3HVx1QUiZkHonpFkQUvcO\na3ZC6pmQujf9tVKFlAmpf5NfvVtImZDmYOJ9kULKhEQkpExI8zLJMQ5CyoQ0J4fzdI08rJAy\nIc1JKaevIxJSJqQZmWrfrJAyIc3I4ROzQuqPkGbkcBj42DMOQsqENCe786/uz8I6HiFlQpqT\n/TEO2wXTiMMKKRPSvExyaqHxQur74+RfEdLMTHFClLFCmvyjwAMIaW4m+LzsaCGdfZ0bIc3N\nBD+2Rwqpl7OPXURI87PQ/UjHhzXLV19IRGOEdPahq1m++kIiGiGkiSb26xESUfuQjh+nn/aj\nwAMIafbabzKNFNK62I9U2VyfzQmMcehd85DKVIfj1iOkmRvjhMYjhHQ8zflcCWnedu/Axj/J\n24fUwQnOBxLSvO0+e944pTG2kSa/5MZAQpq3w0cqmr4LR5lsmHdHQpq707pduzfiSDtkRxil\nHSHN3OnEQu3ODu7zSJmQ5q7sbJ+0Rs+bkDIhzd/pyBohTUZI81cOU3dzCmnu20QfCWkJStvD\npuuHNPe9Rn8T0iKcvzHr/6xvENLZ12UQ0kKcDlerv0emekiz/izs54S0MC2OtRFSJqRlKS3W\nmmqGdJip3/2m4j88MSEtS5Pz7NcL6XjS5d3vqv270xPSsvQd0nF20azdZd/S4RBLdfhZX/PT\nFZVCere/y34kIXXtePaDivN3tUI6Hn2xwNdXSEtzPK/V4TM+w5/MOiEdJhnaHRM4JSEtz+HY\nu3I4x9XQp7NaSPv/X+LLK6Sl2p/bqlRYKFUMaYFbRztCWqpyKmno/Fi1baT1UjMS0oLtAirr\nSUM672Z5c95nhLRYx7mGwZ9Cvzikj+ksdXG0FtKSHTeQhi6SLg/p7OvCCWnBjrN2x2NyLgvq\n0pAWeEjdPwlp0Y4z4LvfrC9r6YKQlnps6j8JaeHeX3fosv1KPw7pdGqj9X7gxRPS4p1fB71c\ndAa8H4Z0PEDJNlLtb+lwiOtT9v8r658ulX4c0mG4RU94vyekq3GsaP3TibyfhXRaEC32MIZP\nCOl6lP3/7T9Z9/1DWgeEdDWEdD0OMw2H07J++5DWn4Z0OBjosrs5T0K6JuU0D3B2kHj8tq9D\n+mu5dvoUxxUR0pU5WyZ9ew3sq5A+Wa5dYUZCukLH3bN/hfSvd/8/QzoeNPGhxqvLSEhX6d0h\nDscGDv/1748wfQzp7EiJT3K8TkK6SuW4v/RUwGE5dZbS/ut/5ewPTod0n3+7kEb5lg6H4MPW\nzTGI0/8Oi5z/zv7uh0XQ6fCjK3/NhHTF3q3Cnboop562f2kX0tmc3/6b99/23bm/ZRMSO3+F\ndDqZ43+HVbp3i6TDL/VO+zVnQmJvv2T5dkjrU0lXn5GQODrO2p02k9afhPTuMxnXdvzCvwmJ\no/LR7j9+3EY69WNRdCQk3ns3+f35rJ1+/iYk/u3T/Uh8RkhELa5qvjRCIhJSJiQiIWVCIhJS\nJiQiIWVCIhJSNmpIfx7utjOod/d/Wg1BA0LKRgzp9eZsn/ltkyFoQkjZiCHdl9Xv5+2tl6dV\nuW8xBE0IKRsxpFV5Pt5+LqsWQ9CEkLIRQ3p3dMnfh5r8dbQkvRBSZolEJKRs3G2kp5ftLdtI\n8yKkbMzp79uzdbeb1yZD0IKQsnH3I91v9yOt7h7sR5oTIWWObCASUiYkIiFlQiISUtZpSPTk\nv6nvwBxc8C6vH86FRrknyxlkQQ9lEc+XkGY6yIIeyiKeLyHNdJAFPZRFPF9CmukgC3ooi3i+\nhDTTQRb0UBbxfAlppoMs6KEs4vkS0kwHWdBDWcTzJaSZDrKgh7KI50tIMx1kQQ9lEc+XkGY6\nyIIeyiKeLyHNdJAFPZRFPF9CmukgC3ooi3i++gkJZkxIUIGQoAIhQQVCggqEBBUICSoQElQg\nJKhASFCBkKACIUEFQoIKhAQVCAkqEBJU0E9Ir79K+fWc/94wjzdldf/lBQbrDNP4eb1fLeJh\nrEd6QUZ4b/UT0mp7GYDGJd1vB1m1fuGeG1/KfXfd0ZumY6zbP4z1WC/ICO+tbkK6L782X+6a\nDvJcfr1ufs7+ajrK+nnV9h34p6yeN4OEa44O1fphrMd6QcZ4b3UT0qpsfio1fuXudv9841Ee\ny23bEe7L09vX3+Wh5SDtH8Z6rBdkjPdWNyHtlNUoo7R91OW+8Qh35WW9+WnedvHd/GGcDTXG\nOG3fW32FdF8eRxjltdw2/fefm//0G+XnePOHcdT6Bdlq/N7qKaTf5e2H4Aget2tGTS0hpFFG\n2BrhBWn+3uoppMe7VePV/q2XVdtVog0h/cAYL0jz91ZPIb351X7d7nU1wnqEkL5vlBdk3fq9\nNXlI768i/dpmi/B8kNtWu1/OB2n7DlwtKqRmL8gHjd5be52F1OilOw3ycnP70mKE9Zgh7Wbt\nXhrP2q1HCanhC/JR00czeUgHu7n+l8a765/GmB9at34HPmw3zp/aT820D2mUF2SM91Y3IW33\nPr/etd1Gehmpo8bvwJGObBghpHFekDHeW92EtD8equ3z+quUD6uSjTQe4WaE52qj+RM10gsy\nwnurn5A2RzTfNJ6zKwsJ6XV79HfTIbaaP1FjvSDt31sdhQTzJSSoQEhQgZCgAiFBBUKCCoQE\nFQgJKhASVCAkqEBIUIGQoAIhQQVCggqEBBUICSoQElQgJKhASFCBkKACIUEFQoIKhAQVCAkq\nEBJUICSoQEhQgZCgAiFBBUKCCoQEFQgJKhASVCAkqEBIUIGQoAIhQQVC6lkp6/vj5csfbw5X\n5n66LeX26f2to9vy5+3rn/Jr1Ht69YTUs1Ieypvbze3bcrj5uL1VHs9vnbyU1dvX1ep1ijt8\nvYTUs1JWz+vnVfm9Xv8+3VyV583vb85vnXksD+uHzd9jRELqWSmb1bancrde3+1v3h7+6+nP\nP7otj5vvYExC6lkpx1/Obt6Xcvf8vPnd6da5l7e1vZcx7yZC6tvnIa0fVm+prF7e3Tp3X+5H\nvZsIqW//COltFe/+Zr9ldLp1ZIk0ASH1rGynsp82U9l3p82lw5/9fWvn7m0b6Xasu8iOkHp2\nmLV7ejdrd7Obxbs5v3Xm99uK3cP7KXGaE1LPStnuPdouhU77kX7v9h79Ob918rra7keycjcu\nIfXsbaXt7nA4w/px9e7Ihj/vbx392h/ZYOVuVELq2cetH7rlleqZkGbDK9Wzb4dUTpreIf7F\n894zIc2G5x0qEBJUICSoQEhQgZCgAiFBBUKCCoQEFQgJKhASVCAkqEBIUIGQoAIhQQVCggqE\nBBUICSoQElQgJKhASFCBkKACIUEFQoIKhAQVCAkqEBJUICSo4H8Pxc/2X445wAAAAABJRU5E\nrkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_pred_vars <- apply(output, 2, var)\n",
    "plot(poss_x, y_pred_vars)\n",
    "abline(v=mean(x),col='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare the confidence band resulting from bootstrap with the $(100-\\alpha)%$ confidence interval for the predicted y $\\hat y_{pred}$, which could be derived from SLR.\\\\\n",
    "The 95% CI is:\n",
    "$$\n",
    "(\\hat y_i + SE(\\hat y_i)\\times t(\\frac{\\alpha}{2};n-2), \\hat y_i+SE(\\hat y_i)\\times t(1-\\frac{\\alpha}{2};n-2)),\n",
    "$$\n",
    "where $MSE = \\frac{\\sum_i (y_i-\\hat y_i)^2}{n-2}$ and $t(q;n-2)$ is defined as $P(Z \\leq t(q;n-2))=q$ for Z~t(n-2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAANlBMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD/AAD///9qpps6AAAACXBIWXMA\nABJ0AAASdAHeZh94AAAfRklEQVR4nO3d4ULiSgxA4SkgchcBef+XvdKKggKWNskkmfP92GXv\n3RJQzsJMUcsRwGyl9g0AMiAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBA\nSIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBA\nSIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBA\nSIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBA\nSIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBA\nSIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBA\nSIAAQgIEEBIggJAAAYQECDAIqQDT/fdfhaETHuXy4VQYgbT++6/CUEJCMlU6IiQkU6cjQkIu\nlToiJGTyX62OCAmJ1OuIkJBHxY4ICWnU7IiQkEXVjggJSVTNiJCQROWOCAkp1O6IkJBB9Y4I\nCQnU74iQEF7d7bpPhITgXHRESAjOR0eEhNicdERICM1LR4SEyLxkREiIzE9HhIS4HHVESAir\nQkf3H5iEhKBcdURICMpXR4SEkGpsez98VBISAqrQ0R/flJiQEE+Njmb+f5lDHI5AXA47IiSE\n4215NPJvSBzicASictkRISEYnx0REmJx2tHxXel65yEk3Oa2I56REIjTjt7feWmHOPydhh2c\nXtYREqLwePro5H3035xw5fMQEn5x3REhIQinHb2/j/+7E65+JkLCD163GZ75yxOufyZCwjXv\nHRESInDfESEhAMfb3l8HTJjx/CEORyAOnx29X70riJDgndOOfhwyYcrzhzgcgSC8bntPOGb+\nIQ5HIIYgHRESXPPZ0fvvL5ogJDjmtKNphwkc4nAEAojTESHBL5+nYW9/LSwhwSuX2943lkfD\noROmPX+IwxHwzmdHd4+dMO75QxyOgHOxOiIk+BRnm2H00QKHOBwB16J1REhwyOW2971thrHH\nSxzicAT88tnR3CuQOOTSZlHKaqs6AqFF7MgypGFTZFl6a5URSMC+o1nbdefrmDD2+UOG404H\nrsv6cDzu12WjMQLxedz2frw8Gq5kwtznDxmOOx3YlcPp8qEsNEYgPJcdjbmWCYOfP2Q4rhy/\nb/XjW09IrYq4PBp5NRKHDMedDnw5h9RpjEBwYTuyDWn1utmWfx8XD+vHuw2E1CaHHY1YHo27\nIpFDhuM+9Re7g8YIhOaxI7FrEjlksNttNqtVv+WwftgRIbXI4+mj8T+Hj3c2wAePp4+e+HmW\nfkIql3RGwK/gHZmGdHgpZfn55iC2v3HJYUdjtxk+r27CLXj+kN6h659sVsOVEBK+OTwN++SP\nKTcMqX9b0GHTLfsrISR8Cb3NMPIKRQ7pdcOB+26xJyRcSNCR/bu/P56UlktCwjd/HT23PBp1\nlTKH9BblfPJosSQknDnsSOE6ZQ7pbcrL56V9WRISBjk6Mt3+Xn/Vs/3jVBEhtcLjtve0qzU5\n5NNudb60fyEkeOxowvJouF6TQxyOgAMOO5p8xSaHOByB+hJ1REioJsk2w7irljnE4QjUlqoj\nQkIl7jqaus0w6sqFDnE4AnX560j12oUOcTgCVSXaZvi8fpNDHI5ATek6IiTYc7ftPW95NEww\nOcThCFTjryOJESaHOByBWip09Ph/S3RESDCWsyNCgi1v294Cy6MRU6QOcTgCVbjryGaM1CEO\nR6CGtB0REgxl3K47TzI5xOEI2EvcESHBjLOOpLYZPmeZHOJwBKx560h4mMkhDkfAlre3Mwh3\nREgwkb0jQoIFZx3JLo+GeSaHOBwBQ85+2IR8RoQEA76ejlQ6IiSoa6EjQoI2Xx0pLI+GoSaH\nOBwBI846UptqcojDEbDRSEeEBFWtdERIUOTt9JHmYJNDHI6APl8daW0zfE42OcThCKgz7+jh\nj69TzYiQoMbV05F2R4QEJW11REjQ4aoj3eXRMN7kEIcjoMpXRxbzTQ5xOAKamuuIkKCgvY4I\nCfI8dWSwPBpugskhDkdAjauOzG6DySEOR0BJhdOw9/+fWUeEBFmNdkRIEOXqZZ1hR4QESa46\nMrwZhARJ7XZESJDTcEeEBDEtd0RIkOKoI9NthgEhQYajr+LTy+j+fSQkiHB0+qhGR4QEEa13\nREgQ4OjtDIrLo4f3kZAwm6dtBrWZf9xJQsJcbXT0x/8nJMxERyeEhHn8dFRredQjJMziqCO9\nmSPuIyFhDkfbdXpDx9xHQsIMbXRESNDVQkcjz5EREqbycxpWdZth3J0kJEzkqCO9mWM7sg3p\n7XVVTlbrN60RsNJGR6PvpGFIh0X5tlQZATONdDT6rxqGtC7dv11/ab/tylpjBKy46Ujza/ie\nuY+GIXVl93V5VzqNETDi5qv4NL8U9qk7aRjS1cfi4U+FIiTn6OgXnpHwNC8v63xsMwxs10jb\nfX+JNVJoXjrSXR49eSctt7+XF7t2i4PKCOhz05HizOf3UmzPI63780jd6pXzSFH52a5TnDnh\nPvLOBjyDju7wE1K5pDMCc9HRPZYhHV5KWW4/r4Tt74CsO7r3T6rqd1Kddh8t3yLUDW+0G66E\nkOJx05HizKn30XT7e/NR06br32ZHSPF4OQ3rsSPbE7L9b/tusSekgNpYHk29kxXeInRYLgkp\nHicdOVwe9QxDWpTzSdjFkpCi8dKR5sw599EwpE15+by0L0tCioWO/mC5/b3++uhs/zhVREjO\n0NFfTE/I7lbnS/sXQorDyWlY5eXRvPvo550NxiMwmpPTR8pPRzPvIyHhD3Q0BiHhMfuObv5n\nx8ujHiHhIR/LI/cdERIe8tGR6jaDzH0kJDzgpCPVmTL3kZBwHx2NRki4q4GOxLZSCAn3uNj2\n1l4eSd1HQsJtPk4fab+sE7uPhISb2uhI7roICbfQ0ZMICTe46Eh3eSS8BCQk/OajI92ZwneR\nkPCLi23vMNsMA0LCT3Q0ASHhBzqagpBwzcXyKNI2w4CQcMVFR7ozVe4iIeESHU1ESPjGtvdk\nhIQvHjpSXh6pPeUSEs48fHcG9Yy07iIh4RMdzUFIGHg4fRT0Zd0JIaHH8mgeQsKJh46UZ+re\nRULCkY7mIyTQkQBCQgMd6e9IElLzHJyG1d9mUL+LhNQ6B6ePor+sOyGkxtGRDEJqm3FHGbcZ\nBoTUtPzLI6u7SEgtc9CR9kyru0hIDcvfkd0rV0JqFx0JIqRmVe/IYHlkdxcJqVXVt73TLI96\nhNSm+qePcnVESG2iI2mE1KLqp2HVl0f2X/BrcojDES2r35H6TOuOCKlB9bfrtEdav3I9ElKD\n6EgDIbWmhY60J9xASI2pvV2Xb5thQEhtqd6R+sg6HRFSW+hICyG1hI7UEFI7ap8+slge1eqI\nkNpRvSP1kRU7IqRmtNCR+oj7CKkRtU8fJe+IkBpRuSP95VHljgipDZW36/Qzqt0RITWBjtQR\nUgPoSB8h5dfAdl31jggpv7od5d9mGBBScpVPHzXxsu6EkHKjIyP2IW0Wpay2qiNwRkdWDEMa\nPsjL0lurjMA1846u/miwPHLTkXlI67I+HI/7ddlojMCVutveBhn56cg8pK4cTpcPZaExApfS\nd+Rh2/vMOqTzR/vGT24TGIELdGTJOqSXc0idxgh8q7rN0NTyqGca0up1sy3/Pi4e1o93Gwhp\ntrodGYz01ZFtSIP+YnfQGIEzOjJmeR5pt9tsVqt+y2H9sCNCmouOrPHOhoQqnz4yGOmuI0ch\nlUs6IxpRtSObbQZ3HTkKyXhEXnU7MpjosSNCSqfq6aM2X9adEFIydFRHhe3vEcsgQpqq5nad\nxfLIa0eWIW0ISV3VjixGeu3I9jxSt9Qe0Tg6qsZ0jbT748uQBEY0jY7qsd1s2JSd9oiGVezI\naHnktyN27dKo+VXlrT8dHQkpjZqnYemIkLKgo8oIKYWap2Hp6ISQMsi+zeC/I0LKIPs2Q4CO\nCCmB7B253vY+I6Tw8ndkMWUuQoquXkcsjy4QUnD1tr1ZHl0ipNAqnj6ioyuEFBkduUFIgdXr\nyGZ5FKgjQgqs4jaDzchAHRFSXNk7CnH66AshRZW/I5MxUggpKDryhZBiqtYR2wy3EVJI9bbr\nbCaG64iQQqIjdwgpnoqnj2wmBuyIkOKp1pHZ8ihgR4QUju3j7HKbwWZiyIwIKZx623U2E4N2\nREjB0JFThBRKrY6MlkdxOyKkUKptMxhNjNsRIUVCR34RUhzJO4q57X1GSGHQkWeEFESt00d2\n2wyhOyKkIKp1ZDQxekeEFAMdeUdIEdCRe4QUQKVtBqvlUeht7zNC8q9WR1YTM3RESP4l7yjB\ny7qTGSGVa5VvVVp0FAIhOVenI8PlUY6OZr+0W3Xbj1/fuheh23NjRNPqbNcZLo+SdDQ3pHXZ\n9b/vylrm9vwe0bJK29509LyZIZXy84IIQurl78hqkr6ZIXVfz0idzO35PaJddToyWx6l6mj+\nS7vu7eO3bVdepW7RzxHNsu5o+J2XddPM3WxYfu7ZraRu0O8Rjaq0XWc1MFlH80/I/ludMtoK\n3ZybI5pER7Hwzgaf6CgYQnIp/TZDto7mh7RdnT4Nq73Q7bk1oj11OjKbmC8joc2Gj//WiZbU\neEh0FM/MkDZleTh9IjZF9D1CTYdU6fSR2cSUHc0/IXsYVqq8s0FKldNHLI/mEniLECFJqtOR\n2cCsHc0NafH5jLQrC7GbdGw5JDoKSmaNtO3KRuwmHRsOiY6imv31SJ9vEVpK3aDfIxpS4zSs\n3fIo6TbDQOQ8Uln9E7o5N0c0o8Z2nV1GqTvinQ2O0FFgM0NaiX5h7M0RzcjdUeLlUU/qK2Rl\ntRgSHYUmsP2toL2QamzXmW4zJO9obkiH1fJN7LbcHtGEKh3ZDczf0fyXdnxfOwHpO7KbVYtp\nSG+vw2mn1fqPp7HGQqKj+Ay3vw+Li+wen8BtK6QKHRkuj9royDKkden+Dd+8a7/tHn9DyaZC\nMu1oeOnA8kiczDc/eRnzzU/O3wPv5I/vg9dSSBW2velInuG34yrl3h8kblVUdJTD/G8QeXoy\nGvXub56RbrDvyHZ51EpHct+y+O+vRzpFN3xnB9ZIZxU6MhzYTka230R/ebFrt3j4johGQqqx\nXWc4sKWO5H6sy5jvWfy27s8jdatXziMd6SiVuZsNr58/aIwv7HsaHWUi986GuW8TUvs5mk7Z\nnz5im0GRn5Bm3qpoKnRkN6+9jvgK2UroKBnJkHhGGi17R4bDnDAM6YmXgdlDMt9msFweNdmR\nZUgbQvpk35HhvDY7Mn1ptxu9SZ47JDpKyHSNtHv8xqARIzJI3VGD2wwD282GzcX7VieNSMC6\nI+PlUaMdsWtnzHy7zvhlXasdEZKt7B1ZTvOFkCzRUVqEZMi4I7brDBGSHeuOTLcZGu9o9rcs\nfhX9aea3RqRh3pHdtKa3GQbz3/2t0VLGkMy3vQ3H0dHs7/3970WjpYQh0VFuAmukt9eFdEv5\nQjLermN5ZE1ms2HXfXzyBH8cc7qQrDuym3ako55ISNvliO/nPW9EbHSU3vyQDq8fT0eL7eGj\npjHfSWjKiOCsT8PSUQVzQ3o7bTash7eich7pJuOOrJdHdNSbex7p48loc/5ej4+/DfHUEdFZ\nd2Q37EhH3+aeR1qN+TkUs0YEZ7ztTUeVzD2PJHZD7o6ILXlHpuNc49txqaKjVhCSpszbDHR0\nhZAUZd5moKNrhKQnc0dsM/xASGroqCWEpMW0I/PlER39QEhKbDuym3VCR78Rko7cHdnOC4GQ\nNNiePqIjBwhJgelDzXh5REe3EZI825d1th8rlkd3EJK43MsjOrqNkKTRUZMISZjlQ81+eURH\n9xCSLNOOjD9MZPQAIYmio1YRkiQ6ahYhCaKjdhGSGMul+Dunj5whJCmmHcl957NR6OhPhCQk\n98s6OvoLIcmgo8YRkojEyyO2GUYhJAmJl0d0NA4hCaAjENJ8dARCms/usfZOR24R0lyGHRVO\nH7lFSDPREU4IaZbcyyM6Go+Q5sjdke284AhpBrYZcEZI0xl2JPjzecehoycR0mR0hG+ENBXb\ndbhASBOZPdbe6SgCQprGrqOjfUem45IgpElMOzIa9YmOJiGkKegIPxDSBHbLIzqKgpCeZvh0\nREdhENKzEnfEdt10hPQkOsIthPQctr1xEyE9xerB9k5HwRDSM8w6OtJRMIT0BNOObEad0dFM\nhDQeHeEuQhrN6MH2TkcREdJYVh0d2faOiJDGMXxZR0cREdIobNfhMUIaw3R5REcR2Ye0WZSy\n2qqOkGa6PKKjkAxDGh4iy9Jbq4zQQUf4k3VI67I+HI/7ddlojFCRd7uOjuRYh9SVw+nyoSw0\nRmjI2xHbdYKsQzo/Vn4/ZsqliSMU2DzY3ukoOOuQXs4hdRoj5Bl1dPqFjiIzDWn1utmWfx8X\nD+vHuw1uQsrckem49ExD+nrZVkp30Bghje06jGR5Hmm322xWq37LYf2wIychmS6P6Cg23tlw\nl+nLOjoKjpDuoSM8gZDuyLzNQEfyCOk2kwfbOx2lQUg32XTU/8rLuhQI6RY6wpMI6Qa7jthm\nyIKQfqMjPI2QfrF4tL3TUTKE9JPl8oiO0iCka6bbDKZ3k21vVYR0hY4wDSFdslsese2dDCFd\nMOlo+I2OkiGkb3Ydsc2QDiF9oSNMR0hnBo+2dzpKi5A+WXQ0/MZ2XUaENKAjzEJIPcvtOuOO\nLKc1jJCOxtvedJQSIVluM3D6KC1CMlwe0VFehJR2m4GOLDUfUtqO2K4z1XpI+o+2dzpqQeMh\nGXT0+Tsd5dZ2SJYdaU+6QkbWmg6JjiCl5ZDMlkd0lF/DIZk9HdFRA5oNye5lnXFHbDNU0WpI\ndARRjYak/mh7p6O2tBmSfkdfl+ioDU2GZNcR2wytaDEkOoK4BkOiI8hrLyTth1utbQY6qqq5\nkNQ7+rrENkNLWguJ7TqoaCwkOoKOpkKyWx6Zd2Q5DTe0FJLd8oiOmtNQSFm3GejIg3ZCoiMo\naiYk5YdbxeURHXnQSkjaHX1fpKMmNRJS3o4sp+G+NkKiIyhrIiS2vaGthZDMthnoqF35Q7J7\nWUdHDUsfUtblER35kj2krB2x7e1M8pB0H241l0d05EvukJQ7urjM8qhxqUOiI1jJHBIdwUzi\nkFQfb+90hEt5Q9Lt6PIPdIS8IZm9rDP9Eexs13mVNKS8yyM68ilnSHQEYylDyrvNQEdeZQyJ\nbQaYSxgSHcFevpDoCBWkC0nz8fZeryOWR85lC0m1o6s/0REuJAvJsCO9Qb/RkXu5Qkq7PKIj\n71KFlHN5xDZDBIlCSro8oqMQTEN6e12Vk9X6TX4EHaEiw5AOi/JtKT0iaUcsj4IwDGldun+7\n/tJ+25W17Ag6QlWGIXVl93V5VzrREYqPt3c6wt8MQ7p6CP5+PJZLT161ZkfXf6Qj3JTiGcms\nI04f4Q7bNdJ2318SXiOl7chyGuax3P5eXrx2WxzERuRcHtFRLLbnkdb9eaRu9Sp4Hinn8oiO\ngon+zoak2wx0FE3wkJJ2xDZDOLFD0nu8vdMRnhE6JMWOfvzZuCPDYZAROSTDjrQG3UJHEQUO\nyawjthnwp7gh0REcCRuS2gOu6jYDHUUVNSS9jn78me06jBE0pKzbdXQUVcyQsnZkOAyyIoZk\ntjyiI4wVMCSz5REdYbR4IeXcZqCj4MKFlHN5REfRRQsp6/KIjoILFlLW5REdRRcsJCV0hJkI\n6Vi/I8NhUEJIdAQBhPRrm4GO8LzmQ/qVER1hgtZDqv10REdJNB4SHUFG0yFVXx7RURoth1R9\neURHeTQcUvWODIdBW7sh0REEtRrS7+URHWGGRkP6nREdYY42Q7r1dGTYEdsM+TQZUvWXdXSU\nToMh1V8e0VE+7YXE8ggKmguJjqChtZDoCCoaC4mOoKOpkG5sM5h2xDZDXi2FdDMjOoKEhkKq\n/7KOjvJqJyQ6gqJWQqq9PGKbIblGQrqRER1BUBsh0RGUNRFS7Y5YHuXXQEi3lkd0BFn5Q7qV\nkXFHdrNQS/qQbj8d0RFkZQ+p9ss6OmpE8pDoCDZSh1R9m4GOmpE5JAfbDHTUisQhOejIbhYq\nyxsSHcFQ1pBuLo/oCFqShnQnIzqCkpwh0RGMpQyp+ss6OmpOwpA8LI/oqDX5QrqZER1BV7qQ\nPHRkNwteZAvpTkeTr+95dNSkZCHREepIFdLtbQY6gr5MIdVfHtFRsxKFREeoJ09I9Tti27th\nWUK6tzyiI5hIEtLtjDh9BCv2IW0Wpay2siPudfTk1cxBR20zDGl4eliW3lpyBB2hNuuQ1mV9\nOB7367IRG3FneURHMGQdUlcOp8uHspAa4WB5REewDun8+H78OH9iBB3BA+uQXs4hdTIjHHTE\ntjeMQ1q9brbl38fFw/rxbsPoES46spsFt0xDGvQXu4PAiLvbDHQEY5bnkXa7zWa16rcc1r87\nKpdGXd+djOgI9gK/s+FuR4I35C90hEHckOgIjtQI6e+XXn+PuLc8oiNUETQkH8sjOsJZzJB8\ndGQ3C+6FDImO4E3AkO4vj+gItcQL6V5GdISKwm1/3+9I44bcQUf4IVpIdASXgoVER/ApVEhs\nM8CrSCGxzQC3AoVER/ArTkh0BMeihHR3eURH8CBISHczoiO4ECOkBx3p3pArdIS7QoTkpCPD\nWYgmQEgsj+Cf/5BYHiEA9yHRESLwHhIdIQTnIdERYnAd0v1tBjqCL55Dup8RHcEZxyE96kj/\nJnwhI4zgNyQfHfF0hFG8hsTyCKE4DYnlEWLxGZKPl3V0hNF8hnR/BB3BpWAh0RF8ChaSITLC\nEwjpNp6O8BRCuomO8BxCuoWO8CRCuoGO8CxC+o2O8DRC+oWO8DxC+omMMAEh/UBHmIKQrvCy\nDtMQ0iU6wkSEdIGOMBUhfaMjTEZIX+gI0xHSGR1hBkL6REaYg5AGdIRZCOmEl3WYiZCOdIT5\nCImOIICQ6AgCCImOIMBpSHb+M5yFxCY8yuXDmcjkluQZkuiupPh4EVLQIYnuSoqPFyEFHZLo\nrqT4eBFS0CGJ7kqKjxchBR2S6K6k+HgRUtAhie5Kio8XIQUdkuiupPh4EVLQIYnuSoqPFyEF\nHZLorqT4eBFS0CGJ7kqKjxchBR2S6K6k+HgRUtAhie5Kio+Xn5CAwAgJEEBIgABCAgQQEiCA\nkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAgJ+QDi+lvOy0p2wWpVsftKcc\nN8of13WX4m4cjT4hBo8tPyF1/Y8BUC5p3Q/ptD9xuyk/z+AJy/5uLFRnHPXvxtHqE2Lw2HIT\n0rq8nH5ZqQ7ZlZfD6d/ZF9Upx12n+wh8K93uNORNc4j+3ThafUIsHltuQurK6V8l5c/carh6\n5SmbstSdsC7bj1//lVfNIfp342j1CbF4bLkJaVA6kym697qslSesyv54+tdc9+lb/W5cjLKY\no/vY8hXSumwMphzKUvX6d+r/+pn8O65+N75of0J6yo8tTyH9Kx//CBrY9K+MVGUIyWRCz+AT\nov7Y8hTSZtUpv+zv7Tvdl0QnhPQEi0+I+mPLU0gfXvRf2x06g9cRhDSeySfkqP3Yqh7S9U+R\nPuisCC+HLLVOv1wO0X0EdqlCUvuE/KD02PrkLCSlT933kP1iudeYcLQMadi12yvv2h1NQlL8\nhPykem+qh3Q27PXvlU/Xby32h47aj8DXfnG+1d+a0Q/J5BNi8dhyE1J/9vmw0l0j7Y06Un4E\nGr2zwSAkm0+IxWPLTUif74fS/bi+lPLjpaQS5QkLg4/VifoHyugTYvDY8hPS6R3NC+U9u5Ik\npEP/7m/VET31D5TVJ0T/seUoJCAuQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIg\ngJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIg\ngJAAAYQECCAkQAAhAQIICRBASIAAQgIEEBIggJAAAYQECCAkz0o5rr9+fPlmcf7J3NtlKcvt\n9aUvy/L28etbeTG9pc0jJM9KeS0flqfLy3K+uOkvlc3lpW/70n382nWHGje4XYTkWSnd7rjr\nyr/j8d/3xa7sTn9eXF66sCmvx9fT34MhQvKslNPLtm1ZHY+rz4vL83/9/v8/LcvmdAQsEZJn\npXz9dnFxXcpqtzv96fvSpf3Hq7295c0EIfl2O6Tja/eRSre/unRpXdamNxOE5NudkD5e4q0X\nnyuj70tfeEaqgJA8K/1W9va0lb36Xi6d/9/vS4PVxxppaXUTMSAkz867dturXbvFsIu3uLx0\n4d/HC7vX6y1xqCMkz0rpzx71z0Lf55H+DWeP3i4vfTt0/XkkXtzZIiTPPl60rc5vZzhuuqt3\nNrxdX/ry8vnOBl7cmSIkz36ufuAWnynPCCkMPlOejQ6pfFO9QbiHj7tnhBQGH3dAACEBAggJ\nEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJ\nEEBIgABCAgQQEiCAkAABhAQIICRAACEBAggJEPA/2SZWJjlMilcAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "error <- rnorm(n, sd=3)\n",
    "y <- b0+b1*x+error\n",
    "slr <- lm(y~x)\n",
    "y_pred <- slr$coefficients[1]+slr$coefficients[2]*poss_x\n",
    "MSE <- sum(slr$residuals^2)/(n-2)\n",
    "y_pred_sigma2 <- MSE*((n+1)/n+(poss_x-mean(x)^2)/sum((x-mean(x))^2))\n",
    "pred_lower <- y_pred+sqrt(y_pred_sigma2)*qt(0.025,n-2)\n",
    "pred_upper <- y_pred+sqrt(y_pred_sigma2)*qt(0.975,n-2)\n",
    "\n",
    "plot(poss_x,y_pred,type='l')\n",
    "abline(a=b0,b=b1,col='red')\n",
    "lines(poss_x, pred_lower,col='grey')\n",
    "lines(poss_x, pred_upper, col='grey')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}