Approximate Bayesian Computation for Finite Mixture Models
Umberto Simola Jessi Cisewski-Kehe Robert L. Wolpert








Approximate Bayesian Computation for Finite Mixture Models
    
Umberto SimolaaCONTACT A. N. Author. Email: umberto.simola@helsinki.fi, Jessi Cisewski–Keheb and Robert L. Wolpertc
aDepartment of Mathematics and Statistics, University of Helsinki, Finland; bDepartment of Statistics, Yale University, New Haven, USA; cDepartment of Statistical Science, Duke University, Durham, USA

    December 30, 2023
==================================================================================================================================================================================================================================================================================================================================




Finite mixture models are used in statistics and other disciplines, but inference for mixture models is challenging due, in part, to the multimodality of the likelihood function and the so-called label switching problem. We propose extensions of the  Approximate Bayesian Computation–Population Monte Carlo (ABC–PMC) algorithm as an alternative framework for inference on finite mixture models. There are several decisions to make when implementing an ABC–PMC algorithm for finite mixture models, including the selection of the kernels used for moving the particles through the iterations, how to address the label switching problem, and the choice of informative summary statistics. Examples are presented to demonstrate the performance of the proposed ABC–PMC algorithm for mixture modeling.  The performance of the proposed method is evaluated in a simulation study and for the popular recessional velocity galaxy data.



Approximate Bayesian computation, Label switching, Finite mixture models, Perturbation kernels, Summary statistics.





§ INTRODUCTION
 

Mixture models have been used in statistics since the late 1800s when Karl Pearson introduced them in an analysis of crab morphometry <cit.>. Subsequently mixture modeling has grown in popularity in the statistical community as a powerful framework for modeling clustered data; McLachlan and Peel <cit.> provides a general overview of mixture modeling, and
Marin et al. <cit.> or Schnatter <cit.> provide a Bayesian perspective. 
In recent decades mixture models have become routinely applied in various applications <cit.>. 
One reason for the general success of mixture models is the ability to specify the number of possibly different component distributions, allowing for flexibility in describing complex systems <cit.>.

The general definition for a finite mixture model with fixed integer K >1 components is:

    ∑_i=1^Kf_i· p_i(y |θ_i),

with mixture weights f_i∈ (0,1) such that ∑_i=1^Kf_i=1 and p_i(y |θ_i) are the component distributions of the mixture, often parametrically specified with a vector of the unknown parameters, θ_i.

Finite mixture models present computational and methodological challenges due, at least in part, to the complex and possibly multimodal likelihood function along with the invariance under permutation of the component indices. The Expectation–Maximization (EM) algorithm <cit.> provides a method for numerically obtaining the maximum likelihood estimates, although the possible multimodality of the likelihood function makes finding the global maximum challenging <cit.>. Extensions of the EM algorithm have been proposed for improving its speed of convergence and avoiding local optima <cit.>. 


Bayesian methodology and applications for mixture modeling have developed in the last decades <cit.>. Bayesian inference for mixture models often relies on Markov Chain Monte–Carlo (MCMC) which can lead to the so-called label switching problem <cit.>, because the likelihood function is invariant to the re-labeling of the mixture components. Additionally the resulting posterior distribution is multimodal and asymmetric, which makes summarizing the posterior distribution using common statistics such as the posterior mean or the Highest Posterior Density (HPD) interval unhelpful <cit.>.


A different framework for inference that can be considered for addressing the issues related to the use of mixture models is Approximate Bayesian Computation (ABC), whose basic algorithm was first presented conceptually by Rubin <cit.> and computationally by Tavaré et al. <cit.> and Pritchard et al. <cit.>. ABC is often used in situations where the likelihood function is complex or not available, but simulation of data through a forward model () is possible. The FM is a simulation model that can take parameters as inputs and returns as output simulated data. With mixture models, though the likelihood function can be available, working with the likelihood is not always straightforward.  
Though ABC has its own challenges, we propose an ABC algorithm for estimating the posterior distributions of the parameters providing an alternative to the MCMC algorithms.


Assuming θ is the possibly multidimensional inferential target, the basic ABC algorithm consists of the following four steps <cit.>: (i) sample from the prior distribution, θ∼π(θ); (ii) produce a generated sample of the data by using θ in the , y∼(Y |θ); (iii) compare the real data, y, with the generated sample, y, using a distance function, ρ(·, ·), letting d = ρ(s(y), s(y)) where s(·) is some (possibly multidimensional) summary statistic of the data; (iv) if the distance, d, is smaller than a fixed tolerance ϵ then θ is retained, otherwise it is discarded.  Repeat until the desired particle sample size N is achieved.

