# Bayesian Finite Mixture Models

Bettina Grün
WU Vienna University
of Economics and Business &Gertraud Malsiner-Walli
WU Vienna University
of Economics and Business<EMAIL_ADDRESS>gertraud.malsiner-
<EMAIL_ADDRESS>
###### doi:

10.1002/9781118445112.stat08373

Bettina Grün, Gertraud Malsiner-Walli Bettina Grün
Institute for Statistics and Mathematics
WU Vienna University of Economics and Business
Welthandelsplatz 1, 1020 Wien, Austria
E-mail:
Gertraud Malsiner-Walli
Institute for Statistics and Mathematics
WU Vienna University of Economics and Business
Welthandelsplatz 1, 1020 Wien, Austria
E-mail: Finite mixture models are a useful statistical model class for
clustering and density approximation. In the Bayesian framework finite mixture
models require the specification of suitable priors in addition to the data
model. These priors allow to avoid spurious results and provide a principled
way to define cluster shapes and a preference for specific cluster solutions.
A generic model estimation scheme for finite mixtures with a fixed number of
components is available using Markov chain Monte Carlo (MCMC) sampling with
data augmentation. The posterior allows to assess uncertainty in a
comprehensive way, but component-specific posterior inference requires
resolving the label switching issue. In this paper we focus on the application
of Bayesian finite mixture models for clustering. We start with discussing
suitable specification, estimation and inference of the model if the number of
components is assumed to be known. We then continue to explain suitable
strategies for fitting Bayesian finite mixture models when the number of
components is not known. In addition, all steps required to perform Bayesian
finite mixture modeling are illustrated on a data example where a finite
mixture model of multivariate Gaussian distributions is fitted. Suitable prior
specification, estimation using MCMC and posterior inference are discussed for
this example assuming the number of components to be known as well as unknown.
cluster analysis, finite mixture, label switching, Markov chain Monte Carlo,
mixture of finite mixtures model, number of components, prior specification,
telescoping sampler This is a preprint of an article published in _Wiley
StatsRef: Statistics Reference Online_. Please cite as: Bettina Grün and
Gertraud Malsiner-Walli (2022) Bayesian finite mixture models. In N.
Balakrishnan, Theodore Colton, Brian Everitt, Walter Piegorsch, Fabrizio
Ruggeri, and Jef L. Teugels, editors, _Wiley StatsRef: Statistics Reference
Online_ .

##  1 Introduction

Finite mixture models have a long history in statistics with their first use
proposed by Newcomb (1886) and Pearson (1894) more than hundred years ago.
They represent a useful model class for cluster analysis as well as semi-
parametric density estimation. In addition the model class enjoys great
flexibility and extensibility because arbitrary statistical models can be
incorporated in a mixture which corresponds to fitting a convex combination of
these models to the data.

McLachlan and Peel (2000) and Frühwirth-Schnatter (2006) provide a thorough
and detailed introduction into finite mixture models discussing different
extensions and their application. McLachlan and Peel (2000) focus on maximum
likelihood estimation, while Frühwirth-Schnatter (2006) discusses the Bayesian
estimation in detail. A more recent overview and discussion of new
developments in mixture analysis are compiled in Frühwirth-Schnatter _et al._
(2019). This contribution focuses on Bayesian estimation of finite mixture
models in a model-based clustering context.

The application of finite mixture models in standard data analysis has been
impeded due to the difficulties faced when estimating and performing
inference. Their uptake has been abetted by the increase in available
computing power. Routine estimation of different kinds of mixture models has
become feasible using the expectation-maximization algorithm (Dempster _et
al._ , 1977) for maximum likelihood estimation and Markov chain Monte Carlo
methods in combination with data augmentation (Diebolt and Robert, 1994) for
Bayesian estimation. Both these estimation methods represent a general
estimation framework equally allowing for general extensibility and inclusion
of arbitrary statistical models in the mixture.

Using the Bayesian framework to specify and estimate a finite mixture model
has several advantages:

  * •

Domain knowledge can be included in a principled way using suitable prior
specifications. When using mixture models for cluster analysis, the priors
allow to specify the prototypical shape of the clusters (Hennig, 2015). As an
explanatory analysis tool, cluster analysis aims at detecting interesting
structures in the data. The priors may guide the estimation procedure to focus
on mixture models with certain characteristics. E.g., the priors specified can
reflect and support aims such as obtaining a cluster solution with equally
sized clusters or clusters which differ primarily in their centroids.

  * •

The mixture likelihood is highly irregular due to identifiability issues and
multimodality as well as the presence of spurious modes. The specification of
priors in Bayesian estimation allows to obtain a regularized and smoothed
version of the mixture likelihood as posterior. In this way problems in
maximum likelihood estimation can be avoided where the application of
optimization methods suffers from this irregularity.

  * •

Parameter uncertainty can be easily assessed using the whole posterior
distribution. No reliance on asymptotic normality is required, allowing for
valid inference in cases where regularity conditions are violated, e.g., for
small data sets and mixtures with small component weights.

Despite these advantages, Bayesian finite mixture models are not routinely
applied yet because of uncertainty regarding the specification of suitable
priors and lack of understanding of their impact (Aitkin, 2001). Also issues
with estimation as well as inference are reported due to computational
challenges (Celeux _et al._ , 2019b). In the following all steps of a Bayesian
finite mixture modeling application are explained in detail, including data
model and prior specification, estimation and inference. This is done for both
cases: if the number of components is assumed to be known as well as if this
number is unknown. In addition to outlining the methodology and discussing
theoretical considerations, a real data set serves as an example to illustrate
the application of each of the steps and as a basis to facilitate the uptake
of Bayesian estimation of finite mixture models in applied research.

##  2 Known Number of Components

###  2.1 Model Specification

####  2.1.1 Data Model

When using finite mixture models for clustering the usual assumption is that a
uni- or multivariate continuous or discrete outcome variable $\bm{y}$ is
sampled from a heterogeneous population where $K$ groups are present. Each
group follows its own group-specific parametric distribution. Each observation
in the sample belongs to one of the $K$ groups, but the group memberships are
not observed. It is further assumed that the group sizes in the sample
correspond to the prevalence of the groups in the population.

Marginally, the finite mixture density of an observation $\bm{y}_{i}$ is given
by a convex combination of the group-specific parametric densities with the
weights corresponding to the group sizes:

$\displaystyle h(\bm{y}_{i}|\bm{\vartheta}_{K})$
$\displaystyle=\sum_{k=1}^{K}\eta_{k}f(\bm{y}_{i}|\bm{\theta}_{k}),$

with the group sizes also referred to as component weights

$\displaystyle\eta_{k}$ $\displaystyle>0,\,k=1,\ldots,K,$
$\displaystyle\sum_{k=1}^{K}\eta_{k}=1$

and the group-specific or component density $f(\cdot|\bm{\theta})$ being a
parametric density for $\bm{y}$ with component-specific parameters
$\bm{\theta}$. $\bm{\vartheta}_{K}$ denotes the complete set of parameters of
a finite mixture model with $K$ components consisting of the component weights
$\bm{\eta}_{K}=(\eta_{1},\ldots,\eta_{K})$ and the component-specific
parameters $\bm{\theta}_{k}$, $k=1,\ldots,K$ which are combined to
$\bm{\Theta}_{K}=(\bm{\theta}_{k})_{k=1,\ldots,K}$.

Assuming that the observations $\bm{y}_{i}$, $i=1,\ldots,N$, are independently
identically distributed, the likelihood of the data
$\bm{y}=(\bm{y}_{i})_{i=1,\ldots,N}$ is given by

$\displaystyle p(\bm{y}|\bm{\vartheta}_{K})$
$\displaystyle=\prod_{i=1}^{N}\left[\sum_{k=1}^{K}\eta_{k}f(\bm{y}_{i}|\bm{\theta}_{k})\right].$

