lecture03_demos (1).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# HIDDEN\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plots\n",
    "plots.style.use('fivethirtyeight')\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_units(x):\n",
    "    \"Convert any array of numbers to standard units.\"\n",
    "    return (x - np.average(x))/np.std(x,ddof=1)  \n",
    "\n",
    "def correlation(x, y):\n",
    "    x_in_standard_units = standard_units(x)\n",
    "    y_in_standard_units = standard_units(y)\n",
    "    if len(x)!=len(y):\n",
    "        raise ValueError('arrays are of different lengths')\n",
    "    return sum(x_in_standard_units * y_in_standard_units)/(len(x)-1)\n",
    "\n",
    "# This function generates a scatter plot with a correlation approximately r\n",
    "def r_scatter(r):\n",
    "    plots.figure(figsize=(5,5))\n",
    "    x = np.random.normal(0, 1, 1000)\n",
    "    z = np.random.normal(0, 1, 1000)\n",
    "    y = r*x + (np.sqrt(1-r**2))*z\n",
    "    plots.scatter(x, y)\n",
    "    plots.xlim(-4, 4)\n",
    "    plots.ylim(-4, 4)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the following functions use the analytical derivations for least squares\n",
    "# formulas can be found in the Lecture slides\n",
    "def slope(x, y):\n",
    "    if len(x)!=len(y):\n",
    "        raise ValueError('arrays are of different lengths')\n",
    "    return  correlation(x, y)* np.std(y,ddof=1)/np.std(x,ddof=1)\n",
    "\n",
    "def intercept(x, y):\n",
    "    b1 = slope(x, y)\n",
    "    return np.average(y) - b1 * np.average(x)\n",
    "\n",
    "def fitted_values(x, y):\n",
    "    \"\"\"Return an array of the regressions estimates at all the x values\"\"\"\n",
    "    b1 = slope(x, y)\n",
    "    b0 = intercept(x, y)\n",
    "    return b1*x + b0\n",
    "\n",
    "def residuals(x, y):\n",
    "    return y - fitted_values(x, y)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] File b'galton.csv' does not exist: b'galton.csv'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-302f19bb1953>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgalton_df\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"galton.csv\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mgalton_df\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhead\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/pandas/io/parsers.py\u001b[0m in \u001b[0;36mparser_f\u001b[0;34m(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, tupleize_cols, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision)\u001b[0m\n\u001b[1;32m    700\u001b[0m                     skip_blank_lines=skip_blank_lines)\n\u001b[1;32m    701\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 702\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    703\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    704\u001b[0m     \u001b[0mparser_f\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/pandas/io/parsers.py\u001b[0m in \u001b[0;36m_read\u001b[0;34m(filepath_or_buffer, kwds)\u001b[0m\n\u001b[1;32m    427\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    428\u001b[0m     \u001b[0;31m# Create the parser.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 429\u001b[0;31m     \u001b[0mparser\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTextFileReader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    431\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mchunksize\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0miterator\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/pandas/io/parsers.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, f, engine, **kwds)\u001b[0m\n\u001b[1;32m    893\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'has_index_names'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'has_index_names'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    894\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 895\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_make_engine\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mengine\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    896\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    897\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/pandas/io/parsers.py\u001b[0m in \u001b[0;36m_make_engine\u001b[0;34m(self, engine)\u001b[0m\n\u001b[1;32m   1120\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_make_engine\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mengine\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'c'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1121\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mengine\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'c'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1122\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mCParserWrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1123\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1124\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mengine\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'python'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/pandas/io/parsers.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, src, **kwds)\u001b[0m\n\u001b[1;32m   1851\u001b[0m         \u001b[0mkwds\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'usecols'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0musecols\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1852\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1853\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reader\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparsers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTextReader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1854\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munnamed_cols\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munnamed_cols\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1855\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/parsers.pyx\u001b[0m in \u001b[0;36mpandas._libs.parsers.TextReader.__cinit__\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/parsers.pyx\u001b[0m in \u001b[0;36mpandas._libs.parsers.TextReader._setup_parser_source\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] File b'galton.csv' does not exist: b'galton.csv'"
     ]
    }
   ],
   "source": [
    "galton_df=pd.read_csv(\"galton.csv\")\n",
    "\n",
    "galton_df.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a data frame with only two columns\n",
    "height_df=galton_df[['midparentHeight','childHeight']]\n",
    "height_df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[slope(height_df['midparentHeight'],height_df['childHeight']),\n",
    " intercept(height_df['midparentHeight'],height_df['childHeight'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_mse(x,y):\n",
    "    if len(x)!=len(y):\n",
    "        raise ValueError('arrays are of different lengths')\n",
    "    def mse(any_slope, any_intercept):\n",
    "        estimate = any_slope*x + any_intercept\n",
    "        return (np.mean((y - estimate) ** 2)) \n",
    "    return minimize(mse) \n",
    "\n",
    "line_mse(height_df['midparentHeight'],height_df['childHeight'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def height_mse(any_slope, any_intercept):\n",
    "    estimate = any_slope*height_df['midparentHeight'] + any_intercept\n",
    "    return (np.mean((height_df['childHeight'] - estimate) ** 2)) \n",
    "\n",
    "[height_mse(0.96399942,0.03064211),height_mse(0.63736,22.63624)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An optimization warning\n",
    "\n",
    "It is clear that the minimize function does a poor job for this dataset (it was performing as expected for the housing data).\n",
    "\n",
    "The documentation for the datascience *minimize* function is here:\n",
    "\n",
    "http://data8.org/datascience/_modules/datascience/util.html#minimize\n",
    "\n",
    "and it shows that the function is based on the following scipy.optimize function using \"Powell\" method as default:\n",
    "\n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html\n",
    "\n",
    "Changing the method to \"CG\" leads to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line_mse_CG(x,y):\n",
    "    if len(x)!=len(y):\n",
    "        raise ValueError('arrays are of different lengths')\n",
    "    def mse(any_slope, any_intercept):\n",
    "        estimate = any_slope*x + any_intercept\n",
    "        return (np.mean((y - estimate) ** 2)) \n",
    "    return minimize(mse,method=\"CG\") \n",
    "\n",
    "line_mse_CG(height_df['midparentHeight'],height_df['childHeight'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "height_df['Fitted']=fitted_values(height_df['midparentHeight'],height_df['childHeight'])\n",
    "height_df.head(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "height_df['Residuals']=height_df['childHeight']-height_df['Fitted']\n",
    "height_df.head(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# plot of residuals - it is a model diagnostic plot\n",
    "plots.scatter(height_df['midparentHeight'],height_df['Residuals'])\n",
    "plots.plot([64,76],[0,0],color=\"black\",lw=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Asymmetry of the regression (least squares) line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulated data with correlation r\n",
    "# side note - generates a standard bivariate normal with correlation 0.5\n",
    "r = 0.5\n",
    "x_demo = np.random.normal(0, 1, 10000)\n",
    "z_demo = np.random.normal(0, 1, 10000)\n",
    "y_demo = r*x_demo + np.sqrt(1 - r**2)*z_demo\n",
    "\n",
    "[slope(x_demo,y_demo),slope(y_demo,x_demo)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scatter plot showing the regression line (black) and the symmetry line (red)\n",
    "plots.figure(figsize=(7,7))\n",
    "plots.xlim(-4, 4)\n",
    "plots.ylim(-4, 4)\n",
    "plots.scatter(x_demo, y_demo, s=10)\n",
    "plots.plot([-4,4],[-4,4], color='red', lw=2)\n",
    "xlims = np.array([-4,4])\n",
    "plots.plot(xlims, 0.5 * xlims + 0.0, lw=2,color='black');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Model ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_and_compare(true_slope, true_int, sample_size):\n",
    "    x = np.random.normal(30, 5, sample_size)\n",
    "    xlims = np.array([np.min(x), np.max(x)])\n",
    "    errors = np.random.normal(0, 6, sample_size)\n",
    "    y = (true_slope * x + true_int) + errors\n",
    "    sample = Table().with_columns('x', x, 'y', y)\n",
    "\n",
    "    plots.figure(figsize=(8, 18))\n",
    "\n",
    "    plots.subplot(4, 1, 1)\n",
    "    plots.scatter(x, y)\n",
    "    plots.plot(xlims, true_slope*xlims + true_int, lw=2, color='green')\n",
    "    plots.title('True Line, and Points Created')\n",
    "\n",
    "    plots.subplot(4, 1, 2)\n",
    "    plots.scatter(x, y)\n",
    "    plots.title('What We Get to See')\n",
    "\n",
    "    plots.subplot(4, 1, 3)\n",
    "    plots.scatter(x, y)\n",
    "    plots.plot(xlims, slope(x,y)*xlims + intercept(x,y), lw=2, color='black')\n",
    "    plots.title('Regression Line: Estimate of True Line')\n",
    "\n",
    "    plots.subplot(4, 1, 4)\n",
    "    plots.scatter(x, y)\n",
    "    plots.plot(xlims, slope(x,y)*xlims + intercept(x,y), lw=2, color='black')\n",
    "    plots.plot(xlims, true_slope*xlims + true_int, lw=2, color='green')\n",
    "    plots.title(\"Regression Line and True Line\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_and_compare(2,-5,10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference for the Slope ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# the default for sample in pandas data frame is \"replace=False\"\n",
    "# To sample with replacement\n",
    "height_df.sample(5,replace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the number of rows in a data frame\n",
    "height_df.shape[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plots.figure(figsize=(8, 18))\n",
    "plots.subplot(5, 1, 1)\n",
    "plots.scatter(height_df['midparentHeight'], height_df['childHeight'], s=10)\n",
    "plots.xlim([64,76])\n",
    "plots.title('Original sample')\n",
    "\n",
    "for i in np.arange(1, 5, 1):\n",
    "    plots.subplot(5,1,i+1)\n",
    "    rep_df = height_df.sample(height_df.shape[0],replace=True)\n",
    "    plots.scatter(rep_df['midparentHeight'], rep_df['childHeight'], s=10)\n",
    "    plots.xlim([64,76])\n",
    "    plots.title('Bootstrap sample '+str(i))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# plot of 10 bootstrapped lines\n",
    "plots.xlim([64,76])\n",
    "plots.ylim([62,72])\n",
    "xlims = np.array([64,76])\n",
    "\n",
    "for i in np.arange(1, 10, 1):\n",
    "    rep_df = height_df.sample(height_df.shape[0],replace=True)\n",
    "    plots.plot(xlims,slope(rep_df['midparentHeight'], rep_df['childHeight'])*xlims+\n",
    "               intercept(rep_df['midparentHeight'], rep_df['childHeight']))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# this is a function that outputs a 95% CI for the slope\n",
    "# input is a data frame, the name of the two columns, and nr of bootstraps\n",
    "\n",
    "def bootstrap_slope(df,x, y, repetitions):\n",
    "    # the number of observations\n",
    "    n=df.shape[0]\n",
    "    # Bootstrap the scatter, find the slope, collect\n",
    "    slopes = make_array()\n",
    "    for i in np.arange(repetitions):\n",
    "        bootstrap_sample = df.sample(n,replace=True)\n",
    "        bootstrap_slope = slope(bootstrap_sample[x], bootstrap_sample[y])\n",
    "        slopes = np.append(slopes, bootstrap_slope)\n",
    "    \n",
    "    # Find the endpoints of the 95% confidence interval for the true slope\n",
    "    left = percentile(2.5, slopes)\n",
    "    right = percentile(97.5, slopes)\n",
    "    \n",
    "    # Slope of the regression line from the original sample\n",
    "    observed_slope = slope(df[x],df[y])\n",
    "    \n",
    "    # Display results\n",
    "    print('Slope of regression line:', observed_slope)\n",
    "    print('Approximate 95%-confidence interval for the slope of the true line:')\n",
    "    print(left, 'to', right)\n",
    " \n",
    "\n",
    "bootstrap_slope(height_df,'midparentHeight','childHeight',1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the confidence interval for the slope can be used for a decision in testing if the true slope is zero.\n",
    "\n",
    "## Prediction of mean response: ##\n",
    "\n",
    "We want a confidence interval for: $\\mu_{y|x_0}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for the height data - predict mean height for children with midparentHeight=68\n",
    "0.63736*68+ 22.63624"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prediction_at(df, x, y, value_of_x):\n",
    "    b1 = slope(df[x],df[y])\n",
    "    b0 = intercept(df[x],df[y])\n",
    "    return b1 * value_of_x + b0\n",
    "\n",
    "prediction_at(height_df,'midparentHeight','childHeight',68)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_mean_response(df, x, y, new_x, repetitions):\n",
    "    # the number of observations\n",
    "    n=df.shape[0]\n",
    "    # Bootstrap the scatter, predict, collect\n",
    "    predictions = make_array()\n",
    "    for i in np.arange(repetitions):\n",
    "        resample = df.sample(n,replace=True)\n",
    "        predicted_y = prediction_at(resample, x, y, new_x)\n",
    "        predictions= np.append(predictions, predicted_y)\n",
    "\n",
    "    # Find the ends of the approximate 95% prediction interval\n",
    "    left = percentile(2.5, predictions)\n",
    "    right = percentile(97.5, predictions)\n",
    "\n",
    "    # Display results\n",
    "    plots.hist(predictions,bins=20)\n",
    "    plots.xlabel('predictions at x='+str(new_x))\n",
    "    plots.plot([left, right], [0, 0], color='yellow', lw=8);\n",
    "    print('Approximate 95%-confidence interval for height of true line:')\n",
    "    print(left, right, '(width =', right - left, ')') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bootstrap_mean_response(height_df,'midparentHeight','childHeight',68,1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the version of the function that outputs only the CI\n",
    "def bootstrap_mean_responseCI(df, x, y, new_x, repetitions):\n",
    "    # the number of observations\n",
    "    n=df.shape[0]\n",
    "    # Bootstrap the scatter, predict, collect\n",
    "    predictions = make_array()\n",
    "    for i in np.arange(repetitions):\n",
    "        resample = df.sample(n,replace=True)\n",
    "        predicted_y = prediction_at(resample, x, y, new_x)\n",
    "        predictions= np.append(predictions, predicted_y)\n",
    "\n",
    "    # Find the ends of the approximate 95% prediction interval\n",
    "    left = percentile(2.5, predictions)\n",
    "    right = percentile(97.5, predictions)\n",
    "\n",
    "    return [left,right]\n",
    "\n",
    "CI68=bootstrap_mean_responseCI(height_df,'midparentHeight','childHeight',68,1000)\n",
    "CI68"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CIs for various values of x\n",
    "CI64=bootstrap_mean_responseCI(height_df,'midparentHeight','childHeight',64,1000)\n",
    "CI70=bootstrap_mean_responseCI(height_df,'midparentHeight','childHeight',70,1000)\n",
    "CI76=bootstrap_mean_responseCI(height_df,'midparentHeight','childHeight',76,1000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plots.scatter(height_df['midparentHeight'],height_df['childHeight'])\n",
    "xlims = np.array([64,76])\n",
    "plots.plot([64, 64], [CI64[0],CI64[1]], color='red', lw=2)\n",
    "plots.plot([68, 68], [CI68[0],CI68[1]], color='red', lw=2)\n",
    "plots.plot([70, 70], [CI70[0],CI70[1]], color='red', lw=2)\n",
    "plots.plot([76, 76], [CI76[0],CI76[1]], color='red', lw=2)\n",
    "plots.plot(xlims, 0.63736 * xlims + 22.636, lw=2,color='black');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction of a future observation##\n",
    "\n",
    "For a future pair $(x_0,y_0)$ for which we only observe $x_0$: we want a CI for $y_0$.\n",
    "\n",
    "The justification for the CI is in the lecture slides."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bootstrap_predictionCI(df, x, y, new_x, repetitions):\n",
    "    # Slope and intercept of the regression line from the original sample\n",
    "    observed_slope = slope(df[x],df[y])\n",
    "    observed_intercept = intercept(df[x],df[y])\n",
    "\n",
    "    # the number of observations\n",
    "    n=df.shape[0]\n",
    "\n",
    "    # Bootstrap the scatter, predict, collect\n",
    "    prediction_error = make_array()\n",
    "    for i in np.arange(repetitions):\n",
    "        resample = df.sample(n,replace=True)\n",
    "        boot_slope=slope(resample[x],resample[y])\n",
    "        boot_intercept=intercept(resample[x],resample[y])\n",
    "        predicted_y = prediction_at(resample, x, y, new_x)\n",
    "        res=np.random.choice(residuals(resample[x],resample[y]),1)\n",
    "        er=(observed_intercept-boot_intercept)+(observed_slope-boot_slope)*new_x+res\n",
    "        prediction_error= np.append(prediction_error, er)\n",
    "\n",
    "    # Find the ends of the approximate 95% prediction interval\n",
    "    left = percentile(2.5, prediction_error)\n",
    "    right = percentile(97.5, prediction_error)\n",
    " \n",
    "    predicted_y = prediction_at(df, x, y, new_x)\n",
    "\n",
    "    return [predicted_y+left,predicted_y+right]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pCI68=bootstrap_predictionCI(height_df,'midparentHeight','childHeight',68,1000)\n",
    "pCI68"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prediction CIs for various values of x\n",
    "pCI64=bootstrap_predictionCI(height_df,'midparentHeight','childHeight',64,1000)\n",
    "pCI70=bootstrap_predictionCI(height_df,'midparentHeight','childHeight',70,1000)\n",
    "pCI76=bootstrap_predictionCI(height_df,'midparentHeight','childHeight',76,1000)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plots.scatter(height_df['midparentHeight'],height_df['childHeight'])\n",
    "xlims = np.array([64,76])\n",
    "plots.plot([64, 64], [pCI64[0],pCI64[1]], color='red', lw=2)\n",
    "plots.plot([68, 68], [pCI68[0],pCI68[1]], color='red', lw=2)\n",
    "plots.plot([70, 70], [pCI70[0],pCI70[1]], color='red', lw=2)\n",
    "plots.plot([76, 76], [pCI76[0],pCI76[1]], color='red', lw=2)\n",
    "plots.plot(xlims, 0.63736 * xlims + 22.636, lw=2,color='black');"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}