# Machine learning for econometricians # Course for the ECB # by Ewen Gallic # December 2021 # Hands-on Session 3 # Deep Learning library(tidyverse) # 1 Neural Networks ------------------------------------------------------- sigmoid_f <- function(z) 1/(1+exp(-z)) hyperbolic_tan_f <- function(z) (exp(2*z) - 1)/(exp(2*z)+1) hard_threshold_f <- function(z, beta) as.numeric(z >= beta) relu_f <- function(z) pmax(0, z) x <- seq(-2, 2, length.out = 1000) df <- tibble( x = x, identity = x, sigmoid = sigmoid_f(x), tanh = hyperbolic_tan_f(x), threshold = hard_threshold_f(x, beta = 1), ReLU = relu_f(x)) df <- df %>% gather(type, value, -x) %>% mutate(type = factor( type, levels = c("identity", "sigmoid", "tanh", "threshold", "ReLU"), labels = c("Identity", "Sigmoid", "Tanh", "Threshold (beta=1)", "ReLU"))) ggplot(data = df, aes(x = x, y = value, colour = type)) + geom_line(size = 1.1) + labs(x = "z", y = "f(z)") + scale_colour_manual( NULL, values = c("Identity" = "#648FFF", "Sigmoid" = "#785EF0", "Tanh" = "#DC267F", "Threshold (beta=1)" = "#FE6100", "ReLU" = "#FFB000")) ## install.packages("tensorflow") ## install.packages("keras") library(keras) # First time: ## tensorflow::install_tensorflow() ## tensorflow::tf_config() ## install_keras() # If it does not work: ## unlink(reticulate::miniconda_path(), recursive = TRUE) ## reticulate::install_miniconda(path = reticulate::miniconda_path(), ## update = TRUE, force = FALSE) ## keras::install_keras() # Run once ## install_keras(tensorflow = "gpu") # load the MNIST dataset mnist <- dataset_mnist() x_train <- mnist$train$x dim(x_train) x_test <- mnist$test$x dim(x_test) y_train <- mnist$train$y head(y_train) y_test <- mnist$test$y head(y_test) # Other way to write it: # c(c(x_train, y_train), c(x_test, y_test)) %<-% mnist # dim(x_train) # dim(x_test) # head(y_train) # head(y_test) # There is approximatively 10% of each of the # 10 different classes (0, 1, 2, …, 9): round(100*prop.table(table(y_train)),1) # First training example y_train[1] # Let us plot an example rotate <- function(x) t(apply(x, 2, rev)) convert_example_tibble <- function(example){ rotate(example) %>% reshape2::melt(c("xx", "yy"), value.name = "grey_scale") %>% as_tibble() } plot_example <- function(example){ ggplot(data = convert_example_tibble(example), mapping = aes(x = xx, y = yy, fill = grey_scale))+ geom_tile() + coord_equal() + scale_fill_gradient("Grey Scale", low = "white", high = "black") } plot_example(x_train[1,,]) # The first 10 observations from MNIST training set : map_df(1:10, ~convert_example_tibble(x_train[.,,]), .id = "ind") %>% mutate(ind = as.numeric(ind)) %>% left_join(tibble(y = y_train[1:10], ind = 1:10), by = "ind") %>% mutate(lab = str_c("ex. ", ind, ": ", y)) %>% ggplot(data = ., mapping = aes(x = xx, y = yy, fill = grey_scale))+ geom_tile() + coord_equal() + scale_fill_gradient("Grey Scale", low = "white", high = "black") + facet_wrap(~lab, ncol = 5) + theme(axis.text = element_blank(), axis.ticks = element_blank()) + labs(x = NULL, y = NULL) # 1.2.1 Practice With Keras: classifier ----------------------------------- # Reshaping data to get a matrix x_train <- array_reshape(x_train, dim = c(nrow(x_train), 784)) x_test <- array_reshape(x_test, dim = c(nrow(x_test), 784)) # Scaling x_train <- x_train/255 # One-hot encoding of the response variable y_train <- to_categorical(y_train) # Scaling on the test set as well x_test <- x_test/255 y_test <- to_categorical(y_test) # Splitting the training in two parts, the first one to train # and the second to observe the predictive abilities of the model n_valid <- round(.2*nrow(x_train)) ind_partial_valid <- sample(1:nrow(x_train), size = n_valid, replace = FALSE) x_train_partial <- x_train[-ind_partial_valid, ] x_validation <- x_train[ind_partial_valid, ] y_train_partial <- y_train[-ind_partial_valid,] y_validation <- y_train[ind_partial_valid,] dim(x_train_partial) dim(x_validation) dim(y_train_partial) dim(y_validation) # sequential model: model_architecture <- keras_model_sequential(name = "MNIST_Model") %>% # First hidden layer layer_dense(units = 256, activation = "relu", input_shape = 784) %>% # Second hidden layer layer_dense(units = 128, activation = "relu") %>% # Output layer layer_dense(units = 10, activation = 'softmax') model_architecture # Loss function and metrics computed during training model_architecture %>% compile( loss = "categorical_crossentropy", optimizer = "rmsprop", metrics = c("accuracy") ) # Training the model # Computing the metrics on both training and validation sets mnist_model_hist <- model_architecture %>% fit(x_train_partial, y_train_partial, epochs = 20, batch_size = 128, validation_data = list(x_validation, y_validation)) # Saving the results if(!dir.exists("deep_learning_models")) dir.create("deep_learning_models") model_architecture %>% save_model_hdf5("deep_learning_models/mnist_model_1.h5") save(mnist_model_hist, file = "deep_learning_models/mnist_model_1_hist.rda") # To load the results: # model_architecture <- load_model_hdf5("deep_learning_models/mnist_model_1.h5") # load("deep_learning_models/mnist_model_1_hist.rda") # Metrics computed during training mnist_model_hist$metrics # Accuracy and loss after each epoch, both on the training and validation sets. plot(mnist_model_hist) # Same kind of graph, custom made (just for fun): as_tibble(mnist_model_hist$metrics) %>% mutate(epoch = row_number()) %>% pivot_longer(cols = -epoch) %>% mutate(data = ifelse(str_detect(name, "^val_"), "validation", "training"), name = str_remove(name, "^val_")) %>% ggplot(data = ., mapping = aes(x = epoch, y = value, colour = data)) + geom_line() + geom_point() + scale_colour_manual("Sample", values = c("validation" = "#E69F00", "training" = "#0072B2")) + facet_wrap(~name, scales = "free_y", ncol = 1) # performances of the model on the test set performances_test <- model_architecture %>% evaluate(x_test, y_test) performances_test # pedictions (on unseen data, here): mnist_pred_test <- model_architecture %>% predict(x_test) # Prdicted probabilities mnist_pred_test[1,] # Observed value y_test[1, ] # make sure that the sum of the predicted probabilities, # thanks to the softmax function used in the output layer, sum up to 1. sum(mnist_pred_test[1,]) # Let us have a look at the predictive capacities for each categories. mnist_pred_test_class <- as_tibble(mnist_pred_test, .name_repair = "unique") %>% mutate(pred_class = max.col(.), pred_class = seq(0,9)[pred_class]) %>% magrittr::extract2("pred_class") mnist_test_class <- y_test %>% as_tibble(.name_repair = "unique") %>% mutate(y_class = max.col(.), y_class = seq(0,9)[y_class]) %>% magrittr::extract2("y_class") # Confusion table table(observed = mnist_test_class, predicted = mnist_pred_test_class) # A graph showing the percentage of correct predictions can be made. # To that end, let us obtain the percentage of correct predictions for each class. table_pct_pred <- tibble( observed = factor(mnist_test_class, levels = 0:9), predicted = factor(mnist_pred_test_class, levels = 0:9) ) %>% group_by(observed) %>% mutate(n_k = sum(n())) %>% group_by(observed, predicted, n_k) %>% count() %>% mutate(pred_pct = n / n_k) %>% ungroup() %>% select(observed, predicted, pred_pct) table_pct_pred # If we use this table directly, there will be missing cases. # Let us complete it with 0 values when there are 0 predictions # made for a given class when the true class is of class k. classes_mnist <- levels(factor(mnist_test_class, levels = 0:9)) df_plot_pred_mnist <- expand_grid( observed = classes_mnist, predicted = classes_mnist) %>% left_join( table_pct_pred ) %>% mutate(pred_pct = replace_na(pred_pct, 0)) ggplot(data = df_plot_pred_mnist, mapping = aes(x = predicted, y = observed, fill = pred_pct, label = round(pred_pct, 2))) + geom_tile() + geom_text(color = "white") + scale_fill_gradient("% of\n correct pred.", low = "#1E88E5", high = "#D81B60") + labs(x = "Predicted", y = "Observed", title = "Percentage of predicted classes for each observed class.") + theme(panel.background = element_blank(), plot.title.position = "plot") # Another way of looking at the predictive capacities consists in # looking at the percentage of correctly predicted observation per class. table_pct_pred_2 <- tibble( observed = factor(mnist_test_class, levels = 0:9), predicted = factor(mnist_pred_test_class, levels = 0:9) ) %>% mutate(correct = observed == predicted) %>% group_by(observed) %>% summarise(pct_correct = sum(correct) / n()) %>% mutate(pct_incorrect = 1-pct_correct) %>% pivot_longer(cols = -observed) %>% mutate(name = factor(name, levels = c("pct_correct", "pct_incorrect"), labels = c("Correctly predicted", "Incorrect predicted"))) ggplot(data = table_pct_pred_2, mapping = aes(x = name, y = observed, fill = value, label = round(value, 2))) + geom_tile() + geom_text(color = "white") + scale_fill_gradient("% of\n correct pred.", low = "#1E88E5", high = "#D81B60") + labs(x = "Accuracy of prediction", y = "Observed", title = "% of correctly or incorrectly predicted obs by observed class.") + theme(panel.background = element_blank(), plot.title.position = "plot") # Increasing the number of units model_architecture_2 <- keras_model_sequential(name = "MNIST_Model_2") %>% # First hidden layer layer_dense(units = 256, activation = "relu", input_shape = 784) %>% # Second hidden layer layer_dense(units = 2, activation = "relu") %>% # Output layer layer_dense(units = 10, activation = 'softmax') model_architecture_2 %>% compile( loss = "categorical_crossentropy", optimizer = "rmsprop", metrics = c("accuracy") ) mnist_model_2_hist <- model_architecture_2 %>% fit(x_train_partial, y_train_partial, epochs = 20, batch_size = 128, validation_data = list(x_validation, y_validation)) mnist_model_2_hist # The accuracy of the model dropped. # This is because the information was too compressed # when going through the second layer. as_tibble(mnist_model_hist$metrics) %>% mutate(model = "Initial model", epoch = row_number()) %>% bind_rows(as_tibble(mnist_model_2_hist$metrics) %>% mutate(model = "Second model", epoch = row_number())) %>% pivot_longer(cols = -c(epoch, model)) %>% mutate(data = ifelse(str_detect(name, "^val_"), "validation", "training"), name = str_remove(name, "^val_")) %>% ggplot(data = ., mapping = aes(x = epoch, y = value, colour = data, linetype = model)) + geom_line() + geom_point() + scale_colour_manual("Sample", values = c("validation" = "#E69F00", "training" = "#0072B2")) + facet_wrap(~name, scales = "free_y", ncol = 1) + labs(y = NULL, title = "Accuracy and loss after each epoch, custom made graph.") + theme(plot.title.position = "plot") # Saving the results model_architecture_2 %>% save_model_hdf5("deep_learning_models/mnist_model.h5") # Loading the results # model_architecture_2 <- # load_model_hdf5("deep_learning_models/mnist_model.h5") # To destroy the current TensorFlow graph and create a new one, # the function k_clear_session() can be used. # This is useful to avoid clutter from previous models or layers. k_clear_session() # 1.2.2 Practice With Keras: Regression ----------------------------------- # Load data load(url("http://egallic.fr/Enseignement/ML/ECB/data/bike.rda")) # transform categorical variables to dummy variables # Example with toy data: library(fastDummies) tibble(x = c("A", "B", "C"), y = factor(c("mon", "tue", "wed")), z = 1:3) %>% # Creating dummies for chr and fct dummy_cols( remove_first_dummy = TRUE, remove_selected_columns = TRUE) # Let us do it with the bike dataset: x_data <- bike %>% select(-date, -rented_bike_count, -y_binary) %>% dummy_cols( remove_first_dummy = TRUE, remove_selected_columns = TRUE) # The variable to predict (the target): y_data <- bike$rented_bike_count # training and test sets n_train <- round(.8*nrow(x_data)) set.seed(123) ind_train <- sample(1:nrow(x_data), size = n_train, replace = FALSE) # Train x_train <- x_data[ind_train, ] y_train <- y_data[ind_train] dim(x_train) length(y_train) # Test x_test <- x_data[-ind_train, ] y_test <- y_data[-ind_train] dim(x_test) length(y_test) # Preparing scaling: mean and standard-deviation # computed on training set, for all variables mean_train <- apply(x_test, 2, mean) std_train <- apply(x_test, 2, sd) # Scaling x_train <- scale(x_train, center = mean_train, scale = std_train) x_test <- scale(x_test, center = mean_train, scale = std_train) # function that defines the model depending on the train data model_structure_bike <- function(train_data) { model <- keras_model_sequential() %>% # First hidden layer layer_dense(units = 128, activation = "relu", input_shape = dim(train_data)[[2]]) %>% # Second hidden layer layer_dense(units = 128, activation = "relu") %>% # Third hidden layer layer_dense(units = 64, activation = "relu") %>% # Output layer layer_dense(units = 1) model %>% compile( optimizer = "rmsprop", loss = "mse", metrics = c("mae") ) } # Number of folds for CV nb_folds <- 5 set.seed(123) folds <- sample(rep(1:nb_folds, length=nrow(x_train))) # Number of epochs and batch size num_epochs <- 10 batch_size <- 32 # Init object that will contain the metrics loss_and_metrics <- NULL for(k in 1:nb_folds){ # cat(str_c("\nFold number: ", k, "\n")) ind_current <- folds != k # Train set x_train_partial <- x_train[which(ind_current),] y_train_partial <- y_train[which(ind_current)] # Validation set x_train_valid <- x_train[which(!ind_current),] y_train_valid <- y_train[which(!ind_current)] # Building the model model <- model_structure_bike(x_train_partial) # Training the model on the training set model %>% fit(x_train_partial, y_train_partial, epochs = num_epochs, batch_size = batch_size, verbose = 0) # Evaluating its performances on test data results_current <- model %>% evaluate(x_train_valid, y_train_valid, verbose = 0) # Store loss and metrics loss_and_metrics <- loss_and_metrics %>% bind_rows( as_tibble(t(results_current)) %>% mutate(fold = k) ) } # Saving the model model %>% save_model_hdf5("deep_learning_models/bike_model.h5") save(loss_and_metrics, file = "deep_learning_models/bike_loss_and_metrics_neural_net.rda") # model <- load_model_hdf5("deep_learning_models/bike_model.h5") # load("deep_learning_models/bike_loss_and_metrics_neural_net.rda") loss_and_metrics # average of the the Mean Average Error over the folds: mean(loss_and_metrics$mae) # compared to the distribution of the target variable: summary(y_train) # Let us train the model over more epochs to see if it yields better results. num_epochs <- 800 batch_size <- 32 loss_and_metrics_2 <- vector(mode = "list", length = nb_folds) library(tictoc) results <- vector(mode = "list", length = nb_folds) tic() for(k in 1:nb_folds){ cat(str_c("\nFold number: ", k, "\n")) ind_current <- folds != k # Train set x_train_partial <- x_train[which(ind_current),] y_train_partial <- y_train[which(ind_current)] # Validation set x_train_valid <- x_train[which(!ind_current),] y_train_valid <- y_train[which(!ind_current)] # Building the model model <- model_structure_bike(x_train_partial) # Training the model on the training set history <- model %>% fit(x_train_partial, y_train_partial, validation_data = list(x_train_valid, y_train_valid), epochs = num_epochs, batch_size = batch_size, verbose = 0) # Evaluating its performances on test data results_current <- model %>% evaluate(x_train_valid, y_train_valid, verbose = 0) results[[k]] <- results_current # Store loss and metrics loss_and_metrics_2[[k]] <- history } toc() model %>% save_model_hdf5("deep_learning_models/bike_model_2.h5") save(loss_and_metrics_2, file = "deep_learning_models/bike_loss_and_metrics_neural_net_2.rda") # model <- load_model_hdf5("deep_learning_models/bike_model_2.h5") # load("deep_learning_models/bike_loss_and_metrics_neural_net_2.rda") # Let us compute, for each epoch, the average of the 5-folds MAE both on # the train and on the validation samples: mae_mean_folds <- map(loss_and_metrics_2, "metrics") %>% map("val_mae") %>% map_df(~as_tibble(.) %>% mutate(epoch = row_number()), .id = "fold") %>% mutate(sample = "validation") %>% bind_rows( map(loss_and_metrics_2, "metrics") %>% map("mae") %>% map_df(~as_tibble(.) %>% mutate(epoch = row_number()), .id = "fold") %>% mutate(sample = "train") ) %>% group_by(epoch, sample) %>% summarise(mae = mean(value), .groups = "drop") mae_mean_folds ggplot(data = mae_mean_folds, mapping = aes(x = epoch, y = mae, colour = sample)) + geom_line() + scale_colour_manual("Sample", values = c("validation" = "#E69F00", "train" = "#0072B2")) + labs(x = "Epoch", y = "MAE (validation sample)", title = "MAE computed on the validation samples, over the epochs.") + theme(plot.title.position = "plot") # Let us say our final model is the following one: num_epochs <- 400 batch_size <- 32 model <- model_structure_bike(x_train) model %>% fit(x_train, y_train, epochs = num_epochs, batch_size = batch_size, verbose = 0) model %>% save_model_hdf5("deep_learning_models/bike_model_3.h5") # model <- load_model_hdf5("deep_learning_models/bike_model_3.h5") # predictive abilities of the model on unseen data: results_test <- model %>% evaluate(x_test, y_test) results_test # 1.2.3 Regularisation techniques ----------------------------------------- x_data <- bike %>% select(-date, -rented_bike_count, -y_binary) %>% dummy_cols( remove_first_dummy = TRUE, remove_selected_columns = TRUE) # 1 if "Low", 0 otherwise # y_data <- to_categorical(ifelse(bike$y_binary == "Low", yes = 1, no = 0)) y_data <- ifelse(bike$y_binary == "Low", yes = 1, no = 0) # Train and test sets n_train <- round(.8*nrow(x_data)) set.seed(123) ind_train <- sample(1:nrow(x_data), size = n_train, replace = FALSE) x_train <- x_data[ind_train, ] y_train <- y_data[ind_train] dim(x_train) length(y_train) x_test <- x_data[-ind_train, ] y_test <- y_data[-ind_train] dim(x_test) length(y_test) # Scaling mean_train <- apply(x_test, 2, mean) std_train <- apply(x_test, 2, sd) x_train <- scale(x_train, center = mean_train, scale = std_train) x_test <- scale(x_test, center = mean_train, scale = std_train) # K-fold Cross-validation nb_folds <- 3 set.seed(123) folds <- sample(rep(1:nb_folds, length=nrow(x_train))) # architecture of the model model_structure_bike_l2 <- function(train_data, l=0.01) { model <- keras_model_sequential() %>% # First hidden layer layer_dense(units = 32, activation = "relu", input_shape = dim(train_data)[[2]], kernel_regularizer = regularizer_l2(l)) %>% # Second hidden layer layer_dense(units = 64, activation = "relu", kernel_regularizer = regularizer_l2(l)) %>% # Third hidden layer layer_dense(units = 32, activation = "relu", kernel_regularizer = regularizer_l2(l)) %>% # Output layer layer_dense(units = 1, activation = "sigmoid") model %>% compile( optimizer = "rmsprop", loss = "binary_crossentropy", metrics = c("accuracy", metric_true_positives(name = "tp"), metric_true_negatives(name = "tn"), metric_false_positives(name = "fp"), metric_false_negatives(name = "fn")) ) } # Let us train a model on 300 epochs, with a batch size of 512 observations. num_epochs <- 300 batch_size <- 512 metrics_regul <- vector(mode = "list", length = nb_folds) for(k in 1:nb_folds){ cat(str_c("\nFold number: ", k, "\n")) ind_current <- folds != k # Train set x_train_partial <- x_train[which(ind_current),] y_train_partial <- y_train[which(ind_current)] # Validation set x_train_valid <- x_train[which(!ind_current),] y_train_valid <- y_train[which(!ind_current)] # Building the model model_with_regul <- model_structure_bike_l2(x_train_partial, l = 0.01) # Training the model on the training set history <- model_with_regul %>% fit(x_train_partial, y_train_partial, validation_data = list(x_train_valid, y_train_valid), epochs = num_epochs, batch_size = batch_size, verbose = 0) # Store loss and metrics metrics_regul[[k]] <- history } save(metrics_regul, file = "deep_learning_models/metrics_regul.rda") # load("deep_learning_models/metrics_regul.rda") # For comparison, a model without regularisation can be trained. model_structure_bike <- function(train_data) { model <- keras_model_sequential() %>% # First hidden layer layer_dense(units = 32, activation = "relu", input_shape = dim(train_data)[[2]]) %>% # Second hidden layer layer_dense(units = 64, activation = "relu") %>% # Third hidden layer layer_dense(units = 32, activation = "relu") %>% # Output layer layer_dense(units = 1, activation = "sigmoid") model %>% compile( optimizer = "rmsprop", loss = "binary_crossentropy", metrics = c("accuracy", metric_true_positives(name = "tp"), metric_true_negatives(name = "tn"), metric_false_positives(name = "fp"), metric_false_negatives(name = "fn")) ) } # Then, the predictive abilities of the model can be assessed # on the validation sets during the 3-fold cross-validation. num_epochs <- 500 batch_size <- 512 metrics_without_regul <- vector(mode = "list", length = nb_folds) for(k in 1:nb_folds){ cat(str_c("\nFold number: ", k, "\n")) ind_current <- folds != k # Train set x_train_partial <- x_train[which(ind_current),] y_train_partial <- y_train[which(ind_current)] # Validation set x_train_valid <- x_train[which(!ind_current),] y_train_valid <- y_train[which(!ind_current)] # Building the model model <- model_structure_bike(x_train_partial) # Training the model on the training set history <- model %>% fit(x_train_partial, y_train_partial, validation_data = list(x_train_valid, y_train_valid), epochs = num_epochs, batch_size = batch_size, verbose = 0) # Store loss and metrics metrics_without_regul[[k]] <- history } save(metrics_without_regul, file = "deep_learning_models/metrics_without_regul.rda") # load("deep_learning_models/metrics_without_regul.rda") # Let us define a function to extract the metric values from # the history stored during the k-fold cross-validation process. # This function will return the average over the k folds for each epoch, # in a tibble. #' @param k_fold_history list with the training history during the k-fold CV #' @param metric (string) name of the metric to extract from the history history_to_tibble <- function(k_fold_history, metric){ map(k_fold_history, "metrics") %>% # Metrics on the validation sets map(str_c("val_",metric)) %>% map_df(~as_tibble(.) %>% mutate(epoch = row_number()), .id = "fold") %>% mutate(sample = "validation") %>% bind_rows( map(k_fold_history, "metrics") %>% map(metric) %>% map_df(~as_tibble(.) %>% mutate(epoch = row_number()), .id = "fold") %>% mutate(sample = "train") ) %>% group_by(epoch, sample) %>% summarise(!!sym(metric) := mean(value), .groups = "drop") %>% mutate(sample = factor(sample, levels = c("validation", "train"))) } df_metrics_regul <- history_to_tibble(metrics_without_regul, "accuracy") %>% mutate(type = "Without regularisation") %>% bind_rows( history_to_tibble(metrics_regul, "accuracy") %>% mutate(type = "L2 regularisation (l=0.01)") ) df_metrics_regul # Then we can look at the results: ggplot(data = df_metrics_regul, mapping = aes(x = epoch, y = accuracy, linetype = sample, colour = type)) + geom_line() + scale_colour_manual( "Regularisation", values = c("Without regularisation" = "#000000", "L2 regularisation (l=0.01)" = "#E69F00")) + scale_linetype_discrete("Sample") + labs(x = "Epoch", y = "Accuracy", title = "Accuracy with or without regularisation.") + theme(legend.position = "bottom", legend.box="vertical", plot.title.position = "plot") # 1.2.3.2 Dropout --------------------------------------------------------- # Let us simulate some data. set.seed(123) x <- rnorm(4) W_1 <- matrix(round(rnorm(20), 2), ncol = 5) colnames(W_1) <- str_c("L1 Unit ", 1:5) W_1 W_2 <- matrix(round(rnorm(20), 2), ncol = 4) colnames(W_2) <- str_c("L2 Unit ", 1:4) W_2 beta <- c(0.23, -0.10, 0.13, 0.02) bias_0 <- rnorm(5) bias_1 <- rnorm(4) bias_2 <- rnorm(1) bias_0 bias_1 bias_2 # Each unit of the first hidden layer makes a linear combination # of the inputs, using the current weights. linear_comb_layer_1 <- x %*% W_1 + bias_0 linear_comb_layer_1 # Let us say we use a ReLU function for the first hidden layer A_k_layer_1 <- relu_f(linear_comb_layer_1) A_k_layer_1 # The linear combination made in the second layer, by each unit: linear_comb_layer_2 <- A_k_layer_1 %*% W_2 + bias_1 linear_comb_layer_2 # The output value of each unit is then, using a ReLU as well as the # activation function A_k_layer_2 <- relu_f(linear_comb_layer_2) A_k_layer_2 # Assuming the identity function as the activation function # in the output layer, the prediction is: pred <- A_k_layer_2 %*% beta + bias_2 pred # Let us apply inverted dropout regularisation # dropout rate of .2 keep_rate <- .8 # Randomly selecting the kept units kept_units_L1 <- matrix(rbinom(n=length(W_1), size = 1, prob = keep_rate), ncol = ncol(W_1)) kept_units_L1 # In the second hidden layer as well kept_units_L2 <- matrix(rbinom(n=length(W_2), size = 1, prob = keep_rate), ncol = ncol(W_2)) kept_units_L2 # The matrix of weights W1 then becomes: W_1_dropout <- W_1 * kept_units_L1 W_1_dropout # Scaling W_1_dropout <- W_1_dropout/keep_rate W_1_dropout # The activations then become: A_k_layer_1_dropout <- relu_f(x %*% W_1_dropout + bias_0) A_k_layer_1_dropout # The matrix of weights W2 becomes: W_2_dropout <- W_2 * kept_units_L2 / keep_rate W_2_dropout # The activations of the second layer are then: A_k_layer_2_dropout <- relu_f(A_k_layer_1_dropout %*% W_2_dropout + bias_1) A_k_layer_2_dropout # Then the prediction becomes: pred_dropout <- A_k_layer_2_dropout %*% beta + bias_2 pred_dropout # Training a model with dropout model_structure_bike_dropout <- function(train_data) { model <- keras_model_sequential() %>% # First hidden layer layer_dense(units = 32, activation = "relu", input_shape = dim(train_data)[[2]]) %>% # Second hidden layer layer_dropout(rate = .4) %>% layer_dense(units = 64, activation = "relu") %>% # Third hidden layer layer_dropout(rate = .4) %>% layer_dense(units = 32, activation = "relu") %>% # Output layer layer_dense(units = 1, activation = "sigmoid") model %>% compile( optimizer = "rmsprop", loss = "binary_crossentropy", metrics = c("accuracy", metric_true_positives(name = "tp"), metric_true_negatives(name = "tn"), metric_false_positives(name = "fp"), metric_false_negatives(name = "fn")) ) } # The model can be trained using 3-fold cross validation and a dropout rate of .2 # applied to the second and third hidden layers. # As earlier, the model is trained over 300 epochs, # using a batch size of 512 observations. num_epochs <- 300 batch_size <- 512 metrics_dropout <- vector(mode = "list", length = nb_folds) for(k in 1:nb_folds){ cat(str_c("\nFold number: ", k, "\n")) ind_current <- folds != k # Train set x_train_partial <- x_train[which(ind_current),] y_train_partial <- y_train[which(ind_current)] # Validation set x_train_valid <- x_train[which(!ind_current),] y_train_valid <- y_train[which(!ind_current)] # Building the model model_with_dropout <- model_structure_bike_dropout(x_train_partial) # Training the model on the training set history <- model_with_dropout %>% fit(x_train_partial, y_train_partial, validation_data = list(x_train_valid, y_train_valid), epochs = num_epochs, batch_size = batch_size, verbose = 0) # Store loss and metrics metrics_dropout[[k]] <- history } save(metrics_dropout, file = "deep_learning_models/metrics_dropout.rda") # load("deep_learning_models/metrics_dropout.rda") df_metrics_regul <- history_to_tibble(metrics_without_regul, "accuracy") %>% mutate(type = "Without regularisation") %>% bind_rows( history_to_tibble(metrics_regul, "accuracy") %>% mutate(type = "L2 regularisation (l=0.01)") ) %>% bind_rows( history_to_tibble(metrics_dropout, "accuracy") %>% mutate(type = "Dropout regularisation (rate=.4)") ) ggplot(data = df_metrics_regul, mapping = aes(x = epoch, y = accuracy, linetype = sample, colour = type)) + geom_line() + scale_colour_manual( "Regularisation", values = c("Without regularisation" = "#000000", "L2 regularisation (l=0.01)" = "#E69F00", "Dropout regularisation (rate=.4)" = "#009E73")) + scale_linetype_discrete("Sample") + labs(x = "Epoch", y = "Accuracy (validation sample)", title = "Accuracy with or without regularisation (L2 or dropout).") + theme(legend.position = "bottom", legend.box="vertical", plot.title.position = "plot") # 2 Recurrent Neural Networks --------------------------------------------- # 2.1 Practice with Keras ------------------------------------------------- # Creating other variables from `hour` bike <- bike %>% mutate(day_moment = case_when( hour > 18 & hour <= 22 ~ "evening", hour > 22 | hour <= 6 ~ "night", hour > 6 & hour <= 18 ~ "day", TRUE ~ "problem" )) %>% mutate(day_peak = case_when( hour == 18 ~ "evening", hour == 8 ~ "morning", TRUE ~ "none" )) prop.table(table(bike$day_moment)) prop.table(table(bike$day_peak)) range(bike$date) # removing the unnecessary variables from the sample, # and creating dummy variables for all qualitative variables. data_bikes <- bike %>% select(-date, -y_binary, -month, -seasons, -hour) %>% dummy_cols( remove_first_dummy = TRUE, remove_selected_columns = TRUE) # 2.1.1 Preparing Training, Validation and Test Datasets ------------------ # Delimitations for train/valid/test samples ind_train <- 1:5000 ind_valid <- 5001:6500 # Position of the target variable ind_target <- which(colnames(data_bikes) == "rented_bike_count") ind_target # Train set train_data <- data_bikes[ind_train, -ind_target] # Preparing scaling (mean <- apply(train_data, 2, mean)) (std <- apply(train_data, 2, sd)) # Scaling all predictors data_bikes_x <- data_bikes[, -ind_target] data_bikes_x <- scale(data_bikes_x, center = mean, scale = std) # Scaling the target differently (count data) max_target <- max(data_bikes[[ind_target]]) target_bikes <- data_bikes[[ind_target]]/max_target # Reunite target and predictors data_bikes <- cbind(target_bikes, data_bikes_x) ind_target <- 1 # Let us say that we want to: # # - make predictions 24 hours ahead: we will set a parameter # named delay to 24, as the unit of time in the bike dataset # is the hour # - rely on the previous 3-days data: we will set a parameter named # lookback to 3*24 # - sample the observations every steps=1 time steps. lookback <- 3*24 step <- 1 delay <- 24 batch_size <- 32 #' WARNING: as variable `i_gen` is updated in the #' global environment you should make sure that no variable #' with the same name is used in your code. #' Otherwise, some conflicts may appear. #' Function from Chollet & Allaire (2018) #' @param data Data frame with all the observations in rows #' @param lookback Number of lags #' @param delay Horizon for the prediction #' @param min_index,max_index Index of the first and last #' examples in `data` in which sampling #' is to take place #' @param shuffle Should the rows of the dataset from which #' samples are drawn be randomly swapped before #' sampling? #' @param batch_size Number of samples per batch #' @param step Period, in time steps, at which samples are drawn #' for daily data, if `step=3`, samples are drawn #' every 3 days #' @param ind_target col index of the target variable in `data` generator <- function(data, lookback, delay, min_index, max_index, shuffle = FALSE, batch_size = 128, step = 6, ind_target) { if (is.null(max_index)) max_index <- nrow(data) - delay - 1 i_gen <- min_index + lookback function() { if (shuffle) { rows <- sample(c((min_index+lookback):max_index), size = batch_size) } else { if (i_gen + batch_size >= max_index) i_gen <<- min_index + lookback rows <- c(i_gen:min(i_gen+batch_size, max_index)) i_gen <<- i_gen + length(rows) } samples <- array(0, dim = c(length(rows), lookback / step, dim(data)[[-1]])) targets <- array(0, dim = c(length(rows))) # targets_j <- array(0, dim = c(length(rows))) for (j in 1:length(rows)) { indices <- seq(rows[[j]] - lookback, rows[[j]], length.out = dim(samples)[[2]]) samples[j,,] <- data[indices,] targets[[j]] <- data[rows[[j]] + delay, ind_target] # targets_j[[j]] <- rows[[j]] + delay } list(samples, targets) } } # Let us define an object that will draw observations # from the training set when it will be called during # the training of the network: train_gen <- generator( data = data_bikes, lookback = lookback, delay = delay, min_index = 1, max_index = max(ind_train), shuffle = TRUE, step = step, batch_size = batch_size, ind_target = ind_target ) val_gen <- generator( data = data_bikes, lookback = lookback, delay = delay, min_index = min(ind_valid), max_index = max(ind_valid), shuffle = FALSE, step = step, batch_size = batch_size, ind_target = ind_target ) test_gen <- generator( data = data_bikes, lookback = lookback, delay = delay, min_index = max(ind_valid)+1, max_index = NULL, shuffle = FALSE, step = step, batch_size = batch_size, ind_target = ind_target ) # how many steps are needed to draw from these functions # in order to see the entire validation and entire set, respectively. val_steps <- (max(ind_valid)-min(ind_valid)-lookback)/batch_size val_steps test_steps <- (nrow(data_bikes)-(max(ind_valid)+1)-lookback)/batch_size test_steps # function to create the complete set of tensors from the training sample: generate_data_full <- function(data, lookback, delay, min_index, max_index, step = 1, ind_target) { rows_min <- max(min_index, lookback) rows_max <- max_index-delay rows <- rows_min:rows_max function(){ samples <- array(0, dim = c(length(rows), lookback / step, dim(data)[[-1]])) targets <- array(0, dim = c(length(rows))) for (j in 1:length(rows)) { indices <- seq(rows[[j]] - (lookback-1), rows[[j]], length.out = dim(samples)[[2]]) samples[j,,] <- data[indices,] targets[[j]] <- data[rows[[j]]+delay, ind_target] } list(samples, targets) } } # Indices of observations in the test sample ind_test <- (last(ind_valid)+1):nrow(data_bikes) data_test <- generate_data_full( data = data_bikes, lookback = lookback, delay = delay, min_index = ind_test[1], max_index = last(ind_test), step = 1, ind_target = ind_target ) c(x_data_test, y_data_test) %<-% data_test() dim(x_data_test) dim(y_data_test) # 2.1.2 Naive Benchmark --------------------------------------------------- #' Computes the average of MAE obtained on batches drawn #' from the validation sample. #' The number of batches drawn is `val_steps` #' Function from Chollet & Allaire (2018) evaluate_naive_method <- function() { batch_maes <- rep(NA, val_steps) for (step in 1:val_steps) { # Gen. batches of data from the validation sample c(samples, targets) %<-% val_gen() # Take the last observed values for each example # in the batch, i.e., target_{t-1} preds <- samples[, dim(samples)[[2]], ind_target] # CMAE using target_{t-1} as the predicted value mae <- mean(abs(preds - targets)) batch_maes[step] <- mae } # Avg of MAE obtained over the different batches mean(batch_maes) } batch_mae <- evaluate_naive_method() batch_mae # Do not forget the target ws scaled: batch_mae*max_target # 2.1.3 A Densely Connected Model ----------------------------------------- model_benchmark <- keras_model_sequential() %>% layer_flatten( input_shape = c(lookback / step, dim(data_bikes)[-1])) %>% layer_dense(units = 16, activation = "relu") %>% layer_dense(units = 1, activation = "relu") model_benchmark %>% compile( optimizer = "rmsprop", loss = "mae" ) history_benchmark <- model_benchmark %>% fit( train_gen, steps_per_epoch = 60, epochs = 50, validation_data = val_gen, validation_steps = val_steps ) model_benchmark %>% save_model_hdf5( "deep_learning_models/bikes_model_benchmark.h5") save(history_benchmark, file = "deep_learning_models/history_benchmark.rda") # model_benchmark <- # load_model_hdf5("deep_learning_models/bikes_model_benchmark.h5") # load("deep_learning_models/history_benchmark.rda") # MAE computed both on the training and the validation sets # during the training. plot(history_benchmark) # MAE on the test set: mae_nn <- evaluate(model_benchmark, x=x_data_test, y=y_data_test) mae_nn*max_target # Predictions on test data preds_nn <- predict(model_benchmark, x = x_data_test) # Let us have a look at the predicted values hour by hour over the test sample. ind_in_data <- ind_test[(delay+1):length(ind_test)] df_test_pred_nn <- tibble(obs = bike$rented_bike_count[ind_in_data], pred = preds_nn[,1]*max_target) %>% mutate(date = bike$date[ind_in_data], hour = bike$hour[ind_in_data], hm = lubridate::ymd_h(str_c(date, " ", hour)) ) %>% mutate(month = lubridate::month(date)) df_test_pred_nn %>% pivot_longer(cols = c(obs, pred), names_to = "type", values_to = "bike") %>% ggplot(data = ., mapping = aes(x = hm, y = bike, colour = type, linetype = type)) + geom_line() + facet_wrap(~month, scale = "free", ncol=1) + labs(x = "Date", y = "Number of rented bikes", title = "Predictions for the densely connected network, by month.") + scale_colour_manual( NULL, values = c("obs" = "#648FFF", "pred" = "#DC267F")) + scale_linetype_discrete(NULL) + theme(plot.title.position = "plot") # Residuals: df_test_pred_nn %>% mutate(residuals = obs-pred) %>% ggplot(data = ., mapping = aes(x=hm, y = residuals)) + geom_line() + geom_hline(yintercept = 0, colour = "#648FFF", linetype = "dashed") + labs(x = "Date", y = "Residuals") # 2.1.4 Recurrent Network without Dropout --------------------------------- model_rnn <- keras_model_sequential() %>% layer_gru( units = 16, input_shape = list(NULL, dim(data_bikes)[[-1]])) %>% layer_dense(units = 1, activation = "relu") model_rnn %>% compile( optimizer = "rmsprop", loss = "mae" ) history_rnn <- model_rnn %>% fit( train_gen, steps_per_epoch = 60, epochs = 50, validation_data = val_gen, validation_steps = val_steps ) model_rnn %>% save_model_hdf5("deep_learning_models/model_rnn.h5") save(history_rnn, file = "deep_learning_models/history_rnn.rda") # model_rnn <- load_model_hdf5("deep_learning_models/model_rnn.h5") # load("deep_learning_models/history_rnn.rda") plot(history_rnn) + labs(title = "Training and validation MAE on bike data (first RNN).") + theme(plot.title.position = "plot") # performances of the model on the test data: mae_rnn <- evaluate(model_rnn, x=x_data_test, y=y_data_test) mae_rnn*max_target # Predictions on test data: preds_rnn <- predict(model_rnn, x=x_data_test) ind_in_data <- ind_test[(delay+1):length(ind_test)] df_test_pred_rnn <- tibble(obs = bike$rented_bike_count[ind_in_data], pred = preds_rnn[,1]*max_target) %>% mutate(date = bike$date[ind_in_data], hour = bike$hour[ind_in_data], hm = lubridate::ymd_h(str_c(date, " ", hour)) ) %>% mutate(month = lubridate::month(date)) df_test_pred_rnn %>% pivot_longer(cols = c(obs, pred), names_to = "type", values_to = "bike") %>% ggplot(data = ., mapping = aes(x = hm, y = bike, colour = type, linetype = type)) + geom_line() + facet_wrap(~month, scale = "free", ncol=1) + labs(x = "Date", y = "Number of rented bikes", title = "Predictions for the first recurrent network, by month.") + scale_colour_manual( NULL, values = c("obs" = "#648FFF", "pred" = "#DC267F")) + scale_linetype_discrete(NULL) + theme(plot.title.position = "plot") # Residuals df_test_pred_rnn %>% mutate(residuals = obs-pred) %>% ggplot(data = ., mapping = aes(x=hm, y = residuals)) + geom_line() + geom_hline(yintercept = 0, colour = "#648FFF", linetype = "dashed") + labs(x = "Date", y = "Residuals", title = "Residuals with the first recurrent network.") + theme(plot.title.position = "plot") # 2.1.5 Recurrent Network with Dropout ------------------------------------ model_rnn_drop <- keras_model_sequential() %>% layer_lstm( units = 16, dropout = 0.3, recurrent_dropout = 0.4, input_shape = list(NULL, dim(data_bikes)[[-1]])) %>% layer_dense(units = 1, activation = "relu") model_rnn_drop %>% compile( optimizer = "rmsprop", loss = "mae" ) history_rnn_drop <- model_rnn_drop %>% fit( train_gen, steps_per_epoch = 60, epochs = 100, validation_data = val_gen, validation_steps = val_steps ) model_rnn_drop %>% save_model_hdf5("deep_learning_models/model_rnn_drop.h5") save(history_rnn_drop, file = "deep_learning_models/history_rnn_drop.rda") # model_rnn_drop <- load_model_hdf5("deep_learning_models/model_rnn_drop.h5") # load("deep_learning_models/history_rnn_drop.rda") plot(history_rnn_drop) + labs(titl = "Training and validation MAE on bike data, (regularised RNN).") + theme(plot.title.position = "plot") # MAE on the test set: mae_rnn_drop <- evaluate(model_rnn_drop, x=x_data_test, y=y_data_test) mae_rnn_drop*max_target # Predictions on test data preds_rnn_drop <- predict(model_rnn_drop, x = x_data_test) ind_in_data <- ind_test[(delay+1):length(ind_test)] df_test_pred_rnn_drop <- tibble(obs = bike$rented_bike_count[ind_in_data], pred = preds_rnn_drop[,1]*max_target) %>% mutate(date = bike$date[ind_in_data], hour = bike$hour[ind_in_data], hm = lubridate::ymd_h(str_c(date, " ", hour)) ) %>% mutate(month = lubridate::month(date)) df_test_pred_rnn_drop %>% pivot_longer(cols = c(obs, pred), names_to = "type", values_to = "bike") %>% ggplot(data = ., mapping = aes(x = hm, y = bike, colour = type, linetype = type)) + geom_line() + facet_wrap(~month, scale = "free", ncol=1) + labs(x = "Date", y = "Number of rented bikes", title = "Predictions for the regularised RNN, by month.") + scale_colour_manual( NULL, values = c("obs" = "#648FFF", "pred" = "#DC267F")) + scale_linetype_discrete(NULL) + theme(plot.title.position = "plot") df_test_pred_rnn_drop %>% mutate(residuals = obs-pred) %>% ggplot(data = ., mapping = aes(x=hm, y = residuals)) + geom_line() + geom_hline(yintercept = 0, colour = "#648FFF", linetype = "dashed") + labs(x = "Date", y = "Residuals", title = "Residuals with the regularised recurrent network.") + theme(plot.title.position = "plot")