In contrast to other model likelihoods the mixture likelihood has the property
that the product over the observations is taken over sums which implies for
the log-likelihood that a sum of the log of sums is determined. This
complicates maximum likelihood estimation.

Other issues with the mixture likelihood are the multimodality of its surface
and the potential presence of spurious modes due to unboundedness for
parameters at the border of the parameter space. Multimodality is caused by
the invariance of the likelihood to relabeling the components i.e., the
likelihood value does not change if the components are re-ordered. As a
consequence, the mixture likelihood function has $K!$ different, but
equivalent modes corresponding to all different ways of labeling. Unbounded
values occur for example for multivariate mixtures of Gaussians in case the
mixture distribution contains a component where the mean value is equal to an
observed value $\bm{y}_{i}$ and the variance-covariance matrix is the zero
matrix.

The structure of the mixture likelihood is complicated because the group
memberships are not observed. If each observation would not only consist of
$\bm{y}_{i}$ but also the group membership $S_{i}\in\\{1,\ldots,K\\}$, the
likelihood for these complete observations
$(\bm{y}_{i},S_{i})_{i=1,\ldots,N}$, which consist of the observed data and
the missing or unobserved data, would be simpler in structure. The complete-
data likelihood is given by

$\displaystyle p(\bm{y},\bm{S}|\bm{\vartheta}_{K})$
$\displaystyle=\prod_{i=1}^{N}\left[\eta_{S_{i}}f(\bm{y}_{i}|\bm{\theta}_{S_{i}})\right].$

The complete-data likelihood is easier to maximize. Dempster _et al._ (1977)
propose to use the expectation-maximization algorithm for maximum likelihood
or maximum a-posteriori estimation exploiting the simplicity of the complete-
data log-likelihood compared to the observed-data log-likelihood.

######  Example 1 (Data set and mixture model specification).

The diabetes data set (Reaven and Miller, 1979) is used in Scrucca _et al._
(2016) to illustrate the use of finite mixture models in a maximum likelihood
context. The data set contains observations for 145 non-obese adults who were
classified into three clinical groups (class) of diabetes ("Normal", "Overt",
"Chemical"). The following measurements are available for the patients:

  * •

glucose: the area under the plasma glucose curve after a three hour oral
glucose tolerance test (OGTT),

  * •

insulin: the area under the plasma insulin curve after a three hour OGTT,

  * •

sspg: the steady state plasma glucose level.

Figure 1: Bivariate scatter plots of the clinical measurements in the data set
diabetes in the off-diagonal with coloring according to the clinical
classification in the lower triangle and kernel density estimates grouped by
clinical classification in the diagonal.

The data is visualized in Figure 1. The upper triangle contains scatter plots
of the observed data $\bm{y}$ which consists of three continuous variables
glucose, insulin and sspg. The unobserved clinical group memberships are
indicated in the pairwise scatter plots in the lower triangle by color and
shape. The diagonal contains for each variable kernel density estimates of the
univariate distributions separately for each clinical group.

In the following it is assumed that the data generating process for this data
set is a finite mixture distribution with three components, each corresponding
to one of the three clinical groups. We assume further that the observations
within each group follow a multivariate normal distribution, i.e., the
component distribution is a three-dimensional multivariate Gaussian
distribution. We will fit a normal mixture model with three components to
assess if the clinical groups can be inferred when observing only the three
clinical measurements. The plot clearly indicates that, in order to capture
the groups, the component distributions need to differ in their mean as well
as their variance-covariance parameters and that using diagonal variance-
covariance matrices for the component distributions would not provide a
suitable fit. Thus in this section, where the number of components $K$ are
assumed to be known, the following data generating process corresponding to a
normal mixture model is imposed:

$\displaystyle
p(\bm{y}|\bm{\eta}_{3},\bm{\mu}_{1},\bm{\mu}_{2},\bm{\mu}_{3},\bm{\Sigma}_{1},\bm{\Sigma}_{2},\bm{\Sigma}_{3})$
$\displaystyle=\sum_{k=1}^{3}\eta_{k}f_{N}(\bm{y}|\bm{\mu}_{k},\bm{\Sigma}_{k}),$

where $f_{N}(\cdot|\bm{\mu},\bm{\Sigma})$ is the multivariate Gaussian
distribution density with mean $\bm{\mu}$ and variance-covariance matrix
$\bm{\Sigma}$.

####  2.1.2 Prior Specification

In the Bayesian framework, all information contained in the data $\bm{y}$
about the parameters $\bm{\vartheta}_{K}$ is summarized in terms of the
posterior density $p(\bm{\vartheta}_{K}|\bm{y})$, which is derived by
combining the likelihood with the prior on the parameters
$p(\bm{\vartheta}_{K})$ using Bayes’ theorem:

$\displaystyle p(\bm{\vartheta}_{K}|\bm{y})$ $\displaystyle\propto
p(\bm{y}|\bm{\vartheta}_{K})p(\bm{\vartheta}_{K}).$

The parameters $\bm{\vartheta}_{K}$ in a finite mixture model with $K$
components consist of the component weights
$\bm{\eta}_{K}=(\eta_{1},\ldots,\eta_{K})$ and the component-specific
parameters $\bm{\Theta}_{K}=(\bm{\theta}_{k})_{k=1,\ldots,K}$.

The prior specification constitutes a crucial part of the finite mixture model
definition as only the joint specification of the data model and the prior
determines the statistical model. No conjugate prior exists for the mixture
likelihood. This implies that no prior specification is possible which would
allow to obtain the posterior in closed form. Thus, usually Markov chain Monte
Carlo (MCMC) methods are used to approximate the posterior.

Despite the need to use MCMC methods for posterior inference anyway,
simplifying assumptions are in general made when specifying priors for finite
mixture models. In particular, in the following we will only consider priors
which have the following characteristics:

  * •

The priors on the component weights and the component-specific parameters are
assumed to be independent. The component-specific parameters may also depend
on some hyper-parameters $\bm{\phi}$ but are conditionally independent. This
implies the following prior structure:

$\displaystyle p(\bm{\eta}_{K},\bm{\Theta}_{K},\bm{\phi})$
$\displaystyle=p(\bm{\eta}_{K})p(\bm{\Theta}_{K}|\bm{\phi})p(\bm{\phi})=p(\bm{\eta}_{K})\prod_{k=1}^{K}p(\bm{\theta}_{k}|\bm{\phi})p(\bm{\phi}).$

Assuming independence simplifies the specification of the prior, in particular
given that usually no prior knowledge is available which would suggest a
specific dependency structure. Hyper-parameters on the component-specific
parameters allow to dilute the influence of the specific priors imposed on the
component-specific parameters and adapt them to the data by combining
information across the components. Again this is a useful strategy in case no
clear prior knowledge is available which can be exploited to specify
particular priors on the component-specific parameters.

  * •

Proper priors are imposed on the component-specific parameters. Improper
priors on the component-specific parameters lead to improper posteriors. This
can easily be seen through the following considerations: The posterior is
proportional to the likelihood times prior and the likelihood of $\bm{y}$
given $\bm{\vartheta}_{K}$ can also be determined by summing the complete-data
likelihood of $\bm{y}$ and $\bm{S}$ given $\bm{\vartheta}_{K}$ over all
possible assignment vectors $\bm{S}$. In this sum, assignment vectors $\bm{S}$
which do not assign observations to all components lead to empty components
where the posterior is identical to the prior. Hence an improper prior leads
to an improper posterior.

  * •

The priors are label-invariant. This implies that a-priori the distribution of
each component is the same. Having label-invariant priors in combination with
the label-invariant likelihood leads to a posterior which is also label-
invariant. This is a sensible approach in case no clear prior information is
available which would allow to differently characterize the components of the
mixture.

  * •

