options(warn=-1)
library(dslabs)
library(tidyverse)
library(caret)
library(broom)
library(matrixStats)
library(purrr)
library(gridExtra)
library(randomForest)
library(Rborist)
options(repr.plot.width=12, repr.plot.height=7)

mnist <- read_mnist()
data("mnist_27")

The MNIST digits data (Modified National Institute of Standards and Technology database digits) is composed of training and testing sets. The training set has 60,000 images and the test set has 10,000 images, each one of them with 784 pixels (28 x 28).

names(mnist)
names(mnist$train)
dim(mnist$train$images)
dim(mnist$test$images)

1. 'train'
2. 'test'

1. 'images'
2. 'labels'

1. 60000
2. 784

1. 10000
2. 784

We start with a general analysis based on the entire train set.
The boxplot below shows the distribution of the total pixels darkness and how it varies from digit to digit:

tibble(labels = as.factor(mnist$train$labels),
       row_averages = rowMeans(mnist$train$images)) %>%
  qplot(labels, row_averages, data = ., geom = "boxplot") +
    theme(axis.text = element_text(size=14))

From the boxplot we can see that the digit 1 has significant less pixels compared to the other digits.

Two labels and two predictors

We will start with only two predictors, to showcase and analyze different machine learning techniques.

To select the two predictors, we group all 784 predictors into two informations:
The proportions of dark pixels in the upper left quadrant and in the lower right quadrant. To calculate the proportions we binarize the data in the preprocessing step.
We show first, from the entire data set, those with largest and smallest proportion of dark pixels in the two quadrants.

# Preprocessing
row_column <- expand.grid(row=1:28, col=1:28)
upper_left_ind <- which(row_column$col <= 14 & row_column$row <= 14)
lower_right_ind <- which(row_column$col > 14 & row_column$row > 14)

xs <- (mnist$train$images > 255/2)*1
xs <- cbind(rowSums(xs[ ,upper_left_ind])/784,
            rowSums(xs[ ,lower_right_ind])/784)
colnames(xs) <- c("x_1", "x_2")

Showing the largest and smallest images ("3" and "1", respectively) with dark pixels proportion in the upper left quadrant:

is <- c(which.min(xs[,1]), which.max(xs[,1]))
titles <- c("smallest","largest")
tmp <- lapply(1:2, function(i){
  expand.grid(Row=1:28, Column=1:28) %>%
    mutate(label=titles[i],
           value = mnist$train$images[is[i],])
})
tmp <- Reduce(rbind, tmp)
tmp %>% ggplot(aes(Row, Column, fill=value)) +
  geom_raster() +
  scale_y_reverse() +
  scale_fill_gradient(low="white", high="black") +
  facet_grid(.~label) +
  geom_vline(xintercept = 14.5) +
  geom_hline(yintercept = 14.5)

Showing the largest and smallest images ("0" and "5", respectively) with dark pixels proportion in the lower right quadrant

is <- c(which.min(xs[,2]), which.max(xs[,2]))
titles <- c("smallest","largest")
tmp <- lapply(1:2, function(i){
  expand.grid(Row=1:28, Column=1:28) %>%
    mutate(label=titles[i],
           value = mnist$train$images[is[i],])
})
tmp <- Reduce(rbind, tmp)
tmp %>% ggplot(aes(Row, Column, fill=value)) +
  geom_raster() +
  scale_y_reverse() +
  scale_fill_gradient(low="white", high="black") +
  facet_grid(.~label) +
  geom_vline(xintercept = 14.5) +
  geom_hline(yintercept = 14.5)

Since only two predictors are unsufficient to classify all ten outcomes, we work first with only two labels: 2 and 7. Later we will work on the entire dataset.

The largest and smallest images ("7" and "2", respectively) with dark pixels proportion in the upper left quadrant:

is <- mnist_27$index_train[c(which.min(mnist_27$train$x_1), which.max(mnist_27$train$x_1))]

tmp <- lapply(1:2, function(i){
  expand.grid(Row=1:28, Column=1:28) %>%
    mutate(label=titles[i],
           value = mnist$train$images[is[i],])
})
tmp <- Reduce(rbind, tmp)
tmp %>% ggplot(aes(Row, Column, fill=value)) +
  geom_raster() +
  scale_y_reverse() +
  scale_fill_gradient(low="white", high="black") +
  facet_grid(.~label) +
  geom_vline(xintercept = 14.5) +
  geom_hline(yintercept = 14.5)

The largest and smallest images (both "7") with dark pixels proportion in the lower right quadrant:

is <- mnist_27$index_train[c(which.min(mnist_27$train$x_2), which.max(mnist_27$train$x_2))]

tmp <- lapply(1:2, function(i){
  expand.grid(Row=1:28, Column=1:28) %>%
    mutate(label=titles[i],
           value = mnist$train$images[is[i],])
})
tmp <- Reduce(rbind, tmp)
tmp %>% ggplot(aes(Row, Column, fill=value)) +
  geom_raster() +
  scale_y_reverse() +
  scale_fill_gradient(low="white", high="black") +
  facet_grid(.~label) +
  geom_vline(xintercept = 14.5) +
  geom_hline(yintercept = 14.5)

All training set is presented here, with the two predictors and labels:

mnist_27$train %>%
  ggplot(aes(x_1, x_2, color = y)) +
  geom_point()

Since this is a controlled study, we will use the true conditional probability, built from both training and test data, to compare each machine learning technique outcome. On real cases we do not have access to this true conditional probability, but here it is used as illustrative purpose.
It can be seen, in this case, the non-linear nature of the true conditional probability.

plot_cond_prob <- function(p_hat=NULL, legend=FALSE, title = "True conditional probability"){
  data <- mnist_27$true_p
  if(!is.null(p_hat)){
    data <- mutate(data, p=p_hat)
  }
  data %>% ggplot(aes(x_1, x_2, z=p, fill=p)) +
    geom_raster(show.legend = legend) +
    scale_fill_gradientn(colors=c("#F8766D","white","#00BFC4")) +
    stat_contour(breaks=c(0.5),color="black") +
    ggtitle(title)
}

options(repr.plot.width=7, repr.plot.height=7)
plot_cond_prob(legend=TRUE)