A clustering algorithm for overlapping Gaussian mixtures

Abstract

Gaussian mixture models (GMM) are widely used as a probabilistic model for density estimation for multivariate data and as an unsupervised clustering algorithm to provide a soft clustering to the available data. The GMM relies on the expectation-maximization algorithm for maximizing the likelihood. A new approach is proposed in the present work, which depends on Approximate Bayesian Computation and aims not only to estimate the population parameters but also to assign each observation to a specific subpopulation. The performance, in terms of mean estimates for the parameters of each subpopulation, number of observations assigned to each cluster, accuracy, and other performance metrics, of the new approach is compared with the expectation-maximization algorithm for the GMM and the K-means algorithm under several challenging simulation scenarios and is applied to a real data set.

KEYWORDS:

• Overlapping subpopulations
• soft clustering
• Gaussian mixture models
• approximate Bayesian computation
• kernel density

1 Introduction

Cluster analysis is one of the most frequently used processes in unsupervised learning. Clustering aims to determine homogeneous subgroups in a data set by identifying similarities between data points’ characteristics (Ezugwu et al. 2022) and has many applications in different fields, like computer science, data science and statistics, healthcare, urban development and transportation, environmental and technological impact and innovations, and pattern recognition and image analysis (Liu, Sun, and Wang 2006; Jyoti and Singh 2011; Bandyopadhyay et al. 2014; Ding et al. 2023).

Clustering algorithms may be classified into two main categories, namely (a) as hard and exclusive membership methods and (b) as soft (fuzzy), overlapping techniques. In the latter, probability model-based approaches play an important role (Bouveyron and Brunet-Saumard 2014; Bouveyron et al. 2019). Among the most popular probability, model-based approaches are the mixture models, especially the Gaussian Mixture Models (GMM), which not only provide a useful tool to model populations consisting of two or more subpopulations but also a powerful approach to revealing underlying structures in bi- and multi-modal data (Quiñones-Grueiro et al. 2019). For fitting a mixture model, the expectation and maximization (EM) algorithm (Dempster, Laird, and Rubin 1977a) is usually used due to the complexity of the likelihood function.

GMMs and mixture models in general, estimate the probability that a data point belongs to each cluster and thus can be used as soft classifiers whenever clusters overlap. However, the phenomenon of cluster overlapping is not mathematically well characterized, especially in multivariate cases (Sun and Wang 2011), and can be, in many cases, a really challenging task to model correctly.

In this paper, we aim at clustering based on GMM in the case of not-well-separated, overlapped clusters, and we propose a likelihood-free approach based on the Approximate Bayesian Computation (ABC) algorithm to identify and estimate the parameters of the underlying subpopulations using a discrepancy measure focusing on the kernel density estimations. Kernel density approaches are often used in classification and clustering problems (see for example Taylor 1997; Hinneburg and Gabriel 2007; Anderson 2009) since a cluster (or a class, in the classification problems) can be easily identified/defined by a local maximum of the estimated density function. Moreover, we focus on Gaussian mixtures since they are not only suitable for many clustering algorithms but also sone of the most frequently used models in many applications (Yang, Lai, and Lin 2012; McLachlan and Rathnayake 2014; Ahmed et al. 2019). Moreover, we compare the proposed method with the classical EM approach for the GMM and the K-means clustering algorithm, which is one of the most popular nonparametric, hard, and exclusive membership methods.

The rest of the paper is organized as follows: In Section 2, we briefly review the related work and present some motivating examples. In Section 3, the proposed method is presented in detail. In Sections 4 and 5 a simulation study is conducted and the proposed method is applied to a real data set, respectively. Finally, some concluding remarks are given in the final section.

2 Preliminaries and related work

2.1 Gaussian mixture models as a clustering algorithm

Finite mixture models are often adopted to describe heterogeneity among individuals or groups in different applications and data (Economou and Caroni 2009; Economou 2013; Punzo 2019; McLachlan, Lee, and Rathnayake 2019). Under a finite mixture model, the population is assumed to consist of k subpopulations or groups with (probability) density function given by (1) $$f(\mathit{x}|\lambda )=\sum _{j=1}^k{p}_j{f}_j(\mathit{x}|{\lambda }_{\mathit{j}})$$(1) where x is a d-dimensional vector of real value random variables, and $p_j,j=1,2,\dots ,k$ are the subpopulation proportions with ${\sum }_{j=1}^k{p}_j=1,f_j(\mathbf{x}|{\lambda }_{\mathbf{j}})$ the (probability) density function of the jth subpopulation, and ${\lambda }_{\mathbf{j}}$ the corresponding parameter vector. Usually, the f_js have the same functional form, as in the rest of the paper, but different parameters.

One of the most frequently used mixture models is the so-called GMM. A GMM is a weighted sum of k component Gaussian densities given by $$f(\mathit{x}|\lambda )=\sum _{j=1}^k{p}_{j}g(\mathit{x}|{\mu }_j, {\mathbf{\Sigma }}_j),$$where $g(\mathbf{x}|{\mu }_j, {\mathbf{\Sigma }}_j),j=1,\dots ,k$ are the k component Gaussian densities given by $$g(\mathit{x}|{\mu }_j, {\mathbf{\Sigma }}_j)=\frac{1}{{(2\pi )}^{d/2}|{\Sigma }_j{|}^{1/2}} \hspace{0.17em}\exp \hspace{0.17em}\left\{-\frac{1}{2}(\mathit{x}-{\mu }_j)\prime {\Sigma }_j^{-1}(\mathit{x}-{\mu }_j)\right\},$$ with mean vector μ_j and covariance matrix Σ_j and λ is the vector of all the parameters. The jth d-dimensional Gaussian distribution is denoted in the following as $MN({\mu }_j,{\mathbf{\Sigma }}_j)$.