The priors are conditionally conjugate. This implies that, conditional on the
data, the latent variables and all other parameters, the marginal posteriors
of the parameters are available in closed form. This facilitates the
application of a Gibbs sampling scheme in MCMC sampling for inference. In
addition to the computational advantages this also leads to the use of well
known parametric priors where the characteristics are well understood. While
this specification restricts the parametric families considered for the
priors, these priors still offer great flexibility by varying the parameter
values or adding a hyper-prior. For an empirical example illustrating the
influence of different parameter values on the results obtained see Grün _et
al._ (2021).

Using priors with these characteristics implies that the prior for the
component weights is specified independently of the component-specific
distribution and hence also is the same regardless of the specific mixture
model. To select the prior on the component weights $\bm{\eta}_{K}$, one uses
that the unobserved group memberships or assignment vectors $\bm{S}$ are drawn
from a multinomial distribution. The conjugate prior is the Dirichlet
distribution. To obtain a label-invariant prior a symmetric $K$-dimensional
Dirichlet distribution is specified which depends only on a single scalar
parameter $\gamma_{K}$:

$\displaystyle\bm{\eta}_{K}$ $\displaystyle\sim\mathcal{D}(\gamma_{K}).$

The prior mean of the component weights is thus label-invariant and given by
$1/K$. The parameter $\gamma_{K}$ influences the prior variability of the
component weights with a prior variance of $(K-1)/(K^{2}(K\gamma_{K}+1))$.

The priors for the component-specific parameters $\bm{\theta}_{k}$,
$k=1,\ldots,K$ depend on the component-specific distribution of the mixture.
They can be selected as the usual conjugate priors of the component-specific
distribution ignoring the mixture structure. For the prior on the hyper-
parameters $\bm{\phi}$ the conjugate prior for the set of component-specific
parameters is an obvious choice.

######  Example 2 (Prior specification).

To fit a finite mixture of multivariate Gaussian distributions with three
components to the diabetes data set, we need to specify $\gamma_{K}$ and
priors for the component-specific means and the component-specific variance-
covariance matrices. We choose $\gamma_{K}=1$ which corresponds to imposing a
uniform prior on the weights. If the number of components is known, this
parameter is usually not influential.

For the component-specific parameters we will use conditionally conjugate
priors, i.e., a normal prior on the component means and an inverse Wishart
prior on the variance-covariance matrices with a hierarchical prior on the
second parameter. This specification is similar to the specification proposed
in Richardson and Green (1997) for univariate Gaussian mixtures. This implies:

$\displaystyle\bm{\mu}_{k}$
$\displaystyle\sim\mathcal{N}(\bm{b}_{0},\bm{B}_{0}),k=1,\ldots,K,$
$\displaystyle\bm{\Sigma}_{k}$
$\displaystyle\sim\mathcal{W}^{-1}(c_{0},\bm{C}_{0}),k=1,\ldots,K,$
$\displaystyle\bm{C}_{0}$ $\displaystyle\sim\mathcal{W}(g_{0},\bm{G}_{0}).$

Note that we use the parameterization of the Wishart and inverse Wishart
distribution as employed by Frühwirth-Schnatter (2006). This means that the
density of $\bm{Y}\sim\mathcal{W}(\alpha,\bm{V})$ is given by

$\displaystyle f(\bm{Y}|\alpha,\bm{V})$
$\displaystyle=\frac{|\bm{V}|^{\alpha}}{\Gamma(\alpha)}|\bm{Y}|^{\alpha-(r+1)/2}\exp\\{-\text{tr}(\bm{V}\bm{Y})\\},$

and of $\bm{Y}\sim\mathcal{W}^{-1}(\alpha,\bm{V})$ by

$\displaystyle f(\bm{Y}|\alpha,\bm{V})$
$\displaystyle=\frac{|\bm{V}|^{\alpha}}{\Gamma(\alpha)}|\bm{Y}^{-1}|^{\alpha+(r+1)/2}\exp\\{-\text{tr}(\bm{V}\bm{Y}^{-1})\\},$

where
$\Gamma(\alpha)=\pi^{r(r-1)/4}\prod_{j=1}^{r}\Gamma(\frac{2\alpha+1-j}{2})$
and $r$ is the dimension of the observations, i.e.,
$\bm{y}_{i}\in\mathbb{R}^{r}$.

We select the prior parameters $\bm{b}_{0}$ and $\bm{B}_{0}$ to ensure that
the analysis is invariant to the scaling of the data and to have proper priors
which are barely informative to induce no shrinkage on the component-specific
means to a common mean value. In particular we set
$\bm{b}_{0}=\text{median}(\bm{y})$, i.e., equal to the component-wise median,
and $\bm{B}_{0}=\bm{R}^{2}$, i.e., equal to the squared data range, in order
to be very uninformative in regard to $\bm{b}_{0}$.

Further, we define the degree of freedom parameters of the Wishart and inverse
Wishart distribution to depend on the dimension $r$ of the data
$\bm{y}_{i}\in\mathbb{R}^{r}$:

$\displaystyle c_{0}$ $\displaystyle=c+\frac{r+1}{2},$ $\displaystyle g_{0}$
$\displaystyle=1+\frac{r-1}{2},$

where $c$ is a constant. $\bm{C}_{0}$ is modeled in order to incorporate
knowledge about the a-priori expected component variance-covariance matrix
$\bm{\Sigma}_{k}$. To be scale-invariant the prior mean of $\bm{\Sigma}_{k}$
is assumed to be a fraction $\phi$ of the diagonal matrix $\bm{S}$, which
contains the diagonal values of the empirical variance-covariance matrix of
the data $\textsf{COV}(\bm{y})$ in the diagonal. Specifically, for the
following specifications

$\displaystyle\bm{C}_{0}$ $\displaystyle=c\phi\bm{S},$
$\displaystyle\bm{G}_{0}$ $\displaystyle=g_{0}\bm{C}_{0}^{-1},$

the prior mean value of $\bm{\Sigma}_{k}$ is $\phi\bm{S}$. A sensible choice
for $\phi$ is, e.g., $\phi=0.75$. The values of $c$ and $\phi\bm{S}$ can be
interpreted as adding $c$ virtual observations with variance-covariance matrix
$\phi\bm{S}$ to each group $k$. If $c=2.5$, the eigenvalues of
$\bm{\Sigma}_{k}\bm{\Sigma}_{j}^{-1}$ are bounded away from zero avoiding
spurious modes in the posterior (Frühwirth-Schnatter, 2006, p. 192).

####  2.1.3 Hierarchical Generative Model

Combining the data model with the prior specifications, the hierarchical
generative model for the parameters $\bm{\vartheta}_{K}$, the unobserved
assignment vectors $\bm{S}$ and the observed data $\bm{y}$ is given by:

$\displaystyle\bm{\eta}_{K}|K,\gamma_{K}$
$\displaystyle\sim\mathcal{D}(\gamma_{K}),$ (1) $\displaystyle\bm{\phi}$
$\displaystyle\sim p(\bm{\phi}),$ (2) $\displaystyle\bm{\theta}_{k}|\bm{\phi}$
$\displaystyle\sim p(\bm{\theta}_{k}|\bm{\phi}),\text{ independently for
}k=1,\ldots,K,$ (3) $\displaystyle S_{i}|K,\bm{\eta}_{K}$
$\displaystyle\sim\mathcal{M}(1;\bm{\eta}_{K}),\text{ independently for
}i=1,\ldots,N,$ (4) $\displaystyle\bm{y}_{i}|K,S_{i}=k,\bm{\theta}_{k}$
$\displaystyle\sim f(\bm{y}_{i}|\bm{\theta}_{k}),\text{ independently for
}i=1,\ldots,N,$ (5)

where $\mathcal{D}(\gamma_{K})$ is the symmetric $K$-dimensional Dirichlet
distribution with scalar parameter $\gamma_{K}$ and
$\mathcal{M}(1;\bm{\eta}_{K})$ is the multinomial distribution with one trial
and vector of success probabilities $\bm{\eta}_{K}$ for the $K$ possible
outcomes.

