{'fraction_non_alphanumeric': 0.029914651090254346, 'fraction_numerical': 0.013752232965661951, 'mean_word_length': 4.692978208232446, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian model averaging in model-based clustering and density estimation

We propose Bayesian model averaging (BMA) as a method for postprocessing the results of model-based clustering. Given a number of competing models, appropriate model summaries are averaged, using the posterior model probabilities, instead of being taken from a single"best"model. We demonstrate the use of BMA in model-based clustering for a number of datasets. We show that BMA provides a useful summary of the clustering of observations while taking model uncertainty into account. Further, we show that BMA in conjunction with model-based clustering gives a competitive method for density estimation in a multivariate setting. Applying BMA in the model-based context is fast and can give enhanced modeling performance.


Introduction
Model-based clustering methods are based on the assumption that the population can be modeled using the finite mixture model (e.g. Banfield and Raftery, 1993;Celeux and Govaert, 1995;Fraley and Raftery, 2002). Within this paradigm, it is assumed that the data come from G subpopulations, which correspond to the mixture components, and within each subpopulation the data are modeled using a single parametric component distribution. The most common finite mixture model that is used for model-based clustering is the finite normal mixture, but many alternatives exist (e.g. McLachlan and Peel, 2000;Lee and McLachlan, 2013b).
Model selection is an intrinsic part of model-based clustering. In particular, the number of clusters (component densities), G, is unknown and a number of competing choices for component densities may also be under consideration. Each combination of component density and number of clusters can be viewed as a separate model, and a model selection approach can be used to select both at the same time (e.g. Fraley and Raftery, 2002). In most implementations of model-based clustering, the "best" model is chosen by using some criterion and clustering is based on the "best" single model. A number of methods have been proposed for selecting the "best" model including choosing the model with the highest Bayesian Information Criterion (BIC) (Schwarz, 1978) or the highest Integrated Complete Likelihood (ICL) (Biernacki et al., 2000).
The approach of reporting the results of model-based clustering based on a single model ignores the uncertainty that arises from the model selection. Consequently, the uncertainty about quantities of interest may be underestimated. We propose basing the results of model-based clustering on a combination of the results from all candidate models rather than on those from a single model. We propose taking a weighted average of the model summaries, where the weights are approximate posterior model probabilities. Thus we propose using Bayesian Model Averaging (BMA) (Hoeting et al., 1999) within the model-based clustering paradigm.
To obtain valid inference using BMA, the quantity of interest should have the same meaning in each model under consideration. In model-based clustering, the quantities of interest must have the same meaning for all values of G and must be invariant to the labeling of the clusters in the finite mixture model. This is because finite mixture models are identifiable only up to permutations of the cluster labels. Here we focus on inference about the clustering consensus matrix. This has the same meaning for all values of G and is invariant to the cluster labeling.
We also consider using model-based clustering as a method for multivariate density estimation, following Fraley and Raftery (2002). In this case, the estimated density has the same meaning in all models, so we again use the posterior model probabilities as weights to offer an alternative density estimation procedure to multivariate kernel density estimation (e.g., Scott, 1992;Duong, 2007) or to model-based clustering density estimation methods based on a single model (Fraley and Raftery, 2002).
The paper is organized as follows. In Section 2, we briefly review the modelbased clustering paradigm. In Section 3, we outline how the BMA approach can be used to deal with model uncertainty. In Section 4, we describe how we can create a matrix for each model that is invariant to the number and labeling of the clusters. In Section 5, we provide an assessment of clustering performance for BMA in conjunction with model-based clustering. In Section 6, we describe the background for density estimation and introduce a framework for combining BMA and model-based clustering for multivariate density estimation. We illustrate this with a number of simulations. We conclude, in Section 7, with a discussion of related work and suggest future directions.

Model-based clustering
Model-based clustering (Banfield and Raftery, 1993;Celeux and Govaert, 1995;Raftery, 2002, 2007) is used for clustering data into groups, where the number of groups G is typically unknown. Here we focus on model-based clustering based on the finite normal mixture model as described in Fraley and Raftery (2002).
We assume that there are G clusters, where each cluster g arises with probability τ g . Data within each cluster follow a normal distribution with cluster-specific mean µ g and covariance Σ g , so that the data are characterized by a finite mixture of normal distributions. The density of each observation where φ(·|·, ·) is a multivariate normal density. In practice, the number of clusters, G, is usually unknown and needs to be determined as part of the model inference.
The assumption of multivariate normal distributed clusters implies that the clusters are elliptical in shape. Banfield and Raftery (1993) proposed that constraints be placed on the covariance matrices to gain parsimony. These are specified using a modified eigenvalue decomposition of Σ g , namely where λ g is a constant which controls the cluster volume, D g is an orthogonal matrix of eigenvectors which control the orientation of the clusters, and A g a diagonal matrix, with entries proportional to the eigenvalues, which control the shape of the cluster.
We can restrict each property of the covariance Σ g (volume, shape, orientation) in different ways, resulting in fourteen different possible models (Biernacki et al., 2006). Throughout this paper, we will consider the ten covariance structures implemented in the mclust software (Fraley et al., 2012), as displayed in Figure 1. Each letter in the name of a model corresponds to the constraint placed on the volume, shape and orientation respectively. The constraint can be equal (E), variable (V) or identity (I). For example, in the EEV model, each cluster has the same volume and the same shape but the orientations of the clusters can differ. Other parametriza- tions of covariance matrices are useful in the context of model-based cluster analysis (e.g. Murphy, 2008, 2010;Biernacki and Lourme, 2014), but we do not consider them further here.
Based on choosing a covariance constraint and the number of groups G, we can fit the finite mixture of normal distributions using the EM algorithm (Dempster et al., 1977) to yield the maximum likelihood estimates of the model parameters and an estimate of the posterior probability that each observation x i belongs to group g. The resulting N × G matrix of probabilities, denoted by Z, has typical entry z ig , which is the estimated conditional posterior probability that observation i belongs to group g. This is where {(τ g ,μ g ,Σ g ) : g = 1, 2 . . . , G} are the maximum likelihood estimates of the model parameters. The probability matrix, Z, can be converted into group labels for the observations and thus the observations can be clustered. In a clustering problem where several possible values for G are considered (e.g. G = 1, 2, . . . , 9) and when all ten model types are also considered, we have to choose between about 83 different models. However, we would expect only a small number of these models to be strongly supported by the data. A frequently-used model selection approach is to choose the "best" model using the Bayesian Information Criterion (BIC), (Schwarz, 1978) where the BIC for model M m is given by BIC m = 2 log(L) − κ m log(N ). (2) In Equation 2, L is the maximized likelihood of the data, κ m is the number of estimated model parameters for model M m , and N is the number of observations. This approach was proposed for clustering by Dasgupta and Raftery (1998), and has been found to perform well (e.g. Steele and Raftery, 2010).
Once the model has been chosen, the cluster membership matrix is based on the parameters for that model alone. The other competing models are then discarded and model uncertainty is ignored.
3 Model Uncertainty and Bayesian model averaging Basing inferences on a single "best" model ignores uncertainty about what the best model is. This can result in underestimating uncertainty about quantities of interest such as cluster membership or the estimated density.
We address this using Bayesian model averaging (BMA) (Leamer, 1978;Madigan and Raftery, 1994;Hoeting et al., 1999), which proceeds as follows. If {M 1 , ..., M M } denotes the set of all models being considered and if ∆ is the quantity of interest, then the posterior distribution of ∆ given the data is All probabilities are conditional on the set of all models being considered. Madigan and Raftery (1994) argued that averaging over all of the models, as in Equation 3, provides better predictive ability than using any single model. One difficulty is that the integral in Equation 4 is intractable in most cases and so an approximation to the posterior model probability is required. We use the Bayesian Information Criterion (BIC) to derive approximate posterior model probabilities for Bayesian model averaging for model-based clustering (Dasgupta and Raftery, 1998). If all the models are equally likely a priori, so that P(M 1 ) = · · · = P(M M ) = 1/M , this yields  Strehl and Ghosh (2003) introduced the concept of a binary similarity matrix in the context of clustering. This matrix is N × N and for any pair of observations (i, j), the (i, j)th entry in the matrix is 1 if the ith and jth observations are in the same cluster in a given model, and 0 if they are not. They used this structure to form what they call pairwise cluster ensembles; Monti et al. (2003) and Kuncheva and Hadjitodorov (2004) also exploited this idea in a clustering context. In these articles, weights were given to various clustering methods and the similarity matrices were combined using these weights. The resulting matrix was described by Kuncheva and Hadjitodorov (2004) as a consensus matrix. Fern and Brodley (2003) extended the binary similarity matrix to soft-clustering techniques, which return a probability vector P(g|i, M m ), g = 1, . . . , G for an observation i, representing the probability that i belongs to each cluster under model M m . The values P(g|i, M m ) are analogous to the z ig values, as defined in Equation 1. Fern and Brodley (2003) used model-based clustering to cluster high-dimensional data by randomly projecting the data into a low-dimensional space and clustering the data in the low-dimensional space. The data are projected multiple times and the consensus matrices are averaged across all projections. Fern and Brodley (2003) defined S m ij as the probability that observations i and j belong to the same cluster under model M m . This can be calculated as

Similarity matrix
because z m ig is an estimate of the probability that observation i belongs to cluster g in model M m .
We can construct the similarity matrix S m for each model M m as follows: We propose using the matrix S m of probabilities that any pair of observations belong to the same cluster, when implementing BMA for model-based clustering. This ensures that we are averaging a quantity that has the same meaning across differing number of clusters and is invariant to cluster labeling.

Properties of the similarity matrix
The matrix S m for any model M m will be N × N where N is the number of observations. It is invariant to label switching, and its dimension is invariant to the number of clusters in the model. It can therefore be used to combine models with different numbers of clusters. It can be viewed as a similarity matrix between the data points (x 1 , . . . , x N ). The element s k ij of the matrix S m represents the probability that i and j belong to the same cluster. We can also define a quantity d ij = 1 − s ij as the probability that they are not in the same cluster. So, S = 1 − D, where D can be thought of as a dissimilarity matrix. Therefore, D can be used with any clustering algorithm which operates directly on a dissimilarity matrix.
We define s A as the minimum probability that two observations i and j in a set A belong together, which we can think of as the minimum probability that two elements in A belong to the same cluster. If we define d A as 1 − s A , we obtain the following result: Thus the maximum probability that two elements of A do not appear in the same cluster is d A . It follows that D can be used in hierarchical clustering with complete linkage (Sokal and Sneath, 1963) and the results will have an intuitive interpretation. In complete linkage, the dissimilarity between two groups G and H is the largest dissimilarity between opposite groups, that is, the maximum distance between an element of G and an element of H, namely Hierarchical complete-linkage clustering will then merge the groups with the smallest d complete . This continues until all the groups are merged.
The results can be shown on a dendrogram which can be cut at any level. There is an intuitive interpretation of the level of cut. If we cut the dendrogram at a particular probability level, any observations clustered together at that level have a probability of at least that value of all being in the same cluster.
We illustrate this with a toy example. Suppose we have six data points A, . . . , F and two equally likely clustering results. The first clustering attempt puts A, B, C in one cluster and D, E, F in the other, while the second attempt puts A, C, E in one cluster and B, D, F in the other. For simplicity we will use hard clustering where the probability of cluster membership is either 0 or 1 for each observation. If both models have equal probability, the S matrix will be If we proceed to implement hierarchical clustering on the corresponding dissimilarity matrix, with complete linkage, we get the dendrogram shown in Figure 2.
Averaging the two models and cutting the resulting dendogram at any probability level above 0.5 leads to a four-cluster solution. This makes sense. Both clusterings agree that A and C belong in the same cluster as do D and F . However, they disagree on B and E. This dendrogram gives a probabilistic representation of this situation. If we cut the dendogram in Figure 2 at 0.2, we see that A, B and C have at least a probability 0.2 of being in the same group, but cutting at 0.8 shows that A and C have at least a probability of 0.8 of being in the same group.

Summary of method
We can now carry out Bayesian model averaging by using the posterior model probabilities estimated from Equation 5 as weights and the similarity matrices for each . For each pair of observations i and j in our dataset, we propose assigning the probability vector
When using the default settings in mclust, a total of 83 candidate models were fitted. These include three one-cluster models (G = 1) and the ten possible covariance structures (Figure 1) for each number of clusters G = 2, 3, . . . , 9. One can of course fit a larger or smaller set of models if desired.

Clustering Fisher's iris data
The iris data (Fisher, 1936) gives the measurements in centimetres of the variables sepal length and width and petal length and width for 150 observations with 50 of each of three species of iris: versicolor, setosa and virginica. Table 1 shows the model-based clustering models with the highest BIC values for the iris data. These results show that the VEV model with two clusters has the highest BIC, but that there is considerable uncertainty about whether this is the best model. In a single model scenario, the VEV model with two clusters would be selected and clustering would be based solely on this model. However, it can be seen from Table 1, that the BIC value for the three-cluster VEV solution is similar to that for the selected model. The posterior model probabilities in the table are estimated using Equation 5 and we see that the posterior probability for this second best model is almost 40%.
In order to visualize the observations that are likely to be in the same cluster, we use a heatmap to represent the similarity matrix. The ordering of the observations in the heatmap is important. Here we present the data in species order, which yields an intuitive heatmap, as shown in Figure 3. However, in many applications this will not be the case, but the structure may be easier to see if the data are ordered using some seriation method. Figure 3(a) shows the heatmap for the similarity matrix for the two-cluster VEV model, while Figure 3(b) shows the three-cluster model; these are the two models with the highest BIC values. The red sections of the figure represent clusters of (i, j) pairs that have a high probability of being in the same cluster, which the white areas show (i, j) pairs that are unlikely to be in the same cluster. The clusters are isolated along the diagonal of the heatmap. The model that would be selected by BIC separates the setosa species very well from the other two species but it merges the two other species into one cluster. It can be seen also that the probabilities assigned to the co-clustering pairs are very high. The fact that there is large uncertainty associated with the two-cluster solution is not clear in the single model results.
The second most likely model shows separation between the three clusters with high probability, as we see in Figure 3(b). Thus the argument can be made that this model should contribute to the final clustering result. Figure 3(c) shows the heatmap for the similarity matrix produced from the Bayesian model averaging process. The method separates the large cluster into two clusters which approximate the known species groups quite well, but also reflect the uncertainty about whether there are really three species.
We can also present these results using dendrograms, as detailed in Section 4.2, and these are shown in Figure 4. The dendrograms show that the single model results cluster the observations into two clear groups, whereas the BMA results show the possibility of both a two-cluster and/or three-cluster solution, depending on the cutoff used.

Clustering Italian wines data
We now apply the methodology to a well-known dataset consisting of 27 chemical measurements on each of 178 wine samples belonging to three cultivars of wine (Barolo, Grignolino and Barbera) produced in the Piedmont region of Italy (Forina et al., 1986). Table 2 shows the number of samples collected in each vintage year for the three cultivars.   Table 3 shows the candidate models with the highest BIC and hence the highest approximate posterior probability; all other models had negligible posterior probability. Although the data consist of samples from three different cultivars, a sevencluster model was preferred to any of the three-cluster models. The seven clusters successfully separate the three cultivars, and also partly reflect the different vintage years shown in Table 2 (McNicholas and Murphy, 2008).
The heatmap of the similarity matrix of the optimal model is presented in Figure 5(a). Here the observations are shown in order of vintage year within cultivar. It appears that the clustering results partly reflect the vintage years as well as the cultivars.
We seriated the data to reorder the observations, using the order of the leaf nodes in a dendrogram produced by hierarchical clustering with complete linkage. In hierarchical displays, a decision is needed at each merge to specify which subtree should go to the left and which to the right (Hurley, 2012). We used the order  suggested in Gruvaeus and Wainer (1972), which ensures that objects at the boundaries of each class were located next to objects outside the class which they most resembled (Gordon, 1987). At a merge of clusters A and B, the new cluster is one of (A, B), (A , B), (A, B ), (A , B ), where A denotes A in reverse order. The new cluster is chosen to minimize the distance between the object in A placed adjacent to an object from B. The reordered similarity matrix, shown in Figure 5(b), has seven clear clusters, with uncertainty about the group membership of some of the observations. From Table 3 we see that the combined posterior probability of the two threecluster models (EVI and VVI) is high, at just under 40%. We use BMA to average the similarity matrices across the candidate models and show the results in Figure 6. Figure 6(a) shows the clustering in the originally presented order, after BMA has been carried out. It gives a clearer picture of the cultivar clustering, although there is some uncertainty, shown by the yellow colors. Figure 6(b) shows the seriated version, where the smaller clusters group naturally into three larger clusters, but with one small cluster that could belong to either of two cultivars.   The two three-cluster solutions give similar clustering results. In the sevencluster solution, three of the clusters are for the Grignolino observations, and two each for the Barbera and Barolo observations. The Barolo and Barbera clusters break down by vintage year, whereas for the Grignolino observations the clusters correspond less clearly to vintage years.

Bayesian model averaging for density estimation.
The fitted mixture model produced by model-based clustering provides an overall density estimate for the data generation process. Ferguson (1983) showed that finite mixtures of normal distributions can approximate any distribution on the real line to within any given accuracy. Thus the density estimate produced by model-based clustering can be used as a density estimation method. The performance of densities fitted in this way was assessed by Roeder and Wasserman (1997) for the univariate case and by Fraley and Raftery (2002) for the bivariate case. Both of these studies showed model-based clustering to be competitive with state of the art kernel density estimation methods.
We propose using Bayesian model averaging of the density estimates for each model using the posterior model probabilities from Equation 5. So, given modelbased density estimatesf Mm (x), for m = 1, 2, . . . , M , we havê We compare the BMA results with kernel density estimation methods. Kernel density estimation has long been used for density estimation and visualization of univariate data (Rosenblatt, 1956;Parzen, 1962;Jones et al., 1996). The kernel density estimatef h of a univariate density f based on a random sample X 1 , . . . , X N of size N isf where K h (·) = (1/h)K(·/h) for a kernel function K, assumed to be a symmetric probability density and a bandwidth h (the smoothing parameter). We also comparê f BMA with the single-model estimatef SM , wherê where M m is the model with the highest value of BIC. (SM stands for "single model".) Wand and Jones (1994) described the extension of the method to the bivariate case. The multivariate kernel density estimate iŝ where H is now a d × d positive definite matrix and K is a d-variate spherically symmetric density function; a common simplification is to use a diagonal H (e.g., Wand and Jones, 1994, Chapter 4.2). The performance of density estimation procedures can be compared using mean integrated squared error (MISE) or expected Kullback-Leibler divergence (KL), defined as follows: where the expected value is taken with respect to the resulting density estimate,f , computed from a random sample, X 1 , X 2 , . . . , X n , drawn from f . For both criteria, smaller is better. MISE (Equation 7) is the most commonly used measure of performance, while Kullback-Leibler divergence (Equation 8) provides an alternative which considers the differences in densities on a logarithmic scale; this places more emphasis on differences in regions of low density. The asymptotic performance of kernel density estimation under the MISE criterion was described by Silverman (1986) and Scott (1992) and the form of the optimal bandwidth for the kernel was derived. Hall (1987) studied the asymptotic performance of kernel density estimation under Kullback-Leibler divergence and showed that the optimal bandwidth in this case leads to a smoother density estimate than under MISE. The optimal bandwidths differ because Kullback-Leibler puts a larger penalty on regions where a density estimate has low density but the true density has high density. Wand and Jones (1993) and Duong and Hazelton (2003) investigated the use of unconstrained parametrisations of the H matrix as opposed to diagonal ones on simulated target densities, and concluded that this can improve efficiency.
We use the extended version of multivariate kernel density estimation as described in Duong and Hazelton (2003) and Duong (2007) as a benchmark for comparison with the model-based clustering density estimates. This approach uses a two-stage estimation procedure, where a pilot density estimate is used to get an improved estimate of the optimal bandwidth compared to the standard plug-in estimates (e.g. Silverman, 1986, Section 3.4.2). Duong (2007) used kernels with non-diagonal H matrices.
6.1 Density estimation results 6.1.1 Simulation study for bivariate density estimation Fraley and Raftery (2002) defined bivariate analogs of the first ten of the univariate datasets defined in Marron and Wand (1992) to compare the performance of the fitted model-based clustering density with multivariate kernel density estimation. Contour plots of the densities are shown in Appendix A as well as the parameter settings for the actual densities we used.
We produced 250 simulations of 250 observations from each of the same ten distributions and compared density estimation using model-based clustering with the Bayesian model averaged result in terms of the mean integrated squared error and the Kullback-Leibler distance. We compared the performance of the model-based clustering approaches to the extended kernel density estimaton method described in Duong and Hazelton (2003) as implemented in the ks R package (Duong, 2007(Duong, , 2014. The results of the simulation study are shown in Table 4 where the numbers in the MISE columns are the MISE for kernel density estimation and for the single model respectively divided by the MISE for the BMA method. Similarly, for the Kullback-Leibler column, we divide the Kullback-Leibler distance for kernel density estimation and for the single model by the Kullback-Leibler distance for the BMA method. Density estimation using model-based clustering with BMA compares very favorably with kernel density estimation according to the KL criterion, while it compares favorably in the majority of cases by the MISE criterion. It also compares favorably with density estimation using model-based clustering with a single model in all cases and by both criteria, although the gains are more modest. Fraley and Raftery (2002) gave further comparisons of single model density estimation with other kernel density estimation methods, namely Gaussian kernel density estimation using both the normal optimal bandwidth and cross-validated bandwidth (Bowman and Azzalini, 1997).

Simulation study for higher dimensional density estimation
We also investigated density estimation for higher dimensional data. We first added dimensions consisting of observations from the standard normal distribution to the existing distributions used in Section 6.1.1. The first three rows in Table 5 refer to bivariate distributions with one extra such dimension added, while the next three rows refer to the same distributions with four extra standard normal dimensions. We also implemented three-and six-dimensional versions of the bimodal distribution. We allowed for variable separation of the modes as in the bimodal case, and the results are in rows seven to twelve of Table 5. The settings used for all these simulations are in Appendix A.2.
Table 5 MISE and Kullback-Leibler (KL) distance ratios for density estimation via model-based clustering with Bayesian model averaging (BMA) as against kernel density estimation with the ks package (KS) and single model model-based clustering (SM). The top six lines comprise simulations of bivariate data as before with either one or four extra dimensions of standard normal. The bottom six lines comprise extensions to three and six dimensions of the bimodal density of Fraley and Raftery (2002)  Model-based clustering with BMA strongly outperformed kernel density estimation in all cases except one for both three-dimensional and six-dimensional data, according to both criteria, MISE and KL. Model-based clustering with BMA also uniformly outperformed single-model model-based clustering, although the gain was more modest. Density estimation using model-based clustering can be carried out for higher dimensions, and performs well. However, we limited ourselves to six dimensions because the ks R package provides results for at most six dimensions. We conjecture that model-based clustering's gain in performance would be even greater for higher dimensions.

Lansing maples
We now consider the Lansing trees data set (Gerrard, 1969) where we are interested in the density of the maple trees in a forest. The models with the highest BIC are shown in Table 6 and we can see that three models have non-negligible posterior probability. A six-cluster solution has approximately 16% posterior model probability with two seven-cluster solutions having approximately 83% between them. There are very small probabilities associated with five and eight cluster solutions respectively. However, the best model that would be chosen according to BIC is VII with seven clusters.
We can compare graphically the contour plots for the models with the highest BIC and the plot for the density estmation after BIC (Figure 8 Figure 9 shows the differences between the density estimates produced by BMA and the density estimate for the VII model with seven clusters. Again, the difference in probability density around the cluster at (0.4, 0.1) is evident. There are some other small areas of smoothing denoted in blue. The brown peaks are higher densities to compensate for the lower density at (0.4, 0.1).

Discussion
We have implemented Bayesian model averaging for model-based clustering and density estimation. The results show that the BMA assesses the model uncertainty better than the single best model approach.
We have proposed similarity matrices as a way of representing ensembles of clustering solutions. They provide an intuitive way of visualising clustering solutions and are consistent across models with different numbers of clusters. They are also invariant to label switching between models. Similar ideas were put forward by Strehl and Ghosh (2003) and Monti et al. (2003) using binary classification matrices, while Fern and Brodley (2003) extended the idea to soft classifications when proposing their random projection clustering method. All of these papers use different ways of combining the matrices than BMA. Kuncheva and Hadjitodorov (2004) suggested that "cutting" the matrix at a certain threshold is equivalent to running the single link algorithm and cutting the dendrogram at that level. We argue that using complete linkage gives a more intuitive interpretation.
We have shown how to apply Bayesian model averaging to clustering solutions using the similarity matrix. This provides a statistical postprocessing method which takes account of model uncertainty. When used on datasets with well-documented structure or a known ground truth, results achieved with Bayesian model averaging are consistent with the ground truth. We have shown that carrying out BMA on datasets that have no ground truth might give additional information than using model-based clustering alone and at little computational cost.
The model-based density estimation method described in Fraley and Raftery (2002) gives better results in many cases than those given by well-known kernel density estimation methods for the simulated datasets analysed here. With the addition of the Bayesian model averaging described in this paper, even more improvement can be seen, and this becomes more pronounced in higher dimensions.
In the case where the derivative of the density estimate is of interest (e.g. Jones (1994) and others), a separate bandwidth is needed for each order derivative (Chacón and Duong, 2013). An advantage of the finite mixture density estimate and the BMA mixture density estimate is that a single estimate is produced which can be used for estimating derivatives of any order.
The proposed methodology could also be used when more than one family of component distributions is under consideration. For example, the normal mixture models with eigendecomposed covariances used in mclust (Fraley and Raftery, 2002), the normal mixtures with factor analytic covariance structure used in pgmm (McNi-cholas and Murphy, 2008), the multivariate t mixtures used in EMMIX (McLachlan and Peel, 1999) and the skew-t mixtures used in EMMIXuskew (Lee and McLachlan, 2013a) could be combined using the BMA procedure.
One previous approach to Bayesian model averaging within model-based clustering (Wei and McNicholas, 2014) considers noninvariance of model-based clustering results to cluster labeling by matching the clusterings for competing models using an cluster agreement criterion. Further, they combine the clusters in the larger model so that the models being combined have the same number of clusters, G. We have proposed a very different approach for Bayesian model averaging for model-based clustering by choosing a quantity of interest, ∆, that is invariant to cluster labeling and that has a common meaning for all values of G.
Overall, our results suggest that Bayesian model averaging is a useful postprocessing tool for model-based clustering and density estimation. It often helps and seldom disimproves the results, so it could be used routinely as part of model-based clustering.    Figure 12 Strongly skewed distribution: Contour map denoting the density of a strongly skewed distribution.        Figure 19 Claw distribution: Contour map denoting the density of the claw distribution. For these results we merely added standard normal observations to the simulated bivariate clusters. We did not differentiate between the clusters as cluster membership is not important in this setting.
A.2.2 Three-and six-dimensional extensions Table 17 3D Bimodal Data: Settings used to simulate the density of a 3D bimodal distribution.Note that displacement means the distance from the origin -in this case along the line y = x. Separation as described in the results tables is the separation of the group centres.