{"cells":[{"metadata":{},"cell_type":"markdown","source":"Objective: Develop a multiclass classifier using the entire MNIST data set for input data. The classifier will be used to predict which of 10 digits is being written, evaluated with and without dimension reduction with principal components analysis (PCA). \n\n**Kaggle ID: Claire Boetticher**<br>\n**Kaggle username: clairence**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import dependencies\nimport numpy as np\nimport pandas as pd\nimport os\nimport random\nimport time\n\n# ML imports\nfrom sklearn.model_selection import KFold, train_test_split, cross_val_score, cross_val_predict, GridSearchCV \nfrom sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score, classification_report\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.decomposition import PCA, IncrementalPCA\n\n# Visualizations\nimport matplotlib.pyplot as plt  # static plotting\nimport seaborn as sns  # pretty plotting, including heat map\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Initialize process time list for code timing across entire program\nprocess_time = []","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Import and inspect training data"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import train and test datasets\n\n# The training data set, (train.csv), has 785 columns. The first column, called \"label\", is the \n# digit that was drawn by the user. The rest of the columns contain the pixel-values of the associated image.\nsample_submission = pd.read_csv(\"../input/digit-recognizer/sample_submission.csv\")\ntest = pd.read_csv(\"../input/digit-recognizer/test.csv\")\ntrain = pd.read_csv(\"../input/digit-recognizer/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 42,000 images\n# Each image has 784 features (each image is 28 x 28 pixels) and 1 label feature (785 total)\n# Each feature represents one pixel's intensity from 0 (white) to 255 (black)\ntrain.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 28,000 images, contains 784 features and no label\ntest.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Inspect class balances for train set, seems relatively balanced\nprint(train['label'].value_counts(ascending=False))\nprint('--------------------')\nprint(train['label'].value_counts(normalize=True))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Figure 1: Label Distribution in Train Set**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot total distribution of labels in Kaggle train set\nmn_plt_total = sns.countplot(train['label'], palette=\"muted\").set_title('Total Digit Distribution')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Save the labels to a Pandas series target\ny = train['label']\n\n# Drop the label feature\nX = train.drop(\"label\",axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Confirm that all 256 values between the min-max of 0-255 exist in the train set\nlen(np.unique(X))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# View images\nimages_to_plot = 9\nrandom_indices = random.sample(range(42000), images_to_plot)\n\nsample_images = X.loc[random_indices, :]\nsample_labels = y.loc[random_indices]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot examples\nplt.clf()\nplt.style.use('seaborn-muted')\n\nfig, axes = plt.subplots(3,3, \n                         figsize=(5,5),\n                         sharex=True, sharey=True,\n                         subplot_kw=dict(adjustable='box', aspect='equal')) #https://stackoverflow.com/q/44703433/1870832\n\nfor i in range(images_to_plot):\n    \n    # axes (subplot) objects are stored in 2d array, accessed with axes[row,col]\n    subplot_row = i//3 \n    subplot_col = i%3  \n    ax = axes[subplot_row, subplot_col]\n\n    # plot image on subplot\n    plottable_image = np.reshape(sample_images.iloc[i,:].values, (28,28))\n    ax.imshow(plottable_image, cmap='gray_r')\n    \n    ax.set_title('Digit Label: {}'.format(sample_labels.iloc[i]))\n    ax.set_xbound([0,28])\n\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step 1: Fit a Random Forest classifier and evaluate on test set\n- OOB model\n- Train and evaluate\n- Grid search\n- Train and evaluate\n- Predict\n- Confusion matrix"},{"metadata":{},"cell_type":"markdown","source":"### 1. Train test split and distribution of labels"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Assign classifier\n# Use default criterion (gini)\ncls_rf = RandomForestClassifier(random_state=1, n_jobs=-1, n_estimators=90)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Split train and test set 90/10\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.20, random_state=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check split\nprint(X_train.shape)\nprint(y_train.shape)\nprint(X_test.shape)\nprint(y_test.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Figure 2: Label Distribution in Training Set**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Evaluate distribution of labels across train set\nmn_plt_trn = sns.countplot(y_train, palette=\"muted\").set_title('Train Digit')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Figure 3: Label Distribution in Test Set**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Evaluate distribution of labels across test set\nmn_plt_trn = sns.countplot(y_test, palette=\"muted\").set_title('Test Digit')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. Initial classification with Random Forest"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fit model and run on test set with timing\nstart_time = time.time()\ncls_rf.fit(X_train, y_train)\nprint('Accuracy: ', cls_rf.score(X_test, y_test))\n\n# Timing\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3. Grid search for fine-tuning\nTry n_estimator hyperparameter"},{"metadata":{"trusted":true},"cell_type":"code","source":"# cls_rf = RandomForestClassifier()\n\nparam_grid = {\n    'bootstrap': [True],\n    'n_estimators': [100,150,200,250,300,400]\n}\n\ngrid_search_rf = GridSearchCV(estimator=cls_rf, param_grid=param_grid, cv=3, \n                              return_train_score=True, n_jobs=-1, verbose=2)\n\ngrid_search_rf.fit(X_train, y_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Display parameter recommendations\nprint('Test set score: ', grid_search_rf.score(X_test, y_test))\nprint('Best parameters: ', grid_search_rf.best_params_)\nprint('Best cross-validation score: ', grid_search_rf.best_score_)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Adjust model parameters, not using max depth\ncls_rf = RandomForestClassifier(random_state=1, n_jobs=-1, n_estimators=300)\n\n# Original\n# cls_rf = RandomForestClassifier(random_state=1, n_jobs=-1, n_estimators=90)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fit the tuned parameters with timing\nstart_time = time.time()\ncls_rf.fit(X_train, y_train)\n\n# Results with tuning\nprint('Accuracy: ', cls_rf.score(X_test, y_test))\n\n# Timing\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# F1 score\npred = cls_rf.predict(X_test)\nprint('F1 accuracy: ', f1_score(pred, y_test, average='macro'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Print classification report\ntarget_names = ['class 0', 'class 1', 'class 2', 'class 3', 'class 4', 'class 5', 'class 6', 'class 7',\n                   'class 8', 'class 9']\nprint(classification_report(y_test, pred, target_names=target_names))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3. Confusion matrix"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot confusion matrix of actual versus predicted labels\nrf_cm = confusion_matrix(y_test, pred)\nrf_cm_plt=sns.heatmap(rf_cm.T, square=True, annot=True, fmt='d', cbar=False, cmap=\"Blues\")\nplt.xlabel('Actual label')\nplt.ylabel('Predicted label')\nplt.title(\"Valid\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4. Evaluate model on test set for Kaggle"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Submission code\ntestData = pd.read_csv(\"test.csv\")\nstart_time = time.time()\npred = cls_rf.predict(testData)\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)\n\n# # Create Dataframe\n# data = pred\n# df_1 = pd.DataFrame(pred)\n# df_1['ImageID'] = df_1.index + 1\n# df_1.columns = ['Label', 'ImageID']\n# submission = df_1[['ImageID', 'Label']]\n\n# # Output to csv\n# submission.to_csv('Boetticher_RF_predictions.csv',header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step 2: Principal Components Analysis"},{"metadata":{},"cell_type":"markdown","source":"### 1. Data preparation and PC identification"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Combine training and test set for PCA\nx = np.concatenate((X_train, X_test), axis=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fit model using PCA, generating principal components that represent 95 percent of the variability in \n# the explanatory features\nstart_time = time.time()\npca = PCA(.95)\npca.fit(x)\ntotimages = pca.transform(x)\n# pca.n_components_\n\n# Timing\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Output number of principal components explaining 95% of variability in the explanatory features: 154\nprint('Principal components count: ', pca.n_components_)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. Explained variance ratio and visualization"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Explained variance ratio for 154 PCs\n# 10% of the train dataset's variance lies along the first PC\n# 7% along the second\n# 6% along the third, etc.\npca.explained_variance_ratio_","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Explained variance plot\nplt.plot(np.cumsum(pca.explained_variance_ratio_))\nplt.xlabel('number of components')\nplt.ylabel('cumulative explained variance')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Reform the training and testing images\n# Training: up to 33,600\n# Testing: 33,600 to 42,000\ntrainimages = totimages[0:33600, :]\ntestimages = totimages[33600:42000, :]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Convert to integers\ntrainimages = trainimages.astype(int)\ntestimages = testimages.astype(int)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3. Test Incremental PCA"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Time component identification\nstart_time = time.time()\n\nn_batches = 100\ninc_pca = IncrementalPCA(n_components=154)\nfor X_batch in np.array_split(X_train, n_batches):\n    inc_pca.partial_fit(X_batch)\n\nX_reduced = inc_pca.transform(X_train)\n\n# Timing\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Explained variance ratio for 154 PCs with incremental PCA\n# 10% of the train dataset's variance lies along the first PC\n# 7% along the second\n# 6% along the third, etc.\n# Not drastically different than PCA and took longer\ninc_pca.explained_variance_ratio_","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step 3: Random Forest Classifier with Principal Components"},{"metadata":{},"cell_type":"markdown","source":"### 1. Fit model with principal components"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Fit the model\nstart_time = time.time()\ncls_rf2 = RandomForestClassifier(random_state=1, n_jobs=-1, n_estimators=300)\ncls_rf2.fit(trainimages, y_train.values.ravel())\nprint('Accuracy: ', cls_rf2.score(testimages, y_test))\n\n# Timing\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Predict a second time with principal components\npred2 = cls_rf2.predict(testimages)\nprint('F1 accuracy: ', f1_score(pred2, y_test, average='macro'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. Confusion matrix"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Plot confusion matrix of actual versus predicted labels\nrf_cm = confusion_matrix(y_test, pred2)\nrf_cm_plt=sns.heatmap(rf_cm.T, square=True, annot=True, fmt='d', cbar=False, cmap=\"Blues\")\nplt.xlabel('Actual label')\nplt.ylabel('Predicted label')\nplt.title(\"Valid\");","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3. Evaluate model on test set for Kaggle"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Submission code with timing\nx_test = pd.read_csv(\"test.csv\")\nstart_time = time.time()\npred_PCA = cls_rf2.predict(pca.transform(x_test))\nelapsed_time = time.time() - start_time\nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\nprint('Time in seconds: ', elapsed_time)\nprocess_time.append(elapsed_time)\n\n# # Create Dataframe\n# data = pred_PCA\n# df_1 = pd.DataFrame(pred_PCA)\n# df_1['ImageID'] = df_1.index + 1\n# df_1.columns = ['Label', 'ImageID']\n# submission = df_1[['ImageID', 'Label']]\n\n# # Output to csv\n# submission.to_csv('Boetticher_RF__PCA_predictions_2.csv',header=True, index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Processing time for entire study"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Final process_time addition for entire program\n# Find sum of elements in process_time list \ntotal = 0\n\n# Iterate over each element in process_time list and add them in variable total \nfor ele in range(0, len(process_time)): \n    total = total + process_time[ele] \n    \n# Conversion\ndef convert(seconds): \n    min, sec = divmod(seconds, 60) \n    hour, min = divmod(min, 60) \n    return \"%d:%02d:%02d\" % (hour, min, sec) \n      \n# Time output\nn = total\nconvert(n)\n\n# Results\nprint(\"Total time of model fitting and evalation for study (in seconds): \", total) \nprint('Formatted time: ', time.strftime(\"%H:%M:%S\", time.gmtime(total)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.10"}},"nbformat":4,"nbformat_minor":1}