LAB resampling methods ISLR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will work through the LAB for resampling methods chapter in ISLR.\n",
    "\n",
    "Basic functions that perform least squares linear regression and other simple analyses come standard with the base distribution named MASS. So we use the library() function to run it. The cv.glm() function is part of the boot library\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "loadlibraries=function(){\n",
    "    library('MASS')\n",
    "    library('class')\n",
    "    library('ISLR')\n",
    "    library('boot')\n",
    "    print('Libraries Loaded')\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "\"package 'ISLR' was built under R version 3.5.2\""
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Libraries Loaded\"\n"
     ]
    }
   ],
   "source": [
    "loadlibraries()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The following objects are masked from Default (pos = 3):\n",
      "\n",
      "    balance, default, income, student\n",
      "\n",
      "The following objects are masked from Default (pos = 4):\n",
      "\n",
      "    balance, default, income, student\n",
      "\n",
      "The following objects are masked from Default (pos = 5):\n",
      "\n",
      "    balance, default, income, student\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>default</th><th scope=col>student</th><th scope=col>balance</th><th scope=col>income</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>No       </td><td>No       </td><td> 729.5265</td><td>44361.625</td></tr>\n",
       "\t<tr><td>No       </td><td>Yes      </td><td> 817.1804</td><td>12106.135</td></tr>\n",
       "\t<tr><td>No       </td><td>No       </td><td>1073.5492</td><td>31767.139</td></tr>\n",
       "\t<tr><td>No       </td><td>No       </td><td> 529.2506</td><td>35704.494</td></tr>\n",
       "\t<tr><td>No       </td><td>No       </td><td> 785.6559</td><td>38463.496</td></tr>\n",
       "\t<tr><td>No       </td><td>Yes      </td><td> 919.5885</td><td> 7491.559</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       " default & student & balance & income\\\\\n",
       "\\hline\n",
       "\t No        & No        &  729.5265 & 44361.625\\\\\n",
       "\t No        & Yes       &  817.1804 & 12106.135\\\\\n",
       "\t No        & No        & 1073.5492 & 31767.139\\\\\n",
       "\t No        & No        &  529.2506 & 35704.494\\\\\n",
       "\t No        & No        &  785.6559 & 38463.496\\\\\n",
       "\t No        & Yes       &  919.5885 &  7491.559\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "default | student | balance | income | \n",
       "|---|---|---|---|---|---|\n",
       "| No        | No        |  729.5265 | 44361.625 | \n",
       "| No        | Yes       |  817.1804 | 12106.135 | \n",
       "| No        | No        | 1073.5492 | 31767.139 | \n",
       "| No        | No        |  529.2506 | 35704.494 | \n",
       "| No        | No        |  785.6559 | 38463.496 | \n",
       "| No        | Yes       |  919.5885 |  7491.559 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  default student balance   income   \n",
       "1 No      No       729.5265 44361.625\n",
       "2 No      Yes      817.1804 12106.135\n",
       "3 No      No      1073.5492 31767.139\n",
       "4 No      No       529.2506 35704.494\n",
       "5 No      No       785.6559 38463.496\n",
       "6 No      Yes      919.5885  7491.559"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "attach(Default)\n",
    "head(Default)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = default ~ balance + income, family = binomial, \n",
       "    data = Default)\n",
       "\n",
       "Deviance Residuals: \n",
       "    Min       1Q   Median       3Q      Max  \n",
       "-2.4725  -0.1444  -0.0574  -0.0211   3.7245  \n",
       "\n",
       "Coefficients:\n",
       "              Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***\n",
       "balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***\n",
       "income       2.081e-05  4.985e-06   4.174 2.99e-05 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 2920.6  on 9999  degrees of freedom\n",
       "Residual deviance: 1579.0  on 9997  degrees of freedom\n",
       "AIC: 1585\n",
       "\n",
       "Number of Fisher Scoring iterations: 8\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lorf=glm(default~balance+income,data=Default,family = binomial)\n",
    "summary(lorf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>split1</th><th scope=col>split2</th><th scope=col>split3</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>0.0286</td><td>0.0276</td><td>0.0248</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " split1 & split2 & split3\\\\\n",
       "\\hline\n",
       "\t 0.0286 & 0.0276 & 0.0248\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "split1 | split2 | split3 | \n",
       "|---|\n",
       "| 0.0286 | 0.0276 | 0.0248 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  split1 split2 split3\n",
       "1 0.0286 0.0276 0.0248"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "validation_set=function(seed){\n",
    "    set.seed(seed)\n",
    "    train=sample(nrow(Default),nrow(Default)/2)\n",
    "    #We know fit the model on the training set and evaluate it using the validation set.\n",
    "    lorfv=glm(default~balance+income,data=Default[train,],family=binomial)\n",
    "    predictions=predict(lorfv,Default[-train,],type = 'response')\n",
    "    #contrasts(default) We discover the response that is superior to our choosen threshold is corrsponding to defaulting!\n",
    "    response=rep('No',nrow(Default)/2)\n",
    "    response[predictions>0.5]='Yes'\n",
    "    return(mean(response!=Default[-train,1]))\n",
    "}\n",
    "results=data.frame('split1'=validation_set(1),'split2'=validation_set(2),'split3'=validation_set(3))\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we would expect changing the split modifies the error rate slightly around 0.026, this is due to the fact we only using half of the available data to fit the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>split1</th><th scope=col>split2</th><th scope=col>split3</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>0.0288</td><td>0.0286</td><td>0.0248</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " split1 & split2 & split3\\\\\n",
       "\\hline\n",
       "\t 0.0288 & 0.0286 & 0.0248\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "split1 | split2 | split3 | \n",
       "|---|\n",
       "| 0.0288 | 0.0286 | 0.0248 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  split1 split2 split3\n",
       "1 0.0288 0.0286 0.0248"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "validation_set2=function(seed){\n",
    "    set.seed(seed) #We use same split so that the comparaison make sense!\n",
    "    train=sample(nrow(Default),nrow(Default)/2)\n",
    "    lorfv=glm(default~.,data=Default,subset=train,family=binomial)\n",
    "    predictions=predict(lorfv,Default[-train,],type = 'response')\n",
    "    response=rep('No',nrow(Default)/2)\n",
    "    response[predictions>0.5]='Yes'\n",
    "    return(mean(response!=Default[-train,1]))\n",
    "}\n",
    "results=data.frame('split1'=validation_set2(1),'split2'=validation_set2(2),'split3'=validation_set2(3))\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding the student dummy variable doesn't reduce the value of test error. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = default ~ balance + income, family = binomial, \n",
       "    data = Default)\n",
       "\n",
       "Deviance Residuals: \n",
       "    Min       1Q   Median       3Q      Max  \n",
       "-2.4725  -0.1444  -0.0574  -0.0211   3.7245  \n",
       "\n",
       "Coefficients:\n",
       "              Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***\n",
       "balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***\n",
       "income       2.081e-05  4.985e-06   4.174 2.99e-05 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 2920.6  on 9999  degrees of freedom\n",
       "Residual deviance: 1579.0  on 9997  degrees of freedom\n",
       "AIC: 1585\n",
       "\n",
       "Number of Fisher Scoring iterations: 8\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>balance</th><th scope=col>income</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>0.0002273731</td><td>4.985167e-06</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       " balance & income\\\\\n",
       "\\hline\n",
       "\t 0.0002273731 & 4.985167e-06\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "balance | income | \n",
       "|---|\n",
       "| 0.0002273731 | 4.985167e-06 | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  balance      income      \n",
       "1 0.0002273731 4.985167e-06"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "set.seed(1)\n",
    "lorf=glm(default~balance+income,data=Default,family = binomial)\n",
    "summary(lorf)\n",
    "data.frame('balance'=summary(lorf)$coefficients[2,2],'income'=summary(lorf)$coefficients[3,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see above for balance glm returns a std error of: 0.0002273731 while for income it is: 4.985167e-06."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "#b-\n",
    "set.seed(1)\n",
    "boot.fn=function(data,index){\n",
    "    return(coef(glm(default~balance+income,data=Default,subset=index,family = binomial)))  \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "ORDINARY NONPARAMETRIC BOOTSTRAP\n",
       "\n",
       "\n",
       "Call:\n",
       "boot(data = Default, statistic = boot.fn, R = 1000)\n",
       "\n",
       "\n",
       "Bootstrap Statistics :\n",
       "         original        bias     std. error\n",
       "t1* -1.154047e+01 -8.008379e-03 4.239273e-01\n",
       "t2*  5.647103e-03  2.299970e-06 2.267955e-04\n",
       "t3*  2.080898e-05  5.870933e-08 4.582525e-06"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot(Default,boot.fn,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are slightly different!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>Year</th><th scope=col>Lag1</th><th scope=col>Lag2</th><th scope=col>Lag3</th><th scope=col>Lag4</th><th scope=col>Lag5</th><th scope=col>Volume</th><th scope=col>Today</th><th scope=col>Direction</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>1990     </td><td> 0.816   </td><td> 1.572   </td><td>-3.936   </td><td>-0.229   </td><td>-3.484   </td><td>0.1549760</td><td>-0.270   </td><td>Down     </td></tr>\n",
       "\t<tr><td>1990     </td><td>-0.270   </td><td> 0.816   </td><td> 1.572   </td><td>-3.936   </td><td>-0.229   </td><td>0.1485740</td><td>-2.576   </td><td>Down     </td></tr>\n",
       "\t<tr><td>1990     </td><td>-2.576   </td><td>-0.270   </td><td> 0.816   </td><td> 1.572   </td><td>-3.936   </td><td>0.1598375</td><td> 3.514   </td><td>Up       </td></tr>\n",
       "\t<tr><td>1990     </td><td> 3.514   </td><td>-2.576   </td><td>-0.270   </td><td> 0.816   </td><td> 1.572   </td><td>0.1616300</td><td> 0.712   </td><td>Up       </td></tr>\n",
       "\t<tr><td>1990     </td><td> 0.712   </td><td> 3.514   </td><td>-2.576   </td><td>-0.270   </td><td> 0.816   </td><td>0.1537280</td><td> 1.178   </td><td>Up       </td></tr>\n",
       "\t<tr><td>1990     </td><td> 1.178   </td><td> 0.712   </td><td> 3.514   </td><td>-2.576   </td><td>-0.270   </td><td>0.1544440</td><td>-1.372   </td><td>Down     </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lllllllll}\n",
       " Year & Lag1 & Lag2 & Lag3 & Lag4 & Lag5 & Volume & Today & Direction\\\\\n",
       "\\hline\n",
       "\t 1990      &  0.816    &  1.572    & -3.936    & -0.229    & -3.484    & 0.1549760 & -0.270    & Down     \\\\\n",
       "\t 1990      & -0.270    &  0.816    &  1.572    & -3.936    & -0.229    & 0.1485740 & -2.576    & Down     \\\\\n",
       "\t 1990      & -2.576    & -0.270    &  0.816    &  1.572    & -3.936    & 0.1598375 &  3.514    & Up       \\\\\n",
       "\t 1990      &  3.514    & -2.576    & -0.270    &  0.816    &  1.572    & 0.1616300 &  0.712    & Up       \\\\\n",
       "\t 1990      &  0.712    &  3.514    & -2.576    & -0.270    &  0.816    & 0.1537280 &  1.178    & Up       \\\\\n",
       "\t 1990      &  1.178    &  0.712    &  3.514    & -2.576    & -0.270    & 0.1544440 & -1.372    & Down     \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "Year | Lag1 | Lag2 | Lag3 | Lag4 | Lag5 | Volume | Today | Direction | \n",
       "|---|---|---|---|---|---|\n",
       "| 1990      |  0.816    |  1.572    | -3.936    | -0.229    | -3.484    | 0.1549760 | -0.270    | Down      | \n",
       "| 1990      | -0.270    |  0.816    |  1.572    | -3.936    | -0.229    | 0.1485740 | -2.576    | Down      | \n",
       "| 1990      | -2.576    | -0.270    |  0.816    |  1.572    | -3.936    | 0.1598375 |  3.514    | Up        | \n",
       "| 1990      |  3.514    | -2.576    | -0.270    |  0.816    |  1.572    | 0.1616300 |  0.712    | Up        | \n",
       "| 1990      |  0.712    |  3.514    | -2.576    | -0.270    |  0.816    | 0.1537280 |  1.178    | Up        | \n",
       "| 1990      |  1.178    |  0.712    |  3.514    | -2.576    | -0.270    | 0.1544440 | -1.372    | Down      | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  Year Lag1   Lag2   Lag3   Lag4   Lag5   Volume    Today  Direction\n",
       "1 1990  0.816  1.572 -3.936 -0.229 -3.484 0.1549760 -0.270 Down     \n",
       "2 1990 -0.270  0.816  1.572 -3.936 -0.229 0.1485740 -2.576 Down     \n",
       "3 1990 -2.576 -0.270  0.816  1.572 -3.936 0.1598375  3.514 Up       \n",
       "4 1990  3.514 -2.576 -0.270  0.816  1.572 0.1616300  0.712 Up       \n",
       "5 1990  0.712  3.514 -2.576 -0.270  0.816 0.1537280  1.178 Up       \n",
       "6 1990  1.178  0.712  3.514 -2.576 -0.270 0.1544440 -1.372 Down     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head(Weekly)\n",
    "attach(Weekly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly)\n",
       "\n",
       "Deviance Residuals: \n",
       "   Min      1Q  Median      3Q     Max  \n",
       "-1.623  -1.261   1.001   1.083   1.506  \n",
       "\n",
       "Coefficients:\n",
       "            Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)  0.22122    0.06147   3.599 0.000319 ***\n",
       "Lag1        -0.03872    0.02622  -1.477 0.139672    \n",
       "Lag2         0.06025    0.02655   2.270 0.023232 *  \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 1496.2  on 1088  degrees of freedom\n",
       "Residual deviance: 1488.2  on 1086  degrees of freedom\n",
       "AIC: 1494.2\n",
       "\n",
       "Number of Fisher Scoring iterations: 4\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lorf1=glm(Direction~Lag1+Lag2,data=Weekly,family=binomial)\n",
    "summary(lorf1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "glm(formula = Direction ~ Lag1 + Lag2, family = binomial, data = Weekly[-1, \n",
       "    ])\n",
       "\n",
       "Deviance Residuals: \n",
       "    Min       1Q   Median       3Q      Max  \n",
       "-1.6258  -1.2617   0.9999   1.0819   1.5071  \n",
       "\n",
       "Coefficients:\n",
       "            Estimate Std. Error z value Pr(>|z|)    \n",
       "(Intercept)  0.22324    0.06150   3.630 0.000283 ***\n",
       "Lag1        -0.03843    0.02622  -1.466 0.142683    \n",
       "Lag2         0.06085    0.02656   2.291 0.021971 *  \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "(Dispersion parameter for binomial family taken to be 1)\n",
       "\n",
       "    Null deviance: 1494.6  on 1087  degrees of freedom\n",
       "Residual deviance: 1486.5  on 1085  degrees of freedom\n",
       "AIC: 1492.5\n",
       "\n",
       "Number of Fisher Scoring iterations: 4\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lorf2=glm(Direction~Lag1+Lag2,data=Weekly[-1,],family=binomial)\n",
    "summary(lorf2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "Down\n",
       "<details>\n",
       "\t<summary style=display:list-item;cursor:pointer>\n",
       "\t\t<strong>Levels</strong>:\n",
       "\t</summary>\n",
       "\t<ol class=list-inline>\n",
       "\t\t<li>'Down'</li>\n",
       "\t\t<li>'Up'</li>\n",
       "\t</ol>\n",
       "</details>"
      ],
      "text/latex": [
       "Down\n",
       "\\emph{Levels}: \\begin{enumerate*}\n",
       "\\item 'Down'\n",
       "\\item 'Up'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "Down\n",
       "**Levels**: 1. 'Down'\n",
       "2. 'Up'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] Down\n",
       "Levels: Down Up"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>1:</strong> TRUE"
      ],
      "text/latex": [
       "\\textbf{1:} TRUE"
      ],
      "text/markdown": [
       "**1:** TRUE"
      ],
      "text/plain": [
       "   1 \n",
       "TRUE "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Weekly[1,]$Direction\n",
    "predict(lorf2,Weekly[1,],type='response')>0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the classification is inccorrect. Prediction is up while the real Directionis down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Error rate is:\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "0.449954086317723"
      ],
      "text/latex": [
       "0.449954086317723"
      ],
      "text/markdown": [
       "0.449954086317723"
      ],
      "text/plain": [
       "[1] 0.4499541"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "error=rep(NA,nrow(Weekly))\n",
    "for (i in 1:nrow(Weekly)){\n",
    "    lorf=glm(Direction~Lag1+Lag2,data=Weekly[-i,],family=binomial)\n",
    "    error[i]=(predict(lorf,Weekly[i,],type='response')>0.5 && Weekly[i,]$Direction=='Up') || (predict(lorf,Weekly[i,],type='response')<0.5 && Weekly[i,]$Direction=='Down')\n",
    "}\n",
    "print('Error rate is:')\n",
    "1-mean(error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well the model does slightly better than random guessing, nothing surprising since this is a complicated problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(1)\n",
    "x=rnorm(100)\n",
    "y=x-2*x^2+rnorm(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$n=100$ and $p=1$\n",
    "\n",
    "$y=x-2x^2+\\epsilon$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1kle1oaGh8\nfHyMjIyampqnp6eysrK9vb3Hx8fQ0NDZ2dnh4eHp6enw8PD/////JqYrAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAUPklEQVR4nO3dbVubWBSGURJ11LZq//+vHY22tZW8wQPsA2t9mMtOVbZ4\nbiGEmu4nMFq39ACwBkKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIA\nIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQ\nICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJ\nAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAh\nQYCQIEBIECAkCBASBAgJAoQEATOE1EFjBqzyfDgLbAKShAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFt136ffLeNE9Jm7fcXJXLhu22dkLZq\nv78okQvfbfOEtFVCihLSZjm1SxLSdh0L5J//r6NLCIl/OAQNIST+5kHRIELib0IaREj8Q0dD\nCGntrq9CRwMIaeUcX+YhpHXziGcmQlo3Ic1ESCuno3kIae10NAshQYCQIEBIECAkCBASBAgJ\nAoQEAUKCACGtimdflyKkNXE/0GKEtCLuUF2OkFZESMsR0proaDFCWhUdLUVIECAkCBASBAgJ\nAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQmqSu7yrEVKL/LujcoTUIP8Sth4hNUhI9Qip\nRcM60t6EhNSkYR0paTpC2grng5MS0lYIaVJC2gwdTUlI26GjCQkJAoS0AbMdijZ8zBPS+s32\n4GjLj8KEtHqzXa7b9HVBIa2ekOYgpFXpXceTLu/Pn3rDHQmpaf+u2yMredKO/i5psg1VJ6SG\n/dvN/OdWmz6b+4uQ2vVlFU+1rI9/TiH9IqR2fV3FU3V0uqT8FhskpIZ9XcVTdXSqpPwWWySk\nlp1axbEV7vTtEkJaqXGL/8tDr9HzrJ2Q1mncYeTLRYzARCsnpHUaFZKTuesJaaXGHpAu/+jD\ne26+OyGt1biHSJeHdHhXRzAh0ePykvZ/TDxTbUKix+VlCOmdkOhxRRlO7Q6ERJ8rbmZwseGN\nkOjl9rrrCKkpBZawx0O9hNSSCktYSL2E1JAaa7jEEOUIqSEVQvq4RKekfwipJcuv4MizRkt/\nEVMQUlOWXoKRp1+X/3EwASFxhVRH6ytJSFzjvYGRJ3ZCGvohBTfBsMUcSGCNHQlpsxZczivs\nSEhbtc4TrOUIaaO+hKSqUYRUxswr+UtHShpj1pB+PNx1b+7uf0y1iXbNvpJn+W3HmzFjSC83\n3R+3k2yiYQuvZCGNNGNI993u29Phrefvu+5+ik00bN6V/HVLua1vs8cZQ9p1T7/ffup2U2yi\nZcdXcn5l9m3r/f+M39ZGj2wzhtR1x/7w8X8+GbiJph37J9v5lXn86Dd+W1s9R3REKqVnEU6w\nMo9+ysC2hDTlhxy8Pkb6/nx4y2OkI/pW4RQr89QBKXFIGvcZmjTn5e/bT+duNy+TbKJxvQs5\ntjI/fZrjD8cSD5LGfoYWzfs80v3heaTd3YPnkfr1LuRYRxd8om1WEODOhlqmW8gXnrYpaRgh\nbcVlIV14bie3fwlpMy7t6Nd7nXjvjV5QOEVI63HulZkvOrH73cipJ4g3eon7FCGtxr+Le9Bi\n/6ujU9fIhfQ3Ia3Fv6t74HL/9SsZzjxpq6N/CGktjoQ0YMHv96c+Vke9hDSvCVdg36ndsGPS\n6Y88+jebrktIs5r0Z3n/zUXDQ7p0O58+7tptrYeQ5jTrWdHgF9Ibdijb+BmfkOY052I7bGno\n764bMKaQZviQgptYxrwdjdiWjq4kpHnNeGI3xcI+9cBJSNN/SMFNrN6gdX3mQ05fgthyR0Ja\nrzH3NRz/66sv5m2EkFZq+OW64X+/ZUJap4EndhcckobPtGpCWqWBx46zH6SjY4S0SkNPwoQy\nlJBW5+P342liVkJam4+EdDQvIa2MK2vLENLKCGkZQlqbKe5o4Cwhrc4Ul705R0g4HQwQEkIK\nEBJO7QKExE8XG8YTEgQIqRmOGpUJqRUex5QmpEa4slabkBohpNqE1IoSHRUYoSghNaPAIvav\nBY8SEhcbdnpZ4lA6OSFxMb/K+DghcbmhByQhhT6k4CYYYuhvy8tPUo2QmtTU0mxq2KGE1KKG\nfsg3M+hIQqroz+o78XrIM84z3MWDNvL1HCWkgv6svv512E5IF0/ayhd0lJDq+bP6jq3DZpbd\npSG186PhGCEV8rGUzofUzonQRX3s90KaxjZD+r2Wzp3anfks2alG2khHQqrj02o6c7Hh3GdJ\nTjW9VXQkpDoiy2nyNTnF69K2n5GQKkmsp6lX5RSffQ0dCamSxHqavqM5Xyu9HUJam4lP7FZx\n+JiAkNZqmvWuoyOEtKjpVuVUK15H/YS0pLOrffCybeMcrP6EFxPSgs6u9uExNBFSAyNeTEgL\nOrfax9TQwCJtIvZLCWlJFxyQRpzcDfzA2QhpclsJ6dxqX9FC67OmL09Ipa1nofVa0ZcnJAgQ\nEgQICQKEBAFCggAhQYCQGGdFl7DHEBKjDH5SdWUBCokxBt/ms6a7Gt4IiTGGhrSq++zeCKmM\nNtfViANSm1/wEUKqotWFNfghUptf7jFCKmJ1P6LPWdkXK6Qixoe0spXZGCFVMb4jJS1ISGWM\n70hJyxHSSghpWUJaCx0tSkjVrPi3nfRqdOx/CamYrR1YTn69De0KIdVy7KFOQ0vqt9EvHtvS\nDxUh1XJkYbW0pH65aOZTITV1/URIxRztqJkl9eHCmc8ckJr5qoVUTf+JXUNL6sOlM598iNTO\nFy2kJrS0pH4JzNzQFy2kcnpXT0NL6rcWZx5MSEs5LLOVXFdASEs59NITTYsPhxDSUvZ/9P3F\nMkMxnJCWcTQkp3ZtEtJCjp3abewx+moIaSnHLjbQJCFBgJAgQEgQICQIEBIECGla5y/LuXC3\nCkKa1PlnVz3/ug5CmtL5+33cEbQSQpqSkDZDSJNyarcVQpqWiw0bISQIEBIECAkChAQBQoIA\nIUGAkCBASFfytA99hHQdNyLQS0hXGXlrnAhXS0hXGReSw9l6Cek6YztS0koJ6UqjTuyEtFpC\nmpGO1ktIc9LRagkJAoREw+oc4oVEuwo96BQSzap0GVRINEtIBTbBCtTpSEi0rExHQoIEIUGA\nkCBASBAgJAgQEgQICQKEBAFCojF1noT9TEi0pdBtQZ8JiaZUulH1MyHRFCEV2wSNqtmRkAYq\n+c3chpq7XkiDHP+xWPPbzNSENMTxE/WiJx5MTUhDHA2p6kNhpiakQU4dkIS0RUIa5vhDJB1t\nkpDCdLRNQoIAIUGAkCBASBAgJAgQEgQIiRaUf1ZBSDSg/vPcQqK+Bu68EhL1fYRUuSUh0YCP\njgqXJCRa8NFR3ZKERCOEVHITNKdyR0KiHYU7EhIkCAkChAQBQoIAIUGAkCBASBAgJAgQUkDl\nJwqZh5DGK33rCvMQ0mi1b6ZctUK7XUijCWkplfa7kMar9P3cklI/wYR0vS/fvDLfzW0RUoVN\nDFfou7dxlb4TQrpWqZ+DG1fo+yCkawmJHkK6mo74av6QHm+67u77pJvIcm2B82YMqTt84G13\ncD/JJqbgAMQF5g7pvrt/+fnz+b57nGITE/CQiEvMHdKue3l7+6W7mWITExASl5g7pK779If4\nJqagIy4wd0j//QppN8UmJqEjzps1pLuHx+/dt9c3X+5PX20oFRKcN2tI7w5v7l6m2AQsZM7n\nkZ6eHh/v7g6XHO6/dtR9NnQTsIyRId08PMdGObIJaMDIkF4PHlO0JCQaMzKkl2//Xd/S+TM3\nIdGYwGOkHw8317UkJFYnc7Hhafd6XDp5089fn0BIrE0kpO/vd6LeXvoJhMTajA/p5eH1cHTz\n/eW1prsLP4GQWJuxIf14u9hw//T+Fxd+MiGxOmOfR3o9GD3+enL19P1zQzcBDRj7PNK5f+s6\njJBozNjnkWKDHN0ENMAvP4EAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACH18ZJIXElIPbxIH9cS0ldeNparCekrIXE1IfXQEdcSUh8dcSUhQYCQIEBI\nECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBJbM8lLbQuJjdnv\npyhJSGzLfj9JSUJiW4QECU7tIOFER8MTExL8MuJgJST4MObhk5Dgg5AgwakdJLjYAIsSEgQI\nCQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQIKTPJvlltmyBkD6Z5jee\nsQVC+mOi38HJFgjpDyExmJA+0RFDCekzHTGQkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBI\nECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQI\nCQKEBAEbCcnrtTCtbYTkFcSY2CZC8pqWTE1IELCJkJzaMbVthORiAwknltFGQoLxTp3YCAku\nc/KhtpDgMkKCBKd2kOBiA0xLSBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQE\nAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQ\nIGDWkH483HVv7u5/TLUJWMSMIb3cdH/cTrIJWMiMId13u29Ph7eev++6+yk2AQuZMaRd9/T7\n7aduN8UmYCEzhtR1x/4Q2wQsxBEJAuZ9jPT9+fCWx0iszZyXv28/XbW7eZlkE7CMeZ9Huj88\nj7S7e/A8EuvizgYIEBIECAkChAQBdULqPhv4Ofb76EhwqVnvbLi4lYGb2O+VxDJmDOlx6pD2\neyWxkDlP7Z52p//xxNhNCInFzPoY6en0jUGjN6EjljLvxYbHT/etTrEJHbGQOlftZt4EJAkJ\nAoQEAUuEdP75ViHRGCFBgJAgQEgQICQIEBIEuPwNAUKCACFBgJAgQEgQICQIEBIECAkChAQB\nQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAg\nQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBAS\nBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFC\nggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCCg\nrZD2++m3DQM0FdJ+ryRqaimk/V5JFCUkCGgpJKd2lNVUSC42UFVbIUFRQoIAIUGAkCBASBAg\nJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgoGhI0JgBqzwfTkih\nyYzSo84kJUapMEO/QpMZpUedSUqMUmGGfoUmM0qPOpOUGKXCDP0KTWaUHnUmKTFKhRn6FZrM\nKD3qTFJilAoz9Cs0mVF61JmkxCgVZuhXaDKj9KgzSYlRKszQr9BkRulRZ5ISo1SYoV+hyYzS\no84kJUapMEO/QpMZpUedSUqMUmGGfoUmM0qPOpOUGKXCDP0KTWaUHnUmKTFKhRmgeUKCACFB\ngJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQIKh/R40+3uX5ae4sNjhR11\nv6uzR0rskDdFlkmV3fHV/eFlAXbL76I3T0NeoCDt9rBHbpYe46DEDnlTZZkU2R1fPXX/vbz9\n3Ptv6UHePO0KrJsf3e7pbZIfSw/ys8gOeVNmmdTYHT3u3icr8f167G4LzHHffX/977fuYelB\nquyQN2WWyeIDnLH8HnrV3VeY4657/vn2I/hu6UGq7JBPlp9m8QFOe+lulx7h1VOF79SvEQpM\nUmSH/FFgmVTaHT0eD6czBRRYN4VC+llmjHcFlkml3fHV867AecxBgXUjpGMqLJNCu+Orl93i\nR+xfCqwbIR1RYpnU2R0fPr+q9O2yz5l8HqXAutkJqd/Cy+Rdnd3x4c/qfb65fS4ySol1837V\n7rnCVbufJXbIu8WXybsqu+Or78tfifmkwLp5ODyi/t7dLz3IQYEdclBlmRTZHV89F9lBHwqs\nm0p3NpTYIW/KLJMau6PHf133+dRqaRXmuDnsjyoLp8AO+VlomSw+wDFdlT30rsIcL4e7v5ee\n4kOFHfKz0DJZfABYAyFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGAkCBASBAgJAgQEgQICQKE\nBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkChAQBQoIAIUGA\nkCBASBAgJAgQUpNuux+v//3R/bf0IHwQUpOeu93rf3e7l6UH4YOQ2vTYPfx86L4tPQa/CKlR\nt91jd7f0EPwmpEY9d133vPQQ/CakVt1390uPwB9CapQjUi1CatTd62Ok26WH4Dchtenb64nd\nQ/e49Bj8IqQmvewOzyM5uStDSE367+POBid3VQgJAoQEAUKCACFBgJAgQEgQICQIEBIECAkC\nhAQBQoIAIUGAkCBASBAgJAgQEgQICQKEBAFCggAhQYCQIEBIECAkCBASBAgJAoQEAUKCACFB\ngJAgQEgQICQI+B+kmP11kFGXTQAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(x,y,col='cornflowerblue',pch=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>7.28816160667281</li>\n",
       "\t<li>0.937423637615551</li>\n",
       "\t<li>0.95662183010894</li>\n",
       "\t<li>0.953904892744803</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 7.28816160667281\n",
       "\\item 0.937423637615551\n",
       "\\item 0.95662183010894\n",
       "\\item 0.953904892744803\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 7.28816160667281\n",
       "2. 0.937423637615551\n",
       "3. 0.95662183010894\n",
       "4. 0.953904892744803\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 7.2881616 0.9374236 0.9566218 0.9539049"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "set.seed(1)\n",
    "cvError=rep(0,4)\n",
    "D=data.frame('x'=x,'y'=y)\n",
    "for (i in 1:4){\n",
    "    fit=glm(y~poly(x,i))\n",
    "    cvError[i]=cv.glm(D,fit)$delta[1]\n",
    "}\n",
    "cvError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             Estimate Std. Error   t value     Pr(>|t|)\n",
      "(Intercept) -1.550023  0.2600138 -5.961308 3.953542e-08\n",
      "poly(x, i)   6.188826  2.6001382  2.380191 1.923846e-02\n",
      "              Estimate Std. Error    t value     Pr(>|t|)\n",
      "(Intercept)  -1.550023 0.09580323 -16.179231 2.656229e-29\n",
      "poly(x, i)1   6.188826 0.95803228   6.459934 4.184810e-09\n",
      "poly(x, i)2 -23.948305 0.95803228 -24.997388 4.584330e-44\n",
      "               Estimate Std. Error    t value     Pr(>|t|)\n",
      "(Intercept)  -1.5500226 0.09626318 -16.101926 4.995066e-29\n",
      "poly(x, i)1   6.1888256 0.96263178   6.429068 4.971565e-09\n",
      "poly(x, i)2 -23.9483049 0.96263178 -24.877950 1.216703e-43\n",
      "poly(x, i)3   0.2641057 0.96263178   0.274358 7.843990e-01\n",
      "               Estimate Std. Error     t value     Pr(>|t|)\n",
      "(Intercept)  -1.5500226 0.09590514 -16.1620379 5.169227e-29\n",
      "poly(x, i)1   6.1888256 0.95905143   6.4530695 4.590732e-09\n",
      "poly(x, i)2 -23.9483049 0.95905143 -24.9708243 1.593826e-43\n",
      "poly(x, i)3   0.2641057 0.95905143   0.2753822 7.836207e-01\n",
      "poly(x, i)4   1.2570950 0.95905143   1.3107691 1.930956e-01\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>7.28816160667281</li>\n",
       "\t<li>0.937423637615552</li>\n",
       "\t<li>0.956621830108939</li>\n",
       "\t<li>0.953904892744804</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 7.28816160667281\n",
       "\\item 0.937423637615552\n",
       "\\item 0.956621830108939\n",
       "\\item 0.953904892744804\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 7.28816160667281\n",
       "2. 0.937423637615552\n",
       "3. 0.956621830108939\n",
       "4. 0.953904892744804\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 7.2881616 0.9374236 0.9566218 0.9539049"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "set.seed(2)\n",
    "cvError=rep(0,4)\n",
    "D=data.frame('x'=x,'y'=y)\n",
    "for (i in 1:4){\n",
    "    fit=glm(y~poly(x,i))\n",
    "    cvError[i]=cv.glm(D,fit)$delta[1]\n",
    "    print(summary(fit)$coefficients)\n",
    "}\n",
    "cvError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same result since we are using a LOOCV, cuz from a set of n observations we can extract n sets with n-1 as size $C^{n-1}_n=n$ and that's exactly the number of models we looking to train so, the random mechanism that chooses the n-1 folds won't affect results of the LOOCV, this cannot be said about other k-folds CVs with $k<n-1$. To put in simple words for babies to understand it: the random mechanism of choice doesn't have any choice compared when using lower tier CV when it has more possibilities to choose from than the acctual number of fitting operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "e- The smallest corresponds to the quadratic fit, which we would expect since the real function is as such. \n",
    "\n",
    "f- The linear model is insignificant, while the quadratic term is very significan in all models it is involved the same is for the first power term. Although beyond quadratic terms all terms are insignificant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6- "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>crim</th><th scope=col>zn</th><th scope=col>indus</th><th scope=col>chas</th><th scope=col>nox</th><th scope=col>rm</th><th scope=col>age</th><th scope=col>dis</th><th scope=col>rad</th><th scope=col>tax</th><th scope=col>ptratio</th><th scope=col>black</th><th scope=col>lstat</th><th scope=col>medv</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>0.00632</td><td>18     </td><td>2.31   </td><td>0      </td><td>0.538  </td><td>6.575  </td><td>65.2   </td><td>4.0900 </td><td>1      </td><td>296    </td><td>15.3   </td><td>396.90 </td><td>4.98   </td><td>24.0   </td></tr>\n",
       "\t<tr><td>0.02731</td><td> 0     </td><td>7.07   </td><td>0      </td><td>0.469  </td><td>6.421  </td><td>78.9   </td><td>4.9671 </td><td>2      </td><td>242    </td><td>17.8   </td><td>396.90 </td><td>9.14   </td><td>21.6   </td></tr>\n",
       "\t<tr><td>0.02729</td><td> 0     </td><td>7.07   </td><td>0      </td><td>0.469  </td><td>7.185  </td><td>61.1   </td><td>4.9671 </td><td>2      </td><td>242    </td><td>17.8   </td><td>392.83 </td><td>4.03   </td><td>34.7   </td></tr>\n",
       "\t<tr><td>0.03237</td><td> 0     </td><td>2.18   </td><td>0      </td><td>0.458  </td><td>6.998  </td><td>45.8   </td><td>6.0622 </td><td>3      </td><td>222    </td><td>18.7   </td><td>394.63 </td><td>2.94   </td><td>33.4   </td></tr>\n",
       "\t<tr><td>0.06905</td><td> 0     </td><td>2.18   </td><td>0      </td><td>0.458  </td><td>7.147  </td><td>54.2   </td><td>6.0622 </td><td>3      </td><td>222    </td><td>18.7   </td><td>396.90 </td><td>5.33   </td><td>36.2   </td></tr>\n",
       "\t<tr><td>0.02985</td><td> 0     </td><td>2.18   </td><td>0      </td><td>0.458  </td><td>6.430  </td><td>58.7   </td><td>6.0622 </td><td>3      </td><td>222    </td><td>18.7   </td><td>394.12 </td><td>5.21   </td><td>28.7   </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llllllllllllll}\n",
       " crim & zn & indus & chas & nox & rm & age & dis & rad & tax & ptratio & black & lstat & medv\\\\\n",
       "\\hline\n",
       "\t 0.00632 & 18      & 2.31    & 0       & 0.538   & 6.575   & 65.2    & 4.0900  & 1       & 296     & 15.3    & 396.90  & 4.98    & 24.0   \\\\\n",
       "\t 0.02731 &  0      & 7.07    & 0       & 0.469   & 6.421   & 78.9    & 4.9671  & 2       & 242     & 17.8    & 396.90  & 9.14    & 21.6   \\\\\n",
       "\t 0.02729 &  0      & 7.07    & 0       & 0.469   & 7.185   & 61.1    & 4.9671  & 2       & 242     & 17.8    & 392.83  & 4.03    & 34.7   \\\\\n",
       "\t 0.03237 &  0      & 2.18    & 0       & 0.458   & 6.998   & 45.8    & 6.0622  & 3       & 222     & 18.7    & 394.63  & 2.94    & 33.4   \\\\\n",
       "\t 0.06905 &  0      & 2.18    & 0       & 0.458   & 7.147   & 54.2    & 6.0622  & 3       & 222     & 18.7    & 396.90  & 5.33    & 36.2   \\\\\n",
       "\t 0.02985 &  0      & 2.18    & 0       & 0.458   & 6.430   & 58.7    & 6.0622  & 3       & 222     & 18.7    & 394.12  & 5.21    & 28.7   \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "crim | zn | indus | chas | nox | rm | age | dis | rad | tax | ptratio | black | lstat | medv | \n",
       "|---|---|---|---|---|---|\n",
       "| 0.00632 | 18      | 2.31    | 0       | 0.538   | 6.575   | 65.2    | 4.0900  | 1       | 296     | 15.3    | 396.90  | 4.98    | 24.0    | \n",
       "| 0.02731 |  0      | 7.07    | 0       | 0.469   | 6.421   | 78.9    | 4.9671  | 2       | 242     | 17.8    | 396.90  | 9.14    | 21.6    | \n",
       "| 0.02729 |  0      | 7.07    | 0       | 0.469   | 7.185   | 61.1    | 4.9671  | 2       | 242     | 17.8    | 392.83  | 4.03    | 34.7    | \n",
       "| 0.03237 |  0      | 2.18    | 0       | 0.458   | 6.998   | 45.8    | 6.0622  | 3       | 222     | 18.7    | 394.63  | 2.94    | 33.4    | \n",
       "| 0.06905 |  0      | 2.18    | 0       | 0.458   | 7.147   | 54.2    | 6.0622  | 3       | 222     | 18.7    | 396.90  | 5.33    | 36.2    | \n",
       "| 0.02985 |  0      | 2.18    | 0       | 0.458   | 6.430   | 58.7    | 6.0622  | 3       | 222     | 18.7    | 394.12  | 5.21    | 28.7    | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  crim    zn indus chas nox   rm    age  dis    rad tax ptratio black  lstat\n",
       "1 0.00632 18 2.31  0    0.538 6.575 65.2 4.0900 1   296 15.3    396.90 4.98 \n",
       "2 0.02731  0 7.07  0    0.469 6.421 78.9 4.9671 2   242 17.8    396.90 9.14 \n",
       "3 0.02729  0 7.07  0    0.469 7.185 61.1 4.9671 2   242 17.8    392.83 4.03 \n",
       "4 0.03237  0 2.18  0    0.458 6.998 45.8 6.0622 3   222 18.7    394.63 2.94 \n",
       "5 0.06905  0 2.18  0    0.458 7.147 54.2 6.0622 3   222 18.7    396.90 5.33 \n",
       "6 0.02985  0 2.18  0    0.458 6.430 58.7 6.0622 3   222 18.7    394.12 5.21 \n",
       "  medv\n",
       "1 24.0\n",
       "2 21.6\n",
       "3 34.7\n",
       "4 33.4\n",
       "5 36.2\n",
       "6 28.7"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "22.5328063241107"
      ],
      "text/latex": [
       "22.5328063241107"
      ],
      "text/markdown": [
       "22.5328063241107"
      ],
      "text/plain": [
       "[1] 22.53281"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "attach(Boston)\n",
    "head(Boston)\n",
    "mu=mean(medv)\n",
    "mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Var(\\bar{X})=Var(\\sum_i^n{\\frac{1}{n} x_i})$\n",
    "$=\\frac{1}{n^2}\\sum_i^n{Var(x_i)}$\n",
    "$=\\frac{1}{n^2}\\sum_i^n{\\sigma^2}$\n",
    "$=\\frac{\\sigma^2}{n^2}\\sum_i^n{1}$\n",
    "$=\\frac{\\sigma^2}{n}$\n",
    "Now that we can't have $\\sigma$ of the population we will use $\\hat{\\sigma}^2$ from data set.\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.408861147497535"
      ],
      "text/latex": [
       "0.408861147497535"
      ],
      "text/markdown": [
       "0.408861147497535"
      ],
      "text/plain": [
       "[1] 0.4088611"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "stdErr=sd(medv)/sqrt(nrow(Boston))\n",
    "stdErr"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Now let's estimate the std error using the bootstrap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "ORDINARY NONPARAMETRIC BOOTSTRAP\n",
       "\n",
       "\n",
       "Call:\n",
       "boot(data = medv, statistic = boot.fn, R = 1000)\n",
       "\n",
       "\n",
       "Bootstrap Statistics :\n",
       "    original     bias    std. error\n",
       "t1* 22.53281 -0.0148002   0.4153876"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.fn=function(data,index){\n",
    "    return(mean(data[index]))\n",
    "}\n",
    "bootstrap=boot(medv,boot.fn,1000)\n",
    "bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two values we get are very close. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>21.7020311241107</li>\n",
       "\t<li>23.3635815241107</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 21.7020311241107\n",
       "\\item 23.3635815241107\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 21.7020311241107\n",
       "2. 23.3635815241107\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 21.70203 23.36358"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "\tOne Sample t-test\n",
       "\n",
       "data:  medv\n",
       "t = 55.111, df = 505, p-value < 2.2e-16\n",
       "alternative hypothesis: true mean is not equal to 0\n",
       "95 percent confidence interval:\n",
       " 21.72953 23.33608\n",
       "sample estimates:\n",
       "mean of x \n",
       " 22.53281 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#d-\n",
    "confinter=c(bootstrap$t0-2*0.4153876,bootstrap$t0+2*0.4153876)\n",
    "confinter\n",
    "t.test(medv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The confidence interval obtained from the t.test is narrower than the one given by the bootstrap, but not by a big margin tho.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "21.2"
      ],
      "text/latex": [
       "21.2"
      ],
      "text/markdown": [
       "21.2"
      ],
      "text/plain": [
       "[1] 21.2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "median=median(medv)\n",
    "median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "ORDINARY NONPARAMETRIC BOOTSTRAP\n",
       "\n",
       "\n",
       "Call:\n",
       "boot(data = medv, statistic = boot.fn, R = 1000)\n",
       "\n",
       "\n",
       "Bootstrap Statistics :\n",
       "    original   bias    std. error\n",
       "t1*     21.2 -0.01975   0.3725066"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.fn=function(data,index){\n",
    "    return(median(data[index]))\n",
    "}\n",
    "bootstrap=boot(medv,boot.fn,1000)\n",
    "bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value predicted by the bootstrap for the median is exactly the one we estimated by median(). And the standard error is somewhat low."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>10%:</strong> 12.75"
      ],
      "text/latex": [
       "\\textbf{10\\textbackslash{}\\%:} 12.75"
      ],
      "text/markdown": [
       "**10%:** 12.75"
      ],
      "text/plain": [
       "  10% \n",
       "12.75 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "quantile(medv,0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constructing percentiles consists on deviding the dataset into 100 part. The first percentile per say is the value for which size of data*1% of the observation is below its value, the 2nd percentile is the value for which size of data*2% of the observations is below it and so on. It is usefull to know that the n^th percentile is the n% quantile. So the tenth percentile is the 0,1 quantile thus the result above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "ORDINARY NONPARAMETRIC BOOTSTRAP\n",
       "\n",
       "\n",
       "Call:\n",
       "boot(data = medv, statistic = boot.fn, R = 1000)\n",
       "\n",
       "\n",
       "Bootstrap Statistics :\n",
       "    original  bias    std. error\n",
       "t1*    12.75 0.02945   0.4870446"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "boot.fn=function(data,index){\n",
    "    return(quantile(data[index],0.1))\n",
    "}\n",
    "bootstrap=boot(medv,boot.fn,1000)\n",
    "bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value is exactly the same and the standard error is somezhat low."
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}