{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Tabular Playground Series - July 2022\n\nBy Michael Mortenson\n\nThe dataset for this challenge is simulated manufacturing control data. The goal is to use unsupervised (clustering) to identify different control states. We are not told the number of control states, the units, time dependencies, or any other information about the data.\n\n## Domain Insights\n\nManufacturing products requires raw materials to be processed through a series of steps (and often many machines). Each step will not be perfectly exact, so the engineering design will have tolerances for the manufacture of each piece of a product. The job of manufacturing control to ensure that production stays efficient by detecting problems as machine pieces begin to wear, devices lose their calibration, or other problems arise.\n\nSince we are not told what each column in the data represents, we will need to analyze to decide how we should treat it. ","metadata":{}},{"cell_type":"code","source":"# Useful Packages\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom matplotlib import pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:23.965310Z","iopub.execute_input":"2022-07-28T06:18:23.965766Z","iopub.status.idle":"2022-07-28T06:18:23.972991Z","shell.execute_reply.started":"2022-07-28T06:18:23.965733Z","shell.execute_reply":"2022-07-28T06:18:23.971683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in data\nX = pd.read_csv('/kaggle/input/tabular-playground-series-jul-2022/data.csv')","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:23.990612Z","iopub.execute_input":"2022-07-28T06:18:23.991082Z","iopub.status.idle":"2022-07-28T06:18:24.738445Z","shell.execute_reply.started":"2022-07-28T06:18:23.991048Z","shell.execute_reply":"2022-07-28T06:18:24.737552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## General Visualization\n\nLet's take a look at the general spread of the data in each column of our dataset.","metadata":{}},{"cell_type":"code","source":"X1 = X.drop(columns='id')\nplt.figure(dpi=100, figsize=(15, 5))\nsns.boxplot(data=X1)\nplt.grid(True)\nplt.title(\"Distribution of each column\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:24.740049Z","iopub.execute_input":"2022-07-28T06:18:24.740647Z","iopub.status.idle":"2022-07-28T06:18:26.054338Z","shell.execute_reply.started":"2022-07-28T06:18:24.740612Z","shell.execute_reply":"2022-07-28T06:18:26.053366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I visually see that there are 4 main groupings here:\n* 00-06, float data, normally distributed. Probably dimension tolerance data.\n* 07-14. integer data, categorical counts, skewed (Poisson?). Probably counting some kind of feature.\n* 15-21, float data, normally distributed. Probably more dimension tolerance data. Very similar to 00-06.\n* 22-28, float data, normally distributed. Mean drifts and spread is wider than the other float data categories.","metadata":{}},{"cell_type":"markdown","source":"### PCA - Principle Component Analysis\n\nFrom what I can tell, PCA requires that you center your data and normalize it so that the variance is consistent across features. The Yeo-Johnson tranformation (default for sklearn's PowerTransformer) improves the normality of data, so I will go with that, but I am open to other ideas.","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import PowerTransformer\npt = PowerTransformer(method='yeo-johnson')\nX_tran = pt.fit_transform(X1)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:26.055473Z","iopub.execute_input":"2022-07-28T06:18:26.055762Z","iopub.status.idle":"2022-07-28T06:18:29.711741Z","shell.execute_reply.started":"2022-07-28T06:18:26.055734Z","shell.execute_reply":"2022-07-28T06:18:29.710676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualize the tranformed data\nplt.figure(dpi=100, figsize=(15, 5))\nsns.boxplot(data=X_tran)\nplt.grid(True)\nplt.title(\"Distribution of each column, Transformed\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:29.715131Z","iopub.execute_input":"2022-07-28T06:18:29.715561Z","iopub.status.idle":"2022-07-28T06:18:30.476609Z","shell.execute_reply.started":"2022-07-28T06:18:29.715519Z","shell.execute_reply":"2022-07-28T06:18:30.475352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Verify that the standard deviation of the tranformed data is consistent.\nstd = np.std(X_tran, axis=0)\nfig, ax = plt.subplots(figsize=[15,5])\nax.plot(std)\nax.set_title(\"Standard Deviation\")\nax.set_xlabel(\"Feature Number\")\nax.grid(True)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:30.478089Z","iopub.execute_input":"2022-07-28T06:18:30.478476Z","iopub.status.idle":"2022-07-28T06:18:30.675547Z","shell.execute_reply.started":"2022-07-28T06:18:30.478434Z","shell.execute_reply":"2022-07-28T06:18:30.674598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.decomposition import PCA\npca = PCA()\npca.fit(X_tran)\nvar_ratio = pca.explained_variance_ratio_\nexp_var = pca.explained_variance_\nsing_vals = pca.singular_values_\ncomponents = pca.components_","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:30.677044Z","iopub.execute_input":"2022-07-28T06:18:30.677569Z","iopub.status.idle":"2022-07-28T06:18:30.793167Z","shell.execute_reply.started":"2022-07-28T06:18:30.677524Z","shell.execute_reply":"2022-07-28T06:18:30.791821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The singular values show the relative 'strength' of each component (see info on how the singular value decomposition relates to PCA).","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=[15,5])\nax.plot(sing_vals)\nax.set_title(\"Singular Values\")\nax.set_xlabel(\"Component Number\")\nax.grid(True)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:30.795512Z","iopub.execute_input":"2022-07-28T06:18:30.796478Z","iopub.status.idle":"2022-07-28T06:18:31.024299Z","shell.execute_reply.started":"2022-07-28T06:18:30.796414Z","shell.execute_reply":"2022-07-28T06:18:31.023073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The main \"elbow\" occurs from 5 or 6 on, which would indicate the first 6 or 7 components are the most informative.","metadata":{}},{"cell_type":"code","source":"# Cumulative proportion of variance (from PC1 to PC6)   \npercent_explained = np.cumsum(var_ratio)\n\nfig, ax = plt.subplots(figsize=[15,5])\nax.plot(percent_explained)\nax.set_title(\"Percent Explained\")\nax.set_xlabel(\"Component Number\")\nax.set_xlim([0,28])\nax.set_ylim([0,1])\nax.grid(True)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:31.025334Z","iopub.execute_input":"2022-07-28T06:18:31.025637Z","iopub.status.idle":"2022-07-28T06:18:31.224726Z","shell.execute_reply.started":"2022-07-28T06:18:31.025609Z","shell.execute_reply":"2022-07-28T06:18:31.223930Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I find it a bit odd that the components look to be accounting for about equal variance in the data. Would this suggest they are all important, or did I mess up in the computation?\n\n* Kaggle hero suggested that the equal variance may be due to the Yeo-Johnson transformation (which makes features have similar variance).\n\nWhat would you suggest for how to move forward? Are there other visualizations or analyses I could do? Is PCA good enough, or is there another (nonlinear?) transformation that would be better?\n\n* One Kaggle hero suggested that I try PCA on each other \"groups\" without any preprocessing, since they seem to have similar behaviors within each group. I will try this next.","metadata":{}},{"cell_type":"markdown","source":"## PCA on each group","metadata":{}},{"cell_type":"code","source":"# The group labels\ng1 = [f\"f_{x:02d}\" for x in range(7)]\ng2 = [f\"f_{x:02d}\" for x in range(7, 14)]\ng3 = [f\"f_{x:02d}\" for x in range(14, 22)]\ng4 = [f\"f_{x:02d}\" for x in range(22, 29)]\ng_label_list = [g1, g2, g3, g4]\n\n# The group data\nX_g1 = X1[g1].copy()\nX_g2 = X1[g2].copy()\nX_g3 = X[g3].copy()\nX_g4 = X1[g4].copy()\ng_data = [X_g1, X_g2, X_g3, X_g4]","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:31.226082Z","iopub.execute_input":"2022-07-28T06:18:31.226608Z","iopub.status.idle":"2022-07-28T06:18:31.247153Z","shell.execute_reply.started":"2022-07-28T06:18:31.226576Z","shell.execute_reply":"2022-07-28T06:18:31.245805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initialize Storage\nvar_ratios = []\nexp_vars = []\nsing_val_list = []\ncomp_list = []\nprojections = []\n\n# Loop through groups\nfor data in g_data:\n    pca = PCA()\n    projection = pca.fit_transform(data)\n    projections.append(projection)\n    var_ratios.append(pca.explained_variance_ratio_)\n    exp_vars.append(pca.explained_variance_)\n    sing_val_list.append(pca.singular_values_)\n    comp_list.append(pca.components_)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:31.251150Z","iopub.execute_input":"2022-07-28T06:18:31.251514Z","iopub.status.idle":"2022-07-28T06:18:31.367186Z","shell.execute_reply.started":"2022-07-28T06:18:31.251483Z","shell.execute_reply":"2022-07-28T06:18:31.365800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Singular Value Plots\nfig, ax = plt.subplots(nrows=len(g_label_list), ncols=2, figsize=[15,15])\ngroup_names = ['g1', 'g2', 'g3', 'g4']\n\nfor idx, label in enumerate(group_names):\n    ax[idx,0].plot(sing_val_list[idx])\n    ax[idx,0].set_title(f\"Singular Values {label}\")\n    ax[idx,0].grid(True)\n    \n    percent_explained = np.cumsum(var_ratios[idx])\n    ax[idx,1].plot(percent_explained)\n    ax[idx,1].set_title(f\"Percent Explained {label}\")\n    ax[idx,1].set_ylim([0,1])\n    ax[idx,1].grid(True)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:31.369138Z","iopub.execute_input":"2022-07-28T06:18:31.369907Z","iopub.status.idle":"2022-07-28T06:18:32.402068Z","shell.execute_reply.started":"2022-07-28T06:18:31.369858Z","shell.execute_reply":"2022-07-28T06:18:32.401179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The Gaussian distributed groups all have similar looking plots to the total PCA plots, but the skewed, integer count group seems to have some variation. \n\nI wonder what would happen if I did t-SNE or UMAP (see below) on the new PCA stuff.","metadata":{}},{"cell_type":"markdown","source":"### t-SNE : t-Distributed Stochastic Neighbor Embedding\n\nOne method for visualizing the data that was recommened to me following my original posting was t-SNE. \n\nIn doing some basic googling, it uses the KL-divergence to create a low-dim embedding of  high-dimensional data.\n\nThings to remember:\n* t-SNE is not convex (intialization-dependent results)\n* perplexity needs to be less than the number of points (easy here)\n* Cluster sizes in a t-SNE plot mean nothing\n* Distances btween clusters is also not informative\n* Lower perplexity captures the local structures better and higher perplexity captures the global structures better.\n","metadata":{}},{"cell_type":"code","source":"from sklearn.manifold import TSNE\nX2 = X1.sample(frac=.10, random_state=0)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:32.403327Z","iopub.execute_input":"2022-07-28T06:18:32.403913Z","iopub.status.idle":"2022-07-28T06:18:32.416314Z","shell.execute_reply.started":"2022-07-28T06:18:32.403868Z","shell.execute_reply":"2022-07-28T06:18:32.415226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"perplexities = [5.0, 10.0, 20.0, 30.0, 40.0, 50.0]\nembeddings = []\nfor perplexity in perplexities:\n    tsne = TSNE(n_components=2, learning_rate='auto', perplexity=perplexity, verbose=1)\n    X_embedded = tsne.fit_transform(X2)\n    embeddings.append(X_embedded)","metadata":{"execution":{"iopub.status.busy":"2022-07-28T06:18:32.417955Z","iopub.execute_input":"2022-07-28T06:18:32.418460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(ncols=1, nrows=len(perplexities), figsize=[15,45])\n\nfor idx, perplexity in enumerate(perplexities):\n    x = embeddings[idx][:,0]\n    y = embeddings[idx][:,1]\n    ax[idx].scatter(x,y)\n    ax[idx].set_ylabel(f\"Peplexity: {perplexity}\")\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So, t-sne doesn't look like it's giving great results. Even varying perplexity, the data is not really splitting into different clusters.","metadata":{}},{"cell_type":"markdown","source":"# t-SNE on each PCA group","metadata":{}},{"cell_type":"code","source":"perplexities = [5.0, 10.0, 20.0, 30.0, 40.0, 50.0]\n\nfig, ax = plt.subplots(ncols=len(projections), nrows=len(perplexities), figsize=[15,45])\n\nfor ridx, perplexity in enumerate(perplexities):\n    for cidx, projection in enumerate(projections):\n        \n        tsne = TSNE(n_components=2, learning_rate='auto', perplexity=perplexity, verbose=1)\n        projection_df = pd.DataFrame(projection)\n        X_samp = projection_df.sample(frac=.10, random_state=0)\n        X_embedded = tsne.fit_transform(X_samp)\n        \n        x = X_embedded[:,0]\n        y = X_embedded[:,1]\n        ax[ridx, cidx].scatter(x,y, alpha=0.1)\n        ax[ridx, cidx].set_ylabel(f\"Peplexity: {perplexity}\")\n        \n\nplt.tight_layout()\nplt.show()   ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Running t-SNE on the PCA of each of the four groups doesn't give us much to go on either. Still blobs.","metadata":{}},{"cell_type":"markdown","source":"### UMAP\n\nSome people have suggested using UMAP, so I will try it out and see how it goes.","metadata":{}},{"cell_type":"code","source":"import umap","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"neighbor_size = [5.0, 50.0, 100.0, 500.0]\numap_embeds = []\nfor n_neighbors in neighbor_size:\n    reducer = umap.UMAP(random_state=0, n_neighbors=n_neighbors, n_components=2, verbose=True)\n    \n    # Fit the umap enbedding with the PowerTransformer normalized data\n    X_embed = reducer.fit_transform(X_tran)\n    umap_embeds.append(X_embed)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(ncols=1, nrows=len(neighbor_size), figsize=[15,45])\n\nfor idx, n_neighbors in enumerate(neighbor_size):\n    x = umap_embeds[idx][:,0]\n    y = umap_embeds[idx][:,1]\n    ax[idx].scatter(x,y)\n    ax[idx].set_ylabel(f\"Num Neighbors: {n_neighbors}\")\n\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The UMAP is just as single-cluster-y as t-sne, so it looks like PCA gives us the most info (even though it's not much). \n\nIf you have tips for how to improve the umap or t-SNE visualization, let me know. :)","metadata":{}},{"cell_type":"markdown","source":"### Clustering - Gaussian Mixture Model\n\nnTry out clustering with the first 7 components of PCA. People have indicated that this is a good place to start, so let's start there. Also, I tried k-means and it wasn't great, so let's try the more flexible GMM.","metadata":{}},{"cell_type":"code","source":"# Set number of components\nnum_comps = 6\n\n# Grab components from PCA\npca = PCA(n_components=num_comps)\npca_array = pca.fit_transform(X_tran)\npca_df = pd.DataFrame(pca_array)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.mixture import GaussianMixture\n\n# Intialize the GMM model\ngmm = GaussianMixture(n_components = num_comps, random_state=0)\n\n# Fit the GMM model for the dataset\ngmm.fit(pca_df)\n\n# Assign a label to each sample\nlabels = gmm.predict(pca_df)\npca_df['Clusters']= labels","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualize the clusters and their distributions\ng = sns.PairGrid(pca_df, vars=list(range(num_comps)), hue=\"Clusters\", palette=\"tab10\")\ng.map_diag(sns.histplot)\ng.map_offdiag(sns.scatterplot)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Wow! Some of these clusters look very distinct, which is great. This is especially true of the first couple columns (which is consistent with the idea of principle components). Let's try submitting and see how it does. My previous best scroe was 0.24602.","metadata":{}},{"cell_type":"markdown","source":"### 4-part PCA Clustering\n\nI want to try doing clustering using the top 3 components of each of the 4 subgroups' individual PCAs.","metadata":{}},{"cell_type":"code","source":"# Combine the top j components of each group into one\nj = 3\n\ntop_g1 = projections[0][:,:j]\ntop_g2 = projections[1][:,:j]\ntop_g3 = projections[2][:,:j]\ntop_g4 = projections[3][:,:j]\n\nnew_data = np.concatenate((top_g1, top_g2, top_g3, top_g4), axis=1)\npca = PCA()\nnew_data = pca.fit_transform(new_data)\n\n# pt = PowerTransformer(method='yeo-johnson')\n# new_data = pt.fit_transform(new_data)\n# pca_df = pd.DataFrame(new_data)\n\n# Visualize the tranformed data\nplt.figure(dpi=100, figsize=(15, 5))\nsns.boxplot(data=pca_df)\nplt.grid(True)\nplt.title(\"Distribution of each column, Transformed\")\nplt.show()\n\n# Intialize the GMM model\ngmm = GaussianMixture(n_components = num_comps, random_state=0)\n\n# Fit the GMM model for the dataset\ngmm.fit(pca_df)\n\n# Assign a label to each sample\nlabels = gmm.predict(pca_df)\npca_df['Clusters']= labels","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualize the clusters and their distributions\ng = sns.PairGrid(pca_df, vars=list(range(num_comps)), hue=\"Clusters\", palette=\"tab10\")\ng.map_diag(sns.histplot)\ng.map_offdiag(sns.scatterplot)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prep for submission\nsubmission = pd.DataFrame()\nsubmission['id'] = pca_df.index\nsubmission['Predicted'] = pca_df.loc[:,['Clusters']]\nsubmission.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Save for submitting\nsubmission.to_csv('submission.csv',index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cool! With 7 GMM components, that gave me a score of 0.27272. An improvement of 0.03 over using my previous method (which was just straight k-means).\n\nI wonder how I would do with different numbers of components?\n\n* With 6 GMM components, I get a score of 0.28254, an improvment of another 0.1.\n* With 5 GMM components, I get a score of 0.22158, a backslide in performance. \n\nSplitting the data into the 4 groups I observed visually, and doing PCA in each one, I thought it might be interesting to use the top couple components with the GMM model to try and cluster the data.\n* For 5 GMM components with the top 3 of each group's PCA components I get a score of 0.25911\n* For 6 GMM components with the top 2 of each group's PCA components I get a score of 0.17945 (Eek!)\n* For 6 GMM components with the top 3 of each group's PCA components fed through a PCA, I get 0.27467 (Bettah!)\n* For 6 GMM components with the top 3 of each group's PCA components fed through a PCA and PowerTransformer, I get 0.17736 (not better.)\n\nTrying a slightly different tactic, I used the top 3 PCA components to create a dataset that I fed into PCA. I only trained the GMM on the top 3 overal components and got a score of 0.24391.\n","metadata":{}}]}