There are several techniques available for estimating the parameters of a GMM but the most popular is the expectation-maximization (EM) algorithm (Dempster, Laird, and Rubin 1977b). During the EM algorithm, the Expectation (E-step) and the Maximization (M-step) steps are repeated until convergence is achieved. In the bth iteration, the E-step basically updates the probabilities $P({\mathbf{x}}_i\in j|{\mathbf{x}}_i,{\mu }_j^{(b-1)}, {\mathbf{\Sigma }}_j^{(b-1)},p_j^{(b-1)},{\lambda }^{(b-1)})$, i.e., the probabilities that a data point ${\mathbf{x}}_i,i=1,\dots ,n$ belongs to component $j=1,\dots ,k$, using the Bayes’ rule $$\begin{array}{c}P({\mathit{x}}_i\in j|{\mathit{x}}_i,{\mu }_j^{(b-1)}, {\mathbf{\Sigma }}_j^{(b-1)},p_j^{(b-1)},{\lambda }^{(b-1)})\\ \hspace{1em}=\frac{g\left({\mathit{x}}_i|{\mu }_j^{(b-1)}, {\mathbf{\Sigma }}_j^{(b-1)}\right)p_j^{(b-1)}}{f({\mathit{x}}_i|{\lambda }^{(b-1)})}\end{array}$$where the parameters ${\mu }_j^{(b-1)},{\mathbf{\Sigma }}_j^{(b-1)}$, $p_j^{(b-1)}$, and ${\lambda }^{(b-1)}$ are the current values, obtained by the previous iteration, of the parameters. The probabilities $P({\mathbf{x}}_i\in j|{\mathbf{x}}_i,{\mu }_j^{(b-1)}, {\mathbf{\Sigma }}_j^{(b-1)},p_j^{(b-1)},{\lambda }^{(b-1)})$ serve as soft labels that a data point ${\mathbf{x}}_i,i=1,\dots ,n$ belongs to component $j=1,\dots ,k$.

In the M-step, the new values of the parameters are determined by maximizing with respect to the parameter and by using the ${\gamma }_{{\mathbf{x}}_i}(j)=P({\mathbf{x}}_i\in j|{\mathbf{x}}_i,{\mu }_j^{(b-1)}, {\mathbf{\Sigma }}_j^{(b-1)},p_j^{(b-1)},{\lambda }^{(b-1)})$, the expected value of the complete log-likelihood given by $$\sum _{i=1}^n\sum _{j=1}^k{\gamma }_{{\mathit{x}}_i}(j)\left(\hspace{0.17em}\log \hspace{0.17em}(p_j)+\hspace{0.17em}\log \hspace{0.17em}\left(g\left({\mathit{x}}_i|{\mu }_j, {\mathbf{\Sigma }}_j\right)\right)\right).$$

Remark 1.

As already mentioned, the EM algorithm provides soft labels for each data point. Although, the final, fitted model can also be used as a hard clustering procedure by assigning each data point to the component yielding the highest posterior probability. Moreover, some modified approaches have been suggested based on penalized maximum likelihood estimation to obtain improved and more robust clustering algorithms based on GMM (Yang, Lai, and Lin 2012; He and Ho 2019).

2.2 K-means

The K-means algorithm is often considered as one of the most frequently used clustering algorithms (Berkhin 2006). K-means is an iterative algorithm that divides the available observations into k clusters by (step a) assigning each observation to the cluster with the nearest centroid (mean), and (step b) revising each cluster centroid (mean) based on the points in that cluster, and (step c) repeating steps (a) and (b) until convergence, i.e. the labels of all the points do not change, is achieved.

Although K-means is a hard clustering method and does not make any assumptions for the distribution of the underlying populations, as the GMM algorithm does, it can still be viewed as a limiting case of a GMM (Barber 2012) or as a special case of them when truncated variational EM approximations are applied to Gaussian Mixture Models (GMM) with isotropic Gaussians (Lücke and Forster 2019). As a result, K-means is good at capturing (hyper)ball-shaped structures (Kłopotek, Wierzchoń, and Kłopotek 2020) and for that reason, K-means is also considered, for comparison purposes, along with the GMM in the rest of the paper.

Remark 2.

An important extension of the K-means algorithm is the Fuzzy K-means. The Fuzzy K-means algorithm is a soft clustering algorithm that relies on the same principles as the K-means algorithm. In the Fuzzy K-means algorithm, each data point is not assigned to a single cluster, but membership scores are assigned to each data point for each cluster based on the distance between the cluster centroid and the data point (Bezdek 1981; Ferraro 2021).

2.3 Simulated examples

GMM and K-means are powerful clustering algorithms but can still underperform in certain situations in which overlapping subpopulations are presented. More specifically, the following examples present cases in which every subpopulation is relatively close to at least one of the other subpopulations. The closeness of the subpopulation is measured based on Matusita’s measure of overlapping coefficient (Matusita 1955) and in all cases, all subpopulations present a small to moderate degree of similarity with at least one of the rest of the subpopulations. In all the following simulated examples, the weakness of GMM and K-means is apparent. The K-means usually works well for unimodal clusters and well-separated groups, so it is expected to underperform in such situations. Nevertheless, K-means is also used for comparison purposes, since it is one of the most frequently used algorithms for clustering. On the other hand, GMM is more suitable for overlapping subpopulations, even if it is not expected to work perfectly whenever the subpopulations are not well separated (Celeux 1985). Of course, as the sample size increases, the GMM is expected to present more acceptable behavior as long as the mixing weights are not too extreme. For small sample sizes, such as the ones presented next, whenever overlapping subpopulations are presented, both GMM and K-means, as will be demonstrated, underperform.

Example 1.

The first data set was generated from a three-component Gaussian mixture with $$\begin{array}{c}{\mu }_1=(-2,-2{)}^{\prime },{\Sigma }_1=(\begin{array}{cc}2& 1\\ 1& 2\end{array})\\ {\mu }_2=(0,0{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}2& -1\\ -1& 2\end{array})\\ {\mu }_3=(2,2{)}^{\prime },{\Sigma }_3=(\begin{array}{cc}2& 1\\ 1& 2\end{array})\end{array}$$and $p_i=1/3,i=1,2,3$. The Matusita’s measure of overlapping coefficients ${\rho }_{i,j}$ between the ith and the jth distribution are ${\rho }_{1,2}={\rho }_{2,3}=0.53$ and ${\rho }_{1,3}=0.26$. These values indicate that the second component presents moderate overlap with each of the other distributions. The simulated sample consisted of 33 observations from the first subpopulation, 46 observations from the second subpopulation, and 41 observations from the third subpopulation. In Figure 1(a) the contour plots of the three subpopulations given the label of the observations are given. Plotting lines in each contour plot define points with common values of the corresponding probability density function (pdf), estimated using the observed sample given the label of each data point. Lines closer to the center of each contour plot indicate high values of the corresponding pdf, i.e. ranges with higher probability. The number of lines in each contour plot and their closeness also represent how rapidly these values change; fewer and closer lines indicate a distribution with smaller standard deviations. In Figure 1(b) the clusters obtained by the K-means are presented, while in Figure 1(c) the clusters identified by GMM are depicted. Figure 1(d) presents the contour plots of the three clusters identified by the proposed method, presented in the following section. To avoid unnecessary repetition, we defer the discussion of the details of this plot to the next section.

