{"cells":[{"metadata":{"_uuid":"91205fc1b01313867750410210588a2763508d53"},"cell_type":"markdown","source":"MNIST Hand written digits:\n* Playing with some dimensionality reduction techniques.\n* Visualization of reduced datasets. \n* Testing performance of some classifiers on reduced datasets.\n* Scoring 0.982 with SVC polynomial kernel on PCA-redused dataset with 200ish components. \n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import accuracy_score\n\nfrom sklearn.decomposition import PCA\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\nfrom sklearn.manifold import TSNE #Only visualization, sklearn does noe support .transform - method\n\nfrom sklearn import svm\nfrom sklearn.linear_model import SGDClassifier as SGD\nfrom sklearn.neighbors import KNeighborsClassifier as KNN\n\nimport time\nimport matplotlib.pyplot as plt\nimport matplotlib\n%matplotlib inline\n\nimport seaborn as sns\nimport warnings\nwarnings.filterwarnings('ignore')\n\ndata = pd.read_csv(\"../input/train.csv\")\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"635de3896db1ff559bd3609235f26c5518818cf7"},"cell_type":"markdown","source":"Only using 2000 datapoits when visualizing and testing multiple classifiers. "},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"y_train = data[:2000]['label'].values\nx_train_raw = data[:2000].drop(labels = 'label',axis = 1).values\n\ny_test = data[40000:42000]['label'].values\nx_test_raw = data[40000:42000].drop(labels = 'label',axis = 1).values\n\nscaler = StandardScaler()\n\ntrain_input = scaler.fit_transform(x_train_raw)\ntest_input = scaler.transform(x_test_raw)\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"95778e41099664c84f41e53c7a82c49fea676d67"},"cell_type":"markdown","source":"The MNIST dataset is arrays of 784 pixel values, representing 28x28 pixels hand written digits. Let's look at some:"},{"metadata":{"trusted":true,"_uuid":"b15ed823fec4bf8cd1a3b4aa3d4781d0fd7a244c","scrolled":true},"cell_type":"code","source":"\nfor i in range(25):\n    plt.subplot(5,5,i+1)\n    plt.imshow((x_train_raw[i*42]).reshape(28,28)) \n    plt.xticks([])\n    plt.yticks([])  \n    \n\n\n    ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a908d7261835c98d9074fb01613f8d73a313a05d"},"cell_type":"markdown","source":"Calling pca without specifying number of components gives Principal Comonent Analysis on all components. This allows us to see the variance of the different components to select a reasonable amount of components for representing the dataset. "},{"metadata":{"_uuid":"18c616ab95bf794a6c75da62c32d95dbc1b4e0d5","trusted":true},"cell_type":"code","source":"    \npca = PCA().fit(x_train_raw)\nplt.plot(pca.explained_variance_)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"66fb78ad12bb4909d2c580d47b88c5a2e5eb40f2"},"cell_type":"markdown","source":"Plotting some data from a reduced dataset. Last plot is the eigenvalues seen as images. "},{"metadata":{"trusted":true,"_uuid":"f929d6e7ca781d3576a9084e02b0decdcf82cb3c"},"cell_type":"code","source":"\n\npca = PCA(n_components = 100).fit(train_input)\neigenvalues = pca.components_[:5]\nfig1, ax1 = plt.subplots() \nax1.plot(eigenvalues[0])\nax1.plot(eigenvalues[1])\nax1.plot(eigenvalues[2])\nax1.set_title('Three first eigenvalues')\nplt.show()\n\nfig2,ax2 = plt.subplots()\nax2.plot(pca.explained_variance_)\nax2.set_title('Variance in components')\nplt.show()\n\nplt.figure(figsize = (10,5))\nfor i in range(len(eigenvalues)):\n    plt.subplot(1,5,i+1)\n    plt.imshow((eigenvalues[i]).reshape(28,28)) \n    plt.xticks([])\n    plt.yticks([])  \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7ba3ea9493172fe689dc18c43a42f06e4f548176"},"cell_type":"markdown","source":"Lets use some different dimensionality reduction methods to visualize data. \n\n**PCA:** Principal component analysis: Linear method that orthogonally projects data into lower dimensionality. Free to choose how many components to keep by calling n_ components. \n\n**LDA:** Linear discriminant analysis: Seeks to find a set of dimensions to maximize the variance along the axis of the new dimensions. n_components is maximum the number of classes - 1 in the dataset. So with 10 different digits, it gives maximum 9 dimensions. \n\n**TSNE:** t-distributed stochastic neighbour embedding. Non-linear visualization on high dimensional data. There is no transform-method for new data in sklearn, so it's not easily applicable for classification on unseen data.\n\n"},{"metadata":{"_uuid":"fe6708a5aee098bd8150a86e2ba123ad518fdf2c"},"cell_type":"markdown","source":"Lets see how they do with only two components each. Using unscaled training data as TSNE performed better without scaling."},{"metadata":{"trusted":true,"_uuid":"9af5c929e428f3aaf8a2fc1fcc10f1555b2bd8a1","scrolled":true},"cell_type":"code","source":"\npca = PCA(n_components = 2).fit(x_train_raw)\ntrain_pca = pca.transform(x_train_raw)\n\nlda = LDA(n_components =2).fit(x_train_raw,y_train)\ntrain_lda = lda.transform(x_train_raw)\n\ntrain_tsne = TSNE(n_components = 2).fit_transform(x_train_raw)\n\n\n\n ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"7a59d4fd3ec001055bec617a5ff6757fb425c50c"},"cell_type":"code","source":"ax = sns.scatterplot(\n                train_pca[:,0],train_pca[:,1],\n                hue = y_train, \n                palette = sns.color_palette( \"hls\",n_colors = 10)\n                )    \nax.set_title('PCA 2-comp')\nplt.show()\n\nax = sns.scatterplot(\n                train_lda[:,0],train_lda[:,1],\n                hue = y_train, \n                palette = sns.color_palette( \"hls\",n_colors = 10)\n                )    \nax.set_title('LDA 2-comp')\nplt.show()\n\nax = sns.scatterplot(\n                train_tsne[:,0],train_tsne[:,1],\n                hue = y_train, \n                palette = sns.color_palette( \"hls\",n_colors = 10)\n                )    \nax.set_title('TSNE')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"055f9c6df221d65787142187c758a6d223afca08"},"cell_type":"markdown","source":"Now, lets see how several components in PCA/LCA relates. "},{"metadata":{"trusted":true,"_uuid":"2861198a4bce79fcb974dcd0bff6bd67ed53c4eb"},"cell_type":"code","source":"pca = PCA(n_components = 9).fit(x_train_raw)\ntrain_pca = pca.transform(x_train_raw)\n\nlda = LDA(n_components =9).fit(x_train_raw,y_train)\ntrain_lda = lda.transform(x_train_raw)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"950745bc931ae1c23ee9523d336296f655444595","_kg_hide-input":false},"cell_type":"code","source":"fig, axs = plt.subplots(2,2, figsize =(12,12))\naxs[0,0].scatter(\n    train_pca[:,0],train_pca[:,1],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[0,0].set_title('PCA 1st and 2nd component')\n\naxs[0,1].scatter(\n    train_pca[:,1],train_pca[:,2],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[0,1].set_title('PCA 2nd and 3rd component')\n\naxs[1,0].scatter(\n    train_pca[:,2],train_pca[:,3],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[1,0].set_title('PCA 3rd and 4th component')\n\naxs[1,1].scatter(\n    train_pca[:,3],train_pca[:,4],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[1,1].set_title('PCA 4th and 5th component')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"90e99aa70a65cdbb40a8e4af917bae9c45573bfe"},"cell_type":"code","source":"fig, axs = plt.subplots(2,2, figsize =(12,12))\naxs[0,0].scatter(\n    train_lda[:,0],train_lda[:,1],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[0,0].set_title('LDA 1st and 2nd component')\n\naxs[0,1].scatter(\n    train_lda[:,1],train_lda[:,2],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[0,1].set_title('LDA 2nd and 3rd component')\n\naxs[1,0].scatter(\n    train_lda[:,2],train_lda[:,3],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[1,0].set_title('LDA 3rd and 4th component')\n\naxs[1,1].scatter(\n    train_lda[:,3],train_lda[:,4],\n    c = y_train,\n    cmap = plt.get_cmap('gist_rainbow' ),\n    alpha = 0.5) \naxs[1,1].set_title('LDA 4th and 5th component')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"918f76162372ccd5796cae8f11244586af31bede"},"cell_type":"markdown","source":"We see that LDA is more able to map the variance down to few dimensions. Now: lets see how some classifiers perform on the redused datasets. Using SVM ploynomial and rbf-kernels, k-NN and logistic regression.  Parameters for the rbf-kernel is stolen from Saba Siddiqi's kernel: https://www.kaggle.com/sabasiddiqi/svm-classification-parameter-selection-0-968"},{"metadata":{"trusted":true,"_uuid":"a9ad6685b59fba5c6a71699afe77f72b3241b6ce"},"cell_type":"code","source":"def testClassifiers(x_train_scaled, y_train, x_test_scaled, y_test,n_components,use_lda = False ):\n    start = time.time()\n    pca = PCA(n_components = n_components)      \n    models = (\n        svm.SVC(kernel = 'rbf', C = 10, gamma = 0.001),\n        svm.SVC(kernel = 'poly' , degree = 2),\n        KNN(n_neighbors = 5),\n        #SGD(loss = 'log') #logistic regression    \n       )    \n    train_pca = pca.fit_transform(x_train_scaled)\n    test_pca = pca.transform(x_test_scaled)    \n    predictions = []\n    \n    for model in models:\n        model.fit(train_pca,y_train)\n        predictions.append(model.predict(test_pca))\n    \n    if use_lda == True:\n        lda = LDA(n_components = n_components)\n        train_lda = lda.fit_transform(x_train_scaled,y_train)\n        test_lda = lda.transform(x_test_scaled)\n        for model in models:\n            model.fit(train_lda,y_train)\n            predictions.append(model.predict(test_lda))\n        labels = ['PCA_SVC_rbf','PCA_SVC_poly','PCA_KNN','LDA_SVC_rbf','LDA_SVC_poly','LDA_KNN']\n    else:\n        \n        labels = ['SVC_rbf','SVC_poly','KNN']\n\n    scores = pd.DataFrame()\n    \n    for i,label in enumerate(labels):\n        scores[label] = [accuracy_score(y_test,predictions[i])]\n    scores['time'] = [time.time()-start]\n    scores['components'] = n_components\n    return scores\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"932cf9b6e627dfc7427f81e4f1134ea30ae5df58"},"cell_type":"markdown","source":"Calculating scores of the \"first\" few components where LDA can be compared to PCA:"},{"metadata":{"trusted":true,"_uuid":"052c53207c01eac2be7924c54c1ffb21ebe850bd","scrolled":true},"cell_type":"code","source":"results = pd.DataFrame()\n\nfor components in range(2,18,2):\n    df = testClassifiers(train_input, y_train, test_input, y_test,components,use_lda=True)\n    results = pd.concat([results,df])\n     \nprint(results.head())\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5e6d46d851a9e7fdc913d66eb11a53de3e25aef3"},"cell_type":"markdown","source":"Plotting accuracy of different classifiers as a function of redused components.  LDA beats PCA for n_components < 5.  LDA has maximum 9 components."},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"783273f6ce006727bcfa8cb9c75e9b77704b681d"},"cell_type":"code","source":"labels = ['PCA_SVC_rbf','PCA_SVC_poly','PCA_KNN','LDA_SVC_rbf','LDA_SVC_poly','LDA_KNN']\nx = results['components'].values\n\nf,axs = plt.subplots(1,2, figsize = (10,5))\n\nfor label in labels:\n    axs[0].scatter(x , results[label].values)\naxs[0].legend(labels)\n\nfor label in labels:\n    axs[1].plot(x , results[label].values)\naxs[1].legend(labels)\n\n    \n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"46e7fb41e4ab55edc9bafb44d46c4b143521d7ce"},"cell_type":"markdown","source":""},{"metadata":{"_uuid":"ac5fac43a065569e014016a7fdece47de4a9a4e7"},"cell_type":"markdown","source":"So, lets plot the classifiers for all components of PCA.  "},{"metadata":{"trusted":true,"_uuid":"49c715d3ec69c574f9a35140b66f117c50dbf87f"},"cell_type":"code","source":"results = pd.DataFrame()\n\nfor components in range(2,100,4):\n    df = testClassifiers(train_input, y_train, test_input, y_test,components)\n    results = pd.concat([results,df])\n\nfor components in range(100,200,10):\n    df = testClassifiers(train_input, y_train, test_input, y_test,components)\n    results = pd.concat([results,df])\nfor components in range(200,784,40):\n    df = testClassifiers(train_input, y_train, test_input, y_test,components)\n    results = pd.concat([results,df])\n     \n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2d3442073fa7121764a3b57ff4d0e8ef0d000aed","scrolled":true},"cell_type":"code","source":"labels = ['SVC_rbf','SVC_poly','KNN']\n\nx = results['components'][5:].values\n\nf,axs = plt.subplots(1,2, figsize = (10,5))\n\nfor label in labels:\n    axs[0].scatter(x , results[label][5:].values)\naxs[0].legend(labels)\n\nfor label in labels:\n    axs[1].plot(x , results[label][5:].values)\naxs[1].legend(labels)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"fc64168a6ac9f60296e8bc321989ba5f6b97a4cf"},"cell_type":"markdown","source":"Most interesting is that the SVM's have a \"peak performance\" at different amount of components. Lets see how they can perform with more training data:\n(still keeping the last 2000 as a test set)"},{"metadata":{"trusted":true,"_uuid":"011ac7efa2a38b4a195098c8c58f3b71c425820a"},"cell_type":"code","source":"y_train = data[:40000]['label'].values\nx_train_raw = data[:40000].drop(labels = 'label',axis = 1).values\n\ny_test = data[40000:42000]['label'].values\nx_test_raw = data[40000:42000].drop(labels = 'label',axis = 1).values\n\nscaler = StandardScaler()\n\ntrain_input = scaler.fit_transform(x_train_raw)\ntest_input = scaler.transform(x_test_raw)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"bc7f6f58a2bd45d94fbb11483d863f7e76fe865a"},"cell_type":"markdown","source":"The Support Vector Machines scores both scores around 0.97 on the remaining 2000 elements of training data when using rbf-kernel on 90 components and poly-kernel on 200 components.. On the 21000 images in test.csv, it performs slightly better, 0.982. "},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"69a1c367166699d84c89f99478d199ce02bf0075"},"cell_type":"code","source":"models = (\n    svm.SVC(kernel = 'rbf', C = 10, gamma = 0.001),\n    svm.SVC(kernel = 'poly' , degree = 3) \n    )\nt_rbf = time.time()\npca = PCA(n_components = 90).fit(train_input)\ntrain_pca = pca.transform(train_input)\ntest_pca = pca.transform(test_input)\n    \nmodels[0].fit(train_pca,y_train)\npredictions = models[0].predict(test_pca)\nrbf_score =  accuracy_score(y_test, predictions)\nt_rbf = time.time()-t_rbf\n\nt_poly = time.time()\npca = PCA(n_components = 200).fit(train_input)\ntrain_pca = pca.transform(train_input)\ntest_pca = pca.transform(test_input)\n    \nmodels[1].fit(train_pca,y_train)\npredictions = models[1].predict(test_pca)\npoly_score =  accuracy_score(y_test, predictions)  \nt_poly = time.time()-t_poly\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"c1bacc4d16729344a3f442cc8b8f93ec719b3c41"},"cell_type":"code","source":"scores_final = pd.DataFrame([[rbf_score,t_rbf],[poly_score,t_poly]],columns = ['score','time'],index = ['rbf','poly'])\nprint(scores_final)\n","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}