ProblemSets_PartII/Math227C20Sp_P10_Bootstrap.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem set 10\n",
    "\n",
    "## Part I\n",
    "\n",
    "The concept that all our estimators have their own distribution, and therefore their own standard deviation, seems like a new concept but we are already familiar with one example:\n",
    "\n",
    "Consider $N_{\\rm obs}$ normal random variables with mean $\\mu$ and standard deviation $\\sigma$. According to textbooks, our estimate of the mean, $\\hat{\\mu}$, has a distribution with standard deviation $\\sigma/\\sqrt{N_{\\rm obs}}$. In reality, we do not know $\\mu$ or $\\sigma$ (if we did, why would we be estimating them!), so we approximate the standard deviation of $\\hat{\\mu}$ with $\\hat{\\sigma}/\\sqrt{N_{\\rm obs}}$. This is called the __standard error of the mean__.\n",
    "\n",
    "In this problem, we will explore whether the nonparametric bootstrap estimate of the standard deviation of $\\hat{\\mu}$ agrees with the __textbook formula of the standard error of the mean__.\n",
    "\n",
    "\n",
    "### 1.1 A single nonparametric bootstrap\n",
    "\n",
    "The following code generates a sample of $N_{\\rm obs}=20$ normal random variables.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>182.841279496126</li>\n",
       "\t<li>191.771695547986</li>\n",
       "\t<li>183.042142127269</li>\n",
       "\t<li>171.155066765154</li>\n",
       "\t<li>180.384972188731</li>\n",
       "\t<li>177.130487648841</li>\n",
       "\t<li>160.499777431114</li>\n",
       "\t<li>163.303138306464</li>\n",
       "\t<li>162.79187993073</li>\n",
       "\t<li>170.898101016939</li>\n",
       "\t<li>182.502058004729</li>\n",
       "\t<li>201.534283993282</li>\n",
       "\t<li>168.880009932178</li>\n",
       "\t<li>173.314032394432</li>\n",
       "\t<li>164.542897453057</li>\n",
       "\t<li>163.274456327346</li>\n",
       "\t<li>184.540365134401</li>\n",
       "\t<li>204.292940305087</li>\n",
       "\t<li>176.091428919027</li>\n",
       "\t<li>172.027995003157</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 182.841279496126\n",
       "\\item 191.771695547986\n",
       "\\item 183.042142127269\n",
       "\\item 171.155066765154\n",
       "\\item 180.384972188731\n",
       "\\item 177.130487648841\n",
       "\\item 160.499777431114\n",
       "\\item 163.303138306464\n",
       "\\item 162.79187993073\n",
       "\\item 170.898101016939\n",
       "\\item 182.502058004729\n",
       "\\item 201.534283993282\n",
       "\\item 168.880009932178\n",
       "\\item 173.314032394432\n",
       "\\item 164.542897453057\n",
       "\\item 163.274456327346\n",
       "\\item 184.540365134401\n",
       "\\item 204.292940305087\n",
       "\\item 176.091428919027\n",
       "\\item 172.027995003157\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 182.841279496126\n",
       "2. 191.771695547986\n",
       "3. 183.042142127269\n",
       "4. 171.155066765154\n",
       "5. 180.384972188731\n",
       "6. 177.130487648841\n",
       "7. 160.499777431114\n",
       "8. 163.303138306464\n",
       "9. 162.79187993073\n",
       "10. 170.898101016939\n",
       "11. 182.502058004729\n",
       "12. 201.534283993282\n",
       "13. 168.880009932178\n",
       "14. 173.314032394432\n",
       "15. 164.542897453057\n",
       "16. 163.274456327346\n",
       "17. 184.540365134401\n",
       "18. 204.292940305087\n",
       "19. 176.091428919027\n",
       "20. 172.027995003157\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] 182.8413 191.7717 183.0421 171.1551 180.3850 177.1305 160.4998 163.3031\n",
       " [9] 162.7919 170.8981 182.5021 201.5343 168.8800 173.3140 164.5429 163.2745\n",
       "[17] 184.5404 204.2929 176.0914 172.0280"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "NObs <- 20\n",
    "rnorm(NObs,178.1,12.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__(a)__ Perform a non-parametric boostrap, with 1000 resamplings, on the above data and estimate the mean $\\hat{\\mu}$ and the standard deviation $\\hat{\\sigma}$ each time. Output the standard deviation of $\\hat{\\mu}$. (In this block, the sample size is fixed at $N_{\\rm obs}=20$.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# the (sqrt)variance in the estimate of the mean -- nonparemetric bootstrap\n",
    "\n",
    "# the ground truth standard deviation of the normal distrubtion\n",
    "# the standard deviation of the sample\n",
    "\n",
    "# the (sqrt)variance in the estimate of the mean -- the textbook way, called \"SEM\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__(b)__ Repeat the above, but in a loop, with $N_{\\rm obs}$ ranging from 10 to 10000. (This is a big loop. If it is slow on your computational machine, feel free to skip, e.g., every tenth or 20th value of $N_{\\rm obs}$.) Compute the standard deviation of $\\hat{\\mu}$ from non-parametric bootstrap for each value of $N_{\\rm obs}$. Store the value of $\\hat{\\sigma}$ for each value of $N_{\\rm obs}$ (because you will need it below).\n",
    "\n",
    "Plot the standard deviation of $\\hat{\\mu}$ from non-parametric bootstrap versus $N_{\\rm obs}$.\n",
    "\n",
    "On the same axis, plot $\\hat{\\sigma}/\\sqrt{N_{\\rm obs}}$ versus $N_{\\rm obs}$.\n",
    "\n",
    "Do they agree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part II: Bootstrap on coefficient of linear regression\n",
    "\n",
    "__(Out of class)__\n",
    "\n",
    "__(EXAMPLE IS IDENTICAL TO IN-CLASS NOTEBOOK)__\n",
    "\n",
    "The code below simulates some linear data with noise. Linear regression using lm will yield estimates for the parameters beta0 and beta1 (and sigma, which we will not focus on in this problem). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Simulate some linear data with noise\n",
    "NObs <- 20\n",
    "beta0 <- 10.5\n",
    "beta1 <- 2.6\n",
    "sigma <- 1.4\n",
    "\n",
    "X <- seq(0,10,length=NObs)\n",
    "eps <- rnorm(N,0,sigma)\n",
    "\n",
    "Y <- beta0 + beta1*X +eps\n",
    "\n",
    "plot(X,Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__(a)__ In a loop, repeat the simulation and perform lm estimation, 1000 times. This will yield a distribution of estimates for beta_0 and beta_1. Plot histograms of these. Output the quantiles of these. (Note in reality one would not have this ability.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__(b)__ In reality, assume we only have 1 instance (in this example, with NObs=20 data points). \n",
    "\n",
    "Perform the simulation once to generate one data set with 20 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__(c)__ Write code to perform non-parametric bootstrap on this data set. For a linear model like this, the resampling should be random pairs of X,Y data (so, to be clear, the resampling does not mix an X-value with a Y-value it did not belong with). \n",
    "\n",
    "Using the distribution of estimates for beta_0 and beta_1 from this non-parametric bootstrap, plot histograms and output quantiles. Compare with the \"true\" distributions from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}