Fig. 1 (a) The contour plots of the three subpopulations given the label of the observations generated under the distribution scheme of Example 1. (b) The clusters obtained by the K-means. (c) The clusters obtained by the GMM. (d) The clusters obtained by the proposed algorithm.

From Figures 1(b) and (c) it is clear that both the K-means and the GMM fail to capture the true underlying distribution of the three clusters (as demonstrated, despite the underrepresentation in that particular sample of the first subpopulation, in Figure 1(a)). More specifically, K-means seems to manage to correctly identify the centroid of each cluster but fails to express correctly the variability (i.e., the variance) of each cluster, mainly due to the misclassification of some points (see Figure 1(b)). On the other hand, GMM fails completely to identify the centroids and the variance-covariance structure of all three clusters (see Figure 1(c)).

Example 2.

The second data set was again generated from a three-component Gaussian mixture, but this time with $$\begin{array}{c}{\mu }_1=(-4,0{)}^{\prime },{\Sigma }_1=(\begin{array}{cc}1& 0\\ 0& 1\end{array})\\ {\mu }_2=(0,0{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}4& 0\\ 0& 4\end{array})\\ {\mu }_3=(4,0{)}^{\prime },{\Sigma }_3=(\begin{array}{cc}1& 0\\ 0& 1\end{array})\end{array}$$and $p_i=1/3,i=1,2,3$. In this case the Matusita’s measure of overlapping coefficients between these distributions are ${\rho }_{1,2}={\rho }_{2,3}=0.36$ and ${\rho }_{1,3}=3.35·{10}^{-4}$ indicating that the second component presents moderate to small overlapping with each of the other distributions since their mean vectors are further away, compared to the first example, from the mean vector of the second subpopulation.

Thirty-nine (39) observations were generated from the first subpopulation, 34 observations from the second subpopulation, and 46 observations from the third subpopulation. Their scatterplot, along with the contour plots of the three subpopulations estimated by this sample, given the label of the observations, is presented in Figure 2(a). Clearly, there is significant overlap of the three clusters. The first (left, red) and the third (right, blue) subpopulations seem to be almost complete in the volume covered by the second subpopulation (middle, green). This feature is missed by both the K-means (which was expected since K-means do not allow overlapping clusters) and the GMM, see Figure 2(b) and (c), respectively. Especially, the clusters identified by the GMM clustering algorithm seem to behave more or less as the ones identified by the K-means, i.e., they present a very small to no overlap. This is mainly due to the noticeable elongation that the middle cluster presents.

Fig. 2 (a) The contour plots of the three subpopulations given the label of the observations generated under the distribution scheme of Example 2. (b) The clusters obtained by the K-means. (c) The clusters obtained by the GMM. (d) The clusters obtained by the proposed algorithm.

Example 3.

The third data set was generated from a four-component Gaussian mixture with $$\begin{array}{c}{\mu }_1=(-3,0{)}^{\prime },{\Sigma }_1=(\begin{array}{cc}2& -1\\ -1& 2\end{array})\\ {\mu }_2=(0,1{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}2& 0\\ 0& 2\end{array})\\ {\mu }_3=(3,-2{)}^{\prime },{\Sigma }_3={\Sigma }_1\\ {\mu }_4=(5,1{)}^{\prime },{\Sigma }_4={\Sigma }_1\end{array}$$

and $p_i=1/4,i=1,2,3,4$. The Matusita’s measure of overlapping coefficients between these distributions are $$\begin{array}{c}{\rho }_{1,2}=0.45{\rho }_{1,3}=9.69·{10}^{-2}\\ {\rho }_{1,4}=2.28·{10}^{-2}{\rho }_{2,3}=0.39\\ {\rho }_{2,4}=0.18{\rho }_{3,4}=0.21\end{array}$$

Thirty-six (36), 46, 43, and 35 observations were generated from the first, second, third, and fourth subpopulations, respectively. It is clear that three of the four clusters, corresponding to the last three subpopulations, were not identified correctly from both the K-means and the GMM algorithm (see Figure 3(b) and (c)). Especially, the cluster of the second subpopulation (green) was completely missed by the GMM since it shrank and covered a very small volume. Moreover, the 4th cluster of the last subpopulation seems to cover a smaller volume based on the K-means algorithm (see Figure 3(b)) compared with that indicated in Figure 3(a). The same cluster based on the GMM seems to cover a larger volume with a small positive correlation (see Figure 3(c)).

Fig. 3 (a) The contour plots of the three subpopulations given the label of the observations generated under the distribution scheme of Example 3. (b) The clusters obtained by the K-means. (c) The clusters obtained by the GMM. (d) The clusters obtained by the proposed algorithm.

3 Proposed methodology

The proposed clustering method relies heavily on the so-called Approximate Bayesian computation (ABC) method. ABC is a likelihood-free technique that is used to perform inference in cases in which the evaluation of the likelihood function is computationally costly or even infeasible (Rubin 1984; Diggle and Gratton 1984). It is used in many different applications and fields (see for example Economou et al. 2020; Park and Kwon 2022). The key idea behind ABC is the production of artificial data sets using some prior distribution and comparing them with the observed sample. If a produced data set is “close enough” to the observed data set, then the parameters used to generate this set are not-rejected. The procedure is repeated until M, a relatively large integer, and not-rejected artificial data sets are obtained. Then, the M not-rejected parameters are summarized to describe the posterior distribution of the parameters.

The ABC method for the proposed clustering algorithm is described in the following steps:

1. (a) Adopt a prior distribution for the cluster membership score for each data point.

(b) Generate a data partition by assigning each data point to one of the k clusters by generating an integer from 1 to k. Each integer represents the k clusters, and the integers are sampled using the prior distribution adopted for each data point.

1. Using the nj members allocated to the jth cluster, estimate the sample mean, ${\widehat{\mu }}_j$, and the sample variance-covariance matrix, ${\widehat{\mathbf{\Sigma }}}_j$.

2. For $j=1,\dots ,k$ generate nj observations from the $MN({\widehat{\mu }}_j,{\widehat{\mathbf{\Sigma }}}_j)$.

