Inference in model-based cluster analysis.
A new approach to cluster analysis has been introduced based on parsimonious geometric modelling of the within-group covariance matrices in a mixture of multivariate normal distributions, using hierarchical agglomeration and iterative relocation. It works well and is widely used via the MCLUST software available in S-PLUS and StatLib. However, it has several limitations: there is no assessment of the uncertainty about the classification, the partition can be suboptimal, parameter estimates are biased, the shape matrix has to be specified by the user, prior group probabilities are assumed to be equal, the method for choosing the number of groups is based on a crude approximation, and no formal way of choosing between the various possible models is included. Here, we propose a new approach which overcomes all these difficulties. It consists of exact Bayesian inference via Gibbs sampling, and the calculation of Bayes factors (for choosing the model and the number of groups) from the output using the LaplaceMetropolis estimator. It works well in several real and simulated examples.
Introduction
Banfield and Raftery (1993) - hereafter BR - building on work of Murtagh and Raftery
(1984), introduced a new approach to cluster analysis based on a mixture of multivariate
normal distributions, where the covariance matrices are modeled parsimoniously in a geometrically
interpretable way. The general finite normal mixture distribution for n data
points dimensions with K groups is
Y
where OE(\Deltaj-; \Sigma) is the multivariate normal density with mean - and covariance matrix \Sigma,
a vector of group mixing proportions such that - k - 0 and
The BR approach is based on a variant of the standard spectral decomposition of \Sigma k ,
namely
where - k is a scalar, A
an orthogonal matrix for each Each factor in equation (2) has a geometric
controls the volume of the kth group, A k its shape and D k its orientation.
By imposing constraints such as - (i.e. each group has the
same volume and shape, but they have different orientations), one obtains different models,
which lead in turn to different clustering algorithms. The models considered here, including
parsimonious spherically shaped ones, are listed in Table 1.
BR developed algorithms aimed at maximizing the classification likelihood
Y
as a function of both ' and -, where - is the vector of group memberships, namely -
belongs to the kth group. These algorithms use hierarchical agglomeration and iterative
relocation. They worked well on several real and simulated data sets and are now fairly
widely used. They are implemented in the software MCLUST, which is both an S-PLUS
function and a Fortran program available from StatLib. 1
However, the BR algorithms have limitations, several of which are common to all agglomerative
hierarchical clustering methods:
To obtain the Fortran version, send the e-mail message "send mclust from general" to the address
statlib@stat.cmu.edu.

Table

