---
title: "Neural Networks - Boston Housing Price - Densily connected NN"
subtitle: "TMA4268 Statistical Learning V2018"
author: "Thiago G. Martins, Department of Mathematical Sciences, NTNU"
date: "Spring 2018"
output: #3rd letter intentation hierarchy
#  beamer_presentation:
###    incremental: true # or >* for one at a time
#  slidy_presentation:
#    font_adjustment: +1  
   prettydoc::html_pretty:
    theme: cayman
    highlight: github
    toc: true
    toc_depth: 4
   # prettydoc::html_pretty:
   #  theme: architect
   #  highlight: github
#  pdf_document:
  #  toc: true
  #  toc_depth: 2

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, include=FALSE}
library(keras)
```

## The Boston housing price dataset

The objective here is to predict the median price of homes in a given Boston suburb in the mid-1970s, given data points about the suburb at the time, such as the crime rate, the local property tax rate, and so on. 

### Preparing the data

* Few data points: only 506, split between 404 training samples and 102 test samples. 

* Each feature in the input data (for example, the crime rate) has a different scale. 

* For instance, some values are proportions, which take values between 0 and 1; others take values between 1 and 12, others between 0 and 100, and so on.

```{r, message=FALSE, warning=FALSE}
dataset <- dataset_boston_housing()
c(c(train_data, train_targets), c(test_data, test_targets)) %<-% dataset
```

```{r, message=FALSE, warning=FALSE}
str(train_data)
str(test_data)
```

The targets are the median values of owner-occupied homes, in thousands of dollars:

```{r, message=FALSE, warning=FALSE}
str(train_targets)
```

A widespread best practice to deal with such data is to do feature-wise normalization.

```{r, message=FALSE, warning=FALSE}
mean <- apply(train_data, 2, mean)                                
std <- apply(train_data, 2, sd)
train_data <- scale(train_data, center = mean, scale = std)       
test_data <- scale(test_data, center = mean, scale = std)
```

Note that the quantities used for normalizing the test data are computed using the training data. You should never use in your workflow any quantity computed on the test data, even for something as simple as data normalization.

### Model building

Because so few samples are available, you’ll use a very small network with two hidden layers, each with 64 units.

Because you’ll need to instantiate the same model multiple times, you use a function to construct it.

```{r, message=FALSE, warning=FALSE}
build_model <- function() {                               
  model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation = "relu",
                input_shape = dim(train_data)[[2]]) %>%
    layer_dense(units = 64, activation = "relu") %>%
    layer_dense(units = 1)
  model %>% compile(
    optimizer = "rmsprop",
    loss = "mse",
    metrics = c("mae")
  )
}
```

The network ends with a single unit and no activation (it will be a linear layer). 

You're also monitoring a new metric during training: mean absolute error (MAE). It’s the absolute value of the difference between the predictions and the targets. For instance, an MAE of 0.5 on this problem would mean your predictions are off by $500 on average.

### Model validation

K-fold CV: Small dataset

```{r, message=FALSE, warning=FALSE}
# randomly allocate each sample to a particular fold
k <- 4
indices <- sample(1:nrow(train_data))
folds <- cut(indices, breaks = k, labels = FALSE)
```

```{r, message=FALSE, warning=FALSE}
# Run CV and keep the MAE scores
num_epochs <- 100
all_scores <- c()
for (i in 1:k) {
  cat("processing fold #", i, "\n")

  val_indices <- which(folds == i, arr.ind = TRUE)                     
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  partial_train_data <- train_data[-val_indices,]                      
  partial_train_targets <- train_targets[-val_indices]

  model <- build_model()                                               

  model %>% fit(partial_train_data, partial_train_targets,             
                epochs = num_epochs, batch_size = 1, verbose = 0)

  results <- model %>% evaluate(val_data, val_targets, verbose = 0)    
  all_scores <- c(all_scores, results$mean_absolute_error)
}
```

Running this with num_epochs = 100 yields the following results:

```{r, message=FALSE, warning=FALSE}
all_scores
mean(all_scores)
```

Let’s try training the network a bit longer: 500 epochs. To keep a record of how well the model does at each epoch, you’ll modify the training loop to save the per-epoch validation score log.

```{r, message=FALSE, warning=FALSE}
num_epochs <- 500
all_mae_histories <- NULL
for (i in 1:k) {
  cat("processing fold #", i, "\n")

  val_indices <- which(folds == i, arr.ind = TRUE)              
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]

  partial_train_data <- train_data[-val_indices,]               
  partial_train_targets <- train_targets[-val_indices]

  model <- build_model()                                        

  history <- model %>% fit(                                     
    partial_train_data, partial_train_targets,
    validation_data = list(val_data, val_targets),
    epochs = num_epochs, batch_size = 1, verbose = 0
  )
  mae_history <- history$metrics$val_mean_absolute_error
  all_mae_histories <- rbind(all_mae_histories, mae_history)
}
```

You can then compute the average of the per-epoch MAE scores for all folds.

```{r, message=FALSE, warning=FALSE}
average_mae_history <- data.frame(
  epoch = seq(1:ncol(all_mae_histories)),
  validation_mae = apply(all_mae_histories, 2, mean)
)
```

Let’s plot this:

```{r, message=FALSE, warning=FALSE}
library(ggplot2)
ggplot(average_mae_history, aes(x = epoch, y = validation_mae)) + geom_line()
```

Let's use `geom_smooth()` to try to get a clearer picture:

```{r, message=FALSE, warning=FALSE}
ggplot(average_mae_history, aes(x = epoch, y = validation_mae)) + geom_smooth()
```

According to this plot, validation MAE stops improving significantly after 125 epochs. Past that point, you start overfitting.

Once you're finished tuning other parameters of the model (in addition to the number of epochs, you could also adjust the size of the hidden layers), you can train a final production model on all of the training data, with the best parameters, and then look at its performance on the test data.

```{r, message=FALSE, warning=FALSE}
model <- build_model()
model %>% fit(train_data, train_targets,                  
          epochs = 80, batch_size = 16, verbose = 0)
result <- model %>% evaluate(test_data, test_targets)

result
```

You’re still off by about $2,540.

### Recommended exercise 3

Try to improve the MAE on the test data by changing model parameters such as number of layers and number of hidden units in each layer.