3. Compute the discrepancy between the artificial $D^{*}$ and the observed data set D using a suitable defined discrepancy measure d.

4. Retain (i) the partition, (ii) the sample mean, ${\widehat{\mu }}_j$ and (iii) the sample variance-covariance matrix, ${\widehat{\mathbf{\Sigma }}}_j$ only if $d<ϵ$, for some $ϵ>0$. Go to step 1 (b) and repeat until M not-rejected samples are obtained.

The above algorithm describes the procedure to obtain artificial data sets that are similar to the observed data set. In the sequel, a few essential details for applying steps 1(a), 1(b), 4, and 5 of the previous algorithm in practice are discussed.

1. (a) A flat, prior categorical distribution for the cluster membership score can be adopted for each data point. If prior knowledge is available for the cluster sizes, a categorical distribution with different parameters that reflects the mixing proportions of the underlying subpopulations can be used. In this case, special attention should be given to the labeling of each cluster. A simple approach is to label the clusters with respect to the order of their mean values with respect to one variable. For example, if one believes that the middle cluster in Figure 2 is twice the size of the other two clusters, then the categorical distribution with parameters 1/4, 2/4, and 1/4 should be used.

Another approach, that is actually adopted in the rest of the paper, is to use a data-driven prior distribution by initially applying a soft, fuzzy clustering algorithm to obtain a fuzzy n × k matrix with the membership score of each data point with respect to the k clusters. From the examples, presented in Section 2.3 it seems that the GMM fails, at least in some cases, to correctly identify the centroid of the clusters. On the other hand, the K-means algorithm seems to present more robust behavior regarding the centroids of the clusters. Given that Fuzzy K-means relies on the same principles as the K-means algorithm, using the membership score matrix obtained by the Fuzzy K-means as a data-driven prior distribution is recommended when applying the proposed algorithm.

(b) A similar procedure can be found in the Stochastic EM (SEM) algorithms proposed by (Celeux 1985). The SEM algorithm incorporates a stochastic step, the S-step, between the E- and M-steps of the classical EM algorithm in which each observation is assigned to each subpopulation according to a multinomial distribution with parameters equal to the a posteriori probabilities of the E-step. This procedure does not allow the EM algorithm to be trapped in local optima (Novais and Faria 2021). Although the idea seems similar the proposed algorithm does not make use of this partition to compute the M-step of the EM algorithm but to generate an artificial data set and compare it with the observed one.

4. A discrepancy measure can be defined based on the integrated squared error given by $$T={\int }_x{\int }_y{\left(f_D(x,y)-f_{D^{*}}(x,y)\right)}^{2}d\mathit{ydx}$$between the kernel density estimates $f_D(x,y)$ and $f_{D^{*}}(x,y)$ of the observed D and the artificial $D^{*}$ sample. T is asymptotically normally distributed with a known mean and standard deviation under the null hypothesis, i.e. under the hypothesis that both samples share the same density (for more details see Duong, Goud, and Schauer 2012) and so an approximate p-value can be calculated based on the standardized version of T, d.

5. Small absolute values of d, or equivalently large p-values, assure that the retained partitions of the data points are more likely to resemble the true, unknown labeling of the data points, and as a consequence, the estimated sample mean, ${\widehat{\mu }}_j$, and sample variance-covariance matrix, ${\widehat{\mathbf{\Sigma }}}_j$ to be closer to their true values. To ensure that only highly similar artificial data sets are retained, we propose the adoption of a very large threshold, at least for the two dimension case and for relatively small sample sizes, of p-value, for example, 0.75.

It is of note that the M not-rejected parameters can be used not only to summarize the posterior distribution of the parameters of the underlying Gaussian distributions but also to assign each data point to a specific subpopulation, i.e., cluster. This can be achieved by simply adopting the majority voting method.

Remark 3.

In mixture models and clustering, “label switching” is always a problem that needs to be resolved before a majority voting method can be applied. A common approach to dealing with this problem is to set an identifiability constraint on the parameter space regarding the order of the mean values in one of the dimensions of the identified clusters. In cases where these mean values are relatively close, the order of the mean values of the identified cluster in a second dimension can also be used to label the identified clusters. This approach is also used in the simulations and in the application presented in Sections 4 and 5, respectively.

To make the proposed algorithm clearer, the first of the two following examples demonstrate analytically the procedure, while the second reports the results of the applications of the proposed algorithm to the data sets of Examples 1–3 presented in Section 2.3.

Example 4.

In Example 1, a sample consisting of 33 observations from the first subpopulation, 46 observations from the second subpopulation, and 41 observations from the third subpopulation was generated and analyzed with K-means and GMM. The scatterplot of this sample, using the labels, which in what follows are assumed unknown, is presented in Figure 4(a). The estimated kernel density plot for the observed data is presented in Figure 4(d). Regarding the label switching problem mentioned in Remark 3, the three identified clusters in Figure 2(d) can be labeled as 1, 2, and 3 based on the order of their mean values with respect to the horizontal axis (from left to right).

Fig. 4 (a) The scatterplot of the three subpopulations given the label of the observations generated under the distribution scheme of Example 1. (b) A rejected data partition. (c) A not-rejected data partition. (d) The kernel density plot for the observed data. (e) The kernel density plot for the rejected artificial data set generated using the partition presented in subfigure (b). (f) The kernel density plot for the not-rejected artificial data set generated using the partition presented in subfigure (c).

From the analysis in Example 1 it is clear that the GMM underperforms with respect to the K-means which justifies the use of the Fuzzy K-means algorithm to obtain a fuzzy 120 × 3 matrix with the membership score of each one of the 120 data points with respect to the k = 3 clusters (step 1(a) of the algorithm). Based on these matrices, several different data partitions can be generated (step 1(b) of the algorithm). Two of them are presented in Figure 4(b) and (c). Obviously, none of them is correct 100% but the scatterplot presented in Figure 4(b), which, as we will see, corresponds to the rejected partition, is less similar to the partition presented in scatterplot of Figure 4(a) (see, for example, the two red points (marked in red cycles) in the up-right part of the plot).

