<jupyter_start><jupyter_text>Chapter 5: Resampling Methods- **Chapter 5 from the book [An Introduction to Statistical Learning](https://www.statlearning.com/).**- **By <NAME>, <NAME>, <NAME> and <NAME>.**- **Pages from $198$ to $201$**- **By [<NAME>](https://www.linkedin.com/in/mosta-ashour/)****Exercises:**- **[5.](5) [(a)](5a) [(b)](5b) [(c)](5c) [(d)](5d)**- **[6.](6) [(a)](6a) [(b)](6b) [(c)](6c) [(d)](6d)**- **[7.](7) [(a)](7a) [(b)](7b) [(c)](7c) [(d)](7d)**- **[8.](8) [(a)](8a) [(b)](8b) [(c)](8c) [(d)](8d) [(e)](8e) [(f)](8f)**- **[9.](9) [(a)](9a) [(b)](9b) [(c)](9c) [(d)](9d) [(e)](9e) [(f)](9f) [(g)](9g) [(h)](9h)**    5.4 Exercises   Applied<jupyter_code>import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scripts.ch4 import *

import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split

from IPython.display import display
from sklearn import datasets
from scipy import stats

# import warnings filter
from warnings import simplefilter
# ignore all future warnings
simplefilter(action='ignore', category=FutureWarning)
def confusion_table(confusion_mtx):
    """Renders a nice confusion table with labels"""
    confusion_df = pd.DataFrame({'y_pred=0': np.append(confusion_mtx[:, 0], confusion_mtx.sum(axis=0)[0]),
                                 'y_pred=1': np.append(confusion_mtx[:, 1], confusion_mtx.sum(axis=0)[1]),
                                 'Total': np.append(confusion_mtx.sum(axis=1), ''),
                                 '': ['y=0', 'y=1', 'Total']}).set_index('')
    return confusion_df

def total_error_rate(confusion_mtx):
    """Derive total error rate from confusion matrix"""
    return 1 - np.trace(confusion_mtx) / np.sum(confusion_mtx)<jupyter_output><empty_output><jupyter_text>$5.$ In Chapter $4$, we used logistic regression to predict the probability of $\color{brown}{\text{default}}$ using $\color{brown}{\text{income}}$ and $\color{brown}{\text{balance}}$ on the $\color{brown}{\text{Default}}$ data set. We will now estimate the test error of this logistic regression model using the validation set approach. Do not forget to set a random seed before beginning your analysis.<jupyter_code># Read the Default dataset
default_df = pd.read_csv('datasets/Default.csv').drop('Unnamed: 0', axis=1)

# Get all dummies in the df
default_df = pd.get_dummies(default_df, dtype=np.float64, drop_first=True)
default_df.head()<jupyter_output><empty_output><jupyter_text>**$(a)$ Fit a logistic regression model that uses $\color{brown}{\text{income}}$ and $\color{brown}{\text{balance}}$ to predict $\color{brown}{\text{default}}$.**<jupyter_code># Logistic regresssion
X = default_df[['income', 'balance']]
y = default_df.default_Yes
log_r = LogisticRegression().fit(X, y)<jupyter_output><empty_output><jupyter_text>**$(b)$ Using the validation set approach, estimate the test error of this model. In order to do this, you must perform the following steps:**- **$i.$ Split the sample set into a training set and a validation set.**<jupyter_code># Splitting the sample set
X_train, X_test, y_train, y_test = train_test_split(default_df[['income', 'balance']],
                                                    default_df.default_Yes,
                                                    test_size=0.5,
                                                    random_state=0)<jupyter_output><empty_output><jupyter_text>- **$ii.$ Fit a multiple logistic regression model using only the training observations.**<jupyter_code># Model
log_r = LogisticRegression().fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>- **$iii.$ Obtain a prediction of default status for each individual in the validation set by computing the posterior probability of default for that individual, and classifying the individual to the default category if the posterior probability is greater than 0.5.**<jupyter_code># Predict
y_pred = log_r.predict(X_test)<jupyter_output><empty_output><jupyter_text>- **$iv.$ Compute the validation set error, which is the fraction of the observations in the validation set that are misclassified.**<jupyter_code>val_err = 1 - sum(sum([y_test == y_pred]))/len(y_test)
print("Validation error rate =", round(val_err*100, 4),"%")<jupyter_output>Validation error rate = 2.9 %<jupyter_text>**$(c)$ Repeat the process in $(b)$ three times, using three different splits of the observations into a training set and a validation set. Comment on the results obtained.**<jupyter_code>for i in range(1,4):
    # Splitting the sample set
    X_train, X_test, y_train, y_test = train_test_split(default_df[['income', 'balance']],
                                                        default_df.default_Yes,
                                                        test_size=0.5,
                                                        random_state=i)
    
    # Logistic regresssion
    log_r = LogisticRegression().fit(X_train, y_train)
    
    # Predict
    y_pred = log_r.predict(X_test)

    val_err = 1 - sum(sum([y_test == y_pred]))/len(y_test)
    print(f"With random_state = {i}")
    print("Validation error rate =", round(val_err*100, 4),"%")
    
    # confusion matrix
    confusion_mtx = confusion_matrix(y_test, y_pred)
    plot_confusion_matrix(confusion_mtx)
    
    # overall stats
#    log_stats = conf_mtx_stats(confusion_matrix(y_test, y_pred))
#    display(log_stats)
    print("---------------------------------------")<jupyter_output>With random_state = 1
Validation error rate = 3.22 %<jupyter_text>>- **Comment:**>  - it seems to average around 3.16% test error rate.>  - As we face **unbalanced problem** in the **Default** dataset, the model tends to predict almost all as **Not Defaulted**. **$(d)$ Now consider a logistic regression model that predicts the probability of default using income, balance, and a dummy variable for student. Estimate the test error for this model using the validation set approach. Comment on whether or not including a dummy variable for student leads to a reduction in the test error rate.**<jupyter_code>for i in range(1,5):
    # Splitting the sample set
    X_train, X_test, y_train, y_test = train_test_split(default_df[['income', 'balance', 'student_Yes']],
                                                        default_df.default_Yes,
                                                        test_size=0.5,
                                                        random_state=i)
    
    # Logistic regresssion
    log_r = LogisticRegression().fit(X_train, y_train)
    
    # Predict
    y_pred = log_r.predict(X_test)

    val_err = 1 - sum(sum([y_test == y_pred]))/len(y_test)
    print(f"With random_state = {i}")
    print("Validation error rate =", round(val_err*100, 4),"%")
    
    # confusion matrix
    confusion_mtx = confusion_matrix(y_test, y_pred)
    plot_confusion_matrix(confusion_mtx)
    
    # overall stats
    log_stats = conf_mtx_stats(confusion_matrix(y_test, y_pred))
    display(log_stats)
    print("---------------------------------------")<jupyter_output>With random_state = 1
Validation error rate = 3.22 %<jupyter_text>>- **Comment:**>  - Now it seems to average around 3.07% test error rate, and that is a slightly reduction in the test error rate. >  - And we can see that including a dummy variable for student leads to an increasement in the Recall! which is the True Positive Rate. Thus, we can say including a dummy variable for student fixed a little in the **unbalanced problem**.  $6.$ We continue to consider the use of a logistic regression model to predict the probability of $\color{brown}{\text{default}}$ using $\color{brown}{\text{income}}$ and $\color{brown}{\text{balance}}$ on the $\color{brown}{\text{Default}}$ data set. In particular, we will now compute estimates for the standard errors of the $\color{brown}{\text{income}}$ and $\color{brown}{\text{balance}}$ logistic regression coefficients in two different ways: $(1)$ using the bootstrap, and $(2)$ using the standard formula for computing the standard errors in the $\color{brown}{\text{glm()}}$ function. Do not forget to set a random seed before beginning your analysis.**$(a)$ Using the $\color{brown}{\text{summary()}}$ and $\color{brown}{\text{glm()}}$ functions, determine the estimated standard errors for the coefficients associated with income and balance in a multiple logistic regression model that uses both predictors.**<jupyter_code># Model
f = 'default_Yes ~ income + balance'
log_r = smf.logit(formula=f, data=default_df).fit()

# Summary
display(log_r.summary())

statsmodels_est = pd.DataFrame({'coef_sm': log_r.params, 'SE_sm': log_r.bse})
display(statsmodels_est)<jupyter_output>Optimization terminated successfully.
         Current function value: 0.078948
         Iterations 10<jupyter_text>**$(b)$ Write a function, $\color{brown}{\text{boot.fn()}}$, that takes as input the $\color{brown}{\text{Default}}$ data set as well as an index of the observations, and that outputs the coefficient estimates for $\color{brown}{\text{income}}$ and $\color{brown}{\text{balance}}$ in the multiple logistic regression model.**<jupyter_code>def boot_fn(df):
    # Model
    f = 'default_Yes ~ income + balance'
    model_logit = smf.logit(formula=f, data=df).fit();
    return model_logit.params<jupyter_output><empty_output><jupyter_text>**$(c)$ Use the $\color{brown}{\text{boot()}}$ function together with your $\color{brown}{\text{boot.fn() }}$function to estimate the standard errors of the logistic regression coefficients for $\color{brown}{\text{income}}$ and $\color{brown}{\text{balance}}$.**<jupyter_code>coef_preds = []

for i in range(2000):
    # bootstrap sample
    sample = default_df.sample(len(default_df), replace=True)
    coef_preds.append(boot_fn(sample))
coef_pred     = np.mean(coef_preds, axis=0)
standard_errs = np.std(coef_preds, axis=0)

bootstrap_est = pd.DataFrame({'coef_boot': coef_pred,
                              'SE_boot': standard_errs})
display(bootstrap_est)<jupyter_output><empty_output><jupyter_text>**$(d)$ Comment on the estimated standard errors obtained using the $\color{brown}{\text{glm()}}$ function and using your bootstrap function.**<jupyter_code>pd.concat([statsmodels_est, bootstrap_est.set_index(statsmodels_est.index)], axis=1)<jupyter_output><empty_output><jupyter_text>>- We can see that the standard error estimates for **income** are equal to 6 decimal places. and the estimates for **balance** are equal to 4 decimal places. This could be an indication that the assumptions of the logistic $SE$ estimates are well-satisfied.>- We can see that the standard errors estimated by the statsmodels approximately equal to the estimates obtained by bootstrap, and with just 10 bootstrap samples we got a slightly more significant than with 2,000 bootstrap samples!  $7.$ In Sections $5.3.2$ and $5.3.3$, we saw that the $\color{brown}{\text{cv.glm()}}$ function can be used in order to compute the LOOCV test error estimate. Alternatively, one could compute those quantities using just the $\color{brown}{\text{glm()}}$ and $\color{brown}{\text{predict.glm()}}$ functions, and a for loop. You will now take this approach in order to compute the LOOCV error for a simple logistic regression model on the $\color{brown}{\text{Weekly}}$ data set. Recall that in the context of classification problems, the LOOCV error is given in $(5.4)$.<jupyter_code># Read the Weekly dataset
weekly_df = pd.read_csv('datasets/Weekly.csv')

# Get all dummies in the df
weekly_df = pd.get_dummies(weekly_df, drop_first=True)
weekly_df.head()<jupyter_output><empty_output><jupyter_text>**$(a)$ Fit a logistic regression model that predicts $\color{brown}{\text{Direction}}$ using $\color{brown}{\text{Lag1}}$ and $\color{brown}{\text{Lag2}}$.**<jupyter_code># Splitting the sample set
X_train, X_test, y_train, y_test = train_test_split(weekly_df[['Lag1', 'Lag2']],
                                                    weekly_df.Direction_Up,
                                                    test_size=0.5,
                                                    random_state=42)

# Logistic regresssion
log_r = LogisticRegression().fit(X_train, y_train)

# Predict
y_pred = log_r.predict(X_test)

test_err = 1 - sum(sum([y_test == y_pred]))/len(y_test)
print("test error rate =", round(test_err*100, 4),"%")

# confusion matrix
confusion_mtx = confusion_matrix(y_test, y_pred)
plot_confusion_matrix(confusion_mtx)

# overall stats
log_stats = conf_mtx_stats(confusion_matrix(y_test, y_pred))
display(log_stats)<jupyter_output>test error rate = 43.6697 %<jupyter_text>**$(b)$ Fit a logistic regression model that predicts $\color{brown}{\text{Direction}}$ using $\color{brown}{\text{Lag1}}$ and $\color{brown}{\text{Lag2}}$ using all but the first observation.**<jupyter_code># Model
X = weekly_df[['Lag1', 'Lag2']].iloc[1:]
y = weekly_df.Direction_Up.iloc[1:]
log_r = LogisticRegression().fit(X, y)<jupyter_output><empty_output><jupyter_text>**$(c)$ Use the model from $(b)$ to predict the direction of the first observation. You can do this by predicting that the first observation will go up if $P\color{brown}{\text{(Direction="Up"|Lag1, Lag2)}} > 0.5$. Was this observation correctly classified?**<jupyter_code>y = [weekly_df.Direction_Up.iloc[0]]
y_pred = log_r.predict([weekly_df[['Lag1', 'Lag2']].iloc[0]])

test_err = 1 - sum(sum([y == y_pred]))/len(y)
print("test error rate =", round(test_err*100, 4),"%")

# confusion matrix
confusion_mtx = confusion_matrix(y, y_pred)
plot_confusion_matrix(confusion_mtx)<jupyter_output>test error rate = 100.0 %<jupyter_text>**$(d)$ Write a for loop from $i = 1$ to $i = n$, where $n$ is the number of observations in the data set, that performs each of the following steps:** - **$i.$ Fit a logistic regression model using all but the $i$th observation to predict $\color{brown}{\text{Direction}}$ using $\color{brown}{\text{Lag1}}$ and $\color{brown}{\text{Lag2}}$.** - **$ii.$ Compute the posterior probability of the market moving up for the ith observation.** - **$iii.$ Use the posterior probability for the $i$th observation in order to predict whether or not the market moves up.** - **$iv.$ Determine whether or not an error was made in predicting the direction for the $i$th observation. If an error was made,then indicate this as a 1, and otherwise indicate it as a 0.**<jupyter_code>y_pred = []

for i in range(len(weekly_df)):
    
    idx = weekly_df.index[weekly_df.index != i]
    
    X_train = weekly_df[['Lag1', 'Lag2']].iloc[idx, :]
    X_test  = [weekly_df[['Lag1', 'Lag2']].iloc[i, :]]
    
    y_train = weekly_df.Direction_Up.iloc[idx]
    y_test  = [weekly_df.Direction_Up.iloc[i]]
    
    # Model
    log_r = LogisticRegression().fit(X_train, y_train)
    
    # Predict
    y_pred.append(log_r.predict(X_test)[0])<jupyter_output><empty_output><jupyter_text>**$(e)$ Take the average of the $n$ numbers obtained in $(d)iv$ in order to obtain the LOOCV estimate for the test error. Comment on the results.**<jupyter_code>X_test
y_test = np.array(weekly_df.Direction_Up)

test_err = 1 - sum(y_test == y_pred)/len(y_pred)
print("test error rate =", round(test_err*100, 4),"%")

# confusion matrix
confusion_mtx = confusion_matrix(y_test, np.array(y_pred))
plot_confusion_matrix(confusion_mtx)<jupyter_output>test error rate = 44.9954 %<jupyter_text>>- **Comments:**>  - `LOOCV` produces a test error rate of $45.0\%$, higher than the $43.67\%$ error rate observed with a $50\%$ holdout set.>  - Still our classifier is not that much good.>  - However, The `LOOCV` approach allows the model to train on a $n - 1$ observations larger than 50% holdout which train on $n/2$, which could be significant if we were to experiment with more flexible models that require more observations to mitigate over-fitting.   $8.$ We will now perform cross-validation on a simulated data set.**$(a)$ Generate a simulated data set as follows:**```> set .seed (1)> y=rnorm (100)> x=rnorm (100)> y=x-2* x^2+ rnorm (100)```**In this data set, what is $n$ and what is $p$? Write out the model used to generate the data in equation form.**<jupyter_code>np.random.seed(1)
x = np.random.normal(size=100)
y = x - (2*x**2) + np.random.normal(size=100)<jupyter_output><empty_output><jupyter_text>>- $n = 100$>- $p = 2$>- $Y = X - 2 X^2 + \epsilon$. **$(b)$ Create a scatterplot of $X$ against $Y$ . Comment on what you find.**<jupyter_code>ax = sns.scatterplot(x=x, y=y)
plt.xlabel('x')
plt.ylabel('y');<jupyter_output><empty_output><jupyter_text>>- The above graph shows $X$ plotted against $Y$.>- From the above graph, it is clear that there is a linear relationship between $X$ and $Y$ from the shape of the curve with some variance.>- The scope of the graph has 2 maximum and -2 minimum with a quadratic shape. **$(c)$ Set a random seed, and then compute the LOOCV errors that result from fitting the following four models using least squares:**- $i.$   $Y = \beta_0 + \beta_1 X + \epsilon$  - $ii.$ $Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \epsilon$  - $iii.$ $Y = \beta_0 + β_1 X + \beta_2 X^2 + \beta_3 X^3 + \epsilon$  - $iv.$ $Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 X^3 + \beta_4 X^4 + \epsilon$  **Note you may find it helpful to use the $\color{brown}{\text{data.frame()}}$ function to create a single data set containing both $X$ and $Y$ .**<jupyter_code>def mse(y_true, y_pred):
    """Calculate mean squared error"""
    return np.sum(np.square(y_pred - y_true)) / y_true.size

np.random.seed(1)
x = np.random.normal(size=100)
y = x - (2*x**2) + np.random.normal(size=100)
df = pd.DataFrame({'x': x, 'y': y})

formulae = {'x'   : 'y ~ x', 
            'x^2' : 'y ~ x + np.power(x, 2)',
            'x^3' : 'y ~ x + np.power(x, 2) + np.power(x, 3)',
            'x^4' : 'y ~ x + np.power(x, 2) + np.power(x, 3) + np.power(x, 4)'}

errors = {}
for f in formulae:
    # predictions state
    y_pred = []
    
    for i in range(df.shape[0]):

        # Model
        train = df.drop(i)
        l_r = smf.ols(formula=formulae[f], data=df.drop(i)).fit()

        ## Predict
        y_pred  += [float(l_r.predict(exog=dict(x=df.iloc[i]['x'])))]
    errors[f] = mse(y, np.array(y_pred))
errors<jupyter_output><empty_output><jupyter_text>**$(d)$ Repeat $(c)$ using another random seed, and report your results.Are your results the same as what you got in $(c)$? Why?**<jupyter_code>np.random.seed(2)
x = np.random.normal(size=100)
y = x - (2*x**2) + np.random.normal(size=100)
df = pd.DataFrame({'x': x, 'y': y})

formulae = {'x'   : 'y ~ x', 
            'x^2' : 'y ~ x + np.power(x, 2)',
            'x^3' : 'y ~ x + np.power(x, 2) + np.power(x, 3)',
            'x^4' : 'y ~ x + np.power(x, 2) + np.power(x, 3) + np.power(x, 4)'}

errors = {}
for f in formulae:
    # predictions state
    y_pred = []
    
    for i in range(df.shape[0]):

        # Model
        train = df.drop(i)
        l_r = smf.ols(formula=formulae[f], data=df.drop(i)).fit()

        ## Predict
        y_pred  += [float(l_r.predict(exog=dict(x=df.iloc[i]['x'])))]
    errors[f] = mse(y, np.array(y_pred))
errors<jupyter_output><empty_output><jupyter_text>>- Changing the random seed that is used to simulate the observations has a large effect on the observed mean squared error.>- In changing the random seed we have changed the sample of observations taken from the population.>- This change in error is due to the small sample size with $n=100$, and so high variance in our data. **$(e)$ Which of the models in $(c)$ had the smallest LOOCV error? Is this what you expected? Explain your answer.** >- The $iv.$ model has the smallest LOOCV error.>- Although the second, third and fourth order polynomial fits perform almost identically.>- This is as we would expect, as the more flexible models "higher-order polynomials" tend to overfit the training set, and as the relationship between $X$ and $Y$ in $(a)$ already quadratic. **$(f)$ Comment on the statistical significance of the coefficient estimates that results from fitting each of the models in $(c)$ using least squares. Do these results agree with the conclusions drawn based on the cross-validation results?**<jupyter_code>np.random.seed(1)
x = np.random.normal(size=100)
y = x - (2*x**2) + np.random.normal(size=100)
df = pd.DataFrame({'x': x, 'y': y})

formula = {'x'    : 'y ~ x', 
            'x^2' : 'y ~ x + np.power(x, 2)',
            'x^3' : 'y ~ x + np.power(x, 2) + np.power(x, 3)',
            'x^4' : 'y ~ x + np.power(x, 2) + np.power(x, 3) + np.power(x, 4)'}

p_values = {}
for f in formulae:
    p_val = [0]
    for i in range(df.shape[0]):

        # Model
        train = df.drop(i)
        l_r = smf.ols(formula=formula[f], data=df.drop(i)).fit()

        ## p-values
        try:
            p_val += l_r.pvalues
        except:
            p_val = l_r.pvalues

    p_values[f] = round(p_val/100, 5)
for k, v in p_values.items():
    print(k)
    print(v)
    print("-------------------------")<jupyter_output>x
Intercept    0.00000
x            0.00277
dtype: float64
-------------------------
x^2
Intercept         0.24966
x                 0.00000
np.power(x, 2)    0.00000
dtype: float64
-------------------------
x^3
Intercept         0.27376
x                 0.00001
np.power(x, 2)    0.00000
np.power(x, 3)    0.23250
dtype: float64
-------------------------
x^4
Intercept         0.02446
x                 0.00000
np.power(x, 2)    0.00000
np.power(x, 3)    0.12781
np.power(x, 4)    0.01814
dtype: float64
-------------------------<jupyter_text>>- The $\text{p-value}$ in $i.$ with $x$ as the only predictor, seems to be statistical significance with a $\text{p-value}$ under our threshold of $0.05$.>- In $ii.$ with the quadratic term, still statistical significance.>- In $iii.$ the quadratic term still statistical significance, but still there's no evidence that the cubic term is not a statistical significance, which concluded that there's no reason to add the cubic term in the model.>- in $iv.$ the quadratic and the Quartic terms are statistical significance, but still there's no evidence that the cubic term is statistical significance.  $9.$ We will now consider the $\color{brown}{\text{Boston}}$ housing data set, from the $\color{brown}{\text{MASS}}$ library.<jupyter_code># read boston dataset 
boston_df = pd.read_csv('datasets/Boston.csv').drop('Unnamed: 0', axis=1)
boston_df.head()<jupyter_output><empty_output><jupyter_text>**$(a)$ Based on this data set, provide an estimate for the population mean of $\color{brown}{\text{medv}}$. Call this estimate $\hat{\mu}$.**<jupyter_code>mu_hat = boston_df.medv.mean()
mu_hat<jupyter_output><empty_output><jupyter_text>Here! **$(b)$ Provide an estimate of the standard error of $\hat{\mu}$. Interpret this result.*****Hint: We can compute the standard error of the sample mean by dividing the sample standard deviation by the square root of the number of observations.***<jupyter_code>se = np.std(boston_df.medv) / np.sqrt(len(boston_df))
se<jupyter_output><empty_output><jupyter_text>**$(c)$ Now estimate the standard error of $\hat{\mu}$ using the bootstrap. How does this compare to your answer from $(b)$?**<jupyter_code>st_err = []
n = 2000

for i in range(n):
    # bootstrap sample
    sample = boston_df.medv.sample(len(boston_df), replace=True)
    
    # standard error of u^
    s_err = np.std(sample) / np.sqrt(len(boston_df))
    st_err.append(s_err)

boot_se = np.mean(st_err)
boot_se<jupyter_output><empty_output><jupyter_text>>- The bootstrap method gives a approximately similar $SE$ as the answer from $(b)$, and that's a remarkably good estimate of the standard error, when compared to the same estimate derived analytically in $(b)$.>- The bootstrap approach is computationally more expensive, but has the advantage that no analytic derivation of the standard error for the statistic is required, making it much more general for application to other statistics. **$(d)$ Based on your bootstrap estimate from $(c)$, provide a $95\%$ confidence interval for the mean of $\color{brown}{\text{medv}}$. Compare it to the results obtained using $\color{brown}{\text{t.test(Boston\$medv)}}$. Hint: You can approximate a $95\%$ confidence interval using the formula $[\hat{\mu} - 2SE(\hat{\mu}), \hat{\mu} + 2SE(\hat{\mu})]$.**<jupyter_code>mu_hat = round(np.mean(boston_df.medv), 6)
conf_l = round(mu_hat - (2*boot_se), 6)
conf_h = round(mu_hat + (2*boot_se), 6)

print('mu^:', mu_hat)
print('SE: ', round(boot_se, 6))
print(f'95% confidence interval: [{conf_l}, {conf_h}]')<jupyter_output>mu^: 22.532806
SE:  0.407535
95% confidence interval: [21.717736, 23.347876]<jupyter_text>**$(e)$ Based on this data set, provide an estimate, ${\hat{\mu}_{med}}$, for the median value of $\color{brown}{\text{medv}}$ in the population.**<jupyter_code>median_hat = np.median(boston_df['medv'])
median_hat<jupyter_output><empty_output><jupyter_text>**$(f)$ We now would like to estimate the standard error of ${\hat{\mu}_{med}}$. Unfortunately, there is no simple formula for computing the standard error of the median. Instead, estimate the standard error of the median using the bootstrap. Comment on your findings.**<jupyter_code>medians = []
n = 2000

for i in range(n):
    # bootstrap sample
    sample = boston_df.medv.sample(len(boston_df), replace=True)
    
    # the sample median
    median = np.median(sample)
    medians.append(median)

med_se = np.std(medians)
med_se<jupyter_output><empty_output><jupyter_text>>- As with the mean, the standard error is quite small relative to the median value. **$(g)$ Based on this data set, provide an estimate for the tenth percentile of $\color{brown}{\text{medv}}$ in Boston suburbs. Call this quantity $\hat{\mu}_{0.1}$. (You can use the $\color{brown}{\text{quantile()}}$ function.)**<jupyter_code>ten_perc = np.percentile(boston_df.medv, 10)
ten_perc<jupyter_output><empty_output><jupyter_text>**$(h)$ Use the bootstrap to estimate the standard error of $\hat{\mu}_{0.1}$. Comment on your findings.**<jupyter_code>ten_percs = []
n = 2000

for i in range(n):
    # bootstrap sample
    sample = boston_df.medv.sample(len(boston_df), replace=True)
    
    # tenth percentile
    ten_perc = np.percentile(sample, 10)
    ten_percs.append(ten_perc)

ten_perc_se = np.std(ten_percs)
ten_perc_se<jupyter_output><empty_output>