In order to improve the computational efficiency of the basic ABC algorithm, there have been many developments  <cit.>. In this work, we extend a sequential version of ABC, the ABC–PMC of Beaumont et al. <cit.>, which is introduced in the next Section. Interestingly, the original PMC sampler <cit.> was introduced for solving problems related to the multimodality of the likelihood function for settings in which mixture models were used.  This suggests that the ABC–PMC version of the ABC algorithm may be more suitable for mixture models compared to other available algorithms, such as the Sequential Monte Carlo version (ABC–SMC) by <cit.>. ABC may be considered preferable over MCMC if the likelihood function is intractable. However, in the finite mixture model setting, the likelihood function is tractable, but there are computational challenges due to the multimodality of the likelihood function as well as the label switching problem that can lead to issues for MCMC <cit.>.

In any ABC analysis the definition of both suitable summary statistics, s(·), and a distance function, ρ,  for comparing the real data  y to the generated sample y is needed and is crucial for getting useful inference results <cit.>. The definition of summary statistics is necessary because ABC suffers of the curse of dimensionality and using the entire dataset is not computationally feasible <cit.>. 
In the proposed method, we generally use the Hellinger distance for comparing the real data, y, to the generated sample, y. 

While the proposed method is valid for y∈ℝ^dim, for illustration purposes, we focus most of the discussion on the one dimensional case where dim=1.

The paper is organized as follow. Section <ref> begins with background on the ABC–PMC algorithm by Beaumont et al. <cit.>, and then details our proposed extensions for the finite mixture model setting.
Sections <ref> and <ref> are dedicated, respectively, to a simulation study and a popular real data example using recessional velocities of galaxy data. 


Concluding remarks are presented in Section <ref>.


 

§ METHODOLOGY
 
 
Sequential versions of the basic ABC algorithm have been proposed to improve the computational efficiency of the algorithm. This is typically accomplished by considering a sequence of decreasing tolerances, ϵ_1:T= (ϵ_1, …, ϵ_T) such that ϵ_1 ≥ϵ_2 ≥⋯≥ϵ_T, where T is the total number of iterations, or time steps, of the algorithm. 
Rather than sampling from the prior at each iteration after the initial step, improved proposals are drawn from the ABC posterior sample of the previous step. 
The  ABC–PMC algorithm of Beaumont et al. <cit.> relies on ideas of importance sampling  in order to improve the efficiency of the algorithm by constructing a series of intermediate proposal distributions using the sequential ABC posteriors. 
The first iteration of the ABC–PMC algorithm uses tolerance ϵ_1 and draws proposals from the specified prior distributions; from the second iteration through the T^th, proposals are selected from the previous iteration's ABC posterior and then moved according to some perturbation kernel (e.g., a Gaussian kernel).  
Since the proposals are drawn from the previous iteration's ABC posterior and moved according to a kernel, rather than being sampled directly from the prior distributions, importance weights are used. The derivation of the importance weight for particle J = 1, …, N in iteration t, where N is the desired particle sample size, is discussed in Beaumont et al. <cit.> and leads to the following importance weights:

    W_t^(J)∝π(θ_t^(J))/∑_L = 1^N W_t-1^(L)ϕ(θ_t^(J)|θ_t-1^(L), τ_θ,t-1^2),


where π(·) is the prior distribution, ϕ(·| a, b) is a Gaussian perturbation kernel with mean a and variance b, where  τ_θ,t-1^2 is twice the (weighted) sample variance of the particles for θ at iteration t-1 (as recommended <cit.>).  
The importance weights are then standardized so that ∑_J = 1^N W_t^(J) = 1.
Other perturbation kernels can be used and the importance weights would need to be updated accordingly. 
Further details on the original ABC–PMC algorithm can be found in Beaumont et al. <cit.>, though there are
several elements that require additional discussion.  In particular, an implementation of ABC–PMC requires the user to select the tolerance sequence ϵ_1:T, the perturbation kernels, the summary statistic(s), and the distance function(s).

First, to select the sequence ϵ_1:T we implement an adaptive approach based on the accepted distances from the previous time step. That is, ϵ_t is set to a particular quantile of the accepted distances from time step t-1 <cit.>.
Selecting too high of a quantile at which to shrink the tolerance could result in a decrease in computational efficiency because more iterations would be needed to shrink the tolerance to a small enough value.
On the other hand, selecting a quantile that is too low can also contribute to computational inefficiencies because more draws from the proposal are needed within each iteration to find datasets that achieve the small tolerance. 
The rule used to shrink the tolerance is important as imprudent selection can lead to the ABC–PMC algorithm getting stuck in local modes <cit.>. 
This latter potential issue is particularly crucial when working with mixture models, since the multimodality of the likelihood function is a well known problem that needs to be addressed in order to produce reliable statistical inference results.   


