{"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":"Following along with Samson Zhang's neural network from scratch tutorial (but doing it all in R instead of Python).\n\nCheck out his YouTube walkthrough (https://www.youtube.com/watch?v=w8yWXqWQYmU) and python notebook (https://www.kaggle.com/code/wwsalmon/simple-mnist-nn-from-scratch-numpy-no-tf-keras/notebook)","metadata":{}},{"cell_type":"code","source":"#load data\ndata <- read.csv(\"../input/digit-recognizer/train.csv\")","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2022-08-11T14:35:04.606003Z","iopub.execute_input":"2022-08-11T14:35:04.607891Z","iopub.status.idle":"2022-08-11T14:35:13.074574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#turn into a matrix and get number of rows and cols\ndata <- as.matrix(data)\nm = nrow(data)\nn = ncol(data)\n\n#shuffle rows\ndata_shuffled <- data[sample(nrow(data)), ]\n\n#test-train split\ndata_dev <- t(data_shuffled[0:1000, ])\nY_dev <- data_dev[1,]\nX_dev <- data_dev[2:n, ]\n\ndata_train <- t(data_shuffled[1000:m, ])\nY_train <- data_train[1,]\nX_train <- data_train[2:n, ]\n\n#scale pixel data (to range from 0 to 1)\nX_dev <- X_dev / 255\nX_train <- X_train / 255","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2022-08-11T14:35:13.077427Z","iopub.execute_input":"2022-08-11T14:35:13.079142Z","iopub.status.idle":"2022-08-11T14:35:14.528898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#initialize parameters\ninit_params <- function(){\n  W1 <- matrix(runif(10 * 784, -.5, .5) , nrow = 10, ncol = 784)\n  b1 <- matrix(runif(10, -.5, .5), nrow = 10, ncol = 1)\n  W2 <- matrix(runif(10 * 10, -.5, .5) , nrow = 10, ncol = 10)\n  b2 <- matrix(runif(10, -.5, .5), nrow = 10, ncol = 1)\n  \n  return(list(W1 = W1,\n              b1 = b1,\n              W2 = W2,\n              b2 = b2))\n}\n\n#ReLU function\nReLU <- function(x){\n  return(max(0, x))\n}\n\n#softmax function\nsoftmax <- function(x){\n  stopifnot(is.matrix(x))\n  \n  temp <- x\n  for(i in 1:nrow(x)){\n    for(j in 1:ncol(x)){\n      temp[i,j] <- exp(x[i, j]) / sum(exp(x[, j]))\n    }\n  }\n  \n  return(temp)\n}\n\n#forward propagation\nforward_prop <- function(params, X){\n  Z1 <- params$W1 %*% X + rep(params$b1, each = ncol(X))\n  A1 <- apply(Z1, c(1, 2), ReLU)\n  Z2 <- params$W2 %*% Z1 + rep(params$b2, each = ncol(X))\n  A2 <- softmax(Z2)\n  \n  return(list(Z1 = Z1,\n              A1 = A1,\n              Z2 = Z2,\n              A2 = A2))\n}\n\n#one hot encode function\none_hot <- function(Y){\n  temp <- matrix(0, 10, length(Y))\n  for(i in 1:length(Y)){\n    temp[Y[i] + 1, i] <- 1 #must add plus one because the lowest Y digit is 0 instead of 1\n  }\n  row.names(temp) <- 0:9\n  return(temp)\n}\n\n#derivative of ReLU\nderiv_ReLU <- function(Z){\n  return(+(Z > 0))\n}\n\n#backward propagation\nback_prop <- function(fparams, params, X, Y){\n  m = length(Y)\n  one_hot_Y <- one_hot(Y)\n  dZ2 = fparams$A2 - one_hot_Y\n  dW2 = 1 / m * dZ2 %*% t(fparams$A1)\n  db2 = 1 / m * rowSums(dZ2)\n  dZ1 = params$W2 %*% dZ2 * deriv_ReLU(fparams$Z1)\n  dW1 = 1 / m * dZ1 %*% t(X)\n  db1 = 1 / m * rowSums(dZ1)\n  \n  return(list(dW1 = dW1,\n              db1 = db1,\n              dW2 = dW2,\n              db2 = db2))\n}\n\n#update parameters\nupdate_params <- function(params, bparams, alpha){\n  W1 = params$W1 - alpha * bparams$dW1\n  b1 = params$b1 - alpha * bparams$db1\n  W2 = params$W2 - alpha * bparams$dW2\n  b2 = params$b2 - alpha * bparams$db2\n  \n  return(list(W1 = W1,\n              b1 = b1,\n              W2 = W2,\n              b2 = b2))\n}\n\n#get model prediction\nget_predictions <- function(fparams){\n  pred <- numeric()\n  for(i in 1:ncol(fparams$A2)){\n    pred[i] <- which(fparams$A2[,i] == max(fparams$A2[,i])) - 1\n  }\n  return(pred)\n}\n\n#get accuracy\nget_accuracy <- function(predictions, Y){\n  stopifnot(length(predictions) == length(Y))\n  return(sum(+(predictions == Y)) / length(Y))\n}\n\n#gradient descent\ngradient_descent <- function(X, Y, iterations, alpha, X_test, Y_test){\n  params <- init_params()\n  W1 <- params$W1\n  b1 <- params$b1\n  W2 <- params$W2\n  b2 <- params$b2\n  \n  for(i in 1:iterations){\n    fparams <- forward_prop(params, X)\n    bparams <- back_prop(fparams, params, X, Y)\n    params <- update_params(params, bparams, alpha)\n    \n    if(i %% 10 == 0){ #run this every 10 iterations\n      print(paste(\"Iteration: \", i))\n      print(paste(\"Train Accuracy: \", round(get_accuracy(get_predictions(fparams), Y), 3)))\n      print(paste(\"Test Accuracy: \", round(get_accuracy(get_predictions(forward_prop(params, X_test)), Y_test), 3)))\n    }\n  }\n  \n  return(params)\n}","metadata":{"execution":{"iopub.status.busy":"2022-08-11T14:35:14.533430Z","iopub.execute_input":"2022-08-11T14:35:14.535307Z","iopub.status.idle":"2022-08-11T14:35:14.571722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#run gradient descent\nparams <- gradient_descent(X_train, Y_train, 500, .1, X_dev, Y_dev)","metadata":{"execution":{"iopub.status.busy":"2022-08-11T14:35:14.575167Z","iopub.execute_input":"2022-08-11T14:35:14.577460Z"},"trusted":true},"execution_count":null,"outputs":[]}]}