{"cells":[{"metadata":{"_cell_guid":"d7a71f50-2d95-5715-7c97-8e4baf788bfb","_uuid":"7b26a4c84039a38b0305918838b6d3bdf2c9da6c"},"cell_type":"markdown","source":"# Analysis of the MNIST dataset\n\n\n\nIn this analysis we will apply the following methods to the MNIST classification problem:\n\n1) Random forest classification\n\n2) Principal component analysis (PCA) + k-nearest neighbours (kNN)"},{"metadata":{"_cell_guid":"b79fc28f-94bc-43f8-8aa6-71315abc8adf","_uuid":"678f9e68029199b7271475dc6db2467611fac40f","trusted":false,"collapsed":true},"cell_type":"code","source":"# load the modules\n\nimport numpy as np\nimport pandas as pd\nimport seaborn as sb\nsb.set_style(\"dark\")\nimport matplotlib.pyplot as plt\n\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.neighbors import KNeighborsClassifier\nfrom sklearn.decomposition import PCA\n\n%pylab inline","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"8bcd0340-39fa-6b3c-9b8d-3e9a18f079d1","_uuid":"085f83f624e137c319b56b4a65eb065d4a581da1","trusted":false,"collapsed":true},"cell_type":"code","source":"# We use this function in order to evaulate a classifier. It trains on a fraction of the data corresponding to \n# aplit_ratio, and evaulates on the rest of the data\n\ndef evaluate_classifier(clf, data, target, split_ratio):\n    trainX, testX, trainY, testY = train_test_split(data, target, train_size=split_ratio, random_state=0)\n    clf.fit(trainX, trainY)\n    return clf.score(testX,testY)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a747b6e0-cb46-7b34-8f29-092ef30733cb","_uuid":"fdb0058ff125e51f0419132efeabb063da5cc919","trusted":false,"collapsed":true},"cell_type":"code","source":"# read in the data\n\ntrain = pd.read_csv('../input/train.csv')\ntest  = pd.read_csv('../input/test.csv')\ntarget = train[\"label\"]\ntrain = train.drop(\"label\",1)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f121719e-4dfc-1dec-f071-639c7a22a944","_uuid":"b1cff8bdfaf7d8a0df2b1037f284692aab5db7cc","trusted":false,"collapsed":true},"cell_type":"code","source":"# plot some of the numbers\n\nfigure(figsize(5,5))\nfor digit_num in range(0,64):\n    subplot(8,8,digit_num+1)\n    grid_data = train.iloc[digit_num].as_matrix().reshape(28,28)  # reshape from 1d to 2d pixel array\n    plt.imshow(grid_data, interpolation = \"none\", cmap = \"bone_r\")\n    xticks([])\n    yticks([])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"bdc541b8-abb8-3507-483c-ae9940dfb9b2","_uuid":"b5855b3618540b8e2009c78a10ce2626d590d9cf","trusted":false,"collapsed":true},"cell_type":"code","source":"# check performance of random forest classifier, as function of number of estimators \n# here we only take 1000 data points to train\n\nn_estimators_array = np.array([1,5,10,50,100,200,500])\nn_samples = 10\nn_grid = len(n_estimators_array)\nscore_array_mu =np.zeros(n_grid)\nscore_array_sigma = np.zeros(n_grid)\nj=0\nfor n_estimators in n_estimators_array:\n    score_array=np.zeros(n_samples)\n    for i in range(0,n_samples):\n        clf = RandomForestClassifier(n_estimators = n_estimators, n_jobs=1, criterion=\"gini\")\n        score_array[i] = evaluate_classifier(clf, train.iloc[0:1000], target.iloc[0:1000], 0.8)\n    score_array_mu[j], score_array_sigma[j] = mean(score_array), std(score_array)\n    j=j+1","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"26b82ddd-10cc-47d5-0c10-bfef4b60b9f1","_uuid":"0c4de549c3ab140ce81ea4baf9f2d0d06993552d","trusted":false,"collapsed":true},"cell_type":"code","source":"# it looks like the performace saturates around 50-100 estimators\n\nfigure(figsize(7,3))\nerrorbar(n_estimators_array, score_array_mu, yerr=score_array_sigma, fmt='k.-')\nxscale(\"log\")\nxlabel(\"number of estimators\",size = 20)\nylabel(\"accuracy\",size = 20)\nxlim(0.9,600)\ngrid(which=\"both\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"d22b7510-3c48-5b5e-90bc-6e5d8294fb2d","_uuid":"88c17b2d8fda74fe180198520eba64989c6c109f"},"cell_type":"markdown","source":"Are there any feature that are particularly important? We can check this using clf.feature_importances:"},{"metadata":{"_cell_guid":"d4801307-c5fe-f96f-66a3-fc2db8ac7e1e","_uuid":"2d51ce9121238a23b7a5fb0fed6a03f6c39b49a5","trusted":false,"collapsed":true},"cell_type":"code","source":"importances = clf.feature_importances_\nindices = np.argsort(importances)[::-1]\n\nprint(\"Feature ranking:\")\nfor f in range(0,10):\n    print(\"%d. feature %d (%f)\" % (f + 1, indices[f], importances[indices[f]]))\n\n# Plot the feature importances\n\nfigure(figsize(7,3))\nplot(indices[:],importances[indices[:]],'k.')\nyscale(\"log\")\nxlabel(\"feature\",size=20)\nylabel(\"importance\",size=20)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3427c56a-7620-22dc-909a-2cf2fa821ff5","_uuid":"e98846e7b5eb6d5dc492654ee870252cf2fa7837"},"cell_type":"markdown","source":"It looks like there are no significantly important features (i.e., pixels) in the original data. Next, let us try to decompose the data using a principal component analysis (PCA):"},{"metadata":{"_cell_guid":"ee4f8e1c-ccc3-193c-1763-e1fd1d47852b","_uuid":"f845319da999fe0b2f0491b920a16df09b50cb77","trusted":false,"collapsed":true},"cell_type":"code","source":"pca = PCA(n_components=2)\npca.fit(train)\ntransform = pca.transform(train)\n\nfigure(figsize(6,5))\nplt.scatter(transform[:,0],transform[:,1], s=20, c = target, cmap = \"nipy_spectral\", edgecolor = \"None\")\nplt.colorbar()\nclim(0,9)\n\nxlabel(\"PC1\")\nylabel(\"PC2\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"6c302666-9ba8-18d9-1129-cc45ec68f715","_uuid":"c862811f20c445bcd5cad5e03dd3c425d3517dd7"},"cell_type":"markdown","source":"It is interesting to see how well PCA separates the feature space into visible clusters already for 2 components. Next, let's look at what happens if we increase the number of components in PCA. In particular, we would like to know how many components are needed to capture most of the variance in the data. For this we will use the pca.explained_variance_ratio function."},{"metadata":{"_cell_guid":"627aeeab-eddd-1aba-d19a-537e74c453c6","_uuid":"806dfab39154f0b211b79dd70fbed70c4a81fe89","trusted":false,"collapsed":true},"cell_type":"code","source":"n_components_array=([1,2,3,4,5,10,20,50,100,200,500])\nvr = np.zeros(len(n_components_array))\ni=0;\nfor n_components in n_components_array:\n    pca = PCA(n_components=n_components)\n    pca.fit(train)\n    vr[i] = sum(pca.explained_variance_ratio_)\n    i=i+1    ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b737cb19-bcaa-cd24-37b8-6ddb3e3ee684","_uuid":"a759fd31d1ed0d752c38ec63c28e3444b49e71cd","trusted":false,"collapsed":true},"cell_type":"code","source":"figure(figsize(8,4))\nplot(n_components_array,vr,'k.-')\nxscale(\"log\")\nylim(9e-2,1.1)\nyticks(linspace(0.2,1.0,9))\nxlim(0.9)\ngrid(which=\"both\")\nxlabel(\"number of PCA components\",size=20)\nylabel(\"variance ratio\",size=20)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"695af730-5d79-7c8c-97cb-a1fb2f8068bc","_uuid":"e526129c9198c4672db217a4eab89043ee7b5aed"},"cell_type":"markdown","source":"We see that ~100 PCA components are needed to capture ~90% of the variance in the data. This seems a lot of components. Maybe the more important question is: How good is our prediction as a function of number of components? Let's look at this next. We will train a kNN classifier on the PCA output."},{"metadata":{"_cell_guid":"cc65a69f-7073-b937-d975-88c98b5c66d2","_uuid":"d674314493d80e976145cb1fb86450552160e411","trusted":false,"collapsed":true},"cell_type":"code","source":"clf = KNeighborsClassifier()\nn_components_array=([1,2,3,4,5,10,20,50,100,200,500])\nscore_array = np.zeros(len(n_components_array))\ni=0\n\nfor n_components in n_components_array:\n    pca = PCA(n_components=n_components)\n    pca.fit(train)\n    transform = pca.transform(train.iloc[0:1000])\n    score_array[i] = evaluate_classifier(clf, transform, target.iloc[0:1000], 0.8)\n    i=i+1","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"68e04bb4-c783-4cba-cd6e-b3db3a5a279a","_uuid":"560b7d3e37460e99a5993de671d950be04385920","trusted":false,"collapsed":true},"cell_type":"code","source":"figure(figsize(8,4))\nplot(n_components_array,score_array,'k.-')\nxscale('log')\nxlabel(\"number of PCA components\", size=20)\nylabel(\"accuracy\", size=20)\ngrid(which=\"both\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"73193bdb-b039-e0d6-6db4-689799382f53","_uuid":"7ce87e2f407861c10b721561f5717cfe255d026a"},"cell_type":"markdown","source":"The accuracy seems to saturate at ~90% (roughly matching the performance of the random forest classifier) for >~20 PCA components. In fact, the accuracy even seems to drop for much larger numbers, even though a larger number of PCA components captures more of the variance in the data, as seen in the plot above. The drop in accuracy is probably due to overfitting.\n\nFinally, we will train on the whole training set and prepare a submit file for the Kaggle competition."},{"metadata":{"_cell_guid":"ed495d33-8254-4058-ef27-dae03ccf4d30","_uuid":"928172a69df2aa446c3537b39d3e3b6288ba306a","trusted":false,"collapsed":true},"cell_type":"code","source":"# PCA + kNN\n    \npca = PCA(n_components=50)\npca.fit(train)\ntransform_train = pca.transform(train)\ntransform_test = pca.transform(test)\n\nclf = KNeighborsClassifier()\nclf.fit(transform_train, target)\nresults=clf.predict(transform_test)\n\n# prepare submit file\n\nnp.savetxt('results.csv', \n           np.c_[range(1,len(test)+1),results], \n           delimiter=',', \n           header = 'ImageId,Label', \n           comments = '', \n           fmt='%d')\n\n# Kaggle score 0.97343","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b2aea948-1727-e33c-040a-57e1e749ccc8","_uuid":"e286af6126f2885ff1637558d2e115b4d913be0e"},"cell_type":"markdown","source":"This gives a score on Kaggle of 0.97343 - not too bad! Using the random forest classifier:"},{"metadata":{"_cell_guid":"d441f1ae-7ba7-9001-23b4-c72533568aec","_uuid":"75a45119bc956e0953f7daf03e341694feb38d60","trusted":false,"collapsed":true},"cell_type":"code","source":"# random forest classification\n\nclf = RandomForestClassifier(n_estimators = 100, n_jobs=1, criterion=\"gini\")\nclf.fit(train, target)\nresults=clf.predict(test)\n\n# prepare submit file\n\nnp.savetxt('results.csv', \n           np.c_[range(1,len(test)+1),results], \n           delimiter=',', \n           header = 'ImageId,Label', \n           comments = '', \n           fmt='%d')\n\n# Kaggle score ~0.96","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"6d476c2d-e855-0ca5-c444-649ba49af635","_uuid":"419118fe2567b6c2e9a02476ebc37eeba40c40e7"},"cell_type":"markdown","source":"This gives a slightly worse score (0.96).\n\nAny feedback on my analysis is more than welcome!"}],"metadata":{"_change_revision":0,"_is_fork":false,"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.5.2"}},"nbformat":4,"nbformat_minor":1}