{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# MNIST digit recognition","metadata":{}},{"cell_type":"markdown","source":"things to do:\n- [done] flesh out some commentary on accuracy vs loss and different scoring procedures\n- need to explain what the layers of code in the CNN model do, with hyperlinks, etc.\n- [done] make tileplot or similar graphic for confusion matrix\n- [done] add data frame with precision, recall, f1 by class\n- add explanation of precision vs recall\n- [done] add \"rogue's gallery\" of incorrectly predicted images\n","metadata":{}},{"cell_type":"markdown","source":"### Setup: install libraries and load data\nIf we ran this code in the R console, we would need to rely on a library called miniconda that basically passes the neural network commands to Python (keras itself was written in Python.) By working in a Kaggle notebook, this dependency is seamless.","metadata":{"execution":{"iopub.status.busy":"2022-12-08T00:58:54.952823Z","iopub.execute_input":"2022-12-08T00:58:54.954779Z","iopub.status.idle":"2022-12-08T00:58:54.969497Z"}}},{"cell_type":"code","source":"library(dplyr)       # tidyverse logic for data frames\nlibrary(ggplot2)     # graphics\nlibrary(keras)       # neural network\nlibrary(tensorflow)  # neural network\nlibrary(readr)       # load csv files\nlibrary(caret)       # Classification And REgression Training - includes confusionMatrix\n\ntrain <- read_csv('/kaggle/input/digit-recognizer/train.csv')\ntest  <- read_csv('/kaggle/input/digit-recognizer/test.csv')\n\noptions(warn=-1)     # turn off some warnings\nverbose = 0\n\nnrow(train)     # our training data set has 42,000 characters with their labels\nnrow(test)      # we are going to train a model to classify these 28,000\n\n# Define some constants here instead of hunting through the code for them later.\nIm_R     <- 28      # of rows in the image\nIm_C     <- 28      # of columns in the image\ntr_num   <- 36000   # size of our training set (rest will be validation)\n\nset.seed(142857)    # seed for random number generation (train / validation split, etc.)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-12-10T19:49:26.390659Z","iopub.execute_input":"2022-12-10T19:49:26.392090Z","iopub.status.idle":"2022-12-10T19:49:31.979446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The label column tells the true value of the digit (0-9) and the remaining 784 columns (pixel0 through pixel783) encode the image. Each image is 28x28 grayscale, so there are 784 pixels, and an intensity (0 = white, 255 = black) for each pixel. The image is encoded left to right and top to bottom, so that pixel0 through pixel27 encode the top row, pixel28 through pixel55 encode the second row, etc. The test data set looks the same, except that it is missing the label column (i.e. you aren't told what digit each label encodes - you're supposed to build a model to figure it out.)","metadata":{}},{"cell_type":"code","source":"head(train)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:31.981969Z","iopub.execute_input":"2022-12-10T19:49:31.983331Z","iopub.status.idle":"2022-12-10T19:49:32.060249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Distribution of values for the label column. It's not exactly 4200 of each, but the labels in the training data set are pretty evenly distributed.\n\nOur model makes probabilistic predictions: e.g. it will look at an image and say \"there is a 90% chance that this is a 9 and a 10% chance that this is an 8.\" There are two main ways that models are evaluated:\n- accuracy: the \"dominant\" (largest value) prediction is 9, so score either 100% if true label is 9 or 0% if true label is not 9\n- loss: measures confidence of prediction, so score 90% if the true label is 9, 10% if the true label is 8, or 0% if the true label is anything else\n\nFor a well-balanced distribution of labels, accuracy is a more appropriate scoring metric, and the Kaggle competition is scored on accuracy.","metadata":{}},{"cell_type":"code","source":"train$label %>% table() %>% barplot(col = \"#47C7FC\", main = \"Distribution of label in Training Set\")","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:32.062010Z","iopub.execute_input":"2022-12-10T19:49:32.063009Z","iopub.status.idle":"2022-12-10T19:49:32.158743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample image\nLet's decode one of these rows of data to see what these images look like. \nFirst, let's define the function \"draw\" that converts a set values into a 28x28 image and then plots the image in grayscale. We can re-use this function in other graphics. Then, we can apply that function to the fourth row of the train dataset, which according to the label in the head above should be a \"4\".","metadata":{}},{"cell_type":"code","source":"draw <- function(x) {                                         # where x is vector of length 784\n    mat <- matrix(as.numeric(x), Im_R, Im_C, byrow = TRUE)    # convert to 28x28\n    image(t(mat)[,Im_R:1], axes = FALSE,                      # plot a 28x28 grayscale image\n          col = grey(seq(1, 0, length = 256)))                \n}\n\nt <- train[4,-1]\ndraw(t)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:32.160429Z","iopub.execute_input":"2022-12-10T19:49:32.161396Z","iopub.status.idle":"2022-12-10T19:49:32.225827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Data Exploration","metadata":{}},{"cell_type":"markdown","source":"Here are 10 of each true label. You can see some of the challenges related to classification here:\n - differences in line thickness, pen pressure, steady hand, slant, etc.\n - variations such as open vs closed \"4\", slashed \"0\" and \"7\", self-intersecting \"2\", etc.\n - schmutz, stray pixels, etc. in some images\n - some of these images are unclear even to the naked eye (e.g. the second \"5\" from the right.)","metadata":{}},{"cell_type":"code","source":"par(mfrow=c(10, 10), pty='s', mai=c(0.0, 0, 0, 0.0))\nfor (row in 0:9) {\n  samp <- train %>% filter(label == row)\n  for (col in 1:10) {\n    t <- samp[col, -1]\n    draw(t)\n    box(lty = 'solid')\n  }\n}","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:32.228197Z","iopub.execute_input":"2022-12-10T19:49:32.229240Z","iopub.status.idle":"2022-12-10T19:49:33.191289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Below: for each row of data, I took the average intensity of all the observations and then made a boxplot by label. In other words, we are looking at how much \"ink\" it takes to write each numeral. Note that the numeral \"1\" takes the least amount of space, while \"0\" seems to take the most. The violin plot is sort of a combination of boxplot and density plot - it has been in the Python Seaborn graphics package for a while and was a recent addition to ggplot2.","metadata":{}},{"cell_type":"code","source":"train$AvgRow <-  rowMeans(train[,-1])   # new column is the average of all the pixels in a row\np <- ggplot(train, aes(group = label, x=as.factor(label), y=AvgRow)) + geom_violin()\np + stat_summary(fun=median, geom=\"point\", size=2, color=\"red\")\n\ntrain <- select(train, -\"AvgRow\")       # drop this column when we're done making the plot","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:33.194411Z","iopub.execute_input":"2022-12-10T19:49:33.195972Z","iopub.status.idle":"2022-12-10T19:49:34.007826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"What if we took the average intensity of every pixel in the 28x28 grid for each true label? The resulting graphics are \"fuzzy\" but your brain really wants to see the digits 0-9 in there! This gives us hope that maybe we can build a classification model.","metadata":{}},{"cell_type":"code","source":"par(mfrow=c(3,4), pty='s', mai=c(0, 0, 0, 0))\n\nfor (lab in 0:9) {\n  subs <- train %>% filter(label == lab)\n  t <- colMeans(subs)[2:785]\n  draw(t)  \n  box(lty = 'solid')\n}","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:34.009844Z","iopub.execute_input":"2022-12-10T19:49:34.011004Z","iopub.status.idle":"2022-12-10T19:49:34.671482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Building the CNN model","metadata":{}},{"cell_type":"code","source":"# The MNIST data is known to be \"clean\" in the sense that there are no missing values to impute or anything \n# like that - but we still need to do some preprocessing: rescale the data, convert to data.matrix, split our \n# training observations into \"train\" and \"holdout\" subsets and reshape the 1x784 arrays into 28x28 arrays.\n\n# Best practice is to rescale all numbers from 0 to 1 - so let's divide our pixel values by 255.\ntrain[, c(2:ncol(train))] <- train[, c(2:ncol(train))]/255\ntest      <- test/255\n\n# Data Matrix is a more efficient way to encode sparse data.\ntrain.mtx <- data.matrix(train)\ntest.mtx  <- data.matrix(test)\n\n# Next step is to create a \"validation\" (or \"holdout\") group.\n\nHO_rows   <- sample(nrow(train.mtx), size = tr_num)\nx_train   <- train.mtx[ HO_rows, -1]\nx_holdout <- train.mtx[-HO_rows, -1]\ny_train   <- train.mtx[ HO_rows,  1] %>% to_categorical(10)\ny_holdout <- train.mtx[-HO_rows,  1] %>% to_categorical(10)\n\n# Reshape the training data into 28x28 array, instead of 1x784 columns.\nx_train   <- array_reshape(x_train,   c(nrow(x_train),   Im_R, Im_C, 1))\nx_holdout <- array_reshape(x_holdout, c(nrow(x_holdout), Im_R, Im_C, 1))\nx_test    <- array_reshape(test.mtx,  c(nrow(test.mtx),  Im_R, Im_C, 1))\n\ndim(x_train)      # each observation should now take the form of a 28x28 matrix.\n","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:34.673501Z","iopub.execute_input":"2022-12-10T19:49:34.674604Z","iopub.status.idle":"2022-12-10T19:49:37.316403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This is the \"guts\" of the model. Obviously we need to flesh out what all these lines of code are doing.\n\nmodel <- keras_model_sequential() \nmodel %>% \n  layer_conv_2d(filters = 64, kernel_size = c(3,3), activation = 'relu', input_shape = c(Im_R, Im_C, 1)) %>% \n  layer_max_pooling_2d(pool_size = c(2, 2)) %>% \n  layer_dropout(rate = .4) %>% \n  layer_flatten() %>% \n  layer_dense(units = 128, activation = 'relu') %>% \n  layer_dropout(rate = .4) %>% \n  layer_dense(units = 10, activation = 'softmax')\n\nsummary(model)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:37.318130Z","iopub.execute_input":"2022-12-10T19:49:37.319092Z","iopub.status.idle":"2022-12-10T19:49:37.421767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compile the model. The \"history\" feature lets us plot the evolution of the error and loss as the model \"learns.\"\n\nmodel %>% compile(loss='categorical_crossentropy', \n                  optimizer=optimizer_adam(learning_rate = 0.002),\n                  metrics=c('accuracy')\n                 )\n\nhistory <- model %>% fit(x_train, y_train, epochs=15, batch_size=32, validation_split=0.2)\n\nplot(history)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:49:37.423523Z","iopub.execute_input":"2022-12-10T19:49:37.424590Z","iopub.status.idle":"2022-12-10T19:53:41.593392Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Evaluating results: confusion matrix, precision, recall, f1","metadata":{}},{"cell_type":"code","source":"model %>% evaluate(x_train, y_train)\nmodel %>% evaluate(x_holdout, y_holdout)\n\n# Confusion matrix on holdout data set\nholdout_pred <- model %>% predict(x_holdout) %>% k_argmax() %>% as.vector()\nholdout_lab  <- y_holdout %>% k_argmax() %>% as.vector()\n\nconfusionMatrix(as.factor(holdout_pred), as.factor(holdout_lab))","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:53:41.596295Z","iopub.execute_input":"2022-12-10T19:53:41.597317Z","iopub.status.idle":"2022-12-10T19:53:50.198638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here is the confusion matrix plotted as a tile plot in ggplot2.","metadata":{}},{"cell_type":"code","source":"label      <- as.factor(holdout_lab)\npredicted  <- as.factor(holdout_pred)\ncM         <- data.frame(label, predicted) %>% group_by(label, predicted) %>% summarise(count = n()) \n\nggplot(cM, aes(x = label, y = ordered(predicted, seq(9,0)), fill = count)) +\n  geom_tile(color = \"black\") + geom_text(aes(label = count)) +\n  scale_fill_gradient(low = \"white\", high = \"red\") +\n  ylab(\"Predicted Class\") + \n  xlab(\"True label\") + \n  ggtitle(\"Confusion Matrix for Holdout Data\")\n\n","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:53:50.200817Z","iopub.execute_input":"2022-12-10T19:53:50.202361Z","iopub.status.idle":"2022-12-10T19:53:50.481600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"c     <- as.matrix(table(Actual = holdout_lab, Predicted = holdout_pred))\nn     <- sum(c)\nr_sum <- apply(c, 1, sum) # number of instances per class\nc_sum <- apply(c, 2, sum) # number of predictions per class\ndiag  <- diag(c)\n\naccuracy  <- sum(diag)/n\nprecision <- diag / c_sum   # precision by class\nrecall    <- diag / r_sum   # recall by class\nf1        <- 2 * precision * recall / (precision + recall)\n\naccuracy\ndata.frame(precision, recall, f1) ","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:53:50.483413Z","iopub.execute_input":"2022-12-10T19:53:50.484429Z","iopub.status.idle":"2022-12-10T19:53:50.529797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, here is a rogue's gallery of images in the holdout set that were predicted incorrectly. The caption in the top left of each cell tells you the true label, followed by our prediction. For some of these (e.g. top left \"9\") you can clearly see both why the prediction was wrong, and why it guessed incorrectly. For others (e.g. the \"9\" in the second row) the value looks ambiguous anyway - is that a 9 or a 4? And some of these images are almost impossible to parse.","metadata":{}},{"cell_type":"code","source":"label <- y_holdout %>% k_argmax() %>% as.vector()     # labels for holdout set\npred  <- holdout_pred                                 # predictions for holdout set\nt_row <- setdiff(seq(1,42000), HO_rows)               # row index in train for holdout set\ndf <- data.frame(t_row,label,pred) %>% filter(label != pred)\n\npar(mfrow=c(5,5), pty='s', mai=c(0.0, 0, 0, 0.0))     # plot the first 25 errors\nfor (x in 1:25) {\n    t <- train[df[x,1], -1]\n    draw(t)\n    text(0.2,0.9, paste(df[x,2],\"/\",df[x,3]), cex=2)\n    box(lty = 'solid')\n}","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:53:50.531605Z","iopub.execute_input":"2022-12-10T19:53:50.532629Z","iopub.status.idle":"2022-12-10T19:53:50.780453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Kaggle Competition Submission","metadata":{}},{"cell_type":"code","source":"# Run our final model to generate predictions for the 28,000 test observations.\n# Store the results to a data frame that we can save, download and submit for scoring.\n\ntest_predictions <- model %>% predict(x_test) %>% k_argmax() %>% as.vector()\npreds_df <- data.frame(ImageId=c(1:nrow(test)), Label=test_predictions)\nhead(preds_df)\nwrite.csv(preds_df, \"MNIST_cjp_v1.csv\", row.names=FALSE)","metadata":{"execution":{"iopub.status.busy":"2022-12-10T19:53:50.783663Z","iopub.execute_input":"2022-12-10T19:53:50.784883Z","iopub.status.idle":"2022-12-10T19:53:55.585705Z"},"trusted":true},"execution_count":null,"outputs":[]}]}