The decreasing tolerance sequence's impact on the achievement of the global maximum is addressed in two ways in the proposed algorithm.
First, we initialize the tolerance sequence by oversampling in the first iteration of the ABC–PMC algorithm <cit.>. 
Let N be the desired number of particles used to approximate the posterior distribution, then the initial tolerance, ϵ_1, can be adaptively selected by sampling N = lN draws from the prior, for some l ∈ℤ^+.  
Then the N particles of the lN total particles with the smallest distances are retained, and ϵ_1 = max(d_1^(1*), …, d_1^(N*)) where d_1^(1*), …, d_1^(N*) are the N smallest distances of the lN particles sampled.
For all examples considered in this work, we use l=5. 
Using this adaptive initialization, the parameter space is better sampled than if we only considered N draws from the prior distribution. 

Secondly, we use lower quantiles for the first several iterations of the ABC–PMC algorithm per the suggestion in Silk et al. <cit.> and Simola et al. <cit.>, which can help the algorithm avoid local modes.



One of the challenges in using ABC–PMC for mixture models is related to selecting an appropriate perturbation kernel for the mixture weights because the individual weights must be between 0 and 1, and the weights must sum to 1.  Hence, the usual Gaussian perturbation kernel is not a viable option because this kernel can lead to proposed mixture weights that are not within the noted constraints. 
An additional challenge in using ABC–PMC for mixture models is due to the label switching problem, and for reasons discussed below, this has to be addressed at the end of each iteration.  Finally, the selection of an appropriate and informative summary statistic for multimodal data requires a summary that can capture the overall shape of the distribution of the data.
 





For these reasons the original version of the ABC–PMC algorithm is not directly able to accommodate mixture models. Our proposed extensions are discussed in more detail later, but can be summarized as follows. In Algorithm <ref>, we adapt the ABC–PMC algorithm to work with finite Gaussian mixture models. In particular, we propose a procedure for resampling the mixture weights in a way that capitalizes on distributional properties of the prior distribution for the mixture weights (i.e., the Dirichlet distribution), which is presented in Algorithm <ref>. The mixture weights could be resampled using other methods, such as resampling the particles according to a Gaussian perturbation kernel, but the resampled weights would need to be adjusted in order to preserve properties of weights such as being between or equal to 0 and 1 and summing to 1. Additionally, in order to address the label switching problem, we designed a procedure that is executed at the end of each iteration of the ABC–PMC algorithm, which is presented in Algorithm <ref>. An overview of different available algorithms designed for addressing the label switching problem and more details on Algorithm <ref> are presented in Section <ref>.  If not properly addressed, the label switching problem can hinder the algorithm's capability of finding the true posterior distribution.








 §.§ Finite Gaussian Mixture Models

A common choice for the p_i(·|·) of Eq. (<ref>) is the Gaussian distribution. This particular class of models, called Gaussian Mixture Models (GMMs), is  flexible and widely used in various applications (e.g., <cit.>). Maintaining the notation of Eq. (<ref>), a GMM is defined as:


    ∑_i=1^Kf_i·ϕ(y |θ_i),

where ϕ is the density function of the Gaussian distribution, θ_i = (μ_i, σ_i^2) is the vector of unknown means and variances for each of the K groups with the unknown parameters  μ = (μ_1, …, μ_K) ∈ℝ^K, σ^2 = (σ_1^2, …, σ_K^2) ∈ℝ_+^K, and the mixture weights, f =(f_1, …, f_K). Note that the f =(f_1, …, f_K) lie in the unit simplex, Δ^K-1≡{x∈ℝ^K_+: ∑_j x_j=1} inside the unit cube [0,1]^K. 

A common choice of the prior distribution for f=(f_1, …, f_K) is the Dirichlet distribution of order K with hyperparameter δ=(δ_1,…δ_K), where often δ≡ (1, …, 1), as proposed by Wasserman <cit.>. Another common choice has been proposed by Rousseau and Mengersen <cit.>, who defined the hyperparameter δ≡ (1/2, …, 1/2); in this way the prior is marginally a Jeffreys prior distribution. The priors for the mean and the variance of the GMM from Eq. (<ref>) can be defined as follows:


    μ_i∼ϕ(μ|ξ,κ),     
    σ_i^2∼Inv–Gamma(α,β),

