ch04soln.Rmd

---
title: "Chapter 4: Classification"
author: "Solutions to Exercises"
date: "January 12, 2016"
output: 
  html_document: 
    keep_md: no
---

***
## CONCEPTUAL
***

<a id="ex01"></a>

>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)}$

***

<a id="ex02"></a>

>EXERCISE 2:

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^2 - 2x \mu_k))$
  
  __Step 3:__ Simplify and rearrange $log(p_k(x)) =  (\frac {2x \mu_k} {2 \sigma^2} -\frac {\mu_k^2} {2 \sigma^2}) + log(\pi_k) + log(C)$

***

<a id="ex03"></a>

>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)$
  
  __Step 2:__ Take log of both sides $log(p_k(x)) = log(C') + log(\pi_k) - log(\sigma_k) + (- \frac {1} {2 \sigma_k^2} (x - \mu_k)^2)$
  
  __Step 3:__ Simplify and rearrange $log(p_k(x)) = (- \frac {1} {2 \sigma_k^2} (x^2 + \mu_k^2 - 2x\mu_k)) + log(\pi_k) - log(\sigma_k) + log(C')$
  
There's the $x^2$.

***

<a id="ex04"></a>

>EXERCISE 4:

__Part a)__

If $X$ is uniformly distributed, then (0.65-0.55)/(1-0) = 10%

__Part b)__

For two features, $10\% \times 10\% = 1\%$

__Part c)__

For 100 features, $10\%^{100}=$ _a very small number_

__Part d)__

When there are a large number of dimensions, the percentage of observations that can be used to predict with KNN becomes very small. This means that for a set sample size, more features leads to fewer neighbors.

__Part e)__

* For p=1, side = 0.1
* For p=2, side = 0.1^(1/2) = 0.316
* For p=100, side = 0.1^(1/100) = 0.977

This is saying that when the number of features is high (i.e. p=100), to use on average 10% of the training observations would mean that we would need to include almost the entire range of each individual feature.

***

<a id="ex05"></a>

>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.

***

<a id="ex06"></a>

>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_1 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}} =$

```{r}
exp(-6+0.05*40+1*3.5)/(1+exp(-6+0.05*40+1*3.5))  #0.38
```

__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$

```{r}
(log(0.5/(1-0.5)) + 6 - 3.5*1)/0.05  #50
```

Student needs to study for 50 hours.

***

<a id="ex07"></a>

>EXERCISE 7:

For constant variance, $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) }}$

Evaluating this becomes $p_{yes}(4) = \frac {0.8 \exp(- \frac {1} {2 \times 36} (4 - 10)^2)} {0.8 \exp(- \frac {1} {2 \times 36} (4 - 10)^2) + (1-0.8) \exp(- \frac {1} {2 \times 36} (4 - 0)^2)}$

```{r}
(0.8*exp(-1/(2*36)*(4-10)^2))/(0.8*exp(-1/(2*36)*(4-10)^2)+(1-0.8)*exp(-1/(2*36)*(4-0)^2))
```

Probability is 75.2%

***

<a id="ex08"></a>

>EXERCISE 8:

There's not enough information to say which method is better. With such a high error rate for the logistic regression, it's possible that the true decision boundary is not linear, so KNN=1 might have a better fit. On the other hand, KNN=1 has a high propensity to overfit. With KNN=1 having an average error of 18%, it's possible that the training error is close to 0% and the test error is more than 30%. If we are selecting the model with only error rate data, then we want to know which model has the lower __test__ error rate.

***

<a id="ex09"></a>

>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)__

```{r}
0.16/(1-0.16)
```

Odds of defaulting is 0.19

***
## APPLIED
***

<a id="ex10"></a>

>EXERCISE 10:

__Part a)__

```{r, warning=FALSE, message=FALSE}
require(ISLR)
data(Weekly)
summary(Weekly)
pairs(Weekly)
```

`Year` and `Volume` are positively correlated similar to the `Smarket` data set.

__Part b)__

```{r, warning=FALSE, message=FALSE}
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)__

```{r, warning=FALSE, message=FALSE}
logit.prob <- predict(fit.logit, Weekly, type="response")
logit.pred <- ifelse(logit.prob > 0.5, "Up", "Down")
table(logit.pred, Weekly$Direction)
(54+557)/nrow(Weekly)  # Accuracy=0.56
```

* When prediction is "Down", model is right 54/(54+48)=52.9%.
* When prediction is "Up", model is right 557/(430+557)=56.4%

Model is has higher accuracy when the prediction is "Up"

__Part d)__

```{r, warning=FALSE, message=FALSE}
train.yrs <- Weekly$Year %in% (1990:2008)
train <- Weekly[train.yrs,]
test <- Weekly[!train.yrs,]
fit2 <- glm(Direction~Lag2, data=train, family=binomial)
fit2.prob <- predict(fit2, test, type="response")
fit2.pred <- ifelse(fit2.prob > 0.5, "Up", "Down")
table(fit2.pred, test$Direction)
mean(fit2.pred == test$Direction)  # Accuracy=0.625
```

__Part e)__

```{r, warning=FALSE, message=FALSE}
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.625
```

__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.587
```

