Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Classification - Exercise R code as soutution manual ISLR Introduction to Statistical Learning James, Witten, Hastie, Tibshirani
Typology: Exercises
1 / 19
This page cannot be seen from the preview
Don't miss anything!
title: "Chapter 4: Classification" author: "Solutions to Exercises" date: "January 12, 2016" output: html_document: keep_md: no
>EXERCISE 1: Let $Z=e^{\beta_0+\beta_1X}$, Equation (4.2) becomes Step 1: $p(X) = \frac{Z}{1+Z}$ Step 2: $\frac{1}{p(X)} = \frac{1+Z}{Z} = 1+\frac{1}{Z}$ Step 3: $Z = \frac{1}{\frac{1}{p(X)}-1} = \frac{1}{\frac{1-p(X)}{p(X)}} = \frac{p(X)}{1-p(X)}$
Equation (4.12): $p_k(x) = \frac {\pi_k \frac {1} {\sqrt{2 \pi} \sigma} \exp(- \frac {1} {2 \sigma^2} (x - mu_k)^2) } {\sum { \pi_l \frac {1} {\sqrt{2 \pi} \sigma} \exp(- \frac {1} {2 \sigma^2} (x - \mu_l)^2) }}$ Substitute $C = \frac { \frac {1} {\sqrt{2 \pi} \sigma} \exp(- \frac {1} {2 \sigma^2} (x^2)) } {\sum { \pi_l frac {1} {\sqrt{2 \pi} \sigma} \exp(- \frac {1} {2 \sigma^2} (x - \mu_l)^2) }}$ as this term does not vary across $k$ Step 1: Equation becomes $p_k(x) = C \pi_k \exp(- \frac {1} {2 \sigma^2} (\mu_k^2 - 2x \mu_k))$ Step 2: Take log of both sides $log(p_k(x)) = log(C) + log(\pi_k) + (- \frac {1} {2 \sigma^2} (\mu_k^
>EXERCISE 3: If $\sigma$ varies by $k$ then Equation (4.12) becomes: $p_k(x) = \frac {\pi_k \frac {1} {\sqrt{2 \pi} sigma_k} \exp(- \frac {1} {2 \sigma_k^2} (x - \mu_k)^2) } {\sum { \pi_l \frac {1} {\sqrt{2 \pi} \sigma_k} exp(- \frac {1} {2 \sigma_k^2} (x - \mu_l)^2) }}$ The constant term that does not vary by $k$ becomes $C' = \frac { \frac {1} {\sqrt{2 \pi}}} {\sum { \pi_l frac {1} {\sqrt{2 \pi} \sigma_k} \exp(- \frac {1} {2 \sigma_k^2} (x - \mu_l)^2) }}$ Step 1: Equation becomes $p_k(x) = C' \frac{\pi_k}{\sigma_k} \exp(- \frac {1} {2 \sigma_k^2} (x - mu_k)^2)$
Part e)
>EXERCISE 5: Part a) If the actual decision boundary is linear, then we would expect LDA to perform better on the test set. For the training set, QDA has a chance of performing better if it overfits. Part b) QDA would likely perform better on both the training set and the test set. Part c) In general a large sample size is more beneficial for QDA so would expect QDA accuracy to increase more than LDA.
Part d) FALSE: We might achieve a better error rate on the training set but not on the test set because if the true decision boundary is linear then the QDA is not flexible in any predictive way.
>EXERCISE 6: Part a) For logistic regression, $p(X) = \frac{e^{\beta_0+\beta_1 X_1+\beta_2 X_2}}{1+e^{\beta_0+\beta_ X_1+\beta_2 X_2}}$ Plugging in the values $p(X) = \frac{e^{-6 + 0.05 \times 40 + 1 \times 3.5}}{1+e^{-6+0.05 \times 40 + 1 times 3.5}} =$
exp(-6+0.05*40+1*3.5)/(1+exp(-6+0.05*40+1*3.5)) #0.Part b) Solve this equation $0.5 = \frac{e^{-6 + 0.05 X_1 + 1 \times 3.5}}{1+e^{-6+0.05 X_1 + 1 \times 3.5}}$ Which equates to solving the logit equation $log(\frac{0.5}{1-0.5}) = -6 + 0.05 X_1 + 1 \times 3.5$
(log(0.5/(1-0.5)) + 6 - 3.5*1)/0.05 # are selecting the model with only error rate data, then we want to know which model has the lower __test__ error rate. *** >EXERCISE 9: __Part a)__ We want to solve $0.37 = \frac{p_{default}}{1-p_{default}}$ Rearranging, this becomes $\frac{1}{0.37} = \frac{1-p_{default}}{p_{default}} = \frac{1}{p_{default}}-1$ Finally $p_{default} = \frac{1}{\frac{1}{0.37}+1}$ ```{r} 1/(1/0.37+1)Probability of default is 27.0% Part b)
0.16/(1-0.16)Odds of defaulting is 0.
>EXERCISE 10: Part a)
require(ISLR) data(Weekly) summary(Weekly) pairs(Weekly)Year and Volume are positively correlated similar to the Smarket data set. Part b)
fit.logit <- glm(Direction~., data=Weekly[,c(2:7,9)], family=binomial) summary(fit.logit)Lag2 seems to have statistically significant predictive value Part c)
require(MASS) fit.lda <- lda(Direction~Lag2, data=train) fit.lda.pred <- predict(fit.lda, test)$class table(fit.lda.pred, test$Direction) mean(fit.lda.pred == test$Direction) # Accuracy=0.
__Part f)__ ```{r, warning=FALSE, message=FALSE} fit.qda <- qda(Direction~Lag2, data=train) fit.qda.pred <- predict(fit.qda, test)$class table(fit.qda.pred, test$Direction) mean(fit.qda.pred == test$Direction) # Accuracy=0.Part g)
require(class) set.seed(1) train.X <- as.matrix(train$Lag2) test.X <- as.matrix(test$Lag2) knn.pred <- knn(train.X, test.X, train$Direction, k=1) table(knn.pred, test$Direction) mean(knn.pred == test$Direction) # Accuracy=0.Part h)
The Logistic Regression and LDA models produced the best results Part i)
knn.pred <- knn(train.X, test.X, train$Direction, k=5) table(knn.pred, test$Direction) mean(knn.pred == test$Direction) knn.pred <- knn(train.X, test.X, train$Direction, k=10) table(knn.pred, test$Direction) mean(knn.pred == test$Direction) knn.pred <- knn(train.X, test.X, train$Direction, k=20) table(knn.pred, test$Direction) mean(knn.pred == test$Direction) knn.pred <- knn(train.X, test.X, train$Direction, k=30) table(knn.pred, test$Direction) mean(knn.pred == test$Direction)Higher k values for KNN (around 20) seemed to produce the best results when using only Lag2 as predictor.
fit.lda <- lda(Direction~Lag2+I(Lag1^2), data=train) fit.lda.pred <- predict(fit.lda, test)$class table(fit.lda.pred, test$Direction) mean(fit.lda.pred == test$Direction) # Accuracy=0.test <- mydf[-trainid,]
__Part d)__ ```{r, warning=FALSE, message=FALSE} fit.lda <- lda(mpg01~displacement+horsepower+weight+acceleration, data=train) fit.lda.pred <- predict(fit.lda, test)$class table(fit.lda.pred, test$mpg01) mean(fit.lda.pred != test$mpg01) # error ratePart e)
fit.qda <- qda(mpg01~displacement+horsepower+weight+acceleration, data=train) fit.qda.pred <- predict(fit.qda, test)$class table(fit.qda.pred, test$mpg01) mean(fit.qda.pred != test$mpg01) # error ratePart f)
fit.logit <- glm(mpg01~displacement+horsepower+weight+acceleration, data=train, family=binomial) logit.prob <- predict(fit.logit, test, type="response") logit.pred <- ifelse(logit.prob > 0.5, 1, 0) table(logit.pred, test$mpg01) mean(logit.pred != test$mpg01) # error rate ## ``` __Part g)__ ```{r, warning=FALSE, message=FALSE} train.X <- cbind(train$displacement, train$horsepower, train$weight, train$acceleration) test.X <- cbind(test$displacement, test$horsepower, test$weight, test$acceleration) knn.pred <- knn(train.X, test.X, train$mpg01, k=1) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01) knn.pred <- knn(train.X, test.X, train$mpg01, k=10) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01) knn.pred <- knn(train.X, test.X, train$mpg01, k=20) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01) knn.pred <- knn(train.X, test.X, train$mpg01, k=30) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01) knn.pred <- knn(train.X, test.X, train$mpg01, k=50) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01) knn.pred <- knn(train.X, test.X, train$mpg01, k=100) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01) knn.pred <- knn(train.X, test.X, train$mpg01, k=200) table(knn.pred, test$mpg01) mean(knn.pred != test$mpg01)Power2(10,3) Power2(8,17) Power2(131,3)
__Part d)__ ```{r} Power3 <- function(x, a) { return(x^a) } Power3(3,8)Part e)
x <- 1: plot(x, Power3(x,2), log="y", main="log(x^2) vs. x", xlab="x", ylab="log(x^2)")Part f)
PlotPower <- function(x, a) { plot(x, Power3(x,2), main="x^a versus x", xlab="x", ylab=paste0("x^",a)) } PlotPower(1:10,3)>EXERCISE 13:
data(Boston) summary(Boston) crim01 <- ifelse(Boston$crim > median(Boston$crim), 1, 0) mydf <- data.frame(Boston, crim01) pairs(mydf) # pred1 = age, dis, lstat, medv sort(cor(mydf)[1,]) # pred2 = tax, rad (highest correlations with crim) set.seed(1) trainid <- sample(1:nrow(mydf), nrow(mydf)*0.7 , replace=F) # 70% train, 30% test train <- mydf[trainid,] test <- mydf[-trainid,] train.X1 <- cbind(train$age, train$dis, train$lstat, train$medv) test.X1 <- cbind(test$age, test$dis, test$lstat, test$medv) train.X2 <- cbind(train$tax, train$rad) test.X2 <- cbind(test$tax, test$rad) # Logistic Regression models fit.logit1 <- glm(crim01~age+dis+lstat+medv, data=train, family=binomial) logit1.prob <- predict(fit.logit1, test, type="response") logit1.pred <- ifelse(logit1.prob > 0.5, 1, 0) mean(logit1.pred != test$crim01) # error rate fit.logit2 <- glm(crim01~tax+rad, data=train, family=binomial) knn1.pred <- knn(train.X1, test.X1, train$crim01, k=100) mean(knn1.pred != test$crim01) knn1.pred <- knn(train.X1, test.X1, train$crim01, k=200) mean(knn1.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=1) mean(knn2.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=5) mean(knn2.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=10) mean(knn2.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=20) mean(knn2.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=50) mean(knn2.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=100) mean(knn2.pred != test$crim01) knn2.pred <- knn(train.X2, test.X2, train$crim01, k=200) mean(knn2.pred != test$crim01)Surprisingly, the KNN model with two predictors tax and rad and k=1 had the best error rate