###  2.2 Model Estimation

Different computational solutions for estimating Bayesian finite mixture
models with a fixed number of components $K$ and obtaining an approximation of
the posterior of $\bm{\vartheta}_{K}$ are considered in Celeux _et al._
(2019b). A generic approach builds on data augmentation and includes the
component assignment vectors $\bm{S}$ in the sampling scheme (Diebolt and
Robert, 1994). Using conditionally conjugate priors, the resulting MCMC
sampling scheme only requires Gibbs steps. More details on MCMC methods and
the Metropolis–Hastings algorithm in general are given in [stat08285.
stat07834].

Given some initial parameter values $\bm{\eta}_{K}$ and $\bm{\Theta}_{K}$,
iterate Steps 1 and 2:

  1. 1.

Classification of each observation $\bm{y}_{i}$, $i=1,\ldots,N$ conditional on
$\bm{\eta}_{K}$ and $\bm{\Theta}_{K}$:

    1. (a)

For $i=1,\ldots,N$: sample $S_{i}$ from

$\displaystyle
P(S_{i}=k|\bm{y}_{i},\bm{\vartheta}_{K})\propto\eta_{k}f(\bm{y}_{i}|\bm{\theta}_{k}).$

  2. 2.

Parameter simulation conditional on the classification
$\bm{S}=(S_{1},\ldots,S_{N})$:

    1. (a)

Sample the mixture weights $\bm{\eta}_{K}$ from

$\displaystyle\bm{\eta}_{K}|\bm{S}\sim\mathcal{D}(e_{1},\ldots,e_{K}),$

where

$\displaystyle e_{k}=\gamma_{K}+N_{k}$

with $N_{k}=\\#\\{i:S_{i}=k\\}$ the number of observations assigned to
component $k$.

    2. (b)

Sample the component-specific parameters $\bm{\theta}_{k}$, $k=1,\ldots,K$
conditional on $\bm{y}$, $\bm{S}$ and $\bm{\phi}$ from their conditional
conjugate posteriors.

    3. (c)

Sample hyper-parameters $\bm{\phi}$ conditional on $\bm{\Theta}_{K}$.

######  Example 3 (MCMC sampling scheme).

To fit the finite mixture of multivariate Gaussian distributions to the
diabetes data set, Steps 2(b) and 2(c) are specifically given by:

  1. 2(b)

For all $k=1,\ldots,K$, sample

$\displaystyle\bm{\mu}_{k}|\bm{y},\bm{S},\bm{\Sigma}_{k}$
$\displaystyle\sim\mathcal{N}(\bm{b}_{k},\bm{B}_{k}),$
$\displaystyle\bm{\Sigma}_{k}^{-1}|\bm{y},\bm{S},\bm{\mu}_{k}$
$\displaystyle\sim\mathcal{W}(c_{k},\bm{C}_{k}),$

with

$\displaystyle\bm{b}_{k}$
$\displaystyle=\bm{B}_{k}(\bm{B}_{0}^{-1}\bm{b}_{0}+\bm{\Sigma}_{k}^{-1}N_{k}\bar{\bm{y}}_{k}),$
$\displaystyle\bm{B}_{k}$
$\displaystyle=(\bm{B}_{0}^{-1}+N_{k}\bm{\Sigma}_{k}^{-1})^{-1},$
$\displaystyle c_{k}$ $\displaystyle=c_{0}+\frac{N_{k}}{2},$
$\displaystyle\bm{C}_{k}$
$\displaystyle=\bm{C}_{0}+\frac{1}{2}\sum_{i:S_{i}=k}(\bm{y}_{i}-\bm{\mu}_{k})(\bm{y}_{i}-\bm{\mu}_{k})^{\top},$

where $\bar{\bm{y}}_{k}=\frac{1}{N_{k}}\sum_{i:S_{i}=k}\bm{y}_{i}$.

  2. 2(c)

Sample

$\displaystyle\bm{C}_{0}|(\bm{\Sigma}_{k})_{k=1,\ldots,K}\sim\mathcal{W}(g_{0}+Kc_{0},\bm{G}_{0}+\sum_{k=1}^{K}\bm{\Sigma}_{k}^{-1}).$

See Frühwirth-Schnatter (2006, Chapter 6).

####  2.2.1 Initialization

A standard default initialization strategy for MCMC sampling is to draw the
initial parameter values from their priors. For Bayesian finite mixture models
this strategy is not advisable because such an initial parameter vector
$\bm{\vartheta}_{K}$ might have a very low posterior probability leading to
many burn-in iterations being required to arrive at a parameter region with
higher posterior probability. In addition one usually also is at risk to be
trapped in a local optimum.

To obtain a better initial parameter vector, in particular in the context of
cluster analysis, the initial parameter values $\bm{\eta}_{K}$ and
$\bm{\Theta}_{K}$ may be obtained based on a partition of the data set which
is derived using some cluster analysis method. E.g., $k$-means clustering can
be applied to the data $\bm{y}$ with $k$-means clustering implicitly assuming
equally sized clusters and multivariate normal cluster distributions with
spherical variance-covariance matrices which are identical across clusters.

######  Example 4 (Initialization).

For the diabetes data set, the initial partition of the data into three groups
is obtained by applying $k$-means clustering. The initial values for the
component means are set to the cluster means of the initial partition, whereas
the initial variance-covariance matrices are set to $0.75\bm{S}$ with $\bm{S}$
being a diagonal matrix containing the diagonal elements of
$\textsf{COV}(\bm{y})$ and the mixture weights set equal to $1/K$.

####  2.2.2 Mixing and Convergence Diagnostics

Using MCMC with data augmentation for fitting Bayesian mixture models is known
to be susceptible to poor mixing (Celeux _et al._ , 2000). In fact, the known
multimodality of the mixture posterior due to label invariance may be
exploited to assess convergence. Using label-invariant priors, the posterior
takes the same values regardless of how the labels are permuted. If a mixture
with $K$ components is fitted, this means that the posterior has $K!$
identical modes induced by all possible permutations of the labels. Thus a
well mixing Markov chain should equally visit all these modes and all
posteriors of component-specific parameters should be identical across $K$. To
resolve the problem of partial label switching and ensure that all these modes
are equally visited, Frühwirth-Schnatter (2001) proposed the permutation
sampler which ensures balanced label-switching.

Visual inspection of the trace plots of the posterior draws of the parameters
as well as other quantities of interest such as the mixture likelihood allow
to assess convergence of the chain as well as the presence of label switching
in case non-label-invariant quantities are inspected. Insights from these
plots can be used to determine a suitable number of burn-in iterations to be
omitted and decide if more iterations are required for suitable posterior
inference. If the visual inspection indicates that during MCMC sampling only
one mode of the mixture posterior is explored, there is in fact no need to
resolve label switching. However, as label switching cannot be always assessed
visually, it is common practice to assume the presence of label switching and
resolve it in any case before performing component-specific inference.

######  Example 5 (Trace plots).

The Gibbs sampling scheme outlined in Section 2.2 and Example 3 is used to
draw $M=30000$ samples from the finite Gaussian mixture model specified in
Examples 1 and 2 to analyze the diabetes data set. All draws are retained
without omitting any burn-in iterations in order to study the convergence
behavior of the chain. Figure 2 shows trace plots of the sampled component
means of the first variable and the number of allocated observations to each
of the components for the first 10000 iterations. The figure indicates that
the chain converged rather fast after about 1000 iterations. Also, it seems
that no label switching occurred after convergence, at least for the plotted
range of iterations. In the following the first 5000 iterations are omitted as
burn-in and the remaining 25000 iterations used for posterior inference.

Figure 2: Trace plots of the $\mu_{k1}$ draws of the first variable (left) and
the number of observations allocated to a component $N_{k}$ (right), in
different colors for the three components.

###  2.3 Post-Processing and Inference