with mean ξ, variance κ, and shape parameter α and rate parameter β. 
There are several methods for selecting the hyperparameters, η = (ξ, κ, α, β), such as the Empirical Bayes approach <cit.> and the `weakly informative principle' <cit.>. Both of these options are considered in examples later so as to be consistent with the original analysis of the example.

Our proposed algorithm is specified in Algorithm <ref> for the model defined in Eq. (<ref>), and the details of important steps are presented next.







 §.§ Perturbation kernel functions
 
One of the advantages of ABC–PMC over the basic ABC algorithm is that, starting from the second iteration, rather than drawing proposals from the prior distribution, proposed particles are drawn from the previous step's ABC posterior according to their importance weights. 
Then, instead of using the actual proposed value that was drawn, it is perturbed according to some kernel. Constructing computationally efficient perturbation kernels is crucial to improve the computational performance of the ABC algorithm  <cit.>.
There are a number of kernel functions, K(·|·), available for perturbing the proposed particles. Beaumont et al. <cit.> suggest a Gaussian kernel centered on the selected particle from the previous iteration and a variance equal to twice the empirical variance of the previous iteration's ABC posterior. This is a reasonable choice if there are no constraints on the support of the parameter of interest. 
For the means and the variances of the mixture model, a multivariate Gaussian distribution and a multivariate Gaussian distribution truncated to ℝ_+  were used, respectively. Both perturbation kernels are assigned variances per the suggestion from Beaumont et al. <cit.> (i.e., variances set to be twice the empirical variances of the previous iteration's ABC posteriors).
When constraints on the parameter support are present, such as for variances or mixture weights, a perturbation kernel should be selected so that it does not propose values outside the parameter's support. 


When moving the selected values for the mixture weights, not only is there the constraint that each mixture weight component must be in [0,1], but it is also required that ∑_i = 1^K f_i = 1, making the Gaussian kernel inappropriate.


In the first iteration of the proposed ABC–PMC algorithm and for J = 1, …, N, the mixture weights f^(J)_1 =  {f_1,1, …, f_1,K}^(J)
are directly sampled from the prior distribution, which is a Dirichlet(δ), where δ = (δ_1, …, δ_K). 
For t>1 and for J = 1, …, N, proposals are drawn from the previous iterations's particle system according to their importance weights, as defined in Algorithm <ref>. 
After randomly selecting some particle L ∈{1, …, N} from the t-1 iteration's particle system we want to “jitter” or move the selected particle's mixture weight, 
f_t^*, in manner that preserves some information coming from the selected particle, but not let it be an identical copy, in order to obtain the resampled mixture weights, 
f_t^(L). This is carried out using Algorithm <ref> and the mathematical justification can be found in Appendix <ref>.


In Algorithm <ref>, for each resampled mixture weight for iteration t = 2, …, T and mixture component i = 1, …, K, rely on the random variable denoted ξ^*_t,i (the particle index J is suppressed for notational convenience in what follows), which is designed to be the sum of two independent random variables, Z_t f_t^*B_t,i∼Gamma(pδ_i,1) and η_t,i∼Gamma((1-p)δ_i,1).
The resampled mixture weight is f_t =  ξ_t,i^*/ξ_t,+^*, with ξ_t,+^*=∑_i=1^Kξ_t,i^*,
 so that f_t ∼Dirichlet(δ).  
The key benefits of this procedure are that the resampled weight vector, f_t, retains some information from the f_t^* (i.e., the previous iterations particle that was selected for resampling) through the parameter p (discussed next), but also follows the same distribution as  the prior, Dirichlet(δ), so that the importance weights due to this resampling are simply unity.

The parameter p is a fixed real number with range [0,1] that determines how much information to retain from f_t^*. The choice of p has an impact on both the allowed variability of the marginal ABC posterior distributions for the mixture weights and the efficiency of the entire procedure. In particular, fixing a p close to 1 leads to a Dirichlet resampling in which the new set of mixture weights f_t is close to the previous set f_t^* (if p=1, then f_t= f_t^*). On the other hand a choice for p close to 0 implies that the information coming from f_t^* is weakly considered (for p=0 a new set of particles is drawn directly from the prior distribution and no information about f_t^* has been retained). We found p=0.5 to be a good choice for balancing efficiency and variability (i.e., it allows for some retention of information of the selected particle, but does not lead to nearly identical copies of it). 










 §.§ Algorithm for addressing the label switching problem
 

As noted earlier, a common problem arising in the Bayesian mixture models framework is label switching. When drawing a sample from a posterior (for both MCMC and ABC), the sampled values are not necessarily ordered according to their mixture component assignments because the likelihood is exchangeable. 



For example, suppose a particle 
{(f_1^(J), …, f_K^(J)), (μ_1^(J), …, μ_K^(J)), (σ^2^(J)_1, …, σ^2^(J)_K)},  where J ∈{1, …, N} , is accepted for a K component GMM.  (Note that in this Section, we do not include the algorithm step index, t, to simply the notation.)  This particle was selected with a particular ordering of the 1, …, K components with f_i^(J), μ_i^(J), and σ^2^(J)_i coming from the same mixture component, i = 1, …, K; however, a new particle that is accepted will not necessarily follow that same ordering of the i = 1, …, K components.
Somehow the particles have to be ordered in such a way that aligns different realizations of the i = 1, …, K components in order to eliminate the ambiguity. 



Several approaches have been proposed to address the label switching problem and are known as relabeling algorithms.  A first group of relabeling algorithms consists of imposing an artificial identifiability constraint in order to arbitrarily pick a parameter (e.g., the mixture weights) and sort all the parameters for each accepted particle according to that parameter's order <cit.>. However, many of the algorithms proposed for addressing the label switching are deterministic (e.g., Stephen's method <cit.> and the pivotal reordering algorithm <cit.>). A third class of strategies, called probabilistic relabeling algorithms, uses a probabilistic approach for addressing the label switching problem <cit.>. A detailed overview of methods for addressing label switching is presented in Papastamoulis <cit.>.
In Section <ref>, we provide an example that illustrates a problem with the artificial identifiability constraint approach.
Instead, we propose a deterministic strategy for addressing the label switching by selecting a parameter that has at least two well-separated components.  



Addressing the label switching problem is especially important for the proposed sequential ABC algorithm because each time step's ABC posterior is used as the proposal in the subsequent step of the algorithm so the label switching has to be resolved before using it as a proposal distribution.
Algorithm <ref> outlines the proposed strategy, and is carried out at the end of each iteration. 
The key aspect of Algorithm <ref> is to select a parameter for reordering that has at least two well-separated components. 
To determine this, each set of parameters (e.g., the means, the variances, the mixture weights), is arranged in increasing order. For example, for each particle J, J=1,…,N, μ^(J) would be ordered so that μ_(1)^(J)≤μ_(2)^(J)≤⋯≤μ_(K)^(J), with μ_(i)^(J) as the i^th order statistic. This is carried out for each set of parameters with analogous notation.


The next step is to determine which set of parameters has the best separated components values. We propose first shifting and scaling each set of parameters to be supported within the range [0,1] so that scaling issues are mitigated and the different parameter set values are comparable. 

One option for this adjustment is to use some distribution function, such as a Gaussian distribution with a mean and variance equal to the sample mean and the sample variance of the considered parameter set (e.g., the sample mean for the μ's is μ̅=(NK)^-1∑_k=1^K∑_J=1^Nμ_(k)^(J), and the sample variance for the μ's is var(μ)=(NK)^-1∑_k=1^K∑_J=1^N (μ_(k)^(J)-μ̅)^2). 

The resulting k^th largest standardized value for the mixture mean is μ^(J)_(k). This is carried out for each set of parameters with analogous notation. 


For each component of each ordered and standardized particle a representative value, such as a mean, is computed (e.g., μ̅_(k) = N^-1∑_J = 1^N μ^(J)_(k) is the representative value of the k^th component of the mean parameter). This is carried out for each set of parameters with analogous notation. 
The pairwise distances (pdist) between the representative values within each parameter set is determined.  
The parameter that has the largest separation between any two of its representative values is selected for the overall ordering of the particle system.



Variations of Algorithm <ref> were considered.  
For example, rather than sorting based on the parameter set with largest separation between any two of its representative values, we considered basing the sorting on the parameter with the largest separation between its two closest representative values (i.e., the maximum of the minimum separations); however, this sorting did not perform well empirically.  
The issue seemed to be that  parameter with the largest separation between its two closest representative values may actually have all of its components relatively close; after multiple iterations, none of the components were well separated  from the other components.  This lead to iteration after iteration of components that remained a blend of components rather than separating clearly into K components.
Overall, from empirical experiments, the approach outlined in Algorithm <ref> performed the best and thus is our recommendation.  However, we emphasize that alterations to this procedure may be necessary for mixture models that have additional structure or correlations among the parameters of the component distributions or to deal with multi–dimensional mixture models.



 §.§ Summary statistics
 

Comparing the entire simulated dataset, y, with the entire set of observations, y, in an ABC procedure is not computationally feasible. For this reason the definition of a lower dimensional summary statistic is necessary. 
For mixture models, due to the multimodality of the data, common summaries such as means or higher order moments (e.g., variance, skewness or kurtosis) do not capture relevant aspects of the distribution which can lead to incorrect ABC posterior distributions. However, an estimate of the density of the data can better account for its key features (e.g., the shape of
each component of the mixture). For the examples presented in Sections <ref> and <ref>, we suggest using kernel density estimates (KDEs, <cit.>) of the simulated data, f̂_y_prop, n, and the observations, f̂_y_obs, n, to summarize the data, and then the Hellinger distance, H, to carry out the comparison. We considered other distances, such as the Kolmogorov–Smirnov (KS) distance <cit.>, but found the Hellinger distance to perform better than any other tested distance function. The Hellinger distance quantifies the similarity between two density functions, f and g, and is defined as:

    H(f,g) = (∫(√(f(y))-√(g(y)))^2 dy )^1/2.

At each iteration t of the proposed ABC–PMC procedure, a proposed θ is accepted if H(f̂_y_obs, n,f̂_y_prop, n) < ϵ_t, where ϵ_t is the tolerance.












§ SIMULATION STUDY
 

In this Section a simulation study is presented to evaluate the behavior of the proposed ABC–PMC algorithm, Algorithm <ref>, using popular GMM examples. 
In particular we are interested in evaluating the success of the procedure with respect to the label switching problem and the reliability of the Hellinger distance as a summary statistic. 
To determine the number of iterations, T, a stopping rule was defined based on the Hellinger distance between sequential ABC posteriors; once the sequential Hellinger distance dropped below a threshold of 0.05 for each of the marginal ABC posteriors, the algorithm was stopped.



 §.§ Mixture Model with equal and unequal group sizes
 

The first example is taken from Mena and Walker <cit.>, which considered a GMM obtained by simulating data coming from K=2 groups of equal size which was designed for evaluating performance with respect to the label switching problem.
The 40 observations were simulated as follows:

    Y_i=1,…,20∼ϕ(y | -20,1),   Y_i=21,…,40∼ϕ(y | 20,1),

where ϕ is the Normal density function as defined previously. 

The variance is assumed known for both groups and hence the parameters of interest are the mixture weights, f = (f_1, f_2) (with f_2 = 1 - f_1), and the means μ = (μ_1, μ_2). 
The prior distributions selected for this analysis are



μ_i∼ϕ(μ| 0,100), 









and
f=(f_1, f_2) ∼Dirichlet(1,1).


The number of particles is set to N=5000 and the quantile used for shrinking the tolerance is q=0.5. The algorithm was stopped after T=20 iterations, since further reduction of the tolerance did not lead to a notable improvement of the ABC posterior distributions (evaluated by calculating the Hellinger distance between the sequential ABC posterior distributions as noted previously).  

Figure <ref> displays the resulting ABC posteriors (`ABC Posterior Alg. 3 LS' and `ABC Posterior naive LS' are discussed later for illustrating the label switching issue) 
and the corresponding MCMC posteriors, which are used as a benchmark to access the performance of the proposed method. The ABC posteriors closely match the MCMC posteriors. The summary of the results presented in Table <ref> demonstrate that the ABC posterior distributions are a suitable approximation of the MCMC posteriors. 
The Hellinger distances between the marginal ABC and MCMC posteriors are also displayed in the last column of Table <ref>; the Hellinger distance between the MCMC and the ABC posterior is 0.032 for the mixture weights and 0.21 for the means. 


































