Chapter 5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ed2facd6",
   "metadata": {},
   "source": [
    "# Chapter 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "226d1833",
   "metadata": {},
   "source": [
    "## Exercise 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71a89d5a",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7ca53f0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "from ISLP import load_data\n",
    "from ISLP.models import (ModelSpec as MS,\n",
    "                         summarize,\n",
    "                         poly)\n",
    "from sklearn.model_selection import train_test_split\n",
    "from functools import partial\n",
    "from sklearn.model_selection import \\\n",
    "     (cross_validate,\n",
    "      KFold,\n",
    "      ShuffleSplit)\n",
    "from sklearn.base import clone\n",
    "from ISLP.models import sklearn_sm\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import confusion_matrix, accuracy_score\n",
    "import statsmodels.formula.api as smf "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7ac69cdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "Default=load_data(\"Default\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6c5a0bc6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0       No\n",
       "1       No\n",
       "2       No\n",
       "3       No\n",
       "4       No\n",
       "        ..\n",
       "9995    No\n",
       "9996    No\n",
       "9997    No\n",
       "9998    No\n",
       "9999    No\n",
       "Name: default, Length: 10000, dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Default.default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ee257fde",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>        <td>default</td>     <th>  No. Observations:  </th>  <td> 10000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>  9997</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     2</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>         <td>Logit</td>      <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -789.48</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Mon, 11 Sep 2023</td> <th>  Deviance:          </th> <td>  1579.0</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>18:10:52</td>     <th>  Pearson chi2:      </th> <td>6.95e+03</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>9</td>        <th>  Pseudo R-squ. (CS):</th>  <td>0.1256</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>intercept</th> <td>  -11.5405</td> <td>    0.435</td> <td>  -26.544</td> <td> 0.000</td> <td>  -12.393</td> <td>  -10.688</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>balance</th>   <td>    0.0056</td> <td>    0.000</td> <td>   24.835</td> <td> 0.000</td> <td>    0.005</td> <td>    0.006</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>income</th>    <td> 2.081e-05</td> <td> 4.99e-06</td> <td>    4.174</td> <td> 0.000</td> <td>  1.1e-05</td> <td> 3.06e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                default   No. Observations:                10000\n",
       "Model:                            GLM   Df Residuals:                     9997\n",
       "Model Family:                Binomial   Df Model:                            2\n",
       "Link Function:                  Logit   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -789.48\n",
       "Date:                Mon, 11 Sep 2023   Deviance:                       1579.0\n",
       "Time:                        18:10:52   Pearson chi2:                 6.95e+03\n",
       "No. Iterations:                     9   Pseudo R-squ. (CS):             0.1256\n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "intercept    -11.5405      0.435    -26.544      0.000     -12.393     -10.688\n",
       "balance        0.0056      0.000     24.835      0.000       0.005       0.006\n",
       "income      2.081e-05   4.99e-06      4.174      0.000     1.1e-05    3.06e-05\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "allvars = Default.columns.drop(['default','student'])\n",
    "design=MS(allvars)\n",
    "X=design.fit_transform(Default)\n",
    "y=Default.default=='Yes'\n",
    "glm=sm.GLM(y,X,family=sm.families.Binomial())\n",
    "res=glm.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f854de41",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce81c01",
   "metadata": {},
   "source": [
    "#### (i) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0e3492eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>default</th>\n",
       "      <th>student</th>\n",
       "      <th>balance</th>\n",
       "      <th>income</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9394</th>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>13911.441282</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>898</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>396.985985</td>\n",
       "      <td>55454.631040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2398</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>1046.416673</td>\n",
       "      <td>47598.307369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5906</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>836.343137</td>\n",
       "      <td>34559.158405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2343</th>\n",
       "      <td>No</td>\n",
       "      <td>Yes</td>\n",
       "      <td>534.692907</td>\n",
       "      <td>18729.566240</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3996</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>794.176213</td>\n",
       "      <td>43335.985041</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5889</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>81.531865</td>\n",
       "      <td>40847.811311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4577</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>1137.791157</td>\n",
       "      <td>21103.429621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8600</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>646.081240</td>\n",
       "      <td>36686.959447</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>847</th>\n",
       "      <td>No</td>\n",
       "      <td>No</td>\n",
       "      <td>388.548350</td>\n",
       "      <td>18007.832960</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5000 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     default student      balance        income\n",
       "9394      No     Yes     0.000000  13911.441282\n",
       "898       No      No   396.985985  55454.631040\n",
       "2398      No      No  1046.416673  47598.307369\n",
       "5906      No      No   836.343137  34559.158405\n",
       "2343      No     Yes   534.692907  18729.566240\n",
       "...      ...     ...          ...           ...\n",
       "3996      No      No   794.176213  43335.985041\n",
       "5889      No      No    81.531865  40847.811311\n",
       "4577      No      No  1137.791157  21103.429621\n",
       "8600      No      No   646.081240  36686.959447\n",
       "847       No      No   388.548350  18007.832960\n",
       "\n",
       "[5000 rows x 4 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Default_train, Default_valid = train_test_split(Default, test_size=5000, random_state=0)\n",
    "Default_valid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "946f6a94",
   "metadata": {},
   "source": [
    "#### (ii) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6b660a23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>        <td>default</td>     <th>  No. Observations:  </th>  <td>  5000</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>  4997</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     2</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>         <td>Logit</td>      <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -364.78</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Mon, 11 Sep 2023</td> <th>  Deviance:          </th> <td>  729.56</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>18:10:52</td>     <th>  Pearson chi2:      </th> <td>2.57e+03</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>9</td>        <th>  Pseudo R-squ. (CS):</th>  <td>0.1125</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>intercept</th> <td>  -11.3896</td> <td>    0.635</td> <td>  -17.935</td> <td> 0.000</td> <td>  -12.634</td> <td>  -10.145</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>balance</th>   <td>    0.0056</td> <td>    0.000</td> <td>   16.792</td> <td> 0.000</td> <td>    0.005</td> <td>    0.006</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>income</th>    <td>  1.58e-05</td> <td> 7.34e-06</td> <td>    2.151</td> <td> 0.031</td> <td>  1.4e-06</td> <td> 3.02e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                default   No. Observations:                 5000\n",
       "Model:                            GLM   Df Residuals:                     4997\n",
       "Model Family:                Binomial   Df Model:                            2\n",
       "Link Function:                  Logit   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -364.78\n",
       "Date:                Mon, 11 Sep 2023   Deviance:                       729.56\n",
       "Time:                        18:10:52   Pearson chi2:                 2.57e+03\n",
       "No. Iterations:                     9   Pseudo R-squ. (CS):             0.1125\n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "intercept    -11.3896      0.635    -17.935      0.000     -12.634     -10.145\n",
       "balance        0.0056      0.000     16.792      0.000       0.005       0.006\n",
       "income       1.58e-05   7.34e-06      2.151      0.031     1.4e-06    3.02e-05\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "allvars_train = Default_train.columns.drop(['default','student'])\n",
    "design_train=MS(allvars_train)\n",
    "X_train=design_train.fit_transform(Default_train)\n",
    "y_train=Default_train.default=='Yes'\n",
    "glm=sm.GLM(y_train,X_train,family=sm.families.Binomial())\n",
    "res=glm.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e08be31",
   "metadata": {},
   "source": [
    "#### (iii)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "55371cfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_valid=design_train.transform(Default_valid)\n",
    "valid_pred = res.predict(X_valid)\n",
    "labels=np.array(['No']*5000)\n",
    "labels[valid_pred>0.5]='Yes'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "072b5955",
   "metadata": {},
   "source": [
    "#### (iv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a6fbb08c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.03979999999999995"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_valid = Default_valid['default']\n",
    "1-accuracy_score(y_valid, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d469894",
   "metadata": {},
   "source": [
    "## Exercise 6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfc7811a",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "982ea4c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#We should add a column of dummy variables related to 'default' because without this manipulation\n",
    "#GLM used hereafter won't succed in converting categorical variables into dummy ones.\n",
    "Default['default_yes'] = (Default['default'] == 'Yes').astype('int')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "33022df1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>-11.540500</td>\n",
       "      <td>0.435000</td>\n",
       "      <td>-26.544</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>balance</th>\n",
       "      <td>0.005600</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>24.835</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>income</th>\n",
       "      <td>0.000021</td>\n",
       "      <td>0.000005</td>\n",
       "      <td>4.174</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                coef   std err       z  P>|z|\n",
       "intercept -11.540500  0.435000 -26.544    0.0\n",
       "balance     0.005600  0.000000  24.835    0.0\n",
       "income      0.000021  0.000005   4.174    0.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "allvars=Default.columns.drop(['default','student','default_yes'])\n",
    "design=MS(allvars)\n",
    "X=design.fit_transform(Default)\n",
    "y=Default.default=='Yes'\n",
    "model=sm.GLM(y,X,family=sm.families.Binomial())\n",
    "result=model.fit()\n",
    "summarize(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "caf2d4c4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "intercept    0.434772\n",
       "balance      0.000227\n",
       "income       0.000005\n",
       "dtype: float64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.bse#Generating the standard errors ofthe coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e87bc73b",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "90011fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def boot_SE(func,D,n=None,B=1000,seed=0):\n",
    "    rng = np.random.default_rng(seed)\n",
    "    first_, second_ = 0, 0\n",
    "    n = n or D.shape[0]\n",
    "    for _ in range(B):\n",
    "        idx = rng.choice(D.index,n,replace=True)\n",
    "        value = func(D,idx)\n",
    "        first_ += value\n",
    "        second_ += value**2\n",
    "    return np.sqrt(second_ / B - (first_ / B)**2)\n",
    "def boot_fn(model_matrix,response,D,idx):\n",
    "    D_=D.loc[idx]\n",
    "    y=D_[response]\n",
    "    X=clone(model_matrix).fit_transform(D_)\n",
    "    model=sm.GLM(y,X,family=sm.families.Binomial())\n",
    "    result=model.fit()\n",
    "    coef_income=result.params[1]\n",
    "    coef_balance=result.params[2]\n",
    "    return result.params"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c4609af",
   "metadata": {},
   "source": [
    "### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f67408ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.37418699e+01,  4.46108654e-05,  6.54456450e-03],\n",
       "       [-9.73343593e+00,  1.29161125e-05,  4.36172555e-03],\n",
       "       [-1.09652934e+01,  6.88991070e-06,  5.60127348e-03],\n",
       "       [-1.29994735e+01,  5.80207964e-05,  5.99463447e-03],\n",
       "       [-1.76482871e+01,  3.85717484e-05,  8.94706073e-03],\n",
       "       [-1.37799054e+01,  7.20354708e-05,  5.64428958e-03],\n",
       "       [-1.00286504e+01,  1.25804073e-05,  5.16837631e-03],\n",
       "       [-1.36627122e+01,  3.52359713e-05,  6.88246994e-03],\n",
       "       [-1.30954565e+01,  8.56134990e-06,  7.38734450e-03],\n",
       "       [-1.41538872e+01,  3.45878793e-05,  7.22152035e-03]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hp_func=partial(boot_fn,MS(['income','balance']),'default_yes')\n",
    "#Below is a test for hp_func instantiation\n",
    "rng = np.random.default_rng(0)\n",
    "np.array([hp_func(Default,\n",
    "          rng.choice(392,\n",
    "                     392,\n",
    "                     replace=True)) for _ in range(10)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c338caa3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "intercept    0.425280\n",
       "income       0.000005\n",
       "balance      0.000227\n",
       "dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boot_SE(hp_func,Default,B=1000,seed=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f1fd9ae",
   "metadata": {},
   "source": [
    "### (d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab58f5ff",
   "metadata": {},
   "source": [
    "Clearly we obtain approximate values by both methods."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b274aba",
   "metadata": {},
   "source": [
    "## Exercise 7"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01a4a61e",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3b8b9772",
   "metadata": {},
   "outputs": [],
   "source": [
    "Weekly=load_data('Weekly')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "28f5bb42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Year</th>\n",
       "      <th>Lag1</th>\n",
       "      <th>Lag2</th>\n",
       "      <th>Lag3</th>\n",
       "      <th>Lag4</th>\n",
       "      <th>Lag5</th>\n",
       "      <th>Volume</th>\n",
       "      <th>Today</th>\n",
       "      <th>Direction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1990</td>\n",
       "      <td>0.816</td>\n",
       "      <td>1.572</td>\n",
       "      <td>-3.936</td>\n",
       "      <td>-0.229</td>\n",
       "      <td>-3.484</td>\n",
       "      <td>0.154976</td>\n",
       "      <td>-0.270</td>\n",
       "      <td>Down</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1990</td>\n",
       "      <td>-0.270</td>\n",
       "      <td>0.816</td>\n",
       "      <td>1.572</td>\n",
       "      <td>-3.936</td>\n",
       "      <td>-0.229</td>\n",
       "      <td>0.148574</td>\n",
       "      <td>-2.576</td>\n",
       "      <td>Down</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1990</td>\n",
       "      <td>-2.576</td>\n",
       "      <td>-0.270</td>\n",
       "      <td>0.816</td>\n",
       "      <td>1.572</td>\n",
       "      <td>-3.936</td>\n",
       "      <td>0.159837</td>\n",
       "      <td>3.514</td>\n",
       "      <td>Up</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1990</td>\n",
       "      <td>3.514</td>\n",
       "      <td>-2.576</td>\n",
       "      <td>-0.270</td>\n",
       "      <td>0.816</td>\n",
       "      <td>1.572</td>\n",
       "      <td>0.161630</td>\n",
       "      <td>0.712</td>\n",
       "      <td>Up</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1990</td>\n",
       "      <td>0.712</td>\n",
       "      <td>3.514</td>\n",
       "      <td>-2.576</td>\n",
       "      <td>-0.270</td>\n",
       "      <td>0.816</td>\n",
       "      <td>0.153728</td>\n",
       "      <td>1.178</td>\n",
       "      <td>Up</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Year   Lag1   Lag2   Lag3   Lag4   Lag5    Volume  Today Direction\n",
       "0  1990  0.816  1.572 -3.936 -0.229 -3.484  0.154976 -0.270      Down\n",
       "1  1990 -0.270  0.816  1.572 -3.936 -0.229  0.148574 -2.576      Down\n",
       "2  1990 -2.576 -0.270  0.816  1.572 -3.936  0.159837  3.514        Up\n",
       "3  1990  3.514 -2.576 -0.270  0.816  1.572  0.161630  0.712        Up\n",
       "4  1990  0.712  3.514 -2.576 -0.270  0.816  0.153728  1.178        Up"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Weekly.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "80626c70",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "RangeIndex: 1089 entries, 0 to 1088\n",
      "Data columns (total 9 columns):\n",
      " #   Column     Non-Null Count  Dtype  \n",
      "---  ------     --------------  -----  \n",
      " 0   Year       1089 non-null   int64  \n",
      " 1   Lag1       1089 non-null   float64\n",
      " 2   Lag2       1089 non-null   float64\n",
      " 3   Lag3       1089 non-null   float64\n",
      " 4   Lag4       1089 non-null   float64\n",
      " 5   Lag5       1089 non-null   float64\n",
      " 6   Volume     1089 non-null   float64\n",
      " 7   Today      1089 non-null   float64\n",
      " 8   Direction  1089 non-null   object \n",
      "dtypes: float64(7), int64(1), object(1)\n",
      "memory usage: 76.7+ KB\n"
     ]
    }
   ],
   "source": [
    "Weekly.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b66b9cfb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>0.2212</td>\n",
       "      <td>0.061</td>\n",
       "      <td>3.599</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lag1</th>\n",
       "      <td>-0.0387</td>\n",
       "      <td>0.026</td>\n",
       "      <td>-1.477</td>\n",
       "      <td>0.140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lag2</th>\n",
       "      <td>0.0602</td>\n",
       "      <td>0.027</td>\n",
       "      <td>2.270</td>\n",
       "      <td>0.023</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             coef  std err      z  P>|z|\n",
       "intercept  0.2212    0.061  3.599  0.000\n",
       "Lag1      -0.0387    0.026 -1.477  0.140\n",
       "Lag2       0.0602    0.027  2.270  0.023"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Fitting a logistic regression of Lag1 and Lag2 as exog and Direction as endog\n",
    "allvars=Weekly.columns.drop(['Direction','Year','Volume','Today','Lag3','Lag4','Lag5'])\n",
    "design=MS(allvars)\n",
    "X=design.fit_transform(Weekly)\n",
    "y=Weekly.Direction=='Up'\n",
    "model=sm.GLM(y,X,family=sm.families.Binomial())\n",
    "result=model.fit()\n",
    "summarize(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3f443af",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "65c4fa21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>z</th>\n",
       "      <th>P&gt;|z|</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>intercept</th>\n",
       "      <td>0.2232</td>\n",
       "      <td>0.061</td>\n",
       "      <td>3.630</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lag1</th>\n",
       "      <td>-0.0384</td>\n",
       "      <td>0.026</td>\n",
       "      <td>-1.466</td>\n",
       "      <td>0.143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Lag2</th>\n",
       "      <td>0.0608</td>\n",
       "      <td>0.027</td>\n",
       "      <td>2.291</td>\n",
       "      <td>0.022</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             coef  std err      z  P>|z|\n",
       "intercept  0.2232    0.061  3.630  0.000\n",
       "Lag1      -0.0384    0.026 -1.466  0.143\n",
       "Lag2       0.0608    0.027  2.291  0.022"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_1=sm.GLM(y.iloc[1:],X.iloc[1:],family=sm.families.Binomial())\n",
    "result_1=model_1.fit()\n",
    "summarize(result_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10cb8b22",
   "metadata": {},
   "source": [
    "### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3c46a189",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['Up'], dtype='<U4')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probs = result_1.predict([X.iloc[0]])\n",
    "labels = np.array(['Down']*1)\n",
    "labels[probs>0.5] = \"Up\"\n",
    "labels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07485ed1",
   "metadata": {},
   "source": [
    "### (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c78722c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "n=len(Weekly)\n",
    "errors=np.zeros(n)\n",
    "for i in range(n):\n",
    "    x_train=X.drop([i])\n",
    "    y_train=y.drop([i])\n",
    "    model_i=sm.GLM(y_train,x_train,family=sm.families.Binomial()).fit()\n",
    "    x_test=X.iloc[i]\n",
    "    pred=model_i.predict([x_test])\n",
    "    labels = np.array(['Down']*1)\n",
    "    labels[pred>0.5] = \"Up\"\n",
    "    if labels[0]!=Weekly.Direction[i]:\n",
    "        errors[i]=1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ab14c4",
   "metadata": {},
   "source": [
    "### (e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e80c8428",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.44995408631772266"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "errors.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f6e9ae2",
   "metadata": {},
   "source": [
    "The test error estimate is 45%. The logistic regression model performs poorly on average (It is slightly better than random guessing). This must be a consequence of the model not enough flexible relatively to the Bayes decision boundary.<br>\n",
    "A possible remedy is to raise the power degrees of the predictors to increase flexibility."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "803495fa",
   "metadata": {},
   "source": [
    "## Exercise 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e04a98c1",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "a59f4975",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng=np.random.default_rng(1)\n",
    "x=rng.normal(size=100)\n",
    "y=x-2*x**2+rng.normal(size=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ce75bea",
   "metadata": {},
   "source": [
    "n is the number of observations and p is the number of predictors. In this case there are two predictors which are $x$ and $x^2$.<br>\n",
    "Hence, n=100 and p=2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a40f4b24",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "f5b81447",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6KElEQVR4nO3de3yU5Z3///ckkglgMhxGmLBGiIoVjFsOHjhYOSgYq3jqA8WWLvwq9KECKyJfLesqIIvU6lZcraAtX1BRsa2ywleXFRRw1VDkENfIlhUaGgxJMYATDpJgMr8/cMYc5nDfk3vmnrnn9Xw88tidyT2Ta4aaeee6PtfncgUCgYAAAAAcJMvuAQAAAFiNgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAABznDLsHYIempiYdOHBAeXl5crlcdg8HAAAYEAgEdPToUfXq1UtZWdHnaDIy4Bw4cECFhYV2DwMAAMRh//79Ovvss6Nek5EBJy8vT9LpNyg/P9/m0QAAACPq6upUWFgY+hyPJiMDTnBZKj8/n4ADAECaMVJeQpExAABwHNsDzqJFi3TppZcqLy9PPXr00E033aTdu3dHfcymTZvkcrnafP35z39O0qgBAEAqsz3gbN68WdOmTdOWLVu0fv16ffPNNxo7dqyOHz8e87G7d+9WdXV16Ktv375JGDEAAEh1ttfgrFu3rsXt5cuXq0ePHtq+fbuuvPLKqI/t0aOHunTpksDRAQCAdGT7DE5rfr9fktStW7eY1w4cOFAFBQW66qqrtHHjxojX1dfXq66ursUXAABwrpQKOIFAQLNmzdIVV1yh4uLiiNcVFBTo+eef1+uvv6433nhD3/ve93TVVVfp/fffD3v9okWL5PF4Ql/0wAEAwNlcgUAgYPcggqZNm6a33npLH3zwQcwGPq2NGzdOLpdLa9asafO9+vp61dfXh24H99H7/X62iQMAkCbq6urk8XgMfX6nzAzOjBkztGbNGm3cuNF0uJGkIUOG6PPPPw/7PbfbHep5Q+8bAACcz/Yi40AgoBkzZmj16tXatGmTioqK4nqenTt3qqCgwOLRAUDyNTYFtLXisA4ePakeebm6rKibsrM4Nw8ww/aAM23aNL3yyit68803lZeXp5qaGkmSx+NRx44dJUlz5sxRVVWVXnzxRUnS4sWL1adPH1100UVqaGjQypUr9frrr+v111+37XUAgBXWlVdr/tpdqvafDN1X4MnV3HH9VVLMH3GAUbYHnCVLlkiSRo4c2eL+5cuXa/LkyZKk6upqVVZWhr7X0NCg2bNnq6qqSh07dtRFF12kt956Sz/84Q+TNWwAsNy68mrdtXKHWhdG1vhP6q6VO7Rk4iBCDmBQShUZJ4uZIiUASIbGpoCueOy9FjM3zbkk+Ty5+uCB0SxXIWOlZZExAGSyrRWHI4YbSQpIqvaf1NaKw8kbFJDGbF+iAgBIB49GDjexrqMoGWiLgAMAKaBHXm5c11GUDITHEhUApIDLirqpwJOrSPMuLp0OLpcVfXeMTbAoufXSVrAoeV15deIGDKQ4Ag4ApIDsLJfmjusvSW1CTvD23HH9Q0tPjU0BzV+7q82OK0mh++av3aXGpozbRwJIIuAAQMooKS7QkomD5PO0XIbyeXLbbBGnKBmIjhocAEghJcUFGtPfF7NouD1FyUAmIOAAQIrJznJp6Hndo14Tb1EykCkIOADQDnZt0Q4WJdf4T4atwwk2BmxelAxkEgIOAMTJzi3awaLku1bukEtqEXLCFSUDmYYiYwCIQyps0TZTlAxkGmZwAMCkWFu0XTq9RXtMf1+7ZlCMLH8ZLUoGMg0BBwBMMrNFO1axcCRmlr+MFCUDmYYlKgAwKdFbtFNh+QtIdwQcADApkVu06VAMWIOAAwAmxXNulFF0KAasQcABAJPMnhtlBh2KAWsQcAAgDonaok2HYsAa7KICgDglYos2HYoBaxBwAKAdWm/RbmwKqHTvobgDj5EOxQ9d15++N0AMBBwAkDVnSll1dENw+av1c/k8ubrh+wVa8JY9x0MA6cQVCAQybq9hXV2dPB6P/H6/8vPz7R4OAJuFCya+/Fzdftk56uPtZCjwBHvXtP6FGnxEPHU5rUPXkeP1mvbKTkt/BpBOzHx+E3AIOEBGixRMWos2S9LYFNAVj70XcXt3sG7mgwdGRwxJsWaQrPgZQLoz8/nNEhWAjBWtqV5rwS7C4WZJ2nt0g5GlrWQcDwE4CdvEAWSsWKGhuWhdhNvTu8bosQz0xwHMIeAAyFhmw0CkLsLx9q4xcywD/XEAcwg4ADLWvtrjcT2udTCK9+gGM8tOiTweAnAiAg6AjLSuvFpPbvg8rse2niWJ9+gGM8tOiTweAnAiAg6AjBNcGjIr2ixJPEc3mF12StTxEKkk2CjxzbIqle49xKnpiBu7qACkHCua7kVjpri4tWizJGaPbojnWIZ4jodI9PtpFasaJQISAQdAiknGh1y8O41+fmVRzDG0Proh1rWxjmUIF6jM/Ix0CQ2R+hFF254PRMMSFYCUYXTLdHvFu9Po99u+sHzJJJHLTsl6P80ItwRlZjcZYBQzOABSQqwPOZdOf8iN6e9r9/JKrKWhSI6cOKVn3vtc91x9Qbt+fmtj+vuUl9tBpXsPSQpo6LleDTmve7teZzLfT6MizSZNuLSQJoawHDM4AFKCmS3TRkUqWI22IymW5R/us3QmYV15ta547D395Hd/0jMb9+iZjXs1+4+faP2umnY9byLez/aINptkdDcbTQxhBjM4AFKC1Z16Y9WeRDqxO5avvj4VdiYhnkLeRNadpFLnYyNLUEbQxBBmEHAApAQrO/UaDQ7NdyTV+L/WQ2+W61h9Y8znr/F/rdK9h5qd8t2gBW+ZK+RN9BKS0fdzX+2JiOOzaudVe3atSeF3kwGxEHAApIR4tkyHYzY4NN+RVHn4az254X9jjnXBW/+jw8cbol4TaxYm0YdnGq0zWrzhf/U935kttp7vqz2hV7dWqqbOmp1XVswS0cQQZlGDAyAlWNWptz21J9NHn68unTrEHGuscBP8OVLk3T8bDNbYxBsOgu+nkSWgOW98quG/fFe3/3aL7llVpic3/G+LcCPFt/MqWAP1+d+OmRz9d1wu6Tc/HsgWcZhGwAGQMqzYMt36gzmScMEhO8ulX95ysbHBGhApTK0rr9ayD/cZeo7ao/WhIumGb5pCRdMf7qnVh5/XRu34W1JcoHuv7htzjEdOnFJNXX3M6yTj27WDxdO3/3aLntm4J+b1EX9uQPr8YPwBCZkrJZaonn32WT3++OOqrq7WRRddpMWLF+sHP/hBxOs3b96sWbNm6bPPPlOvXr10//33684770ziiAEkSjydeoPWlVdrwf/7zNDPiVSjUlJcoKVhio+7de6gw8dPGXsRrTQPU2aOichynV4Oa347UraItITUx9vZ/IAjMLpsFqkGKl7LP9yn6aP7skQFU2wPOK+99ppmzpypZ599VsOHD9dzzz2na6+9Vrt27dI555zT5vqKigr98Ic/1NSpU7Vy5Up9+OGHuvvuu3XWWWfpRz/6kQ2vAIDVzHTqDTL6oWqklidcyKqpO6l7XyszNaag5mHKTMFt6zATbeIkUs1PInYeRVs2i1YDFa9IO9eAaGxfovr1r3+tO+64Q1OmTFG/fv20ePFiFRYWasmSJWGvX7p0qc455xwtXrxY/fr105QpU/Szn/1MTzzxRJJHDiBVDkY0+qFqppYnGLKu//tekqS9B4+aHle4wzkTtS070hJSsNjYyrmPaKHJaICbPuo8vXzH5fJ0jF3zJNEDB+bZOoPT0NCg7du36xe/+EWL+8eOHauPPvoo7GNKS0s1duzYFvddc801WrZsmU6dOqUOHdr+x1JfX6/6+u/Wl+vq6iwYPZDZknXGkZHtykY/VLt1ztHCm4sNjy/cazSrdZhKZC+XcEtI0c67iteR45HrdYwGkb498zS8r1c/G15kaOcaPXBglq0zOLW1tWpsbFTPnj1b3N+zZ0/V1ITfYVBTUxP2+m+++Ua1tbVhH7No0SJ5PJ7QV2FhoTUvAMhQiTrjqPWM0Nv/fSBUqHrPqjLd/tstuuKx99o8v9EP1X++rp+pcBPuNRqV5Qp/OGciZlRaa/1+RCreLvDkqkunDqbHsuCt/4k4W7ev9rih5wgGllg718LNggFG2F6DI0kuV8v/vAKBQJv7Yl0f7v6gOXPmaNasWaHbdXV1hBwgTolqUGd0tiRcrYnRv+59no6GrrOijqQpID3/foUGntO1RciJdYK4FTMs4d6PSMXb63fVmJ7diVRovK682tCxC80DS3Dn2p0rd7S5zsySItCarTM4Xq9X2dnZbWZrDh482GaWJsjn84W9/owzzlD37uEL0Nxut/Lz81t8AYhPIs44MjNbEq7WJDgrEonZWQCjS16d3dkxr5m/dleb7dzRtsM/++OBUV9LNLFeZ7Cu6MYBf6eh3x7mGWkssXy458sWsziNTQHNW2Nsd9hD17UMLMGda61ftxUnqiNz2TqDk5OTo8GDB2v9+vW6+eabQ/evX79eN954Y9jHDB06VGvXrm1x3zvvvKNLLrkkbP0NAGtZfcZRPLMlrWtNsrNcuuH7BXru/YqIjzEzC2B07JOG9tazm/4Sc5w/Wfan0H3N65Raz6gM7t1V2/96RNcW+/R/DfbJCWrPbEfzsXy450s9s3FvzMc8s3GvXt9RFXotWysOG+5B1LVzTtQxxGoPYOUxEnAu25eoZs2apZ/+9Ke65JJLNHToUD3//POqrKwM9bWZM2eOqqqq9OKLL0qS7rzzTj3zzDOaNWuWpk6dqtLSUi1btkyvvvqqnS8DyBhWnhklte+comAQWVdereejhJtwtTDRGB17lsv8JHjrJbbgMs+68mqNeHxj3O+Fr50F3sHZncuKuun1HVUxj3iQWr6W+m+aDP+sSAHSSHuAZBW3I/3Zvk38tttu0+LFi/XII49owIABev/99/X222+rd+/ekqTq6mpVVlaGri8qKtLbb7+tTZs2acCAAVqwYIH+7d/+jR44QJLEKpI1uxzUnu2/PfJyY84AuSSt+aTa1BZ2o68xnr4s4ZbY4i1onnlVXz01YYBenTpEHzww2pIP+GhHZrTW/LV4z3Qb/hnx7ohKVHE7nMn2gCNJd999t/bt26f6+npt375dV155Zeh7K1as0KZNm1pcP2LECO3YsUP19fWqqKigizGQRFadGRUUz4dd8xCViJogo69xyLnd49oR1XxMRpboWr+VBZ5cLZ04SDOu6qseebk6ePS757KCmbqc4GtRQPLlx77el++Oa0dUrOJ2yfgxEsgMKRFwAKQXK86MChrcu6u6hanJiKR1iErUoZVGXqOZ2Y5IYzKyRNcUkB66rl+L2RpJhrbQx6ukuEAfPDBa00edb+j62uP1mndD/5jXzbvhorjqZRIRZOFsttfgAEhP7TkzKihYT2HkdO6g5rUmjU0BrS6rMvS4eGaKjLzGYBCat2aX4SLb5mMyGry8eW7dOODvJEU+liLScQ3xys5yafj5XkOHZfbIO71kt3TiIP3ijU/11YmW53Z16dRBv7zl4rjHZXVxO5yPgAMgbvGcGRVk5OyoAk+uHrquv7p2zgkbMLZWHDZ0AGb3zjlxN4oz+hq/PvWN4edsfh6W0RmHYEBLVB+iSIL1SJGKjluf7RUMhVv+ckilew9JCmjouV4N+Xa3W7ysLm6H8xFwACSdkbqTbp07aPP/GaWcMyKvpBv9a/3GAb0Sto14XXl12CZ1kbReYgsGiGjLL82Lts0s1VhxOGWsxoTNX0vzxww/36vh53vb/fODzAYtgBocAElnpO7k8PFT2v7XI1GvMfrX+pj+PsNjM+N0c7vPTD2mdZ1SsIdPNDd8vyAUIOxYqrGy5ipeVhe3w/mYwQGQdFZ9SMf6q15K7DlGp5vbRT54srnpo87X8PO9bWp4GpsCWvNJ9MLgNZ9U6/6SfsrOctm2VGNFzZUVY1gycVCbPjjt7QEEZyLgAEg6qz6k41k+sZKZWZK+Pc8Mu2RkZDar+ZKTnUs17am5skoqBC2kB5aoACSdkRO1s1zSkeOxZ0fsXD4xM0sS6Vqzs1nBUBcp3EjOX6oJd6YW0BozOACSrvnMSyRNAWnaKzu15NsDIaOx66/6y4q6yZfvjrlMFW2ZLN7ZrC6dOrTZiu1p51ZswEmYwQFgi5LiAv3mx4PadOltrXV32samgEr3HmpxOrdkz1/12VkuzbvhopjXPXRdP22tONxmzJL5oy+C2+tbhxtJ8oe5LxEi/RukyvMBEjM4AGzUtXOOon2Wtd7ynIoHLZYUF0Rsbte1UwfdesnZWvDW/0Qcs5k6IiPb663sgROO1f8G68qrNW/NZy1mwXz5bs274SJTz8cJ42jNFQgEMi4q19XVyePxyO/3Kz8/3+7hIINl+i/lN8uqdM+qspjXPTVhgNxnZEVtDHjv1X01fXRf296/xqaAtuw9pNK/1Eo6PZvkP3FK015pO+bgCJvXCBkJDqV7D+n2326JOZZXpw5JSDFwpOaM4V6P0eeL1kNoqcHnS8Xgi8Qw8/nNDA5gE34pG68/8XZ2a/YfP4k6c/Hkhs/16tb9mneDPe9fdpZLw/t6Nbzv6eZ2jU0BXfHYe4Y7DhupI7LzuAKrOyg3NgX0izc+jXrNL974NObzJevYCqQfanAAGwR/KbfeHhz8pWzVgYmpzmj9iVyKuZVakmrqUuf9i+dwyFh1RHYeV2D1YZdb9h4KW0fU3FcnTmnL3kMRv88J44iGgAMkGb+Uv2O0O23tMWPN9IJS4f1LxGyL2YJkK1n9ek4v5bXvOk4YRzQEHCDJ+KX8ncamgDwdc/Sz4X3UtXOHFt9r3sfGzIxEqrx/iZhtsfO4Autfj9ExRr6OE8YRDTU4QJLxS/m0cDVI3Trn6KYBvTSmv69F/YmRIxlas/v9S1THYbuOKzDyb5Dlkv70l0OGiuaHntddz2zcE/PnRiuW5oRxREPAAZKMX8qRC0OPHG/Q8g/3tflgNNIYsDUj718id7El8hgJOxobGm3OuPjdz0O3oxXNDzm3e9hmhc117dRBQ86NHHA4YRzRsEQFJJmddRSpIN4apNCRDPnuqM9v9P1bV16tKx57T7f/dovuWVWm23+7RVc89p6lBcqJPEbCTGNDqxrpnW7OODBmc8agaEXz2Vku/fKWi6M+ftEtF0d9XZwwjmjog0MfHNggOIMhhf/L3slbW9vby6WxKaBn3tujJzf8b5vvGX3/rO7nEku8M0VWzDBZ3Y7A6L9fUHAW5YMHRocd++lGf7tUUxf/+Gi5kDnogwOkOLvqKFJBe2uQsrNcuufqvvqe78y43j+r+7kYEc8p3FZ8aMfbIyZasDJb29S6G3VrViy3ccI4wiHgADbJ1F/KVtUgxfv+mdnFlohuwEZY0bwu3iAXK1jFWxsWLRjFEwAT8RxwFgIOYKNM/KVsZWFoPO9fqu9is2qGKZ4gZyRYjenvM72jTXJ20TxSE0XGQBpxwqnLdheGpvouNqv6JJkNckaLvyVF/PcLx+lF80hdzOAAacJJhZR21iCl+tZiq2aYzAY5M8Eq0r9fa+xkgp0IOEAacOKBgnbVICWyP40VrJphMhvkzAar1v9++2qP69Wtlaqp++5YDTOBNZE9iVL5ZyNxCDhAirNj10+yZGe5dFlRt9CHy9aKw0n5cEnlXWxWzTCZDXLxBKvWNVDTR/eNKyiEm53s0rGD/r/hfTR9dN+E/u/BSTOjaIk+OPTBQYprb9+YVBT8i3n9rhr9e9kBHT7eEPqeFR8uRv8iT9W/3K3sk2T0A7yxKaArHnsvZrCK1M8mXpFmJ4O6dOqgX95ycULCRrL7IaH9zHx+E3AIOEhxb5ZV6Z5VZTGve2rCAN044O8SNg6rwkC4D9zm2vvh4pS/yK18HUb/7ZLdgDIYqqLV8QR/frJ/dqICHdqHRn+Ag6TCrh+rPmxj/bUutW/ZzUm1SlbWKBndTp/spbtYhc1BAVm/DJsO/ZDQPgQcIMXZvevHqtAQrZaotXg+XJxYq2RHn6RkFn+b6TVkddhI9X5IaD/64AApLpl9Y1r32Wn4pimugzHDMfrXenNmPlys6h8Dcwd5tofZWUcrw0YqzIwisZjBAdJAMpYOwi1Ddeuc06IAuLXWMy1WnmEkmftw4S/y9BOcnTQafK0MG3bPjCLxCDhAmkjk0kGkZaho4aa5g0dPWnqGUTwfLvxFnn6Cs5N3flvYHEkiwkaq90NC+7FEBaSRRCwdmKmNiWRf7QndtXJHm7/Eg3U668qrQ38xxxpxvB8usZ6fIwNSU0lxgZZOHKQunTqE/X4iw0ZwZtTnaRl6fZ7ctCpIR3hsE2ebODKc0T474bgk9cx3S3Kppi72dtv1u2rCbkNurj1bupO9zRnWaWwK6Jn39mj5hxX66utTofuTscU/VfshoS364MRAwAG+Y7TPTmvBX/8zr+6rJzd8HvP6YCPC8LU+HXTzgL/T1f197f5wcUofnExF2EA09MEBYJjRmpRunTvo8PHv/rIOFjjXf9Nk6PGRzjCy+kPMrjOuYA07tsbDmWwLOPv27dOCBQv03nvvqaamRr169dLEiRP14IMPKicnJ+LjJk+erBdeeKHFfZdffrm2bIlvih3IdEZ3k2z+P6O0/a9H2oSG0r2HDP2caGcYWY0PSQC2BZw///nPampq0nPPPafzzz9f5eXlmjp1qo4fP64nnngi6mNLSkq0fPny0O1ogQhAdEZ3k+SckRU2NLDdFkAqsi3glJSUqKSkJHT73HPP1e7du7VkyZKYAcftdsvn8yV6iEDGaE+fHbbbAkhFKVWD4/f71a1b7L/yNm3apB49eqhLly4aMWKEFi5cqB49eiRhhIBztad2JdlnGAFALCkTcPbu3aunn35a//qv/xr1umuvvVbjx49X7969VVFRoYceekijR4/W9u3b5Xa7wz6mvr5e9fX1odt1dXWWjh0Axb0AUovl28TnzZun+fPnR73m448/1iWXXBK6feDAAY0YMUIjRozQ7373O1M/r7q6Wr1799aqVat0yy23mBoT28ThdGa23LK9GkCqs7UPTm1trWpra6Ne06dPH+Xmnt5RceDAAY0aNUqXX365VqxYoaws882V+/btqylTpuiBBx4I+/1wMziFhYUEHDiamcAS6agGGuQBSCW29sHxer3yer2Grq2qqtKoUaM0ePBgLV++PK5wc+jQIe3fv18FBZF/+brd7ojLV4ATRQoswaMTmgeWaEc1BHQ65Mxfu0tj+vtYbgKQNmw7i+rAgQMaOXKkCgsL9cQTT+jLL79UTU2NampqWlx34YUXavXq1ZKkY8eOafbs2SotLdW+ffu0adMmjRs3Tl6vVzfffLMdLwNIObECi3Q6sDQ2nb61teJw1NOcm58YDgDpwrYi43feeUd79uzRnj17dPbZZ7f4XvNVs927d8vv90uSsrOz9emnn+rFF1/UV199pYKCAo0aNUqvvfaa8vLykjp+IFWZCSxDz+se6jAci9HrACAV2BZwJk+erMmTJ8e8rnnY6dixo/7zP/8zgaMC0p/ZwGL0qAbvmSzzAkgfti1RAYitsSmg0r2H9GZZlUr3HgotK0VjNLAErwt2Io5VXXPf78u0rrza0HMDgN1Spg8OgJbi3bZt9uiEaJ2Im/tbXX2bAuVUx8nUQOayfJt4OjCzzQywQ3u3bQcfL4U/OiHc49eVV2veml2qqYu8xBUMRx88MDrlgwJ9fQDnMfP5zRIVkGLM7oIKJ3h0gs/TcrnK58mNGI5Kigv0r+O/H3Vs6bKjKhjwWhdbB7fJs9QGOB9LVECKMbsLKpJ4jk6oPV4f8XvNpfKOKvr6AJAIOEDKsXLbdnaWK2oIas1sgXIqsiogAkhvLFEBKcbOkBFrR5VLp+tYggXKqYi+PgAkAg6QcoIhI5Yjxxss/9nBHVWS2oSc4O254/qn9NKOE2ahALQfAQdIMdlZLj10Xb+Y1y14K3qhcbziKVBOJU6YhQLQftTgACmoa+fYXYPD1ZFY1fclngLlVBGtr0+6zEIBaD8CDpCC4qkjsbrvi9kC5VQSnIVq/X746IMDZAwCDpCCzNaRRGoMGOz7kg5LS1ZL51koAO1HwAFSkJnjFuj7Elk6z0IBaB+KjIEUZGY3k5m+LwCQKQg4QIoyupuJvi8A0BZLVEAKM1JHQt8XAGiLgAOkuFh1JGbqdQAgU7BEBaQ5J3QfBgCrEXAAB0j37sMAYDWWqIAEsKqjsBn0fQGA7xBwAItZ3VHYDPq+AMBpLFEBFgp2FG7dlybYUXhdebVNIwOAzELAASwSq6OwdLqjcCJOAAcAtETAASxCR2EASB0EHMAidBQGgNRBwAEsQkdhAEgdBBzAIsGOwpE2Zbt0ejcVHYUBIPEIOIBF6CgMAKmDgANYiI7CAJAaaPQHWIyOwoBz2NGVHNYg4AAJQEdhIP3Z2ZUc7ccSFQAArdCVPP0RcAAAaIau5M5AwAEAoBm6kjsDAQcAgGboSu4MBBwAAJqhK7kzEHAAAGiGruTOQMABDGpsCqh07yG9WVal0r2HKDAEHIqu5M5AHxzAAPphAJkl2JW89X/3Pv67Txu2zuD06dNHLperxdcvfvGLqI8JBAKaN2+eevXqpY4dO2rkyJH67LPPkjRiZCL6YQCZqaS4QB88MFqvTh2ipyYM0KtTh+iDB0YTbtKE7TM4jzzyiKZOnRq6feaZZ0a9/le/+pV+/etfa8WKFbrgggv0L//yLxozZox2796tvLy8RA8XGSZWPwyXTvfDGNPfx3Q14EB0JU9fttfg5OXlyefzhb6iBZxAIKDFixfrwQcf1C233KLi4mK98MILOnHihF555ZUkjhqZgn4YAJCebA84jz32mLp3764BAwZo4cKFamhoiHhtRUWFampqNHbs2NB9brdbI0aM0EcffRTxcfX19aqrq2vxBRhBPwwASE+2LlHdc889GjRokLp27aqtW7dqzpw5qqio0O9+97uw19fU1EiSevbs2eL+nj176q9//WvEn7No0SLNnz/fuoEjY9APAwDSk+UzOPPmzWtTONz6a9u2bZKke++9VyNGjNDf//3fa8qUKVq6dKmWLVumQ4cORf0ZLlfLWodAINDmvubmzJkjv98f+tq/f3/7XygyAv0wACA9WT6DM336dE2YMCHqNX369Al7/5AhQyRJe/bsUffubYu6fD6fpNMzOQUF31WxHzx4sM2sTnNut1tutzvW0IE2gv0w7lq5Qy6pRbFxvP0wGpsC2lpxWAePnlSPvNPhiAJlALCW5QHH6/XK6/XG9didO3dKUovw0lxRUZF8Pp/Wr1+vgQMHSpIaGhq0efNmPfbYY/ENGIjByn4Y9NMBgOSwrQantLRUW7Zs0ahRo+TxePTxxx/r3nvv1Q033KBzzjkndN2FF16oRYsW6eabb5bL5dLMmTP16KOPqm/fvurbt68effRRderUST/+8Y/teinIACXFBRrT39eumZdgP53WW86D/XSWTBxEyAEAi9gWcNxut1577TXNnz9f9fX16t27t6ZOnar777+/xXW7d++W3+8P3b7//vv19ddf6+6779aRI0d0+eWX65133qEHDhKuPf0w6KcDAMnlCgQCGXegTl1dnTwej/x+v/Lz8+0eDjJA6d5Duv23W2Je9+rUITQVA4AIzHx+294HB8gE9NMBgOQi4ABJQD8dAEguAg6QBPTTAYDkIuAASRDspyOpTciJt58OACAyAg6QJMF+Oj5Py2UonyeXLeIAYDFbz6ICMo0V/XQAALERcIAka08/HQCAMSxRAQAAx2EGBwAAmJIOhwYTcAAAgGHpcmgwS1QAAMCQ4KHBzcON9N2hwevKq20aWVsEHAAAEFOsQ4Ol04cGNzalxhGXBBwAABDT1orDbWZumgtIqvaf1NaKw8kbVBQEHAAAEFO6HRpMwAEAADGl26HBBBwAABBTuh0aTMABAAAxpduhwQQcAABgSDodGkyjPyBJ0qHzJwDEki6HBhNwgCRIl86fAGBEOhwazBIVkGDp1PkTAJyCgAMkULp1/gQApyDgAAmUbp0/AcApCDhAAqVb508AcAoCDpBA6db5EwCcgoADJFC6df4EAKcg4AAJlG6dPwHAKQg4QIKlU+dPAHAKGv0BSZAunT8BwCkIOEAYiThWIR06fwKAUxBwgFY4VgEA0h81OEAzHKsAAM5AwAG+xbEKANB+jU0Ble49pDfLqlS695BtvzNZogK+ZeZYBWppAKCtVFriZwYH+BbHKgBA/FJtiZ+AA3yLYxUAID6puMRPwAG+xbEKABAfM0v8yULAAb7FsQoAEJ9UXOK3LeBs2rRJLpcr7NfHH38c8XGTJ09uc/2QIUOSOHI4GccqAIB5qbjEb9suqmHDhqm6umXB0UMPPaQNGzbokksuifrYkpISLV++PHQ7JycnIWNEZuJYBQAwJ7jEX+M/GbYOx6XTfygmc4nftoCTk5Mjn88Xun3q1CmtWbNG06dPl8sV/YPE7Xa3eCxgNY5VAADjgkv8d63cIZfUIuTYtcSfMjU4a9asUW1trSZPnhzz2k2bNqlHjx664IILNHXqVB08eDDxA0RaSZVGUwCQKVJtid8VCARS4jf/D3/4Q0nS22+/HfW61157TWeeeaZ69+6tiooKPfTQQ/rmm2+0fft2ud3usI+pr69XfX196HZdXZ0KCwvl9/uVn59v3YtASkilRlMAkGkScVhxUF1dnTwej6HPb8sDzrx58zR//vyo13z88cct6my++OIL9e7dW7///e/1ox/9yNTPq66uVu/evbVq1SrdcsstpsZEwHGeYKOp1v+jDv6nRaEwAKQvWwNObW2tamtro17Tp08f5eZ+N4W1YMECPf3006qqqlKHDh1M/8y+fftqypQpeuCBB8J+nxmczNDYFNAVj70XsRdDsMjtgwdGUzAMAGnITMCxvMjY6/XK6/Uavj4QCGj58uX6h3/4h7jCzaFDh7R//34VFET+q9ztdkdcvoJzcJYUACDI9iLj9957TxUVFbrjjjvCfv/CCy/U6tWrJUnHjh3T7NmzVVpaqn379mnTpk0aN26cvF6vbr755mQOGynIaAOp9btqEjwSAIDdbA84y5Yt07Bhw9SvX7+w39+9e7f8fr8kKTs7W59++qluvPFGXXDBBZo0aZIuuOAClZaWKi8vL5nDRgoy2kDq/364L+mHvgEAkitldlElk5k1PKSPWDU4QdTiAEB6MvP5bfsMDmCV5mdJRWPHoW8AgOQi4CBjJfPQNwBAchFw4BiNTQHNX7vL8PXJPPQNAJBcBBw4Rqxt4s0VJPnQNwBAchFw4BhmlpySfegbACC5CDhwDKNLTvdefQHHNQCAwxFw4BiXFXVTgSdX0eZlfPluTR99ftLGBACwBwEHjtF8m3jrkOP69mveDRexNAUAGYCAA0cpKS7QkomD5PO0XK7yeXI5SRwAMojlh20CdispLtCY/j5trTisg0dPqkfe6R1TzNwAQOYg4MCRsrNcnBgOABmMJSoAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4BBwAAOA4CQ04Cxcu1LBhw9SpUyd16dIl7DWVlZUaN26cOnfuLK/Xq3/8x39UQ0ND1Oetr6/XjBkz5PV61blzZ91www364osvEvAKkGiNTQGV7j2kN8uqVLr3kBqbAnYPCQDgAGck8skbGho0fvx4DR06VMuWLWvz/cbGRl133XU666yz9MEHH+jQoUOaNGmSAoGAnn766YjPO3PmTK1du1arVq1S9+7ddd999+n666/X9u3blZ2dnciXBAutK6/W/LW7VO0/GbqvwJOrueP6q6S4wMaRAQDSnSsQCCT8T+YVK1Zo5syZ+uqrr1rc/x//8R+6/vrrtX//fvXq1UuStGrVKk2ePFkHDx5Ufn5+m+fy+/0666yz9NJLL+m2226TJB04cECFhYV6++23dc0118QcT11dnTwej/x+f9ifgcRbV16tu1buUOv/8bm+/b9LJg4i5AAAWjDz+W1rDU5paamKi4tD4UaSrrnmGtXX12v79u1hH7N9+3adOnVKY8eODd3Xq1cvFRcX66OPPgr7mPr6etXV1bX4gn0amwKav3ZXm3AjKXTf/LW7WK4CAMTN1oBTU1Ojnj17triva9euysnJUU1NTcTH5OTkqGvXri3u79mzZ8THLFq0SB6PJ/RVWFhozQtAXLZWHG6xLNVaQFK1/6S2VhxO3qAAAI5iOuDMmzdPLpcr6te2bdsMP5/L5WpzXyAQCHt/NNEeM2fOHPn9/tDX/v37TT03rHXwaORwE891AAC0ZrrIePr06ZowYULUa/r06WPouXw+n/70pz+1uO/IkSM6depUm5md5o9paGjQkSNHWsziHDx4UMOGDQv7GLfbLbfbbWhMSLweebmWXgcAQGumA47X65XX67Xkhw8dOlQLFy5UdXW1CgpOF5S+8847crvdGjx4cNjHDB48WB06dND69et16623SpKqq6tVXl6uX/3qV5aMC4l1WVE3FXhyVeM/GbYOxyXJ58nVZUXdkj00AIBDJLQGp7KyUmVlZaqsrFRjY6PKyspUVlamY8eOSZLGjh2r/v3766c//al27typd999V7Nnz9bUqVND1dFVVVW68MILtXXrVkmSx+PRHXfcofvuu0/vvvuudu7cqYkTJ+riiy/W1VdfnciXA4tkZ7k0d1x/Sd/tmgoK3p47rr+ys8wtUwIAEJTQPjgPP/ywXnjhhdDtgQMHSpI2btyokSNHKjs7W2+99ZbuvvtuDR8+XB07dtSPf/xjPfHEE6HHnDp1Srt379aJEydC9z355JM644wzdOutt+rrr7/WVVddpRUrVtADJ42UFBdoycRBbfrg+OiDAwCwQFL64KQa+uCkjsamgLZWHNbBoyfVI+/0shQzNwCAcMx8fid0BgeIJTvLpaHndbd7GAAAh+GwTQAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgEHAAA4DgJDTgLFy7UsGHD1KlTJ3Xp0qXN9z/55BPdfvvtKiwsVMeOHdWvXz899dRTMZ935MiRcrlcLb4mTJiQgFcAAADS0RmJfPKGhgaNHz9eQ4cO1bJly9p8f/v27TrrrLO0cuVKFRYW6qOPPtLPf/5zZWdna/r06VGfe+rUqXrkkUdCtzt27Gj5+AEAQHpKaMCZP3++JGnFihVhv/+zn/2sxe1zzz1XpaWleuONN2IGnE6dOsnn81kyTgAA4CwpV4Pj9/vVrVu3mNe9/PLL8nq9uuiiizR79mwdPXo0CaMDAADpIKEzOGaVlpbq97//vd56662o1/3kJz9RUVGRfD6fysvLNWfOHH3yySdav3592Ovr6+tVX18ful1XV2fpuAEAQGoxPYMzb968NgW+rb+2bdtmeiCfffaZbrzxRj388MMaM2ZM1GunTp2qq6++WsXFxZowYYL++Mc/asOGDdqxY0fY6xctWiSPxxP6KiwsND0+AACQPlyBQCBg5gG1tbWqra2Nek2fPn2Um5sbur1ixQrNnDlTX331Vdjrd+3apVGjRmnKlClauHChmeFIkgKBgNxut1566SXddtttbb4fbgansLBQfr9f+fn5pn8eAABIvrq6Onk8HkOf36aXqLxer7xeb9yDa+2zzz7T6NGjNWnSpLjCTfA5Tp06pYKCgrDfd7vdcrvd7RkmAABIIwktMq6srFRZWZkqKyvV2NiosrIylZWV6dixY5JOB5NRo0ZpzJgxmjVrlmpqalRTU6Mvv/wy9BxVVVW68MILtXXrVknS3r179cgjj2jbtm3at2+f3n77bY0fP14DBw7U8OHDE/lyAABAmkhokfHDDz+sF154IXR74MCBkqSNGzdq5MiR+sMf/qAvv/xSL7/8sl5++eXQdb1799a+ffskSadOndLu3bt14sQJSVJOTo7effddPfXUUzp27JgKCwt13XXXae7cucrOzk7kywEAAGnCdA2OE5hZwwMAAKnBzOd3yvXBAQAAaC8CDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcJwz7B4A7NHYFNDWisM6ePSkeuTl6rKibsrOctk9LAAALEHAyUDryqs1f+0uVftPhu4r8ORq7rj+KikusHFkAABYgyWqDLOuvFp3rdzRItxIUo3/pO5auUPryqttGhkAANYh4GSQxqaA5q/dpUCY7wXvm792lxqbwl0BAED6IOBkkK0Vh9vM3DQXkFTtP6mtFYeTNygAABKAgJNBDh6NHG7iuQ4AgFRFwMkgPfJyLb0OAIBURcDJIJcVdVOBJ1eRNoO7dHo31WVF3ZI5LAAALEfAySDZWS7NHddfktqEnODtueP60w8HAJD2CDgZpqS4QEsmDpLP03IZyufJ1ZKJg+iDAwBwBBr9ZaCS4gKN6e+jkzEAwLEIOBkqO8uloed1t3sYAAAkBEtUAADAcQg4AADAcQg4AADAcQg4AADAcQg4AADAcdhFZaHGpgBbrwEASAEEHIusK6/W/LW7WpzWXeDJ1dxx/WmeBwBAkrFEZYF15dW6a+WOFuFGkmr8J3XXyh1aV15t08gAAMhMBJx2amwKaP7aXQqE+V7wvvlrd6mxKdwVAAAgERIacBYuXKhhw4apU6dO6tKlS9hrXC5Xm6+lS5dGfd76+nrNmDFDXq9XnTt31g033KAvvvgiAa8gtq0Vh9vM3DQXkFTtP6mtFYeTNygAADJcQgNOQ0ODxo8fr7vuuivqdcuXL1d1dXXoa9KkSVGvnzlzplavXq1Vq1bpgw8+0LFjx3T99dersbHRyuEbcvBo5HATz3UAAKD9ElpkPH/+fEnSihUrol7XpUsX+Xw+Q8/p9/u1bNkyvfTSS7r66qslSStXrlRhYaE2bNiga665pl1jNqtHXm7si0xcBwAA2i8lanCmT58ur9erSy+9VEuXLlVTU1PEa7dv365Tp05p7Nixoft69eql4uJiffTRR2EfU19fr7q6uhZfVrmsqJsKPLmKtBncpdO7qS4r6mbZzwQAANHZHnAWLFigP/zhD9qwYYMmTJig++67T48++mjE62tqapSTk6OuXbu2uL9nz56qqakJ+5hFixbJ4/GEvgoLCy0bf3aWS3PH9ZekNiEneHvuuP70wwEAIIlMB5x58+aFLQxu/rVt2zbDz/fP//zPGjp0qAYMGKD77rtPjzzyiB5//HGzw1IgEJDLFT5EzJkzR36/P/S1f/9+088fTUlxgZZMHCSfp+UylM+TqyUTB9EHBwCAJDNdgzN9+nRNmDAh6jV9+vSJdzwaMmSI6urq9Le//U09e/Zs832fz6eGhgYdOXKkxSzOwYMHNWzYsLDP6Xa75Xa74x6TESXFBRrT30cnYwAAUoDpgOP1euX1ehMxFknSzp07lZubG3Fb+eDBg9WhQwetX79et956qySpurpa5eXl+tWvfpWwcRmRneXS0PO62zoGAACQ4F1UlZWVOnz4sCorK9XY2KiysjJJ0vnnn68zzzxTa9euVU1NjYYOHaqOHTtq48aNevDBB/Xzn/88NONSVVWlq666Si+++KIuu+wyeTwe3XHHHbrvvvvUvXt3devWTbNnz9bFF18c2lUFAAAyW0IDzsMPP6wXXnghdHvgwIGSpI0bN2rkyJHq0KGDnn32Wc2aNUtNTU0699xz9cgjj2jatGmhx5w6dUq7d+/WiRMnQvc9+eSTOuOMM3Trrbfq66+/1lVXXaUVK1YoOzs7kS8HAACkCVcgEMi4MwTq6urk8Xjk9/uVn59v93AAAIABZj6/bd8mDgAAYDUCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcJyE9sFJVcGd8VaeKg4AABIr+LltpMNNRgaco0ePSpKlp4oDAIDkOHr0qDweT9RrMrLRX1NTkw4cOKC8vLyIJ5Cjferq6lRYWKj9+/fTTDGJeN+Tj/c8+XjP7ZEK73sgENDRo0fVq1cvZWVFr7LJyBmcrKwsnX322XYPIyPk5+fzC8gGvO/Jx3uefLzn9rD7fY81cxNEkTEAAHAcAg4AAHAcAg4Swu12a+7cuXK73XYPJaPwvicf73ny8Z7bI93e94wsMgYAAM7GDA4AAHAcAg4AAHAcAg4AAHAcAg4AAHAcAg4Sbt++fbrjjjtUVFSkjh076rzzztPcuXPV0NBg99AcbeHChRo2bJg6deqkLl262D0cR3r22WdVVFSk3NxcDR48WP/1X/9l95Ac7/3339e4cePUq1cvuVwu/fu//7vdQ3K0RYsW6dJLL1VeXp569Oihm266Sbt377Z7WIYQcJBwf/7zn9XU1KTnnntOn332mZ588kktXbpU//RP/2T30BytoaFB48eP11133WX3UBzptdde08yZM/Xggw9q586d+sEPfqBrr71WlZWVdg/N0Y4fP67vf//7euaZZ+weSkbYvHmzpk2bpi1btmj9+vX65ptvNHbsWB0/ftzuocXENnHY4vHHH9eSJUv0l7/8xe6hON6KFSs0c+ZMffXVV3YPxVEuv/xyDRo0SEuWLAnd169fP910001atGiRjSPLHC6XS6tXr9ZNN91k91AyxpdffqkePXpo8+bNuvLKK+0eTlTM4MAWfr9f3bp1s3sYQFwaGhq0fft2jR07tsX9Y8eO1UcffWTTqIDE8/v9kpQWv78JOEi6vXv36umnn9add95p91CAuNTW1qqxsVE9e/ZscX/Pnj1VU1Nj06iAxAoEApo1a5auuOIKFRcX2z2cmAg4iNu8efPkcrmifm3btq3FYw4cOKCSkhKNHz9eU6ZMsWnk6Sue9xyJ43K5WtwOBAJt7gOcYvr06frv//5vvfrqq3YPxZAz7B4A0tf06dM1YcKEqNf06dMn9P8fOHBAo0aN0tChQ/X8888neHTOZPY9R2J4vV5lZ2e3ma05ePBgm1kdwAlmzJihNWvW6P3339fZZ59t93AMIeAgbl6vV16v19C1VVVVGjVqlAYPHqzly5crK4vJw3iYec+RODk5ORo8eLDWr1+vm2++OXT/+vXrdeONN9o4MsBagUBAM2bM0OrVq7Vp0yYVFRXZPSTDCDhIuAMHDmjkyJE655xz9MQTT+jLL78Mfc/n89k4MmerrKzU4cOHVVlZqcbGRpWVlUmSzj//fJ155pn2Ds4BZs2apZ/+9Ke65JJLQrOSlZWV1JYl2LFjx7Rnz57Q7YqKCpWVlalbt24655xzbByZM02bNk2vvPKK3nzzTeXl5YVmLT0ejzp27Gjz6GIIAAm2fPnygKSwX0icSZMmhX3PN27caPfQHOM3v/lNoHfv3oGcnJzAoEGDAps3b7Z7SI63cePGsP+7njRpkt1Dc6RIv7uXL19u99Biog8OAABwHAohAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4/z/mSjGHVobXDUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x,y);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eff9d1f",
   "metadata": {},
   "source": [
    "The scatter shows a quadratic plot. \n",
    "The two points in the bottom left corner could be outsiders."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97675a2c",
   "metadata": {},
   "source": [
    "### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0e1c3de4",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "219662b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.63302984, 1.12293686, 1.30179655, 1.33239427])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n=len(x)\n",
    "x.reshape(1,-1)\n",
    "y.reshape(1,-1)\n",
    "cv_error = np.zeros(4)\n",
    "M = sklearn_sm(sm.OLS)\n",
    "for i, d in enumerate(range(1,5)):\n",
    "    X = np.power.outer(x, np.arange(d+1))\n",
    "    M_CV = cross_validate(M,X,y,cv=n)\n",
    "    cv_error[i] = np.mean(M_CV['test_score'])\n",
    "cv_error"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75348424",
   "metadata": {},
   "source": [
    "The model reported in ii returns the smallest CV error."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6664a9a",
   "metadata": {},
   "source": [
    "### (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "358947fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "cc131adf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.63302984, 1.12293686, 1.30179655, 1.33239427])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cv_error2 = np.zeros(4)\n",
    "x.reshape(1,-1)\n",
    "y.reshape(1,-1)\n",
    "M = sklearn_sm(sm.OLS)\n",
    "for i, d in enumerate(range(1,5)):\n",
    "    X = np.power.outer(x, np.arange(d+1))\n",
    "    M_CV = cross_validate(M,X,y,cv=n)\n",
    "    cv_error2[i] = np.mean(M_CV['test_score'])\n",
    "cv_error2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a2e015",
   "metadata": {},
   "source": [
    "The results are exactly the same because we only remove one observation from the training set. Thus, there is no random effect resulting from the observations used for the test set. LOOCV will always be the same, no matter the random seed is."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee20b77b",
   "metadata": {},
   "source": [
    "### (e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a5e6091",
   "metadata": {},
   "source": [
    "The smallest error is found for the quadratic model ii. This result was expected since that model has the same form of the true one defined in question a."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6870e8d5",
   "metadata": {},
   "source": [
    "### (f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "8ef4ff14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>std err</th>\n",
       "      <th>t</th>\n",
       "      <th>P&gt;|t|</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>poly(x, degree=4)[0]</th>\n",
       "      <td>16.6000</td>\n",
       "      <td>1.911</td>\n",
       "      <td>8.685</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>poly(x, degree=4)[1]</th>\n",
       "      <td>-22.2244</td>\n",
       "      <td>1.911</td>\n",
       "      <td>-11.627</td>\n",
       "      <td>0.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>poly(x, degree=4)[2]</th>\n",
       "      <td>-1.0768</td>\n",
       "      <td>1.911</td>\n",
       "      <td>-0.563</td>\n",
       "      <td>0.575</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>poly(x, degree=4)[3]</th>\n",
       "      <td>2.2726</td>\n",
       "      <td>1.911</td>\n",
       "      <td>1.189</td>\n",
       "      <td>0.237</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                         coef  std err       t  P>|t|\n",
       "poly(x, degree=4)[0]  16.6000    1.911   8.685  0.000\n",
       "poly(x, degree=4)[1] -22.2244    1.911 -11.627  0.000\n",
       "poly(x, degree=4)[2]  -1.0768    1.911  -0.563  0.575\n",
       "poly(x, degree=4)[3]   2.2726    1.911   1.189  0.237"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df=pd.DataFrame({'x':x,'y':y})\n",
    "X= MS([poly('x', degree=4)],intercept=False).fit_transform(df)\n",
    "model=sm.OLS(y,X)\n",
    "result=model.fit()\n",
    "summarize(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144ea991",
   "metadata": {},
   "source": [
    "As it was expected, the pvalues show that the only significant coefficients are those related to $x$ and $x^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9193ebc2",
   "metadata": {},
   "source": [
    "## Exercise 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "7ec727fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "Boston=load_data('Boston')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2c0fe51",
   "metadata": {},
   "source": [
    "### (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "5d009aa7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22.532806324110677"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_hat=Boston['medv'].mean()\n",
    "mu_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b884f01e",
   "metadata": {},
   "source": [
    "The population mean estimate for medv is 22.53."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b30db23",
   "metadata": {},
   "source": [
    "### (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "97c69dfa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4088611474975351"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SE_mu_hat=Boston['medv'].std()/(Boston.shape[0])**0.5\n",
    "SE_mu_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50a53349",
   "metadata": {},
   "source": [
    "The obtained estimation of the standard deviation is 0.41."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "058ba61e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RangeIndex(start=0, stop=506, step=1)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Boston['medv'].index"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb62b399",
   "metadata": {},
   "source": [
    "### (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "9e5c774f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4047759328720298"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "means=[]\n",
    "for _ in range(1000):\n",
    "    means.append(Boston['medv'].sample(n=Boston.shape[0],replace=True,random_state=None).mean())\n",
    "boot_std=np.std(means)\n",
    "boot_std"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d58b9ea",
   "metadata": {},
   "source": [
    "The answer is approximately the same as the one obtained in (b)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb47d6d4",
   "metadata": {},
   "source": [
    "### (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "f89a1c5e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The 95% confidence interval for the mean of medv is 21.72325445836662 23.342358189854735\n"
     ]
    }
   ],
   "source": [
    "print('The 95% confidence interval for the mean of medv is', mu_hat-2*boot_std,mu_hat+2*boot_std)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff10c71",
   "metadata": {},
   "source": [
    "### (e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "b6590552",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.2"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_med=Boston['medv'].median()\n",
    "mu_med"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7482bb05",
   "metadata": {},
   "source": [
    "### (f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "39ae668b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3806817673267263"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "medians=[]\n",
    "for _ in range(10000):\n",
    "    medians.append(Boston['medv'].sample(n=Boston.shape[0],replace=True).median())\n",
    "boot_med_std=np.std(medians)\n",
    "boot_med_std"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "562c3e3d",
   "metadata": {},
   "source": [
    "The standard error of the median is less than the standard error of the mean. This a sign of non normality of the distribution of feature 'medv'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "b5dcf88a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.75"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_01=np.percentile(Boston['medv'],10)\n",
    "mu_01"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14368a54",
   "metadata": {},
   "source": [
    "### (g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "edefd3ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.500366615593007"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "percentiles=[]\n",
    "for _ in range(5000):\n",
    "    sample=Boston['medv'].sample(n=Boston.shape[0],replace=True)\n",
    "    percentiles.append(np.percentile(sample,10))\n",
    "boot_mu01_std=np.std(percentiles)\n",
    "boot_mu01_std"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a405e2a",
   "metadata": {},
   "source": [
    "The standard error is small relatively to the tenth percentile value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e66072f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}