####  2.3.1 Resolving label switching

The approximations of the posterior distributions of parameters and latent
variables should be identical across the $K$ components in case the chain has
converged and equally visited all modes induced by the $K!$ permutations of
the labels. However, these approximations do not help to characterize the
component sizes and distributions or assign observations to components. Hence,
parameters which are label-dependent need to be relabeled in order to focus on
a single mode of the posterior. Many different methods have been proposed to
obtain a unique labeling (for an overview, see for example Papastamoulis,
2016). For example, a unique labeling can be achieved by imposing ordering
constraints on the draws or by solving an optimization problem to minimize the
dissimilarity between the classification matrices and some average of them.

Alternatively, a straightforward approach is to apply $k$-means clustering to
the point process representation of the MCMC draws to identify the mixture
model (Frühwirth-Schnatter, 2006; Malsiner-Walli _et al._ , 2016). We will
refer to this approach as the ppr approach in the following. The component-
specific parameters $\bm{\theta}_{k}$, $k=1,\ldots,K$, or some (low-
dimensional) functional, are clustered in the point process representation
into $K$ clusters using $k$-means clustering. The information which draws are
from which component or from which MCMC iteration is not included in the
clustering procedure. Hence, even though the draws are clustered into $K$
groups, not necessarily all $K$ sampled component-specific parameter vectors
of one MCMC sweep are assigned to $K$ different clusters. Only for those MCMC
sweeps where the $K$ sampled parameter vectors are clustered into exactly $K$
groups a unique labeling of the draws can be obtained. Thus, only these MCMC
draws are retained and the other MCMC draws are omitted before relabeling and
performing a component-specific analysis. The proportion of MCMC sweeps which
need to be omitted serves as indicator how well the components can be
distinguished. A high proportion of omitted sweeps usually results if draws of
different components are very similar and therefore allocated to the same
cluster by $k$-means. This is in particular in a cluster analysis context an
indication that the fitted mixture model is overfitting and has too many
components.

######  Example 6 (Resolving label switching).

In order to obtain a unique labeling of the sampled draws, the sampled
component means $\bm{\mu}_{k}$ are clustered into three groups using
$k$-means. The bivariate scatter plots of the sampled component means
$\bm{\mu}_{k}$ are shown in Figure 3 in the first row. After applying
$k$-means clustering the classification is shown in the second row. The
sampled component means $\bm{\mu}_{k}$ are clearly allocated to three separate
and distinct clusters. The non-permutation rate, i.e., the proportion of MCMC
sweeps where the draws are clustered into less than three different clusters,
is smaller than 0.01. Thus, for more than 99% of the MCMC sweeps a unique
labeling could be obtained and the remaining less than 1% sweeps were excluded
from further analysis. The obtained unique labeling is used to reorder not
only the component mean draws $\bm{\mu}_{k}$, but also the draws of the other
parameters, e.g., the mixture weights $\bm{\eta}_{k}$ and the assignment
variables $\bm{S}$. After relabeling, component-specific inference is possible
by, e.g., computing the posterior means. In Table 1 on the left the posterior
means of the mixture weights and the component means are reported for the
model identified using the ppr approach.

Figure 3: Bivariate scatter plots of the $\mu$ draws for all three components and for the three variables (on the top) and the retained $\mu$ draws after model identification with the ppr approach where the colors indicate the identified component labels (on the bottom). Only the first 1000 iterations are shown after discarding the first 5000 iterations as burn-in.  Table 1: Component- or cluster-specific inference of the identified mixture models estimated with known $K=3$ (left), using a sparse finite mixture (SFM) model with $K=10$ and $\gamma=0.01$ (middle) and a dynamic mixture of finite mixtures (MFM) model with BNB$(1,4,3)$ prior on $K$ and $\gamma_{K}=0.5/K$ (right): Number of observations allocated to the components or clusters using the posterior mode of the assignment probabilities ($N_{k}$) and the posterior mean estimates of the component or cluster weights $\hat{\eta}_{k}$ and component or cluster means $\hat{\bm{\mu}}_{k}$, $k=1,2,3$.  |  Known $K=3$ |  |  SFM with $K=10$ |  | Dynamic MFM
---|---|---|---|---|---
| 1 | 2 | 3 |  | 1 | 2 | 3 |  | 1 | 2 | 3
$N_{k}$ | 28 | 33 | 84 |  | 28 | 33 | 84 |  | 28 | 33 | 84
$\hat{\eta}_{k}$ | 0.20 | 0.25 | 0.55 |  | 0.20 | 0.24 | 0.56 |  | 0.20 | 0.24 | 0.56
$\hat{\mu}_{k1}$ | 229.41 | 104.37 | 91.41 |  | 229.39 | 104.49 | 91.44 |  | 229.41 | 104.49 | 91.45
$\hat{\mu}_{k2}$ | 1098.04 | 496.87 | 361.43 |  | 1097.89 | 497.94 | 361.73 |  | 1097.97 | 497.78 | 361.89
$\hat{\mu}_{k3}$ | 82.66 | 319.27 | 165.19 |  | 82.72 | 321.17 | 165.47 |  | 82.71 | 321.89 | 165.44

####  2.3.2 Final partition

After having resolved the label switching issue, the posterior of the
assignment vectors $\bm{S}$ can be used to assign the observations to
components and derive a partition of the data. For example, the partition
based on the maximum a-posteriori estimate is obtained if each observation
$\bm{y}_{i}$ is assigned to the component to which it is allocated most often
during MCMC sampling, i.e., based on the posterior mode. Note that the
partition derived in this way may contain less than $K$ components.

Alternatively, if interest only lies in obtaining a final partition, loss-
based approaches can be used where the loss is label-invariant and hence no
unique labeling is required and thus resolving the label switching can be
avoided. This approach requires to solve an optimization problem which
consists of determining the partition which minimizes the posterior loss
(Rastelli and Friel, 2018). Yet another possibility is to estimate a
similarity matrix between observations based on their posterior co-allocation
probabilities and then use clustering methods to obtain a suitable partition
based on the similarity matrix. Molitor _et al._ (2010) apply partitioning
around medoids (PAM) to such a similarity matrix and select a suitable number
of clusters based on the silhouette width.

######  Example 7 (Final partition).

The final partition of the data is obtained by relabeling the assignment
vectors $\bm{S}$ according to the unique labeling obtained in the ppr approach
and then assigning each observation to the component to which it was assigned
most frequently during MCMC sampling. The resulting group sizes are reported
in Table 1 on the left. Table 2 on the left reports the confusion table of the
estimated and true classification. The estimated cluster labels seem to be in
line with the true class labels if both are sorted by group size and hence,
the observations on the diagonal are assumed to be correctly classified. Based
on the number of observations in the off-diagonal a misclassification rate
(MCR) of 14% is obtained (see Table 2 on the right). An additional measure to
assess congruence between two partitions is the adjusted Rand index (ARI;
Hubert and Arabie, 1985). The adjusted Rand index has the advantage that it is
label invariant. The maximum value of one indicates perfect congruence between
two partitions, the adjustment is used to obtain a value of zero if the number
of co-allocations in the pair of partitions equals the number expected by
chance. For the example data set the final partition estimated has an adjusted
Rand index value of 0.65 (see Table 2 on the right). Determining the final
partition based on the minimization of the posterior loss defined by the
variation of information criterion (Rastelli and Friel, 2018) yields almost
the same result with respect to the performance measures ARI and MCR. The
final partition of the ppr approach is also visualized and compared to the
true classification using bivariate scatter plots of the measurements in
Figure 4.

Table 2: Left: Confusion table of the true and the estimated classification
obtained with the ppr approach. Right: Partition congruence measures
consisting of the adjusted Rand index (ARI) and the misclassification rate
(MCR) for final partitions, obtained by the ppr method (ppr) or minimizing the
posterior loss based on the variation of information method (VI).