As noted in Section <ref>, the label switching problem has to be carefully addressed when using the ABC–PMC algorithm. For each time step following the initial step, the previous step's ABC posterior is used as the proposal rather than the prior distribution so the procedure for addressing the label switching proposed in Section <ref>, Algorithm <ref>, is used at the end of each time step.
In order to illustrate the consequences of incorrectly addressing the label switching, we ran the proposed ABC algorithm on the example proposed by Mena and Walker <cit.>, except rather than using Algorithm <ref>, the ordering of the particle system is carried out using the ordering of the mixture weights; the mixture weights are both equal to 0.5 for this example, making them a poor choice for attempting to separate out the mixture components.  
The resulting ABC posteriors are displayed in Figure <ref> (blue dashed lines).  The means, μ = (μ_1, μ_2), of the mixture components are shuffled and not close to the MCMC posterior, while the mixture weights are sorted such a way all the elements of f_1 are smaller than 0.5 and all the elements of f_2 are larger than 0.5. 

Mena and Walker <cit.> added a third component to the mixture in Eq. (<ref>), simulating five additional observations from a standard Gaussian distribution, obtaining a three-component GMM with known variances. The ABC–PMC algorithm was run with the same specifications as the first part of the example, but required T=25 steps to achieve our stopping rule.

Figure <ref> shows the MCMC and the ABC posteriors for the weights and the means of the mixture components. 
The behavior of the ABC posterior distributions is consistent with their MCMC benchmarks. 
A summary of the results presented in Table <ref> shows that the posterior means and the posterior standard deviations for the ABC posterior distributions are consistent with the ones retrieved using MCMC. Finally, in the third column, the Hellinger distances between the ABC and MCMC posteriors are provided. 





































 §.§ Mixture Model with unequal group size
 