Next (step 2 of the algorithm) the ${\widehat{\mu }}_j$, and the ${\widehat{\mathbf{\Sigma }}}_j$, j = 1, 2, 3 were calculated for each of the two partitions. Using the corresponding sample size n_j of each partition, random observations were generated from the $MN({\widehat{\mu }}_j,{\widehat{\mathbf{\Sigma }}}_j)$ distributions in order to produce artificial data sets (step 3 of the algorithm). The kernel density plots for these artificial data sets are presented in Figure 4(e) and (f). From these plots, it is clear that the kernel density plot in Figure 4(f) resembles more that obtained by the observed data. This is also demonstrated by the p-value obtained by testing the hypothesis that these samples share the same density as the observed sample. More specifically, the p-values are 0.05384673 and 0.7563502 for the artificial data sets of Figure 4(e) and (f), respectively (step 5 of the algorithm). Since the p-value of the first artificial data set is smaller than the threshold 0.75 set in the 5th step of the proposed plot, the used data partition and the corresponding ${\widehat{\mu }}_j$, and ${\widehat{\mathbf{\Sigma }}}_j$, j = 1, 2, 3 are dropped. The opposite occurs for the second artificial data set.

Examples 1–3 continued. In Section 2.3 we have deferred the discussion of the details regarding the last plots, labeled as plots (d), in Figures 1–3. These plots present the contour plots of the clusters identified by the proposed method using 101 not-rejected samples. In Examples 1 and 3 it is clear that the proposed method managed to capture the true underlying distribution of the subpopulations. In the second example, the performance of the proposed method is not as great as before, but it is still acceptable and outperforms that of GMM and K-means.

4 Simulation study

Examples 1–3 indicate that the performance of the proposed algorithm outperforms that of GMM and K-means. Although, to evaluate the performance of the proposed algorithm in cases in which the K-Means and GMM are likely to underperform, and compare its performance with the GMM and K-means clustering algorithms a simulation study was carried out. More specifically, 1000 simulated samples were generated for each of the following scenarios: $$\begin{array}{c}S\mathit{cenario}1:{\mu }_1=(-4,0{)}^{\prime },{\Sigma }_1=(\begin{array}{cc}1& 0\\ 0& 1\end{array})\\ {\mu }_2=(0,0{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}2& 0\\ 0& 2\end{array})\\ {\mu }_3=(4,0{)}^{\prime },{\Sigma }_3=(\begin{array}{cc}1& 0\\ 1& 1\end{array})\end{array}$$ $$\begin{array}{c}S\mathit{cenario}2:{\mu }_1=(-4,0{)}^{\prime },{\Sigma }_1=(\begin{array}{cc}1& 0\\ 0& 1\end{array})\\ {\mu }_2=(0,0{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}4& 0\\ 0& 4\end{array})\\ {\mu }_3=(4,0{)}^{\prime },{\Sigma }_3=(\begin{array}{cc}1& 0\\ 1& 1\end{array})\end{array}$$ $$\begin{array}{c}S\mathit{cenario}3:{\mu }_1=(-4,0{)}^{\prime },{\Sigma }_1=(\begin{array}{cc}1& 0\\ 0& 1\end{array})\\ {\mu }_2=(-1,-1{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}2& 0\\ 0& 2\end{array})\\ {\mu }_3=(1,1{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}2& 0\\ 0& 2\end{array})\\ {\mu }_4=(4,0{)}^{\prime },{\Sigma }_2=(\begin{array}{cc}1& 0\\ 0& 1\end{array})\end{array}$$

To make the comparison more straightforward, the sample size of each subpopulation was fixed at 40 observations. Additionally, to investigate the performance of the proposed algorithm with non-comparable cluster sizes, the following sample sizes $(n_1,n_2,n_3,n_4)=(40,50,30,40)$ were also used for the four subpopulations in the third scenario.

For each data set, the K-means, the GMM (as a hard clustering algorithm), and the proposed algorithm (with M = 101) were applied. The value of M was selected through a small, initial simulation study in which the membership status of each point at each cluster was carefully monitored. From the simulation, it was clear that for the majority of the cases (and points), the membership status of each point was unchanged from M larger than 80–90 not-rejected samples. Furthermore, the value 101 was selected for M since 101 is a prime number, which makes it highly improbable that a point is assigned to 2 or more clusters the exactly same maximum number of times. However, in case this unlikely event occurs, the membership status of such a point may be determined by simply increasing the value of M The evaluation of the performance and the comparison of the three approaches were made in terms of the mean (based on the 1000 simulated samples for each scenario)

• estimates for the parameters of each subpopulation,

• number of observations assigned to each cluster

• accuracy

• Cohen’s Kappa statistic and

• sensitivity, specificity, recall, and F1 score per cluster.

In the upper tabular in Tables 1–4 the mean values of the posterior distributions of the parameters of each subpopulation (Note that μ_x, μ_y in the tables denote the mean values of X and Y, respectively, while σ₁₁, σ₁₂, and σ₂₂ denote the (1,1)th, (1,2)th and (2,2)th elements of the corresponding variance-covariance matrix) and the mean number of observations assigned to each cluster for the 1000 simulated samples in each case are reported. Their standard deviations are reported in the second tabular, while the bottom two tabulars report the performance evaluation metrics for each scenario. In all cases, the best values (closer to the true value, minimum standard deviation, or higher value of the performance evaluation metric) are indicated in bold.

Table 1 Summary results (Mean values-upper tabular, standard deviations-second tabular, and performance evaluation metrics-two bottom tabular) for the 1000 simulated samples from Scenario 1.

Mean values
	1st cluster						2nd cluster						3rd cluster
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
True	40	−4	0	1	0	1	40	0	0	4	0	4	40	4	0	1	0	1
K-means	42.38	–3.93	–0.00	0.99	–0.00	1.05	35.06	–0.01	0.01	1.04	0.01	1.86	42.56	3.93	–0.00	1.00	0.00	1.05
GMM	44.47	–3.81	0.00	1.11	0.01	1.01	31.18	0.01	0.01	0.73	0.01	1.93	44.35	3.82	0.00	1.11	–0.01	1.00
New	42.38	–3.94	–0.00	0.97	0.01	1.04	35.47	0.02	0.00	1.07	0.01	1.91	42.14	3.96	0.00	0.97	–0.01	1.03
St. Dev. of the estimates
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
K-means	3.53	0.21	0.21	0.28	0.17	0.24	3.16	0.49	0.31	0.25	0.23	0.48	3.75	0.24	0.18	0.30	0.18	0.25
GMM	5.72	0.34	0.18	0.55	0.22	0.25	8.25	0.53	0.37	0.41	0.39	0.66	5.99	0.34	0.17	0.53	0.22	0.25
New	2.38	0.19	0.17	0.19	0.16	0.24	2.93	0.29	0.28	0.19	0.24	0.47	2.41	0.19	0.16	0.19	0.17	0.24

