10-deep-learning.Rmd

# Deep Learning

## Conceptual

### Question 1

> Consider a neural network with two hidden layers: $p = 4$ input units, 2 units
> in the first hidden layer, 3 units in the second hidden layer, and a single
> output.
>
> a. Draw a picture of the network, similar to Figures 10.1 or 10.4.

```{r, echo=FALSE, out.width="80%"}
knitr::include_graphics("images/nn.png")
```

> b. Write out an expression for $f(X)$, assuming ReLU activation functions. Be
> as explicit as you can!

The three layers (from our final output layer back to the start of our network)
can be described as:

\begin{align*}
f(X) &= g(w_{0}^{(3)} + \sum^{K_2}_{l=1} w_{l}^{(3)} A_l^{(2)}) \\
A_l^{(2)} &= h_l^{(2)}(X) = g(w_{l0}^{(2)} + \sum_{k=1}^{K_1} w_{lk}^{(2)} A_k^{(1)})\\
A_k^{(1)} &= h_k^{(1)}(X) = g(w_{k0}^{(1)} + \sum_{j=1}^p w_{kj}^{(1)} X_j) \\
\end{align*}

for $l = 1, ..., K_2 = 3$ and $k = 1, ..., K_1 = 2$ and $p = 4$, where,

$$
g(z) = (z)_+ = \begin{cases}
  0, & \text{if } z < 0 \\
  z, & \text{otherwise}
\end{cases}
$$

> c. Now plug in some values for the coefficients and write out the value of
> $f(X)$.

We can perhaps achieve this most easily by fitting a real model. Note,
in the plot shown here, we also include the "bias" or intercept terms.

```{r}
library(ISLR2)
library(neuralnet)
library(sigmoid)
set.seed(5)
train <- sample(seq_len(nrow(ISLR2::Boston)), nrow(ISLR2::Boston) * 2/3)

net <- neuralnet(crim ~ lstat + medv + ptratio + rm,
    data = ISLR2::Boston[train, ],
    act.fct = relu,
    hidden = c(2, 3)
)
plot(net)
```

We can make a prediction for a given observation using this object.

Firstly, let's find an "ambiguous" test sample

```{r}
p <- predict(net, ISLR2::Boston[-train, ])
x <- ISLR2::Boston[-train, ][which.min(abs(p - mean(c(max(p), min(p))))), ]
x <- x[, c("lstat", "medv", "ptratio", "rm")]
predict(net, x)
```

Or, repeating by "hand":

```{r}
g <- function(x) ifelse(x > 0, x, 0) # relu activation function
w <- net$weights[[1]] # the estimated weights for each layer
v <- as.numeric(x) # our input predictors

# to calculate our prediction we can take the dot product of our predictors
# (with 1 at the start for the bias term) and our layer weights, lw)
for (lw in w) v <- g(c(1, v) %*% lw)
v
```

> d. How many parameters are there?

```{r} 
length(unlist(net$weights))
```

There are $4*2+2 + 2*3+3 + 3*1+1 = 23$ parameters.

### Question 2

> Consider the _softmax_ function in (10.13) (see also (4.13) on page 141)
> for modeling multinomial probabilities.
>
> a. In (10.13), show that if we add a constant $c$ to each of the $z_l$, then
> the probability is unchanged.

If we add a constant $c$ to each $Z_l$ in equation 10.13 we get:

\begin{align*}
Pr(Y=m|X) 
 &= \frac{e^{Z_m+c}}{\sum_{l=0}^9e^{Z_l+c}} \\
 &= \frac{e^{Z_m}e^c}{\sum_{l=0}^9e^{Z_l}e^c} \\
 &= \frac{e^{Z_m}e^c}{e^c\sum_{l=0}^9e^{Z_l}} \\
 &= \frac{e^{Z_m}}{\sum_{l=0}^9e^{Z_l}} \\
\end{align*}

which is just equation 10.13.