| Estimated
---|---
True | 1 | 2 | 3
Overt | 27 | 6 | 0
Chemical | 1 | 24 | 11
Normal | 0 | 3 | 73

Method | Groups | ARI | MCR
---|---|---|---
ppr | 3 | 0.65 | 0.14
VI | 3 | 0.64 | 0.15

Figure 4: Bivariate scatter plots of the clinical measurements with the colors
indicating the component membership estimated using the posterior modes of the
component assignments after identification with the ppr approach and the
shapes the true classes. The size of misclassified observations is increased.

##  3 Unknown Number of Components

In many applications the number of groups is not known. In this case data
analysis does not only consist of fitting a mixture model with a fixed number
of components $K$, but usually also needs to provide insights into how many
data clusters are present in the data set. Celeux _et al._ (2019a) point out
that determining the number of components to include in a finite mixture model
is not only an important question because inference about the fitted model is
highly sensitive to the $K$ selected. Determining $K$ is also most difficult
due to the identifiability issues of the mixture model and due to the fact
that a suitable selection of $K$ might also depend on the modeling purpose of
the mixture model.

An important distinction to make when estimating the number of clusters in the
data set at hand, is the conceptual difference between the number of
components $K$ in the finite mixture model and the number of clusters in the
data set, henceforth called $K_{+}$. McCullagh and Yang (2008) discuss this
distinction in the context of species sampling where $K$ represents the number
of species in the population whereas $K_{+}$ corresponds to the number of
species observed in the available data set. These numbers need not to be the
same as not all species of the population might also be observed in the
sample. Correspondingly, not all available mixture components might be used to
generate the data at hand. Thus, it might be that $K_{+}<K$.

Maximum likelihood estimation of finite mixture models does not make this
distinction between $K$ and $K_{+}$, as the usual assumption is that the
number of clusters $K_{+}$ in the data exactly corresponds to the number of
mixture components $K$. Problems with empty components are usually only
acknowledged if, for example, uncertainty assessment based on bootstrapping is
performed, see, e.g., O’Hagan _et al._ (2019) who propose to use a weighted
likelihood bootstrap to ensure that observations from all components present
in the original data set are also present in the bootstrap samples.

In contrast, in Bayesian estimation of finite mixture models the distinction
between $K$ and $K_{+}$ is always present. This distinction is of course made
in Bayesian nonparametrics [stat07850, stat07905] where an infinite number of
mixture components is assumed while the data at hand only consists of a finite
number of clusters. However, also in finite mixture modeling the distinction
between $K$ and $K_{+}$ is natural. When assigning observations to the
components by drawing from the multinomial distribution with success
probabilities equal to the component weights, there is a-priori a positive
probability that not to all components observations are assigned. This
probability is the larger, the more unbalanced the weights are. A-priori the
balancedness of the weights is influenced by their prior, i.e., the Dirichlet
parameter $\gamma_{K}$.

Unfortunately, the handling of an overfitting mixture model where the number
of mixture components $K$ exceeds the number of clusters $K_{+}$ in the data
set is not straightforward as identifiability problems arise. Additional
components can be added to any mixture model without changing the induced
mixture distribution: either a component weight of zero is assigned to the new
components or the new components are obtained by duplicating an already
existing component and splitting the component weight between them. Rousseau
and Mengersen (2011) provide theoretic insights into how the value
$\gamma_{K}$ of the Dirichlet prior on the mixture weights impacts on the
asymptotic behavior of the posterior of an overfitting finite mixture model.
Frühwirth-Schnatter (2011) points further out how these theoretic results can
be used for practical data analysis. In fact, if $K$ is not known and has to
be estimated, $\gamma_{K}$ needs to be specified in dependence of the model
selection strategy pursued, as will be explained in the following sections.
The choice of $\gamma_{K}$ determines how large the gap between $K$ and
$K_{+}$ is a-priori, i.e., the prior assumption on how many of the species in
the population are also observed in the data set.

In the following we discuss three strategies for Bayesian finite mixture
modeling when the number of components $K$ is unknown and one aims at
estimating the number of clusters $K_{+}$ in the data set. Note that the
problem of estimating the number of groups in a data set is usually referred
to as “estimating $K$”. However, we explicitly pursue a distinction between
$K$ and $K_{+}$ and aim at estimating the number of data clusters $K_{+}$,
which may be different from $K$, the number of mixture components or number of
clusters in the population.

###  3.1 Model Selection

If it is assumed that $K_{+}=K$, estimating $K_{+}$ can be viewed as a model
selection problem, where the most appropriate mixture model within the class
of mixture models with varying $K$ is selected. The posterior probability of
each model is determined by calculating the marginal likelihood in combination
with a suitable prior on the models. Such an approach has the advantage to not
only allow comparisons across mixture models with different number of
components, but also mixture models with different component distributions and
different prior specifications. Using equal prior model probabilities the
ratio between two posterior model probabilities corresponds to the Bayes
factor which also is often used as basis for model selection and to compare
different models.

If the model selection strategy aims at fitting finite mixture models with
different number of components $K$ and comparing their fit using, for example,
the posterior probabilities and the marginal likelihoods, it is important to
ensure that all components are filled and hence to choose the Dirichlet
parameter $\gamma_{K}$ large (Frühwirth-Schnatter, 2011). Determining the
marginal likelihood for finite mixture models is again complicated by the
multimodality due to label switching. Frühwirth-Schnatter (2019) suggests to
estimate the marginal likelihood with the full permutation or the double
random permutation bridge sampling estimator, however only for $K$ up to
seven.

###  3.2 Sparse Finite Mixture Models

A computational less intensive approach to estimate the number of clusters in
the data set is the “sparse finite mixture” (SFM) approach. Inspired by the
theoretic insights provided by Rousseau and Mengersen (2011), Malsiner-Walli
_et al._ (2016) propose the SFM model where $K$ is selected large to induce on
purpose an overfitting finite mixture model. This model specification is
combined with a (very) small value of $\gamma_{K}$ such that for the
superfluous components the component weights a-priori and a-posteriori tend to
zero. This implies that a-priori a large gap between $K$ and $K_{+}$ is
imposed and the posterior of $K_{+}$ may flexibly adapt to the observed
clusters in the data. The number of data clusters $K_{+}$ can be estimated
using the posterior mode of the number of filled components induced by the
sampled assignment vectors. In contrast to pursuing a model selection strategy
as explained in Section 3.1, fitting a sparse finite mixture model does not
require to fit finite mixture models with increasing number of components $K$,
but a single MCMC run of such an overfitting mixture model is sufficient to
estimate the number of data clusters $K_{+}$.

If the posterior mode of the number of filled components is used to select a
suitable number of data clusters $K_{+}$, cluster-specific posterior inference
can be performed after omitting all MCMC sweeps where more or less components
were filled and retaining only MCMC sweeps where exactly $K_{+}$ components
were filled. Using the draws of these MCMC sweeps first only the draws of
component-specific parameters associated with filled components are retained
and those of empty components omitted. These sweeps and draws are then used to
resolve label-switching and perform inference for the parameters of the
cluster-specific distributions. The link between sparse finite mixture models
and Dirichlet process mixtures where $K$ is in fact equal to infinity is
empirically investigated in Frühwirth-Schnatter and Malsiner-Walli (2019).

###  3.3 Prior on the Number of Components

In Bayesian analysis a natural approach to account for the unknown number of
components $K$ is to introduce a prior on $K$ and include the determination of
the posterior of $K$ as well as the posterior of $K_{+}$ in the Bayesian
inference. The hierarchical generative model given in (1) then has the
additional layer at the top:

$\displaystyle K$ $\displaystyle\sim p(K),$

