MNIST ("Modified National Institute of Standards and Technology") is computer vision dataset released in 1999. It contains data of handwritten images and it is the "de facto" benchmark for classification algorithms. The goal is to correctly identify digits from a dataset of tens of thousands of handwritten images.

The data

The data description can be found at Kaggle:

The data files train.csv and test.csv contain gray-scale images of hand-drawn digits, from zero through nine.

Each image is 28 pixels in height and 28 pixels in width, for a total of 784 pixels in total. Each pixel has a single pixel-value associated with it, indicating the lightness or darkness of that pixel, with higher numbers meaning darker. This pixel-value is an integer between 0 and 255, inclusive.

The training data set, (train.csv), has 785 columns. The first column, called "label", is the digit that was drawn by the user. The rest of the columns contain the pixel-values of the associated image.

Each pixel column in the training set has a name like pixelx, where x is an integer between 0 and 783, inclusive. To locate this pixel on the image, suppose that we have decomposed x as x = i * 28 + j, where i and j are integers between 0 and 27, inclusive. Then pixelx is located on row i and column j of a 28 x 28 matrix, (indexing by zero).

Exploration


In [29]:
library(tidyverse)
library(purrrlyr)
library(forcats) # factors munging
library(ggthemes) # visualization
library(scales) # visualization
library(caret) # ML
library(Rtsne)
library(viridis)
library(fastknn)

theme_update(plot.title = element_text(hjust = 0.5))

In [2]:
train <- read_csv("data/train.csv")
test <- read_csv("data/test.csv")


Parsed with column specification:
cols(
  .default = col_integer()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  .default = col_integer()
)
See spec(...) for full column specifications.

In [3]:
head(train)


labelpixel0pixel1pixel2pixel3pixel4pixel5pixel6pixel7pixel8pixel774pixel775pixel776pixel777pixel778pixel779pixel780pixel781pixel782pixel783
10000000000000000000
00000000000000000000
10000000000000000000
40000000000000000000
00000000000000000000
00000000000000000000

In [4]:
dim(train)


  1. 42000
  2. 785

In [5]:
dim(test)


  1. 28000
  2. 784

In [6]:
train <- train %>%
    mutate(label = factor(label))

In [7]:
train <- train %>%
    mutate(intensity = 
           select(., starts_with("pixel")) %>%
               rowMeans()
          ) %>%
    select(intensity, everything())

In [8]:
options(repr.plot.width=7, repr.plot.height=5)

ggplot(data=train, aes(x=label)) + 
    geom_bar(aes(y = (..count..)/sum(..count..))) +
    ylab("percentage of train data") +
    xlab("digit") +
    ggtitle("Distribution of digits") +
    scale_y_continuous(labels = percent, limits=c(0, 0.12)) +
    theme_hc()



In [9]:
options(repr.plot.width=7, repr.plot.height=5)

ggplot(train, aes(x = intensity, y = ..density..)) +
    geom_density(aes(fill = label), alpha = 0.3) +
    ggtitle("Distribution of intensity per digit") +
    theme_hc()



In [10]:
ggplot(train, aes(x = label, y = intensity)) +
    geom_boxplot(aes(fill = label), alpha = 0.5) +
    xlab("digit") +
    theme_hc()



In [11]:
options(repr.plot.width=5, repr.plot.height=5)

flip <- function(matrix) {
      apply(matrix, 2, rev)
}

plot_digit <- function (obs) {
    m <- flip(matrix(rev(as.numeric(obs[-1])), 28, 28))
    image(m, axes = FALSE, col = grey(seq(0, 1, length = 256)))
    title(main = obs$label)
}

train[4, ] %>%
    select(-intensity) %>%
    plot_digit



In [12]:
digit_groups <- train %>%
    group_by(label)

Mean


In [13]:
options(repr.plot.width=7, repr.plot.height=10)

par(mfrow=c(4,3), mar=c(1.5, 1.5, 1.5, 1.5))

r <- digit_groups %>%
    summarise_all(funs(mean)) %>%
    select(-intensity) %>%
    by_row(plot_digit)


Median


In [14]:
par(mfrow=c(4,3), mar=c(1.5, 1.5, 1.5, 1.5))

r <- digit_groups %>%
    summarise_all(funs(median)) %>%
    select(-intensity) %>%
    by_row(plot_digit)


Standard deviation


In [15]:
par(mfrow=c(4,3), mar=c(1.5, 1.5, 1.5, 1.5))

r <- digit_groups %>%
    summarise_all(funs(sd)) %>%
    select(-intensity) %>%
    by_row(plot_digit)


Missing data?


In [16]:
anyNA(train)


FALSE

Visualizing using t-SNE

t-SNE (t-Distributed Stochastic Neighbour) allows to map high-dimensional data to a 2D or 3D plane. It is most commonly used for visualization purposes. However, the original algorithm is computationally expensive $O(N^4)$ with the number of samples in our data.

We will use Barnes-Hut-SNE algorithm which is an improvement over the standard t-SNE and runs in O(N log N) time. The algorithm is implemented in the Rtsne package.


In [17]:
nrow(train)


42000

In [18]:
train_sample <- sample_n(train, 5000)

In [19]:
tsne_data <- as.matrix(select(train_sample, -intensity, -label))
tsne_out <- Rtsne(tsne_data, check_duplicates = FALSE, pca = TRUE, dims = 2)

In [60]:
options(repr.plot.width=7, repr.plot.height=6)

