An effective strategy for initializing the EM algorithm in finite mixture models

Abstract

Finite mixture models represent one of the most popular tools for modeling heterogeneous data. The traditional approach for parameter estimation is based on maximizing the likelihood function. Direct optimization is often troublesome due to the complex likelihood structure. The expectation–maximization algorithm proves to be an effective remedy that alleviates this issue. The solution obtained by this procedure is entirely driven by the choice of starting parameter values. This highlights the importance of an effective initialization strategy. Despite efforts undertaken in this area, there is no uniform winner found and practitioners tend to ignore the issue, often finding misleading or erroneous results. In this paper, we propose a simple yet effective tool for initializing the expectation–maximization algorithm in the mixture modeling setting. The idea is based on model averaging and proves to be efficient in detecting correct solutions even in those cases when competitors perform poorly. The utility of the proposed methodology is shown through comprehensive simulation study and applied to a well-known classification dataset with good results.

Appendix

1.1 Assessment of the choice of weights

Fig. 10

Results of the simulation study comparing weights to be used in emaEM algorithm. Gaussian mixtures are simulated from different (K, p). The y-axis is the mean and standard deviation of the adjusted Rand index calculated from 100 replication of each case. Solid line is for equal weights and dashed line is for weights based on BIC

As shown in Sect. 2, we need to assign weights to different models when aggregating. The weights can either be based on BIC or we can assign equal weights to all models. In this section, we conducted a small simulation study to see which one gives better performance. We considered a Gaussian mixtures with $K = 3,15,30$ components and $p = 2,6$ dimensions. All mixtures are simulated using MixSim package. From each mixture we simulated 100 datasets of size $n = 50 \times K$. We run the emaEM initialization algorithm on each dataset with the two approaches to weights (BIC and equal). In all experiments, short EM is run 100 times. Although the weights based on BIC have shown a good performance in a variety of settings (Hoeting et al. 1999), in this case, assigning equal weights to all models gives better output. This occurs because, almost all the weight is assigned to the model with the lowest BIC and all other models are assigned almost zero weight when BIC-based weights are used. Therefore, this is essentially running the emEM algorithm. On the other hand, assigning equal weights to all models, the information from all models is gathered. Therefore, we will get some information for each pair of points from all models. Figure 10 shows results for each mixture combination considered. We report the mean and standard deviation of $\mathcal A\mathcal R$ for equal weights as well as those for weights based on BIC. As we can see, there is a better performance observed when we use equal weights. Hence, for all experimental and real data analysis we used equal weights.

1.2 Performance assessment for unequal cluster sizes

Fig. 11

Results of the simulation study for unequal cluster sizes. The x-axis represents the competing initialization methods and y-axis represents the mean adjusted Rand index value achieved by each method. Each sub-plot represents different ($K,p,\check{\omega },\tau _{(1)}$) combinations

In this section, we assess the performance of the emaEM initialization method when groups have unequal sizes. This interesting point was raised by one of the anonymous reviewers. That is, whether or not the new initialization methodology (emaEM) will be affected by unequal cluster sizes. To check this, we simulated data from a Gaussian mixture with the smallest cluster having certain fixed size. This can be specified by using the PiLow option of the function MixSim in the MixSim package. Similarly to previous simulation studies presented, Gaussian mixtures varied in the number of components ($K \in \{3,15,30\}$), dimensions ($p\in \{2,6\}$), and amount of maximum overlap between components ($\check{\omega }\in \{0, 0.05, 0.1\}$). The smallest cluster is specified to have a certain proportion of the total observations. To be consistent with our notation, we will use $\tau _{(1)}$ to represent the mixing proportion of the smallest cluster. In our simulations, we used different values of $\tau _{(1)}$ to see if there will be a clear trend in performance as we increase the size of the smallest cluster. The values considered are $\tau _{(1)} \in \{0.1,0.15, 0.2\}$ for mixtures with $K = 3$, $\tau _{(1)} \in \{0.02,0.03, 0.04\}$ for mixtures with $K = 15$, and $\tau _{(1)} \in \{0.01,0.015, 0.02\}$ for mixtures with $K = 30$. Note that the clusters with similar sizes would have $\tau _{(k)}\simeq 1/K$ for a given value of K. Therefore, $\tau _{(k)} \simeq 0.33, 0.067, 0.033$ for $K = 3,15,30$, respectively. We simulated 100 mixtures from each combination of ($K, p, \check{\omega }, \tau _{(1)}$) and one dataset of size $n=100 \times K$ from each mixture. The EM algorithm is run initialized by the competing methods M-EM, RndEM, emEM, SEM and emaEM. Adjusted Rand index values are calculated for each of the methods for all 100 datasets in each combination. Mean adjusted Rand index values for each combination ($K,p,\check{\omega },\tau _{(1)}$) are given in Fig. 11. Mean ranks (${\bar{\mathcal R}}$) based on the adjusted Rand index are presented in Fig. 12. Overall, the results are consistent with the ones found previously with approximately equal clusters. It can be noted that emaEM performed slightly worse than other competing methods when $K = 3$ and $\tau _{(1)}$ is small. This is shown in the first row of Figs. 11 and 12. However, when the number of clusters increases emaEM is superior to other methods. To make emaEM superior under all scenarios, the algorithm can be modified slightly to pick the maximum between the aggregated result and the best of random restarts based on likelihood values. Hence, this will guarantee that at least as good as an emEM solution to be found.

Fig. 12

Results of the simulation study for unequal cluster sizes. The x-axis represents the competing initialization methods and y-axis represents the mean rank achieved by each method. Each sub-plot has a description as in Fig. 11