Table

Performance evaluation metrics
	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1
K-means	0.97	0.82	0.97	0.96	0.97	0.95	0.97	0.82	0.97	0.94	0.87	0.94
GMM	0.98	0.73	0.98	0.93	0.98	0.93	0.98	0.73	0.98	0.93	0.81	0.93
New	0.98	0.84	0.97	0.96	0.97	0.96	0.98	0.84	0.97	0.95	0.89	0.95

Table 2 Summary results (Mean values-upper tabular, standard deviations-second tabular, and performance evaluation metrics-two bottom tabular) for the 1000 simulated samples from Scenario 2.

Mean values
	1st cluster						2nd cluster						3rd cluster
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
True	40	–4	0	1	0	1	40	0	0	2	0	2	40	4	0	1	0	1
K-means	46.18	–3.84	0.01	1.15	0.01	1.41	27.51	–0.00	–0.01	1.28	0.01	2.98	46.31	3.83	0.00	1.18	–0.01	1.41
GMM	47.79	–3.73	-0.00	1.23	–0.00	1.31	24.58	0.01	0.02	0.87	–0.02	3.50	47.63	3.73	0.00	1.26	0.00	1.31
New	45.52	–3.90	0.00	1.01	0.00	1.30	29.36	0.04	–0.00	1.30	0.02	3.76	45.12	3.92	–0.00	1.01	–0.00	1.26
St. Dev. of the estimates
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
K-means	5.27	0.26	0.28	0.40	0.35	0.41	4.88	0.85	0.92	0.46	0.38	1.03	5.35	0.28	0.29	0.44	0.37	0.41
GMM	8.07	0.41	0.33	0.66	0.36	0.46	9.72	1.01	0.89	0.60	0.73	1.69	8.79	0.46	0.32	0.69	0.39	0.50
New	3.07	0.19	0.19	0.20	0.21	0.34	3.54	0.45	0.58	0.30	0.47	1.08	3.08	0.20	0.20	0.20	0.21	0.33

Table

Performance evaluation metrics
	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1
K-means	0.96	0.61	0.96	0.90	0.96	0.90	0.96	0.61	0.96	0.89	0.72	0.89
GMM	0.97	0.54	0.96	0.89	0.96	0.89	0.97	0.54	0.96	0.88	0.66	0.88
New	0.97	0.68	0.97	0.92	0.97	0.92	0.97	0.68	0.97	0.91	0.78	0.91

Table 3 Summary results (Mean values-upper tabular, standard deviations-second tabular, and performance evaluation metrics-two bottom tabular) for the 1000 simulated samples from Scenario 3 with equal cluster sizes.

Mean values
	1st cluster						2nd cluster						3rd cluster						4th cluster
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
True	40	–4	0	1	0	1	40	–1	–1	2	0	2	40	1	1	2	0	2	40	4	0	1	0	1
K-means	45.21	–3.90	–0.04	0.97	0.02	1.06	35.03	–0.82	–1.19	1.10	0.02	1.21	34.93	0.85	1.21	1.09	0.01	1.22	44.82	3.92	0.03	0.97	0.03	1.05
GMM	50.46	–3.68	–0.12	1.17	–0.09	1.09	29.29	–0.78	–1.03	0.66	0.07	1.46	29.52	0.74	1.02	0.66	0.06	1.52	50.73	3.67	0.11	1.20	–0.08	1.09
New	43.82	–3.97	0.01	0.89	0.07	1.04	36.34	–0.88	–1.28	1.09	–0.07	1.34	36.29	0.91	1.29	1.09	–0.07	1.33	43.56	3.98	–0.02	0.90	0.08	1.03
St. Dev. of the estimates
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
K-means	6.36	0.24	0.27	0.25	0.22	0.28	7.39	0.68	0.78	0.28	0.22	0.41	7.38	0.70	0.76	0.28	0.22	0.42	6.71	0.26	0.27	0.24	0.21	0.28
GMM	12.12	0.46	0.29	0.63	0.34	0.37	18.12	0.88	0.97	0.61	0.50	1.10	17.85	0.89	0.98	0.59	0.53	1.27	13.17	0.49	0.29	0.66	0.33	0.38
New	3.64	0.20	0.19	0.18	0.16	0.24	5.11	0.37	0.45	0.26	0.22	0.39	4.99	0.40	0.42	0.26	0.22	0.36	3.97	0.20	0.20	0.18	0.15	0.24

Table

Performance evaluation metrics
	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1
K-means	0.94	0.63	0.63	0.93	0.94	0.92	0.92	0.94	0.94	0.63	0.63	0.93	0.88	0.67	0.67	0.88
GMM	0.95	0.47	0.48	0.94	0.89	0.91	0.91	0.89	0.95	0.47	0.48	0.94	0.84	0.51	0.52	0.83
New	0.94	0.69	0.69	0.94	0.95	0.93	0.93	0.95	0.94	0.69	0.69	0.94	0.90	0.72	0.72	0.90

Table 4 Summary results (Mean values-upper tabular, standard deviations-second tabular, and performance evaluation metrics-two bottom tabular) for the 1000 simulated samples from Scenario 3 with non-equal cluster sizes.

Mean values
	1st cluster						2nd cluster						3rd cluster						4th cluster
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
True	40	–4	0	1	0	1	50	–1	–1	2	0	2	30	1	1	2	0	2	40	4	0	1	0	1
K-means	46.01	–3.90	–0.03	0.95	0.04	1.04	37.19	–0.92	–1.38	1.09	0.01	1.14	32.51	0.63	1.01	1.10	0.03	1.23	44.29	3.92	0.02	0.99	0.02	1.08
GMM	51.39	–3.68	–0.11	1.13	–0.06	1.08	32.37	–0.88	–1.14	0.70	0.10	1.40	26.99	0.58	0.80	0.59	0.08	1.50	49.26	3.68	0.10	1.22	–0.10	1.12
New	44.68	–3.95	0.01	0.90	0.09	1.03	38.47	–1.00	–1.43	1.08	–0.07	1.25	34.11	0.73	1.02	1.09	–0.05	1.40	42.75	3.99	0.00	0.89	0.05	1.06
St. Dev. of the estimates
	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22	n	μx	μy	σ11	σ12	σ22