plot_data <- tibble(
    x = tsne_out$Y[,1],
    y = tsne_out$Y[,2],
    digit = train_sample$label
)

ggplot(plot_data) + 
    geom_point(aes(x=x, y=y, color=digit)) +
    ggtitle("t-SNE clustering of MNIST digits") +
    theme_hc()


Building a model


In [21]:
train_index <- createDataPartition(train$label, p=0.8, list=FALSE)
cv_train <- train[train_index, ]
cv_validation <- train[-train_index, ]

In [22]:
train_data <- select(cv_train, -label, -intensity)
train_labels <- cv_train$label

validation_data <- select(cv_validation, -label, -intensity)
validation_labels <- cv_validation$label

In [32]:
knn_model <- fastknn(
    as.matrix(train_data), 
    train_labels,
    as.matrix(validation_data), 
    k=10,
    method = "dist"
)

Evaluation


In [24]:
confusion_matrix <- caret::confusionMatrix(
    table(predicted=knn_model$class, actual=validation_labels)
)

In [25]:
confusion_matrix


Confusion Matrix and Statistics

         actual
predicted   0   1   2   3   4   5   6   7   8   9
        0 821   0   5   3   1   1   3   0   4   2
        1   0 934  11   3  12   1   2  12  12   2
        2   0   0 788   3   0   1   0   2   4   0
        3   0   0   2 833   0  15   0   0  20   5
        4   0   0   2   1 776   2   1   4   7   8
        5   2   0   1  10   0 718   2   0  13   1
        6   2   1   1   2   3  12 818   0   3   2
        7   1   0  21   7   1   0   0 847   4  14
        8   0   0   3   4   1   0   1   0 737   1
        9   0   1   1   4  20   9   0  15   8 802

Overall Statistics
                                          
               Accuracy : 0.9616          
                 95% CI : (0.9573, 0.9657)
    No Information Rate : 0.1115          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.9574          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: 0 Class: 1 Class: 2 Class: 3 Class: 4 Class: 5
Sensitivity           0.99395   0.9979  0.94371  0.95747  0.95332  0.94598
Specificity           0.99749   0.9926  0.99868  0.99442  0.99670  0.99620
Pos Pred Value        0.97738   0.9444  0.98747  0.95200  0.96879  0.96118
Neg Pred Value        0.99934   0.9997  0.99381  0.99508  0.99500  0.99464
Prevalence            0.09838   0.1115  0.09945  0.10362  0.09695  0.09040
Detection Rate        0.09778   0.1112  0.09385  0.09921  0.09242  0.08552
Detection Prevalence  0.10005   0.1178  0.09505  0.10422  0.09540  0.08897
Balanced Accuracy     0.99572   0.9952  0.97119  0.97595  0.97501  0.97109
                     Class: 6 Class: 7 Class: 8 Class: 9
Sensitivity           0.98912   0.9625  0.90764  0.95818
Specificity           0.99656   0.9936  0.99868  0.99233
Pos Pred Value        0.96919   0.9464  0.98661  0.93256
Neg Pred Value        0.99881   0.9956  0.99019  0.99536
Prevalence            0.09850   0.1048  0.09671  0.09969
Detection Rate        0.09743   0.1009  0.08778  0.09552
Detection Prevalence  0.10052   0.1066  0.08897  0.10243
Balanced Accuracy     0.99284   0.9781  0.95316  0.97526

In [49]:
confusion_data <- as.tibble(confusion_matrix$table) %>%
    mutate(rate = n/sum(n)) %>%
    mutate(error_rate = ifelse(actual == predicted,
                             0, rate))

In [61]:
options(repr.plot.width=8, repr.plot.height=6)

ggplot(confusion_data) + 
    geom_tile(aes(x=actual, y=predicted, fill=n)) +
    ggtitle("Absolute predicted vs actual values") +
    scale_x_discrete(name="Actual") + 
    scale_y_discrete(name="Predicted") +
    scale_color_viridis() +
    scale_fill_viridis(option="viridis") +
    theme_hc()



In [83]:
max_scale = max(confusion_data$rate)
mid_scale = max_scale / 2

ggplot(confusion_data) + 
    geom_tile(aes(x=actual, y=predicted, fill=rate)) +
    ggtitle("Weighted predicted vs actual values") +
    scale_x_discrete(name="Actual") + 
    scale_y_discrete(name="Predicted") +
    scale_color_viridis() +
    scale_fill_viridis(
        option="viridis",
        breaks=c(0, mid_scale, max_scale),
        labels=c(0, round(mid_scale, 5), round(max_scale, 5))
    ) + theme_hc()



In [82]:
max_scale = max(confusion_data$error_rate)
mid_scale = max_scale / 2

ggplot(confusion_data) + 
    geom_tile(aes(x=actual, y=predicted, fill=error_rate)) +
    ggtitle("Error rates") +
    scale_x_discrete(name="Actual") + 
    scale_y_discrete(name="Predicted") +
    scale_fill_viridis(
        option="viridis", 
        breaks=c(0, mid_scale, max_scale),
        labels=c(0, round(mid_scale, 5), round(max_scale, 5))
    ) + theme_hc()


Sending our solution


In [27]:
train_data <- select(train, -label, -intensity)
train_labels <- train$label

knn_model <- fastknn(
    as.matrix(train_data), 
    train_labels,
    as.matrix(test), 
    k=10,
    method = "dist"
)

In [28]:
predictions <- tibble(
    ImageId=1:length(knn_model$class), 
    Label=knn_model$class
)
write_csv(predictions, "results/knn.csv")