__Part g)__

```{r, warning=FALSE, message=FALSE}
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.500
```

__Part h)__

The Logistic Regression and LDA models produced the best results

__Part i)__

```{r, warning=FALSE, message=FALSE}
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.

```{r, warning=FALSE, message=FALSE}
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.644
```

***

<a id="ex11"></a>

>EXERCISE 11:

__Part a)__

```{r, warning=FALSE, message=FALSE}
require(ISLR)
data(Auto)
mpg01 <- ifelse(Auto$mpg > median(Auto$mpg), 1, 0)
mydf <- data.frame(Auto, mpg01)
```

__Part b)__

```{r, warning=FALSE, message=FALSE}
pairs(mydf)
```

`displacement`, `horsepower`, `weight` and `acceleration` seem to be highly correlated

__Part c)__

```{r, warning=FALSE, message=FALSE}
set.seed(1)
trainid <- sample(1:nrow(mydf), nrow(mydf)*0.7 , replace=F)  # 70% train, 30% test
train <- mydf[trainid,]
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 rate
```

__Part e)__

```{r, warning=FALSE, message=FALSE}
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 rate
```

__Part f)__

```{r, warning=FALSE, message=FALSE}
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)
```

KNN performs best around k=30 and k=100

***

<a id="ex12"></a>

>EXERCISE 12:

__Part a)__

```{r}
Power <- function() {
  print(2^3)
}
Power()
```

__Part b)__

```{r}
Power2 <- function(x, a) {
  print(x^a)
}
Power2(3,8)
```

__Part c)__

```{r}
Power2(10,3)
Power2(8,17)
Power2(131,3)
```

__Part d)__

```{r}
Power3 <- function(x, a) {
  return(x^a)
}
Power3(3,8)
```

__Part e)__

```{r}
x <- 1:10
plot(x, Power3(x,2), log="y", main="log(x^2) vs. x",
     xlab="x", ylab="log(x^2)")
```

__Part f)__

```{r}
PlotPower <- function(x, a) {
  plot(x, Power3(x,2), main="x^a versus x",
       xlab="x", ylab=paste0("x^",a))
}
PlotPower(1:10,3)
```

***

<a id="ex13"></a>

>EXERCISE 13:

```{r, warning=FALSE, message=FALSE}
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)
logit2.prob <- predict(fit.logit2, test, type="response")
logit2.pred <- ifelse(logit2.prob > 0.5, 1, 0)
mean(logit2.pred != test$crim01)  # error rate

# LDA models
fit.lda1 <- lda(crim01~age+dis+lstat+medv, data=train)
fit.lda1.pred <- predict(fit.lda1, test)$class
mean(fit.lda1.pred != test$crim01)  # error rate
fit.lda2 <- lda(crim01~tax+rad, data=train)
fit.lda2.pred <- predict(fit.lda2, test)$class
mean(fit.lda2.pred != test$crim01)  # error rate

# QDA models
fit.qda1 <- qda(crim01~age+dis+lstat+medv, data=train)
fit.qda1.pred <- predict(fit.qda1, test)$class
mean(fit.qda1.pred != test$crim01)  # error rate
fit.qda2 <- qda(crim01~tax+rad, data=train)
fit.qda2.pred <- predict(fit.qda2, test)$class
mean(fit.qda2.pred != test$crim01)  # error rate

# KNN models
set.seed(1)
knn1.pred <- knn(train.X1, test.X1, train$crim01, k=1)
mean(knn1.pred != test$crim01)
knn1.pred <- knn(train.X1, test.X1, train$crim01, k=5)
mean(knn1.pred != test$crim01)
knn1.pred <- knn(train.X1, test.X1, train$crim01, k=10)
mean(knn1.pred != test$crim01)
knn1.pred <- knn(train.X1, test.X1, train$crim01, k=20)
mean(knn1.pred != test$crim01)
knn1.pred <- knn(train.X1, test.X1, train$crim01, k=50)
mean(knn1.pred != test$crim01)
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