K-means	6.27	0.23	0.26	0.23	0.19	0.28	7.40	0.63	0.73	0.26	0.20	0.39	7.30	0.70	0.80	0.28	0.23	0.41	5.97	0.24	0.27	0.25	0.22	0.28
GMM	14.15	0.50	0.29	0.62	0.31	0.39	19.69	0.87	0.97	0.59	0.48	1.06	17.39	0.96	1.06	0.50	0.50	1.27	10.99	0.45	0.28	0.62	0.33	0.36
New	3.85	0.20	0.19	0.18	0.16	0.24	5.28	0.37	0.51	0.27	0.21	0.36	4.94	0.46	0.54	0.26	0.24	0.43	3.95	0.20	0.19	0.18	0.17	0.26

Table

Performance evaluation metrics
	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1	Sensitivity	Specificity	Recall	F1
K-means	0.94	0.61	0.64	0.94	0.93	0.94	0.90	0.95	0.94	0.61	0.64	0.94	0.87	0.69	0.62	0.89
GMM	0.93	0.48	0.46	0.95	0.88	0.92	0.90	0.91	0.93	0.48	0.46	0.95	0.82	0.54	0.46	0.86
New	0.94	0.65	0.71	0.94	0.94	0.94	0.90	0.96	0.94	0.65	0.71	0.94	0.89	0.73	0.67	0.91

In all the scenarios studied, the proposed method presents on average a slightly better behavior with respect to the number of observations assigned to each cluster and the overall performance evaluation metrics (Accuracy & Kappa). Moreover, the mean values of the estimates for the parameters of each subpopulation for the proposed method were, in the majority of the cases, closer to their true values. In the cases in which the mean values of the estimates for the parameters obtained by the proposed method were not the best, the mean estimates of the parameters were still close to the true ones. Furthermore, it seems that the estimates for the parameters obtained by the proposed method seem to be less variable since they present, in almost all the cases, the smallest standard deviation.

Regarding the performance evaluation metrics for each cluster in every scenario, it seems that, with the exception of sensitivity, the proposed method presents the best performance with respect to all metrics in all the simulated scenarios among the three methods tested. Comparing Scenarios 1 and 2, which differ only in the variance of the middle subpopulation, it is clear that all the methods present better behavior in the first scenario. This indicates that when there is a significant overlap of the subpopulations, all the methods face larger problems in assigning the observations to the correct cluster. Scenario 3 differs from scenarios 1 and 2 in the number of subpopulations. From the results, it seems, as may be expected, that as the number of overlapped subpopulations increases, the values of the performance evaluation metrics decrease, demonstrating the more complex setup of the population. Moreover, the case of non-equal clusters (Table 4) seems to suffer more compared with the equal size clusters (Table 3), mainly due to the mean number of observations assigned to the most extreme, with respect to their number of observations, clusters.

Remark 4.

The proposed algorithm seems to perform better compared to K-means and GMM in the above-studied cases, i.e. in cases of relatively small sample sizes and overlapping subpopulations. This is probably due to the fact that the proposed algorithm not only permits the overlapping of the subpopulations (which K-means do not) but also assigns the observed points to clusters that produce similar data sets as the observed one and not simply to the multivariate normal components that maximize the component posterior probability, given the data, as GMM do. Very often, in cases of relatively small sample sizes and overlapping subpopulations, the log-likelihood surface is littered with saddle points, plateaux, ridges, and sub-optimal maxima which often lead the GMM to wrong assignments and incorrect inference. On the other hand, the ABC does not make use of the underlying likelihood function and focuses only on the distance between the observed and the artificial data sets generated under the assumed model. As a result, the proposed algorithm, which relies on the ABC algorithm, does not suffer as much as the EM in these cases and can be used to carry out not only the Bayesian posterior inference, through the lens of density estimation but also the assignment of each data point to a cluster.

Remark 5.

The proposed algorithm assumes a known number of clusters in the data. If the number of clusters is misspecified, then the performance of the algorithm may be seriously affected. From a parallel simulation study that assumed an incorrect number of clusters for the above-studied scenarios, it seems that whenever the number of clusters is underestimated, the acceptance rate may be small, and a large number of artificial data sets need to be generated before M not-rejected samples are obtained. On the other hand, if the number of clusters is overestimated, the number of rejected samples becomes very small. While this may be an indicator that the acceptance rate can be used as a measure of selecting the optimal number of clusters, the simulations did not highlight a clear, solid rule that can be applied in all cases. Moreover, it is worth mentioning that due to the overlapping of the subpopulations studied, the algorithm does not manage to assign zero (or almost zero) observations to the extra cluster(s), but splits the overlapping region into two or more clusters.

Remark 6.

The average acceptance rates (standard deviation, minimum, and maximum values) for the four simulation scenarios were 0.2014 (0.1263, 0.0002, 0.6392), 0.0928 (0.0770, 0.0002, 0.4529), 0.1503 (0.1015, 0.0005, 0.5050), and 0.2560 (0.1819, 0.0103, 0.6159), respectively. This is equivalent, for example, to the fact that 501.49 samples were required on average in the first scenario to obtain the desired 101 not-rejected samples. From the results, it is clear that there is such a large heterogeneity in the acceptance rate across the simulated samples that no clear pattern can be identified.

5 Application

For illustration purposes, the proposed method is applied to a flea beetle data set that can be found in Lubischew (1962) and Hand et al. (1993). The data consists of 74 beetles in the genus Flea Beetle. Chaetocnema, which contains three species, namely the heikertingeri (number of observations: 31, subpopulation 1), the heptapotamica (number of observations: 22, subpopulation 2), and the concinna (number of observations: 21, subpopulation 3). These three subpopulations and the corresponding label of each observation are used to evaluate the performance of the proposed algorithm. For each observation, two characteristics of the beetle were measured: the maximal width (x) of the aedeagus in the forepart (in microns) and the front angle (y) of the aedeagus (1 unit = 7.5 degrees).

For comparison purposes, the K-means and the GMM are also applied.

5.1 Clustering algorithms

In Figure 5(a) the contour plots of the three subpopulations given their true labels are given. The red contour corresponds to the heikertingeri, the green to heptapotamica, and the blue to the concinna species. In Figure 1(b) the clusters obtained by the K-means are presented, while in Figure 1(c) the clusters identified by GMM are depicted. Figure 1(d) presents the contour plots of the three clusters identified by the proposed method. Again, the proposed approach to resolve the label-switching problem presented in Remark 3 was adopted and the three identified clusters were labeled as 1, 2, and 3 based on the order of their mean values with respect to the horizontal axis, i.e., the clusters were labeled with respect to their mean value of the maximal width (x) of the aedeagus in the forepart (in microns). In Table 5 the confusion tables of the three approaches are also presented.