> b. In (4.13), show that if we add constants $c_j$, $j = 0,1,...,p$, to each of
> the corresponding coefficients for each of the classes, then the predictions
> at any new point $x$ are unchanged.

4.13 is 

$$
Pr(Y=k|X=x) = \frac
{e^{\beta_{K0} + \beta_{K1}x_1 + ... + \beta_{Kp}x_p}}
{\sum_{l=1}^K e^{\beta_{l0} + \beta_{l1}x1 + ... + \beta_{lp}x_p}}
$$

adding constants $c_j$ to each class gives:

\begin{align*}
Pr(Y=k|X=x) 
&= \frac
  {e^{\beta_{K0} + \beta_{K1}x_1 + c_1 + ... + \beta_{Kp}x_p + c_p}}
  {\sum_{l=1}^K e^{\beta_{l0} + \beta_{l1}x1 + c_1 + ... + \beta_{lp}x_p + c_p}} \\
&= \frac
  {e^{c1 + ... + c_p}e^{\beta_{K0} + \beta_{K1}x_1 + ... + \beta_{Kp}x_p}}
  {\sum_{l=1}^K e^{c1 + ... + c_p}e^{\beta_{l0} + \beta_{l1}x1 + ... + \beta_{lp}x_p}} \\
&= \frac
  {e^{c1 + ... + c_p}e^{\beta_{K0} + \beta_{K1}x_1 + ... + \beta_{Kp}x_p}}
  {e^{c1 + ... + c_p}\sum_{l=1}^K e^{\beta_{l0} + \beta_{l1}x1 + ... + \beta_{lp}x_p}} \\
&= \frac
  {e^{\beta_{K0} + \beta_{K1}x_1 + ... + \beta_{Kp}x_p}}
  {\sum_{l=1}^K e^{\beta_{l0} + \beta_{l1}x1 + ... + \beta_{lp}x_p}} \\
\end{align*}

which collapses to 4.13 (with the same argument as above).

> This shows that the softmax function is _over-parametrized_. However,
> regularization and SGD typically constrain the solutions so that this is not a
> problem.

### Question 3

> Show that the negative multinomial log-likelihood (10.14) is equivalent to
> the negative log of the likelihood expression (4.5) when there are $M = 2$
> classes.

Equation 10.14 is 

$$
-\sum_{i=1}^n \sum_{m=0}^9 y_{im}\log(f_m(x_i))
$$

Equation 4.5 is:

