{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "015a3955", "metadata": { "lines_to_next_cell": 0, "name": "setup", "tags": [ "remove_input" ] }, "outputs": [], "source": [ "knitr::opts_chunk$set(error = TRUE)" ] }, { "cell_type": "markdown", "id": "06f0bdb7", "metadata": {}, "source": [ "# Deep Learning\n", "In this section, we show how to fit the examples discussed in the text. We use the `keras` package, which interfaces to the\n", "`tensorflow` package which in turn links to efficient\n", "`python` code. This code is\n", "impressively fast, and the package is well-structured. A good companion\n", "is the text *Deep Learning with R* (F. Chollet and J.J. Allaire, *Deep Learning with R* (2018), Manning Publications.), and most of our code is adapted from\n", "there.\n", "\n", "Getting `keras` up and running on your computer can be a challenge.\n", " The book website gives step-by-step\n", " instructions on how to achieve this. (Many thanks to Balasubramanian Narasimhan for preparing the `keras` installation instructions.)\n", "Guidance can also be found at .\n", "\n", "Since the first printing of this book, the `torch` package has become available as an alternative to the `keras` package for deep learning. While `torch` does not require a `python` installation, the current implementation appears to be a fair bit slower than `keras`. A version of this lab that makes use of `torch` is available on the book website. (Many thanks to Daniel Falbel and Sigrid Keydana for preparing the `torch` version of this lab.)\n", "\n" ] }, { "cell_type": "markdown", "id": "fea5b3a6", "metadata": {}, "source": [ "## A Single Layer Network on the Hitters Data\n", "We start by fitting the models in Section 10.6.\n", "We set up the data, and separate out a training and\n", "test set." ] }, { "cell_type": "code", "execution_count": null, "id": "dd2615a4", "metadata": { "lines_to_next_cell": 2, "name": "chunk1" }, "outputs": [], "source": [ "library(ISLR2)\n", "Gitters <- na.omit(Hitters)\n", "n <- nrow(Gitters)\n", "set.seed(13)\n", "ntest <- trunc(n / 3)\n", "testid <- sample(1:n, ntest)" ] }, { "cell_type": "markdown", "id": "4fd8f578", "metadata": {}, "source": [ "The linear model should be familiar, but we present it anyway." ] }, { "cell_type": "code", "execution_count": null, "id": "a9896c02", "metadata": { "name": "chunk2" }, "outputs": [], "source": [ "lfit <- lm(Salary ~ ., data = Gitters[-testid, ])\n", "lpred <- predict(lfit, Gitters[testid, ])\n", "with(Gitters[testid, ], mean(abs(lpred - Salary)))" ] }, { "cell_type": "markdown", "id": "7e68da0e", "metadata": {}, "source": [ "Notice the use of the `with()` command: the first argument is a\n", "dataframe, and the second an expression that can refer to\n", "elements of the dataframe by name. In this instance the dataframe corresponds to the test data and the expression computes the mean absolute prediction error on this data.\n", "\n", "Next we fit the lasso using `glmnet`. Since this package does\n", "not use formulas, we create `x` and `y` first." ] }, { "cell_type": "code", "execution_count": null, "id": "626584a2", "metadata": { "name": "chunk3" }, "outputs": [], "source": [ "x <- scale(model.matrix(Salary ~ . - 1, data = Gitters))\n", "y <- Gitters$Salary" ] }, { "cell_type": "markdown", "id": "de5c1b07", "metadata": {}, "source": [ "The first line makes a call\n", "to `model.matrix()`, which produces the same matrix that was\n", "used by `lm()` (the `-1` omits the intercept).\n", "This function automatically converts factors to dummy variables.\n", "The `scale()` function standardizes the matrix so each column\n", "has mean zero and variance one." ] }, { "cell_type": "code", "execution_count": null, "id": "4ecbb2d4", "metadata": { "name": "chunk4" }, "outputs": [], "source": [ "library(glmnet)\n", "cvfit <- cv.glmnet(x[-testid, ], y[-testid],\n", " type.measure = \"mae\")\n", "cpred <- predict(cvfit, x[testid, ], s = \"lambda.min\")\n", "mean(abs(y[testid] - cpred))" ] }, { "cell_type": "markdown", "id": "fc6edb73", "metadata": {}, "source": [ "To fit the neural network, we first set up a model structure\n", "that describes the network." ] }, { "cell_type": "code", "execution_count": null, "id": "4567679a", "metadata": { "name": "chunk5" }, "outputs": [], "source": [ "library(keras)\n", "modnn <- keras_model_sequential() %>%\n", " layer_dense(units = 50, activation = \"relu\",\n", " input_shape = ncol(x)) %>%\n", " layer_dropout(rate = 0.4) %>%\n", " layer_dense(units = 1)" ] }, { "cell_type": "markdown", "id": "6eb8f4a9", "metadata": {}, "source": [ "We have created a vanilla model\n", "object called `modnn`, and have added details about the successive\n", "layers in a sequential manner, using the function `keras_model_sequential()`.\n", "The *pipe* operator `\\%>\\%`\n", " passes the previous term as the first argument to the next\n", "function, and returns the result.\n", "It allows us to specify the layers of\n", "a neural network in a readable form.\n", "\n", "We illustrate the use of the pipe operator on a simple example. Earlier, we created `x` using the command" ] }, { "cell_type": "code", "execution_count": null, "id": "bbba68cc", "metadata": { "name": "chunk6" }, "outputs": [], "source": [ "x <- scale(model.matrix(Salary ~ . - 1, data = Gitters))" ] }, { "cell_type": "markdown", "id": "51bd13a2", "metadata": {}, "source": [ "We first make a matrix, and then we center each of the variables.\n", "Compound expressions like this can be difficult to parse. We could have obtained the same result using the pipe operator:" ] }, { "cell_type": "code", "execution_count": null, "id": "f85daf35", "metadata": { "name": "chunk7" }, "outputs": [], "source": [ "x <- model.matrix(Salary ~ . - 1, data = Gitters) %>% scale()" ] }, { "cell_type": "markdown", "id": "205af290", "metadata": {}, "source": [ "Using the pipe operator makes it easier to follow the sequence of operations.\n", "\n", "We now return to our neural network. The object `modnn` has a single hidden layer with 50 hidden units, and\n", "a ReLU activation function. It then has a dropout layer, in which a\n", "random 40\\% of the 50 activations from the previous layer are set to zero\n", "during each iteration of the stochastic gradient descent\n", "algorithm. Finally, the output layer has just one unit with no\n", "activation function, indicating that the model provides a single\n", "quantitative output." ] }, { "cell_type": "markdown", "id": "f5de523f", "metadata": {}, "source": [ "Next we add details to `modnn` that control the fitting\n", "algorithm. Here we have simply followed the examples given in\n", "the Keras book. We minimize squared-error loss as in\n", "(10.23). The algorithm\n", "tracks the mean absolute error on the training data, and\n", "on validation data if it is supplied." ] }, { "cell_type": "code", "execution_count": null, "id": "1bae2bc5", "metadata": { "name": "chunk8" }, "outputs": [], "source": [ "modnn %>% compile(loss = \"mse\",\n", " optimizer = optimizer_rmsprop(),\n", " metrics = list(\"mean_absolute_error\")\n", " )" ] }, { "cell_type": "markdown", "id": "ae2fb03a", "metadata": {}, "source": [ "In the previous line, the pipe operator passes `modnn` as the first argument to `compile()`. The `compile()` function does not actually change the `R`\n", "object `modnn`, but it does communicate these specifications to\n", "the corresponding `python` instance of this model that has been created\n", "along the way." ] }, { "cell_type": "markdown", "id": "6fdef655", "metadata": {}, "source": [ "Now we fit the model. We supply the training data and two\n", "fitting parameters, `epochs` and `batch_size`. Using 32\n", "for the latter means that at each step of SGD, the algorithm randomly\n", "selects 32 training observations for the computation of the gradient. Recall from\n", "Sections 10.4 and 10.7\n", "that an epoch amounts to the number of SGD steps required to process $n$\n", "observations. Since the training set has\n", "$n=176$, an epoch is $176/32=5.5$ SGD steps. The `fit()` function has an argument\n", "`validation_data`; these data are not used in the fitting,\n", "but can be used to track the progress of the model (in this case\n", "reporting the \n", "mean absolute error). Here we\n", "actually supply the test data so we can see the mean absolute error of both the\n", "training data and test data as the epochs proceed. To see more options\n", "for fitting, use `?fit.keras.engine.training.Model`." ] }, { "cell_type": "code", "execution_count": null, "id": "fc5fbbe1", "metadata": { "name": "chunk9" }, "outputs": [], "source": [ "history <- modnn %>% fit(\n", " x[-testid, ], y[-testid], epochs = 1500, batch_size = 32,\n", " validation_data = list(x[testid, ], y[testid])\n", " )" ] }, { "cell_type": "markdown", "id": "3f1c66ab", "metadata": {}, "source": [ "We can plot the `history` to display the mean absolute error for the training and test data. For the best aesthetics, install the `ggplot2` package before calling the `plot()` function. If you have not installed `ggplot2`, then the code below will still run, but the plot will be less attractive." ] }, { "cell_type": "code", "execution_count": null, "id": "6ac54e83", "metadata": { "name": "chunk10" }, "outputs": [], "source": [ "plot(history)" ] }, { "cell_type": "markdown", "id": "9b6b6790", "metadata": {}, "source": [ "It is worth noting that if you run the `fit()` command a second time in the same \\text{R} session, then the fitting process will pick up where it left off. Try re-running the `fit()` command, and then the `plot()` command, to see!\n", "\n", "Finally, we predict from the final model, and\n", "evaluate its performance on the test data. Due to the use of SGD, the results vary slightly with each\n", "fit. Unfortunately the `set.seed()` function does not ensure\n", "identical results (since the fitting is done in `python`), so your results will differ slightly." ] }, { "cell_type": "code", "execution_count": null, "id": "667b354f", "metadata": { "lines_to_next_cell": 2, "name": "chunk11" }, "outputs": [], "source": [ "npred <- predict(modnn, x[testid, ])\n", "mean(abs(y[testid] - npred))" ] }, { "cell_type": "markdown", "id": "436e1d7f", "metadata": {}, "source": [ "## A Multilayer Network on the MNIST Digit Data\n", "The `keras` package comes with a number of example datasets,\n", "including the `MNIST` digit data. Our first step is to load the\n", "`MNIST` data. The `dataset_mnist()` function is provided for this purpose." ] }, { "cell_type": "code", "execution_count": null, "id": "096901ae", "metadata": { "name": "chunk12" }, "outputs": [], "source": [ "mnist <- dataset_mnist()\n", "x_train <- mnist$train$x\n", "g_train <- mnist$train$y\n", "x_test <- mnist$test$x\n", "g_test <- mnist$test$y\n", "dim(x_train)\n", "dim(x_test)" ] }, { "cell_type": "markdown", "id": "0f6d8fc9", "metadata": {}, "source": [ "There are 60,000 images in the training data and 10,000 in the test data. The images are $28\\times 28$, and stored as a three-dimensional array,\n", "so we need to reshape them into a matrix. Also, we need to \"one-hot\"\n", "encode the class label. Luckily `keras` has a lot of built-in\n", "functions that do this for us." ] }, { "cell_type": "code", "execution_count": null, "id": "491f2d5a", "metadata": { "name": "chunk13" }, "outputs": [], "source": [ "x_train <- array_reshape(x_train, c(nrow(x_train), 784))\n", "x_test <- array_reshape(x_test, c(nrow(x_test), 784))\n", "y_train <- to_categorical(g_train, 10)\n", "y_test <- to_categorical(g_test, 10)" ] }, { "cell_type": "markdown", "id": "8de6c21f", "metadata": {}, "source": [ "Neural networks are somewhat sensitive to the scale of the inputs. For example, ridge and\n", "lasso regularization are affected by scaling. Here the inputs are eight-bit (Eight bits means $2^8$, which equals 256. Since the convention is to start at $0$, the possible values range from $0$ to $255$.)\n", "grayscale values between 0 and 255, so we rescale to the unit\n", "interval." ] }, { "cell_type": "code", "execution_count": null, "id": "7ca5b8ed", "metadata": { "lines_to_next_cell": 2, "name": "chunk14" }, "outputs": [], "source": [ "x_train <- x_train / 255\n", "x_test <- x_test / 255" ] }, { "cell_type": "markdown", "id": "d0a7011d", "metadata": {}, "source": [ "Now we are ready to fit our neural network." ] }, { "cell_type": "code", "execution_count": null, "id": "43945f5a", "metadata": { "name": "chunk15" }, "outputs": [], "source": [ "modelnn <- keras_model_sequential()\n", "modelnn %>%\n", " layer_dense(units = 256, activation = \"relu\",\n", " input_shape = c(784)) %>%\n", " layer_dropout(rate = 0.4) %>%\n", " layer_dense(units = 128, activation = \"relu\") %>%\n", " layer_dropout(rate = 0.3) %>%\n", " layer_dense(units = 10, activation = \"softmax\")" ] }, { "cell_type": "markdown", "id": "27cbb05d", "metadata": {}, "source": [ "The first layer goes from $28\\times28=784$ input units to a hidden\n", "layer of $256$ units, which uses the ReLU activation function.\n", "This is specified by a call to `layer_dense()`, which takes as\n", "input a `modelnn` object, and returns a modified `modelnn`\n", "object.\n", "This is then piped through `layer_dropout()` to perform\n", "dropout regularization. The second hidden layer comes\n", "next, with $128$ hidden units, followed by a dropout layer.\n", "The final layer is the output layer, with activation\n", "`\"softmax\"` (10.13) for the 10-class classification\n", "problem, which defines the map from the second hidden layer to class\n", "probabilities.\n", "Finally, we use `summary()` to summarize the model, and to make sure we got it\n", "all right." ] }, { "cell_type": "code", "execution_count": null, "id": "cfd42c26", "metadata": { "name": "chunk16" }, "outputs": [], "source": [ "summary(modelnn)" ] }, { "cell_type": "markdown", "id": "e28828f1", "metadata": {}, "source": [ "The parameters for each layer include a bias term, which results in a\n", "parameter count of 235,146. For example, the first hidden\n", "layer involves $(784+1)\\times256=200{,}960$ parameters.\n", "\n", "Notice that the layer names such as `dropout_1` and\n", "`dense_2` have subscripts. These may appear somewhat random; in\n", "fact, if you fit the same model again, these will change. They are of\n", "no consequence: they vary because the model\n", "specification code is run in `python`, and these subscripts are incremented every time\n", "`keras_model_sequential()` is called.\n", "\n", "Next, we add details to the model to specify the fitting algorithm. We fit the model by minimizing the cross-entropy function given by (10.14)." ] }, { "cell_type": "code", "execution_count": null, "id": "b2dc5034", "metadata": { "lines_to_next_cell": 2, "name": "chunk17" }, "outputs": [], "source": [ "modelnn %>% compile(loss = \"categorical_crossentropy\",\n", " optimizer = optimizer_rmsprop(), metrics = c(\"accuracy\")\n", " )" ] }, { "cell_type": "markdown", "id": "495841a6", "metadata": {}, "source": [ "Now we are ready to go. The final step is to supply training data, and fit the model." ] }, { "cell_type": "code", "execution_count": null, "id": "c0d65a6d", "metadata": { "name": "chunk18" }, "outputs": [], "source": [ "system.time(\n", " history <- modelnn %>%\n", " fit(x_train, y_train, epochs = 30, batch_size = 128,\n", " validation_split = 0.2)\n", ")\n", "plot(history, smooth = FALSE)" ] }, { "cell_type": "markdown", "id": "fe243b6a", "metadata": {}, "source": [ "We have suppressed the output here, which is a progress report on the\n", "fitting of the model, grouped by epoch. This is very useful, since on\n", "large datasets fitting can take time. Fitting this model took 144\n", "seconds on a 2.9 GHz MacBook Pro with 4 cores and 32 GB of RAM.\n", "Here we specified a\n", "validation split of 20\\%, so the training is actually performed on\n", "80\\% of the 60,000 observations in the training set. This is an\n", "alternative to actually supplying validation data, like we did in\n", "Section 10.9.1. See\n", "`?fit.keras.engine.training.Model` for all the optional fitting arguments. SGD uses batches\n", "of 128 observations in computing the gradient, and doing the\n", "arithmetic, we see that an epoch corresponds to 375 gradient steps.\n", "The last `plot()` command produces a figure similar to Figure 10.18.\n", "\n", "To obtain the test error in Table 10.1, we first write\n", "a simple function `accuracy()` that compares predicted and true\n", "class labels, and then use it to evaluate our predictions." ] }, { "cell_type": "code", "execution_count": null, "id": "f565c590", "metadata": { "name": "chunk19" }, "outputs": [], "source": [ "accuracy <- function(pred, truth)\n", " mean(drop(as.numeric(pred)) == drop(truth))\n", "modelnn %>% predict(x_test) %>% k_argmax() %>% accuracy(g_test)" ] }, { "cell_type": "markdown", "id": "9de0ba2b", "metadata": {}, "source": [ "The table also reports LDA (Chapter 4) and multiclass logistic\n", "regression. Although packages such as `glmnet` can handle\n", "multiclass logistic regression, they are quite slow on this large\n", "dataset. It is much faster and quite easy to fit such a model\n", "using the `keras` software. We just have an input layer and output layer, and omit the hidden layers!" ] }, { "cell_type": "code", "execution_count": null, "id": "25e0a201", "metadata": { "name": "chunk20" }, "outputs": [], "source": [ "modellr <- keras_model_sequential() %>%\n", " layer_dense(input_shape = 784, units = 10,\n", " activation = \"softmax\")\n", "summary(modellr)" ] }, { "cell_type": "markdown", "id": "8cd941b2", "metadata": {}, "source": [ "We fit the model just as before." ] }, { "cell_type": "code", "execution_count": null, "id": "7ff9d995", "metadata": { "lines_to_next_cell": 2, "name": "chunk21" }, "outputs": [], "source": [ "modellr %>% compile(loss = \"categorical_crossentropy\",\n", " optimizer = optimizer_rmsprop(), metrics = c(\"accuracy\"))\n", "modellr %>% fit(x_train, y_train, epochs = 30,\n", " batch_size = 128, validation_split = 0.2)\n", "modellr %>% predict(x_test) %>% k_argmax() %>% accuracy(g_test)" ] }, { "cell_type": "markdown", "id": "adaa3905", "metadata": {}, "source": [ "## Convolutional Neural Networks\n", "In this section we fit a CNN to the `CIFAR` data, which is available in the `keras`\n", "package. It is arranged in a similar fashion as the `MNIST` data." ] }, { "cell_type": "code", "execution_count": null, "id": "317363b9", "metadata": { "name": "chunk22" }, "outputs": [], "source": [ "cifar100 <- dataset_cifar100()\n", "names(cifar100)\n", "x_train <- cifar100$train$x\n", "g_train <- cifar100$train$y\n", "x_test <- cifar100$test$x\n", "g_test <- cifar100$test$y\n", "dim(x_train)\n", "range(x_train[1,,, 1])" ] }, { "cell_type": "markdown", "id": "6a69c19a", "metadata": {}, "source": [ "The array of 50,000 training images has four dimensions:\n", " each three-color image is represented as a set of three channels, each of which consists of\n", "$32\\times 32$ eight-bit pixels. We standardize as we did for the\n", "digits, but keep the array structure. We one-hot encode the response\n", "factors to produce a 100-column binary matrix." ] }, { "cell_type": "code", "execution_count": null, "id": "a39d82f5", "metadata": { "lines_to_next_cell": 2, "name": "chunk23" }, "outputs": [], "source": [ "x_train <- x_train / 255\n", "x_test <- x_test / 255\n", "y_train <- to_categorical(g_train, 100)\n", "dim(y_train)" ] }, { "cell_type": "markdown", "id": "6352d97f", "metadata": {}, "source": [ "Before we start, we look at some of the training images using the `jpeg` package; similar code produced\n", "Figure 10.5 on page 411." ] }, { "cell_type": "code", "execution_count": null, "id": "786c0cb2", "metadata": { "name": "chunk24" }, "outputs": [], "source": [ "library(jpeg)\n", "par(mar = c(0, 0, 0, 0), mfrow = c(5, 5))\n", "index <- sample(seq(50000), 25)\n", "for (i in index) plot(as.raster(x_train[i,,, ]))" ] }, { "cell_type": "markdown", "id": "ceb48a7e", "metadata": {}, "source": [ "The `as.raster()` function converts the feature map so that it can be plotted as a color image.\n", "\n", "Here we specify a moderately-sized CNN for\n", "demonstration purposes, similar in structure to Figure 10.8." ] }, { "cell_type": "code", "execution_count": null, "id": "a22c648c", "metadata": { "name": "chunk25" }, "outputs": [], "source": [ "model <- keras_model_sequential() %>%\n", " layer_conv_2d(filters = 32, kernel_size = c(3, 3),\n", " padding = \"same\", activation = \"relu\",\n", " input_shape = c(32, 32, 3)) %>%\n", " layer_max_pooling_2d(pool_size = c(2, 2)) %>%\n", " layer_conv_2d(filters = 64, kernel_size = c(3, 3),\n", " padding = \"same\", activation = \"relu\") %>%\n", " layer_max_pooling_2d(pool_size = c(2, 2)) %>%\n", " layer_conv_2d(filters = 128, kernel_size = c(3, 3),\n", " padding = \"same\", activation = \"relu\") %>%\n", " layer_max_pooling_2d(pool_size = c(2, 2)) %>%\n", " layer_conv_2d(filters = 256, kernel_size = c(3, 3),\n", " padding = \"same\", activation = \"relu\") %>%\n", " layer_max_pooling_2d(pool_size = c(2, 2)) %>%\n", " layer_flatten() %>%\n", " layer_dropout(rate = 0.5) %>%\n", " layer_dense(units = 512, activation = \"relu\") %>%\n", " layer_dense(units = 100, activation = \"softmax\")\n", "summary(model)" ] }, { "cell_type": "markdown", "id": "72c355b7", "metadata": {}, "source": [ "Notice that we used the `padding = \"same\"` argument to\n", "`layer_conv_2D()`, which ensures that the output channels have the\n", "same dimension as the input channels. There are 32 channels in the first\n", "hidden layer, in contrast to the three channels in the input layer. We\n", "use a $3\\times 3$ convolution filter for each channel in all the layers. Each\n", "convolution is followed by a max-pooling layer over $2\\times2$\n", "blocks. By studying the summary, we can see that the channels halve in both\n", "dimensions\n", "after each of these max-pooling operations. After the last of these we\n", "have a layer with 256 channels of dimension $2\\times 2$. These are then\n", "flattened to a dense layer of size 1,024:\n", "in other words, each of the $2\\times 2$ matrices is turned into a\n", "$4$-vector, and put side-by-side in one layer. This is followed by a\n", "dropout regularization layer, then\n", "another dense layer of size 512, which finally reaches the softmax\n", "output layer.\n", "\n", "Finally, we specify the fitting algorithm, and fit the model." ] }, { "cell_type": "code", "execution_count": null, "id": "a8d1ceb3", "metadata": { "name": "chunk26" }, "outputs": [], "source": [ "model %>% compile(loss = \"categorical_crossentropy\",\n", " optimizer = optimizer_rmsprop(), metrics = c(\"accuracy\"))\n", "history <- model %>% fit(x_train, y_train, epochs = 30,\n", " batch_size = 128, validation_split = 0.2)\n", "model %>% predict(x_test) %>% k_argmax() %>% accuracy(g_test)" ] }, { "cell_type": "markdown", "id": "f46924ec", "metadata": {}, "source": [ "This model takes 10 minutes to run and achieves 46\\% accuracy on the test\n", "data. Although this is not terrible for 100-class data (a random\n", "classifier gets 1\\% accuracy), searching the web we see results around\n", "75\\%. Typically it takes a lot of architecture carpentry,\n", "fiddling with regularization, and time to achieve such results.\n", "\n", "## Using Pretrained CNN Models\n", "We now show how to use a CNN pretrained on the `imagenet` database to classify natural\n", "images, and demonstrate how we produced Figure 10.10.\n", "We copied six jpeg images from a digital photo album into the\n", "directory `book_images`. (These images are available from the data section of , the ISL book website. Download **book_images.zip**; when clicked it creates the **book_images** directory.) We first read in the images, and\n", "convert them into the array format expected by the `keras`\n", "software to match the specifications in `imagenet`. Make sure that your working directory in `R` is set to the folder in which the images are stored." ] }, { "cell_type": "code", "execution_count": null, "id": "73d408ff", "metadata": { "name": "chunk27" }, "outputs": [], "source": [ "img_dir <- \"book_images\"\n", "image_names <- list.files(img_dir)\n", "num_images <- length(image_names)\n", "x <- array(dim = c(num_images, 224, 224, 3))\n", "for (i in 1:num_images) {\n", " img_path <- paste(img_dir, image_names[i], sep = \"/\")\n", " img <- image_load(img_path, target_size = c(224, 224))\n", " x[i,,, ] <- image_to_array(img)\n", "}\n", "x <- imagenet_preprocess_input(x)" ] }, { "cell_type": "markdown", "id": "7624e1e4", "metadata": {}, "source": [ "We then load the trained network. The model has 50 layers, with a fair bit of complexity.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c7304261", "metadata": { "lines_to_next_cell": 2, "name": "chunk28" }, "outputs": [], "source": [ "model <- application_resnet50(weights = \"imagenet\")\n", "summary(model)" ] }, { "cell_type": "markdown", "id": "e659def9", "metadata": {}, "source": [ "Finally, we classify our six images, and return the top three class\n", "choices in terms of predicted probability for each." ] }, { "cell_type": "code", "execution_count": null, "id": "a8828a51", "metadata": { "lines_to_next_cell": 0, "name": "chunk29" }, "outputs": [], "source": [ "pred6 <- model %>% predict(x) %>%\n", " imagenet_decode_predictions(top = 3)\n", "names(pred6) <- image_names\n", "print(pred6)" ] }, { "cell_type": "markdown", "id": "a9973647", "metadata": {}, "source": [] }, { "cell_type": "markdown", "id": "5038d522", "metadata": {}, "source": [ "## IMDb Document Classification\n", "Now we perform document classification (Section 10.4) on the `IMDB` dataset, which is available as part of the `keras`\n", "package. We limit the dictionary size to the\n", "10,000 most frequently-used words and tokens." ] }, { "cell_type": "code", "execution_count": null, "id": "29d6b056", "metadata": { "name": "chunk30" }, "outputs": [], "source": [ "max_features <- 10000\n", "imdb <- dataset_imdb(num_words = max_features)\n", "c(c(x_train, y_train), c(x_test, y_test)) %<-% imdb" ] }, { "cell_type": "markdown", "id": "12d4ee2b", "metadata": {}, "source": [ "The third line is a shortcut for unpacking the list of lists.\n", "Each element of `x_train` is a vector of numbers between 0 and\n", "9999 (the document), referring to the words found in the dictionary.\n", "For example, the first training document is the positive review on\n", "page 419. The indices of the first 12 words are given below." ] }, { "cell_type": "code", "execution_count": null, "id": "1f568495", "metadata": { "name": "chunk31" }, "outputs": [], "source": [ "x_train[[1]][1:12]" ] }, { "cell_type": "markdown", "id": "902877c7", "metadata": {}, "source": [ "To see the words, we create a function, `decode_review()`, that provides a simple interface to the dictionary." ] }, { "cell_type": "code", "execution_count": null, "id": "baf4450a", "metadata": { "lines_to_next_cell": 2, "name": "chunk32" }, "outputs": [], "source": [ "word_index <- dataset_imdb_word_index()\n", "decode_review <- function(text, word_index) {\n", " word <- names(word_index)\n", " idx <- unlist(word_index, use.names = FALSE)\n", " word <- c(\"\", \"\", \"\", \"\", word)\n", " idx <- c(0:3, idx + 3)\n", " words <- word[match(text, idx, 2)]\n", " paste(words, collapse = \" \")\n", "}\n", "decode_review(x_train[[1]][1:12], word_index)" ] }, { "cell_type": "markdown", "id": "48bfa88a", "metadata": {}, "source": [ "Next we write a function to \"one-hot\" encode each document in a list\n", "of documents, and return a binary matrix in sparse-matrix format." ] }, { "cell_type": "code", "execution_count": null, "id": "e0bcbe04", "metadata": { "name": "chunk33" }, "outputs": [], "source": [ "library(Matrix)\n", "one_hot <- function(sequences, dimension) {\n", " seqlen <- sapply(sequences, length)\n", " n <- length(seqlen)\n", " rowind <- rep(1:n, seqlen)\n", " colind <- unlist(sequences)\n", " sparseMatrix(i = rowind, j = colind,\n", " dims = c(n, dimension))\n", "}" ] }, { "cell_type": "markdown", "id": "84f250d6", "metadata": {}, "source": [ "To construct the sparse matrix, one supplies just the entries that are\n", "nonzero. In the last line we call the function `sparseMatrix()` and supply the\n", "row indices corresponding to each document and the column indices\n", "corresponding to the words in each document, since we omit the\n", "values they are taken to be all ones.\n", "Words that appear more than once in any given document still get\n", "recorded as a one." ] }, { "cell_type": "code", "execution_count": null, "id": "4c922961", "metadata": { "name": "chunk34" }, "outputs": [], "source": [ "x_train_1h <- one_hot(x_train, 10000)\n", "x_test_1h <- one_hot(x_test, 10000)\n", "dim(x_train_1h)\n", "nnzero(x_train_1h) / (25000 * 10000)" ] }, { "cell_type": "markdown", "id": "9a21090f", "metadata": {}, "source": [ "Only 1.3\\% of the entries are nonzero, so this amounts to considerable\n", "savings in memory.\n", "We create a validation set of size 2,000, leaving 23,000 for training." ] }, { "cell_type": "code", "execution_count": null, "id": "d88a4673", "metadata": { "name": "chunk35" }, "outputs": [], "source": [ "set.seed(3)\n", "ival <- sample(seq(along = y_train), 2000)" ] }, { "cell_type": "markdown", "id": "070b63f5", "metadata": {}, "source": [ "First we fit a lasso logistic regression model using `glmnet()`\n", "on the training data, and evaluate its performance on the validation\n", "data. Finally, we plot the accuracy, `acclmv`, as a function of\n", "the shrinkage parameter, $\\lambda$. Similar expressions compute the\n", "performance on the test data, and were used to produce the left plot\n", "in Figure 10.11.\n", "The code takes advantage of the sparse-matrix format of\n", "`x_train_1h`, and runs in about 5 seconds; in the usual\n", "dense format it would take about 5 minutes." ] }, { "cell_type": "code", "execution_count": null, "id": "50073e95", "metadata": { "name": "chunk36" }, "outputs": [], "source": [ "library(glmnet)\n", "fitlm <- glmnet(x_train_1h[-ival, ], y_train[-ival],\n", " family = \"binomial\", standardize = FALSE)\n", "classlmv <- predict(fitlm, x_train_1h[ival, ]) > 0\n", "acclmv <- apply(classlmv, 2, accuracy, y_train[ival] > 0)" ] }, { "cell_type": "markdown", "id": "03719c63", "metadata": {}, "source": [ "We applied the `accuracy()` function that we wrote in Lab 10.9.2 to every column of the prediction matrix\n", "`classlmv`, and since this is a logical matrix of `TRUE/FALSE` values, we\n", "supply the second argument `truth` as a logical vector as well.\n", "\n", "Before making a plot, we adjust the plotting window." ] }, { "cell_type": "code", "execution_count": null, "id": "bf236b17", "metadata": { "lines_to_next_cell": 0, "name": "chunk37" }, "outputs": [], "source": [ "par(mar = c(4, 4, 4, 4), mfrow = c(1, 1))\n", "plot(-log(fitlm$lambda), acclmv)" ] }, { "cell_type": "markdown", "id": "00ab0cb5", "metadata": {}, "source": [ "\n" ] }, { "cell_type": "markdown", "id": "a43ca7e3", "metadata": {}, "source": [ "Next we fit a fully-connected neural network with two hidden layers,\n", "each with 16 units and ReLU activation.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "309863d8", "metadata": { "name": "chunk38" }, "outputs": [], "source": [ "model <- keras_model_sequential() %>%\n", " layer_dense(units = 16, activation = \"relu\",\n", " input_shape = c(10000)) %>%\n", " layer_dense(units = 16, activation = \"relu\") %>%\n", " layer_dense(units = 1, activation = \"sigmoid\")\n", "model %>% compile(optimizer = \"rmsprop\",\n", " loss = \"binary_crossentropy\", metrics = c(\"accuracy\"))\n", "history <- model %>% fit(x_train_1h[-ival, ], y_train[-ival],\n", " epochs = 20, batch_size = 512,\n", " validation_data = list(x_train_1h[ival, ], y_train[ival]))" ] }, { "cell_type": "markdown", "id": "7b145e30", "metadata": {}, "source": [ "The `history` object has a `metrics` component that\n", "records both the training and validation accuracy at each epoch.\n", "Figure 10.11 includes test accuracy at each epoch as well. To\n", "compute the test accuracy, we\n", "rerun the entire sequence above, replacing the last line with" ] }, { "cell_type": "code", "execution_count": null, "id": "18e63c04", "metadata": { "lines_to_next_cell": 2, "name": "chunk39" }, "outputs": [], "source": [ "history <- model %>% fit(\n", " x_train_1h[-ival, ], y_train[-ival], epochs = 20,\n", " batch_size = 512, validation_data = list(x_test_1h, y_test)\n", " )" ] }, { "cell_type": "markdown", "id": "1e988a9b", "metadata": {}, "source": [ "## Recurrent Neural Networks\n", "In this lab we fit the models illustrated in\n", "Section 10.5.\n", "### Sequential Models for Document Classification\n", "Here we fit a simple LSTM RNN for sentiment analysis with\n", "the `IMDb` movie-review data, as discussed in Section 10.5.1. We showed how to input the data in\n", "10.9.5, so we will not repeat that here.\n", "\n", "We first calculate the lengths of the documents." ] }, { "cell_type": "code", "execution_count": null, "id": "bc7d5157", "metadata": { "name": "chunk40" }, "outputs": [], "source": [ "wc <- sapply(x_train, length)\n", "median(wc)\n", "sum(wc <= 500) / length(wc)" ] }, { "cell_type": "markdown", "id": "217eae93", "metadata": {}, "source": [ "We see that over 91\\% of the documents have fewer than 500 words. Our\n", "RNN requires all the document sequences to have the same length. We hence\n", "restrict the document lengths to the last $L=500$ words, and pad the\n", "beginning of the\n", "shorter ones with blanks." ] }, { "cell_type": "code", "execution_count": null, "id": "3344bc1b", "metadata": { "name": "chunk41" }, "outputs": [], "source": [ "maxlen <- 500\n", "x_train <- pad_sequences(x_train, maxlen = maxlen)\n", "x_test <- pad_sequences(x_test, maxlen = maxlen)\n", "dim(x_train)\n", "dim(x_test)\n", "x_train[1, 490:500]" ] }, { "cell_type": "markdown", "id": "16e0fdf2", "metadata": {}, "source": [ "The last expression shows the last few words in the first document. At this stage, each of the 500 words in the document is represented using an integer\n", "corresponding to the location of that word in the 10,000-word dictionary.\n", "The first layer of the RNN is an embedding layer of size 32, which will be\n", "learned during training. This layer one-hot encodes each document\n", "as a matrix of dimension $500 \\times 10,000$, and then maps these\n", "$10,000$ dimensions down to $32$." ] }, { "cell_type": "code", "execution_count": null, "id": "676aa031", "metadata": { "name": "chunk42" }, "outputs": [], "source": [ "model <- keras_model_sequential() %>%\n", " layer_embedding(input_dim = 10000, output_dim = 32) %>%\n", " layer_lstm(units = 32) %>%\n", " layer_dense(units = 1, activation = \"sigmoid\")" ] }, { "cell_type": "markdown", "id": "019db081", "metadata": {}, "source": [ "The second layer is an LSTM with 32 units, and the output\n", "layer is a single sigmoid for the binary classification task.\n", "\n", "The rest is now similar to other networks we have fit. We\n", "track the test performance as the network is fit, and see that it attains 87\\% accuracy." ] }, { "cell_type": "code", "execution_count": null, "id": "c6fbcd9e", "metadata": { "lines_to_next_cell": 2, "name": "chunk43" }, "outputs": [], "source": [ "model %>% compile(optimizer = \"rmsprop\",\n", " loss = \"binary_crossentropy\", metrics = c(\"acc\"))\n", "history <- model %>% fit(x_train, y_train, epochs = 10,\n", " batch_size = 128, validation_data = list(x_test, y_test))\n", "plot(history)\n", "predy <- predict(model, x_test) > 0.5\n", "mean(abs(y_test == as.numeric(predy)))" ] }, { "cell_type": "markdown", "id": "d2075145", "metadata": {}, "source": [ "### Time Series Prediction\n", "We now show how to fit the models in Section 10.5.2\n", "for time series prediction.\n", "We first set up the data, and standardize each of the variables." ] }, { "cell_type": "code", "execution_count": null, "id": "3c8b0938", "metadata": { "name": "chunk44" }, "outputs": [], "source": [ "library(ISLR2)\n", "xdata <- data.matrix(\n", " NYSE[, c(\"DJ_return\", \"log_volume\",\"log_volatility\")]\n", " )\n", "istrain <- NYSE[, \"train\"]\n", "xdata <- scale(xdata)" ] }, { "cell_type": "markdown", "id": "ede6b449", "metadata": {}, "source": [ "The variable `istrain` contains a `TRUE` for each year that is in the training set, and a `FALSE` for each year\n", " in the test set.\n", "\n", "We first write functions to create lagged versions of the three time series. We start with a function that takes as input a data\n", "matrix and a lag $L$, and returns a lagged version of the matrix. It\n", "simply inserts $L$ rows of `NA` at the top, and truncates the\n", "bottom." ] }, { "cell_type": "code", "execution_count": null, "id": "a09feef2", "metadata": { "name": "chunk45" }, "outputs": [], "source": [ "lagm <- function(x, k = 1) {\n", " n <- nrow(x)\n", " pad <- matrix(NA, k, ncol(x))\n", " rbind(pad, x[1:(n - k), ])\n", "}" ] }, { "cell_type": "markdown", "id": "f47eda44", "metadata": {}, "source": [ "We now use this function to create a data frame with all the required\n", "lags, as well as the response variable." ] }, { "cell_type": "code", "execution_count": null, "id": "00a55161", "metadata": { "name": "chunk46" }, "outputs": [], "source": [ "arframe <- data.frame(log_volume = xdata[, \"log_volume\"],\n", " L1 = lagm(xdata, 1), L2 = lagm(xdata, 2),\n", " L3 = lagm(xdata, 3), L4 = lagm(xdata, 4),\n", " L5 = lagm(xdata, 5)\n", " )" ] }, { "cell_type": "markdown", "id": "7bf8b0d3", "metadata": {}, "source": [ "If we look at the first five rows of this frame, we will see some\n", "missing values in the lagged variables (due to the construction above). We remove these rows, and adjust `istrain`\n", "accordingly." ] }, { "cell_type": "code", "execution_count": null, "id": "ee49f4ef", "metadata": { "name": "chunk47" }, "outputs": [], "source": [ "arframe <- arframe[-(1:5), ]\n", "istrain <- istrain[-(1:5)]" ] }, { "cell_type": "markdown", "id": "26430d53", "metadata": {}, "source": [ "We now fit the linear AR model to the training data using `lm()`, and predict on the\n", "test data." ] }, { "cell_type": "code", "execution_count": null, "id": "6801cf43", "metadata": { "name": "chunk48" }, "outputs": [], "source": [ "arfit <- lm(log_volume ~ ., data = arframe[istrain, ])\n", "arpred <- predict(arfit, arframe[!istrain, ])\n", "V0 <- var(arframe[!istrain, \"log_volume\"])\n", "1 - mean((arpred - arframe[!istrain, \"log_volume\"])^2) / V0" ] }, { "cell_type": "markdown", "id": "a8a5b310", "metadata": {}, "source": [ "The last two lines compute the $R^2$ on the test data, as defined in (3.17).\n", "\n", "We refit this model, including the factor variable `day_of_week`." ] }, { "cell_type": "code", "execution_count": null, "id": "45e8eb7a", "metadata": { "lines_to_next_cell": 2, "name": "chunk49" }, "outputs": [], "source": [ "arframed <- data.frame(day = NYSE[-(1:5), \"day_of_week\"], arframe)\n", "arfitd <- lm(log_volume ~ ., data = arframed[istrain, ])\n", "arpredd <- predict(arfitd, arframed[!istrain, ])\n", "1 - mean((arpredd - arframe[!istrain, \"log_volume\"])^2) / V0" ] }, { "cell_type": "markdown", "id": "8e9ffd1a", "metadata": {}, "source": [ "To fit the RNN, we need to reshape these data, since it expects a\n", "sequence of $L=5$ feature vectors $X={X_\\ell}_1^L$ for each observation, as in (10.20) on\n", "page 428. These are lagged versions of the\n", "time series going back $L$ time points." ] }, { "cell_type": "code", "execution_count": null, "id": "1e3833fc", "metadata": { "name": "chunk50" }, "outputs": [], "source": [ "n <- nrow(arframe)\n", "xrnn <- data.matrix(arframe[, -1])\n", "xrnn <- array(xrnn, c(n, 3, 5))\n", "xrnn <- xrnn[,, 5:1]\n", "xrnn <- aperm(xrnn, c(1, 3, 2))\n", "dim(xrnn)" ] }, { "cell_type": "markdown", "id": "3ea26e7e", "metadata": {}, "source": [ "We have done this in four steps. The first simply extracts the\n", "$n\\times 15$ matrix of lagged versions of the three predictor\n", "variables from `arframe`. The second converts this matrix to an\n", "$n\\times 3\\times 5$ array. We can do this by simply changing the\n", "dimension attribute, since the new array is filled column wise. The\n", "third step reverses the order of lagged variables, so that index $1$\n", "is furthest back in time, and index $5$ closest. The\n", "final step rearranges the coordinates of the array (like a partial\n", "transpose) into the format that the RNN module in `keras`\n", "expects.\n", "\n", "Now we are ready to proceed with the RNN, which uses 12 hidden units." ] }, { "cell_type": "code", "execution_count": null, "id": "b84cd646", "metadata": { "name": "chunk51" }, "outputs": [], "source": [ "model <- keras_model_sequential() %>%\n", " layer_simple_rnn(units = 12,\n", " input_shape = list(5, 3),\n", " dropout = 0.1, recurrent_dropout = 0.1) %>%\n", " layer_dense(units = 1)\n", "model %>% compile(optimizer = optimizer_rmsprop(),\n", " loss = \"mse\")" ] }, { "cell_type": "markdown", "id": "31f4128b", "metadata": {}, "source": [ "We specify two forms of dropout for the units feeding into the hidden\n", " layer. The first is for\n", "the input sequence feeding into this layer, and the second is for the\n", "previous hidden units feeding into the layer.\n", "The output layer has a single unit for the response." ] }, { "cell_type": "markdown", "id": "4fbb2b80", "metadata": {}, "source": [ "We fit the model in a similar fashion to previous networks. We\n", "supply the `fit` function with test data as validation data, so that when\n", "we monitor its progress and plot the history function we can see the\n", "progress on the test data. Of course we should not use this as a basis for\n", "early stopping, since then the test performance would be biased." ] }, { "cell_type": "code", "execution_count": null, "id": "f92366f1", "metadata": { "name": "chunk52" }, "outputs": [], "source": [ "history <- model %>% fit(\n", " xrnn[istrain,, ], arframe[istrain, \"log_volume\"],\n", " batch_size = 64, epochs = 200,\n", " validation_data =\n", " list(xrnn[!istrain,, ], arframe[!istrain, \"log_volume\"])\n", " )\n", "kpred <- predict(model, xrnn[!istrain,, ])\n", "1 - mean((kpred - arframe[!istrain, \"log_volume\"])^2) / V0" ] }, { "cell_type": "markdown", "id": "495255d6", "metadata": {}, "source": [ " This model takes about one minute to train.\n", "\n", "We could replace the `keras_model_sequential()` command above with the following command:" ] }, { "cell_type": "code", "execution_count": null, "id": "9ca991e9", "metadata": { "name": "chunk53" }, "outputs": [], "source": [ "model <- keras_model_sequential() %>%\n", " layer_flatten(input_shape = c(5, 3)) %>%\n", " layer_dense(units = 1)" ] }, { "cell_type": "markdown", "id": "3130c0c8", "metadata": {}, "source": [ "Here, `layer_flatten()` simply takes the input sequence and\n", "turns it into a long vector of predictors. This results in a linear AR model.\n", "To fit a nonlinear AR model, we could add in a hidden layer.\n", "\n", "However, since we already have the matrix of lagged variables from the AR\n", "model that we fit earlier using the `lm()` command, we can actually fit a nonlinear AR model without needing to perform flattening.\n", "We extract the model matrix `x` from `arframed`, which\n", "includes the `day_of_week` variable." ] }, { "cell_type": "code", "execution_count": null, "id": "b8011901", "metadata": { "name": "chunk54" }, "outputs": [], "source": [ "x <- model.matrix(log_volume ~ . - 1, data = arframed)\n", "colnames(x)" ] }, { "cell_type": "markdown", "id": "ba1fdb4a", "metadata": {}, "source": [ "The `-1` in the formula avoids the creation of a column of ones for\n", "the intercept. The variable `day_of_week` is a five-level\n", "factor (there are five trading days), and the\n", " `-1` results in five rather than four dummy variables.\n", "\n", "The rest of the steps to fit a nonlinear AR model should by now be familiar." ] }, { "cell_type": "code", "execution_count": null, "id": "237c7cc9", "metadata": { "lines_to_next_cell": 0, "name": "chunk55" }, "outputs": [], "source": [ "arnnd <- keras_model_sequential() %>%\n", " layer_dense(units = 32, activation = 'relu',\n", " input_shape = ncol(x)) %>%\n", " layer_dropout(rate = 0.5) %>%\n", " layer_dense(units = 1)\n", "arnnd %>% compile(loss = \"mse\",\n", " optimizer = optimizer_rmsprop())\n", "history <- arnnd %>% fit(\n", " x[istrain, ], arframe[istrain, \"log_volume\"], epochs = 100, \n", " batch_size = 32, validation_data =\n", " list(x[!istrain, ], arframe[!istrain, \"log_volume\"])\n", " )\n", "plot(history)\n", "npred <- predict(arnnd, x[!istrain, ])\n", "1 - mean((arframe[!istrain, \"log_volume\"] - npred)^2) / V0" ] }, { "cell_type": "markdown", "id": "e7de3148", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "name,tags,-all" }, "kernelspec": { "display_name": "R", "language": "R", "name": "ir" } }, "nbformat": 4, "nbformat_minor": 5 }