\n
\n\nQ-Q plots, aka **Quantile-Quantile** plots, are used to **visually compare** how similar two distributions are to each other. They consist of plotting the quantiles (i.e. regular intervals) of the **observed** distribution against the quantiles of the **theoretical** distribution. The closer the Q-Q plots are to forming a **straight line**, the more confident you can be that the observed and theoretical distributions are the **same**. \n\n**Normal Q-Q plots**","metadata":{}},{"cell_type":"code","source":"# Normal Q-Q plots\nfigure = plt.figure(figsize = (16,12))\nfor i in range(len(data.columns)):\n \n # Q-Q plot\n ax = plt.subplot(6,5, i+1)\n stats.probplot(data.iloc[:,i], dist='norm', plot=plt)\n \n # Aesthetics\n ax.get_lines()[0].set_markersize(6.0)\n ax.get_lines()[1].set_linewidth(3.0)\n ax.set_xticklabels([])\n ax.set_yticklabels([])\n plt.title(data.columns[i])\n \nfigure.tight_layout(h_pad=1.0, w_pad=0.5)\nplt.suptitle('Normal Q-Q Charts', y=1.02, fontsize=20)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-18T10:58:50.812582Z","iopub.execute_input":"2022-07-18T10:58:50.812992Z","iopub.status.idle":"2022-07-18T10:59:00.270311Z","shell.execute_reply.started":"2022-07-18T10:58:50.812962Z","shell.execute_reply":"2022-07-18T10:59:00.269058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Even though the features *f_22* to *f_28* failed the Shapiro-Wilk test, they still appear to be quite close to being normally distributed. This behaviour could be because these features are made up of a **mixture** of normal distributions. There is not an easy way to verify this however.\n\n3.4 Q-Q plots
\n\n
\n\n**Poisson Q-Q plots**","metadata":{}},{"cell_type":"code","source":"# Poisson Q-Q plots\nfigure = plt.figure(figsize = (16,5))\nfor i, col in enumerate(int_feats):\n \n # Q-Q plot\n ax = plt.subplot(2, 4, i+1)\n mu = data[col].mean()\n stats.probplot(data[col], dist='poisson', sparams=(mu,), plot=plt)\n \n # Aesthetics\n ax.get_lines()[0].set_markersize(6.0)\n ax.get_lines()[1].set_linewidth(3.0)\n ax.set_xticklabels([])\n ax.set_yticklabels([])\n plt.title(col)\n \nfigure.tight_layout(h_pad=1.0, w_pad=0.5)\nplt.suptitle('Poisson Q-Q Charts', y=1.02, fontsize=20)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-18T10:59:00.272197Z","iopub.execute_input":"2022-07-18T10:59:00.272910Z","iopub.status.idle":"2022-07-18T10:59:03.866156Z","shell.execute_reply.started":"2022-07-18T10:59:00.272862Z","shell.execute_reply":"2022-07-18T10:59:03.864867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see more clearly that these features are not distributed according to independent Poisson distributions. However, some of them are quite close. It could be also that these features are made up of a **mixture** of Poisson distributions. Unfortunately, there isn't an easy way to verify this.\n\n
\n
\n
\n\n* Features *f_00* to *f_06* and *f_14* to *f_21* are **independent** of all other features.\n* Discrete features (*f_07* to *f_13*) and features *f_22* to *f_28* are **weakly dependent** of each other.","metadata":{}},{"cell_type":"code","source":"# Heatmap of correlations\nplt.figure(figsize=(7,5))\nsns.heatmap(data.corr().abs(), cmap='Greens', vmin=0, vmax=1)\nplt.title('Absolute correlations')","metadata":{"execution":{"iopub.status.busy":"2022-07-18T10:59:03.868260Z","iopub.execute_input":"2022-07-18T10:59:03.868620Z","iopub.status.idle":"2022-07-18T10:59:04.535596Z","shell.execute_reply.started":"2022-07-18T10:59:03.868586Z","shell.execute_reply":"2022-07-18T10:59:04.534324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Elbow method\n\n3.5 Correlations
\n\n
\n\n\nThe **elbow method** is a practical way to determine the number of clusters in a dataset. It works by plotting the **inertia** (or sometimes distortion) against the **number of clusters**, where **inertia** is defined to be the sum of squared distances of samples to their closest cluster center, i.e. a measure of the models bias. The '**elbow**' (point of sudden flattening) of the curve is then chosen to be the optimal number of clusters in the dataset.\n\n4.1 How it works
\n![](\"https://www.oreilly.com/library/view/statistics-for-machine/9781788295758/assets/995b8b58-06f1-4884-a2a1-f3648428e947.png\")
\n\n
\n\nSee my [discussion post](https://www.kaggle.com/competitions/tabular-playground-series-jul-2022/discussion/335079) where I used 50 clusters and the whole dataset. To save time here, we will just use 30 clusters and 10% of the data.","metadata":{}},{"cell_type":"code","source":"%%time\n\ninertias = []\nfor k in range(1,30):\n km = KMeans(n_clusters=k)\n km.fit(data.iloc[:10000,:])\n inertias.append(km.inertia_)\n\n# Plot inertias\nplt.figure(figsize=(16,6))\nplt.plot(range(1,30), inertias, 'bx-')\nplt.xlabel('Number of clusters, k')\nplt.ylabel('Inertia')\nplt.title('Elbow method')\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-11T15:49:37.554060Z","iopub.execute_input":"2022-07-11T15:49:37.554446Z","iopub.status.idle":"2022-07-11T15:51:19.687720Z","shell.execute_reply.started":"2022-07-11T15:49:37.554413Z","shell.execute_reply":"2022-07-11T15:51:19.686236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is **hard** to tell exactly what the optimal value for the number of clusters should be since the curve is quite smooth. We will go with **k=7** for now, but it might be worth experimenting with different values of k as well. \n\n# 5. Competition metric\n\nIt is worth spending some time trying to understand the competition metric. This is called the **Adjusted Rand Index (ARI)**. But to do this, we first need to look at the **Rand Index (RI)**.\n\n4.2 Applying it
\n\n
\n\nThe Rand Index (named after **William Rand** from 1971) is a measure of **similarity** between the predicted clusters and the ground truth clusters. It looks at whether **pairs** of data points are in the same or different clusters. Let's work through an **example** to see how it works.\n\n$$\\Large RI = \\frac{a+b}{{n \\choose 2}}$$\n\n5.1 Rand Index
\n![](\"https://i.postimg.cc/3Ng0nRwZ/RI-n5.jpg\")
\n\n\nTo work out the Rand Index, we need to calculate **two quantities**:\n* $a$ = # pairs in the **same** cluster in the prediction and the **same** cluster in the ground truth.\n* $b$ = # pairs in **different** clusters in the prediction and **different** clusters in the ground truth.\n\nThis can be a **bit confusing** but for example, the points ${\\color{orange} \\alpha}, {\\color{orange} \\beta}$ are in the same cluster in the prediction (orange) and in the same cluster in the ground truth (orange), so the pair $\\{{\\color{orange} \\alpha}, {\\color{orange} \\beta}\\}$ counts towards $a$. On the other hand, the points ${\\color{orange} \\alpha}, {\\color{green} \\delta}$ are in different clusters in both the prediction and ground truth (orange, green) so the pair $\\{{\\color{orange} \\alpha}, {\\color{green} \\delta}\\}$ counts towards $b$. \n\n
\n\nIf we continue like this (**check this yourself**), you will find that the pairs $\\{{\\color{orange} \\alpha}, {\\color{orange} \\beta}\\}, \\{{\\color{green} \\delta}, {\\color{green} \\epsilon}\\}$ are in the **same** cluster for both prediction and ground truth so $a=2$ and the pairs $\\{{\\color{orange} \\alpha}, {\\color{green} \\delta}\\}, \\{{\\color{orange} \\alpha}, {\\color{green} \\epsilon}\\}, \\{{\\color{orange} \\beta}, {\\color{green} \\delta}\\}, \\{{\\color{orange} \\beta}, {\\color{green} \\epsilon}\\}, \\{{\\color{red} \\gamma}, {\\color{green} \\delta}\\}, \\{{\\color{red} \\gamma}, {\\color{green} \\epsilon}\\}$ are in **different** clusters for both prediction and ground truth so $b=6$.\n\nGreat, so putting the numbers in we find that $RI=\\frac{2+6}{10}=0.8$.\n\n
\n
\n\n* RI lies **between 0 and 1**. The closer to 1 the better.\n* If the prediction is **perfect**, i.e. equal to the ground truth, then **RI = 1**.\n* We say a pair of points is in '**agreement**' if they count towards a or b (above), and in '**disagreement**' otherwise. If we pick two points at **random**, RI gives the **probability** that this pair of points is in agreement. (i.e. the predicted 'state' of the pair is 'correct')\n* RI is equivalent to **accuracy** when viewed from a **binary classification** problem over the **pairs** of data points. In particular, each pair is either in agreement (1) or in disagreement (0), in which case $a=\\text{True Positives} \\, (TP)$ and $b=\\text{True Negatives} \\, (TN)$ so the Rand Index becomes:\n\n$$RI = \\frac{TP + TN}{TP + FP + FN + TN} = \\, \\text{accuracy of pairs}$$\n* The only main **downside** to RI is that the **expected value** of RI, $\\mathbb{E}(RI)$, isn't the same for different clustering problems. That means, some problems are **easier** to get a good RI score than others so we can't really **compare** RI between different problems. This is where the adjusted RI comes in. \n\n5.2 Properties of RI
\n\n
\n\nThe **Adjusted Rand Index** (ARI) (introduced by **Hubert** and **Arabie** in 1985) is a 'corrected-for-chance' version of the Rand Index. It substracts RI by the expected value of RI for the specific clusterting problem. It then scales this number so that it has a maximum value of 1. \n\n$$\\Large ARI = \\frac{RI - \\mathbb{E}(RI)}{max(RI)-\\mathbb{E}(RI)}$$\n\n**Properties of ARI:**\n* It has a **maximum value of 1** but **no minimum** value (it can be negative). \n* If the prediction is **perfect**, i.e. equal to the ground truth, then **ARI = 1**.\n* A score of **0**, means the prediction is as good as picking all the clusters at **random**. \n* It is **comparable** between different clustering problems as its expected value is **constant** (0).\n\n5.3 Adjusted Rand Index
\n\n\nWe know how to work out the RI and also that $max(RI)=1$, but working out the expectation $\\mathbb{E}(RI)$ is much **trickier**. Hubert derived the following (rather complicated) **formula** for the entire ARI. See the **appendix** if you are interested to see the derivation.\n\n$$\n\\large ARI = \\frac{ \\left. \\sum_{ij} \\binom{n_{ij}}{2} - \\left[\\sum_i \\binom{a_i}{2} \\sum_j \\binom{b_j}{2}\\right] \\right/ \\binom{n}{2} }{ \\left. \\frac{1}{2} \\left[\\sum_i \\binom{a_i}{2} + \\sum_j \\binom{b_j}{2}\\right] - \\left[\\sum_i \\binom{a_i}{2} \\sum_j \\binom{b_j}{2}\\right] \\right/ \\binom{n}{2} }\n$$\n\n**Note:** This is equivalent to the ARI formula above. \n\n
\n
\n\nFirst, we start by **formalising** the clustering problem. We denote our dataset with $n$ objects by $S = \\{o_1, o_2, \\ldots, o_n \\}$ (Each element is just a data point). Let's represent the ground truth and predicted clustering by two **partitions**: $X=\\{X_1, \\ldots, X_r\\}$ and $Y=\\{Y_1, \\ldots, Y_s\\}$, respectively.\n\nContinuing from our previous example, $X=\\{X_1, X_2\\}$ with $X_1=\\{{\\color{orange} \\alpha},{\\color{orange} \\beta},{\\color{orange} \\gamma}\\}$, $X_2=\\{{\\color{green} \\delta},{\\color{green} \\epsilon}\\}$ and $Y=\\{Y_1, Y_2, Y_3\\}$ with $Y_1=\\{{\\color{orange} \\alpha},{\\color{orange} \\beta}\\}$, $Y_2=\\{{\\color{red} \\gamma}\\}$, $Y_3=\\{{\\color{green} \\delta},{\\color{green} \\epsilon}\\}$.\n\n5.4 ARI example
\n\n\nThen we draw a **contingency** table. Each entry, $n_{i,j}$, denotes how many data points there are in common between the ground truth cluster $X_i$ and the predicted cluster $Y_j$. Mathematically, the formula is $n_{i,j} = |X_i \\cap Y_j |$, i.e. the **intersection**.\n\n$$ \n\\begin{array}{c|cccc|c}\n{{} \\atop X}\\!\\diagdown\\!^Y &\nY_1&\nY_2&\n\\cdots&\nY_s&\n\\text{sums}\n\\\\\n\\hline\nX_1&\nn_{11}&\nn_{12}&\n\\cdots&\nn_{1s}&\na_1\n\\\\\nX_2&\nn_{21}&\nn_{22}&\n\\cdots&\nn_{2s}&\na_2\n\\\\\n\\vdots&\n\\vdots&\n\\vdots&\n\\ddots&\n\\vdots&\n\\vdots\n\\\\\nX_r&\nn_{r1}&\nn_{r2}&\n\\cdots&\nn_{rs}&\na_r\n\\\\\n\\hline\n\\text{sums}&\nb_1&\nb_2&\n\\cdots&\nb_s&\n\\end{array}\n$$\n\nFor example, to work out $n_{11}$, we look at clusters $X_1=\\{{\\color{orange} \\alpha},{\\color{orange} \\beta},{\\color{orange} \\gamma}\\}$ and $Y_1=\\{{\\color{orange} \\alpha},{\\color{orange} \\beta}\\}$ and find the points which appear in both of them. In this case, there are 2: $X_1 \\cap Y_1 = \\{{\\color{orange} \\alpha},{\\color{orange} \\beta}\\}$ so $n_{11}=2$. If we continue like this (**check this yourself**) we get:\n\n$$ \n\\begin{array}{c|ccc|c}\n{{} \\atop X}\\!\\diagdown\\!^Y &\nY_1&\nY_2&\nY_3&\n\\text{sums}\n\\\\\n\\hline\nX_1&\n2&\n1&\n0&\n3\n\\\\\nX_2&\n0&\n0&\n2&\n2\n\\\\\n\\hline\n\\text{sums}&\n2&\n1&\n2&\n\\end{array}\n$$\n\n
\n\nWe are almost there. To work out ARI, we need to calculate these 3 quantities: $\\sum_{i,j} {n_{ij} \\choose 2}$, $\\sum_{i} {a_{i} \\choose 2}$, $\\sum_{j} {b_{j} \\choose 2}$.\n\nFirst recall the formula ${n \\choose 2} = \\frac{n(n-1)}{2}$, which we will be using a lot. E.g. ${5 \\choose 2} = 10$\n\n1. $\\sum_{i,j} {n_{ij} \\choose 2} = {2 \\choose 2} + {1 \\choose 2} + {0 \\choose 2} + {0 \\choose 2} + {0 \\choose 2} + {2 \\choose 2} = 1 + 0 + 0 + 0 + 0 + 1 = 2$.\n\n2. $\\sum_{i} {a_{i} \\choose 2} = {3 \\choose 2} + {2 \\choose 2} = 3 + 1 = 4$.\n\n3. $\\sum_{j} {b_{j} \\choose 2} = {2 \\choose 2} + {1 \\choose 2} + {2 \\choose 2} = 1 + 0 + 1 = 2$.\n\nSo if we plug everything into the formula, we get $ARI = \\frac{2 - (4 \\times 2) / 10}{(4 + 2)/2 - (4 \\times 2)/10} = 0.55$. \n\nThe ARI score (0.55) is quite a bit smaller than the RI score (0.80) we got earlier. This is somewhat expected though with such a small clustering problem; since the number of points $n$ is small, the problem is relatively easy so the ARI makes a **large adjustment**.\n\n# 6. Modelling\n\n
\n
\n\nIt is always important to scale the data for clustering problems so that it is easier to compare the distance between data points. \n\n6.1 Scaling
\n![](\"https://149695847.v2.pressablecdn.com/wp-content/uploads/2021/09/image-47.png\")
\n\n
\n\nk-Means is an **iterative** clustering algorithm that works as follows:\n1. Choose coordinates (e.g. randomly) for the locations of the k centroids.\n2. Group datapoints together by finding the nearast centroid. (There will always be k goups).\n3. Calculate the new centre of each centroid by taking the mean position of datapoints in each group.\n4. Iterative until the centroids stop moving by a significant amount.\n\n6.2 k-Means
\n![](\"https://upload.wikimedia.org/wikipedia/commons/e/ea/K-means_convergence.gif\")
\n\n\nk-Means is popular because it is a **reliable** and **fast** algorithm. The main downside is that it assumes the clusters are **spherical**, which is not always the case. ","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Baseline\nmodel_km = KMeans(n_clusters=7, random_state=0)\npreds_km = model_km.fit_predict(scaled_data)","metadata":{"execution":{"iopub.status.busy":"2022-07-11T15:51:19.768861Z","iopub.execute_input":"2022-07-11T15:51:19.769814Z","iopub.status.idle":"2022-07-11T15:51:31.761952Z","shell.execute_reply.started":"2022-07-11T15:51:19.769777Z","shell.execute_reply":"2022-07-11T15:51:31.760745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"k-Means doesn't end up doing so well on the leaderboard (it scores around ARI=0.23), so we will need more sophisticated models that allow for non-spherical clusters.\n\n
\n
\n\nA Gaussian Mixture Model (GMM) is a clustering algorithm that assumes the data is made up **combination/mixture** of several (multivariate) **Gaussian/Normal distributions**. In contrast to k-Means, data points are assigned a **probability** of belonging to each cluster as opposed to being assigned a single cluster. These probabilities are worked out using the **Expectation-Maximization** (EM) algorithm, which estimates the **parameters** (mean & covariance matrix) of Gaussian distributions using an **iterative Maximum Likelihood Estimation** (MLE) method. At the end of the learning process, data points are assigned to a cluster by choosing the cluster with the **highest probability** out of all of them. \n\nNote that the sklearn implementation can put restrictions on the type of covariance matrices learned. In particular, they can be **tied**, where all clusters have the same covariance matrix, **diagonal**, where all covariance matrices are diagonal, **spherical**, where the resulting clusters are spherical and **full**, where there are no restrictions on the covariance matrices.\n\n6.3 Gaussian Mixture Model (GMM)
\n![](\"https://c.tenor.com/i1rNMdaKd7MAAAAC/gaussian-mixture-models-em-method-math.gif\")
\n\n\n
\n
\n\nA Bayesian Gaussian Mixture Model (BGMM) is very similar to a GMM. It uses the **same assumptions** on the data and the same approach to finding the clusters. The only difference is in the learning algorithm. Instead of Maximum Likelihood Estimation, BGMMs use **Variational Bayesian Estimation**. \n\nIn contrast to traditional frameworks, the Bayesian approach views parameters as **random variables** rather than fixed unknown quantities. It estimates the distribution of these parameters by sampling from the posterior distribution using methods like **Markov Chain Monte Carlo** (MCMC).\n\n","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Baseline\nmodel_bgmm = BayesianGaussianMixture(n_components=7, covariance_type='full', max_iter=100, n_init=5, init_params='random', random_state=0)\npreds_bgmm = model_bgmm.fit_predict(scaled_data_crop)","metadata":{"execution":{"iopub.status.busy":"2022-07-18T11:00:06.262251Z","iopub.execute_input":"2022-07-18T11:00:06.262700Z","iopub.status.idle":"2022-07-18T11:00:16.718097Z","shell.execute_reply.started":"2022-07-18T11:00:06.262662Z","shell.execute_reply":"2022-07-18T11:00:16.716819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Bayesian GMM performs the **best** out of all of the models we've tried so far (it scores around ARI=0.59 on the public leaderboard). It does take **longer** to run though.\n\n# 7. Visualise predictions\n\n6.4 Bayesian Gaussian Mixture Model (Bayesian GMM)
\n\n
","metadata":{}},{"cell_type":"code","source":"# Countplot\nplt.figure(figsize=(10,4))\nsns.countplot(x=preds_bgmm)\nplt.title('Predicted clusters')","metadata":{"execution":{"iopub.status.busy":"2022-07-18T11:00:20.866069Z","iopub.execute_input":"2022-07-18T11:00:20.866465Z","iopub.status.idle":"2022-07-18T11:00:21.096942Z","shell.execute_reply.started":"2022-07-18T11:00:20.866434Z","shell.execute_reply":"2022-07-18T11:00:21.095818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"7.1 Label distribution
\n\n
\n\nTo get insight into the underlying distributions, we can plot cluster-wise histograms of each of the remaining features.\n\n**Continuous: (f_22 to f_28)**","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/ambrosm/tpsjul22-gaussian-mixture-cluster-analysis\nfig, axs = plt.subplots(2, 4, figsize=(20, 7))\naxs = axs.ravel()\nfloat_columns = ['f_22','f_23','f_24','f_25','f_26','f_27','f_28']\ny=preds_bgmm\nfor ax, f in zip(axs, float_columns):\n for i in range(7):\n h, edges = np.histogram(data[f][y == i], bins=np.linspace(-5, 5, 26))\n ax.plot((edges[:-1] + edges[1:]) / 2, h, label=f\"Cluster {i}\", lw=3)\n ax.set_title(f)\n#axs[-2].axis('off')\naxs[-1].axis('off')\nplt.suptitle('Histograms of continuous features by cluster', y=1.02, fontsize=28)\nfig.tight_layout(h_pad=1.0, w_pad=0.5)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-18T11:00:25.709234Z","iopub.execute_input":"2022-07-18T11:00:25.709660Z","iopub.status.idle":"2022-07-18T11:00:27.185430Z","shell.execute_reply.started":"2022-07-18T11:00:25.709623Z","shell.execute_reply":"2022-07-18T11:00:27.184075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Discrete: (f_07 to f_13)**","metadata":{}},{"cell_type":"code","source":"# From https://www.kaggle.com/code/ambrosm/tpsjul22-gaussian-mixture-cluster-analysis\nprop_cycle = plt.rcParams['axes.prop_cycle']\n\nfig, axs = plt.subplots(2, 4, figsize=(20, 7))\naxs = axs.ravel()\nint_columns = [col for col in data.columns if data[col].dtype == 'int']\nfor ax, f in zip(axs, int_columns):\n for i in range(7):\n uv, uc = np.unique(data[f][y == i], return_counts=True)\n ax.plot(uv, uc, alpha=1, color=prop_cycle.by_key()['color'][i % 10], lw=3)\n ax.set_title(f)\n #ax.legend()\naxs[-1].axis('off')\nplt.suptitle('Histograms of discrete features by cluster', y=1.02, fontsize=28)\nfig.tight_layout(h_pad=1.0, w_pad=0.5)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-18T11:00:34.753815Z","iopub.execute_input":"2022-07-18T11:00:34.754266Z","iopub.status.idle":"2022-07-18T11:00:36.147505Z","shell.execute_reply.started":"2022-07-18T11:00:34.754228Z","shell.execute_reply":"2022-07-18T11:00:36.144824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"7.2 Cluster distributions
\n\n
\n\n**Principal Component Analysis (PCA)** was the first dimensionality reduction technique discovered (by Karl Pearson - yes, the guy from Pearson's correlation coefficient) and dates back to as early as **1901**. It is very popular because it is **fast**, **easy to implement** and **easy to interpret**. \n\nPCA works by finding a low dimensional subspace that **maximises the variance** of the data in that subspace and performing a **linear projection**. This basically means the data will be as **spread out** as possible, without changing the relationship between the data points. This allows us to find patterns in dimensions we can visualised.","metadata":{}},{"cell_type":"code","source":"%%time\n\n# PCA\npca = PCA(n_components=3)\ncomponents = pca.fit_transform(scaled_data_crop)\n\n# 3D scatterplot\nfig = px.scatter_3d(\n components, x=0, y=1, z=2, color=preds_bgmm, size=0.1*np.ones(len(scaled_data_crop)), opacity = 1,\n title='PCA plot in 3D',\n labels={'0': 'PC 1', '1': 'PC 2', '2': 'PC 3'},\n width=650, height=500\n)\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-11T16:02:41.682499Z","iopub.execute_input":"2022-07-11T16:02:41.683321Z","iopub.status.idle":"2022-07-11T16:02:44.189533Z","shell.execute_reply.started":"2022-07-11T16:02:41.683276Z","shell.execute_reply":"2022-07-11T16:02:44.186587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Explained variance** shows how much of the variance/spread of the data is captured in each dimension, i.e. how **important** each additional **principal component** is to the original data representation.","metadata":{}},{"cell_type":"code","source":"# PCA\npca_var = PCA()\npca_var.fit(scaled_data_crop)\n\n# Plot\nplt.figure(figsize=(10,5))\nxi = np.arange(1,1+scaled_data_crop.shape[1], step=1)\nyi = np.cumsum(pca_var.explained_variance_ratio_)\nplt.plot(xi, yi, marker='o', linestyle='--', color='b')\n\n# Aesthetics\nplt.ylim(0.0,1.1)\nplt.xlabel('Number of Components')\nplt.xticks(np.arange(1,1+scaled_data_crop.shape[1], step=1))\nplt.ylabel('Cumulative variance (%)')\nplt.title('Explained variance by each component')\nplt.axhline(y=1, color='r', linestyle='-')\nplt.gca().xaxis.grid(False)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-11T16:02:44.191235Z","iopub.execute_input":"2022-07-11T16:02:44.191704Z","iopub.status.idle":"2022-07-11T16:02:44.740703Z","shell.execute_reply.started":"2022-07-11T16:02:44.191669Z","shell.execute_reply":"2022-07-11T16:02:44.739478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"7.3 Principle Component Analysis (PCA)
\n\n
\n\n**t-SNE** (pronounced tiz-knee) stands for **t-distributed Stochastic Neighbor Embedding** and was proposed much more recently by Laurens van der Maaten and Geoffrey Hinton in their [2008 paper](https://www.jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf). \nThis works in a similar way to PCA but has some key differences:\n* Firstly, this is a **stochastic method**. So if you run multiple t-SNE plots on the same dataset it can look different.\n* Another difference is that this is an **iterative method**. It works by repeatedly moving datapoints closer or further away from each other depending on how 'similar' they are. \n* The new representation is **non-linear**. This makes it harder to interpret but it can be very effective at 'unravelling' highly non-linear data.\n\nThe main downside to t-SNE is that is **very slow** compared to the other dimensionality techniques. This is because it makes calculations on a pair-wise basis, which does not scale well with large datasets.","metadata":{}},{"cell_type":"code","source":"%%time\n\n# PCA\ntsne = TSNE(n_components=3)\ncomponents = tsne.fit_transform(scaled_data_crop.iloc[:5000,:])\n\n# 3D scatterplot\nfig = px.scatter_3d(\n components, x=0, y=1, z=2, color=preds_bgmm[:5000], size=0.1*np.ones(len(scaled_data_crop.iloc[:5000,:])), opacity = 1,\n title='t-SNE plot in 3D',\n labels={'0': 'comp. 1', '1': 'comp. 2', '2': 'comp. 3'},\n width=650, height=500\n)\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-11T16:02:44.741987Z","iopub.execute_input":"2022-07-11T16:02:44.742327Z","iopub.status.idle":"2022-07-11T16:05:41.446053Z","shell.execute_reply.started":"2022-07-11T16:02:44.742297Z","shell.execute_reply":"2022-07-11T16:05:41.444890Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Even with just 5% of the data it still takes several minutes to run.","metadata":{}},{"cell_type":"markdown","source":"7.4 t-SNE
\n\n
\n\n**UMAP**, which stands for **Uniform Manifold Approximation and Projection** was proposed by Leland McInnes, John Healy and James Melville in their [2018 paper](http://gobie.csb.pitt.edu/SML/umap.pdf).\n\nIt is similar to t-SNE in that it learns a non-linear mapping that preserves clusters but its main advantage is that it is **significantly faster**. It also tends to do better at preserving **global structure** of the data compared to t-SNE. \n\nReference: https://pair-code.github.io/understanding-umap/ ","metadata":{}},{"cell_type":"code","source":"%%time\n\n# UMAP\num = umap.UMAP(n_components=3)\ncomponents_umap = um.fit_transform(scaled_data_crop)\n\n# 3D scatterplot\nfig = px.scatter_3d(\n components_umap, x=0, y=1, z=2, color=preds_bgmm, size=0.1*np.ones(len(scaled_data_crop)), opacity = 1,\n title='UMAP plot in 3D',\n labels={'0': 'comp. 1', '1': 'comp. 2', '2': 'comp. 3'},\n width=650, height=500\n)\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-11T16:05:41.447627Z","iopub.execute_input":"2022-07-11T16:05:41.448085Z","iopub.status.idle":"2022-07-11T16:07:18.817464Z","shell.execute_reply.started":"2022-07-11T16:05:41.448050Z","shell.execute_reply":"2022-07-11T16:07:18.815453Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"UMAP's **connectivity plot** is a weighted graph that gives insight into the representation of the embedding. It basically shows which connections were most important when creating the projection. ","metadata":{}},{"cell_type":"code","source":"%%time\n\n# Connectivity plot\num = umap.UMAP()\nX_fit = um.fit(scaled_data_crop)\numap.plot.connectivity(X_fit, show_points=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-07-11T16:07:18.819101Z","iopub.execute_input":"2022-07-11T16:07:18.819484Z","iopub.status.idle":"2022-07-11T16:08:59.424365Z","shell.execute_reply.started":"2022-07-11T16:07:18.819450Z","shell.execute_reply":"2022-07-11T16:08:59.422858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 8. Submission","metadata":{}},{"cell_type":"code","source":"sub = pd.read_csv('../input/tabular-playground-series-jul-2022/sample_submission.csv')\nsub['Predicted'] = preds_bgmm\nsub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-07-11T16:08:59.426176Z","iopub.execute_input":"2022-07-11T16:08:59.426659Z","iopub.status.idle":"2022-07-11T16:08:59.626172Z","shell.execute_reply.started":"2022-07-11T16:08:59.426610Z","shell.execute_reply":"2022-07-11T16:08:59.625299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 9. References\n\n* [Notebok: Understading the competition metric: Adjusted rand Index](https://www.kaggle.com/competitions/tabular-playground-series-jul-2022/discussion/334534) by [Towhidul Tonmoy](https://www.kaggle.com/towhidultonmoy).\n* [Paper: Comparing Partitions](https://link.springer.com/content/pdf/10.1007/BF01908075.pdf) by Hubert and Arabie, 1985.\n* [Paper: MEASURING AGREEMENT WHEN TWO OBSERVERS CLASSIFY PEOPLE INTO CATEGORIES NOT DEFINED IN ADVANCE](https://bpspsychub.onlinelibrary.wiley.com/doi/abs/10.1111/j.2044-8317.1974.tb00535.x?casa_token=n4nz2gru0rYAAAAA:gl9377Ehcq4eoUyIfpxJ5Z6CotPTTk8QTptWEC1rfkM3Wp9MEk16UDfr-4CRWNIECO7otF-Fp00ux1I) by Brennan and Light, 1974.\n* [Notebook: Outliers+EDA+Clustering Tutorial](https://www.kaggle.com/code/javigallego/outliers-eda-clustering-tutorial) by [Francisco Javier Gallego\n](https://www.kaggle.com/javigallego).\n* [Article: Details of the Adjusted Rand index and Clustering algorithms](https://faculty.washington.edu/kayee/pca/supp.pdf) by Yeung and Ruzzo, 2001. \n* [Notebook: TPS Jul 22 ADVANCED + 2% SOL](https://www.kaggle.com/code/kartushovdanil/tps-jul-22-advanced-2-sol) by [Torch me](https://www.kaggle.com/kartushovdanil).\n* [Discussion: Visualizing the seven clusters](https://www.kaggle.com/competitions/tabular-playground-series-jul-2022/discussion/334808) by [AmbrosM](https://www.kaggle.com/ambrosm).\n* [Article: POISSON DISPERSION TEST](https://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/poisdisp.htm) by NIST.\n* [Paper: A Test for the Poisson Distribution](http://www-stat.wharton.upenn.edu/~lbrown/Papers/2002c%20A%20new%20test%20for%20the%20Poisson%20distribution%20(with%20L.%20H.%20Zhao).pdf) by Lawrence D. Brown and Linda H. Zhao, 2002.\n* [Lecture notes: Stats 200 Introduction to Statistical Inference](https://artowen.su.domains/courses/200/lec11.pdf) by Art B. Owen, 2018.\n\n# 10. Appendix\n\n7.5 UMAP
\n\n
\n\nIf you're like me, then you will be **curious** as to where **Hubert's formula** for ARI even comes from. I spent several hours reading the **original** papers (Hubert & Arabie 1985, Brennan & Light 1974) to try to understand the derivation and I wanted to **share** my findings. Feel free to **skip** this section though, as it is going to be **very mathematical**. \n\nARI derivation
\n\n\nWe start by extending the definitions we had before.\n\n* $a$ = # pairs in the **same** cluster in the prediction and the **same** cluster in the ground truth.\n* $b$ = # pairs in **different** clusters in the prediction and **different** clusters in the ground truth.\n* $c$ = # pairs in **different** clusters in the prediction and in the **same** cluster in the ground truth.\n* $d$ = # pairs in the **same** cluster in the prediction and **different** clusters in the ground truth.\n\nNote: since $a+b+c+d= {n \\choose 2}$ (total number of pairs), the RI can also be expressed as $RI=\\frac{a+b}{a+b+c+d}$.\n\n
\n\nTo derive **Hubert's formula**, we need to calcuate $a, b, c, d$ in terms of the **contingency table**. This will allow us to calculate the expected RI later. Recall $n_{ij}=|X_i \\cap Y_j|$, so ${n_{ij} \\choose 2}$ gives the number of pairs of points that are both in $X_i$ and $Y_j$. By definition, these pairs will be in the **same** cluster ($X_i$) in the prediction and the **same** cluster ($Y_j$) in the ground truth they count towards $a$. To actually calculate $a$, all we need to do is **sum** these combinations over all the **intersections**. That is,\n\n$$a = \\sum_{i,j} {n_{ij} \\choose 2}$$\n\n
\n\nWe can calculate $b$ using a **trick** once we have $a,c,d$. For now, we focus on probably the most **difficult** part, i.e. calculating $c$. We need to work out how many pairs there are that are in the **same** cluster in the ground truth and in **different** clusters in the prediction. This will involve a technique called **double counting** (instead of working out $c$, we work out $2c$ then divide the answer by 2). \n\nConsider taking one point from $n_{ij}$, i.e. in $X_i \\cap Y_j$. Where can the other point be so that it counts towards $c$? It needs to be in the **same** cluster as in the ground truth, so it needs to be in $X_i$ but in a **different** cluster to the prediction, so it can't be in $Y_j$. Graphically, the other point needs to be from the **same row** but **different column** in the contingency table. \n\n
![](\"https://i.postimg.cc/k5v5Wb3m/Contingency-table.jpg\")
\n\n\nThe argument for $d$ is almost identical because of **symmetry**. All you have to do is swap the $i$'s with the $j$'s and swap the $a_i$'s with the $b_j$'s. If you want to test your understanding, this would be a really good **exercise** to derive it for yourself.\n\n$$\nd = \\frac12 \\left(\\sum_{j} b_{j}^2 - \\sum_{ij} n_{ij}^2 \\right)\n$$\n\n
\n\nAnd now we just need to work out $b$. There is probably a clever combinatorial argument for this, but I'm going to **cheat** a little bit by using the fact that $a+b+c+d = {n \\choose 2}$ to make things easier.\n\n$$\n\\begin{align*}\nb &= {n \\choose 2} - a - c - d \\newline\n&= {n \\choose 2} - \\sum_{i,j} {n_{ij} \\choose 2} - \\frac12 \\left(\\sum_{i} a_{i}^2 - \\sum_{ij} n_{ij}^2 \\right) - \\frac12 \\left(\\sum_{j} b_{j}^2 - \\sum_{ij} n_{ij}^2 \\right) \\newline\n&= \\frac12 \\left(n^2 - n - \\sum_{i,j} n_{ij}^2 + \\sum_{i,j} n_{ij} - \\left(\\sum_{i} a_{i}^2 + \\sum_{j} b_{j}^2 \\right) + 2 \\sum_{i,j} n_{ij}^2 \\right)\n\\end{align*}\n$$\nUsing the fact that $\\sum_{ij} n_{ij}=n$ (there are n data points in total) and after some simplifying we get\n\n$$\nb= \\frac12 \\left(n^2 + \\sum_{i,j} n_{ij}^2 - \\left(\\sum_{i} a_{i}^2 + \\sum_{j} b_{j}^2 \\right) \\right)\n$$\n\n
\n\nSo now we have **all the pieces** to work out the Rand Index using the contingency table. You will see why this helpful when we work out its expected value. For now, let's work out the number of pairs in **agreement**, i.e. $a+b$. Unfortunately, this is a bit of **messy algebra** to get it into the form we want. I've tried to break it down as much as I can. \n\n$$\n\\begin{align*}\na+b &= \\sum_{i,j} {n_{ij} \\choose 2} + \\frac12 \\left(n^2 + \\sum_{i,j} n_{ij}^2 - \\left(\\sum_{i} a_{i}^2 + \\sum_{j} b_{j}^2 \\right) \\right) \\newline\n&= \\sum_{i,j} {n_{ij} \\choose 2} + {n \\choose 2} + \\frac{n}{2} + \\frac12 \\sum_{ij} n_{ij}^2 - \\frac12 \\left( \\sum_{i} a_i^2 + \\sum_{j} b_j^2 \\right) \\newline\n&= \\sum_{i,j} {n_{ij} \\choose 2} + {n \\choose 2} + \\sum_{i,j} {n_{ij} \\choose 2} + \\sum_{i,j} n_{ij} - \\frac12 \\left( \\sum_{i} a_i^2 + \\sum_{j} b_j^2 \\right) \\newline\n&= {n \\choose 2} + 2 \\sum_{i,j} {n_{ij} \\choose 2} + n - \\frac12 \\left( \\sum_{i} a_i^2 + \\sum_{j} b_j^2 \\right) \\newline\n&= {n \\choose 2} + 2 \\sum_{i,j} {n_{ij} \\choose 2} - \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right)\n\\end{align*}\n$$\n\nTherefore,\n\n$$\nRI = \\frac{a+b}{{n \\choose 2}} = 1 + 2 \\sum_{i,j} {n_{ij} \\choose 2} \\Big/ {n \\choose 2} - \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right) \\Big/ {n \\choose 2}\n$$\n\n
\n\nThis is where the **magic** starts. Assuming a generalised **hypergeometric distribution** (i.e. the **probabilistic** way to model clusters) Hubert in 1977, showed the following very useful result.\n\n$$\n\\mathbb{E} \\left( \\sum_{ij} {n_{ij} \\choose 2} \\right)\n = \\sum_{i} {a_i \\choose 2} \\sum_{j} {b_j \\choose 2} \\Big/ {n \\choose 2}\n$$\n \nUnfortunately, the paper where this was proved is **not open access** so I couldn't verify the details. Regardless, it allows us to calculate the expected value of the Rand Index. Combining the above two formulas gives:\n\n$$\n\\mathbb{E} (RI) = 1 + 2 \\sum_{i} {a_i \\choose 2} \\sum_{j} {b_j \\choose 2} \\Big/ {n \\choose 2}^2 - \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right) \\Big/ {n \\choose 2}\n$$\n\nAlong with the fact that $max(RI)=1$, we have everything we need to calculate the Adjusted Rand Index.\n\n
\n\nRecall the formula,\n\n$$\nARI = \\frac{RI - \\mathbb{E}(RI)}{\\max(RI)-\\mathbb{E}(RI)}\n$$\n\nFor the numerator, we have\n\n$$\n\\begin{align*}\nRI - \\mathbb{E}(RI) &= 1 + 2 \\sum_{i,j} {n_{ij} \\choose 2} \\Big/ {n \\choose 2} - \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right) \\Big/ {n \\choose 2} \\newline\n&- 1 - 2 \\sum_{i} {a_i \\choose 2} \\sum_{j} {b_j \\choose 2} \\Big/ {n \\choose 2}^2 + \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right) \\Big/ {n \\choose 2} \\newline\n&= 2 \\sum_{i,j} {n_{ij} \\choose 2} \\Big/ {n \\choose 2} - 2 \\sum_{i} {a_i \\choose 2} \\sum_{j} {b_j \\choose 2} \\Big/ {n \\choose 2}^2\n\\end{align*}\n$$\n\nFor the denominator, we have\n\n$$\n\\begin{align*}\n\\max(RI)-\\mathbb{E}(RI) &= 1 - 1 - 2 \\sum_{i} {a_i \\choose 2} \\sum_{j} {b_j \\choose 2} \\Big/ {n \\choose 2}^2 + \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right) \\Big/ {n \\choose 2} \\newline\n&= \\left( \\sum_{i} {a_i \\choose 2} + \\sum_{j} {b_j \\choose 2} \\right) \\Big/ {n \\choose 2} - 2 \\sum_{i} {a_i \\choose 2} \\sum_{j} {b_j \\choose 2} \\Big/ {n \\choose 2}^2\n\\end{align*}\n$$\n\nFinally, putting the numerator on top and denominator on the bottom and dividing both by the common factor of $2 / {n \\choose 2}$, we get\n\n$$\n\\large ARI = \\frac{ \\left. \\sum_{ij} \\binom{n_{ij}}{2} - \\left[\\sum_i \\binom{a_i}{2} \\sum_j \\binom{b_j}{2}\\right] \\right/ \\binom{n}{2} }{ \\left. \\frac{1}{2} \\left[\\sum_i \\binom{a_i}{2} + \\sum_j \\binom{b_j}{2}\\right] - \\left[\\sum_i \\binom{a_i}{2} \\sum_j \\binom{b_j}{2}\\right] \\right/ \\binom{n}{2} }\n$$\n\nas required.","metadata":{}}]}