{ "cells": [ { "cell_type": "markdown", "metadata": { "_cell_guid": "6ceb4651-642e-73f2-0fd9-0d9f86468ae2" }, "source": [ "## Introduction ##\n", "\n", "\n", "----------\n", "\n", "\n", "\n", "*The work is based on dataset from Digit Recognition Competition. In the next lines I will 1) draw digits from dataset; 2) train multinomial logistic regression with gradient descent using my handmade script; 3) calculte model performance; 4) draw heatmap of model coefficients; 5) draw observations predicted correctly vs falsely. This work gives understanding how packages and functions calculated model coeffificent.*\n", "\n", "\n", "----------\n", "\n", "\n", "Have you ever had this uncomfortable feeling when you don't understand what is going on? I have. Actually I have it every time when I use **glm** function. I print a line of code, wait a couple of minutes - and that's it. Some magic is going on and I am planning to find out what exactly.\n", "\n", "My plan is the following: to build a multinomial logistic regression with gradient descent using just simple functions and a lot of loops.\n", "\n", "Let's see what I've got!\n", "\n", "CAUTION! The code below may contain examples of bad design. Do not try to repeat it yourself. \n", "\n", "(I have never worried about how my code looks like as long as it works. But I know my code will be read by PEOPLE, and people have feelings, rather than these soulless peaces of glass, plastic and metal called computers. So I've done my best not to make you cry.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "8774f395-7591-0e0c-828f-2d8091497108" }, "outputs": [], "source": [ "#Loading data\n", "train <- read.csv(\"../input/train.csv\")\n", "Y <- train[,1]\n", "X <- train[,-1]\n", "X_ones <- cbind(ones=rep(1,42000),X)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "029cd196-1c10-d058-2bee-a3d68b19da5e" }, "source": [ "## Data Visualization ##\n", "\n", "\n", "----------\n", "\n", "\n", "**Rule #1. First of all find out how your data looks like.**\n", "\n", "When your dataset is just a big table of numbers that used to be pictures, it is kind of hard. But I really really wanted to see these handwritten digits that I am going to recognize. So I wrote a script.\n", "\n", "The idea is the following:\n", "\n", "1) choose which rows from dataset you want to see and store them in variable **data**; \n", "\n", "2) extract these rows one by one, convert them to matrix (actually, I converted data to matrix by columns, but not by rows as it should be done, because of function **image**);\n", "\n", "3) use function to visualize values in matrix - **image** is what I need. But it rotates matrix (why?), so I rotated it too before putting it in function. Now picture looks appropriately.\n", "\n", "Yeah, that's it! I am ready for the next step." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "c0bd561d-88d3-6d8e-30ba-5692e0d0dabb" }, "outputs": [], "source": [ "data <- sample(1:42000,100)\n", "par(mfrow=c(10,10),mar=c(0.1,0.1,0.1,0.1))\n", "\n", "for (k in data)\n", "{\n", " row <- NULL\n", " for (n in 1:784)\n", " row[n] <- X[k,n]\n", " \n", " matrix1 <- matrix(row,28,28,byrow=FALSE)\n", " matrix2 <- matrix(rep(0,784),28,28)\n", "\n", " for (i in 1:28)\n", " for (j in 1:28)\n", " matrix2[i,28-j+1] <- matrix1[i,j]\n", " \n", " image(matrix2, axes=FALSE, col=topo.colors(12))\n", "}\n", "\n", "rm(i,n,j,k,row,matrix1,matrix2)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "60eb67d3-cbf0-2a28-f0b6-f6feb96af9ae" }, "source": [ "## Splitting Data into Train and Test Sets ##\n", "\n", "\n", "----------\n", "\n", "\n", "**Rule #2. The more data is used to train a model, the better the model.**\n", "\n", "It is a good practice to check how model performs before using it on the set for classification. So I leave 10,000 observations for testing (which is about 25% of total dataset). \n", "\n", "I also don't want model coefficients to vary when I recalculate them, so I use the same observations in train and test sets each time. Each set has sufficient amount of observations from all classes, there is no asymmenty. \n", "\n", "Successful split!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "5248797b-3be7-2b21-8fff-2e056be58cb6" }, "outputs": [], "source": [ "set.seed(1234)\n", "intest <- sample(1:42000,10000)\n", "X_train <- X_ones[-intest,]\n", "Y_train <- Y[-intest]\n", "X_test <- X_ones[intest,]\n", "Y_test <- Y[intest]\n", "\n", "table(Y_train)\n", "table(Y_test)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "c1434ed9-6f30-90b5-bbf4-664987085e07" }, "source": [ "## Data Normalization ##\n", "**Rule #3. Gradient descent works badly if data is not normalized.**\n", "\n", "The process is simple: for each observation subtract its column mean and divide by its column standard deviation. However, for train set means and sds calculated as they are, but for test set - means and sds from train set should be used.\n", "\n", "(Yeah-yeah! There is no need to create variables **mean_new** and **sd_new**, but recalculation of this variables in line 12 will slow down script execution, because this line has to be executes 32,000 * 875 = 25,120,000 times. So while working with big data set, you'd better add some extra variables.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "73cb08c1-727c-1910-6d50-9f7420d0599f" }, "outputs": [], "source": [ "#Normalizing data in Train set\n", "means <- 0\n", "sd <- 1 \n", "for (i in 2:785)\n", "{\n", " if (mean(X_train[,i]) != 0)\n", " {\n", " mean_new <- mean(X_train[,i])\n", " sd_new <- sd(X_train[,i])\n", " means <- c(means,mean_new)\n", " sd <- c(sd,sd_new)\n", " X_train[,i] <- (X_train[,i]-mean_new)/sd_new #line 12\n", " }\n", " else\n", " {\n", " means <- c(means,0)\n", " sd <- c(sd,1)\n", " }\n", "}\n", "\n", "#Normalizing data in Test set\n", "z <- dim(X_test)[2]\n", "for (i in 2:z)\n", "{\n", " X_test[,i] <- (X_test[,i]-means[i])/sd[i]\n", "}\n", "\n", "rm(mean_new,sd_new,i,z)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "e8b4acbc-58c7-043c-b460-ebc7cb073d07" }, "source": [ "## Model Training ##\n", "\n", "\n", "----------\n", "\n", "\n", "**Rule #4. Training model on big data sets takes time. Spend this time with pleasure.**\n", "\n", "I assume, you know that building multinomial logit with **n** classes (n>2) requires building **n** different logistic regressions - for each class. That gradients descent is iterative method which updates coefficients on each step. I belive, that you've already seen cost and gradient functions for logit. And, of course, you know what is the meaning of life.\n", "\n", "So, no details here. I will do the following:\n", "\n", "1. Train model for each class - 10 models in total.\n", "\n", "2. For each observation, predict the probability from all 10 models. So I will get the probability that particular observation is 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9.\n", "\n", "3. Choose model with highest probability and assign observation to its class.\n", "\n", "But first of all, useful handmade functions:\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "667fe72d-6865-fa39-2572-21840acf4d8a" }, "outputs": [], "source": [ "#sigmoid function\n", "sigmoid <- function(z)\n", "{\n", " g<-1/(1+exp(-z))\n", " return(g)\n", "}\n", "\n", "#cost function\n", "cost <- function(x,y,theta)\n", "{\n", " m <- dim(x)[1]\n", " J <- -1/m * sum(y*log(sigmoid(x %*% theta)) + (1 - y)*log(1 - sigmoid(x %*% theta))) \n", " return(J)\n", "}\n", "\n", "#gradient\n", "grad <- function(x,y,theta)\n", "{\n", " grad <- t(sigmoid(x %*% theta) - y) %*% x \n", " grad <- t(grad)\n", " return(grad)\n", "}\n", "\n", "#model valuation\n", "accuracy <- function (pred,actual)\n", "{\n", " conf_matrix <- table(predictions = pred, actual = actual)\n", " sum <- 0\n", " for (j in 1:10)\n", " for (i in 1:10)\n", " if (j == i)\n", " sum<-sum+conf_matrix[i,j]\n", " return(list(conf_matrix,sum/length(pred)))\n", "}" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "583d36e3-2feb-1641-72ac-63ed5451d0fc" }, "source": [ "And, of course, here is script for gradient descent. \n", "\n", "\n", "What it does: calculates all thetas for 10 models and stores them in variable **all_thetas**; predicts probabilities from 10 models and stores them in variable **pred_train**; saves output of cost function after each iteration for all 10 models in variable **all_costs** (because cost behavior tells whether gradient descent works OK)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "af05157f-5a64-5d37-2e44-7ad462b84171" }, "outputs": [], "source": [ "all_thetas <- NULL\n", "all_costs <- NULL\n", "X_train <- as.matrix(X_train)\n", "pred_train <- Y_train\n", "alpha <- 0.00002\n", "iterations <- 200\n", "\n", "for (j in 1:10) \n", "{\n", " theta <- matrix(rep(0,dim(X_train)[2]),dim(X_train)[2],1)\n", " costs <- as.data.frame(matrix(rep(0,iterations),iterations,2))\n", " n <- j-1\n", " Y <- as.matrix(ifelse(Y_train == n, 1, 0))\n", "\n", " for (i in 1:iterations)\n", " {\n", " gradient <- grad(X_train,Y,theta)\n", " theta <- theta - alpha*gradient\n", " costs[i,2] <- cost(X_train,Y,theta)\n", " costs[i,1] <- i\n", " }\n", " \n", " all_costs <- cbind(all_costs,costs[,2])\n", " all_thetas <- cbind(all_thetas,theta)\n", " pred_new <- round(sigmoid(X_train %*% theta), digits = 2)\n", " pred_train <- cbind(pred_train, pred_new)\n", "}\n", "\n", "rm(costs,gradient,theta,pred_new,Y,alpha,i,iterations,j,n)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "9e67a152-2177-a6fc-4b50-bdd635d4e0b5" }, "source": [ "Here are graphs of costs vs interations. If gradient descent works ok, costs decrease on each step. But if you see something strange:\n", "\n", "1) costs still rapidly decrease near the 200th iteration - try higher alpha or increase number of iterations;\n", "\n", "2) all costs are **NaN** or **Inf** - decrease alpha.\n", "\n", "My costs look good. That's the result of many trials with different number of iterations and value of alpha." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "98ec7709-588c-1119-3ca6-e1dbb70c4328" }, "outputs": [], "source": [ "all_costs <- as.data.frame(all_costs)\n", "par(mfrow=c(3,4),mar=c(1.5,2.5,2.5,1))\n", "\n", "for (j in 1:10)\n", "{\n", " plot(x=rownames(all_costs),y=all_costs[,j],main=c(j-1),ylab=\"Cost Function\",xlab=\"Iterations\")\n", "}" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "b18a800f-6629-9201-7c39-6b496cf74d93" }, "source": [ "Cool! Move on." ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "999c2448-7a27-8dc9-6929-6cdab50e50d2" }, "source": [ "## Model Testing ##\n", "\n", "----------\n", "**Rule #5. Nobody cares about your train set errors, except your Mom.**\n", "\n", "It is time to find out how good my model is. Well, train set accuracy is 91.5% (Mom, that's for you). Test set accuracy - 90.0%.\n", "\n", "(Yeah-yeah, there is a lot of place for improvements. I am accepting ideas!)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "bcbe65c9-59fd-1062-e104-334950f00b25" }, "outputs": [], "source": [ "pred_train<-as.data.frame(pred_train)\n", "all_thetas<-as.data.frame(all_thetas)\n", "\n", "#Train set predictions\n", "pred_train$predictions <- NA\n", "pred_train$max_prob <- NA\n", "m<-dim(pred_train)[1]\n", "\n", "for (i in 1:m)\n", "{\n", " pred_train$predictions[i] <- which.max(pred_train[i,2:11])[[1]] - 1\n", " pred_train$max_prob[i] <- max(pred_train[i,2:11])[[1]]\n", "}\n", "\"train set\"\n", "accuracy(pred_train$predictions,pred_train$pred_train)\n", "\n", "\n", "#Test set predictions\n", "pred_test <- Y_test\n", "X_test <- as.matrix(X_test)\n", "\n", "for (i in 1:10)\n", "{\n", " pred_new_t <- round(sigmoid(X_test %*% as.matrix(all_thetas[,i])), digits = 2)\n", " pred_test <- cbind(pred_test,pred_new_t)\n", "}\n", "\n", "pred_test <- as.data.frame(pred_test)\n", "pred_test$predictions <- NA\n", "pred_test$max_prob <- NA\n", "m<-dim(pred_test)[1]\n", "\n", "for (i in 1:m)\n", "{\n", " pred_test$predictions[i] <- which.max(pred_test[i,2:11])[[1]] - 1\n", " pred_test$max_prob[i] <- max(pred_test[i,2:11])[[1]]\n", "}\n", "\n", "rm(i,pred_new_t,m)\n", "\"test set\"\n", "accuracy(pred_test$predictions,pred_test$pred_test)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "4f410d04-7fef-0711-2fc8-8f53fe4f489a" }, "source": [ "## Visualization of Results ##\n", "\n", "\n", "----------\n", "**Rule #6. Graphical information is easier for human perception.**\n", "\n", "I wanted to find out how thetas look like, maybe it will give me a clue why model does not predict 100% of observations. So I drew a heatmap of thetas from each model. The script is similar to that at the begginning (when I drew digits). \n", "\n", "Just a minute, I am going to explain how to read it. Close look at this formula:\n", "\n", " predictions <- 1/(1+exp(- X_train %*% theta))\n", "\n", "It means that when **X_train %*% theta** is more than 0, we will get probability more than 0.5. In other words, if **X_train %*% theta** more then 0, observation will be classified as being from the class that this model predicts, otherwise, it is from other class. (To understand this idea, I have looked at graph of sigmoid function for 30 minutes. Recommended!).\n", "\n", "Next. Variables (darkness of the pixels) are zeros or positive numbers. Thetas - either negative or positive number. Colors on the heatmap: red - for low (negative) values, yellow - for average, white - for high. Now you have all the information.\n", "\n", "If color of theta-pixel on the theta heatmap is more *red* than the background (these zero pixels in the corners) and observation has a nonzero value of this variable-pixel, there will be *penaly*, so term **X_train %*% theta** *decreases* and probability of being from this class also *decreases*.\n", " \n", "Read one more time, putting \"yellow\", \"reward\", \"increases\" and \"increases\" instead of \"red\", \"penalty\", \"decreases\" and \"decreases\".\n", "\n", "Let's consider heatmap of thetas from model for predicting ones (the most obvious one). It seems like for all observations, that have nonzero pixels in the middle and zero pixels around, we expect model to predict this observation to be one. (Because model-1 will give higher probability, then model-7, model-0 or model-5)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "61fa6100-6b8f-85ba-ea19-243fe8a59989" }, "outputs": [], "source": [ "par(mfrow=c(3,4),mar=c(0.1,0.1,2,0.1))\n", "\n", "for (k in 1:10)\n", "{\n", " matrix1 <- matrix(t(all_thetas[-1,k]),28,28,byrow=FALSE)\n", " matrix2 <- matrix(rep(0,784),28,28)\n", " \n", " for (i in 1:28)\n", " for (j in 1:28)\n", " matrix2[i,28-j+1] <- matrix1[i,j]\n", " \n", " image(matrix2, axes=FALSE, col=heat.colors(12),main = k-1)\n", "}\n", "\n", "rm(matrix1,matrix2,k,i,j)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "e63fee72-6a17-3539-6afd-ab82556d549e" }, "source": [ "Here are the set of observation from test set, predicted correctly (blue) and incorrectly (red). You can compare them with thetas heatmaps and thinking about reasons, why they were predicted in this way." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "30609e81-6886-f0ba-638d-01164c37141e" }, "outputs": [], "source": [ "pred_test$false <- ifelse(pred_test$pred_test == pred_test$predictions,0,1)\n", "par(mfrow=c(14,14),mar=c(0.1,0.1,0.1,0.1))\n", "data<-c(sample(as.integer(row.names(pred_test[pred_test$false == 0,])),98), \n", " sample(as.integer(row.names(pred_test[pred_test$false == 1,])),98))\n", "\n", "z<-0\n", "for (k in data)\n", "{\n", " row <- NULL\n", " for (n in 1:784)\n", " row[n] <- X[k,n]\n", " \n", " matrix1 <- matrix(row,28,28,byrow=FALSE)\n", " matrix2 <- matrix(rep(0,784),28,28)\n", "\n", " for (i in 1:28)\n", " for (j in 1:28)\n", " matrix2[i,28-j+1] <- matrix1[i,j]\n", " \n", " z<-z+1\n", " if (z