Even in those cases in which the definition of the mixture model does not lead to the label switching problem, a second category of issues related to the multimodality of the likelihood function is present. This behavior has been studied from both the frequentist and the Bayesian perspectives. In particular, Marin et al. <cit.> defined the following simple two-component mixture model to illustrate the multimodality issue,

    f_1 ·ϕ(μ|μ_1,1)+(1-f_1) ·ϕ(μ|μ_2,1),

where the weight f_1 is assumed known and different from 0.5 (avoiding the label switching problem). 

According to the specifications by  Marin et al. <cit.>, n=500 samples were drawn from Eq. (<ref>), with θ=(f_1,μ_1,μ_2)=(0.7,0,2.5). 
The bimodality of the likelihood function (see Figure <ref>) makes the use of both the EM algorithm <cit.> and the Gibbs Sampler <cit.> risky because their success depends on the set of initial values selected for initiating the algorithms.

The PMC sampler <cit.> is used as a benchmark for the proposed ABC–PMC solution. 
Figure <ref> displays the log likelihood function (note the two modes), and the final ABC–PMC posteriors with MCMC posteriors using good and bad starting values along with the posteriors using the PMC algorithm.


Table <ref> lists the means for the final ABC, MCMC, and PMC posteriors for μ_1 and μ_2, along with the Hellinger distance between the final ABC–PMC posteriors and the PMC posterior.  The ABC–PMC posteriors closely match the PMC posteriors.
















