$$
\ell(\beta_0, \beta_1) = \prod_{i:y_i=1}p(x_i) \prod_{i':y_i'=0}(1-p(x_i'))
$$

So, $\log(\ell)$ is:

\begin{align*}
\log(\ell) 
 &= \log \left( \prod_{i:y_i=1}p(x_i) \prod_{i':y_i'=0}(1-p(x_i')) \right ) \\
 &= \sum_{i:y_1=1}\log(p(x_i)) + \sum_{i':y_i'=0}\log(1-p(x_i')) \\
\end{align*}

If we set $y_i$ to be an indicator variable such that $y_{i1}$ and $y_{i0}$ are
1 and 0 (or 0 and 1) when our $i$th observation is 1 (or 0) respectively, then
we can write:

$$
\log(\ell) = \sum_{i}y_{i1}\log(p(x_i)) + \sum_{i}y_{i0}\log(1-p(x_i'))
$$

If we also let $f_1(x) = p(x)$ and $f_0(x) = 1 - p(x)$ then:

\begin{align*}
\log(\ell) 
 &= \sum_i y_{i1}\log(f_1(x_i)) + \sum_{i}y_{i0}\log(f_0(x_i')) \\
 &= \sum_i \sum_{m=0}^1 y_{im}\log(f_m(x_i)) \\
\end{align*}

When we take the negative of this, it is equivalent to 10.14 for two classes 
($m = 0,1$).

### Question 4

> Consider a CNN that takes in $32 \times 32$ grayscale images and has a single
> convolution layer with three $5 \times 5$ convolution filters (without
> boundary padding).
>
> a. Draw a sketch of the input and first hidden layer similar to Figure 10.8.

```{r, echo=FALSE, out.width="50%"}
knitr::include_graphics("images/nn2.png")
```

> b. How many parameters are in this model?

There are 5 convolution matrices each with 5x5 weights (plus 5 bias terms) to
estimate, therefore 130 parameters 

> c. Explain how this model can be thought of as an ordinary feed-forward
> neural network with the individual pixels as inputs, and with constraints on
> the weights in the hidden units. What are the constraints?

We can think of a convolution layer as a regularized fully connected layer.
The regularization in this case is due to not all inputs being connected to
all outputs, and weights being shared between connections.

Each output node in the convolved image can be thought of as taking inputs from
a limited number of input pixels (the neighboring pixels), with a set of
weights specified by the convolution layer which are then shared by the
connections to all other output nodes.

> d. If there were no constraints, then how many weights would there be in the
> ordinary feed-forward neural network in (c)?

With no constraints, we would connect each output pixel in our 5x32x32
convolution layer to each node in the 32x32 original image (plus 5 bias terms),
giving a total of 5,242,885 weights to estimate.

### Question 5

> In Table 10.2 on page 433, we see that the ordering of the three methods with
> respect to mean absolute error is different from the ordering with respect to
> test set $R^2$. How can this be?

Mean absolute error considers _absolute_ differences between predictions and 
observed values, whereas $R^2$ considers the (normalized) sum of _squared_
differences, thus larger errors contribute relatively ore to $R^2$ than mean
absolute error.

## Applied

### Question 6

> Consider the simple function $R(\beta) = sin(\beta) + \beta/10$.
>
> a. Draw a graph of this function over the range $\beta \in [−6, 6]$.

```{r}
r <- function(x) sin(x) + x/10
x <- seq(-6, 6, 0.1)
plot(x, r(x), type = "l")
```

> b. What is the derivative of this function?

$$
cos(x) + 1/10
$$

> c. Given $\beta^0 = 2.3$, run gradient descent to find a local minimum of
> $R(\beta)$ using a learning rate of $\rho = 0.1$. Show each of 
> $\beta^0, \beta^1, ...$ in your plot, as well as the final answer.

The derivative of our function, i.e. $cos(x) + 1/10$ gives us the gradient for
a given $x$. For gradient descent, we move $x$ a little in the _opposite_
direction, for some learning rate $\rho = 0.1$:

$$
x^{m+1} = x^m - \rho (cos(x^m) + 1/10)
$$

```{r}
iter <- function(x, rho) x - rho*(cos(x) + 1/10)
gd <- function(start, rho = 0.1) {
  b <- start
  v <- b
  while(abs(b - iter(b, 0.1)) > 1e-8) {
    b <- iter(b, 0.1)
    v <- c(v, b)
  }
  v
}

res <- gd(2.3)
res[length(res)]
```

```{r}
plot(x, r(x), type = "l")
points(res, r(res), col = "red", pch = 19)
```


> d. Repeat with $\beta^0 = 1.4$.

```{r}
res <- gd(1.4)
res[length(res)]
```

```{r}
plot(x, r(x), type = "l")
points(res, r(res), col = "red", pch = 19)
```

### Question 7

> Fit a neural network to the `Default` data. Use a single hidden layer with 10
> units, and dropout regularization. Have a look at Labs 10.9.1–-10.9.2 for
> guidance. Compare the classification performance of your model with that of
> linear logistic regression.

```{r, cache = TRUE}
library(keras)

dat <- ISLR2::Boston
x <- scale(model.matrix(crim ~ . - 1, data = dat))
n <- nrow(dat)
ntest <- trunc(n / 3)
testid <- sample(1:n, ntest)
y <- dat$crim

# logistic regression
lfit <- lm(crim ~ ., data = dat[-testid, ])
lpred <- predict(lfit, dat[testid, ])
with(dat[testid, ], mean(abs(lpred - crim)))

# keras
nn <- keras_model_sequential() |>
  layer_dense(units = 10, activation = "relu", input_shape = ncol(x)) |>
  layer_dropout(rate = 0.4) |>
  layer_dense(units = 1)

compile(nn, loss = "mse", 
  optimizer = optimizer_rmsprop(), 
  metrics = list("mean_absolute_error") 
)

history <- fit(nn,
  x[-testid, ], y[-testid], 
  epochs = 100, 
  batch_size = 26, 
  validation_data = list(x[testid, ], y[testid]),
  verbose = 0
)
plot(history, smooth = FALSE)
npred <- predict(nn, x[testid, ])
mean(abs(y[testid] - npred))
```

In this case, the neural network outperforms logistic regression having a lower
absolute error rate on the test data.

### Question 8

> From your collection of personal photographs, pick 10 images of animals (such
> as dogs, cats, birds, farm animals, etc.). If the subject does not occupy a
> reasonable part of the image, then crop the image. Now use a pretrained image
> classification CNN as in Lab 10.9.4 to predict the class of each of your
> images, and report the probabilities for the top five predicted classes for
> each image.

```{r, echo=FALSE}
knitr::include_graphics(c(
  "images/animals/bird.jpg",
  "images/animals/bird2.jpg",
  "images/animals/bird3.jpg",
  "images/animals/bug.jpg",
  "images/animals/butterfly.jpg",
  "images/animals/butterfly2.jpg",
  "images/animals/elba.jpg",
  "images/animals/hamish.jpg",
  "images/animals/poodle.jpg",
  "images/animals/tortoise.jpg"
))
```

```{r}
library(keras)
images <- list.files("images/animals")
x <- array(dim = c(length(images), 224, 224, 3))
for (i in seq_len(length(images))) {
  img <- image_load(paste0("images/animals/", images[i]), target_size = c(224, 224))
  x[i,,,] <- image_to_array(img)
}

model <- application_resnet50(weights = "imagenet")

pred <- model |>
  predict(x) |>
  imagenet_decode_predictions(top = 5)
  
names(pred) <- images
print(pred)
```

### Question 9

> Fit a lag-5 autoregressive model to the `NYSE` data, as described in the text
> and Lab 10.9.6. Refit the model with a 12-level factor representing the
> month. Does this factor improve the performance of the model?

Fitting the model as described in the text.

```{r}
library(tidyverse)
library(ISLR2)
xdata <- data.matrix(NYSE[, c("DJ_return", "log_volume","log_volatility")])
istrain <- NYSE[, "train"]
xdata <- scale(xdata)

lagm <- function(x, k = 1) {
  n <- nrow(x)
  pad <- matrix(NA, k, ncol(x))
  rbind(pad, x[1:(n - k), ])
}

arframe <- data.frame(
  log_volume = xdata[, "log_volume"], 
  L1 = lagm(xdata, 1), 
  L2 = lagm(xdata, 2),
  L3 = lagm(xdata, 3),
  L4 = lagm(xdata, 4),
  L5 = lagm(xdata, 5)
)

arframe <- arframe[-(1:5), ]
istrain <- istrain[-(1:5)]

arfit <- lm(log_volume ~ ., data = arframe[istrain, ])
arpred <- predict(arfit, arframe[!istrain, ])
V0 <- var(arframe[!istrain, "log_volume"])
1 - mean((arpred - arframe[!istrain, "log_volume"])^2) / V0
```

Now we add month (and work with tidyverse).

```{r}
arframe$month = as.factor(str_match(NYSE$date, "-(\\d+)-")[,2])[-(1:5)]
arfit2 <- lm(log_volume ~ ., data = arframe[istrain, ])
arpred2 <- predict(arfit2, arframe[!istrain, ])
V0 <- var(arframe[!istrain, "log_volume"])
1 - mean((arpred2 - arframe[!istrain, "log_volume"])^2) / V0
```

Adding month as a factor marginally improves the $R^2$ of our model (from 
0.413223 to 0.4170418). This is a significant improvement in fit and model
2 has a lower AIC.

```{r}
anova(arfit, arfit2)
AIC(arfit, arfit2)
```

### Question 10

> In Section 10.9.6, we showed how to fit a linear AR model to the `NYSE` data
> using the `lm()` function. However, we also mentioned that we can "flatten"
> the short sequences produced for the RNN model in order to fit a linear AR
> model. Use this latter approach to fit a linear AR model to the NYSE data.
> Compare the test $R^2$ of this linear AR model to that of the linear AR model
> that we fit in the lab. What are the advantages/disadvantages of each
> approach?

The `lm` model is the same as that fit above:

```{r}
arfit <- lm(log_volume ~ ., data = arframe[istrain, ])
arpred <- predict(arfit, arframe[!istrain, ])
V0 <- var(arframe[!istrain, "log_volume"])
1 - mean((arpred - arframe[!istrain, "log_volume"])^2) / V0
```

Now we reshape the data for the RNN

```{r}
n <- nrow(arframe)
xrnn <- data.matrix(arframe[, -1])
xrnn <- array(xrnn, c(n, 3, 5))
xrnn <- xrnn[, , 5:1]
xrnn <- aperm(xrnn, c(1, 3, 2))
```

We can add a "flatten" layer to turn the reshaped data into a long vector of
predictors resulting in a linear AR model.

```{r}
model <- keras_model_sequential() |>
  layer_flatten(input_shape = c(5, 3)) |>
  layer_dense(units = 1)
```

Now let's fit this model.

```{r}
model |>
  compile(optimizer = optimizer_rmsprop(), loss = "mse")

history <- model |>
  fit(
    xrnn[istrain,, ],
    arframe[istrain, "log_volume"],
    batch_size = 64,
    epochs = 200,
    validation_data = list(xrnn[!istrain,, ], arframe[!istrain, "log_volume"]),
    verbose = 0
  )

plot(history, smooth = FALSE)
kpred <- predict(model, xrnn[!istrain,, ])
1 - mean((kpred - arframe[!istrain, "log_volume"])^2) / V0
```

Both models estimate the same number of coefficients/weights (16):

```{r}
coef(arfit)
model$get_weights()
```

The flattened RNN has a lower $R^2$ on the test data than our `lm` model
above. The `lm` model is quicker to fit and conceptually simpler also 
giving us the ability to inspect the coefficients for different variables.

The flattened RNN is regularized to some extent as data are processed in
batches.

### Question 11

> Repeat the previous exercise, but now fit a nonlinear AR model by "flattening"
> the short sequences produced for the RNN model.

From the book:

> To fit a nonlinear AR model, we could add in a hidden layer.

```{r, c10q11}
xfun::cache_rds({

  model <- keras_model_sequential() |> 
    layer_flatten(input_shape = c(5, 3)) |>
    layer_dense(units = 32, activation = "relu") |>
    layer_dropout(rate = 0.4) |> 
    layer_dense(units = 1)

  model |> compile(
    loss = "mse", 
    optimizer = optimizer_rmsprop(), 
    metrics = "mse"
  )

  history <- model |>
    fit(
      xrnn[istrain,, ],
      arframe[istrain, "log_volume"],
      batch_size = 64,
      epochs = 200,
      validation_data = list(xrnn[!istrain,, ], arframe[!istrain, "log_volume"]),
      verbose = 0
    )

  plot(history, smooth = FALSE, metrics = "mse")
  kpred <- predict(model, xrnn[!istrain,, ])
  1 - mean((kpred - arframe[!istrain, "log_volume"])^2) / V0

})
```

This approach improves our $R^2$ over the linear model above.

### Question 12

> Consider the RNN fit to the `NYSE` data in Section 10.9.6. Modify the code to
> allow inclusion of the variable `day_of_week`, and fit the RNN. Compute the
> test $R^2$.

To accomplish this, I'll include day of the week as one of the lagged variables
in the RNN. Thus, our input for each observation will be 4 x 5 (rather than
3 x 5).

```{r, c10q12}
xfun::cache_rds({
  xdata <- data.matrix(
    NYSE[, c("day_of_week", "DJ_return", "log_volume","log_volatility")] 
  )
  istrain <- NYSE[, "train"]
  xdata <- scale(xdata)

  arframe <- data.frame(
    log_volume = xdata[, "log_volume"], 
    L1 = lagm(xdata, 1),
    L2 = lagm(xdata, 2),
    L3 = lagm(xdata, 3), 
    L4 = lagm(xdata, 4),
    L5 = lagm(xdata, 5)
  )
  arframe <- arframe[-(1:5), ]
  istrain <- istrain[-(1:5)]

  n <- nrow(arframe)
  xrnn <- data.matrix(arframe[, -1])
  xrnn <- array(xrnn, c(n, 4, 5))
  xrnn <- xrnn[,, 5:1]
  xrnn <- aperm(xrnn, c(1, 3, 2))
  dim(xrnn)

  model <- keras_model_sequential() |>
      layer_simple_rnn(units = 12,
      input_shape = list(5, 4),
      dropout = 0.1, 
      recurrent_dropout = 0.1
    ) |>
    layer_dense(units = 1)

  model |> compile(optimizer = optimizer_rmsprop(), loss = "mse")

  history <- model |> 
    fit(
      xrnn[istrain,, ],
      arframe[istrain, "log_volume"],
      batch_size = 64,
      epochs = 200,
      validation_data = list(xrnn[!istrain,, ], arframe[!istrain, "log_volume"]),
      verbose = 0
  )

  kpred <- predict(model, xrnn[!istrain,, ])
  1 - mean((kpred - arframe[!istrain, "log_volume"])^2) / V0

})
```

### Question 13

> Repeat the analysis of Lab 10.9.5 on the `IMDb` data using a similarly
> structured neural network. There we used a dictionary of size 10,000. Consider
> the effects of varying the dictionary size. Try the values 1000, 3000, 5000,
> and 10,000, and compare the results.

```{r, c10q13}
xfun::cache_rds({
  library(knitr)
  accuracy <- c()
  for(max_features in c(1000, 3000, 5000, 10000)) {
    imdb <- dataset_imdb(num_words = max_features)
    c(c(x_train, y_train), c(x_test, y_test)) %<-% imdb

    maxlen <- 500
    x_train <- pad_sequences(x_train, maxlen = maxlen)
    x_test <- pad_sequences(x_test, maxlen = maxlen)

    model <- keras_model_sequential() |>
      layer_embedding(input_dim = max_features, output_dim = 32) |>
      layer_lstm(units = 32) |>
      layer_dense(units = 1, activation = "sigmoid")

    model |> compile(
      optimizer = "rmsprop",
      loss = "binary_crossentropy", 
      metrics = "acc"
    )

    history <- fit(model, x_train, y_train, 
      epochs = 10, 
      batch_size = 128, 
      validation_data = list(x_test, y_test),
      verbose = 0
    )

    predy <- predict(model, x_test) > 0.5
    accuracy <- c(accuracy, mean(abs(y_test == as.numeric(predy))))
  }

  tibble(
    "Max Features" = c(1000, 3000, 5000, 10000),
    "Accuracy" = accuracy
  ) |>
    kable()

})
```

Varying the dictionary size does not make a substantial impact on our estimates
of accuracy. However, the models do take a substantial amount of time to fit and
it is not clear we are finding the best fitting models in each case. For
example, the model using a dictionary size of 10,000 obtained an accuracy of
0.8721 in the text which is as different from the estimate obtained here as
are the differences between the models with different dictionary sizes.