1: Clustering models. The entries indicate whether the feature of interest (shape,
orientation or volume) is the same for each group or not.
Model \Sigma k Shape Orientation Volume
1. -I Spherical None Same
2. - k I Spherical None Different
3. \Sigma Same Same Same
4. - k \Sigma Same Same Different
5. -D k AD t
k Same Different Same
k Same Different Different
7. - k DA k D t Different Same Different
8. \Sigma k Different Different Different
(a) They give only point classifications of each individual and produce no assessment of
the associated uncertainty.
(b) They tend to yield partitions that are suboptimal (even if often good). This is due to
the use of hierarchical agglomeration.
(c) Estimates of the model parameters ' based on the estimated partition tend to be biased
(Marriott, 1975).
(d) They assume the mixing proportions - k in equation (1) to be equal.
(e) The algorithms based on models 5 and 6 in Table 1 require the shape matrix A to
be specified in advance by the user. This is a drawback in general, although it can
sometimes be useful.
(f) To choose K, the number of groups, BR proposed an approximation to the posterior
probabilities based on a quantity called the AWE (Approximate Weight of Evidence).
While this has worked fairly well in practice, it is quite crude.
(g) BR proposed no formal way of choosing among the possible models; this must be done
by the user.
We know of no way of fully overcoming limitations (a), (f) and (g) other than the fully
Bayesian analysis that we develop here. Other possible ways of overcoming (b)-(e) are
discussed in Section 4.
Here we present a new approach to clustering based on the models in Table 1; it consists
of fully Bayesian inference for these models via Gibbs sampling. This overcomes all the
limitations mentioned. A recently proposed way of calculating Bayes factors from posterior
simulation output, the Laplace-Metropolis estimator (Raftery 1995; Lewis and Raftery 1994),
is used to choose the model and determine the number of groups in one step. In Section 2.1
we describe Bayesian estimation for the models of Table 1 using Markov Chain Monte Carlo
methods. In Section 2.2 we outline how Bayes factors can be calculated from the
MCMC output and used to determine the appropriate model and the number of groups. In
Section 3 we show the methods at work on real and simulated data sets.
Bayesian Inference for the Banfield-Raftery Clustering
Models Using the Gibbs Sampler
2.1 Estimation
We assume that the data to be classified, x , arise from a random
vector with density (1) and that the corresponding classification variables - i are unobserved.
We are concerned with Bayesian inference about the model parameters ', - and the classification
indicators - i . MCMC methods provide a general recipe for the Bayesian analysis of
mixtures. For instance, Lavine and West (1992) and Soubiran, Celeux, Diebolt, and Robert
(1991) have used the Gibbs Sampler for estimating the parameters of a multivariate Gaussian
mixture without assuming any specific characteristics for the component variance matrices.
Diebolt and Robert (1994) have considered the Gibbs sampler and the Data Augmentation
method of Tanner and Wong (1987) for general univariate Gaussian mixtures and proved
that both algorithms converge in distribution to the true posterior distribution of the mixture
parameters.
Like these authors, we use conjugate priors for the parameters - and ' of the mixture
model. The prior distribution of the mixing proportions is a Dirichlet distribution
and the prior distributions of the means - k of the mixture components
conditionally on the variance matrices \Sigma k are Gaussian: - k j\Sigma k - N (-
prior distribution of the variance matrices depends on the model, and will be given for
each model in turn.
We estimate the eight models in Table 1 by simulating from the joint posterior distribution
of -, ' and - using the Gibbs sampler (Smith and Roberts 1993). In our case, this consists
of the following steps:
1 Simulate the classification variables - i according to their posterior probabilities
namely
2 Simulate the vector - of mixing proportions according to its posterior distribution
conditional on the - i 's.
3 Simulate the parameters ' of the model according to their posterior distributions conditional
on the - i 's.
The validity of this procedure, namely the fact that the Markov chain associated with the
algorithm converges in distribution to the true posterior distribution of ', was shown by
Diebolt and Robert (1990) in the context of one-dimensional normal mixtures. Their proof
is based on a duality principle, which uses the finite space nature of the chain associated
with the - i 's. This chain is ergodic with state space Kg and is thus geometrically
convergent and even '-mixing. These properties transfer automatically to the sequence of
values of ' and -, and important properties such the Central Limit Theorem or the Law of
the Iterated Logarithm are then satisfied (Diebolt and Robert 1994; Robert 1993).
The same results also apply here, the only difference being the more complex simulation
structure imposed by the variance assumptions. Steps 1 and 2 do not depend on the considered
model. Step 1 is straightforward, and Step 2 consists of simulating - from its conditional
posterior distribution, namely - D (ff 1
3 is not the same for the different models of Table 1, and is described in the Appendix for
each model in turn.
2.2 Choosing the Number of Groups and the Model by Bayes
Factors
BR left the choice of model to the user, while noting that model 6 of Table 1 (corresponding
to their S   criterion) gives good results in many situations. They based the choice of number
of clusters on the AWE criterion, which is a crude approximation to twice the log Bayes factor
for that number of clusters versus just one cluster.
Here we develop a way of choosing both the model and the number of groups at the
same time, using a more accurate approximation to the Bayes factor than that of BR.
We compute approximate Bayes factors from the Gibbs sampler output using the Laplace-
Metropolis estimator of Raftery (1995). This was shown to give accurate results by Lewis
and Raftery (1994).
In what follows, the word "model" refers to a combination of one of the models in Table
1 with a specified number of clusters. The Bayes against another
model M 0 given data D is the ratio of posterior to prior odds, namely
In equation (4),
Z
is the vector of parameters of M k , and pr(' k jM k ) is its prior density
this is called the integrated likelihood of model M k . For a review of Bayes factors, their
calculation and interpretation, see Kass and Raftery (1995). Bayesian model selection is
based on Bayes factors, whose key ingredient is the integrated likelihood of a model. Our
main computational challenge is thus to approximate the integrated likelihood using the
Gibbs sampler output.
We do this using the Laplace-Metropolis estimator of the integrated likelihood (Raftery
1995; Lewis and Raftery 1994). The Laplace method for integrals is based on a Taylor series
expansion of the real-valued function f(u) of the d-dimensional vector u, and yields the
approximation Z
where u   is the value of u at which f attains its maximum, and A is minus the inverse
Hessian of f evaluated at u   . When applied to equation (5) it yields
where d is the dimension of ', ~
' is the posterior mode of ', and \Psi is minus the inverse
Hessian of logfpr(Dj')pr(')g, evaluated at
'. Arguments similar to those in the
Appendix of Tierney and Kadane (1986) show that in regular statistical models the relative
error in equation (7), and hence in the resulting approximation to
is sample size.
While the Laplace method is often very accurate, it is not directly applicable here because
the derivatives it requires are not easily available. The idea of the Laplace-Metropolis estimator
is to get around the limitations of the Laplace method by using posterior simulation
to estimate the quantities it needs. The Laplace method requires the posterior mode, (~-; ~ '),
and j\Psij. The Laplace-Metropolis estimator estimates these from the Gibbs sampler output
using robust location and scale estimators. The likelihood at the approximate posterior
mode is
pr(Dj ~
Y
These quantities are then substituted into equation (7) to obtain the integrated likelihood,
and Bayes factors are computed by taking ratios of integrated likelihoods, as in equation (4).
3 Examples
We now present three examples to illustrate the ability of our methods to overcome the
limitations (a)-(g) of other methods described in Section 1. The first example uses simulated
data, while the second and third examples are based on real data sets.
For each example, we consider only the models [-I]; [- k I]; [\Sigma]; and [- k \Sigma]. Models [-I] and
are probably the most used Gaussian mixture models for clustering data (e.g. McLachlan
and Basford 1988), and the generalizations of these, [- k I] ands [- k \Sigma], to allow for different
volumes have proved to be powerful in several practical situations (Celeux and Govaert
1994).
Our priors are chosen from among the conjugate priors of Section 2.1 so as to be fairly
flat in the region where the likelihood is substantial and not much greater elsewhere. Thus
they satisfy the "Principle of Stable Estimation" (Edwards, Lindman and Savage 1963), and
so it could be expected that the results would be relatively insensitive to reasonable changes
in the prior; we also checked this empirically for each example.
We used -
x and S are the empirical mean vector and variance matrix of the whole
data set, and -
oe 2 is the greatest eigenvalue of S. (The other notation used is defined in the


Appendix

). The amount of information contained in this prior is similar to that contained
in a typical single observation. Thus the prior may be viewed as comparable to the true
prior of a person with some, but rather little, information. Similar priors have been used
for generalized linear models by Raftery (1993) and for linear regression models by Raftery,
Madigan and Hoeting (1994). In each example we assessed the sensivity of the results to
changes in this prior and found it to be small; some of the sensitivity results for the first
example are included below.
Gauss

Figure

Example 1: Simulated data.
3.1 Example 1: Simulated Data
We simulated 200 points from a bivariate two-component Gaussian mixture with equal pro-
portions, mean vectors - t
the
data are shown in Figure 1. The first 600 iterations from the Gibbs sampler output for the
model [- k I] with two groups are shown in Figure 2. Convergence was almost immediate and
successive draws were almost independent; similar results were obtained from other starting
values. We used 1,500 iterations, estimated by the gibbsit program to be enough to estimate
the cumulative distribution function at the .025 and .975 quantiles to within \Sigma:01 for
all the parameters (Raftery and Lewis 1993, 1995); the first 10 iterations were discarded.
The model comparison results are shown in Table 2. The correct model, [- k I], and the
correct number of groups, 2, are strongly favored. The posterior means of the parameters for
the preferred model are - are close to
the true values. The marginal posterior distribution is summarized in Figure 3, which shows
the posterior distribution of the principal circles of the two groups.
Sensitivity to the prior distribution is investigated in Table 3. A new Gibbs sampler

Figure

2: Example 1: Time series plot of the first 600 Gibbs sampler iterations: (a) volume
parameters; (b) mean for group 1; (c) mean for group 2.

Table

2: Example 1: Approximate Log Integrated Likelihoods.
No. of groups [-I] [- k I] [\Sigma] [- k \Sigma]

Figure

3: Example 1: Posterior Distribution of Principal Circles for the [- k I] model with
two groups. There is one circle for each Gibbs sampler iteration and each group, with center

Table

3: Example 1: Sensitivity of selected results to changes in the prior hyperparameters
for 1500 simulations for the selected model.
Prior parameters log B 23 Pr[-
is the log Bayes factor for two groups against three groups. D denotes the data. -
x
is the overall mean of the data -
denotes the means of the two optimal partitions,
namely -
run was done for each choice of prior parameters, and so the differences in Table 3 are due
to both true sensitivity and Monte Carlo variation; the true sensitivity is thus likely to
be smaller. The estimation results are quite insensitive. The testing results are somewhat
more sensitive, which is to be expected (Kass and Raftery 1995), but the overall conclusions
remain the same over all combinations of prior parameters considered.
Perhaps the greatest advantage of the present approach is that it fully assesses uncertainty
about group membership, rather than merely giving a single "best" partition.
In

Figure

4, this is summarized by showing the uncertainty for each point, measured by
. When it is clear that x i belongs to the k-th group,
then
is small, and so U i is also small. In these data, U i is large for
only one point, no. 55, the one that lies on the boundary of the two groups, for which
3.2 Example 2: Butterfly Classification

Figure

3.2 shows four wing measurements of a butterfly. Here we analyze data on two of
these measurements, z 3 and z 4 , for 23 butterflies, shown in Figure 6, from Celeux and Robert
(1993). The aim is to decide how many species are represented in this group of insects, and
to classify them.

Table

4 shows that model [- k \Sigma] with four groups is favored quite strongly over the
alternatives. The posterior means of the parameters are: -
Figure

4: Example 1: Uncertainty plot. At each point a vertical line of length proportional
to U
is plotted. The longest line is of length 0.5.

Figure

5: Example 2: Butterfly Measurements.
z3
z4
20 22 24 26 282030246810121315171921
22 23

Figure

Example 2: Butterfly Data. Values of (z 3 ; z 4 ) for 23 butterflies.

Table

4: Example 2: Approximate Log Integrated Likelihoods.
No. of groups [-I] [- k I] [\Sigma] [- k \Sigma]
0:13. The most likely group memberships a posteriori are shown in Figure
7, along with the associated uncertainties.
All the butterflies are classified with confidence except numbers 4 and 15 which are close
to the boundary between groups 1 and 3. Group 4 consists of just one butterfly, which is
clearly out on its own. The correct classification is known, and there are indeed four groups.
In addition, the correct classification is exactly equal to the optimal classification found by
our methods (Celeux and Robert, 1993).
Figure

7: Example 2: Estimated Group Memberships and Uncertainty Plot for the Butterfly
3.3 Example 3: Kinematic Stellar Data
Until fairly recently, it was believed that the Galaxy consists of two stellar populations,
the disk and the halo. More recently, it has been hypothesized that there are in fact three
stellar populations, the old (or thin) disk, the thick disk, and the halo, distinguished by their
spatial distributions, their velocities, and their metallicities. These hypotheses have different
implications for theories of the formation of the Galaxy. Some of the evidence for deciding
whether there are two or three populations is shown in Figure 8, which shows radial and
rotational velocities for stars, from Soubiran (1993).

Table

5 shows that model [- k \Sigma] is preferred and that there is strong evidence for three
groups as against two. The balance of astronomical opinion has also tilted towards this
conclusion, but based on much more information than just the velocities used here, including
star positions and metallicities (Soubiran 1993). It is impressive that such a strong conclusion
can be reached with the present methods using only a relatively small part of the total
available information.
The posterior means of the parameters for the preferred model are: -
552. The corresponding partition is shown in Figure 9.
U
Figure

8: Example 3: Kinematic Stellar Data. Radial (U) and rotational (V ) velocities for
2,370 stars in the Galaxy. Source: Soubiran (1993).

Table

5: Example 3: Approximate Log Integrated Likelihoods.
No. of groups [-I] [- k I] [\Sigma] [- k \Sigma]
astr

Figure

9: Example 3: Optimal partition for the model [- k \Sigma].
Figure

10: Example 3: Uncertainty plot.
The uncertainty plot is shown in Figure 10. The areas of high uncertainty are those on
the boundaries between any two of the three groups. The greatest uncertainty is in the two
small areas where all three groups intersect.
We have presented a fully Bayesian analysis of the model-based clustering methodolology of
Banfield and Raftery (1993), which overcomes many of the limitations of that approach. It
appears to work well in several examples.
Alternative frequentist approaches, which might be easier to implement, consist of maximizing
the likelihood using the EM algorithm or of maximizing the classification likelihood
using the Classification EM (CEM) algorithm. Celeux and Govaert (1995) considered those
approaches to the full range of clustering models derived from the eigenvalue decomposition
of the group variance matrices, including those considered here. They have shown in particular
how it is possible to find the maximum likelihood estimate of the shape matrix A.
Both approaches could overcome limitations (b), (d) and (e) of section 1, and the maximum
likelihood approach could also overcome limitation (c). They would ovecome difficulty (a)
only partly: they do provide an estimate of the uncertainty about group membership, but
this assessment is incomplete because it does not take account of uncertainty about - and
'. They do not overcome limitations (f) and (g).
In our examples, we explicitly considered only models 1-4 of Table 1; these were sufficient
for the data we considered. However, more generally it would be useful to consider all eight
models, proceeding in the same way and using the results of Section 2.1.


Appendix

Gibbs Sampling for the Clustering Models
We now describe Step 3 of Gibbs sampling for each of the eight clustering models. Given a classification
vector use the notation
the component-wise statistics of location and scale In what follows, we denote
a model by the eigenvalue decomposition of its variance matrix written between brackets. For
instance, [-D k AD t
denotes the mixture model with equal volumes, equal shapes and different
orientations.
(a) Model [-I ]. Here the scale parameter - is common to all components of the mixture. We
assume that the prior distribution on the parameters is conjugate, namely that
2 ae) means that - has the inverted gamma distribution
\Gamma[ 1r]ae \Gammar=2
The posterior distribution on (- is therefore a convolution of normal distributions on
the - k 's and of an inverse gamma distribution on -.
The Gibbs components of Step 3 are then:
3.1 For
I p
with -
3.2 Simulate
(b) Model [- k I ]. When the variance scales are different, the prior distributions are similar for all
components:
We recover the case treated in Diebolt and Robert(1994), namely that in which the groups (-
are generated separately:
3.1 For
I p
with -
3.2 Simulate
(c) Model [\Sigma]. There is no need to consider the eigenvalue decomposition of the covariance matrix
\Sigma, and the prior distribution is given by
means that \Sigma has the inverse Wishart distribution
pr(\Sigma) / j\Sigmaj \Gamma(m+p+1)=2 expftr(\Psi\Sigma \Gamma1 )=2g:
Step 3 of Gibbs sampling is then decomposed as follows:
3.1 For
with -
3.2 Simulate
ae
(d) Model [- k \Sigma 0 ]. The prior distribution has three components:
We make the model identifiable by setting of Gibbs sampling is then decomposed
as follows:
3.1 For
with -
3.2 Simulate
3.3 Simulate
ae
Model [-D k AD t
Here, the distinction between - and A is entirely geometric (volume versus
shape), and we will therefore consider a single parameter A, where the first term of the diagonal
a 1 is no longer constrained to be equal to 1. The prior distribution of the parameters - k and A of
the model is then
The distribution of the orthogonal matrix D k is more delicate to specify. We use the marginal of
an inverse Wishart distribution W \Gamma1
analogy with the previous cases. The marginal
distribution can then be derived explicitly, because of the choice of I p as a "scale" matrix in the
Z
Y
(D k D t
Y
(D k D t
where (D k D t
denotes the t-th element of the diagonal matrix D k D t
k . The posterior distribution
of the whole set of parameters, fA; (D
Y
\Gamman k =2
Y
Y
a \Gammar t
Y
(D k D t
and leads to the following Gibbs steps:
3.1 For
3.2 Simulate
a t jD
3.3 We use the approximation to the distribution of the D k in (14), i.e. we assume that D k and
A are the direction and shape components of an inverse Wishart random variable and that they
are independent. This is true only asymptotically (i.e. as the number of degrees of freedom goes
to infinity) (Anderson 1984, Theorem 13.5.1) but it considerably simplifies the simulation, with
moderate effects (if any) on the resulting posterior distribution. If we thus assume that the couples
(D are all distributed as W \Gamma1
we derive from the posterior distribution
Y
\Gamma2 tr
A
\Theta
that
exp
A
This is a condensed way of saying that the diagonal elements of A are distributed according to the
inverted gamma distributions,
a t j- Ig@ n
A
The D k 's are then distributed a posteriori as the principal direction vectors from the following
inverse Wishart distribution,
(f) Model [- k D k AD t
k ]. The Gibbs simulations of the - k 's and D k 's are similar to those of Model
replacing A in (15).
The simulation of A is also quite close to the previous version since
exp
\Gammatr
A
If we assume - k - Ig(l k =2; ae k =2), the posterior distribution of - k in the Gibbs sampler is
(g) Model [- k DA k D t ]. There is no need to isolate the volume, - k , and we will consider parameters
A k , where the first term of diagonal a 1k is not constrained to be equal to 1. The prior distribution
on the parameters - k , A k and D of the model is then
a tk - Ig(r tk =2; ae tk =2)
We have
Z
Y
Y
The posterior distribution of the whole set of parameters, (D;
Y
\Gamman k =2
Y
Y
a \Gammar tk
tk
ae tk =2
Y
and leads to the followings Gibbs steps.
3.1 For
3.2 Simulate
a tk jD; - Ig
3.3 Simulate
and use a Metropolis step to recover
exp
\Gamma2 tr
A
Model [\Sigma k ]. This is the standard Gaussian mixture model considered by Lavine and West
(1992) and by Soubiran et al. (1991). In this case, there is no need to make use of the eigenvalue
decomposition of \Sigma k . The prior distribution on (- k ; \Sigma k ) is then
and the corresponding Gibbs step is, for



--R

An Introduction to Multivariate Statistical Analysis
"Model-based Gaussian and non Gaussian Clustering,"
"Gaussian parsimonious clustering models,"

" Bayesian estimation of finite mixture distributions,"
" Bayesian estimation of finite mixture distributions Part II: sampling implementation,"
"Bayesian statistical inference for psychological research,"
"Bayes factors,"
" A Bayesian method for classification and discrimination,"
"Estimating Bayes factors via posterior simulation with the Laplace-Metropolis estimator,"
"Separating mixtures of normal distributions,"
Mixture Models
"Fitting straight lines to point patterns,"
"Approximate Bayes factors and accounting for model uncertainty in generalized linear models,"
"How many iterations in the Gibbs sampler?"
"The number of iterations, convergence diagnostics, and generic Metropolis algorithms."

"Acounting for model uncertainty in linear regression,"
"Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods,"
"Kinematics of the Galaxy's stellar population from a proper motion survey,"
"Analyse de m'elanges gaussiens pour de petits 'echantillons : application `a la cin'ematique stellaire,"
--TR

--CTR
Mario A. T. Figueiredo , Anil K. Jain, Unsupervised Learning of Finite Mixture Models, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.3, p.381-396, March 2002
Zhihua Zhang , Kap Luk Chan , Yiming Wu , Chibiao Chen, Learning a multivariate Gaussian mixture model with the reversible jump MCMC algorithm, Statistics and Computing, v.14 n.4, p.343-355, October 2004
bayes factor;gaussian mixture;eigenvalue decomposition;gibbs sampler