§ APPLICATION TO GALAXY DATA
 
The galaxy dataset was introduced to the statistical community in Roeder <cit.>, and has been commonly used for testing clustering methods. The data contains the recessional velocities of 82 galaxies (km/sec) from six well separated sections of the Corona Borealis region. The galaxy data has been well-studied by the statistical community <cit.>.  
The recessional velocities of the galaxies are typically considered realizations of independent and identically distributed Gaussian random variables, but there is discrepancy in the conclusions about the number of groups in the GMM; estimates vary from three components <cit.> to five or six <cit.>.

In this analysis, we focused on the model with three components <cit.>, in order to be consistent with <cit.> and <cit.>. 
For the hyperparameters, we use the Empirical Bayes approach suggested by Casella et al. <cit.>. 
Additionally, since each recessional velocity was also assigned a measurement error, the ABC forward model has been modified to take into account this information. 
In order to include the measurement errors in the forward model, each simulated recessional velocity is assigned one of the observed measurement errors. The simulated and observed recessional velocities were matched according to their ranks, and the measurement error of the observations were assigned to the simulated data according to this matching. Then, Gaussian noise was added to each simulated recessional velocity with a standard deviation equal to its assigned measurement error.

The posterior means for each component's parameters are listed in Table <ref>.  
The third component was found to have a weight equal to 0.057, and a mean and variance equal to 32.94 km/sec and 1.16 km^2/sec^2, respectively. 
The main difference with the proposed ABC–PMC estimates and those from Marin et al. <cit.>, in which they also fixed K = 3, is with the estimated mean and variance of the third component; however, the proposed ABC–PMC estimates are consistent with those obtained by Mena and Walker <cit.>. 


By using the additional information about the measurement errors, the proposed ABC–PMC algorithm can provide posterior distributions that better reflect the uncertainty in the observations.
Including measurement errors in the forward model affects the resulting ABC posterior as reported in Table <ref> and the relative estimates plotted in orange in Figure <ref>. In particular, the estimated variances are smaller for the proposed ABC–PMC algorithm when measurement errors are accounted for, which results in more clearly defined clusters for the first and third components of the mixture model.















§ CONCLUDING REMARKS
 

The popularity of ABC is, at least in part, due to its capacity to handle complex models.  Extensions of the basic algorithm have lead to improved efficiency of the sampling, such as 
the ABC–PMC algorithm of Beaumont et al. <cit.>.  We propose an ABC–PMC algorithm that can successfully handle finite mixture models.
Some of the challenges with inference for finite mixture models are due to the complexity of the likelihood function including its possible multimodality and the exchangeability of the mixture component labels leading to the label switching problem.
Fortunately, ABC can handle complicated likelihood functions, but the label switching problem needs to be addressed.  
We suggest a procedure for addressing the label switching problem within the proposed ABC algorithm that works well empirically.

Some additional challenges with using ABC for mixture models include the selection of informative summary statistics, and defining a kernel for moving the mixture weights since they are constrained to be between 0 and 1 and must sum to 1.

For the summary statistics, we propose using the Hellinger distance between kernel density estimates of the real and simulated observations, which allows the multimodality of the data to be accounted for and compared between the two sets of data.
We propose a Dirichlet resampling algorithm for moving the mixture component weights that preserves some information from the sampled particle, but also improves the efficiency of the ABC–PMC procedure (by not having to draw from the Dirichlet prior at each time step). 


The proposed ABC algorithm has been explored and tested empirically using popular examples from the literature.  The resulting ABC posteriors were compared to the corresponding MCMC posteriors, and in all cases considered the proposed ABC and MCMC posteriors were similar, as desired.  

We also considered the galaxy velocity data from the Corona Borealis Region <cit.>, which is often used in assessing the performance of mixture models.  
An advantage of ABC over other commonly used methods is that the forward model can be easily expanded to better represent the physical process that is being modeled.  For the galaxy data, measurement errors are available in the original data set, but are not generally used when analyzing the data.  We extend the proposed ABC–PMC forward model used in this example to include the measurement errors, which provides a more accurate assessment of the uncertainty in the data.

Though the presented examples focused on one–dimensional GMMs, the component distributions can be adapted to other distributions in the ABC forward model.  Overall, the proposed ABC–PMC algorithm performs well and is able to recover the benchmark MCMC posteriors (when they were available) suggesting that the proposed ABC–PMC algorithm is a viable approach for carrying out inference for FMMs. 









§ MIXTURE WEIGHT RESAMPLING KERNEL
 

