#####################################################################################################
##
##                                    Deep Learning
##
## E. Flachaire, 20/12/2021
##
#####################################################################################################


############################## Simulation experiment

set.seed(1) ; n=200
x=sort(runif(n))
y=sin(12*(x+.2))/(x+.2) + rnorm(n)/2
df=data.frame(y,x)

# fit a polynomial regression with $M=3$:
ols=lm(y~x+I(x^2)+I(x^3))
plot(x,y, main="Polynomial: M=3")
lines(x,predict(ols), col="blue")

# fit a neural network model with $M=3$:
library(neuralnet)
nn=neuralnet(y~x, data=df, hidden=3, threshold=.05) 
yfit=compute(nn,data.frame(x))$net.result
plot(x,y, main="Neural Networks: M=3")
lines(x,yfit, col="red")



########### Application 1: Mincer equation

library(AER) 
data("CPS1985")
CPS1985$gender=as.numeric(CPS1985$gender)
library(neuralnet)

# Single Hidden layer neural network
nn=neuralnet(log(wage)~education+experience+gender,data=CPS1985,hidden=3,threshold=.05)
plot(nn)

# Multilayer Neural Networks
nn=neuralnet(log(wage)~education+experience+gender,data=CPS1985,hidden=c(3,3),threshold=.05)
plot(nn)
nn=neuralnet(log(wage)~education+experience+gender,data=CPS1985,hidden=c(3,3,3),threshold=.05)
plot(nn)



########### Application 2: MNIST handwritten digit dataset (ISLR section 10.9.2)

library(keras)
# load the MNIST digit data
mnist <- dataset_mnist() 
x_train <- mnist$train$x 
g_train <- mnist$train$y 
x_test <- mnist$test$x
g_test <- mnist$test$y
# reshapei inmages into matrices
x_train <- array_reshape(x_train, c(nrow(x_train), 784)) 
x_test <- array_reshape(x_test, c(nrow(x_test), 784))
y_train <- to_categorical(g_train, 10)
y_test <- to_categorical(g_test, 10)
#rescale to the unit interval
x_train <- x_train / 255 
x_test <- x_test / 255
# define the multilayer NN
modelnn <- keras_model_sequential()
modelnn %>%
  layer_dense(units = 256, activation = "relu",
                input_shape = c(784)) %>%
  layer_dropout(rate = 0.4) %>%
  layer_dense(units = 128, activation = "relu") %>%
  layer_dropout(rate = 0.3) %>%
  layer_dense(units = 10, activation = "softmax")
summary(modelnn)
# Add details to the model
modelnn %>% compile(loss = "categorical_crossentropy", 
                    optimizer = optimizer_rmsprop(), metrics = c("accuracy")
)
# Fit the NN with training data
system.time(
  history <- modelnn %>%
    fit(x_train, y_train, epochs = 30, batch_size = 128,
          validation_split = 0.2) 
)
plot(history, smooth = FALSE)
# Obtain the test error
accuracy <- function(pred, truth)
  mean(drop(pred) == drop(truth))
modelnn %>% predict(x_test) %>% max.col %>% accuracy(g_test+1)

# fit a multinomial logit as a NN without hidden layer
modellr <- keras_model_sequential() %>%
  layer_dense(input_shape = 784, units = 10,
              activation = "softmax") 
summary(modellr)
modellr %>% compile(loss = "categorical_crossentropy", optimizer = optimizer_rmsprop(), metrics = c("accuracy"))
modellr %>% fit(x_train, y_train, epochs = 30, batch_size = 128, validation_split = 0.2)
modellr %>% predict(x_test) %>% max.col %>% accuracy(g_test+1)