Fig. 5 (a) The contour plots of the three subpopulations given the labels of the observations for the Flea Beetle data set. (b) The clusters obtained by the K-means. (c) The clusters obtained by the GMM. (d) The clusters obtained by the proposed algorithm.

Table 5 Confusion tables for the flea beetle data set using the three clustering algorithms.

		cluster	heikertingeri	heptapotamica	concinna
	K-means	1	30	3	0
		2	1	19	9
		3	0	0	12
	GMM	1	20	0	0
		2	10	0	0
		3	1	22	21
Proposed	algorithm	1	30	1	0
		2	1	21	5
		3	0	0	16

From Figure 1(b) and 1(c) and the first two tabulars in Table 5 it is clear that both the K-means and the GMM fail to capture the true underlying distribution of the three clusters/subpopulations. More specifically, K-means failed mainly to correctly identify the concinna species since K-means assigned 9 of them (those with large values of the front angle of the aedeagus) to the second cluster that consists mainly of heptapotamica beetles. On the other hand, GMM allocates all the heptapotamica and concinna species to the third cluster and, at the same time, splits the heikertingeri into two clusters. It is worth noticing that the second cluster identified by the GMM covers a very small volume and seems to be shrunken between the other two clusters, a behavior that was observed in the simulated examples presented in Section 2.3.

The proposed algorithm managed to identify almost perfectly the first two species (heikertingeri and heptapotamica) and seems to suffer less, compared to the K-Means, in identifying the third species (concinna). Although it is of note that the proposed method still fails to capture the true structure of the variance-covariance matrices of the underlying subpopulations, especially that of the heptapotamica as can be seen also in Table 6, in which the descriptive statistics for the three identified clusters by the three clustering algorithms are presented. Although it is worth mentioning that in all the cases, the proposed algorithm identified the centroid of each subpopulation better. Finally, it is worth mentioning that 2534 samples (acceptance rate: 0.03986) were generated until 101 not-rejected samples were obtained.

Table 6 Descriptive statistics for the three identified clusters by the used clustering algorithms.

	Cluster 1
	n	$\overline{x}$	$\overline{y}$	σ11	σ12	σ22
True labels	31	124.65	14.29	21.37	–0.33	1.22
K-means	34	124.88	13.94	20.55	–2.07	2.75
GMM	20	121.86	14.34	7.22	–0.166	1.20
Prop. alg.	31	124.52	14.23	19.02	–0.44	1.47
	Cluster 2
	n	$\overline{x}$	$\overline{y}$	σ11	σ12	σ22
True labels	22	138.27	10.09	17.16	–0.50	0.94
K-means	29	139.79	11.55	9.96	0.51	5.04
GMM	10	129.63	14.33	1.19	0.74	0.97
Prop. alg.	27	138.52	10.96	12.77	–0.54	3.81
	Cluster 3
	n	$\overline{x}$	$\overline{y}$	σ11	σ12	σ22
True labels	21	146.19	14.10	31.66	–0.97	0.79
K-means	12	150.08	13.83	10.08	0.94	0.70
GMM	44	141.87	12.07	39.51	6.82	4.72
Prop. alg.	16	148.50	14.00	14.50	–0.135	0.75

The first line of each tabular represents the corresponding values using the true label of each observation.

6 Conclusions and discussion

Clustering, whether hard or soft, is one of the most frequently used techniques to get an intuition about the structure of the data. Among the techniques appearing in the literature, the most frequently used are the K-means and the GMM. The GMM is also one of the most widely used probabilistic models for density estimation. Although, K-means and GMM are powerful clustering algorithms, they can still underperform in certain situations, as presented in Section 2. Especially in some situations, GMM tends to shrink certain subpopulations/clusters resulting in incorrect partitions of the data.

To overcome these problems, a new algorithm was proposed based on the Approximate Bayesian Computation. This method seems to be able to deal better with overlapping clusters and provide more accurate estimates for the parameters of the underlying subpopulations. More research is required on how the proposed algorithm may be applied in high-dimensional data and/or in cases in which some clusters may be well-separated from the rest of the clusters. For example, for the case of a well-separated cluster from the rest of the clusters other discrepancy measures can also be studied to compare the samples, like the kernel maximum mean discrepancy measure (Gretton et al. 2006). For high-dimensional data, since density estimation is usually a challenging task, (a) the determination of the threshold value with respect to the dimension of the data should be also considered and (b) the use of different concepts, as the Friedman-Rafsky test (Friedman and Rafsky 1979) and the nearest neighbor test (Henze 1988), can be worth studying to better compare the samples.

The future aspects of using other discrepancy measures are mentioned but it remains unclear if this is something to consider for the scenarios of focus in the manuscript or also for the higher dimensional data scenario.

Moreover, a study on hyper-parameters tuning may be necessary to achieve an even better performance of the proposed algorithm. In the present study, the default multivariate plug-in bandwidth selection with unconstrained pilot matrices was implemented in R (for more details, see Wand and Jones 1994; Duong and Hazelton 2003; Chacón and Duong 2010). Additionally, it will be interesting to develop, as future work, a technique, based on the proposed algorithm, to determine also the optimal number of the underlying subpopulations/clusters in a given data set. Another interesting point for further research could be the combination of the proposed algorithm with other clustering algorithms to develop an ensemble clustering algorithm and/or generate posterior samples after the best partition to provide more insights for a given sample (see for example Rastelli and Friel 2018).

Acknowledgments

The author thanks all the reviewers and the editor for their valuable comments and suggestions.

Disclosure statement

No potential conflict of interest was reported by the author.

Table

	Accuracy	Kappa
K-means	0.92	0.88
GMM	0.89	0.84
New	0.93	0.89

The best values for each case are indicated in bold.

Table

	Accuracy	Kappa
K-means	0.85	0.77
GMM	0.82	0.74
New	0.88	0.81

The best values for each case are indicated in bold.

Table

	Accuracy	Kappa
K-means	0.78	0.71
GMM	0.71	0.61
New	0.81	0.75

The best values for each case are indicated in bold.

Table

	Accuracy	Kappa
K-means	0.78	0.71
GMM	0.71	0.61
New	0.80	0.74

The best values for each case are indicated in bold.