The steps for resampling the mixture weights at iteration t, {f_t^(J)}_J=1^N, are outlined in Algorithm <ref> with details presented in Section <ref>.  
In this Section, we present mathematical justification for the proposed resampling procedure.
In particular, we want to show that the resampled weights, f_t = (f_t,1, …, f_t,K) with f_t,i = ξ_t,i^*/ξ_t,+^*, i = 1, …, K, ξ_t,+^*=∑_i = 1^Kξ_t,i^*, and t = 2, …, T, follows a Dirichlet(δ_1, …, δ_K).

Recall that to resample a particle at iteration t, a particle from the previous iteration, t-1, is randomly sampled according to the importance weights, {W_t-1^(J)/∑_L = 1^NW_t-1^(L)}_J=1^N.
Suppose some mixture weight vector f_t-1^(J*) = f_t-1 = {f_t-1,1, …, f_t-1,K}, t ∈{2, …, T} has been selected for resampling with  some finite number of components K>1.  To show that carrying out the resampling steps as outlined in Algorithm <ref> results in f_t∼Dirichlet(δ_1, …, δ_K), consider the following propositions.
Proposition <ref> and <ref> correspond to Step 5 of Algorithm <ref>, and Proposition <ref> and <ref> correspond to Steps 1 and 3, respectively.

 
Let f_t-1 = (ξ_t-1,1/ξ_t-1,+, …, ξ_t-1,K/ξ_t-1,+) where ξ_t-1,iGamma(δ_i, β) and ξ_t-1,+ = ∑_i = 1^K ξ_t-1,i, then f_t-1∼Dirichlet(δ).



This is a well-known result that can be proved by using the change-of-variables approach with y_1 = ξ_t-1,1/y_K, …, y_K-1=ξ_t-1,K-1/y_K, and y_K = ξ_t-1,+.  The determinant of the Jacobian of the transformation is y_K^K-1.  Noting the joint distribution ξ_t-1,1, …, ξ_t-1,K∼  e^-∑_i=1^K βξ_t-1,i∏_i = 1^K β
^δ_iξ_t-1,i^δ_i-1/Γ(δ_i), then the joint distribution using the change-of-variables is

    y_1, …, y_K ∼β^∑_i=1^Kδ_i e^- β y_K(y_K^K-1) (y_K (1 - ∑_i = 1^K-1 y_i))^δ_K-1∏_i = 1^K-1(y_iy_K)^δ_i-1/∏_i = 1^KΓ(δ_i).


Then, integrating out the y_K results in 

    y_1, …, y_K-1∼(1 - ∑_i = 1^K-1 y_i)^δ_K-1Γ(∑_i=1^Kδ_i)/∏_i = 1^KΓ(δ_i)∏_i = 1^K-1y_i^δ_i-1.


Noting that 0≤ y_1, …, y_K-1≤ 1 and 0 ≤∑_i=1^K-1 y_i ≤ 1, then we can let y_K = 1 - ∑_i = 1^K-1 y_i, resulting in

    y_1, …, y_K ∼Γ(∑_i=1^Kδ_i)/∏_i = 1^KΓ(δ_i)∏_i = 1^K y_i^δ_i-1,

which is a Dirichlet(δ_1, …, δ_K) distribution.




 
Let f_t-1∼Dirichlet(δ) with δ = (δ_1, …, δ_K) and Z_t ∼Gamma(δ_+, β) where δ_+ = ∑_i = 1^K δ_i.  Then ξ_t,i = Z_tf_t-1,i∼Gamma(δ_i, β).



This follows from Theorem 1 of Yeo and Milne <cit.> by noting that the marginal distribution of a Dirichlet(δ), that is, the distribution of f_t-1,i, is a Beta(δ_i, ∑_i = 1^K δ_j - δ_i).





 
Given p ∈  [0,1] and independent B_t,i∼Beta(pδ_i, (1-p)δ_i) for i = 1, …, K, then ξ_t,iB_t,i∼Gamma(pδ_i, β).



This follows directly from Theorem 1 of Yeo and Milne <cit.>.



 
Given p ∈  [0,1],  independent B_t,i∼Beta(pδ_i, (1-p)δ_i), independent η_t,iGamma((1-p)δ_i, β), i = 1, …, K, 
and set ξ^*_t,i = Z_tf_t-1,iB_t,i + η_t,i and f_t,i = ξ_t,i^*/ξ_t,+^*, with ξ_t,+^*=∑_i=1^Kξ_t,i^*.  Then f_t = (f_t,1, …, f_t,K) ∼Dirichlet(δ_1, …, δ_K).



By Propositions <ref> and <ref>, Z_tf_t-1,iB_t,i∼Gamma(pδ_i, β), and then ξ^*_t,i∼Gamma(δ_i, β) by noting that the sum of the two independent gamma random variables with the same rate parameter, β, is also a gamma random variable with rate parameter β and with the sum δ_i = pδ_1+(1-p)δ_1 as the shape parameter.  Then Proposition <ref> gives the distribution for f_t.


unsrt