where $p(K)$ is the prior on the number of components. Following Miller and
Harrison (2018) we call the resulting model a mixture of finite mixtures (MFM)
model. Frühwirth-Schnatter _et al._ (2020) review previous suggestions for the
prior on $K$ and propose to use the beta-negative-binomial (BNB) distribution
for $K-1$ to ensure that the support is equal to the positive integer numbers.
The BNB distribution is a three-parameter distribution which allows to not
only calibrate the prior mean and variance, but also take into account the
probability for the homogeneity model $K=1$ and the tail behavior, i.e., the
mass assigned to large values of $K$. The BNB distribution has the Poisson,
geometric and negative-binomial distribution as special cases.

The mixture of finite mixtures model has been investigated by Richardson and
Green (1997) with inference based on the reversible jump MCMC for mixtures of
univariate Gaussians. Specific extensions of this reversible jump MCMC
algorithm to also cover other component distributions were considered in Zhang
_et al._ (2004) and Papastamoulis and Iliopoulos (2009). As an alternative to
reversible jump MCMC Stephens (2000) proposes a sampling scheme based on
creating a Markov birth-death process. Further Miller and Harrison (2018)
suggest to re-use methods from Bayesian nonparametrics and apply the Jain-Neal
sampler developed for Dirichlet process mixtures (Jain and Neal, 2004, 2007).
In this sampling scheme the number of components $K$ is not included but each
observation in turn is re-assigned to an already filled component or a new
component during MCMC sampling, thus effectively sampling partitions of the
data and hence $K_{+}$. To fit the mixture of finite mixtures model,
Frühwirth-Schnatter _et al._ (2020) propose the “telescoping sampler” which
extends the usual MCMC scheme with data augmentation by including also a step
where $K$ is explicitly sampled conditional on the current partition of the
data and another step to obtain parameters for the empty components:

  1. 3.

Sample $K$ conditional on current cluster sizes $(N_{1},\ldots,N_{K_{+}})$ and
$\gamma_{K}$ from

$\displaystyle p(K|(N_{1},\ldots,N_{K_{+}}),\gamma_{K})$
$\displaystyle\propto\frac{K!}{(K-K_{+})!}\frac{\Gamma(K\gamma_{K})}{\Gamma(K\gamma_{K}+N)}\prod_{k=1}^{K_{+}}\frac{\Gamma(N_{k}+\gamma_{K})}{\Gamma(1+\gamma_{K})}p(K).$

  2. 4.

Add $K-K_{+}$ empty components with component-specific parameters sampled from
the priors.

Step 2 is modified to draw the mixture weights only after the empty components
were added and the hyper-parameters are drawn conditional on the component-
specific parameters of the filled components, disregarding the empty
components.

This sampling scheme is very generic allowing for arbitrary component
distributions by just re-using sampling schemes already available for finite
mixtures with a fixed number of components $K$. In contrast to Richardson and
Green (1997) and Miller and Harrison (2018) the telescoping sampler is also
applicable with arbitrary values selected for $\gamma_{K}$. Frühwirth-
Schnatter _et al._ (2020) consider in particular the “dynamic” specification
of a mixture of finite mixtures model where $\gamma_{K}=\alpha/K$, i.e., the
Dirichlet parameter decreases with increasing number of components $K$
implying a larger gap between $K$ and $K_{+}$ for larger values of $K$.

######  Example 8 (Fitting a SFM and MFM model).

A sparse finite mixture model with $K=10$ mixture components and $\gamma=0.01$
is fitted to the data set. The sampler is run for $M=30000$ iterations. The
trace plot of the number of components $K$ (in gray) and filled components
$K_{+}$ (in black) is shown in Figure 5a) on the left. After starting with an
initial partition of 10 filled components, rather fast most of them become
empty and the sampler switches mostly between three and four filled
components. On the right, the relative frequencies of the different $K_{+}$
values sampled after omitting a burn-in of 5000 iterations is plotted. A clear
mode at $K_{+}=3$ is discernible. Thus three clusters are estimated for the
data set and all MCMC sweeps where the number of filled components is not
equal to three are discarded. After discarding also the draws sampled from
empty components the ppr approach is used for the remaining draws to identify
the three clusters. In Table 1 (middle) the posterior means of the cluster
means and cluster weights of the identified model are reported as well as the
group sizes of the final partition. The results correspond to those when
fitting a mixture model with exactly three components and the same final
partition is obtained.

Figure 5: a) Sparse finite mixture with $K=10$, $\gamma=0.01$. b) Dynamic MFM
with $K\sim\text{BNB}(1,4,3)$, $\alpha=0.5$. For each modeling approach, on
the left, trace plot of the number of components (gray) and the number of
filled components (black) and on the right, bar plot of the relative
frequencies of the number of filled components.

Finally, a dynamic MFM model with the prior $K-1\sim\text{BNB}(1,4,3)$ and
$\gamma_{K}=0.5/K$ is fitted to the data using the telescoping sampler. The
trace plots of the sampled $K$ (in gray) and induced $K_{+}$ (in black) are
shown in Figure 5b) on the left. The sampled $K$ values assume large values,
even greater than 20. The sampled $K_{+}$ values again switch between three
and four filled components. The bar plot in Figure 5b) on the right visualizes
the relative frequencies of $K_{+}$ values sampled after omitting 5000 burn-in
iterations. Again the mode of the posterior of $K_{+}$ is at $K_{+}=3$,
although not so clear as when fitting a SFM model. After identification of the
model for $K_{+}=3$, the posterior means of the cluster-specific parameters
and the group sizes of the final partition are reported in Table 1 (right).
Again the results obtained correspond to those for $K$ fixed to three and the
exact same final partition is obtained.

##  4 Concluding Remarks

Recent research has provided insights into the specification and estimation of
Bayesian finite mixture models including the impact of different priors and
improvements in estimation and inference, in particular regarding
computational issues. In the future we would expect the routine application of
Bayesian finite mixture models in model-based clustering in different areas of
application as well as using data with “unusual” structure. This will induce
the need to include various component-specific distributions or models. We
have focused on mixtures of distributions in this paper. However, clearly also
a regression setting can be included where the components differ with respect
to the regression coefficients or mixture-of-expert models where the component
weights also depend on additional covariates. Recent advances might also help
to use Bayesian finite mixture models for semi-parametric density
approximation.

## Acknowledgments

The authors gratefully acknowledge support from the Austrian Science Fund
(FWF): P28740; and through the WU Projects grant scheme: IA-27001574.

## References

  * Aitkin (2001) Aitkin M (2001).  “Likelihood and Bayesian Analysis of Mixtures.”  _Statistical Modelling_ , 1(4), 287–304.  10.1177/1471082x0100100404.
  * Celeux _et al._ (2019a) Celeux G, Frühwirth-Schnatter S, Robert CP (2019a).  “Model Selection for Mixture Models—Perspectives and Strategies.”  In S Frühwirth-Schnatter, G Celeux, CP Robert (eds.), _Handbook of Mixture Analysis_ , chapter 7, pp. 117–154. CRC Press, Boca Raton, FL.
  * Celeux _et al._ (2000) Celeux G, Hurn M, Robert CP (2000).  “Computational and Inferential Difficulties with Mixture Posterior Distributions.”  _Journal of the American Statistical Association_ , 95(451), 957–970.  10.2307/2669477.
  * Celeux _et al._ (2019b) Celeux G, Kamary K, Malsiner-Walli G, Marin JM, Robert CP (2019b).  “Computational Solutions for Bayesian Inference in Mixture Models.”  In S Frühwirth-Schnatter, G Celeux, CP Robert (eds.), _Handbook of Mixture Analysis_ , chapter 5, pp. 73–96. CRC Press, Boca Raton, FL.
  * Dempster _et al._ (1977) Dempster AP, Laird NM, Rubin DB (1977).  “Maximum Likelihood from Incomplete Data Via the EM-Algorithm.”  _Journal of the Royal Statistical Society B_ , 39, 1–38.  10.1111/j.2517-6161.1977.tb01600.x.
  * Diebolt and Robert (1994) Diebolt J, Robert CP (1994).  “Estimation of Finite Mixture Distributions Through Bayesian Sampling.”  _Journal of the Royal Statistical Society B_ , 56, 363–375.  10.1111/j.2517-6161.1994.tb01985.x.
  * Frühwirth-Schnatter (2001) Frühwirth-Schnatter S (2001).  “Markov Chain Monte Carlo Estimation of Classical and Dynamic Switching and Mixture Models.”  _Journal of the American Statistical Association_ , 96(453), 194–209.  10.1198/016214501750333063.
  * Frühwirth-Schnatter (2006) Frühwirth-Schnatter S (2006).  _Finite Mixture and Markov Switching Models_.  Springer Series in Statistics. Springer-Verlag, New York.  ISBN 0-387-32909-9.
  * Frühwirth-Schnatter (2011) Frühwirth-Schnatter S (2011).  “Dealing with Label Switching Under Model Uncertainty.”  In K Mengersen, CP Robert, D Titterington (eds.), _Mixture Estimation and Applications_ , chapter 10, pp. 213–239. Wiley, Chichester.
  * Frühwirth-Schnatter (2019) Frühwirth-Schnatter S (2019).  “Keeping the Balance—Bridge Sampling for Marginal Likelihood Estimation in Finite Mixture, Mixture of Experts and Markov Mixture Models.”  _Brazilian Journal of Probability and Statistics_ , 33(4), 706–733.  10.1214/19-BJPS446.
  * Frühwirth-Schnatter _et al._ (2019) Frühwirth-Schnatter S, Celeux G, Robert CP (eds.) (2019).  _Handbook of Mixture Analysis_.  Chapman & Hall/CRC Handbooks of Modern Statistical Methods. Chapman and Hall/CRC.
  * Frühwirth-Schnatter and Malsiner-Walli (2019) Frühwirth-Schnatter S, Malsiner-Walli G (2019).  “From Here to Infinity: Sparse Finite Versus Dirichlet Process Mixtures in Model-Based Clustering.”  _Advances in Data Analysis and Classification_ , 13, 33–64.  10.1007/s11634-018-0329-y.
  * Frühwirth-Schnatter _et al._ (2020) Frühwirth-Schnatter S, Malsiner-Walli G, Grün B (2020).  “Generalized Mixtures of Finite Mixtures and Telescoping Sampling.”  ArXiv preprint arXiv: 2005.09918.
  * Grün _et al._ (2021) Grün B, Malsiner-Walli G, Frühwirth-Schnatter S (2021).  “How Many Data Clusters Are in the Galaxy Data Set? Bayesian Cluster Analysis in Action.”  ArXiv preprint arXiv:2101.12686.
  * Hennig (2015) Hennig C (2015).  “What Are the True Clusters?”  _Pattern Recognition Letters_ , 64, 53–62.  10.1016/j.patrec.2015.04.009.
  * Hubert and Arabie (1985) Hubert L, Arabie P (1985).  “Comparing Partitions.”  _Journal of Classification_ , 2(1), 193–218.  10.1007/BF01908075.
  * Jain and Neal (2004) Jain S, Neal RM (2004).  “A Split-Merge Markov Chain Monte Carlo Procedure for the Dirichlet Process Mixture Model.”  _Journal of Computational and Graphical Statistics_ , 13, 158–182.
  * Jain and Neal (2007) Jain S, Neal RM (2007).  “Splitting and Merging Components of a Nonconjugate Dirichlet Process Mixture Model.”  _Bayesian Analysis_ , 3, 445–500.
  * Malsiner-Walli _et al._ (2016) Malsiner-Walli G, Frühwirth-Schnatter S, Grün B (2016).  “Model-Based Clustering Based on Sparse Finite Gaussian Mixtures.”  _Statistics and Computing_ , 26(1), 303–324.  10.1007/s11222-014-9500-2.
  * McCullagh and Yang (2008) McCullagh P, Yang J (2008).  “How Many Clusters?”  _Bayesian Analysis_ , 3(1), 101–120.  10.1214/08-BA304.
  * McLachlan and Peel (2000) McLachlan GJ, Peel D (2000).  _Finite Mixture Models_.  John Wiley & Sons.  10.1002/0471721182.
  * Miller and Harrison (2018) Miller JW, Harrison MT (2018).  “Mixture Models with a Prior on the Number of Components.”  _Journal of the American Statistical Association_ , 113(521), 340–356.  10.1080/01621459.2016.1255636.
  * Molitor _et al._ (2010) Molitor J, Papathomas M, Jerrett M, Richardson S (2010).  “Bayesian Profile Regression with an Application to the National Survey of Children’s Health.”  _Biostatistics_ , 11(3), 484–498.  10.1093/biostatistics/kxq013.
  * Newcomb (1886) Newcomb S (1886).  “A Generalized Theory of the Combination of Observations So As To Obtain the Best Result.”  _American Journal of Mathematics_ , 8(4), 343–366.  10.2307/2369392.
  * O’Hagan _et al._ (2019) O’Hagan A, Murphy TB, Scrucca L, Gormley IC (2019).  “Investigation of Parameter Uncertainty in Clustering Using a Gaussian Mixture Model via Jackknife, Bootstrap and Weighted Likelihood Bootstrap.”  _Computational Statistics_ , 34(4), 1779–1813.  10.1007/s00180-019-00897-9.
  * Papastamoulis (2016) Papastamoulis P (2016).  “label.switching: An R Package for Dealing with the Label Switching Problem in MCMC Outputs.”  _Journal of Statistical Software, Code Snippets_ , 69(1), 1–24.  10.18637/jss.v069.c01.
  * Papastamoulis and Iliopoulos (2009) Papastamoulis P, Iliopoulos G (2009).  “Reversible Jump MCMC in Mixtures of Normal Distributions with Same Component Mean.”  _Computational Statistics & Data Analysis_, 53, 900–911.
  * Pearson (1894) Pearson K (1894).  “Contributions to the Mathematical Theory of Evolution.”  _Philosophical Transactions of the Royal Society A_ , 185, 71–110.  10.1098/rsta.1894.0003.
  * Rastelli and Friel (2018) Rastelli R, Friel N (2018).  “Optimal Bayesian Estimators for Latent Variable Cluster Models.”  _Statistics and Computing_ , 28(6), 1169–1186.  10.1007/s11222-017-9786-y.
  * Reaven and Miller (1979) Reaven GM, Miller RG (1979).  “An Attempt to Define the Nature of Chemical Diabetes Using a Multidimensional Analysis.”  _Diabetologia_ , 16(1), 17–24.  10.1007/BF00423145.
  * Richardson and Green (1997) Richardson S, Green PJ (1997).  “On Bayesian Analysis of Mixtures With an Unknown Number of Components.”  _Journal of the Royal Statistical Society B_ , 59(4), 731–792.  10.1111/1467-9868.00095.
  * Rousseau and Mengersen (2011) Rousseau J, Mengersen K (2011).  “Asymptotic Behaviour of the Posterior Distribution in Overfitted Mixture Models.”  _Journal of the Royal Statistical Society B_ , 73(5), 689–710.  10.1111/j.1467-9868.2011.00781.x.
  * Scrucca _et al._ (2016) Scrucca L, Fop M, Murphy TB, Raftery AE (2016).  “mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models.”  _The R Journal_ , 8(1), 289–317.  10.32614/RJ-2016-021.
  * Stephens (2000) Stephens M (2000).  “Bayesian Analysis of Mixture Models with an Unknown Number of Components – An Alternative to Reversible Jump Methods.”  _The Annals of Statistics_ , 28, 40–74.  10.1214/aos/1016120364.
  * Zhang _et al._ (2004) Zhang Z, Chan KL, Wu Y, Chen C (2004).  “Learning a Multivariate Gaussian Mixture Model with the Reversible Jump MCMC Algorithm.”  _Statistics and Computing_ , 14, 343–355.  10.1023/B:STCO.0000039484.36470.41.